CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Another attempt at adding temperature to simpleFoam (http://www.cfd-online.com/Forums/openfoam-programming-development/112376-another-attempt-adding-temperature-simplefoam.html)

sdharmar January 22, 2013 14:27

Another attempt at adding temperature to simpleFoam
 
2 Attachment(s)
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

2 Attachment(s)
Attachment 18451

Attachment 18452

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

the solver
 
5 Attachment(s)
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

My simpleFoam
 
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 (Post 351117)
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 (Post 404356)
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" :rolleyes:

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

2 Attachment(s)
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:
  1. 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.
  2. 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

Some corrections
 
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

1 Attachment(s)
Quote:

Originally Posted by sdharmar (Post 403854)
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


All times are GMT -4. The time now is 21:51.