
[Sponsors] 
kappatJayatillekeWallFunctionFvPatchScalarField changes between OpenFOAM 171 and 201 

LinkBack  Thread Tools  Display Modes 
October 13, 2011, 16:34 
kappatJayatillekeWallFunctionFvPatchScalarField changes between OpenFOAM 171 and 201

#1 
Senior Member

Hi everyone!
I discovered following while simulating buoyant flow with 1.7.0 and 1.7.1 versions (Ubuntu 10.04 LTS packages from apt): 1. Comparable to 1.7.0 in version 1.7.1 kappat field is introduced in TEqn for incompressible SIMPLE and PIMPLE solvers. v1.7.0: Code:
{ volScalarField kappaEff ( "kappaEff", turbulence>nu()/Pr + turbulence>nut()/Prt ); fvScalarMatrix TEqn ( fvm::ddt(T) + fvm::div(phi, T)  fvm::laplacian(kappaEff, T) ); TEqn.relax(); TEqn.solve(mesh.solver(T.select(finalIter))); rhok = 1.0  beta*(T  TRef); } 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(mesh.solver(T.select(finalIter))); rhok = 1.0  beta*(T  TRef); } diff for kappatJayatillekeWallFunctionFvPatchScalarField.C between 1.7.1 and 2.0.1 results in Code:
210c210 < const scalarField& nuw = rasModel.nu().boundaryField()[patchI];  > const scalarField& nuw = rasModel.nu()().boundaryField()[patchI]; 238c238 < scalar kt = nu*(yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P)  1/Pr);  > scalar kt = nu*(yPlus/(Prt_*(log(E_*yPlus)/kappa_ + P))  1/Pr); 264c264,268 < makePatchTypeField(fvPatchScalarField, kappatJayatillekeWallFunctionFvPatchScalarField);  > makePatchTypeField > ( > fvPatchScalarField, > kappatJayatillekeWallFunctionFvPatchScalarField > ); log(E_*yPlusTherm) where yPlusTherm can be equal zero. You cal also notice that parenthesis around (log(E_*yPlus)/kappa_ + P)) result in different expression than Prt_/kappa_*log(E_*yPlusTherm) + P. Hope it would be helpful for someone)))) Cheers, A.
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at Last edited by makaveli_lcf; October 13, 2011 at 16:46. Reason: Parenthesis aroun expression. 

October 17, 2011, 11:52 

#2 
New Member
Howard NJOKU
Join Date: Nov 2010
Location: Nsukka, Nigeria
Posts: 9
Rep Power: 8 
Please can you help me with an explanation of what rhok is?


October 17, 2011, 12:34 

#3 
Senior Member

__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

October 19, 2011, 10:06 

#4 
New Member
Howard NJOKU
Join Date: Nov 2010
Location: Nsukka, Nigeria
Posts: 9
Rep Power: 8 
Thanks for the prompt response. I know about the boussinesque approximation. My question is because of a case I'm setting up in OpenFOAM. I am simulating stratification in a cylindrical tank. I have modified the 'hotroom' example under the bouyantBoussinesqPimpleFoam solver and changed the entries in the 'transportProperties' file to those of water. When the run the case with the original entries (which I suspect to be for air), I get realistic results  the stratification is obvious. But when I change the properties  Pr, beta, nu, to those of water the results when viewed show no signs of stratification. Could you be of help please?


October 19, 2011, 15:13 

#5 
Senior Member

Hello Howard!
It is always better when you ask detailed questions. Because for the obvious ones you can find answers yourself and it is more useful for you to learn finding answers. Would you please post you pictures which do not satisfy you and also the parameters which you used? Cheers, A.
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

October 19, 2011, 15:41 

#6 
Senior Member

Do you mean simulating compressible flow? Of couse for water (almost incompressible liquid) there is no such stratification as for e.g. air.
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

October 24, 2011, 16:58 

#7 
New Member
Howard NJOKU
Join Date: Nov 2010
Location: Nsukka, Nigeria
Posts: 9
Rep Power: 8 
Hello Dr. Alexander,
Please pardon the delay in my response. I travelled out of town without my laptop and was without access to computer facilities. The changes I made in the transportProperties file were nu : changed from 1e05 to 4.5e07 (m2sec1) beta: changed from 3e03 to 5e04 (K1) Pr: changed from 0.9 to 3.0 Prt: changed from 0.7 to 2.3. Tref was left unchanged. The RASModel keyword in the RASProperities file was changed to laminar. In the 0/p file the wall pressures were set with the $internalField entry, while the p_rgh value was set to constant 0 in the 0/p_rgh file. You'll see from the attached picture that there is no evidence of stratification. Stratification is obtainable in incompressible flows (the 'buoyantBoussinesqPimpleFoam' solver is for "Transient solver for buoyant, turbulent flow of incompressible fluids" according to the user guide.) This is the basis for stratified hot water storages which is what I'm trying to study. Howard 

