Adding temperature dependent heat source using fvOptions in chtMultiRegionSimpleFoam

shaileshbg · December 19, 2017, 10:33

Hello Foamers,

I want to add a temperature dependent heat source into the chtMultiRegionSimpleFoam solver using fvOptions.
I was looking at semiImplicitSource option which works on source term linearization using SuSp injection rate treatment where:
S(x) = Su + Sp(x)
S(x): net source for field (x)
Su: explicit source contribution
Sp: linearized implicit contribution

I understand that the solver works by evaluating enthalpy (so my field is H) and that the unit of enthalpy is same as that of internal energy (J/Kg), so if I want to add a constant heat source I need not convert anything I can put it in directly under Su (taking absolute and specific into consideration), but I do not know how to add a temperature dependent source term of this form:

source = constant * (T - T_const)

Problem definition:
Heat conduction through four different solids each having its own temperature-dependent source.

Could anyone please help me to understand how to implement this.
Many thanks in advance,
Shailesh

TobiF · December 19, 2017, 14:50

Hey,
use scalarCodedSource to calculate the temperature-depending heat-source.

If you need the mean-temp of the region you should define the selectionMode type cellZone and in the code something like :

Code:

const volScalarField& Tm = mesh_.lookupObject<volScalarField>("T");    
Tvol = Tm.weightedAverage(mesh_.V()).value();                //averageValue of the volScalarField

to calculate the mean tmep of the cellZone.

After that you can use the Tvol - value for your calculation.

Code:

            const vectorField& C = mesh_.C();                    //List of cellcentres
            scalarField& hSource = eqn.source();                //defining source
            forAll (C, i)
            {
                hSource = (Tvol - Tconst)*yourvalue;
           }

Maybe you have to do some adjustments for your case.
Hope it helps...

shaileshbg · December 20, 2017, 06:24

Hello Tobi,

Thank you very much for your reply, this is exactly what I wanted. But when I implement it I am getting an error

Code:

--> FOAM Warning : 
    From function virtual void Foam::fv::option::checkApplied() const
    in file cfdTools/general/fvOptions/fvOption.C at line 125
    Source blood_perfusion_tumor defined for field T but never used

I have defined two heat sources in my fvOptions file:
1) constant q of 13600 W/m^3 (as a scalarSemiImplicit source)
2) temperature dependent heat source with T_const at 310K

Please find my fvOptions file below:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
metabolic_heat_generation_tumor
{
    type            scalarSemiImplicitSource;
    active          true;
 
    scalarSemiImplicitSourceCoeffs
    {
        selectionMode   cellZone; // all, cellSet, cellZone, points
        cellZone         tumor;
        volumeMode      specific; // absolute;
        injectionRateSuSp
        {
            h   (13600 0);
        }
    }
}

blood_perfusion_tumor
{
    type          scalarCodedSource;
    name    sourceTime;
    active          true;
    

    scalarCodedSourceCoeffs
    {
        fields      (T);
        selectionMode   cellZone;
        cellZone        tumor;

        codeInclude
        #{

        #};

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

        codeAddSup
        #{
            const Time& time = mesh().time();
            const volScalarField& Tm = mesh_.lookupObject<volScalarField>("T"); 
            Tvol = Tm.weightedAverage(mesh_.V()).value();                //averageValue of the volScalarField
            const vectorField& C = mesh_.C();                    //List of cellcentres
            scalarField& hSource = eqn.source();                //defining source
            forAll(C, i)
            {
                hSource[i] = (Tvol - 310)* 48000;
            }
            Pout << "***codeAddSup***" << endl;
        #};

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

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

    sourceTimeCoeffs
    {
        // Dummy entry
    }
}


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

In order to monitor temperature, I have added this to my controlDict:

Code:

functions
{
    probes
    {
        type            probes;
        functionObjectLibs ("libsampling.so");
        writeControl   timeStep;
        writeInterval  1;
        probeLocations
        (
            ( 0 0.0495 0)
            ( 0 0.0595 0 )
            ( 0 0.0395 0)
            (0 0.022 0)
            (0 0.067 0 )
            (0 0.072 0)
        );
        fields
        (
            T
        );
    }
}

I have also attached my checkMesh log for your reference. Could you please let me know if I have implemented the coded source properly.

Thank you in advance
Shailesh

TobiF · December 21, 2017, 03:14

Can you please tell me, how you set up your case?

multiRegion via topoSet or via multiple stl-files?

Usually (thats how I set up my cases - doesn't have to be the best way, but is working for me) you write a single fvOptions for each region where you define the code for the single region.

Tobi

shaileshbg · December 21, 2017, 04:40

Hello Tobi,

Thank you for your reply. I have created multiregion with toposet and defined a fvOptions file separately for each region.

Problem detail: My problem has four solid hemispherical regions (muscle, gland, tissue, tumor), these are the steps that I have followed:
1) Use snappyHexMesh with searchable sphere option to define the hemispherical regions for each solid separately (as I required different mesh concentrations in each region)
2) Use merge and stitch meshes to join each meshed region into one whole meshed domain.
3) Use toposet to create different cellZones.
4) Split mesh regions based on these cellZones.
5) Define initial and boundary conditions for each cellZone use changeDictionary.

