
[Sponsors] 
icoScalarTransportFoam + Buoyancy term (Boussinesq approximation) 

LinkBack  Thread Tools  Display Modes 
August 5, 2012, 17:41 
icoScalarTransportFoam + Buoyancy term (Boussinesq approximation)

#1 
New Member
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 9 
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) ); 

August 12, 2012, 16:19 

#2 
Senior Member

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? 

August 12, 2012, 21:03 

#3 
Senior Member
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 437
Rep Power: 15 
I implemented a turbulent Boussinesq buoyant SIMPLE solver. No pressure problems... did you use buoyantPressure boundary conditions where appropriate?
__________________
~~~ Follow me on twitter @DavidGaden 

August 12, 2012, 21:27 

#4  
Senior Member

Quote:
Quote:
Tnx for answer ~ 

August 15, 2012, 16:25 

#5 
New Member
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 9 
Hi Mojtaba
I do not have a way to validate the pressure values. For the "nobouyancy 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 

August 16, 2012, 03:57 

#6  
Senior Member

Quote:
Well i searched about this issue and I found this post: www.cfdonline.com/Forums/openfoam/74012thermalanalysisbuoyantboussinesqsimplefoam2.html#21 Quote:
http://www.opencae.jp/wiki/OpenFOAMVandVSIG/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 ~ 

August 21, 2012, 18:42 

#7 
New Member
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 9 
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*(TTref)*g. Why can't I use that term instead of (1beta(TTref))*g. Does this mean that FIDAP solves for p_rgh and OF1.6 solves for p? 

August 22, 2012, 08:36 

#8  
Senior Member

Quote:
http://www.cfdonline.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*(TTref)*g instead of (1beta(TTref))*g? 

August 22, 2012, 08:51 

#9 
New Member
Patrick
Join Date: Feb 2010
Posts: 10
Rep Power: 9 
You can find the equation in FIDAP theory Manual (in first few pages). I don't have the page number though.


Thread Tools  
Display Modes  


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 FluentUDF  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 ke model  Catherine Meissner  Phoenics  4  June 26, 2008 07:33 
Problem: buoyancy source term for KEEP model  Andrey V.Ivanov  Phoenics  0  May 23, 2003 03:50 