[Red Devil]

Hi All,

If someone has successfully added temperature to pisoFoam solver, I request you to share your Solver folder (one with *.C and createFields.H files) and any one of your case folders (with '0', 'Constant' and 'System' folders) which you have successfully run.

I have successfully added temperature to icoFoam using the instructions given here:

http://openfoamwiki.net/index.php/How_to_add_temperature_to_icoFoam

.
Using the same procedure, I added temperature to pisoFoam also and the solver runs fine. I just want to confirm the correctness of my solver. Please provide some help.

Thanks in advance.

[sharonyue]

Createfields

    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
);
   
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

dimensionedScalar DT
(
transportProperties.lookup("DT")
);
#   include "createPhi.H"


    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);


    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr<incompressible::turbulenceModel> turbulence
    (
        incompressible::turbulenceModel::New(U, phi, laminarTransport)
    );

solver:

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

                volVectorField HbyA("HbyA", U);
                HbyA = rAU*UEqn.H();
                surfaceScalarField phiHbyA
                (
                    "phiHbyA",
                    (fvc::interpolate(HbyA) & mesh.Sf())
                  + fvc::ddtPhiCorr(rAU, U, phi)
                );

                adjustPhi(phiHbyA, U, p);

                // Non-orthogonal pressure corrector loop
                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                {
                    // Pressure corrector

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

                    pEqn.setReference(pRefCell, pRefValue);

                    if
                    (
                        corr == nCorr-1
                     && nonOrth == nNonOrthCorr
                    )
                    {
                        pEqn.solve(mesh.solver("pFinal"));
                    }
                    else
                    {
                        pEqn.solve();
                    }

                    if (nonOrth == nNonOrthCorr)
                    {
                        phi = phiHbyA - pEqn.flux();
                    }
                }

                #include "continuityErrs.H"

                U = HbyA - rAU*fvc::grad(p);
                U.correctBoundaryConditions();
		fvScalarMatrix TEqn
		(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);


TEqn.solve();
            }
        }

        turbulence->correct();

        runTime.write();

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

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

    return 0;
}

Is this right?

[Red Devil]

Thanks for posting your files sharonyue.

I have one change in my solver file. I kept the TEqn outside the piso loop. Which one is correct?

Also, I've seen many people using kappat while adding temperature to pisoFoam. Please help me to understand kappat and it's significance

[alexeym]

Hi,

would you like to reimplement buoyantBoussinesqPimpleFoam (with thermal expansion coefficient equal to zero, as you did not modify momentum equations)? In that solver temperature equation is solver just after velocity equation.

UPD. kappat is turbulent heat conductivity. If you're planning to simulate heat transfer in turbulent flow, then you should account for it, otherwise, you can forget about it.

[sharonyue]

Originally Posted by alexeym

Hi,

would you like to reimplement buoyantBoussinesqPimpleFoam (with thermal expansion coefficient equal to zero, as you did not modify momentum equations)? In that solver temperature equation is solver just after velocity equation.

UPD. kappat is turbulent heat conductivity. If you're planning to simulate heat transfer in turbulent flow, then you should account for it, otherwise, you can forget about it.

Absolutely rite! I once made a comparison with this TpisoFoam and buoyantBoussinesqPimpleFoam with CFX, the latter is more accurate. Maybe we should consider this kappat into a turbulent solver.

And whether to put it in or out piso loop(after turbulence correction). I just made a comparison too. It is totally the same without any difference... See the picture. So I think put it out of the loop is better to save more calculating time.

Regards,

[Red Devil]

Thanks alexeym and sharonyue for your inputs.
I'm trying to create a turbulent solver to solve forced convection heat transfer problem. What all changes I have to make for that?

[alexeym]

Hi,

if you'd like to create this solver for the sake of creation of the solver, then you should also add buoyancy force to momentum equation.

