# Problem in modified pisoFoam with temperature equation + thermophysical model

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

 August 29, 2012, 09:35 Problem in modified pisoFoam with temperature equation + thermophysical model #1 Member   Join Date: Feb 2012 Posts: 35 Rep Power: 12 Hi guys, I recently modified "pisoFoam" to consider temperature varying in the domain. After some problems, I succeed with compiling the solver. For now, I just entered the temperature equation and let it vary along the time running. It works fine until now. At this point, I'd like to get vary some thermophysical property by temperature. I read dozens of threads in this forum but I think I still don't get very well how to do it. First of all, I would like to vary the density by temperature; in a second time I'd like to use this density to vary nu() and then to affect the momentum equation. At this moment, I just want to see density vary by temperature: I'm not even able to do this right now . I chose "basicRhoThermo" as my thermo model and I'm trying to use "thermo.rho()" which shall be calculated by means of "icoPolynomial"; I see that T is correctly calculated as nonuniform field along the mesh, for every time step; on the contrary, rho is calculated with icoPolynomial only on the boundaries as uniform field and then, for every time step, remains unchanged: I don't really understand how to overcome this....where am I mastaken ???. I post here the fundamental part of my code (in red the new code referring to plain pisoFoam): thermophysicalProperties dictionary: Code: thermoType hRhoThermo>; pRef 101325; mixture { equationOfState { rhoCoeffs<8> ( 1000 -1.1 0 0 0 0 0 0); <------this is just a dummy test } specie { nMoles 1; molWeight 28.9; } thermodynamics { Hf 0; Sf 0; CpCoeffs<8> ( 1000 0 0 0 0 0 0 0); } transport { muCoeffs<8> (0.3 -0.0008 0.0000007 -0.0000000001 0 0 0 0); kappaCoeffs<8> ( 1 1e-5 0 0 0 0 0 0); } } transportProperties dictionary: Code: transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1e-05; // Laminar Prandtl number Pr Pr [0 0 0 0 0 0 0] 0.9; // Turbulent Prandtl number Prt Prt [0 0 0 0 0 0 0] 0.7; pisoFoamT.C: Code: #include "fvCFD.H" #include "singlePhaseTransportModel.H" #include "turbulenceModel.H" #include "basicRhoThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "createFields.H" #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; #include "readPISOControls.H" #include "CourantNo.H" // Pressure-velocity PISO corrector { // Momentum predictor fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence->divDevReff(U) ); UEqn.relax(); if (momentumPredictor) { solve(UEqn == -fvc::grad(p)); } // --- PISO loop for (int corr=0; corrcorrect(); #include "TEqn.H" runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } createFields.H: Code:  Info<< "Reading field p\n" << endl; volScalarField p ( IOobject ( "p", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field T\n" << endl; <----added new T field volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field kappat\n" << endl; <----added new kappat field volScalarField kappat for turbulent T field calculation ( IOobject ( "kappat", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); # include "createPhi.H" label pRefCell = 0; scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); singlePhaseTransportModel laminarTransport(U, phi); dimensionedScalar Pr(laminarTransport.lookup("Pr")); dimensionedScalar Prt(laminarTransport.lookup("Prt")); autoPtr turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Info<< "Reading thermophysical properties\n" << endl; autoPtr pThermo ( basicRhoThermo::New(mesh) ); basicRhoThermo& thermo = pThermo(); volScalarField rho ( IOobject ( "rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), thermo.rho() ); rho.oldTime().write();  TEqn.H: Code:  kappat = turbulence->nut()/Prt; kappat.correctBoundaryConditions(); volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat); fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(); thermo.correct(); rho=thermo.rho(); Any hint would be appreciated!! Thank you. Matteo y_jiang and Thamali like this.

 August 30, 2012, 11:44 #2 Member   Join Date: Feb 2012 Posts: 35 Rep Power: 12 Problem not yet resolved, but I made a test that let me think the problem is in the thermo model. I made the test on T to understand why rho is not calculated properly. I define T as before: Code:  Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); In this way, after calculation of TEqn, T is calculated properly along the whole mesh as nonuniform field, for every time step. But if I define T in this different way, I would expect the same result: Code:  Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), thermo.T() ); ....it is not! Now, trying to exploit the thermo model for T, T behaves like rho, so generates just copy of the 0/T file for the new time steps folders, ignoring the presence of the TEqn.....How to explain this?? Please, I need someone's help!!

 September 5, 2012, 12:14 Issues resolved! #3 Member   Join Date: Feb 2012 Posts: 35 Rep Power: 12 Hi guys, finally I fixed and resolved all problems in my new pisoFoamT solver (for now....). The main error was that I didn't really understand very well how to update and calculate a thermo property, exploiting the thermo type specified; then I realized that I have specified "hRhoThermo" as a model to calculate density which is dependent on enthalpy calculation: so, what I have to do was to calculate somewhere enthalpy by means of an equation definition. Below you can find the already posted files of my solver/case (the modified ones). Hope can help someone to understand better how to manage this issues. Enjoy! thermophysicalProperties dictionary: Code: thermoType hRhoThermo>; mixture { specie { Liquid1; nMoles 0.5; molWeight 28.9; Liquid2; nMoles 0.5; molWeight 50; } equationOfState { Liquid1; rhoCoeffs<8> ( 1110 -0.447 0 0 0 0 0 0); Liquid2; rhoCoeffs<8> ( 500 -0.1 0 0 0 0 0 0); } thermodynamics { Liquid1; Hf 0; Sf 0; CpCoeffs<8> ( 1000 0.05 0 0 0 0 0 0); Liquid2; Hf 0; Sf 0; CpCoeffs<8> ( 800 0.05 0 0 0 0 0 0); } transport { Liquid1; muCoeffs<8> (0.3 -0.0008 0.0000007 -0.0000000001 0 0 0 0); kappaCoeffs<8> ( 1 1e-5 0 0 0 0 0 0); Liquid2; muCoeffs<8> (1 -0.0008 0.0000007 -0.0000000001 0 0 0 0); kappaCoeffs<8> ( 1 1e-5 0 0 0 0 0 0); } } Here I added properties of a second liquid. Mean properties values between Liquid1 and Liquid2 properties are calculated during running. createFields.H: Code:  Info<< "Reading field p\n" << endl; volScalarField p ( IOobject ( "p", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field kappat\n" << endl; volScalarField kappat ( IOobject ( "kappat", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); # include "createPhi.H" label pRefCell = 0; scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); singlePhaseTransportModel laminarTransport(U, phi); dimensionedScalar Pr(laminarTransport.lookup("Pr")); dimensionedScalar Prt(laminarTransport.lookup("Prt")); autoPtr turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Info<< "Reading thermophysical properties\n" << endl; autoPtr pThermo ( basicRhoThermo::New(mesh) ); basicRhoThermo& thermo = pThermo(); volScalarField& h=thermo.h(); volScalarField rho ( IOobject ( "rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), thermo.rho() ); rho.write(); // scrive rho al tempo 0 volScalarField nu ( IOobject ( "nu", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), thermo.mu()/rho ); TEqn.H: Code: { kappat = turbulence->nut()/Prt; kappat.correctBoundaryConditions(); // volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat); volScalarField kappaEff("kappaEff", nu/Pr + kappat); fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(); h=thermo.Cp()*T; //<--------this is enthalpy calculation! thermo.correct(); rho=thermo.rho(); nu=thermo.mu()/rho; } Mojtaba.a, songwukong, DuarteMagalhaes and 1 others like this.

 November 20, 2012, 13:22 #5 Member   Join Date: Feb 2012 Posts: 35 Rep Power: 12 Hi Markus, here below I post the files you asked for, but I think it would be even better if you'd expose which problems/errors you encountered along your running case. fvSchemes: Code: ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes { default none; div(phi,U) Gauss limitedLinearV 1; div(phi,T) Gauss limitedLinear 1; div(phi,k) Gauss limitedLinear 1; div(phi,epsilon) Gauss limitedLinear 1; div(phi,R) Gauss limitedLinear 1; div(R) Gauss linear; div(phi,nuTilda) Gauss limitedLinear 1; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear corrected; laplacian(kappaEff,T) Gauss linear corrected; laplacian((1|A(U)),p) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } fvSolution: Code: solvers { p { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0.1; } pFinal { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; } U { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } T { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } k { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } epsilon { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } R { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } nuTilda { solver PBiCG; preconditioner DILU; tolerance 1e-05; relTol 0; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } controlDict: Code: application pisoFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 10; deltaT 0.005; writeControl timeStep; writeInterval 100; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; See ya! Matteo songwukong likes this.

 November 20, 2012, 13:41 #6 New Member   Markus Trompa Join Date: Nov 2012 Location: Regensburg, Germany Posts: 13 Rep Power: 11 Hi Matteo, Thank you very much for your quick reply. The problem which encoutered concerned the Courant number, I just had to switch adjustableTimeStep to on in the controlDict and define a maxCo. Now the problem is fixed. See ya Markus

 October 10, 2013, 06:36 #7 New Member   Konrad Join Date: Sep 2013 Posts: 4 Rep Power: 10 Hi Matteo, As I pmed you I will describe my case here, I really need your help beacause you have written long time before that u managed to modify pisoFoam so that it compute temperature pool. Right now my case is a turbulent water flow through the tube with concentrical cylinder inside it with lets say radius = about 1/4 diameter of tube. ( the cylinder is in the middle of tube and tooks about 1/3 lengh of all tube in mesh geometry model ) I have succesfully solved this flow with LES model included with pisoFoam and results are reasonable and satisfactory. What I want to do now is adding to this model a temperature so that the cylinder is a heat source ( with constant temp. on the wall ) and recompile flow. Regretfully just having the knowledge how to add T to icoFoam is not enough for me to work this case out, code files are different and although solver starts running actual solving doesn't take place ( without any foam fatal error ) and for example T field isn't even runing. I consider density as constant. The temperature difference will be rather small. To be honest, the simplest way would be if you, Matteo, could remind how you managed to add T to pisoFoam succesfully in the simplest way there is for the first time. Best Regards, Konrad Last edited by Byxon; October 10, 2013 at 08:50.

