
[Sponsors] 
December 23, 2013, 09:56 
BuoyantBoussinesqPimpleFoam Modification for LES Capability

#1 
New Member
Simon Zhang
Join Date: Dec 2013
Posts: 12
Rep Power: 4 
Hi OpenFOAMers
I am currently simulating a turbulent wall jet and I want to investigate its mixing process using temperature as a tracer. Therefore I used BuoyantBoussinesqPimpleFoam solver, neglecting the buoyancy by set the value in the file /constant/g to (0 0 0) and I got reasonable velocity spreading results (see attachment). Then I changed the code for the capability of LES so as to get better results. However the LES results is quite abnormal that there seems to be no spreading and large velocity remains along the centerline (see attachment). It seems that there is no turbulence in LES simulation! I succeeded to compile the new solver named myBuoyantBoussinesqPimpleFoam and I changed the files in 0 directory by removing all the wall functions used for original BuoyantBoussinesqPimpleFoam solve. The mesh is checked ok. Since original solver can succeed but new solver cannot, I am afraid that there should be something wrong with the new solver. So I share all the solver code here and it will be greatly appreciated if someone can help me. creatFields.H Code:
Info<< "Reading thermophysical properties\n" << endl; Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading field p_rgh\n" << endl; volScalarField p_rgh ( IOobject ( "p_rgh", 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 ); #include "createPhi.H" #include "readTransportProperties.H" Info<< "Creating turbulence model\n" << endl; // autoPtr<incompressible::RASModel> turbulence // ( // incompressible::RASModel::New(U, phi, laminarTransport) // ); autoPtr<incompressible::turbulenceModel> turbulence //MJC 12212009 ( //MJC 12212009 incompressible::turbulenceModel::New //MJC 12212009 ( //MJC 12212009 U, //MJC 12212009 phi, //MJC 12212009 laminarTransport //MJC 12212009 ) //MJC 12212009 ); //MJC 12212009 // Kinematic density for buoyancy force volScalarField rhok ( IOobject ( "rhok", runTime.timeName(), mesh ), 1.0  beta*(T  TRef) ); // kinematic turbulent thermal thermal conductivity m2/s Info<< "Reading field kappat\n" << endl; volScalarField kappat ( IOobject ( "kappat", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("ghf", g & mesh.Cf()); volScalarField p ( IOobject ( "p", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), p_rgh + rhok*gh ); label pRefCell = 0; scalar pRefValue = 0.0; setRefCell ( p, p_rgh, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue ); if (p_rgh.needReference()) { p += dimensionedScalar ( "p", p.dimensions(), pRefValue  getRefCellValue(p, pRefCell) ); } Code:
#include "fvCFD.H" #include "singlePhaseTransportModel.H" //#include "RASModel.H" //add here #include "turbulenceModel.H" //stop here #include "pimpleControl.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" pimpleControl pimple(mesh); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; #include "readTimeControls.H" #include "CourantNo.H" #include "setDeltaT.H" //  Pressurevelocity PIMPLE corrector loop while (pimple.loop()) { #include "UEqn.H" #include "TEqn.H" //  Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence>correct(); } } runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } Code:
singlePhaseTransportModel laminarTransport(U, phi); // Thermal expansion coefficient [1/K] dimensionedScalar beta(laminarTransport.lookup("beta")); // Reference temperature [K] dimensionedScalar TRef(laminarTransport.lookup("TRef")); // Laminar Prandtl number dimensionedScalar Pr(laminarTransport.lookup("Pr")); // Turbulent Prandtl number dimensionedScalar Prt(laminarTransport.lookup("Prt")); Code:
{ volScalarField rAU("rAU", 1.0/UEqn.A()); surfaceScalarField rAUf("(1A(U))", fvc::interpolate(rAU)); U = rAU*UEqn.H(); phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, U, phi); surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); phi = buoyancyPhi; while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { // Calculate the conservative fluxes phi = p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure U = rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); } } #include "continuityErrs.H" p = p_rgh + rhok*gh; if (p_rgh.needReference()) { p += dimensionedScalar ( "p", p.dimensions(), pRefValue  getRefCellValue(p, pRefCell) ); p_rgh = p  rhok*gh; } } Code:
{ kappat = turbulence>nut()/Prt; kappat.correctBoundaryConditions(); volScalarField kappaEff("kappaEff", turbulence>nu()/Pr + kappat); fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T)  fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(); rhok = 1.0  beta*(T  TRef); } Code:
// Solve the momentum equation fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence>divDevReff(U) ); UEqn.relax(); if (pimple.momentumPredictor()) { solve ( UEqn == fvc::reconstruct ( (  ghf*fvc::snGrad(rhok)  fvc::snGrad(p_rgh) )*mesh.magSf() ) ); } files Code:
myBuoyantBoussinesqPimpleFoam.C EXE = $(FOAM_USER_APPBIN)/myBuoyantBoussinesqPimpleFoam Code:
EXE_INC = \ I../buoyantBoussinesqSimpleFoam \ I$(LIB_SRC)/finiteVolume/lnInclude \ I$(LIB_SRC)/turbulenceModels \ I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel/lnInclude \ I$(LIB_SRC)/transportModels \ I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel EXE_LIBS = \ lfiniteVolume \ lmeshTools \ lincompressibleTurbulenceModel \ lincompressibleRASModels \ lincompressibleLESModels \ lincompressibleTransportModels Simon 

