CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

icoScalarTransportFoam + Buoyancy term (Boussinesq approximation)

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 5, 2012, 17:41
Default icoScalarTransportFoam + Buoyancy term (Boussinesq approximation)
  #1
New Member
 
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 7
will_avoid_comm_solvers is on a distinguished road
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)
);




will_avoid_comm_solvers is offline   Reply With Quote

Old   August 12, 2012, 16:19
Default
  #2
Senior Member
 
Mojtaba.a's Avatar
 
Mojtaba Amiraslanpour
Join Date: Jun 2011
Location: Zanjan, Iran
Posts: 233
Rep Power: 7
Mojtaba.a is on a distinguished road
Send a message via Yahoo to 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?
Mojtaba.a is offline   Reply With Quote

Old   August 12, 2012, 21:03
Default
  #3
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
I implemented a turbulent Boussinesq buoyant SIMPLE solver. No pressure problems... did you use buoyantPressure boundary conditions where appropriate?
__________________
~~~
Follow me on twitter @DavidGaden
marupio is offline   Reply With Quote

Old   August 12, 2012, 21:27
Default
  #4
Senior Member
 
Mojtaba.a's Avatar
 
Mojtaba Amiraslanpour
Join Date: Jun 2011
Location: Zanjan, Iran
Posts: 233
Rep Power: 7
Mojtaba.a is on a distinguished road
Send a message via Yahoo to Mojtaba.a
Quote:
Originally Posted by marupio View Post
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.

Quote:
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 ~
Mojtaba.a is offline   Reply With Quote

Old   August 15, 2012, 16:25
Default
  #5
New Member
 
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 7
will_avoid_comm_solvers is on a distinguished road
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
will_avoid_comm_solvers is offline   Reply With Quote

Old   August 16, 2012, 03:57
Default
  #6
Senior Member
 
Mojtaba.a's Avatar
 
Mojtaba Amiraslanpour
Join Date: Jun 2011
Location: Zanjan, Iran
Posts: 233
Rep Power: 7
Mojtaba.a is on a distinguished road
Send a message via Yahoo to Mojtaba.a
Quote:
Originally Posted by will_avoid_comm_solvers View Post
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
Quote:
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 ~
Mojtaba.a is offline   Reply With Quote

Old   August 21, 2012, 18:42
Default
  #7
New Member
 
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 7
will_avoid_comm_solvers is on a distinguished road
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?
will_avoid_comm_solvers is offline   Reply With Quote

Old   August 22, 2012, 08:36
Default
  #8
Senior Member
 
Mojtaba.a's Avatar
 
Mojtaba Amiraslanpour
Join Date: Jun 2011
Location: Zanjan, Iran
Posts: 233
Rep Power: 7
Mojtaba.a is on a distinguished road
Send a message via Yahoo to Mojtaba.a
Quote:
Originally Posted by will_avoid_comm_solvers View Post
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?
Mojtaba.a is offline   Reply With Quote

Old   August 22, 2012, 08:51
Default
  #9
New Member
 
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 7
will_avoid_comm_solvers is on a distinguished road
You can find the equation in FIDAP theory Manual (in first few pages). I don't have the page number though.
will_avoid_comm_solvers is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
buoyancy production term in turbulence model braennstroem OpenFOAM Running, Solving & CFD 5 January 13, 2012 06:43
How to realize Boussinesq approximation (buoyancy) in Fluent-UDF dochanxu FLUENT 0 October 14, 2010 14:42
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 15:45
buoyancy source term in k-e model Catherine Meissner Phoenics 4 June 26, 2008 07:33
Problem: buoyancy source term for KE-EP model Andrey V.Ivanov Phoenics 0 May 23, 2003 03:50


All times are GMT -4. The time now is 02:17.