If you'd like to have this solver for practical reasons (running simulations, obtaining results) then you can just use buoyantBousinesqPimpleFoam. And if you like it to behave like "pisoFoam with temperature and buoyancy" set number of outer correctors (nOuterCorrectors in PIMPLE dictionary) to 1.

[Red Devil]

Hi alexeym,
I'm not considering buoyancy in my problem.

[alexeym]

OK. Set beta (thermal expansion coefficient) to zero and there will be no buoyancy in buoyantBoussinesqPimpleFoam.

Or if you'd like to go with "your own" solver, here is how turbulent heat diffusivity is calculated:

    alphat = turbulence->nut()/Prt;
    alphat.correctBoundaryConditions();

    volScalarField alphaEff("alphaEff", turbulence->nu()/Pr + alphat);

then alphaEff is used in laplacian of temperature equation.

[Red Devil]

Thanks a lot alexeym. Will try that.

[Red Devil]

Hi alexeym,
I've successfully added temperature to pisoFoam by following your inputs. Thanks a lot.

I've another question. I want to set a constant heat flux (q w/m2) boundary condition on a wall. I did it this way:

    wall
    {
        type      fixedGradient;
        gradient      uniform 100;
    }

where the uniform value of temperature gradient is calculated as q/k (k-thermal conductivity).
Is this correct?

[alexeym]

Hi,

no, it won't be correct, you need heat transfer coefficient (http://en.wikipedia.org/wiki/Heat_transfer_coefficient).

[Red Devil]

Oh! Then how to implement this constant heat flux bc? I'm totally confused.

[alexeym]

Well...

start with description of the problem, as "...forced convection heat transfer problem..." is not enough (you'd like to use q = hc*(T - Tinf) relation for calculation of heat flux, it OK, but to calculate that hc one needs to know more about problem).

Then

$ cd $FOAM_SRC
$ find . -iname '*wallHeat*'

this will give you possible starting points for implementation of your BC.

[Red Devil]

Problem is flow of water through a square duct. At one wall of the duct, I want to give constant heat flux. Want to run both laminar and turbulent cases. For laminar case i'm using modified icoFoam and for turbulent, modified pisoFoam.

So I have to find h for my problem, then from q=h*(T-Tinf) I can get T.
Then the temperature bc is a uniform fixedValue? Is that the way?

[alexeym]

Originally Posted by Red Devil

Problem is flow of water through a square duct. At one wall of the duct, I want to give constant heat flux. Want to run both laminar and turbulent cases. For laminar case i'm using modified icoFoam and for turbulent, modified pisoFoam.

If you just go with buoyantBoussinesqPimpleFoam (with beta equal zero cause you don't need buoyancy), you can avoid this zoo of solvers, as it is capable of simulating momentum and heat transfer in laminar and turbulent regimes.

So I have to find h for my problem, then from q=h*(T-Tinf) I can get T.
Then the temperature bc is a uniform fixedValue? Is that the way?

Actually everything is described in the wikipedia article. But if you'd like to see the math...

1. You've got heat transfer from fluid into a wall and them heat transfer inside the wall. And three temperatures T_fluid, T_wall (temperature of a wall near the fluid interface), T_inf (temperature of the wall at the infinity). You'd like to estimate hc of the system, so you can write q = hc*(T_fluid - T_inf)

2. 1/hc = 1/hc_flud->wall + 1/hc_wall. As we have got fluid, we can assume ideal contact with the wall, so 1/hc_fluid->wall will be equal 0 (so T_wall = T_fluid).

3. Heat transfer coefficient inside the wall can be estimated as kappa_wall/L_wall, where kappa_wall is heat conductivity of the wall material, L_wall is the thickness of the wall.

And finally if your wall thickness is zero, then yes, BC should be fixed value with T_inf.

[jherb]

I think you are looking for the turbulentHeatFluxTemperature boundary condition:
https://github.com/OpenFOAM/OpenFOAM...hScalarField.H