[strakakl]

hi foamers,

I wonder how to set a transient source term.
I want to simulate an electric heater band around a metal cylinder. Therefor i added a heat source to the chtMultiRegion solver and, in the first step, initalized it with a constant value. Results seems to be ok.
In the next step i want to use transient values for the heat source because I got set values for the heater band from an experiment which i want to use with the heat source term. But i don't know how to write them on the field of the heat source.

I tried to modify one of the timeVarying BCs in a way, that they also modify the internal field in the heater and not only the boundary field. By now, i had no sucess . I'm also not sure if this is the best way to do this.

Does anyone know how to impement this kind of transient source term?

Cheers, Klaus

[robingilbert]

does anybody know how to do this?? transient source term?

[kathrin_kissling]

Hi guys,

what are your exact plans?
Can you post an equation?
Do you want to have a real source term or a specific boundary condition?


Best

Kathrin

[robingilbert]

Hi,

Thank you for your reply. This is what I want:

        fvm::div(phi, T)
      - fvm::Sp(fvc::div(phi), T)
      - fvm::laplacian(kappaEff, T)
    ==
    Q/(rho0*Cp0)

where Q is dependent on temperature at a particular point in the solution domain. I thought of using something like:

    if (T at (x1,y1,z1) > Tref)
    {
        Q=Qref1;
    }
   else
   {
   Q= Qref2
   }

where Tref, Qref1 and Qref2 are read by the solver in the beginning. But my question is how to extract temperature at the position (x1,y1,z1) ?

I have been trying to figure that out from the forum. It would be great help if you can answer this question or just give me some hints.

Thank you,

Robin

[kathrin_kissling]

Hi Robin,

so you're searching for a spatial dependent source term.

To get you values: The easiest way is to grap them at the cell center:

Something like

dimensionedScalar Tref("Tref", dimensionSet(0,0,0,1,0,0,0), 350.);

dimensionedScalar Qref1("Qref1", dimensionSet(yourDimensionsHere), scalarValue1);
dimensionedScalar Qref2("Qref2", dimensionSet(yourDimensionsHere), scalarValue2);

forAll(T, cellI) //This loops over all cell centers
{
    if(T[cellI]>Tref)
    {
       Q = Qref1;
    }
    else
    {
        Q = Qref2;
    }
}

This at least will set your field Q;

Be aware that the convergence can be poor!

Best

Kathrin

[boger]

Actually, I think this earlier post by Kathrin is the one that Robin needs, or the post from Su Junwei that immediately precedes it. Both show methods for returning the cell index for a specified (x,y,z) location. Otherwise, I think Robin's coding was essentially correct, where "T at (x1,y1,z1)" would become T[labelOfClosestCell] or T[nearestCellIndex] according to the referenced links.

[kathrin_kissling]

Oh maybe I got something wrong. I thought that he would like to chance the complete Q field!

Robin: Do you want only one single position to determine you Q value, and that this Q value is then constant for the whole domain or do you want to have a field of different Q values each depending on the actual cell value?

Best

Kathrin

[gschaider]

May I add a bit more entropy to the discussion? Yes? Thanks

Originally Posted by boger

Actually, I think this earlier post by Kathrin is the one that Robin needs, or the post from Su Junwei that immediately precedes it. Both show methods for returning the cell index for a specified (x,y,z) location. Otherwise, I think Robin's coding was essentially correct, where "T at (x1,y1,z1)" would become T[labelOfClosestCell] or T[nearestCellIndex] according to the referenced links.

But be warned that this only works reliable for serial runs. For parallel runs with N CPUs there is a (N-1):N chance that the cell in question is on another processor and therefor the value would have to be sent to the others

For your first implementation I wouldn't bother with that too much but be warned that parallel runs will most likely produce "funny" (depends on your sense of humor) results

[robingilbert]

Originally Posted by kathrin_kissling

Robin: Do you want only one single position to determine you Q value, and that this Q value is then constant for the whole domain or do you want to have a field of different Q values each depending on the actual cell value?

Hi Kathrin, thank you for your reply. I want a single position to determine my Q value and this Q value is not for the whole domain but only for a small volume inside the domain. Now, the problem that I am facing is that I am not able to extract the temperature at that single position at run time. And for making it work in a small volume in the domain I thought of doing something like:

LHS= Q/(rho*Cp0)*N

where N is 1 in the required domain and 0 elsewhere (I can set that using setFields). Do you think that will work? or is it a messy solution?

I will go through the threads suggested by Boger and get back with results.

Originally Posted by gschaider

But be warned that this only works reliable for serial runs. For parallel runs with N CPUs there is a (N-1):N chance that the cell in question is on another processor and therefor the value would have to be sent to the others

