
[Sponsors] 
January 22, 2013, 14:27 
Another attempt at adding temperature to simpleFoam

#1 
Member
Suranga Dharmarathne
Join Date: Jan 2011
Location: TX, USA
Posts: 39
Rep Power: 6 
Dear formers,
I use following fvSchemes for a film cooling problem. As attached plots show laminar case of the same problem converges. But I still have concern about the the peak of the pressure residuals. The turbulent case with standard kepsilon model does not converge. Dose anyone of you can give me any idea. Your help is much appreciated. I am burning weeks for this problem now. Best regards, Suranga. /** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 2.1.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes { default none; div(phi,U) Gauss upwind; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; //div(phi,R) Gauss upwind; //div(R) Gauss linear; //div(phi,nuTilda) Gauss upwind; div((nuEff*dev(T(grad(U))))) Gauss linear; div(phi,T) Gauss upwind; } laplacianSchemes { default none; laplacian(nuEff,U) 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; laplacian(DT,T) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************** *********************** // 

January 24, 2013, 10:59 

#2 
Senior Member
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12 
It is unlikely that you will find the problem in the fvSchemes. Also, which solver are you using? What are your tolerances in fvSolution? And, which residuals are you plotting here? Initial, final or what? What are the boundary conditions?


January 24, 2013, 14:16 

#3 
Member
Suranga Dharmarathne
Join Date: Jan 2011
Location: TX, USA
Posts: 39
Rep Power: 6 
0.tar.gz
system.tar.gz Hi Bernhard, Thanks for your reply. Here with I have attached 0 and system directories of my case. Please be kind to go through these and help me. I uses pyFoamPlotWatcher.py to plot residuals and it plots initial residuals. BR, Suranga. 

January 24, 2013, 14:28 
the solver

#4 
Member
Suranga Dharmarathne
Join Date: Jan 2011
Location: TX, USA
Posts: 39
Rep Power: 6 
I am sorry I couldn't attach the solver in the last post. Here it is. I use modified simpleFoam by adding temperature equation to the general simpleFoam.


January 27, 2013, 14:52 
My simpleFoam

#5 
Member
Suranga Dharmarathne
Join Date: Jan 2011
Location: TX, USA
Posts: 39
Rep Power: 6 
Hi Bruno,
I am very sorry if I interpreted the problem wrong. It seems that you got little upset. I am sorry about that. I am also very inexperience in OpenFoam. And I don't understand most of the things in codes. And I am pretty new to CFD too. That is why I have been seeking help in the forum. I just wanted to see whether I am doing right.Anyway here is what I did. I copied Make directory, pEqn.H, UEqn.H, createFields.H and simpleFoam.C files form /opt/openfoam211/applications/solvers/incompressible/simpleFoam to /home/suranga/OpenFOAM/suranga2.1.1/applications/solvers/tempSimpleFoam and changed the file name simpleFoam.C to tempSimpleFoam.C. Then I modified createFields.H as shown below. Modifications are shown in red. Code:
Info<< "Reading field p\n" << endl; volScalarField p ( IOobject ( "p", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); // temperature field Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); //end temperature field Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", 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); // laminar prandtl number [Pr] dimensionedScalar Pr(laminarTransport.lookup("Pr")); //turbulent Prandtl number [Prt] dimensionedScalar Prt(laminarTransport.lookup("Prt")); autoPtr<incompressible::turbulenceModel> turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Code:
{ volScalarField kappaEff ( "kappaEff", turbulence>nu()/Pr + turbulence>nut()/Prt ); fvScalarMatrix TEqn ( //fvm::ddt(T) // changed this since simpleFoam is steady state. fvm::div(phi, T)  fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(); //rhok = 1.0  beta*(T  TRef); } Code:
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 2011 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Application simpleFoam Description Steadystate solver for incompressible, turbulent flow \**/ #include "fvCFD.H" #include "singlePhaseTransportModel.H" #include "RASModel.H" #include "simpleControl.H" #include "IObasicSourceList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "createFields.H" #include "initContinuityErrs.H" simpleControl simple(mesh); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (simple.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; //  Pressurevelocity SIMPLE corrector { #include "UEqn.H" #include "pEqn.H" #include "TEqn.H" // this is to solve temperature equation. } turbulence>correct(); runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // Code:
tempSimpleFoam.C EXE = $(FOAM_USER_APPBIN)/tempSimpleFoam Code:
wmake /home/suranga/OpenFOAM/suranga2.1.1/applications/solvers/tempSimpleFoam directory. Please don't hesitate to try this and give your comments on this. Your valuable ideas are very welcome. Best regards, Suranga. 