Please find a dropbox link with the case file along with a run and clean script.

https://www.dropbox.com/s/5c1hwvaoka...nFoam.zip?dl=0

Thank you very much for your time.
Shailesh

shaileshbg · January 5, 2018, 02:20

Hello Foamers,

Thanks to Tobi we have managed to add the temperature-dependent heat source as required. Below is the snippet of the fvOptions file for anyone who might face this trouble in future.

Code:

blood_perfusion_gland
{
    type          scalarCodedSource;

    active          true;
    
    name    sourceTime;

  

    scalarCodedSourceCoeffs
    {
         selectionMode   cellZone;
         cellZone        gland;
	 fields      (h);

        codeInclude
        #{

        #};

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

        codeAddSup
        #{
	     scalar Tvol = 0;
            //const Time& time = mesh().time();
            const volScalarField& Tm = mesh_.lookupObject<volScalarField>("T"); 
            Tvol = Tm.weightedAverage(mesh_.V()).value();                //averageValue of the volScalarField
            const vectorField& C = mesh_.C();                    //List of cellcentres
	    const scalarField& V = mesh_.V();
            scalarField& hSource = eqn.source();                //defining source
            forAll(C, i)
            {
                hSource[i] = (Tvol - 310)*2400*V[i] ;
            }
            Pout << "***codeAddSup***" << endl;
        #};

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

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

    sourceTimeCoeffs
    {
         $scalarCodedSourceCoeffs;
    }
}

Regards,
Shailesh

SomeUser · March 22, 2018, 03:20

Hi, I'm working on the similar case.
But I noticed that the scalarCodedCase is always working on the whole domain and not only in cellSet or cellZone.
Did your scalarCodedSource work only in the zone you specified of in the whole domain?

Thanks

TobiF · March 22, 2018, 05:40

If you choose selectionMode cellZone, like in the code stated below the codedSource should only work for this cellZone

(and the fvOptions-file is in the ./system/yourZoneName/ - directory)

Code:

    scalarCodedSourceCoeffs
    {
         selectionMode       cellZone;
         cellZone               yourZoneName;
	         fields                    (h);
 
...

And for your next post, please give some more information. of - Version, case-setup, whatyou want to do, etc.

regards Tobi

SomeUser · March 22, 2018, 10:14

Thanks Tobi, I already tried that. Although I don't understand what you mean by fvOptions should be in the ..., if I create the directory like that and put fvOptions inside than fvOptions doesn't work at all.
Btw I create cellZone from cellSet with commant setsToZones.

The Goal:
I'm using OF 3.0 the solver is reactingMultiphaseEulerFoam, I have 4 phases (one has volume fraction of 0, so is "dormant").
I'm not trying to change the temperature but the volume fraction of one phase to another.
Like the volumetric fraction of phase1 (perfectFluid in thermophysical properties) is reducing while for the phase2 (also perfectFluid) the volumetric fraction is increasing.
Case is quasi 2D (square pipe with only one cell in y), gravity is 0.
So I want to have the change in only one half of the pipe.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
options
{
    massSource1
    {
        type            scalarCodedSource;
		active			false;
        scalarCodedSourceCoeffs
        {
		selectionMode   cellZone;
		cellZone			zone1;
        fieldNames      (thermo:rho.phase2);
		redirectType	scalarCodedSource;
        name    sourceTime;
        codeInclude
        #{
		
        #};
        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};
        codeAddSup
        #{
            const Time& time1 = mesh().time(); 
            const scalarField& V1 = mesh_.V();
			const volScalarField& Am1 = mesh_.lookupObject<volScalarField>("alpha.air");
			const volVectorField& Ui1 = mesh_.lookupObject<volVectorField>("U.air");
            scalarField& deltaPhase2 = eqn.source();
            deltaPhase2 -= V1*mag(Am1)*mag(Ui1);
        #};
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};
        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
        }
    }
	massSource2
    {
        type            scalarCodedSource;
		active			true;
        scalarCodedSourceCoeffs
        {
		selectionMode   cellZone;
		cellZone			zone1;
        fieldNames      (thermo:rho.phase1);
		redirectType	scalarCodedSource;
        name    sourceTime2;
        codeInclude
        #{
		
        #};
        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};
        codeAddSup
        #{
            const Time& time2 = mesh().time();
            const scalarField& V2 = mesh_.V();
			const volScalarField& Am2 = mesh_.lookupObject<volScalarField>("alpha.air");
			const volVectorField& Ui2 = mesh_.lookupObject<volVectorField>("U.air");
            scalarField& deltaPhase1 = eqn.source();
            deltaPhase1 -= -V2*mag(Am2)*mag(Ui2);
        #};
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};
        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
        }
    }
}
// ************************************************************************* //

