# transient source term

 Register Blogs Members List Search Today's Posts Mark Forums Read

 March 13, 2011, 11:53 transient source term #1 New Member   Klaus Join Date: Jul 2010 Location: Linz / Austria Posts: 20 Rep Power: 7 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

 October 23, 2011, 22:47 #2 Member   Robin Gilbert Join Date: Jan 2010 Posts: 66 Rep Power: 7 does anybody know how to do this?? transient source term?

 October 24, 2011, 02:40 #3 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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

 October 24, 2011, 16:42 #4 Member   Robin Gilbert Join Date: Jan 2010 Posts: 66 Rep Power: 7 Hi, Thank you for your reply. This is what I want: Code: ``` 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: Code: ``` 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

 October 25, 2011, 02:50 #5 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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 Code: ```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

 October 25, 2011, 10:10 #6 Senior Member   David Boger Join Date: Mar 2009 Location: Penn State Applied Research Laboratory Posts: 146 Rep Power: 8 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. __________________ David A. Boger

 October 25, 2011, 10:20 #7 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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

October 25, 2011, 15:06
#8
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
May I add a bit more entropy to the discussion? Yes? Thanks

Quote:
 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

October 25, 2011, 16:18
#9
Member

Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 7
Quote:
 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:
Code:
`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.

Quote:
 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.

 October 26, 2011, 14:34 #10 Senior Member   Bernhard Linseisen Join Date: May 2010 Location: Magdeburg/Geneva Posts: 178 Blog Entries: 1 Rep Power: 7 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. ;-)

 October 26, 2011, 15:00 #11 Member   Robin Gilbert Join Date: Jan 2010 Posts: 66 Rep Power: 7 Hi Bernard, Thank you so much. I will look into it when you post the documentation.

 October 27, 2011, 04:01 #12 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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 Code: ```IOdictionary defineQPositionDict ( IOobject ( "defineQPositionDict", runTime.constant(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); PtrList regions(defineQPositionDict.lookup("regions")); Info << regions << endl; forAll(regions, regionI) { const entry& region = regions[regionI]; autoPtr 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 Code: `MUST_READ_IF_MODIFIED` . The dictionary is placed in constant and called Code: `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 Code: ```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 robingilbert, jiez and Tobi like this.

October 27, 2011, 07:20
#13
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
Quote:
 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)
Code:
```
expressionSource<scalar> heaterSource
(
IOdictionary
(
IOobject
(
"heaterSourceDict",
runTime.constant(),
mesh,
IOobject::NO_WRITE
)
),
mesh
);```
then modify the equation like this
Code:
```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)
Code:
```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
Code:
```functions {
theSensor {
type createSampledSet;
outputControl timeStep;
outputInterval 1;
setName TSensor;
setName {
type cloud;
axis x;
points (
(0 0 0) // core center
);
}
}
}

libs (
"libswakFunctionObjects.so"
);```
- 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?)

 October 27, 2011, 09:56 #14 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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

October 27, 2011, 14:52
#15
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
Quote:
 Originally Posted by kathrin_kissling Sorry for provoking you...
You're forgiven
Quote:
 Originally Posted by kathrin_kissling and for getting Aleister in the game...
That was me.
Quote:
 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:

Code:
```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

October 27, 2011, 16:08
#16
Member

Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 7
Quote:
 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.

Last edited by robingilbert; October 27, 2011 at 16:26. Reason: edit question

October 28, 2011, 02:06
#17
Senior Member

Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Quote:
 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.

 October 28, 2011, 02:19 #18 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 8 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

October 28, 2011, 03:40
#19
Member

Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 7
Quote:
 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 :

Code:
```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.

 October 28, 2011, 06:16 HowTos online #20 Senior Member   Bernhard Linseisen Join Date: May 2010 Location: Magdeburg/Geneva Posts: 178 Blog Entries: 1 Rep Power: 7 I guess some people were quicker than me. Fine as well! :-) Just for the sake of completeness: The HowTo mentioned above in transient source term 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... Last edited by Linse; November 3, 2011 at 18:42. Reason: Links changed...

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post wind Fluent UDF and Scheme Programming 5 June 21, 2013 05:39 Rajani Kanth.B Fluent UDF and Scheme Programming 4 May 1, 2013 09:31 BlGene Open Source Meshers: Gmsh, Netgen, CGNS, ... 10 August 6, 2009 04:26 segersson OpenFOAM Running, Solving & CFD 5 March 3, 2006 00:06 Brian FLUENT 1 October 24, 2005 09:15

All times are GMT -4. The time now is 17:28.