Another attempt at adding temperature to simpleFoam

sdharmar · January 22, 2013, 14:27

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 k-epsilon 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((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;
laplacian(DT,T) Gauss linear corrected;
}

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

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}

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

Bernhard · January 24, 2013, 10:59

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?

sdharmar · January 24, 2013, 14:16

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.

sdharmar · January 24, 2013, 14:28

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.

sdharmar · January 27, 2013, 14:52

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/suranga-2.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)
    );

And then changed then created a file TEqn.H in tempSimpleFoam directory. (actually I cloned it from some other code but I don't remember what the original application is)

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);
}

Then tempSimpleFoam.C was changed as shown below.

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
    Steady-state 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;

        // --- Pressure-velocity 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;
}


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

Then in Make/files I change the following.

Code:

tempSimpleFoam.C

EXE = $(FOAM_USER_APPBIN)/tempSimpleFoam

Then I gave the command

Code:

wmake

in

/home/suranga/OpenFOAM/suranga-2.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.

wyldckat · January 27, 2013, 17:39

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: http://www.cfd-online.com/Forums/ope...-get-help.html

Quote:

Originally Posted by niklas

Don't hijack other threads and don't post the same question multiple times.
I see that some people post new questions in active threads cause they think this will attract more readers. Do not do this, since it makes it really hard to search the forum.

About what you wrote in the previous post:

Quote:

Originally Posted by sdharmar

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.

I think I was already too tired yesterday to be upset

But put yourself in my shoes: how would you feel if someone simply asked you what is the specific and concrete value to give to "z", so that the following equation is true: "x + y + z = 1"

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:

wyldckat · January 27, 2013, 18:23

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:

tempSimpleFoam.tar.gz - solver modified somewhat based on your instructions. There is a discrepancy, because you posted a "createFields.H" that refers to "PISO" near the end of the file, but I used the original "SIMPLE".
suranga_temp.tar.gz - based on the tutorial "incompressible/simpleFoam/pitzDaily" and on the files you provided.

The final iteration 1000 had this:

Code:

DILUPBiCG:  Solving for Ux, Initial residual = 9.3952e-006, Final residual = 3.09815e-008, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 9.72644e-005, Final residual = 2.51259e-007, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.000164337, Final residual = 1.4318e-006, No Iterations 14
time step continuity errors : sum local = 7.01076e-006, global = -1.60798e-007, cumulative = -0.0238757
DILUPBiCG:  Solving for T, Initial residual = 0.0313769, Final residual = 8.66675e-005, No Iterations 99
DILUPBiCG:  Solving for epsilon, Initial residual = 1.40552e-005, Final residual = 1.13521e-007, No Iterations 4
DILUPBiCG:  Solving for k, Initial residual = 3.8993e-005, Final residual = 3.73099e-007, No Iterations 3

The residual for T looks rather high, but that's possibly because things were not properly defined somewhere...

So, I suggest that you do the following steps:

If you cannot share the case you're working on, then create a case based on a tutorial that has a similar geometry to yours, similarly to what I did. Then share that case with us here on this thread.
If you do not have convergence problems with that test case, then check the sanity of the mesh on your current case:
Code:
```
checkMesh
checkMesh -allTopology -allGeometry
```

Best regards,
Bruno

sdharmar · January 27, 2013, 21:07

Hi Bruno,

Thank you very much for your valuable reply and your decision to move my post.

Quote:

tempSimpleFoam.tar.gz - solver modified somewhat based on your instructions. There is a discrepancy, because you posted a "createFields.H" that refers to "PISO" near the end of the file, but I used the original "SIMPLE".

I am sorry, It seems I have attached details of the wrong file. This is the correct createFields.H which is created for tempSimpleFoam solver.

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);

And I missed to mention that I defined Pr and Prt in transportProperties dictionary in the case/system directory. It is shown below.

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 ] 1e-05;
DT              DT [ 0 2 -1 0 0 0 0 ] 2e-05;
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 ] 1e-06;
    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
    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 ] 1e-06;
    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
    k               k [ 0 0 1 0 0 0 0 ] 0;
    n               n [ 0 0 0 0 0 0 0 ] 1;
}


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

for the completeness of the question I have given all the other files in

/home/naren/OpenFOAM/naren-2.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);

    // Non-orthogonal 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);
}

tempSimpleFoam/UEqn.H

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));

tempSimpleFoam/TEqn.H

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);
}

tempSimpleFoam/tempSimpleFoam.C

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
    Steady-state 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;

        // --- Pressure-velocity 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;
}


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

And changes in the file Make/file has been done as in previously explained. Please be kind to give your feedback.

BR,

Suranga,

wyldckat · January 28, 2013, 05:51

Hi Suranga,

Please read post #7, if you haven't already

Best regards,
Bruno

Sunxing · March 21, 2013, 21:42

Quote:

Originally Posted by sdharmar

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.

Hi sdharmar,

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();

Am I right to modify it like this?

Best regards

Sunxing · March 24, 2013, 05:16

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

lukasf · March 16, 2020, 08:45

This post explains the whole idea:

https://www.factorai.io/2019/08/16/a...to-pimplefoam/

January 28, 2013, 05:51		#9
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,982 Blog Entries: 45 Rep Power: 128	Hi Suranga, Please read post #7, if you haven't already Best regards, Bruno __________________ OpenFOAM: FAQ \| Getting started Forum: How to get help, to post code/output and forum guide Read this before sending me PM

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
interFoam vs. simpleFoam channel flow comparison	DanM	OpenFOAM Running, Solving & CFD	12	January 31, 2020 16:26
Calculation of the Governing Equations	Mihail	CFX	7	September 7, 2014 07:27
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

January 24, 2013, 10:59		#2
Bernhard Senior Member Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 22	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
sdharmar Member Suranga Dharmarathne Join Date: Jan 2011 Location: TX, USA Posts: 39 Rep Power: 15	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.

March 24, 2013, 05:16		#11
Sunxing Member sqing Join Date: Sep 2012 Location: Dalian Posts: 77 Rep Power: 14	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

March 16, 2020, 08:45		#12
lukasf Senior Member Lukas Fischer Join Date: May 2018 Location: Germany, Munich Posts: 117 Rep Power: 8	This post explains the whole idea: https://www.factorai.io/2019/08/16/a...to-pimplefoam/