December 24, 2013, 22:17 

#2 
Member
sqing
Join Date: Sep 2012
Location: Dalian
Posts: 77
Rep Power: 5 
Hi Simon,
Is your case a 2d case or 3d for RANS/LES? And you said that you have gotten reasonable velocity result with RANS. Is it concluded by comparing it with exp data or just coming out with the outlook? what's more can you share your LES case(0 and system) with us? Best regards Sunxing 

December 25, 2013, 01:21 

#3  
New Member
Simon Zhang
Join Date: Dec 2013
Posts: 12
Rep Power: 4 
Quote:
Thanks for your reply. The case is a 3d one and the OF version is 2.1.1. Actually I have performed RANS/LES simulations of wall jets in FLUENT in the past and it was concluded that LES was better than RANS. This is the first time for me to use OF and I want to test the solver in OF by simulating the same case. The solver code posted is created by following guide named buoyantBoussinesqPisoFoam on the openfoamwiki. The computational domain is a block(up: free water surface, down: bottom wall, lelf/right/back: side wall, hole: inlet encircled by back surface, exit: far field outlet). The system files are as following: controlDict Code:
application myBuoyantBoussinesqPimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 30; deltaT 0.0004; writeControl timeStep; writeInterval 2500; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; adjustTimeStep no; // maxCo 0.5; functions { fieldAverage1 { type fieldAverage; functionObjectLibs ( "libfieldFunctionObjects.so" ); enabled true; outputControl outputTime; fields ( U { mean on; prime2Mean on; base time; } T { mean on; prime2Mean on; base time; } ); } } Code:
ddtSchemes { default backward; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) Gauss upwind; div(phi,T) Gauss upwind; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(phi,R) Gauss upwind; div(R) Gauss linear; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear uncorrected; laplacian((1A(U)),p_rgh) Gauss linear uncorrected; laplacian(kappaEff,T) Gauss linear uncorrected; laplacian(DkEff,k) Gauss linear uncorrected; laplacian(DepsilonEff,epsilon) Gauss linear uncorrected; laplacian(DREff,R) Gauss linear uncorrected; } interpolationSchemes { default linear; } snGradSchemes { default uncorrected; } fluxRequired { default no; p_rgh; } // ************************************************************************* // Code:
solvers { p_rgh { solver PCG; preconditioner DIC; tolerance 1e8; relTol 0.01; } p_rghFinal { $p_rgh; relTol 0; } "(UTkepsilonR)" { solver PBiCG; preconditioner DILU; tolerance 1e6; relTol 0.1; } "(UTkepsilonR)Final" { $U; relTol 0; } } PIMPLE { momentumPredictor no; nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } relaxationFactors { fields { } equations { "(UTkepsilonR)" 1; "(UTkepsilonR)Final" 1; } } // ************************************************************************* // transportProperties Code:
transportModel Newtonian; // Laminar viscosity nu nu [ 0 2 1 0 0 0 0 ] 1e06; // Thermal expansion coefficient beta beta [0 0 0 1 0 0 0] 2.7e04; // Reference temperature TRef TRef [0 0 0 1 0 0 0] 300; // Laminar Prandtl number Pr Pr [0 0 0 0 0 0 0] 7; // Turbulent Prandtl number Prt Prt [0 0 0 0 0 0 0] 0.85; Code:
LESModel oneEqEddy; delta cubeRootVol; printCoeffs on; cubeRootVolCoeffs { deltaCoeff 1; } PrandtlCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } Cdelta 0.158; } vanDriestCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } Aplus 26; Cdelta 0.158; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } Code:
simulationType LESModel; p_rgh Code:
dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { left { type buoyantPressure; rho rhok; value uniform 0; } right { type buoyantPressure; rho rhok; value uniform 0; } down { type buoyantPressure; rho rhok; value uniform 0; } up { type symmetryPlane; } hole { type buoyantPressure; rho rhok; value uniform 0; } back { type buoyantPressure; rho rhok; value uniform 0; } exit { type fixedValue; value uniform 0; } } Code:
dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { left { type fixedValue; value uniform (0 0 0); } right { type fixedValue; value uniform (0 0 0); } down { type fixedValue; value uniform (0 0 0); } up { type symmetryPlane; } hole { type fixedValue; value uniform (2.2288 0 0); } back { type fixedValue; value uniform (0 0 0); } exit { type zeroGradient; } } Code:
dimensions [0 0 0 1 0 0 0]; internalField uniform 300; boundaryField { left { type fixedValue; value uniform 300; } right { type fixedValue; value uniform 300; } down { type zeroGradient; } up { type symmetryPlane; } hole { type fixedValue; value uniform 305; } back { type zeroGradient; } exit { type zeroGradient; } } Code:
dimensions [0 2 1 0 0 0 0]; internalField uniform 0; boundaryField { left { type fixedValue; value uniform 0; } right { type fixedValue; value uniform 0; } down { type fixedValue; value uniform 0; } up { type symmetryPlane; } hole { type fixedValue; value uniform 1.18e5; } back { type fixedValue; value uniform 0; } exit { type zeroGradient; } } Code:
dimensions [0 2 1 0 0 0 0]; internalField uniform 0; boundaryField { left { type zeroGradient; } right { type zeroGradient; } down { type zeroGradient; } up { type symmetryPlane; } hole { type zeroGradient; } back { type zeroGradient; } exit { type zeroGradient; } } Code:
dimensions [0 2 2 0 0 0 0]; internalField uniform 0.018628; boundaryField { left { type fixedValue; value uniform 0; } right { type fixedValue; value uniform 0; } down { type fixedValue; value uniform 0; } up { type symmetryPlane; } hole { type fixedValue; value uniform 0.018628; } back { type fixedValue; value uniform 0; } exit { type zeroGradient; } } 