January 14, 2012, 11:27 
kappatw instead of kappat?

#8 
New Member
chirag patel
Join Date: Feb 2011
Posts: 1
Rep Power: 0 
Dear Dr.Alexander,
In the file "kappatJayatillekeWallFunctionFvPatchScalarField.C" kappatw shouldn't be kappat? I searched much, but didn't find any use of kappatw. With Regards, Chirag Patel 

January 31, 2012, 15:48 

#9 
Senior Member
Aram Amouzandeh
Join Date: Mar 2009
Location: Vienna, Vienna, Austria
Posts: 186
Rep Power: 10 
dear all,
i wonder why the JayatillekeWallFunction is not available for the the calculation of alphat in the compressible RANS version. for now, in the compressible version alphat is evaluated as mut/Prt at the wall. i m thinking of implementing the procedure of the JayatillekeWallFunction for compressible RANS in alphatWallFunctions. does it make sens or is alphatw = mutw/Prt resonable anyway?? @Chirag Patel kappatw is taken from the expression Code:
scalarField& kappatw = *this; cheers, aram 

February 1, 2012, 04:25 

#10 
Senior Member
Aram Amouzandeh
Join Date: Mar 2009
Location: Vienna, Vienna, Austria
Posts: 186
Rep Power: 10 
hi,
I though about something like this: Code:
// Populate boundary values scalarField& alphatw = *this; forAll(alphatw, faceI) { [ ... ] if (yPlus > yPlusTherm) { scalar mu = muw[faceI]; scalar at = mu*(yPlus/(Prt_*(log(E_*yPlus))/kappa_ + P)  1/Pr); alphatw[faceI] = max(0.0, at); } else { alphatw[faceI] = 0.0; } cheers, Aram 

February 1, 2012, 05:53 

#11 
Senior Member