And I have an additional question. For quite some time, I'm searching for the equation thermo:rho.phase1 (probably the true name is thermo:rho()) in my solver and up to date I couldn't find it. Could you also help me with that?

Thanks in advance and best regards
SomeUser

BlnPhoenix · August 13, 2019, 05:15

Quote:

Originally Posted by SomeUser

Thanks Tobi, I already tried that. Although I don't understand what you mean by fvOptions should be in the ..., if I create the directory like that and put fvOptions inside than fvOptions doesn't work at all.
Btw I create cellZone from cellSet with commant setsToZones.

The Goal:
I'm using OF 3.0 the solver is reactingMultiphaseEulerFoam, I have 4 phases (one has volume fraction of 0, so is "dormant").
I'm not trying to change the temperature but the volume fraction of one phase to another.
Like the volumetric fraction of phase1 (perfectFluid in thermophysical properties) is reducing while for the phase2 (also perfectFluid) the volumetric fraction is increasing.
Case is quasi 2D (square pipe with only one cell in y), gravity is 0.
So I want to have the change in only one half of the pipe.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
options
{
    massSource1
    {
        type            scalarCodedSource;
		active			false;
        scalarCodedSourceCoeffs
        {
		selectionMode   cellZone;
		cellZone			zone1;
        fieldNames      (thermo:rho.phase2);
		redirectType	scalarCodedSource;
        name    sourceTime;
        codeInclude
        #{
		
        #};
        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};
        codeAddSup
        #{
            const Time& time1 = mesh().time(); 
            const scalarField& V1 = mesh_.V();
			const volScalarField& Am1 = mesh_.lookupObject<volScalarField>("alpha.air");
			const volVectorField& Ui1 = mesh_.lookupObject<volVectorField>("U.air");
            scalarField& deltaPhase2 = eqn.source();
            deltaPhase2 -= V1*mag(Am1)*mag(Ui1);
        #};
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};
        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
        }
    }
	massSource2
    {
        type            scalarCodedSource;
		active			true;
        scalarCodedSourceCoeffs
        {
		selectionMode   cellZone;
		cellZone			zone1;
        fieldNames      (thermo:rho.phase1);
		redirectType	scalarCodedSource;
        name    sourceTime2;
        codeInclude
        #{
		
        #};
        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};
        codeAddSup
        #{
            const Time& time2 = mesh().time();
            const scalarField& V2 = mesh_.V();
			const volScalarField& Am2 = mesh_.lookupObject<volScalarField>("alpha.air");
			const volVectorField& Ui2 = mesh_.lookupObject<volVectorField>("U.air");
            scalarField& deltaPhase1 = eqn.source();
            deltaPhase1 -= -V2*mag(Am2)*mag(Ui2);
        #};
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};
        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
        }
    }
}
// ************************************************************************* //