January 27, 2013, 17:39 

#6  
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,301
Blog Entries: 34
Rep Power: 84 
Hi Suranga,
I've moved some of your posts related to this problem into a single thread. The reason is rather simple, as explained here: How to give enough info to get help Quote:
Quote:
Now, as for the problem you have... I'm going to have to think a bit on this, gather the information on the whole thread and I'll reply back when I have any news. In the mean time, if anyone else wants to answer to Suranga, feel free to do so! Best regards, Bruno edit: the posts were moved from these two threads:
__________________
Last edited by wyldckat; January 27, 2013 at 17:43. Reason: see "edit:"  and fixed a few typos on the hypothetical 

January 27, 2013, 18:23 

#7 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,301
Blog Entries: 34
Rep Power: 84 
Hi Suranga,
I don't know what is the problem you're getting, because I was not able to reproduce it. Apparently I'm still missing either "x" or "y", in order to tell you what "z" is! Attached are the following files:
Code:
DILUPBiCG: Solving for Ux, Initial residual = 9.3952e006, Final residual = 3.09815e008, No Iterations 2 DILUPBiCG: Solving for Uy, Initial residual = 9.72644e005, Final residual = 2.51259e007, No Iterations 2 GAMG: Solving for p, Initial residual = 0.000164337, Final residual = 1.4318e006, No Iterations 14 time step continuity errors : sum local = 7.01076e006, global = 1.60798e007, cumulative = 0.0238757 DILUPBiCG: Solving for T, Initial residual = 0.0313769, Final residual = 8.66675e005, No Iterations 99 DILUPBiCG: Solving for epsilon, Initial residual = 1.40552e005, Final residual = 1.13521e007, No Iterations 4 DILUPBiCG: Solving for k, Initial residual = 3.8993e005, Final residual = 3.73099e007, No Iterations 3 So, I suggest that you do the following steps:
Bruno
__________________


January 27, 2013, 21:07 
Some corrections

