
[Sponsors] 
Problem in modified pisoFoam with temperature equation + thermophysical model 

LinkBack  Thread Tools  Display Modes 
August 29, 2012, 09:35 
Problem in modified pisoFoam with temperature equation + thermophysical model

#1 
Member
Join Date: Feb 2012
Posts: 35
Rep Power: 6 
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<pureMixture<icoPoly8ThermoPhysics>>; 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 1e5 0 0 0 0 0 0); } } Code:
transportModel Newtonian; nu nu [ 0 2 1 0 0 0 0 ] 1e05; // 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; 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" // Pressurevelocity 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; corr<nCorr; corr++) { volScalarField rAU(1.0/UEqn.A()); U = rAU*UEqn.H(); phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, U, phi); adjustPhi(phi, U, p); // Nonorthogonal pressure corrector loop for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { // Pressure corrector fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); if ( corr == nCorr1 && nonOrth == nNonOrthCorr ) { pEqn.solve(mesh.solver("pFinal")); } else { pEqn.solve(); } if (nonOrth == nNonOrthCorr) { phi = pEqn.flux(); } } #include "continuityErrs.H" U = rAU*fvc::grad(p); U.correctBoundaryConditions(); } } turbulence>correct(); #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<incompressible::turbulenceModel> turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Info<< "Reading thermophysical properties\n" << endl; autoPtr<basicRhoThermo> 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(); 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(); Thank you. Matteo 

August 30, 2012, 11:44 

#2 
Member
Join Date: Feb 2012
Posts: 35
Rep Power: 6 
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 ); Code:
Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), thermo.T() ); Please, I need someone's help!! 

September 5, 2012, 12:14 
Issues resolved!

#3 
Member
Join Date: Feb 2012
Posts: 35
Rep Power: 6 
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<pureMixture<icoPoly8ThermoPhysics>>; 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 1e5 0 0 0 0 0 0); Liquid2; muCoeffs<8> (1 0.0008 0.0000007 0.0000000001 0 0 0 0); kappaCoeffs<8> ( 1 1e5 0 0 0 0 0 0); } } 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<incompressible::turbulenceModel> turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Info<< "Reading thermophysical properties\n" << endl; autoPtr<basicRhoThermo> 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 ); 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; } 

November 20, 2012, 09:47 

#4 
New Member
Markus Trompa
Join Date: Nov 2012
Location: Regensburg, Germany
Posts: 13
Rep Power: 5 
Hello Matt_B,
thank you very much for uploading your files. With your topic I could made a big step ahead with my problem. I want to simulate inlet flow jet in a room with a heated box and so had to implement TEqn in pisoFoam too or otherwise use heatTransfer solver with some changes. But now i have some problem when executing the solver pisoFoam. Could you please upload your fvSchemes, fvSolution and controlDict. You would do my a big favor. Thank you very much in advance Markus 

November 20, 2012, 14:22 

#5 
Member
Join Date: Feb 2012
Posts: 35
Rep Power: 6 
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((1A(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 ; } Code:
solvers { p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.1; } pFinal { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0; } U { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } T { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } k { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } epsilon { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } R { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } nuTilda { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } 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; Matteo 

November 20, 2012, 14:41 

#6 
New Member
Markus Trompa
Join Date: Nov 2012
Location: Regensburg, Germany
Posts: 13
Rep Power: 5 
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: 5 
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: 6 
Quote:
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: 5 
Here is a pic of my case. It is a crosssection 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 heatsource 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
Konrad
Join Date: Sep 2013
Posts: 4
Rep Power: 5 
I'm getting runs like this:
Create time Create mesh for time = 0 Reading field p Reading field U Reading/calculating face flux field phi 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: 2e05 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; } gradSchemes { default Gauss linear; grad(p) Gauss linear; } divSchemes { default none; div(phi,U) Gauss linear; div((nuEff*dev(T(grad(U))))) Gauss linear; div(phi,T) Gauss upwind; } laplacianSchemes { default none; laplacian(nu,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; laplacian(DT,T) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(HbyA) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } Code:
FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.05; } pFinal { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0; } T { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0.01; nSweeps 3; maxIter 100; minIter 10; } U { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } k { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } B { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } nuTilda { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } 

October 23, 2013, 08:46 

#11 
Member
Join Date: Feb 2012
Posts: 35
Rep Power: 6 
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: 4 
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: 4 
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: 4 
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: 6 
Quote:
Greetings. 

September 16, 2014, 09:50 

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


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Use of kepsilon and komega Models  Jade M  Main CFD Forum  14  September 7, 2016 14:14 
Mixture model problem  could someone please advise?  matlab_monkey  FLUENT  2  July 26, 2012 08:20 
Poor convergence with temperature dependent density on modified pisoFOAM  ovie  OpenFOAM  1  March 20, 2011 04:19 
Problem with Joulebs effect source term in the energy equation  galaad  OpenFOAM Running, Solving & CFD  0  January 19, 2006 13:01 
Writing a BCDEFI problem for RSM model  S. Bottenheim  CDadapco  2  January 28, 2005 09:55 