And I have an additional question. For quite some time, I'm searching for the equation thermo:rho.phase1 (probably the true name is thermo:rho()) in my solver and up to date I couldn't find it. Could you also help me with that?

Thanks in advance and best regards
SomeUser

Hi!

I'm trying to do the same as you. Mass transfer from one phase to another in reactingMultiphaseEulerFoam. Have you managed to do so?
Maybe we can figure this out together..

roth · November 3, 2019, 10:47

@SomeUser, you are correct. If you are not careful it will apply the source to the entire domain. There is discussion and work-around here:
https://bugs.openfoam.org/view.php?id=2361

I'm copying the example given by @MattijsJ:

Code:

    energySource
    {
        type       scalarCodedSource;
        active     true;
        name       sourceTime;  

        scalarCodedSourceCoeffs
        {
            selectionMode   points;
            points ( 
                   ( 2.0  2.0 2.0)
                                  ); 
//        selectionMode   all;

//        selectionMode   cellZone; 
//        cellZone        heater; 

            fields          (h);

            codeInclude
            #{
            #};

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

            codeAddSup
            #{
                const Time& time = mesh().time();
                const scalarField& V = mesh_.V();
                scalarField& heSource = eqn.source();

                const labelList& cellIDs = cells();

                forAll(cellIDs, i)
                {
                    label cellI = cellIDs[i];
                    heSource[cellI] -= 0.0001*sqr(time.value())*V[cellI];
                }
            #};

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

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

        sourceTimeCoeffs
        {
            $scalarCodedSourceCoeffs;
        }
    }

julieng · December 30, 2020, 17:06

Hello,

Is it possible to add an heat source on a boundary (using facezone or patch name) ?
I want to add on the interface solid-fluid the heat flux:

qr=emi x cst_boltz x (T(boundary temp)^4 - Tcst^4)

Best regards

Julien

choist31 · March 18, 2021, 14:33

Dear Mr. Julieng

I am woking on the similar job. I wonder if you have any resolution on your approch. It would be great if you can share your experience with me.
Thanks for any help in advance.

Best wishes,

Sung

Sully · January 2, 2023, 13:57

Dear all,

I would like to know if it is possible to set heat source with temperature dependency (in fvOptions or not) but by avoiding the use of scalarCodedSourceCoeffs ? Coded part does not work in my case.

Thank you,

jose.ceballos · January 26, 2024, 19:48

Can i use fvOptions in OpenFOAM 4.x? I want to add source terms to the energy equation of a PEM fuel cell.

December 19, 2017, 10:33	Adding temperature dependent heat source using fvOptions in chtMultiRegionSimpleFoam	#1
shaileshbg Member Shailesh BG Join Date: Aug 2017 Location: Bangalore Posts: 39 Rep Power: 9	Hello Foamers, I want to add a temperature dependent heat source into the chtMultiRegionSimpleFoam solver using fvOptions. I was looking at semiImplicitSource option which works on source term linearization using SuSp injection rate treatment where: S(x) = Su + Sp(x) S(x): net source for field (x) Su: explicit source contribution Sp: linearized implicit contribution I understand that the solver works by evaluating enthalpy (so my field is H) and that the unit of enthalpy is same as that of internal energy (J/Kg), so if I want to add a constant heat source I need not convert anything I can put it in directly under Su (taking absolute and specific into consideration), but I do not know how to add a temperature dependent source term of this form: source = constant * (T - T_const) Problem definition: Heat conduction through four different solids each having its own temperature-dependent source. Could anyone please help me to understand how to implement this. Many thanks in advance, Shailesh

December 19, 2017, 14:50		#2
TobiF New Member Join Date: Apr 2014 Posts: 24 Rep Power: 12	Hey, use scalarCodedSource to calculate the temperature-depending heat-source. If you need the mean-temp of the region you should define the selectionMode type cellZone and in the code something like : Code: const volScalarField& Tm = mesh_.lookupObject<volScalarField>("T"); Tvol = Tm.weightedAverage(mesh_.V()).value(); //averageValue of the volScalarField to calculate the mean tmep of the cellZone. After that you can use the Tvol - value for your calculation. Code: const vectorField& C = mesh_.C(); //List of cellcentres scalarField& hSource = eqn.source(); //defining source forAll (C, i) { hSource = (Tvol - Tconst)*yourvalue; } Maybe you have to do some adjustments for your case. Hope it helps...

December 21, 2017, 03:14		#4
TobiF New Member Join Date: Apr 2014 Posts: 24 Rep Power: 12	Can you please tell me, how you set up your case? multiRegion via topoSet or via multiple stl-files? Usually (thats how I set up my cases - doesn't have to be the best way, but is working for me) you write a single fvOptions for each region where you define the code for the single region. Tobi Tobi likes this.

December 21, 2017, 04:40	multiRegion via topoSet	#5
shaileshbg Member Shailesh BG Join Date: Aug 2017 Location: Bangalore Posts: 39 Rep Power: 9	Hello Tobi, Thank you for your reply. I have created multiregion with toposet and defined a fvOptions file separately for each region. Problem detail: My problem has four solid hemispherical regions (muscle, gland, tissue, tumor), these are the steps that I have followed: 1) Use snappyHexMesh with searchable sphere option to define the hemispherical regions for each solid separately (as I required different mesh concentrations in each region) 2) Use merge and stitch meshes to join each meshed region into one whole meshed domain. 3) Use toposet to create different cellZones. 4) Split mesh regions based on these cellZones. 5) Define initial and boundary conditions for each cellZone use changeDictionary. Please find a dropbox link with the case file along with a run and clean script. https://www.dropbox.com/s/5c1hwvaoka...nFoam.zip?dl=0 Thank you very much for your time. Shailesh