#8  
Member
Suranga Dharmarathne
Join Date: Jan 2011
Location: TX, USA
Posts: 39
Rep Power: 6 
Hi Bruno,
Thank you very much for your valuable reply and your decision to move my post. Quote:
Code:
Info<< "Reading field p\n" << endl; volScalarField p ( IOobject ( "p", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); // temperature field Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); //end temperature field Info<< "Reading field U\n" << endl; volVectorField U ( IOobject ( "U", 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("SIMPLE"), pRefCell, pRefValue); singlePhaseTransportModel laminarTransport(U, phi); // laminar prandtl number [Pr] dimensionedScalar Pr(laminarTransport.lookup("Pr")); //turbulent Prandtl number [Prt] dimensionedScalar Prt(laminarTransport.lookup("Prt")); autoPtr<incompressible::RASModel> turbulence ( incompressible::RASModel::New(U, phi, laminarTransport) ); IObasicSourceList sources(mesh); Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 2.1.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object transportProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // transportModel Newtonian; nu nu [ 0 2 1 0 0 0 0 ] 1e05; DT DT [ 0 2 1 0 0 0 0 ] 2e05; Pr Pr [ 0 0 0 0 0 0 0 ] 0.7; Prt Prt [ 0 0 0 0 0 0 0 ] 0.85; CrossPowerLawCoeffs { nu0 nu0 [ 0 2 1 0 0 0 0 ] 1e06; nuInf nuInf [ 0 2 1 0 0 0 0 ] 1e06; m m [ 0 0 1 0 0 0 0 ] 1; n n [ 0 0 0 0 0 0 0 ] 1; } BirdCarreauCoeffs { nu0 nu0 [ 0 2 1 0 0 0 0 ] 1e06; nuInf nuInf [ 0 2 1 0 0 0 0 ] 1e06; k k [ 0 0 1 0 0 0 0 ] 0; n n [ 0 0 0 0 0 0 0 ] 1; } // ************************************************************************* // /home/naren/OpenFOAM/naren2.1.1/applications/solvers/tempSimpleFoam diractory. tempSimpleFoam/pEqn.H Code:
{ p.boundaryField().updateCoeffs(); volScalarField rAU(1.0/UEqn().A()); U = rAU*UEqn().H(); UEqn.clear(); phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf(); adjustPhi(phi, U, p); // Nonorthogonal pressure corrector loop while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve(); if (simple.finalNonOrthogonalIter()) { phi = pEqn.flux(); } } #include "continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Momentum corrector U = rAU*fvc::grad(p); U.correctBoundaryConditions(); sources.correct(U); } Code:
// Momentum predictor tmp<fvVectorMatrix> UEqn ( fvm::div(phi, U) + turbulence>divDevReff(U) == sources(U) ); UEqn().relax(); sources.constrain(UEqn()); solve(UEqn() == fvc::grad(p)); Code:
{ volScalarField kappaEff ( "kappaEff", turbulence>nu()/Pr + turbulence>nut()/Prt ); fvScalarMatrix TEqn ( //fvm::ddt(T) fvm::div(phi, T)  fvm::laplacian(kappaEff, T) ); TEqn.relax(); //TEqn.solve(); TEqn.solve().initialResidual(); //rhok = 1.0  beta*(T  TRef); } Code:
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 2011 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Application simpleFoam Description Steadystate solver for incompressible, turbulent flow \**/ #include "fvCFD.H" #include "singlePhaseTransportModel.H" #include "RASModel.H" #include "simpleControl.H" #include "IObasicSourceList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "createFields.H" #include "initContinuityErrs.H" simpleControl simple(mesh); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (simple.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; //  Pressurevelocity SIMPLE corrector { #include "UEqn.H" #include "pEqn.H" #include "TEqn.H" } turbulence>correct(); runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // BR, Suranga, 

January 28, 2013, 05:51 

#9 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,301
Blog Entries: 34
Rep Power: 84 
Hi Suranga,
Please read post #7, if you haven't already Best regards, Bruno
__________________


March 21, 2013, 21:42 

#10  
Member
sqing
Join Date: Sep 2012
Location: Dalian
Posts: 77
Rep Power: 4 
Quote:
I have read the TEqn.H that your posted. However I have not understood clearly with your TEqn. What is the kappaEff? Can you give me a mathematical equation to describe it? And I have modified the pisoFoam solver to simulate my film cooling case. This is my T equation: Code:
fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T)  fvm::laplacian(DT, T) ); TEqn.solve(); Best regards 

March 24, 2013, 05:16 

#11 
Member
sqing
Join Date: Sep 2012
Location: Dalian
Posts: 77
Rep Power: 4 
Hi sdharmar,
I have another question to you. How do you set the value of Prt? Do you have a formula to figure it out? Regards, Sunxing 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Calculation of the Governing Equations  Mihail  CFX  7  September 7, 2014 06:27 
interFoam vs. simpleFoam channel flow comparison  DanM  OpenFOAM Running, Solving & CFD  11  January 5, 2013 07:21 
Is wall ajacent temperature equal to conservative temperature of the wall?  shenying0710  CFX  8  January 4, 2013 05:03 
Trying to run a benchmark case with simpleFoam  spsb  OpenFOAM  3  February 24, 2012 10:07 
difference of temperature  HKH  FLUENT  0  January 28, 2010 03:45 