Temperature Dependent Thermal Properties in LaplacianFoam

Akram_ · January 4, 2022, 16:13

I'm solving the heat equation using the finite volume method in OpenFOAM and the way it is set up right now is to average the thermal properties of the temperature field. What I want to achieve is to implement the thermal properties that correspond to the temperature of each cell/mesh. Is there a way to do that through the solver itself or in the transportProperties file.

Attached is the code I'm using. I just replaced the source term with the volumetric power that is calculated by another code.

Code:

Info<< "\nCalculating temperature distribution\n" << endl;

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            solve
            (
                fvm::ddt(T) - fvm::laplacian(DT, T) == volPower/(RHO * CP)
            );

        }

        #include "write.H"

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;

This is the script that reads the thermal properties from the transportProperties file.

Code:

Info<< "Reading field T\n" << endl;

    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );


    volScalarField volPower
    (
        IOobject
        (
            "volPower",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );



    Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );


    Info<< "Reading diffusivity DT\n" << endl;

    dimensionedScalar DT
    (
        transportProperties.lookup("DT")
    );

    Info<< "Reading density RHO\n" << endl;

    dimensionedScalar RHO
    (
        transportProperties.lookup("RHO")
    );

    Info<< "Reading thermal conductivity CP\n" << endl;

    dimensionedScalar CP
    (
        transportProperties.lookup("CP")

This is how the transportProperties file is structured

Code:

DT              DT  [0 2 -1 0 0 0 0] 1.7620248e-5; // Diffusion Coefficent 

RHO             RHO [1 -3 0 0 0 0 0] 1.52082743e3;  // Density 

CP              CP  [0 2 -2 -1 0 0 0] 707.36193;  // Heat Capacity

I want DT and CP to be dependant on the temperature of each cell

LuckyTran · January 4, 2022, 16:36

create a volScalarField for CP and DT. It doesn't have to be an IOobject, but you could copy paste the object for T and rename one each for CP and DT (and set AUTO_READ or NO_READ accordingly). Then add two lines in your function somewhere before the solve line where you calculate (i.e. update) CP and DT based on whatever function you have for each iteration.

Akram_ · January 4, 2022, 19:41

So I did that and removed the definitions of DT and CP from the transportproprties file. But I got this error

Code:

FOAM FATAL IO ERROR: 
keyword DT is undefined in dictionary "/home/..../transportProperties"

When I add the DT and CP definitions the code runs normally. I think it ignores the functions I wrote for DT and CP and uses the values in the transportproprties file.

I suspect there is something wrong with the way I defined these parameters.

This is how I edited createFields.H.

Code:

volScalarField CP
    (
        IOobject
        (
            "CP",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ ,
            IOobject::AUTO_WRITE
        ),
        mesh
	);
	
	
	volScalarField DT
	(
        IOobject
        (
            "DT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
	);
	
	volScalarField K
    (
        IOobject
        (
            "K",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
	);
	
	
    Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    Info<< "Reading density RHO\n" << endl;

    dimensionedScalar RHO
    (
        transportProperties.lookup("RHO")
    );

I kept RHO in the transportProperties file

This is how I defined the functions of DT and CP

Code:

while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {	
			CP= 528.99369+4.03815*T
			K= 28.20519-0.01988*T+0.00000882446*T*T
			DT=K/(RHO*CP)
            solve
            (
                fvm::ddt(T) - fvm::laplacian(DT, T) == volPower/(RHO * CP)
            );

Akram_ · January 4, 2022, 23:27

I posted a reply but apparently, it wasn't received. I did all of that but The code ignores it and reads the constant values in the transportproperties file anyway. Is there a way to read the DT and CP values from a file like it does with the volpower and T?

LuckyTran · January 5, 2022, 00:38

The code is not ignoring anything, it is doing exactly what you tell it to do.

If you want it to read from files then use IOobject::MUST_READ. If you don't want it to read from the transportProperties file, then stop telling it to. Delete the lines that say transportProperties.lookup

If you read the fields from one file and overwrite them immediately after, well that is exactly what happens.

It might help your learning if you use separate names for everything instead of DT CP. Use silly names like banana, apple, pineapple, etc. You need to abstract references to fields and what you think they are in your head.

Akram_ · January 5, 2022, 01:19

Do I need to keep the transportproprties file as it is? I get an error if I remove any value from this file.

This is creatFields.H after removing all reference to the transportproprties file.

Code:

Info<< "Reading field T\n" << endl;

    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );


    volScalarField volPower
    (
        IOobject
        (
            "volPower",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );


	volScalarField CP
    (
        IOobject
        (
            "CP",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ ,
            IOobject::NO_WRITE
        ),
        mesh
	);
	
	
	volScalarField DT
	(
        IOobject
        (
            "DT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh
	);
	
	volScalarField K
    (
        IOobject
        (
            "K",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh
	);
	
	volScalarField RHO
    (
        IOobject
        (
            "RHO",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh
	);

And this is the solver part

Code:

while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {	
			CP= 528.99369+4.03815*T
			K= 28.20519-0.01988*T+0.00000882446*T*T
			RHO =1.52082743e3 
			DT=K/(RHO*CP) 
            solve
            (
                fvm::ddt(T) - fvm::laplacian(DT, T) == volPower/(RHO * CP)
            );

        }

Even though I didn't reference transportproprties anywhere it still reads from it and use its values. I changed these values and tested some case and the results are different. It seems that it doesn't depend on the equations of DT and CP at all.

LuckyTran · January 5, 2022, 02:09

And this is exactly why I recommended you to use silly names. Instead of fvm::laplacian(DT, T) if you had done instead fvm::laplacian(cowabunga, T) then you would be sure that DT is not invoked by the solver.

I don't know what else to tell you. Either you forgot to compile the new code or you have omitted other parts of the code where you do read from the transportProperties dict.

Also, c statements end in semicolons ; so there's no way your solver actually compiled the way it is, unless this is not your code or you failed to copy and paste it.

Akram_ · January 5, 2022, 17:58

Yes, there was an issue with the compiling script. I fixed it and then compiled it correctly. I also changed the names of the variables and it doesn't read from the transportproperties file this time.

I have two options; either the code reads the thermal properties the same way it reads T and volpower or it uses defined equations. I tried the first option and constructed files for DT, CP and RHO as DT_new , C_P and RHO_new. The code reads these files with no problem but it crashes before excusting the fvm solver.

This is the error message in the log.laplacianVolPowerFoam file.

Code:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  ? at ??:?
#5  ? at ??:?
#6  __libc_start_main in "/lib64/libc.so.6"
#7  ? at ??:?

I also tried the other option but I get a dimensions error. I defined the thermal properties as follows:

Code:

	volScalarField RHO_new
 (
	 IOobject
	 (
		 "RHO_new",
		 runTime.timeName(),
		 mesh,
		 IOobject::NO_READ,
		 IOobject::NO_WRITE
	 ),
	 mesh,
	 dimensionSet(1, -3, 0, 0, 0, 0, 0)
 );

 forAll( T , i )
 {
	 RHO_new =1520.757982*Foam::sqrt(T[i]); 
 }


	volScalarField CP_new
 (
	 IOobject
	 (
		 "CP_new",
		 runTime.timeName(),
		 mesh,
		 IOobject::NO_READ,
		 IOobject::NO_WRITE
	 ),
	 mesh,
	 dimensionSet(0, 2, -2, -1, 0, 0, 0)
 );

 forAll( T , i )
 {
	 CP_new= 528.99369*T[i]; 
 }
 
 
 volScalarField DT_new
 (
	 IOobject
	 (
		 "DT_new",
		 runTime.timeName(),
		 mesh,
		 IOobject::NO_READ,
		 IOobject::NO_WRITE
	 ),
	 mesh,
	 dimensionSet(0, 2, -1, 0, 0, 0, 0)
 );

 forAll( T , i )
 {
	 DT_new= 1.76054e-5*Foam::sqrt(T[i]);
 }

Note that I removed the equations before the solve line in the solver.c

Code:

while (simple.correctNonOrthogonal())
        {	
		
            solve
            (
                fvm::ddt(T) - fvm::laplacian(DT_new, T) == volPower/(RHO_new * CP_new)
            );

        }

        #include "write.H"

This is the error that I got

Code:

--> FOAM FATAL ERROR: 
Different dimensions for =
     dimensions : [1 -3 0 0 0 0 0] = [0 0 0 0 0 0 0]


    From function bool Foam::dimensionSet::operator=(const Foam::dimensionSet&) const
    in file dimensionSet/dimensionSet.C at line 171.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::dimensionSet::operator=(Foam::dimensionSet const&) const at ??:?
#3  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::operator=(Foam::dimensioned<double> const&) at ??:?
#4  ? at ??:?
#5  __libc_start_main in "/lib64/libc.so.6"
#6  ? at ??:?

I set their dimensions according to how they were set in the transportproperties file. What do you suggest?

LuckyTran · January 6, 2022, 18:07

Putting the formulas in createFields means that they are calculated only once at the start of the run when the fields are initially created. You still need to update them in the code using formulas at runtime if you want temperature dependent properties to be calculated on the fly.

A floating point error is a floating point error. You just need to work through why you would have that. It could a division by zero error or something else that's silly.

Either you do everything in raw numbers or you need to keep track of all the units everywhere all the time. Every number you input needs to have units with it.

Your coefficients in the formulas need to be dimensionedScalars, i.e. 528.99369 needs units. And so does everything else. Btw, what does your T file look like? Your error code shows you are trying to assign a field with dimensions all zeros to a field with dimensions of k/gm^3 (likely your RHO_new). It looks like your T field might also dimensionless. Make sure your T file has the right dimensions. You have set the dimensions on the fields used to store the numbers you care about, but all arithmetic in OpenFOAM is dimensioned. When you multiple two numbers together, both the inputs need to have the proper dimensions as well.

Akram_ · January 7, 2022, 01:59

Thank you so much for your help. I managed to finally make it work.

I defined the coefficients of the following equations

Code:

while (simple.correctNonOrthogonal())
        {	
			CP_new= Cp_1+Cp_2*T;
			k= k_1-k_2*T+k_3*T*T;
			DT=k/(RHO*CP_new);
			
			
            solve
            (
                fvm::ddt(T) - fvm::laplacian(DT, T) == volPower/(RHO * CP_new)
            );

        }

I defined these coefficients as dimensionedScalar in the createFields.H as follows:

Code:


volScalarField k
    (
        IOobject
        (
            "k",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
			IOobject::NO_WRITE
        ),
        mesh,
		dimensionSet(1,1,-3,-1,0,0,0)
    );

	

	volScalarField CP_new
    (
        IOobject
        (
            "CP_new",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
			IOobject::NO_WRITE
        ),
        mesh,
		dimensionSet(0, 2, -2, -1, 0, 0, 0)
    );
  
    
 
	volScalarField DT
    (
        IOobject
        (
            "DT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
			IOobject::NO_WRITE
        ),
        mesh,
		dimensionSet(0, 2, -1, 0, 0, 0, 0)
    );


Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );


   

    dimensionedScalar Cp_1 
    (
        transportProperties.lookup("Cp_1")
    );
	
	

	dimensionedScalar Cp_2 
    (
        transportProperties.lookup("Cp_2")
    );
	
	
	dimensionedScalar k_1 
    (
        transportProperties.lookup("k_1")
    );
	
	
	dimensionedScalar k_2 
    (
        transportProperties.lookup("k_2")
    );

	dimensionedScalar k_3 
    (
        transportProperties.lookup("k_3")
    );
	
	dimensionedScalar RHO 
    (
        transportProperties.lookup("RHO")
    );

Finally, I edited the transportproperties file to include these coefficients with appropriate units

Code:

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

Cp_1             Cp_1  [0 2 -2 -1 0 0 0] 528.99369;

Cp_2             Cp_2  [0 2 -2 -2 0 0 0] 4.03815;  

k_1              k_1   [1 1 -3 -1 0 0 0] 28.20519;

k_2              k_2   [1 1 -3 -2 0 0 0] 0.01988;  

k_3              k_3   [1 1 -3 -3 0 0 0] 0.00000882446;

RHO              RHO   [1 -3 0 0 0 0 0] 1.52082743e3;

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Temperature dependent thermal conductivity and specific heat Property of mixture of t	rohan@123	FLUENT	0	April 17, 2020 11:34
Temperature dependent thermal conductivity and specific heat Property of DH36 steel	rohan@123	FLUENT	1	April 10, 2020 16:33
Domain Reference Pressure and mass flow inlet boundary	AdidaKK	CFX	75	August 20, 2018 06:37
Temperature dependent material properties	stuart230588	CFX	1	November 21, 2013 17:36
how to incorporate temperature dependent thermophysical properties in fluent.	CANDY	Fluent UDF and Scheme Programming	4	October 22, 2012 05:19

January 4, 2022, 16:36		#2
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,752 Rep Power: 66	create a volScalarField for CP and DT. It doesn't have to be an IOobject, but you could copy paste the object for T and rename one each for CP and DT (and set AUTO_READ or NO_READ accordingly). Then add two lines in your function somewhere before the solve line where you calculate (i.e. update) CP and DT based on whatever function you have for each iteration. Akram_ likes this.

January 4, 2022, 23:27		#4
Akram_ New Member Join Date: Jan 2022 Posts: 6 Rep Power: 4	I posted a reply but apparently, it wasn't received. I did all of that but The code ignores it and reads the constant values in the transportproperties file anyway. Is there a way to read the DT and CP values from a file like it does with the volpower and T?

January 5, 2022, 00:38		#5
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,752 Rep Power: 66	The code is not ignoring anything, it is doing exactly what you tell it to do. If you want it to read from files then use IOobject::MUST_READ. If you don't want it to read from the transportProperties file, then stop telling it to. Delete the lines that say transportProperties.lookup If you read the fields from one file and overwrite them immediately after, well that is exactly what happens. It might help your learning if you use separate names for everything instead of DT CP. Use silly names like banana, apple, pineapple, etc. You need to abstract references to fields and what you think they are in your head.

January 5, 2022, 02:09		#7
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,752 Rep Power: 66	And this is exactly why I recommended you to use silly names. Instead of fvm::laplacian(DT, T) if you had done instead fvm::laplacian(cowabunga, T) then you would be sure that DT is not invoked by the solver. I don't know what else to tell you. Either you forgot to compile the new code or you have omitted other parts of the code where you do read from the transportProperties dict. Also, c statements end in semicolons ; so there's no way your solver actually compiled the way it is, unless this is not your code or you failed to copy and paste it.

January 6, 2022, 18:07		#9
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,752 Rep Power: 66	Putting the formulas in createFields means that they are calculated only once at the start of the run when the fields are initially created. You still need to update them in the code using formulas at runtime if you want temperature dependent properties to be calculated on the fly. A floating point error is a floating point error. You just need to work through why you would have that. It could a division by zero error or something else that's silly. Either you do everything in raw numbers or you need to keep track of all the units everywhere all the time. Every number you input needs to have units with it. Your coefficients in the formulas need to be dimensionedScalars, i.e. 528.99369 needs units. And so does everything else. Btw, what does your T file look like? Your error code shows you are trying to assign a field with dimensions all zeros to a field with dimensions of k/gm^3 (likely your RHO_new). It looks like your T field might also dimensionless. Make sure your T file has the right dimensions. You have set the dimensions on the fields used to store the numbers you care about, but all arithmetic in OpenFOAM is dimensioned. When you multiple two numbers together, both the inputs need to have the proper dimensions as well.