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

Howto use scalarCodedSource in fvOptions

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

Like Tree37Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 18, 2013, 09:59
Default Howto use scalarCodedSource in fvOptions
  #1
Member
 
hannes
Join Date: Mar 2013
Posts: 47
Rep Power: 13
hanness is on a distinguished road
Hi,

currently I'm struggling with setting up a volumetric source. In my case I would like to add a volumetric heat source for a given cellSet and run the simulation with the buoyantBoussinesqSimpleFoam solver, later with the respective Pimple solver.
I tried creating a fvOptions file but unfortunately I don't really know what to put where... I adapted the below code from the CodedSource.H and now my question is where can I put the source term and in which form?
Code:
mySource
{
    type            scalarCodedSource;
    active          on;
    selectionMode   cellSet;
    cellSet         myCellSet;
    rampCoeffs      {};
    scalarCodedSourceCoeffs
    {
        fieldNames      (T);
        redirectType    ramp;
        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};
        codeAddSup
        #{
            Pout<< "**codeAddSup**" << endl;
        #};
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};
        code
        #{
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }
}
Moreover, I see a general problem with this approach because I don't actually want to set a source for the temperature but for the heat (I know the value in terms of W/m³) and the Boussinesq-solvers solve a T-equation and not a h/e-equation.
Does anybody have any suggestions?

Thanks a lot
Hannes
hanness is offline   Reply With Quote

Old   November 29, 2013, 08:49
Default
  #2
New Member
 
RB
Join Date: Aug 2013
Posts: 5
Rep Power: 12
rbaud is on a distinguished road
Hi Hanness,

I'm interesting in your case, do you finally work out how to tackle your problem?

(i.e volumetric "heat source" in buoyantBoussinesqSimple/PimpleFoam)

Why are you trying with a codedSource, an explicit or semi-implicit should be able to do the job, shouldn't it?

I'm currently able to add a source of "temperature", (i.e I set a value with an explicitSource and the temperature raise up in my selected region) but I'm struggling to understand which amount of heat I'm setting in my model.

Is it wrong to try to link it to a proper heat load? (as we are solving only temperature?)

I'm currently trying the same basic "testing" than in this post:
http://www.cfd-online.com/Forums/ope...implefoam.html

And it doesn't make sense. (does the source's unit is really in K/s? --> seems ok for me regarding the equation i.e homogeneous with the dT/dt)

Thanks for sharing any hints,
R
rbaud is offline   Reply With Quote

Old   December 6, 2013, 05:25
Default
  #3
Member
 
hannes
Join Date: Mar 2013
Posts: 47
Rep Power: 13
hanness is on a distinguished road
Hi R,

Eventually I got it working. I don't know whether it is the best approach but it works for me:

Code:
heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        heatSourceCells;
    scalarSemiImplicitSourceCoeffs
    {
        volumeMode      specific;
        injectionRateSuSp
        {
            T           (0.1 0); //DT = PHI''' / (rho * cp) 
        }
    }
}
Maybe my comment behind the source term might help you to calculate the correct value for DT?

Regards,
Hannes
student666 and huadian001 like this.
hanness is offline   Reply With Quote

Old   December 9, 2013, 11:23
Red face volumetric heat source in TEqn, dimension understanding
  #4
New Member
 
RB
Join Date: Aug 2013
Posts: 5
Rep Power: 12
rbaud is on a distinguished road
Hi Hannes,

Many thanks for your reply,
I use the same functionality as you to get it run (well), (i.e. semiimplicitsource)

I've took the week end to sleep on it and reach the following:

