[Matt_B]

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:

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 1e-5 0 0 0 0 0 0);
    }
}

transportProperties dictionary:

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:

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

                // Non-orthogonal 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 == nCorr-1
                     && 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:

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

TEqn.H:

    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

[Matt_B]

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:

    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:

    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!!

[Matt_B]

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:

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 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:

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

TEqn.H:

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

}

[MisterX]

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

[Matt_B]

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:

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:

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:

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

[MisterX]

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

[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

[Matt_B]

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

Hi 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

[Byxon]

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

[Byxon]

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: 2e-05 average: 0
oneEqEddyCoeffs
{
ce 1.048;
ck 0.094;
}


Starting time loop

End

I have attached some of my files.

FvSchemes:

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((1|A(U)),p) Gauss linear corrected;

laplacian(DT,T) Gauss linear corrected;
}

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

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

fvSolutions:

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

[Matt_B]

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

[DuarteMagalhaes]

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!

[DuarteMagalhaes]

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!

[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?

[Matt_B]

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.

[774024952]

Thanks for your reply showed me a way to debug in linux. I need to learn some basic knowledge about the linux and OpenFOAM.

[srv537]

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?

[Haijing]

Hi,

I'm a beginner about OpenFOAM. Recrntly, I tried to do a model about thermal internal boundary layer in a box. But when I run the command 'wmake',there are always error hints below:

UEqn.H(14): error: identifier "fvOptions" is undefined
+ fvOptions(U)
^

In file included from ras_pisoFoam_heating.C(110):
UEqn.H(19): error: expression must have class type
fvOptions.constrain(UEqn);
^

In file included from ras_pisoFoam_heating.C(110):
UEqn.H(25): error: expression must have class type
fvOptions.correct(U);
^

In file included from ras_pisoFoam_heating.C(115):
pEqn.H(41): error: expression must have class type
fvOptions.correct(U);
^

ras_pisoFoam_heating.C(138): error: identifier "fvOptions" is undefined
fvOptions(T)
^

ras_pisoFoam_heating.C(141): error: expression must have class type
fvOptions.constrain(TEqn);
^

ras_pisoFoam_heating.C(143): error: expression must have class type
fvOptions.correct(T);
^

compilation aborted for ras_pisoFoam_heating.C (code 2)
make: *** [Make/linux64IccDPInt64Opt/ras_pisoFoam_heating.o] Error 2


I'll poster my code, and wish you could give me some comments, thank you very much.

Haijing

UEqn.H
// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
extPGrad
+ ((T-T0)/T0)*extBuoyancy
+ fvOptions(U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));

fvOptions.correct(U);
}

ras_pisoFoam_heating.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 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
pisoFoam

Group
grpIncompressibleSolvers

Description
Transient solver for incompressible, turbulent flow, using the PISO
algorithm.

\heading Solver details
The solver uses the PISO algorithm to solve the continuity equation:

\f[
\div \vec{U} = 0
\f]

and momentum equation:

\f[
\ddt{\vec{U}} + \div \left( \vec{U} \vec{U} \right) - \div \gvec{R}
= - \grad p
\f]

Where:
\vartable
\vec{U} | Velocity
p | Pressure
\vec{R} | Stress tensor
\endvartable

Sub-models include:
- turbulence modelling, i.e. laminar, RAS or LES
- run-time selectable MRF and finite volume options, e.g. explicit porosity

\heading Required fields
\plaintable
U | Velocity [m/s]
p | Kinematic pressure, p/rho [m2/s2]
\<turbulence fields\> | As required by user selection
\endplaintable

\*---------------------------------------------------------------------------*/


#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"

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

int main(int argc, char *argv[])
{
#include "postProcess.H"

#include "addCheckCaseOptions.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "initContinuityErrs.H"

turbulence->validate();

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

Info<< "\nStarting time loop\n" << endl;
/*
volVectorField centres = U.mesh().C();
forAll (CellCentre,celli)
{
CellCentre[celli]=centres[celli];
}
*/

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "CourantNo.H"

// Pressure-velocity PISO corrector
{
#include "UEqn.H"

// --- PISO loop
while (piso.correct())
{
#include "pEqn.H"
}

// Reynolds stress in RANS
/*
volTensorField gradU = fvc::grad(U);
Tau=-turbulence->nut()*(gradU+gradU.T());
*/


}

laminarTransport.correct();
turbulence->correct();


//T: Temperature
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(turbulence->nuEff(),T)
==
fvOptions(T)
);
TEqn.relax();
fvOptions.constrain(TEqn);
TEqn.solve();
fvOptions.correct(T);


/*
//T0 transport: pollutants source at surface change position
fvScalarMatrix T0Eqn
(
fvm::ddt(T0)
+ fvm::div(phi, T0)
- fvm::laplacian(turbulence->nuEff(),T0)
- TracerSource0
==
fvOptions(T0)
);
T0Eqn.relax();
fvOptions.constrain(T0Eqn);
T0Eqn.solve();
fvOptions.correct(T0);



//T1 transport: pollutants source at surface change position
fvScalarMatrix T1Eqn
(
fvm::ddt(T1)
+ fvm::div(phi, T1)
- fvm::laplacian(turbulence->nuEff(),T1)
- TracerSource1
==
fvOptions(T1)
);
T1Eqn.relax();
fvOptions.constrain(T1Eqn);
T1Eqn.solve();
fvOptions.correct(T1);

//T2 transport: pollutants source after surface change position
fvScalarMatrix T2Eqn
(
fvm::ddt(T2)
+ fvm::div(phi, T2)
- fvm::laplacian(turbulence->nuEff(),T2)
- TracerSource2
==
fvOptions(T2)
);
T2Eqn.relax();
fvOptions.constrain(T2Eqn);
T2Eqn.solve();
fvOptions.correct(T2);
*/



runTime.write();

// #include "patch.H"

//Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
//<< " ClockTime = " << runTime.elapsedClockTime() << " s"
//<< nl << endl;
runTime.printExecutionTime(Info);
}

Info<< "End\n" << endl;

return 0;
}


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

pEqn.H
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr( U, phi))
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU, MRF);

// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector

fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

pEqn.setReference(pRefCell, pRefValue);

pEqn.solve(mesh.solver(p.select(piso.finalInnerIte r())));

if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}

#include "continuityErrs.H"

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

ras_pisoFoam_heating.C