I was using fluent and i just started using OpenFoam. i have added Energy equation in pimpleFoam, i have compared fluent laminar result with it and it is working fine and soon i will compare turbulent result.
I am new in C++ language and OpenFoam so i have some kind of fear that may be i have modified pimpleFoam inaccurately and fortunately i am getting the correct result.
Here is my code, please check it and if there is any mistake please correct me, it will be very helpful for me.
mypimpleFoam.C
PHP Code:
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
laminarTransport.correct();
turbulence->correct();
}
}
#include "TEqn.H"
// Calculates dT/dy
Info<< "Calculating gradT" << endl;
label patchi = mesh.boundaryMesh().findPatchID("HEATEDWALL");
// calculate (dT/dy)_wall
//gradT.boundaryField()[patchi]=T.boundaryField()[patchi].snGrad();
gradT.boundaryFieldRef()[patchi]=T.boundaryFieldRef()[patchi].snGrad();
scalar area1 = gSum(mesh.magSf().boundaryField()[patchi]);
scalar sumField1 = 0;
if (area1 > 0)
{
sumField1 = gSum
(
mesh.magSf().boundaryField()[patchi]
* gradT.boundaryField()[patchi]
)/ area1;
}
Info<< " Average of gradT over patch 1 = " << sumField1 << endl;
//label patchi = mesh.boundaryMesh().findPatchID("myPatch");
//gradT.boundaryField()[patchi]=T.boundaryField()[patchi].snGrad();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
UEqn.H
PHP Code:
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}
PEqn.H
PHP Code:
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
MRF.makeRelative(phiHbyA);
adjustPhi(phiHbyA, U, p);
tmp<volScalarField> rAtU(rAU);
if (pimple.consistent())
{
rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
phiHbyA +=
fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU())*fvc::grad(p);
}
if (pimple.nCorrPISO() <= 1)
{
tUEqn.clear();
}
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAtU()*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
TEqn.H
PHP Code:
// Solve the ENRRGY (Temp) equation
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;
createFields.H
PHP Code:
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
);
//add Temprature field and
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
);
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
"nu",
dimViscosity,
transportProperties.lookup("nu")
);
dimensionedScalar kappa
(
transportProperties.lookup("kappa")
);
dimensionedScalar Pr
(
transportProperties.lookup("Pr")
);
dimensionedScalar Prt
(
transportProperties.lookup("Prt")
);
Info<< "Reading field gradT\n" << endl;
volScalarField gradT
(
IOobject
(
"gradT",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar
(
"gradT",
kappa.dimensions()/dimLength,
0
)
);
#include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, pimple.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
#include "createMRF.H"
please check it and if there is any mistake please correct me
thank you so much.
Srv