October 20, 2013, 13:51
#8
Member

Join Date: Feb 2012
Posts: 35
Rep Power: 12
Quote:
 Originally Posted by Byxon Hi Matteo, As I pmed you I will describe my case here, I really need your help beacause you have written long time before that u managed to modify pisoFoam so that it compute temperature pool. Right now my case is a turbulent water flow through the tube with concentrical cylinder inside it with lets say radius = about 1/4 diameter of tube. ( the cylinder is in the middle of tube and tooks about 1/3 lengh of all tube in mesh geometry model ) I have succesfully solved this flow with LES model included with pisoFoam and results are reasonable and satisfactory. What I want to do now is adding to this model a temperature so that the cylinder is a heat source ( with constant temp. on the wall ) and recompile flow. Regretfully just having the knowledge how to add T to icoFoam is not enough for me to work this case out, code files are different and although solver starts running actual solving doesn't take place ( without any foam fatal error ) and for example T field isn't even runing. I consider density as constant. The temperature difference will be rather small. To be honest, the simplest way would be if you, Matteo, could remind how you managed to add T to pisoFoam succesfully in the simplest way there is for the first time. Best Regards, Konrad
I read your description but I think there is something missing; first of all, what kind of geometry is yours, is it a cavity or an external flux and so, what kind of fluid simulation is the purpose of your task? You should be much more precise...
Anyway, you don't have just to add and create the T field in createFields.H and the trick is done; most important is to add the new T or Enthalpy equation T field, a thermophisicalProperties dictionary and an equation to calculate enthalpy: did you do any of this steps which I listed in the previous posts?