March 22, 2018, 05:40		#8
TobiF New Member Join Date: Apr 2014 Posts: 24 Rep Power: 12	If you choose selectionMode cellZone, like in the code stated below the codedSource should only work for this cellZone (and the fvOptions-file is in the ./system/yourZoneName/ - directory) Code: scalarCodedSourceCoeffs { selectionMode cellZone; cellZone yourZoneName; fields (h); ... And for your next post, please give some more information. of - Version, case-setup, whatyou want to do, etc. regards Tobi

March 22, 2018, 03:20		#7
SomeUser New Member Join Date: Feb 2018 Posts: 3 Rep Power: 8	Hi, I'm working on the similar case. But I noticed that the scalarCodedCase is always working on the whole domain and not only in cellSet or cellZone. Did your scalarCodedSource work only in the zone you specified of in the whole domain? Thanks

December 30, 2020, 17:06		#12
julieng Senior Member julien Join Date: Dec 2018 Posts: 107 Rep Power: 8	Hello, Is it possible to add an heat source on a boundary (using facezone or patch name) ? I want to add on the interface solid-fluid the heat flux: qr=emi x cst_boltz x (T(boundary temp)^4 - Tcst^4) Best regards Julien

March 18, 2021, 14:33	Questions on heat flux model	#13
choist31 New Member sw choi Join Date: Aug 2020 Location: MD in USA Posts: 20 Rep Power: 6	Dear Mr. Julieng I am woking on the similar job. I wonder if you have any resolution on your approch. It would be great if you can share your experience with me. Thanks for any help in advance. Best wishes, Sung

January 2, 2023, 13:57	Heat source with temperature dependent without scalarCodedSourceCoeffs	#14
Sully New Member État Join Date: Jan 2023 Posts: 1 Rep Power: 0	Dear all, I would like to know if it is possible to set heat source with temperature dependency (in fvOptions or not) but by avoiding the use of scalarCodedSourceCoeffs ? Coded part does not work in my case. Thank you,

January 26, 2024, 19:48	fvOptions in OF4.x	#15
jose.ceballos New Member Jose Join Date: Dec 2020 Posts: 18 Rep Power: 6	Can i use fvOptions in OpenFOAM 4.x? I want to add source terms to the energy equation of a PEM fuel cell.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
OpenFOAM without MPI	kokizzu	OpenFOAM Installation	4	May 26, 2014 10:17
Transient spacial dependent heat source	Kumudu	OpenFOAM Running, Solving & CFD	0	October 28, 2013 12:02
[swak4Foam] build problem swak4Foam OF 2.2.0	mcathela	OpenFOAM Community Contributions	14	April 23, 2013 14:59
OpenFOAM on MinGW crosscompiler hosted on Linux	allenzhao	OpenFOAM Installation	127	January 30, 2009 20:08
DxFoam reader update	hjasak	OpenFOAM Post-Processing	69	April 24, 2008 02:24