For your first implementation I wouldn't bother with that too much but be warned that parallel runs will most likely produce "funny" (depends on your sense of humor) results

In that case I will initially start with serial runs before trying it out in parallel.

Thank you all for your replies.

[Linse]

I do not know if it works, but an approach MIGHT be to use abovementioned forAll()-function first to know if it should switch any value. Then you could implement this either by way of changing the N or by adding another source term, again via the forAll() method.
I did something like that recently for testing something and will document this still this week, so look into this thread by the end of the week and you should find the link to some documentation. ;-)

[robingilbert]

Hi Bernard,

Thank you so much. I will look into it when you post the documentation.

[kathrin_kissling]

Hi Robin,

I thought about a little more Foam like solution and I came up with this.
Why just don't use the setFields functionality within your executable. The following piece of code is an example on how this could be done. It uses cellSets to define your desired box/circle or whatever

IOdictionary defineQPositionDict
(
    IOobject
    (
        "defineQPositionDict",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
);

PtrList<entry> regions(defineQPositionDict.lookup("regions"));
Info << regions << endl;
forAll(regions, regionI)
{
   const entry& region = regions[regionI];
   
   autoPtr<topoSetSource> source = topoSetSource::New(region.keyword(), mesh, region.dict());
   

       cellSet selectedCellsSet
       (
         mesh,
     "cellSet",
     mesh.nCells()/10+1 //size estimate
       );
       
       source->applyToSet
       (
         topoSetSource::NEW,
     selectedCellsSet
       );
       
       labelList selection = selectedCellsSet.toc();
       forAll(selection, cellI)
       {
           alpha1[selection[cellI]] = 1.; //Replace this with your source term calculation
       }

}

Here it is coded to set an alpha value during runTime but you can perfectly put in your temperature dependent source term. Of course you need to define it somewhere in advance. You can put this snippet everywhere, either before the beginning of the time loop or within. Then you can even read the dictionary, if you modified it during your run. Just change the dictionary write option to

MUST_READ_IF_MODIFIED

. The dictionary is placed in constant and called

defineQPositionDict

it has the same entries like the setFieldsDict for the regions and you can use all the possibilities to set fields, which you can use with setFields

regions
(
    boxToCell
    {
        box (0.1 0.1 -1) (0.5 0.5 1);
    }
);

I tried to run it in parallel and it works for me. I can not say for certain it will for you.
@ Bernhard can you take some entropy out here? It would be so great, to know if this is a safer way to programm this!

Best

Kathrin

[gschaider]

Originally Posted by kathrin_kissling

I tried to run it in parallel and it works for me. I can not say for certain it will for you.
@ Bernhard can you take some entropy out here? It would be so great, to know if this is a safer way to programm this!

Best

Kathrin

I tried to avoid this, but I feel provoked. Kathrin, nowady half of my postings here end with the words swak4Foam, so don't act surprised.

If I had to do such a thing (mesh region heated depending on a sensor) I'd do the following:

Add something like this to createFields.H (or some other appropriate location)

    
expressionSource<scalar> heaterSource
    (
        IOdictionary
        (
            IOobject
            (
                "heaterSourceDict",
                runTime.constant(),
                mesh,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        ),
        mesh
    );

then modify the equation like this

fvm::div(phi, T)
      - fvm::Sp(fvc::div(phi), T)
      - fvm::laplacian(kappaEff, T)
    ==
      heaterSource()

Now I'm done with C++ and for the rest the motto of Aleister Crowley (no, I'm not into satanism - although C++ is close enough) applies "Do what Thou Wilt".
The heaterSourceDict would look something like this (all of this may contain syntax-errors. I'm doing it from memory)

variables (
    "TThres=666;"  // here comes Aleister
    "TProbes{set'TSensor}=average(T);"
    "QHeat=42;"
    "cp0=1;"
);
dimensions [0 0 0 0 0 0 0]; // This will have to fit the equation
expression "(zone(heaterZone) && ? TProbes<TThres) ? QHeat/rho : 0";

This assumes that there is a cellZone names heaterZone in your mesh. The last thing that remains to be modified is the controlDict to provide the sensor

functions {
  theSensor {
    type createSampledSet;
    outputControl timeStep;
    outputInterval 1;
    setName TSensor;
    setName {
       type cloud;
       axis x;
       points (
         (0 0 0) // core center
       );
    }
  }
}

libs (
   "libswakFunctionObjects.so"
);

The nice thing about this approach is that
- C++-coding is kept to a minimum
- it works in parallel (not swaks doing: sampledSets in OF)
- further case-specific modifications go where the ought to be: into the case
- you can do all kinds of cool stuff without modifying the solver: make the heater depend on the state of a patch, a sampled surface. Using storedVariables you can do state-based heating (once we're heating: keep heating. Hysterese-like switching. Time-dependence ....). Propagate the switching informations to boundary conditions via global variables and groovyBC

Bernhard

PS: Kathrins method of generating a cellSet looks fine to me and is quite flexible. Of course it depends on your application whether you want to keep the cellSet for later computations. BTW: even there you can sneak swak in through the back-door via the expressionToCell-topoSource (am I getting obnoxious?)

[kathrin_kissling]

Sorry for provoking you...
and for getting Aleister in the game...

Thank you very much Bernhard. I should use swak4Foam a lot more often though.

If swak4Foam would not exist and neither cellSet nor sampledSet and I would like to have a cell based solution (in fact I don't want to but maybe sometimes for something else) I would need Pstream functionality, wouldn't I?

Best
Kathrin

[gschaider]

Originally Posted by kathrin_kissling

Sorry for provoking you...

You're forgiven

Originally Posted by kathrin_kissling

and for getting Aleister in the game...

That was me.

Originally Posted by kathrin_kissling

Thank you very much Bernhard. I should use swak4Foam a lot more often though.

If swak4Foam would not exist and neither cellSet nor sampledSet and I would like to have a cell based solution (in fact I don't want to but maybe sometimes for something else) I would need Pstream functionality, wouldn't I?

You mean for distributing the temperature in a point to all processors? Depends on what you mean with Pstream-functionality: Sending around values through OPStreams/IPStreams or a simple reduce (which is based on that)

I THINK (untested and doing the API-calls from memory) this would do the trick:

scalar Tprobe=-1;
label cellID=mesh.findCell(coordinatesOfTheProbe_PickYourOwnName);
if(cellID>=0) { // is -1 if the point is not in the mesh == all other processors
   Tprobe=T[cellID];
}
reduce(Tprobe,maxOp<scalar>());
// now all the processors are on the same page

Of course this assumes that the temperature at the probe location is always positive which some people might think is a bit restrictive

Bernhard

[robingilbert]

Originally Posted by gschaider

- you can do all kinds of cool stuff without modifying the solver: make the heater depend on the state of a patch, a sampled surface. Using storedVariables you can do state-based heating (once we're heating: keep heating. Hysterese-like switching. Time-dependence ....). Propagate the switching informations to boundary conditions via global variables and groovyBC

Thank you Gschaider and Kathrin for your replies. Thats exactly what I want.

I have one more question: Do you think there is a way I can fix the temperature at a faceZone by removing heat from the room using the "heaterZone" as shown in the figure? Basically what I want is that the air flowing out through the faceZone should be at a fixed temperature.

Please forgive me if you think that I am asking way too many questions.

[Bernhard]

I have one more question: Do you think there is a way I can fix the temperature at a faceZone by removing heat from the room using the "heaterZone" as shown in the figure? Basically what I want is that the air flowing out through the faceZone should be at a fixed temperature.

Did you try a fixedValue boundary condition on that face? However, fixing the temperature of the outlet does not make sense physically, and I suppose fixing it as a boundary condition will give you nothing but troubles. I suppose convective heat transfer is larger than the heat diffusion in your case, so a fixed outlet temperature does not have a meaning in my opinion.

[kathrin_kissling]

Hello Bernhard,

thanks again very very much. Now I have something to dig into! Fun stuff for the weekend! Maybe I find out how to send a negative temperature if desired !

Best

Kathrin

[robingilbert]

Originally Posted by Bernhard

Did you try a fixedValue boundary condition on that face? However, fixing the temperature of the outlet does not make sense physically, and I suppose fixing it as a boundary condition will give you nothing but troubles. I suppose convective heat transfer is larger than the heat diffusion in your case, so a fixed outlet temperature does not have a meaning in my opinion.

Hi Bernhard, it is not a boundary face. Its an internal faceZone. I am sorry that the picture is not so clear. Its a duct with heat removal inside the duct so that the air coming out of the duct is of constant temperature that I can fix. I tried something like :

LHS = - (flowDirection & U)*Area*(T-Treq) 
(something like this. the actual setup is in my lab computer)

where "Area" is the area of the faceZone and "Treq" is the required temperature. But with this setup the temperature keeps on rising after sometime. I want to be able to fix the temperature at the faceZone or atleast the temperature of air coming out of the "heaterZone". I know I am doing something fundamentally wrong with the equation that I am using. Thanks for your help. Sorry to keep bothering you again and again.

[Linse]

I guess some people were quicker than me. Fine as well! :-)

Just for the sake of completeness: The HowTo mentioned above in http://www.cfd-online.com/Forums/ope...tml#post329593 Post #10 and some commented code snippets now are online in the OpenFOAM-Section of http://cern.ch/blinseis .

Got something wrong there, will have to reedit the documents I linked...