Matteo

 October 21, 2013, 04:06 #9 New Member   Konrad Join Date: Sep 2013 Posts: 4 Rep Power: 10 Here is a pic of my case. It is a cross-section of a cylindrical tube with a incompressible flow of water with U on inlet 1,5 m/s. The "white" solid cylinder inside is something what I want to be a heat-source with uniform temperature on the surface. I have tried to follow the tutorial from openFoamWiki which adds temperature equation to the icoFoam solver and use it in cavity tutorial, so that it fits my case but I have failed. The solver "starts" but T is ignored. There is no error though. I'm playing with solvers in OpenFOAM for the first time, that's why I need a support of some1 who succeed before in this field. I have all your files considering the case which vary for ex. density by temperature which is not my goal, so I can't just follow your path. Best Regards Konrad

October 21, 2013, 06:56
#10
New Member

Join Date: Sep 2013
Posts: 4
Rep Power: 10
I'm getting runs like this:

Create time

Create mesh for time = 0

Selecting incompressible transport model Newtonian
Selecting turbulence model type LESModel
Selecting LES turbulence model oneEqEddy
Selecting LES delta type cubeRootVol
bounding k, min: 0 max: 2e-05 average: 0
oneEqEddyCoeffs
{
ce 1.048;
ck 0.094;
}