I like the idea of your: Q/(rho*cp) = Ts (1)
It fits with the TEqn solved, which is, I think an energy equation divided by rho*cp, considered as constant.
With:
Ts : source of temperature (0.1 in your example).
Q : heat source
Is that correct with your notation, I try to avoid any confusion with your PHI'''.

I’m still kind of struggle with the dimension of all of this.
Based on the attached small test case (air flow in a duct with a heat source in the middle of the duct.), I get the relation:
Ts = phi*DT
With:
phi : (openfoam calculated) flow (m3/s as we have an incompressible solver)???
DT variation of temperature along the duct (K)

Coupled with the (1), it fits with the classic Q=phi*rho*cp*DT.

I could stop there and be happy but 2 questions remained:
- Which value are you using for your rho and cp, as in the transport property we don't set them explicitly but hidden under the Pr and nu. The solver probably don't use rho and cp at all, so I'm not confident with this step. (do you just use the value which fit the model environment, without any link with the model?)

- This approach leads for the Ts a dimension in m3.K/s where I think K/s is more correct, regarding the TEqn. I guess somewhere the source module link it to the volume to reduce it in K/s. ???

(I’m aware of the absolute/specific functionality, and I’ve understand what is doing in the fact but I’ve not fully understand the code), and so can’t explain my dimension issue..

I hope my post is clear,
I’m probably doing a mistake somewhere or missing something, certainly trivial…
Thanks for any reply,
Remi
Attached Files
File Type: zip T_Source_2_velocity_duct.zip (63.0 KB, 238 views)
huadian001 likes this.
rbaud is offline   Reply With Quote

Old   March 20, 2014, 13:59
Default
  #5
New Member
 
Join Date: Mar 2012
Posts: 29
Rep Power: 14
yanxiang is on a distinguished road
Hi Hannes,

Thanks for the comments first of all.

It looks like "specific", instead of "absolute", is required for this to work properly. However, it is kind of confusing to me when it is named "specific" and "absolute". To me, it sounds more like "intensive" property versus "extensive" because when I calculated the source for T, I didn't really have a m^3 in the denominator.

yanxiang
yanxiang is offline   Reply With Quote

Old   May 6, 2014, 15:36
Default Volumetric heat source (W/m³)
  #6
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 337
Rep Power: 21
zfaraday will become famous soon enough
Hello everybody!

I'm aware that this thread is maybe dead by now but I have a quick question related to what it has been said here. I hope that someone can answer me.

My case is that I have to define a volumetric heat source in chtMultiRegionSimpleFoam in a solid region defined as, say, heatSource. The approach I thought about is similar to what suggested hanness but I don't want to define a "source of temperature" (K) but a "source of heat" (W/m³). I have tried the following fvOptions settings to do so:

Code:
heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
    selectionMode   all;

    scalarSemiImplicitSourceCoeffs
    {
        volumeMode      specific;//absolute;
        injectionRateSuSp
        {
            h           (X Y);
        }
    }
}
I have tried a few combinations of Su-Sp coefficients but I don't really understand the meaning and use of them. What values should I give to the coeffs if I want to set up a source of, for instance, 50 W/m³?

Another question I have is if it's correct or not to use the h variable (enthalpy, if I am not wrong) to create a heat source in a solid region.

I have also made some other attempts using T instead of h but I noticed that nothing happened, like if no source existed. But, according to what I read in this topic above, it shoud have been because I used absolute mode instead of specific. Am I right?

Thanks in advance. Any suggestion or hint will be welcome!


Alex
__________________
Web site where I present my Master's Thesis: foamingtime.wordpress.com

The case I talk about in this site was solved with chtMultiRegionSimpleFoam solver and involves radiation. Some basic tutorials are also resolved step by step in the web. If you are interested in these matters, you are invited to come in!
zfaraday is offline   Reply With Quote

Old   May 6, 2014, 21:33
Default
  #7
Senior Member
 
Andrew Somorjai
Join Date: May 2013
Posts: 175
Rep Power: 12
massive_turbulence is on a distinguished road
Quote:

My case is that I have to define a volumetric heat source in chtMultiRegionSimpleFoam in a solid region defined as, say, heatSource. The approach I thought about is similar to what suggested hanness but I don't want to define a "source of temperature" (K) but a "source of heat" (W/m³). I have tried the following fvOptions settings to do so:
I'm almost positive this would work with the chtMultiRegionSimpleFoam using the heat flux in W/m³.

Code:
        type            compressible::turbulentHeatFluxTemperature;
        heatSource      flux;
        q               uniform 1000;
        alphaEff        kappaEff;
        K               basicThermo;
        Cp              uniform 1000;
        value           uniform 273.5;//set initially?
        KName             K;
I'm assuming that you know the initial temperature of your entire system which means in this case it would be set to 273.5k.
massive_turbulence is offline   Reply With Quote

Old   May 7, 2014, 13:37
Default
  #8
Senior Member
 
Alex
Join Date: Oct 2013
Posts: 337
Rep Power: 21
zfaraday will become famous soon enough
Quote:
Originally Posted by massive_turbulence View Post
I'm almost positive this would work with the chtMultiRegionSimpleFoam using the heat flux in W/m³.
Thanks for your quick replay Andrew!

Yes, I guess it should work at least this configuration of the fvOptions file makes sense to me. But my hesitation is about the values I have to set as Su-Sp coefficients (assuming the configuration is ok) to define the heat generation (for example, 50 W/m³) inside the region.

Regards

Alex
__________________
Web site where I present my Master's Thesis: foamingtime.wordpress.com

The case I talk about in this site was solved with chtMultiRegionSimpleFoam solver and involves radiation. Some basic tutorials are also resolved step by step in the web. If you are interested in these matters, you are invited to come in!

Last edited by zfaraday; May 7, 2014 at 15:59.
zfaraday is offline   Reply With Quote

Old   December 1, 2014, 16:07
Default
  #9
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
Hi,Alex.
Did you figure out what X and Y stand for? I know that they basically describe a source as s=X+Yx where the small x is coordinate. Is that right?

Thanks
Zech
a19910112a is offline   Reply With Quote

Old   December 1, 2014, 16:32
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

you can find an answer in the SemiImplicitSource.H:

Code:
Description
    Semi-implicit source, described using an input dictionary.  The injection
    rate coefficients are specified as pairs of Su-Sp coefficients, i.e.

        \f[
            S(x) = S_u + S_p x
        \f]

    where
    \vartable
        S(x)    | net source for field 'x'
        S_u     | explicit source contribution
        S_p     | linearised implicit contribution
    \endvartable
or if you prefer code, SemiImplicitSource.C

Code:
    UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI].first()/VDash_;
     
    DimensionedField<scalar, volMesh> Sp
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "Sp",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensioned<scalar>
        (
            "zero",
            Su.dimensions()/psi.dimensions(),
            0.0
        ),
        false
    );

    UIndirectList<scalar>(Sp, cells_) = injectionRate_[fieldI].second()/VDash_;

    eqn += Su + fvm::SuSp(Sp, psi);
I.e. x is a field for which source is calculated.
alexeym is offline   Reply With Quote

Old   December 2, 2014, 06:50
Default
  #11
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
yeah, thanks for answering. I have looked at that some time ago. But that doesn't really make sense to me. As I don't really know how to express my heat resource in an implicit and explicit form.
1. Any suggestions on where can I learn about that?
2. What I want to define is a volumetric heat generation rate that changes along x axis, rather than with a field (heat generation rate q"'=cos (x) where x is the system coordinate). Any idea on how to define that?
3. A more challenging question is: what if the source is not even a regular function. It is just a bunch of random numbers. Say a table like
when x=1 q"'=4
when x=2 q"'=7
when x=3 q"'=1 ..... irregular (x is still system coordinate and as assume the q"' doesn't change along y and z direction. )

wishes
Zech
a19910112a is offline   Reply With Quote

Old   December 2, 2014, 07:01
Default
  #12
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

1. Sources. Tutorial. Google.

2. Lots of. Usually at this point people suggest swak4Foam. Or you can modify solver. Or you can take SemiImplicitSource and modify it.

3. In fvOptions file.
Tushar@cfd likes this.

Last edited by alexeym; December 3, 2014 at 07:01.
alexeym is offline   Reply With Quote

Old   December 2, 2014, 07:40
Default
  #13
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
1. that means nothing. every body doing Foam knows....
3. can you say something more? How? If you can do this, question 2 is not even question !!!
a19910112a is offline   Reply With Quote

Old   December 2, 2014, 09:24
Default
  #14
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Well,

If you need someone to Google for you, there is http://www.esi-cfd.com/content/blogcategory/99/128/, you know.

For your case you can take a look at CodedSource. In this case you can just use C++ to extract coordinates from field and to calculate your source term.
Tobi likes this.
alexeym is offline   Reply With Quote

Old   December 3, 2014, 07:36
Default
  #15
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
Thank again for your reply.

Yes,I'm looking at the codedSource. That is exactly what this thread start with discussing, isn't it? Do you have any examples of this? If you know how to use it can you post it here, which, I believe, will be helpful for a lot of people.
I'm asking because I have exhausted the searching result from Google without finding some answers and could not understand what the source code is doing as many new comers to openFoam cannot.
a19910112a is offline   Reply With Quote

Old   December 3, 2014, 11:12
Default
  #16
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Guess, you're right. Google search leads to the first message of this thread. There were other links though with less information.

So here's an example of the codedSource:

Code:
harmonic
{
    type            scalarCodedSource;

    active          true;
    selectionMode   all;

    scalarCodedSourceCoeffs
    {
        fieldNames      (T);
        redirectType    sourceTime;

        codeInclude
        #{

        #};

        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};

        codeAddSup
        #{
            const Time& time = mesh().time();
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
            scalarField& TSource = eqn.source();
            forAll(C, i)
            {
                TSource[i] -= 10*sin(6.28*C[i].x())*V[i];
            }
            Pout << "***codeAddSup***" << endl;
        #};

        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};

        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }

    sourceTimeCoeffs
    {
        // Dummy entry
    }
}
As the configuration is more or less taken from CodedSource.H file (well, sin addition is mine), I don't know what questions this will arise. So I skip commenting the configuration and wait till further questions.

Also I've attached archive of the case (hotRoom case from buoyantBoussinesqSimpleFoam tutorial with removed gravitation (therefore no convection) and suppressed heat diffusivity) and a screen shot of temperature distribution.
Attached Images
File Type: png T.png (26.7 KB, 236 views)
Attached Files
File Type: gz hotRoom.tar.gz (3.2 KB, 258 views)
alexeym is offline   Reply With Quote

Old   December 3, 2014, 13:18
Default
  #17
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
Thank you so much for sharing the example!

I haven't get a chance to run the case yet, but it might be helpful anyway to explain what kind of heat source you have defined here. what it depends on? I mean I see it depends on C[i], V[i] and x(). What are they? I roughly know C and V are values extract from mesh, but not quite sure what are they. And x() is not even defined in advance, is it recognised as the system coordinate by default?

Comparing with the example in the codedSource.H file, you kept the codeCorrect and codeSetValue as they were. I assume their functions are irrelevant to the source your are defining here. If you happen to know, can you please say something on what they can be used to define?
a19910112a is offline   Reply With Quote

Old   December 4, 2014, 04:30
Default
  #18
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Well, at this point your answers for the questions are can be easily found in documentation (http://openfoam.org/docs/cpp/). I will answer some of them, the rest - it is up to you to look up.

You've deduced that V and C are properties of the mesh, so you enter mesh in search dialog in documentation (well, I will just enter fvMesh), click on fvMesh in drop-down menu and go to fvMesh documentation. There we'll find that V is a list of cell volumes, C is a list of cell centers (or centres if you prefer). Also from documentation you've learned that C is a list of vectors. Enter vector in search box, click on the link in the drop-down list, et voilà, x is a method of vector class and it returns x-component of the vector. As there's no convection, no diffusion, and no heat exchange with outer space, strange steady state with negative temperatures reached in 44 iterations.

As it was just a test, I've taken mesh with sizes 2x1x1, and source term in T-equation that is proportional to sin(2*pi*x), so in x-direction we have 2 periods (and there are two periods on temperature field).

As to other methods. Take a look at TEqn.H, relevant part is:

Code:
...
    fvScalarMatrix TEqn
    (
        fvm::div(phi, T)
      - fvm::laplacian(alphaEff, T)
     ==
        fvOptions(T)
    );

    TEqn.relax();

    fvOptions.constrain(TEqn);

    TEqn.solve();

    fvOptions.correct(T);
...
then take a look at solver's log:

Code:
DILUPBiCG:  Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Uz, Initial residual = 0, Final residual = 0, No Iterations 0
***codeAddSup***
**codeSetValue**
DILUPBiCG:  Solving for T, Initial residual = 0.331723, Final residual = 0.00127055, No Iterations 1
**codeCorrect**
DICPCG:  Solving for p_rgh, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 0.36 s  ClockTime = 7 s
So, AddSup is called during construction of the linear system matrix and it is just a source term in temperature equation, SetValue is called before solution of the linear system and can be imagined as a constraint like BC (well, it is called fvOptions.constrain(TEqn)), and finally correct is called after solution of the linear system and can be used to modify T field (well, in fact you can modify almost anything as you've got access to the volScalarField, from where you can get mesh reference, and then other fields, runtime object etc). Surely the same conclusions can be achieved by analysing of the sources of fvOptions framework.
ptbs, PonchO, sharonyue and 2 others like this.
alexeym is offline   Reply With Quote

Old   December 5, 2014, 07:27
Default
  #19
New Member
 
zech
Join Date: Oct 2014
Location: Cambridge,England
Posts: 22
Rep Power: 11
a19910112a is on a distinguished road
Thanks for your answer Alex. I have run your case and don't understand the flowing things:
1. You have set it to stop running after steady state is arrived (44s) somewhere? I wanted to see what happens as time goes, so I set it to run between 0-1000s in the controlDict. But it still runs from 0-44? What is the problem?

2. Because I cannot play with time, I played with the function in the codeAddSup part. I'm so sorry. Although you have told me so clear what V C x are, I'm still not quite sure about what how it defines the source.
Quote:
forAll(C, i)
{
TSource[i] -= 10*sin(6.28*C[i].x())*V[i];
}
a. One of the problem is, what is the value of C ? ‘C is a list of centers’, what is the value of it? I tried to look at the link you gave me on cpp. But I could not even tell ‘C is a list of centers’ from the page.

b. V is the volume, then the quantity 10*sin(6.28*C[i].x()) gives should be q’’’/(rho*Cp). Is that right?

3. As I cannot figure out the previous questions. I thought that I at least know the example you gave defined a temperature source with a period of 1. So, I applied the source into one that I’m familiar with. That is the planeWall2D case. I changed the bottom air with a heat generation solid, and it worked fine as a constant source. I replaced the source with the 1-period source in your example. I changed the proportion to observe a more obvious phenomenon. But sadly, no period could be observed. I have the case uploaded here.
a19910112a is offline   Reply With Quote

Old   December 5, 2014, 08:58
Default
  #20
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Quote:
Originally Posted by a19910112a View Post
1. You have set it to stop running after steady state is arrived (44s) somewhere? I wanted to see what happens as time goes, so I set it to run between 0-1000s in the controlDict. But it still runs from 0-44? What is the problem?
It's convergence criterion set in SIMPLE dictionary in fvSolution file. So you can set endTime to any value greater the 44, if solution converges is 44 iterations it will run for 44 iterations. Also in SIMPLE it's not time, it's iterations.

Quote:
a. One of the problem is, what is the value of C ? ‘C is a list of centers’, what is the value of it? I tried to look at the link you gave me on cpp. But I could not even tell ‘C is a list of centers’ from the page.
Just of curiosity I've clicked on the link I've posted and here's a screen shot of what I've got. To make it easier I've underlined a line you need to read. See attached image.

Quote:
b. V is the volume, then the quantity 10*sin(6.28*C[i].x()) gives should be q’’’/(rho*Cp). Is that right?
Not quite sure I've counted number of primes correctly but it should be power over rhoCp.

Concerning your case: maybe there's no period, I don't know.
Attached Images
File Type: png cell-centers.png (17.8 KB, 282 views)
alexeym is offline   Reply With Quote

Reply

Tags
fvoptions, heat source, scalarcodedsource

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
OpenFoam 2.2 fvOptions temperature limits fredo490 OpenFOAM Running, Solving & CFD 12 January 17, 2020 00:59
Building a solver with fixedTemperatureConstraint using fvOptions Fluido OpenFOAM Programming & Development 9 February 15, 2018 01:30
How to set fvOptions yurifrey OpenFOAM Pre-Processing 5 February 22, 2016 19:14
[swak4Foam] Setting BC for a passive scalar (groovy vs fvOptions) Tobi OpenFOAM Community Contributions 0 May 23, 2013 15:53
A new Howto on the OpenFOAM Wiki Compiling OpenFOAM under Unix mbeaudoin OpenFOAM Installation 2 April 28, 2006 09:54


All times are GMT -4. The time now is 07:52.