icoScalarTransportFoam + Buoyancy term (Boussinesq approximation)
 

Dear All

I am trying to compile a new solver for simulating unsteady species transport with buoyancy force term (Boussinesq approximation).

I started from the scalarTransportFoam and used icoFOAM files (thanks to wikipage and other online resources) to create the unsteady species transport equation (but no buoyancy term). I was able to compile and run this code. The results were also checked experimentally.

Now, I am hoping to include the buoyancy term in the icoScalarTransportFoam. For this, I used the information from BuoyantBoussinesqPisoFoam. I am pasting the .c and the createFields file for the new solver below. It would be great if anyone can have a look at them and point out to any potential issues. So far, I was able to compile the files without any issue

icoScalarGravityTransportFoam.c file

Application
icoScalarGravityTransportFoam

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

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

#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"


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

Info<< "\nStarting time loop\n" << endl;

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

#include "readTimeControls.H"
#include "readPISOControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"


fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

UEqn.relax();

if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(rhok)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
)
);
}

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
# include "TEqn.H"
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();
}

runTime.write();

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

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

return 0;
}


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





createFields.H file


Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

dimensionedScalar nu
(
transportProperties.lookup("nu")
);
dimensionedScalar beta(transportProperties.lookup("beta"));

dimensionedScalar TRef(transportProperties.lookup("TRef"));


dimensionedScalar Diffu
(
transportProperties.lookup("Diffu")
);
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
);

# include "createPhi.H"


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

// Kinematic density for buoyancy force
volScalarField rhok
(
IOobject
(
"rhok",
runTime.timeName(),
mesh
),
1.0 - beta*(T - TRef)
);

Hi Patrick. I'm very interested in this case while i'm trying to do a similar task. I was trying to use buoyantBoussinesqSimpleFoam in which I have implemented 2 scalar transport equations. everything seems ok except pressure. My pressure varies a lot and doesn't have any agreement with fluent and experimental results. Do u have such a problem on pressure?
Your solver is a laminar solver. Have you ever tried to make a turbulence solver?

I implemented a turbulent Boussinesq buoyant SIMPLE solver. No pressure problems... did you use buoyantPressure boundary conditions where appropriate?

Well in p_rgh file I used zeroGradient for inlet and outlet and for all other patches I used buoyantPressure.

Are they true? My p values are about -17 to 12 pa, but the correct result is something like 1 to 2 pa.
Tnx for answer ~

Hi Mojtaba

I do not have a way to validate the pressure values. For the "no-bouyancy code" I was able to compare the species concentration with experiments. The values match well. I am looking for a way to verify this code (with buoyancy term). Any suggestions will be very useful

Hi Patrick
Well i searched about this issue and I found this post:
www.cfd-online.com/Forums/openfoam/74012-thermal-analysis-buoyantboussinesqsimplefoam-2.html#21
For experimenting this I tried to validate my results based on a tutorial made by Masa, which is airConditionedRoom, here:
http://www.opencae.jp/wiki/OpenFOAM-VandV-SIG/Misc
First I used a tetrahedra mesh to simulate the exact case of Masa, But surprisingly buoyantBoussinesqSimpleFoam gave wrong results. So I changed the grid into a simple structured gird and the results were OK.
I also found out that the pressure which we are searching for is p_rgh not p, due to this fact that the momentum equation is solved for p_rgh while buoyancy term is out of the derivatives.
Regards ~

Thanks Mojtaba. That was helpful. I have OpenFoam 1.6. In OpenFOAM1.6, are we solving for p or p - rgh?

Also, I noticed that in FIDAP, the buoyancy term is just -beta*(T-Tref)*g. Why can't I use that term instead of (1-beta(T-Tref))*g. Does this mean that FIDAP solves for p_rgh and OF1.6 solves for p?

Well these are my questions too, I'm trying to understand it in this thread:

http://www.cfd-online.com/Forums/ope...tml#post378166

I will be very happy if you can join us Patrick.
also I am very interested to know how FIDAP solves for boussinesq approximation. because I am validating a OF case with FIDAP. can you tell me where you found out that FIDAP solves -beta*(T-Tref)*g instead of (1-beta(T-Tref))*g?

You can find the equation in FIDAP theory Manual (in first few pages). I don't have the page number though.