Starting time loop

End

I have attached some of my files.

FvSchemes:

Code:
FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
location    "system";
object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default         Euler;
}

{
default         Gauss linear;
}

divSchemes
{
default         none;
div(phi,U)      Gauss linear;
div(phi,T) Gauss upwind;
}

laplacianSchemes
{
default         none;
laplacian(nu,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;

laplacian(DT,T) Gauss linear corrected;
}

interpolationSchemes
{
default         linear;
interpolate(HbyA) linear;
}

{
default         corrected;
}

fluxRequired
{
default         no;
p               ;
}
fvSolutions:
Code:
FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
location    "system";
object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
p
{
solver          PCG;
preconditioner  DIC;
tolerance       1e-06;
relTol          0.05;
}
pFinal
{
solver         PCG;
preconditioner    DIC;
tolerance         1e-06;
relTol            0;
}
T
{
solver           smoothSolver;
smoother         GaussSeidel;
tolerance        1e-6;
relTol           0.01;
nSweeps          3;
maxIter    100;
minIter 10;
}

U
{
solver          PBiCG;
preconditioner  DILU;
tolerance       1e-05;
relTol          0;
}

k
{
solver            PBiCG;
preconditioner     DILU;
tolerance        1e-05;
relTol            0;
}

B
{
solver            PBiCG;
preconditioner    DILU;
tolerance        1e-05;
relTol            0;
}
nuTilda
{
solver         PBiCG;
preconditioner     DILU;
tolerance        1e-05;
relTol            0;

}
}