December 25, 2013, 02:13 

#4 
Member
sqing
Join Date: Sep 2012
Location: Dalian
Posts: 77
Rep Power: 5 
Hi simon,
How many points in your LES case? Does it satisfy y+<1 in the first level grid along the wall side? Best, sunxing 

December 25, 2013, 05:35 

#5  
New Member
Simon Zhang
Join Date: Dec 2013
Posts: 12
Rep Power: 4 
Quote:
I did not do y+ check and the grid is very coarse in my simulation. However it is the same grid in RANS simulation. This is just a trial to test solver and I will try a finer grid for LES in the future. In your view, what effects will be caused from the large y+ along wall side? As seen in the figure, I do not think it is due to the large y+ because we even do not see any turbulence, like vortex and so on. It just looks like a laminar flow rather than a flow from LES. Regards, Simon 

December 28, 2013, 09:05 

#6 
Member
Peter
Join Date: Nov 2011
Posts: 45
Rep Power: 6 
Hi, Simon!
Are these all of the files in your 0 dir? I am recently trying to use buoyantBoussinesqPimpleFoam for LES too. I am not familiar with this solver, so I just try to mix a buoyantBoussinesqPimpleFoam tutorial case and an LES tutorial case together and turn it into my case. In my 0 dir, I have alphat, B, nuTilda and T.org files besides those you showed here. Now I know that nuTilda is not necessary in LES case and B is useless for dynSmagorinsky model which I use. I also find out that T.org is a backup of T, so the content of T.org can be the same as T (please correct me if I am wrong.) My question is: is alphat not needed at all? Regards, Peter Last edited by palmerlee; December 29, 2013 at 09:52. 

