CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

transient source term

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

Like Tree12Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 13, 2011, 11:53
Default transient source term
  #1
New Member
 
Klaus
Join Date: Jul 2010
Location: Linz / Austria
Posts: 20
Rep Power: 15
strakakl is on a distinguished road
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
strakakl is offline   Reply With Quote

Old   October 23, 2011, 23:47
Default
  #2
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
does anybody know how to do this?? transient source term?
robingilbert is offline   Reply With Quote

Old   October 24, 2011, 03:40
Default
  #3
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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
kathrin_kissling is offline   Reply With Quote

Old   October 24, 2011, 17:42
Default
  #4
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
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
Kummi likes this.
robingilbert is offline   Reply With Quote

Old   October 25, 2011, 03:50
Default
  #5
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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
Kummi likes this.
kathrin_kissling is offline   Reply With Quote

Old   October 25, 2011, 11:10
Default
  #6
Senior Member
 
David Boger
Join Date: Mar 2009
Location: Penn State Applied Research Laboratory
Posts: 146
Rep Power: 17
boger is on a distinguished road
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.
Kummi likes this.
__________________
David A. Boger
boger is offline   Reply With Quote

Old   October 25, 2011, 11:20
Default
  #7
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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
Kummi likes this.
kathrin_kissling is offline   Reply With Quote

Old   October 25, 2011, 16:06
Default
  #8
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
gschaider will become famous soon enoughgschaider will become famous soon enough
May I add a bit more entropy to the discussion? Yes? Thanks

Quote:
Originally Posted by boger View Post
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
gschaider is offline   Reply With Quote

Old   October 25, 2011, 17:18
Default
  #9
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
Quote:
Originally Posted by kathrin_kissling View Post
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 View Post
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.
robingilbert is offline   Reply With Quote

Old   October 26, 2011, 15:34
Default
  #10
Senior Member
 
Bernhard Linseisen
Join Date: May 2010
Location: Heilbronn
Posts: 183
Blog Entries: 1
Rep Power: 15
Linse is on a distinguished road
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. ;-)
Linse is offline   Reply With Quote

Old   October 26, 2011, 16:00
Default
  #11
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
Hi Bernard,

Thank you so much. I will look into it when you post the documentation.
robingilbert is offline   Reply With Quote

Old   October 27, 2011, 05:01
Default
  #12
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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<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
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, Tobi and 1 others like this.
kathrin_kissling is offline   Reply With Quote

Old   October 27, 2011, 08:20
Default
  #13
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by kathrin_kissling View Post
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::MUST_READ,
                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"
);
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?)
gschaider is offline   Reply With Quote

Old   October 27, 2011, 10:56
Default
  #14
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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
kathrin_kissling is offline   Reply With Quote

Old   October 27, 2011, 15:52
Default
  #15
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by kathrin_kissling View Post
Sorry for provoking you...
You're forgiven
Quote:
Originally Posted by kathrin_kissling View Post
and for getting Aleister in the game...
That was me.
Quote:
Originally Posted by kathrin_kissling View Post
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
gschaider is offline   Reply With Quote

Old   October 27, 2011, 17:08
Default
  #16
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
Quote:
Originally Posted by gschaider View Post
- 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 17:26. Reason: edit question
robingilbert is offline   Reply With Quote

Old   October 28, 2011, 03:06
Default
  #17
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 21
Bernhard is on a distinguished road
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.
Bernhard is offline   Reply With Quote

Old   October 28, 2011, 03:19
Default
  #18
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 17
kathrin_kissling is on a distinguished road
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
kathrin_kissling is offline   Reply With Quote

Old   October 28, 2011, 04:40
Default
  #19
Member
 
Robin Gilbert
Join Date: Jan 2010
Posts: 66
Rep Power: 16
robingilbert is on a distinguished road
Quote:
Originally Posted by Bernhard View Post
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.
robingilbert is offline   Reply With Quote

Old   October 28, 2011, 07:16
Default HowTos online
  #20
Senior Member
 
Bernhard Linseisen
Join Date: May 2010
Location: Heilbronn
Posts: 183
Blog Entries: 1
Rep Power: 15
Linse is on a distinguished road
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...

Last edited by Linse; November 3, 2011 at 18:42. Reason: Links changed...
Linse is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem of SOURCE term gradient in UDS wind Fluent UDF and Scheme Programming 6 December 1, 2022 15:21
UDF source term Rajani Kanth.B Fluent UDF and Scheme Programming 4 May 1, 2013 10:31
[Gmsh] Compiling gmshFoam with OpenFOAM-1.5 BlGene OpenFOAM Meshing & Mesh Conversion 10 August 6, 2009 05:26
Adding a momentum source term segersson OpenFOAM Running, Solving & CFD 5 March 3, 2006 00:06
UDF Source Term Units? Brian FLUENT 1 October 24, 2005 10:15


All times are GMT -4. The time now is 15:35.