hi Aram!
you are right) But I found an implementation in version OF 2.10: see http://foam.sourceforge.net/docs/cpp...ce.html#l00195 Code:
00195 void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs() 00196 { 00197 if (updated()) 00198 { 00199 return; 00200 } 00201 00202 const label patchI = patch().index(); 00203 00204 const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties"); 00205 00206 const scalar Cmu25 = pow025(Cmu_); 00207 00208 const scalarField& y = rasModel.y()[patchI]; 00209 00210 const scalarField& muw = rasModel.mu().boundaryField()[patchI]; 00211 00212 const scalarField& alphaw = rasModel.alpha().boundaryField()[patchI]; 00213 scalarField& alphatw = *this; 00214 00215 const tmp<volScalarField> tk = rasModel.k(); 00216 const volScalarField& k = tk(); 00217 00218 const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; 00219 const scalarField magUp(mag(Uw.patchInternalField()  Uw)); 00220 const scalarField magGradUw(mag(Uw.snGrad())); 00221 00222 const scalarField& rhow = rasModel.rho().boundaryField()[patchI]; 00223 const fvPatchScalarField& hw = 00224 rasModel.thermo().h().boundaryField()[patchI]; 00225 00226 // Heat flux [W/m2]  lagging alphatw 00227 const scalarField qDot((alphaw + alphatw)*hw.snGrad()); 00228 00229 // Populate boundary values 00230 forAll(alphatw, faceI) 00231 { 00232 label faceCellI = patch().faceCells()[faceI]; 00233 00234 scalar uTau = Cmu25*sqrt(k[faceCellI]); 00235 00236 scalar yPlus = uTau*y[faceI]/(muw[faceI]/rhow[faceI]); 00237 00238 // Molecular Prandtl number 00239 scalar Pr = muw[faceI]/alphaw[faceI]; 00240 00241 // Moleculartoturbulent Prandtl number ratio 00242 scalar Prat = Pr/Prt_; 00243 00244 // Thermal sublayer thickness 00245 scalar P = Psmooth(Prat); 00246 scalar yPlusTherm = this>yPlusTherm(P, Prat); 00247 00248 // Evaluate new effective thermal diffusivity 00249 scalar alphaEff = 0.0; 00250 if (yPlus < yPlusTherm) 00251 { 00252 scalar A = qDot[faceI]*rhow[faceI]*uTau*y[faceI]; 00253 scalar B = qDot[faceI]*Pr*yPlus; 00254 scalar C = Pr*0.5*rhow[faceI]*uTau*sqr(magUp[faceI]); 00255 alphaEff = A/(B + C + VSMALL); 00256 } 00257 else 00258 { 00259 scalar A = qDot[faceI]*rhow[faceI]*uTau*y[faceI]; 00260 scalar B = qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus) + P); 00261 scalar magUc = uTau/kappa_*log(E_*yPlusTherm)  mag(Uw[faceI]); 00262 scalar C = 00263 0.5*rhow[faceI]*uTau 00264 *(Prt_*sqr(magUp[faceI]) + (Pr  Prt_)*sqr(magUc)); 00265 alphaEff = A/(B + C + VSMALL); 00266 } 00267 00268 // Update turbulent thermal diffusivity 00269 alphatw[faceI] = max(0.0, alphaEff  alphaw[faceI]); 00270 00271 if (debug) 00272 { 00273 Info<< " uTau = " << uTau << nl 00274 << " Pr = " << Pr << nl 00275 << " Prt = " << Prt_ << nl 00276 << " qDot = " << qDot[faceI] << nl 00277 << " yPlus = " << yPlus << nl 00278 << " yPlusTherm = " << yPlusTherm << nl 00279 << " alphaEff = " << alphaEff << nl 00280 << " alphaw = " << alphaw[faceI] << nl 00281 << " alphatw = " << alphatw[faceI] << nl 00282 << endl; 00283 } 00284 } 00285 00286 fixedValueFvPatchField<scalar>::updateCoeffs(); 00287 } 00288 Cheers, Alex
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at Last edited by makaveli_lcf; February 1, 2012 at 06:00. Reason: Doxygen link fix 

February 1, 2012, 09:26 

#12 
Senior Member
Aram Amouzandeh
Join Date: Mar 2009
Location: Vienna, Vienna, Austria
Posts: 186
Rep Power: 10 
hi Alex,
thanks for the hint (still using OF 1.7). I had a look on the implementation of OF 2.1. Three main differences in the why how i proposed it: 1) alphat isnt set to zero for yPlus < yPlusLam, but I couldn t find a reference so far were the equation comes from 2) The BC calculates alphaEff and updates alphat by subtracting alpha (laminar) 3) For yPlus > yPlusLam the calculation of alphaEff differs from the version I proposed by an additional term in the denominator which contains velocities (which I couldn t find out yet where it comes from): Code:
alphaEff = mu*(yPlus/(Pr_t*log(E *yPlus)/kappa_ + 0.5*rhow*uTau*(Pr_t*sqr(mag(Up)) + (Pr_t  Pr)*sqr(mag(Uc))))) cheers, Aram 

February 1, 2012, 09:34 

#13 
Senior Member

Hi,
those additional terms are for compressible flow, thus there are none of them in incompressible BC version. Check for the full version: http://my.fit.edu/itresources/manual...pdf/flu_th.pdf, (page 124) Sorry for the wrong link to the Theory Guide) Cheers!
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at Last edited by makaveli_lcf; February 1, 2012 at 11:36. Reason: Theory Guide link 

February 1, 2012, 11:44 

#14 
Senior Member
Aram Amouzandeh
Join Date: Mar 2009
Location: Vienna, Vienna, Austria
Posts: 186
Rep Power: 10 
That s great! Got it now!!
Thanks again!! Cheers, Aram 

July 29, 2012, 02:01 
help me plz

#15 
Member
vahid
Join Date: Feb 2012
Location: MashhadIran
Posts: 80
Rep Power: 6 
Hi every foamers.
I want to add k(kinetic turbulence energy) in one solver. for this work, added next line in the code of solver: const volScalarField &k=U_.db().lookupObject<volScalarField>("k") and wmake was Successfully . but How can I Understand that this k is the same with k(kinetic turbulence energy)???Is it true?????or not???? Because I replase another word in this : const volScalarField &M=U_.db().lookupObject<volScalarField>("M") but nothing error was not occured!!!!!!?????? 

July 29, 2012, 04:33 

#16 
Senior Member

Could the lookup method just give you a null pointer, then no error is detected
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

July 29, 2012, 04:46 
Thanks