December 29, 2013, 09:59 

#7  
Member
Peter
Join Date: Nov 2011
Posts: 45
Rep Power: 6 
Also, how to set up k for LES? I know that for RANS, k can be calculated using equation k=1.5*(u_avg*I)^2. Is it the same for LES?
Quote:


December 29, 2013, 11:01 

#8  
New Member
Simon Zhang
Join Date: Dec 2013
Posts: 12
Rep Power: 4 
Quote:
Are you using OF version 2.2.x? Because alphat is needed in this version. I am using OF version 2.1.1 and kappat is required instead of alphat. For k, I think the estimated k is independent of what solver used. So I estimated k by the same equation described in chapter 2 in the tutorial. Regards, Simon 

December 29, 2013, 21:47 

#9  
Member
Peter
Join Date: Nov 2011
Posts: 45
Rep Power: 6 
Hi, Simon!
Thank you for your reply! By "chapter 2 in the tutorial", do you mean User Guide, where I didn't find the equation? I agree with you on that k is independent of what solver used. But it is not the same thing in RANS and in LES if I understand correctly. In LES, k is subgridscale kinetic energy. So I guess it is estimated by a different equation from that in RANS. That's only my option. Also, I came up with another question. How to calculate kappat on inlet boundary? Quote:
Peter 

December 29, 2013, 22:54 

#10  
Member
Peter
Join Date: Nov 2011
Posts: 45
Rep Power: 6 
Quote:


December 30, 2013, 01:55 

#11  
New Member
Simon Zhang
Join Date: Dec 2013
Posts: 12
Rep Power: 4 
Quote:
Thanks for your comments. Now I have changed the fvscheme and let's see what will happen. For k, I used the equations in section 2.1.8.1 of official user guide. In my simulation, the inlet k is believed not to affect the flow field much and thus I used the same equation as in RANS. For kappat, I have to say sorry because I do not understand it very clearly either. Please see thread http://www.cfdonline.com/Forums/ope...amkappat.html for help. Regards, Simon 

Tags 
les, modification, solver 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Modification of a mesh during runtime  Scofield  OpenFOAM Meshing & Mesh Conversion  6  October 15, 2013 12:49 
Runtime modification of a volScalarField  Scofield  OpenFOAM Programming & Development  4  October 1, 2013 08:31 
Fluent is not accepting Modification in Material properties  er.mkumar  FLUENT  0  September 1, 2012 14:24 
(Dynamic) Mesh Modification  Mr. Joe  Fluent UDF and Scheme Programming  3  March 24, 2010 21:31 
compressible modification of nearwall turbulence  Quain Tchew  Main CFD Forum  0  March 4, 2002 02:29 