[LidaSa]

Hi experienced Foamers,

I'm an OpenFoam greenhorn and I just cannot find the way how to add generated power into my model.

I have a vertical cylindrical channel filled with molten salt (part of Molten Salt Nuclear Reactor). The initial temperature of the salt is 900K, the velocity of forces steady flow is (0 0 3.6) + turbulences - it works just fine (or at least it looks reasonable). The turbulences change the temperature a little bit.

Now I need to add power generated by nuclear fission in the salt and solve the temperature in each cell (I assumed constant rho).

I have the following file generated for the used polyMesh by an other code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      volPower;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -3 0 0 0 0];

internalField   nonuniform List<scalar>
1000
(
 6.23713E+08
 8.93512E+08
 8.92027E+08
 7.60773E+08
 7.56750E+08
 5.81814E+08
 7.42936E+08
 9.79091E+08
etc...
)
;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

and a map file which maps every cell to a separate zone (???).

I'm using chtMultiRegionSimpleFoam with only one region, which might be an overkill, I don't really know, just glad it works.


I understood that it can be done by fvOptions file, but I didn't find a way how to add different heat source to each cell.

My fvOptions looks like this:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
momentumSource
{
    type            meanVelocityForce;
    active          yes;

    meanVelocityForceCoeffs
    {
        selectionMode   all;

        fields          (U);
        Ubar            (0 0 3.6);
        relaxation      1.0;
    }
}

heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
    timeStart       0;
    duration        200;

    scalarSemiImplicitSourceCoeffs
    {
        selectionMode   all;
        volumeMode      specific;//absolute;
        injectionRateSuSp
        {
            e           ( 5.87000E+08 10);

        }
    }
}

// ************************************************************************* //

Where the heatSource is not what I want and I'm actually not sure what it does.

Can anyone help me with this?
Thank you
Lida

[Bloerb]

The fvOptions are easily found within the source code with some explanations. You might want to take a look at the .H files.
https://github.com/OpenFOAM/OpenFOAM.../src/fvOptions

A quick look into these makes it seem like this is not as easy though in this case. To start things:

heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
    timeStart       0;
    duration        200;

    scalarSemiImplicitSourceCoeffs
    {
        selectionMode   all;
        volumeMode      specific;//absolute;
        injectionRateSuSp
        {
            e           ( 5.87000E+08 10);

        }
    }
}

SelectionMode is where this source is active. Currently this is everywhere. You can limit this to a specific cellZone, which you can create with topoSet for example. The volumeMode means if your source is given as the absolute value or per volume (.../m³)
Now inside the brackets is the explicit part and implicit part of your source.

            e           ( Su Sp);

Or in other words source(e) = Su + Sp*e. This means that the following would be added to your equation:

    EEqn += Su + fvm::SuSp(Sp, e);

Adding a thousand of these for each cell seems excessive.
A quick glance at the options did not yield an easy result, though i haven't used these much. A coded source might be better if you are able to program it. There is an example inside the codedSource directory.

[Swagga5aur]

Hello Lida,
I was wondering what kind of heat source dependence of the geometry you want?

Here is a fvoptions file I used in a previous project as a axially variable heat source:

energySource
{
    type            scalarCodedSource;	//scalarSemiImplicitSource
    active          true;
    name	    sourceTime;

    scalarCodedSourceCoeffs	//scalarSemiImplicitSourceCoeffs    S(x) = Su + Sp*x  //   q in [W]; or in [W/m³] if you use specific mode
    {
        selectionMode   cellZone;
        cellZone        porousity1;
	fields 		(h);

	fieldNames	(h);
	name		sourceTime;

        codeInclude
        #{

        #};

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

        codeAddSup
        #{
//            const Time& time = mesh().time();
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
	    const volVectorField& U = mesh().lookupObject<volVectorField>("U");
	    const volScalarField z = U.mesh().C() & vector(0,0,1);

            scalarField& hSource = eqn.source();
            forAll(C, i)
            {
            	hSource[i] -= 6*((497530 - 530371*exp(-174.834*mag(0.311582803726196-z[i])))*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
    }
}

hope this helps and please let me know if you have any questions!

[LidaSa]

Hi Lasse,

My heat source is different in each cell and it cannot be solved by a physical equation. Or it could, by the point of my calculation is to use the volScalarField solution produced by a Monte Carlo calculation to solve the temperature. I hoped I could just somehow add it to the 0 directory as and additional internal energy field.

And I thing you gave me the solution

I have some questions, probably just because my C++ skills are very limited.

Does this:

const volVectorField& U = mesh().lookupObject<volVectorField>("U");

mean that U (as velocity?) is looked up from the time directory file?


In that case I could use something like this instead of the equation for heatSource definition if it would be stored in volPower file in 0 directory?

scalarField& hSource = mesh().lookupObject<volScalarField>("volPower");


Is this the dictionary of cell volumes?
scalarField& V = mesh_.V();

What does this mean? Is the & like a scalar product of velocity and (0,0,1) vectors?
const volScalarField z = U.mesh().C() & vector(0,0,1);

Thank you so much
Lida

[Swagga5aur]

I'll try to answer the best I can, however, I'm no expert myself when it comes to C++.

const volVectorField& U = mesh().lookupObject<volVectorField>("U");

The volVectorField& U is the velocity in the individual cells coinciding in cellZone porousity1.

scalarField& V = mesh_.V();

Yes it is the volume of an individual cell in the cellZone.

const volScalarField z = U.mesh().C() & vector(0,0,1);

In my code the point was to have an axially varying heat source (exponentially rising until an upper limit), where z is the axial coordinate, as the U.mesh().C() fetches the cell center with & vector(0,0,1) being the scalar product, resulting in the axial cell center coordinate.

I haven't tried adding an additional energy field in the 0 directory so won't be to much help there sorry.

Let me know if there is any other questions or if I missed something.

Additional links I used to develop the code in my previous post:
adding a constant volumetric source term to transport equation in a particular region
Howto use scalarCodedSource in fvOptions
How to Define heat generation rate that changes with the system coordinate.

[LidaSa]

Thank you so much. I just keep running into this warning:

Source energySource defined for field h but never used

Do you have any idea what do I have to do to make it work with chtMultiRegionSimpleFoam?

[Swagga5aur]

I implemented the heatsource code you posted initially in the heatExchanger example from chtMultiRegionSimplefoam and got the same error code, the following code should work:

heatSource
{
     type scalarSemiImplicitSource;
     active true;
     
     scalarSemiImplicitSourceCoeffs
     {
	 volumeMode specific;
         selectionMode all;
	 injectionRateSuSp
	 {
		h (5.87000E+08 10);
	 }
     }
}

I believe the issue is that your thermophysicalProperties uses the energy model: energy sensibleEnthalpy;
This requires the fvOptions to be in enthalpy h and not internal energy e, as you defined it in the first post.

Best regards Lasse.

[LidaSa]

Yes yes yes, I realized that the moment I posted it. Sorry about that, sensibleEnthalpy it is.

But this is not it.

scalarField& volPower = mesh().lookupObject<volScalarField>("volPower");

It runs into this error:

request for volScalarField volPower from objectRegistry salt failed
available objects of type volScalarField are

16
(
thermo:mu
thermosi
nut
h
rho
k
p_rgh
gh
momentumSource:rA
(C&(0,0,1))
alphat
p
T
thermo:rho
epsilon
thermo:alpha
)


I now have fvOptions basically copied from you, defining a heat source depending on z coordinate of each cell.
hSource[i] -= 100000000*(z[i]+1);
Which is a great improvement, at least I can define it separately for each cell and it works.


What I need is something like this:

codeAddSup
        #{

            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
	    
            const volScalarField volPower = mesh().lookupObject<volScalarField>("volPower");


            scalarField& hSource = eqn.source();
            forAll(C, i)
            {
            	hSource[i] -= volPower[i]*V[i];
            }
            //Pout << "***codeAddSup***" << endl;
        #};

In other words, I need to load the volPower file from my first post to a dictionary assigning a volPower to each cell in order to work with forAll(C, i) in the equation. I guess I'm gonna need to refresh my C++ skills.

[Swagga5aur]

I took the liberty to try and implement what it seems you want, into a new chtMultiRegionSimpleFoam solver.
I have attached the code and an example of how to setup the case.
I don't have the time to thoroughly test it, but it should be a decent starting point to your solution.
https://www.dropbox.com/s/g8fz1xmxtf...ission.gz?dl=0

Best regards Lasse.

[LidaSa]

Oh Lasse, thank you so much. You just saved my master thesis (due friday...). I'm so grateful

Lida

[Swagga5aur]

No problem, let me know if any other issues should occur and I'll try to answer as soon as possible.

[shaileshbg]

Hi Lasse,

I am also working on adding heat source to every cell in chtMultiRegionSimpleFoam, could you please share the dropbox link again, it doesn't work now. That would be very helpful to me.

[Swagga5aur]

Hi Shailesh,

Sorry for the inconvenience I have re uploaded at the following link:

https://www.dropbox.com/s/4dxbyvk5a2...on.tar.gz?dl=0

Let me know if there is anything else I can help with.

Regards Lasse

[shaileshbg]

Hi Lasse,

Thank you very much for your reply and for re-uploading the files.

My problem is very similar to the one you described in post #3. I need to add a heat source to each cell based on its temperature at the cell centre, so my lookup object is temperature.
The heat source description is:

Heat Source = constant1* (Temperature (at cell centre) - constant2)

The fvOptions file that I have defined is:

            const vectorField& C = mesh_.C();                    //vector of cellcentres
	   
	    const volScalarField& T = mesh().lookupObject<volScalarField>("T"); // Look up temperature 
	    const volScalarField Tcc = T.mesh().C() & vector(1,1,1);  // Temperature at each cell centre, & vector used to convert from volVectorField to volScalarField
	    scalarField& hSource = eqn.source();                //defining source
	    	    
            forAll(C, i)
	    {
		
		hSource[i] = (Tcc[i] - 310)*2400 ; // -c_b w_b (T - Ta) = h = constant * (T (temp at each cell centre) - 310 (wall temp))
            }

With this I am getting unusually high-temperature values, I have attached the temperature graph alongside. The temperature range that I am expecting is in between 300K to 315K

Could you please help me know where my mistake is.

Thanking You,

[Swagga5aur]

Hello Shailesh,

I would like to see some more information about your case, but it seems to me that there is no limitation on the heat source as it increases with the temperature which increases the source. You could try to implement a simple source to check your definitions.

Feel free to send me your case and I'll give it a look.

Additionally, I noticed that you have hSource[i] = and in my example I use -= being different, but I'm not sure how your source should be defined, just wanted to let you know.

Regards Lasse

[shaileshbg]

Hello Lasse,

Thank you very much for your reply.

Please find the dropbox link to my case.

https://www.dropbox.com/s/pel9rkq2db...reast.zip?dl=0

Thank you for taking the time to look at my case, I look forward to hearing back from you.

[Swagga5aur]

Hello Shaileshbg,

Just to make clear this theory is outside of my field, but I'll give it a go with helping you.

I believe the issue lies in the heat source addition, its independent of the cell size so an infinitely small cell would have the same heat source ie infinite temperature as a very poor mesh. Tried adding the volume as in the following:

            const scalarField& V = mesh_.V();
	    const volScalarField& T = mesh().lookupObject<volScalarField>("T"); // Look up temperature 
            const volScalarField Tcc = T.mesh().C() & vector(1,1,1);  // Temperature at each cell centre, & vector used to convert from volVectorField to volScalarField
            const vectorField& C = mesh_.C();                    //vector of cellcentres

            scalarField& hSource = eqn.source();                //defining source
            forAll(C, i)
            {
                 hSource[i] -= -(Tcc[i] - 310)*48000*V[i] ; // c_b w_b (T - Ta) = h = constant * (T (temp at each cell centre) - 310 (wall temp))
            }
            Pout << "***codeAddSup***" << endl;

Reducing the temperature to about 3000K, still being far too much.
However, I also suspect that the constant value "48000" is the product of the blood density ~constant, blood heat capacity ~constant and the bloodperfusion which seems to be a unit of volume flow / volume which may be considered constant resulting in a unit of [W/m^3].

Saw some litterature defining the above as all constants with bloodperfusion of 0.0005 and I assume that you have a bloodperfusion of 0.011429 for the 48000 constant.

Without knowing more about the physical problem and general theory of the field I won't be of more assistance, and as of now I won't have time to investigate it further but feel free to link me some material regarding the issue and and I'll try find some time to help you in the next week.

Regards Lasse.

[priyankap]

Hello Lasse,


I have a very basic question.



I have a box of dimension say (300m*300m*10m). Now in this box, can I define just the front face of the box as a heat source using the fvOptions file or not?



If yes, how can I do it? Because I know if I want to define this entire box as a heat source using power then my fvOptions file should be something like this:

heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
 
    scalarSemiImplicitSourceCoeffs
    {
        selectionMode   all; // all, cellSet, cellZone, points
       // cellZone        heatSrc;
        //cellSet         c1;
        volumeMode      absolute; // specific;
        injectionRateSuSp
        {
            h     (10 0);
        }
    }
}

but now, if I want to define just one side of the box as heat source what do I do?


Thank you
Priyanka

[shaileshbg]

Hi,

You can use the selectionMode option to specify the source region.

selectionMode   cellSet; // all, cellSet, cellZone, points

But there is no option to set faces or boundary patches. To overcome this, we can use topoSet to select all the faces in the patch (front face) to form a faceSet and then use that faceSet to form a cellSet.

First, define your patch ("front", in this case), easiest is to set it in blockMeshDict.

Then in your topoSetDict, use that patch info to form your faceSet.

{
        name    frontFaces;
        type    faceSet;
        action  new;
        source  patchToFace;
        sourceInfo
        {
            name "front";
        }
    }

And then use this faceSet to define your cellSet

// Select based on faceSet

    name    frontFaceCellSet;
    type      cellSet;
    action    new;
    source   faceToCell;
    sourceInfo
   {
       set frontFaces;             // Name of faceSet

       //option neighbour; // cell with neighbour in faceSet
        //option owner;     //  ,,       owner
        option any;         // cell with any face in faceSet
       //option all;       // cell with all faces in faceSet
    }

Use this cellSet in the selectionMode while defining the source.

selectionMode   cellSet; // all, cellSet, cellZone, points
        cellSet        frontFaceCellSet;

Hope this helps.

[priyankap]

Thanks a lot for your answer. I will try to do it this way. I have one more question.


In my previous post I had shared my fvOptions file which was:

heatSource
{
    type            scalarSemiImplicitSource;
    active          true;
 
    scalarSemiImplicitSourceCoeffs
    {
        selectionMode   all; // all, cellSet, cellZone, points
       // cellZone        heatSrc;
        //cellSet         c1;
        volumeMode      absolute; // specific;
        injectionRateSuSp
        {
            h     (1000 0);
        }
    }
}

Here what is 'h'? I am putting my power value here to define my heat source, which is 1000W. If I give this value to my heat source it generates some temperature on the region which I defined as the heat source.



Now I want to know how is this temperature calculated? I know the formula



Q = (k.A. (delta T))/L

where Q is the power in W, k is thermal conductivity on W/m.K, A is the area if cross section, delta T is temperature gradient and L is the length.



But when I use this formula in my simple case of just two plates (one of which is made heat source and the other is at room temperature) joined with a rod, I cannot verify the power value that I have given in 'h' of fvOptions file.


Now, I don't understand if this 'h' is directly power or it is something else which in turn calculates the temperature of the heat source.