PISO
{
nCorrectors     2;
nNonOrthogonalCorrectors 0;
pRefCell    0;
pRefValue    0;

}
Attached Files
 createFields.H (1.3 KB, 26 views) pisoFoamT.C (4.7 KB, 39 views)

 October 23, 2013, 08:46 #11 Member   Join Date: Feb 2012 Posts: 35 Rep Power: 12 Hi Konrad, just a brief reply to your questions for now. To me seems very strange you cannot verify the T creation and calculation during running of your case. Everything seems to me correct, the T creation in createFields.H and the T equation definition in the pisoFoamT.C file. Then I have a doubt on what you did in preparation of the simulation of your case: did you correctly compile your solver? I'm thinking you forgot to compile the solver after you added new lines regarding T, so that's why new code lines are skipped on the new runnings. Maybe you already know, but anyway, to compile the solver you have to go to solver directory and type on the terminal "wclean" and after that "wmake"; try again now to run the solver. If this is not your ploblem then let me know again. Sorry if I will not answer very soon likely, but I'm very very busy on these days. Matteo

 June 11, 2014, 06:40 error: ‘class Foam::basicRhoThermo’ has no member named #12 New Member   Duarte Magalhães Join Date: Apr 2014 Location: Lisbon, Portugal Posts: 24 Rep Power: 10 Hi everyone! First of all thanks a lot Matteo for this thread, it was really helpful to have a starting point. I am implementing this code to my icoFoam code (already with temperature equation as in http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam and with Courant number control). OpenFOAM version is 2.3.0. While compiling i get the following errors: createFields.H:62:30: error: ‘class Foam::basicRhoThermo’ has no member named ‘h’ createFields.H:91:12: error: ‘class Foam::basicRhoThermo’ has no member named ‘mu’ own_icoFoamPV.C:121:15: error: ‘class Foam::basicRhoThermo’ has no member named ‘mu’ I think i have done everything as explained here, and i can't understand why I am getting this error. I tried to use rhoThermo instead of basicRhoThermo but i still get the error: createFields.H:62:30: error: ‘class Foam::rhoThermo’ has no member named ‘h’ Any help is much appreciated! Thanks in advance! Last edited by DuarteMagalhaes; June 11, 2014 at 10:56.

 June 11, 2014, 12:35 Call functions from rhoThermo and basicThermo in createFields.H file #13 New Member   Duarte Magalhães Join Date: Apr 2014 Location: Lisbon, Portugal Posts: 24 Rep Power: 10 Ok, I fixed this problem. The problem was that i was using class rhoThermo to call the function h( ) but this function is not present in rhoThermo (only basicThermo, for example) and also in OpenFOAM 2.3.0. it is he( ), not h( ). I have another problem now. On my createFields.H file, i will need to create fields from two different classes inside thermophysicalProperties library because basicThermo has all the functions i need (rho, he, kappa) but not mu( ) which is present in class rhoThermo. How can i make my createFields.H file so that it reads from both classes? Thank you in advance! Last edited by DuarteMagalhaes; June 11, 2014 at 14:07.

 September 9, 2014, 11:27 #14 New Member   Wentao Zheng Join Date: Nov 2013 Posts: 7 Rep Power: 10 hi. I have a question. What's the meaning of the “thermo.correct”? I also find "thermo.correct" in rhoSimpleFoam. How does it get the temperature field? Why does it not use " h = Cp*T" to solve the temperature field?

September 15, 2014, 11:06
#15
Member

Join Date: Feb 2012
Posts: 35
Rep Power: 12
Quote:
 Originally Posted by 774024952 hi. I have a question. What's the meaning of the “thermo.correct”? I also find "thermo.correct" in rhoSimpleFoam. How does it get the temperature field? Why does it not use " h = Cp*T" to solve the temperature field?
The reason of this expression is to update all the variables defined into the Thermo Model it has been adopted. In my case, basicRhoThermo, does not have enthalpy as a variable automatically updated with a "thermo.correct" command but others as the density. To better understand how this process happens you could simply just do a "gdb yourCase" to debug your code and see how it works.

Greetings.

 September 16, 2014, 09:50 #16 New Member   Wentao Zheng Join Date: Nov 2013 Posts: 7 Rep Power: 10 Thanks for your reply showed me a way to debug in linux. I need to learn some basic knowledge about the linux and OpenFOAM.

 December 30, 2016, 07:51 #17 Member   Saurav Kumar Join Date: Jul 2016 Posts: 80 Rep Power: 7 Hi, i just want to use pisoFoam for heat transfer in laminar and turbulent flows without thermal effects. Matt_B used energy equation and considered thermal effects. but i dont want to consider thermal effect then what modification i should do in energy equation. Here is MATT_B energy equation with thermal effects. ************************************************** ************** kappat = turbulence->nut()/Prt; kappat.correctBoundaryConditions(); volScalarField kappaEff("kappaEff", turbulence->nu()/Pr + kappat); fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(); thermo.correct(); rho=thermo.rho(); *********************************************** Here is my energy equation without thermal effect. ************************************************** ************* fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T) ); TEqn.relax(); TEqn.solve(); eqnResidual = TEqn.solve().initialResidual(); maxResidual = max(eqnResidual, maxResidual); ************************************************** *********************************8 is it right? will it work for laminar as well as turbulent flow? raj kumar saini likes this.