[hawkeye321]

Hey Foamers

I have added energy equation to icoFoam, and I am solving natural convection in cavity. The problem is that as I am using smaller and smaller time steps, the velocities become smaller and smaller, and I can not get a results which is independent of time step.

So do you have any ideas to help me getting time step independent results?

[niaz]

dear mahdi
do you add temperature truly?
it seams that your solver has a problem.
you should attach your solver. it is hard to make decision without looking to your solver

[hawkeye321]

Dear Niyaz

this is my solver

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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
icoFoam

Description
Transient solver for incompressible, laminar flow of Newtonian fluids.

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

#include "fvCFD.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"

fvVectorMatrix UEqn
(
(Density)*fvm::ddt(U)
+ (Density)*fvm::div(phi, U)
- (Visc)*(fvm::laplacian(U))
);

solve(UEqn == GravityVector*g*Density*(BetaT*(T-TNot)) - (fvc::grad(p)));
//solve(UEqn == -fvc::grad(p));
// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

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

#include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
fvScalarMatrix TEqn
(
Density*Cp*fvm::ddt(T)
+ Density*Cp*fvm::div(phi, T)
- Cond*(fvm::laplacian(T))
);
TEqn.solve();

//F = Gravity*(1-Beta*(T-TemAmb));
///////////
fvScalarMatrix CEqn
(
fvm::ddt(C)
+ fvm::div(phi, C)
);
CEqn.solve();
///////////

// forAll( h, celli)
// {
// T[celli] = (h[celli] - 1)/Ste.value();
//T[celli] = (h[celli] - L.value())/Cp.value();
// h[celli] = 1 + (Ste.value()*T[celli]);
// h[celli] = (Cp.value()*T[celli]) + L.value();

// }
///////////
//F = (GravityVector*Pr*T)-(GravityVector*Le*C);
//F = gl*GravityVector*g*Density*(1-BetaT*(T-TNot));
F = GravityVector*g*Density*(BetaT*(T-TNot));
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************** *********************** //

[niaz]

your changes are strange.
for U:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
==
-beta*(T - T0)*g
);
for T:
solve
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);

where nu is viscosity and DT is conductivity
your equation seems have a problem.

[hawkeye321]

Hi niaz
Thanks for your comments. I tried with the new equations, but still have problem with dt. any more comments?

[niaz]

solve(UEqn == GravityVector*g*Density*(BetaT*(T-TNot)) - (fvc::grad(p)));
I think that before the GravityVector, you should have (-)

[hawkeye321]

Thank you. But it will not fix the problem. It will only change the direction of U and V.

[niaz]

Can you show a contour of two results?
your case is a bit strange.

[niaz]

at first, I suggest you to solve a forced convection to prove the exactness of the solver then add boussinesqe term to equation.

[hawkeye321]

I fixed the problem.
The key is
in order to get a proper behavior of velocity-pressure coupling in PISO algorithm, I should include the source term inside the velocity equation such that

fvVectorMatrix UEqn
(
fvm::ddt(U)
+fvm::div(phi,U)
-fvm::laplacian(ModifVisc,U)
==
GravityVector*betaT*(T-Tnot)
);
solve(UEqn == -fvc::grad(p) );