#17 
Member
vahid
Join Date: Feb 2012
Location: MashhadIran
Posts: 80
Rep Power: 6 
Thanks Dr. Alexander A. Vakhrushev for your answer.
its means this work not true???this k is not kinetic turbulence energy??? if your answer yes!!!Please suggest me one way to Correction this code with k(kinetic turbulence energy)!!! 

July 29, 2012, 04:51 

#18 
Senior Member

Could you show more information, e.g.post your solver (if it is not secret) to understand your problem. Actually it should be ok what yo use to get k field.
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

July 29, 2012, 05:28 
my solver is in Appendix

#19 
Member
vahid
Join Date: Feb 2012
Location: MashhadIran
Posts: 80
Rep Power: 6 
Thanks for your answer Dr.
it is not secret!!! my solver is in Appendix.I want to added kinetic turbulence energy in Directory: OpenFOAM/vahid 2.0.1/run/myinterPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.c in file:SchnerrSauer.c I have add (turbulence kinetic Energy or k ) in a part of this solver. first: added: I$(LIB_SRC)/turbulenceModels \ I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ in option. and then: added //.................................................. .....changed // Construct incompressible turbulence model autoPtr<incompressible::RASModel> turbulence ( incompressible::RASModel::New(U, phi, twoPhaseProperties()) ); //.................................................. .....changed in creatField. and then: added #include "RASModel.H" and: //.................................................. ............. volScalarField turbKinEnergy = turbulence().k(); //turbKinEnergy = turbulence().k(); volScalarField turbDisEnergy = turbulence().epsilon(); in this and turbkinEnergy=k //.................................................. ............. in myinterPhaseChangeFoam.c and wmake is ok!!! but I want in this solver in directory: /run/myinterPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.c I added #include "RASModel.H" in it and: Foam::tmp<Foam::volScalarField> Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff ( const volScalarField& p ) const { volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1))); volScalarField rho ( limitedAlpha1*rho1() + (scalar(1)  limitedAlpha1)*rho2() ); return //......I want to change it( <<k>> turbulence multiple in it): (3*rho1()*rho2())*sqrt(2/(3*rho1()))*k() *rRb(limitedAlpha1)/(rho*sqrt(mag(p  pSat()) + 0.01*pSat())); //.................................................. ...... } dont successful wmake, and seen(was not declared ): I/opt/openfoam201/src/OSspecific/POSIX/lnInclude fPIC c $SOURCE o Make/linux64GccDPOpt/SchnerrSauer.o phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C: In member function â€کFoam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff(const Foam::volScalarField&) constâ€™: phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C:113: error: k() was not declared in this scope make: *** [Make/linux64GccDPOpt/SchnerrSauer.o] Error 1 please help me. another way that I worked: Foam::tmp<Foam::volScalarField> Foam:haseChangeTwoPhaseMixtures::SchnerrSauer: Coeff ( const volScalarField& p ) const { volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1))); volScalarField rho ( limitedAlpha1*rho1() + (scalar(1)  limitedAlpha1)*rho2() ); return const volScalarField &k=U_.db().lookupObject<volScalarField>("k"); (3*rho1()*rho2())*sqrt(2/(3*rho1()))*k() *rRb(limitedAlpha1)/(rho*sqrt(mag(p  pSat()) + 0.01*pSat())); and wmake was Successfully . but How can I Understand that this k is the same with k(kinetic turbulence energy)???Is it true??? Because I replase another world in this : const volScalarField &M=U_.db().lookupObject<volScalarField>("M") but not error occured!!!!!!?????? please help me more about this problem????Thanks. 

February 24, 2014, 09:57 

#20 
Senior Member
Join Date: Mar 2009
Posts: 138
Rep Power: 10 
Hello Everybody,
I am exhuming this thread because it is the only one were somebody talks about yPlustherm. So maybe the Pro's see their old thread poping up and somebody can give me a reply: I am working with the above mentioned wall function. What I dont understand is the calculation of yPlustherm. The docu says it is the value of yPlus at the edge of the thermal laminar sublayer. But where does the formula come from: Code:
00084 scalar f = ypt  (log(E_*ypt)/kappa_ + P)/Prat; 00085 scalar df = 1.0  1.0/(ypt*kappa_*Prat); 00086 scalar yptNew = ypt  f/df; So my question is: What is yPlustherm and whats the background for this formula? So I would be gratefull if somebody could lift the fog. Thanks! Camoesas
__________________
OF  2.0.0 Last edited by camoesas; February 28, 2014 at 03:49. 

Thread Tools  
Display Modes  

