[will_avoid_comm_solvers]

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

[Mojtaba.a]

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?

[marupio]

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

[Mojtaba.a]

Originally Posted by marupio

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.

wall
{
type buoyantPressure;
rho rhok;
value uniform 0;
}

air-inlet
{
type zeroGradient;
}

out_flow
{
type zeroGradient;
}

symmetry_plane
{
type symmetryPlane;
}

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 ~

[will_avoid_comm_solvers]

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

[Mojtaba.a]

Originally Posted by will_avoid_comm_solvers

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

The solvers are applicable to any case where the Boussinesq approximation can be justified, they just do not work very well on poor meshes. So if you have a perfect block structured mesh, then the solvers work very well (in our experience), but if you have tets, polyhedra or deformed hexes, then the results are completely unphysical. It might just be that we did not use the right combination of differencing schemes, but believe me we tried many.
Also, always use 2nd order or near 2nd order schemes like linearUpwind if you are trying to match experiment - you will be disappointed otherwise.

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 ~

[will_avoid_comm_solvers]

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?

[Mojtaba.a]

Originally Posted by will_avoid_comm_solvers

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?

[will_avoid_comm_solvers]

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