
[Sponsors] 
October 29, 2011, 08:08 
Wall functions  questions about implementation

#1 
Senior Member
Robert Sawko
Join Date: Mar 2009
Posts: 117
Rep Power: 16 
Dear Foam Users and Developers,
I am really interested in understanding wall functions well. The main reason for this interest is that I would like to know precisely which RAS model is compatible with each of those functions. Therefore, I embarked on the quest to understand the implantation in OF. To paraphrase one sketch, I would like to emerge one day from this journey, probably a little bruised, but hopefully a little wiser. So here I am after a few days of literature review and perusing the code with a set of compiled questions  not yet wiser, but not too bruised as well. I am hoping to get both with this post. I am looking at OF 2.0 Q1) nutUWallFunction stands for Standard Wall function. yplus is estimated with at most 10 Newton iterations under given U and y and then turbulent viscosity is modified. According to the literature e.g. Wilcox having u_\tau allows to adjust k and epsilon/omega. But in OF I don't see any equivalent wall boundary conditions for k and omega. Instead there is an implementation which changes omega/epsilon and adjusts the production term of k. In fact k and omega/epsilon have only one wall function each. Is nutUWallFunction compatible with epsilon and kqR wall functions? Q2) G is k production i.e. \bar{uv} dU/dy. Please look at the comments, the code is taken from omega wall function. Code:
G[faceCellI] = //This is an estimate of shear stress \bar uv (nutw[faceI] + nuw[faceI]) *magGradUw[faceI] //This is an estimate of dU/dy from the logarithmic velocity profile *Cmu25*sqrt(k[faceCellI]) /(kappa_*y[faceI]); Q2') When is this actually invoked? k equation is always solved second in RAS. omega/epsilon is always solved before and it has this suspiciously looking method boundaryManipulate. Code:
omegaEqn().boundaryManipulate(omega_.boundaryField()); What is db object? It appears from for example here: Code:
volScalarField& G = const_cast<volScalarField&>(db().lookupObject<volScalarField>(GName_)); Q5) Finally, what is the correct boundary condition for nut in kepsilon lowRe models? Should I put calculated or should I put nutLowReWallFunction or... is there no difference? The latter seems to have only constructors and cloning functions. I will appreciate any comments or hints on the above questions. I find OF implementation very neat but as I outlined above I have still some doubts and gaps in my knowledge. 

May 14, 2012, 09:02 

#2  
Senior Member
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 168
Rep Power: 9 
Hey Robert,
I am also investigating in wall functions and how they are implemented in OpemFOAM. Maybe, if we merge our knowledge, we'll become wiser Quote:
This wall function is compatible with both, kqRWallFunction and omegaWallFunction/epsilonWallFunction. It uses the log law to derive the boundary condition for nut, like can be read in this thread http://www.cfdonline.com/Forums/ope...lfunction.html Quote:
\bar uv is equal to (u_tau)^2 = tau_w= (nu + nut) dU/dy = (nutw[faceI] + nuw[faceI])) *magGradUw[faceI] and dU/dy = u_tau/ (kappa*y) = c_mu^(0.25) * sqrt(k) But I ask myself, what happens when using a low Reynolds turbulence model? u_tau is computed according to the log law. This is only valid for a high Reynolds turbulence model using wall functions. Is this ansatz right for a computation without wall functions and y+ ~1? Quote:
Quote:
I am concerning about the production term, especially what happens without wall functions for the boundary condition. According to the implementation, G is computed using the log law, no matter how big y+ is. Isn't this wrong? Best, Anne 

May 14, 2012, 09:49 

#3  
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: London, UK
Posts: 368
Rep Power: 14 
Quote:
__________________
~roman 

May 15, 2012, 06:20 

#4  
Senior Member
Robert Sawko
Join Date: Mar 2009
Posts: 117
Rep Power: 16 
Thanks for the comments Anne.
Quote:
nutUWallFunction nutkWallFunction and their roughwall equivalents. Have you seen this: http://http://www.tfd.chalmers.se/~lada/postscript_files/jonas_report_WF.pdf Please look at Standard wall function (p 11) and LaunderSpalding methodology (incidentally, also p 11. ;)) This is I believe the difference between nutU* and nutk* wall function. Interestingly, the latter were implemented first and this is why I think omegaWallFunction is and kqRWallFunction were also tailed to this methodology. Quote:
Wilcox, Turbulence Modeling for CFD, Chapter 4 section 9 pp 138 Low ReynoldsNumber Effects Have a look if you haven't done so already. An interesting aspect is that some models employ certain blending procedures. I don't think this is at present implemented in OF. Two comments here: 1) I once expanded the article in CFDwiki about this (feel free to correct it you find anything!) http://www.cfdonline.com/Wiki/Near...komega_models 2) Read the description to kOmegaSST in OF. For some reason OF developers do not believe this sort of tricks. Nevertheless, I think we could probably implement these low to high switching techniques in OF. Are you interested in doing that? It's just that they might not be valid, but they should at least provide correct limiting behaviours for very low wall y+ and very high y+. Quote:
Quote:
I believe you're right here. Shall we run some simple case with y+~1 to see what is the difference between omegaWallFunction with and without this altered production term? That should be fairly straightforward and should settle the matter, shouldn't it? 

June 8, 2012, 09:30 

#5  
Senior Member
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 168
Rep Power: 9 
Hey Robert,
Quote:
Quote:
For the implementation of the production term have a look at this thread http://www.cfdonline.com/Forums/ope...tml#post361420 There is an interesting remark, which states that the computation of the production term is correct for low Reynolds number, as it is related to turbulent kinetic energy, which is very small at the viscous sublayer. So the production term also goes to zero in this region. Quote:
Look at the source code of OF2.0.0. HTML Code:
scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI])); scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]); omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog)); 

December 27, 2012, 12:55 

#6  
Member
Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 53
Rep Power: 7 
Quote:
I am also wondering what is hidden behind all thos wall conditions, as I want to know which range of y+ is allowed for a given wall condition. based on what you said romant, the nut file in a case does not define the wall function that will be applied to the case. But where is it defined if you have for example a kepsilon turbulence model ? Indeed kWall and epsilonWall functions are only zeroGradient BC. Thx for your reply. 

January 2, 2013, 04:31 

#7 
Senior Member
Anne Gerdes
Join Date: Aug 2010
Location: Hamburg
Posts: 168
Rep Power: 9 
Hey,
the nutWallFunction implements the wall function via the boundary condition for nut. Look at this thread http://www.cfdonline.com/Forums/ope...lfunction.html Kind Regards Anne 

January 8, 2013, 11:16 

#8 
Member
Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 53
Rep Power: 7 
Thanks a lot, this thread answered clearly my questions and even more !


January 16, 2013, 05:51 

#9  
Member
Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 53
Rep Power: 7 
Quote:
I tried to understand why this wallboundarycondition was used in the k epsilon model, and actually it seems weird to me that the gradient of k near the wall could be zero. Tell me if I am wrong but on the wall : k = utau²/sqrt(Cmu) dk/dn = d[utau²/sqrt(Cmu)]/dn with utau = sqrt(nu*du/dn) we have dk/dn = d[ (nu*du/dn)/sqrt(Cmu) ]/dn dk/dn = [nu/sqrt(Cmu)] * d²u/dn² if dk/dn =0 then this means that on the wall there is no viscosity constraint which is false. I must be wrong somewhere but I can't find where... If you know, let me know ! Thx for all ... to all ! 

January 16, 2013, 06:27 

#10 
Senior Member
Robert Sawko
Join Date: Mar 2009
Posts: 117
Rep Power: 16 
This nearwall modelling business is so exciting. :)
My understanding is the following. From the point of view of physics k=0 at the wall. Why? Because by definition k=1/2 u_i'u_i' and U=0 and therefore u'=0 at the wall. Meaning no fluctuations at the wall. But when you do wall functions we approximate the behaviour of the near wall region by using various experimental and theoretical results accounting for turbulence boundary layer structure. You must look up viscous, buffer and logarithmic regions in the literature. The expression you give k = utau²/sqrt(Cmu) has been established for logarithmic subregion only (300>y+>30). This is also called constant shear layer (correct me if I am wrong!). Experiments show that it is an okay approximation for log region. Note that wall functions impose the value of k in the centre of the cell (rather than the face). dk/dn is just a way of saying  the wall will have no further influence on the value in the cell and will approximate any gradients as if shear in the cell was constant. This is a good for logregion, but if you have a much finer mesh (y+~1) it is simply not true. Hope this description is a good approximation of what is really happening and hope it helps you on your quest! :) 

January 16, 2013, 08:38 

#11 
Member
Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 53
Rep Power: 7 
Hi Robert, thanks for the answer.
I agree with you k = 0 at the wall. For the relation I gave at first, k = utau²/sqrt(Cmu) it seems to be only applicable near the wall because in this region the generation of turbulence is equal to its dissipation see (4.2.18 in the link below): http://staffweb.cms.gre.ac.uk/~ct02/...is/node54.html If I got it, when we say zeroGradient for k on the wall, it is a way to avoid any more calculation as we already use the wall functions, which is enough ? But then we should do it for every variable... Actually, there must be something I still don't understand 

January 23, 2013, 10:18 
function definitions

#12 
Member
Ali Khalifesoltani
Join Date: Mar 2011
Location: Esfahan, Iran
Posts: 54
Rep Power: 9 
Hi everybody,
Is there any way to see the functions implementation in Openfoam? I want to see the definition of some turbulencerelated functions but don't know how?! Thanks in advance 

January 23, 2013, 12:33 

#13 
Member
David GISEN
Join Date: Jul 2009
Location: Germany
Posts: 54
Rep Power: 10 
Hello Ali,
they are hidden at /opt/OpenFOAM2.1.0/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions (your path may differ.) Regards, David 

January 23, 2013, 14:48 

#14 
Member
Ali Khalifesoltani
Join Date: Mar 2011
Location: Esfahan, Iran
Posts: 54
Rep Power: 9 
Vielen Dank David,
But is there any easier way to find the implementation of some functions that are declared in *.C files in the Openfoam website? For instance if I want to see the implementation of the function DepsilonEff() at line 233 of kEpsilon.C file at link below what should I do? http://foam.sourceforge.net/docs/cpp/a08551_source.html 

January 24, 2013, 04:08 

#15 
Member
David GISEN
Join Date: Jul 2009
Location: Germany
Posts: 54
Rep Power: 10 
Just follow the links, I guess:
"Return the effective diffusivity for epsilon. Definition at line 129 of file kEpsilon.H." (Note: better to have this kind of conversations in a separate thread.) Regards, David 

September 20, 2013, 07:24 

#16 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 20 
Hi all,
Could anyone explain where the "k" equation (let's say in the kepsilon model) does the equivalent to Code:
epsilon_.boundaryField().updateCoeffs(); ... epsEqn().boundaryManipulate(epsilon_.boundaryField()); Code:
omega_.boundaryField().updateCoeffs(); ... omegaEqn().boundaryManipulate(omega_.boundaryField());
__________________
The skeleton ran out of shampoo in the shower. 

September 20, 2013, 07:56 

#17  
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: London, UK
Posts: 368
Rep Power: 14 
Quote:
I am not completely sure, but I think the zeroBoundary condition is evaluated while solving kEqn in Code:
kEqn().relax(); solve(kEqn); bound(k_, kMin_); Code:
template<class Type> void kqRWallFunctionFvPatchField<Type>::evaluate ( const Pstream::commsTypes commsType ) { zeroGradientFvPatchField<Type>::evaluate(commsType); }
__________________
~roman 

September 20, 2013, 07:59 

#18 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 20 
Ok, for the "k=0" boundary condition, we don't change the equations at all, since only a zero is added.
Now for dk/dn=0, why is the "evaluate" function for k called. I am no professional programmer, but can read the code. Why is the "omega" b.c. explicitly called by the turbulence model, where "k" b.c (you were showing in your post) isn't?
__________________
The skeleton ran out of shampoo in the shower. 

September 20, 2013, 08:09 

#19  
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: London, UK
Posts: 368
Rep Power: 14 
Quote:
The code reads Code:
forAll(mutw, faceI) { label faceCellI = patch().faceCells()[faceI]; scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI])); scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]); omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog)); G[faceCellI] = (mutw[faceI] + muw[faceI]) *magGradUw[faceI] *Cmu25*sqrt(k[faceCellI]) /(kappa_*y[faceI]); } As you can see, not only is the value for omega in the first cell manipulated, but also for G in the first cell, where G is the average turbulence production term sometimes also called
__________________
~roman 

September 20, 2013, 08:15 

#20 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 20 
So what you are saying is that omega and epsilon wall functions have to change the cell value itself, so thay need to be explitcitly called and "overwrite" the calculated cell value. Whereas the k boundary condition defines the value at the cell's boundary and thus just changes theboundary condition of the wall cell but not the value in the first cell. Thus, k boundary condition changes something that is more hidden in the solving routine.
Is that correct?
__________________
The skeleton ran out of shampoo in the shower. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
UDF Wall functions in Fluent  syler3321  Fluent UDF and Scheme Programming  2  September 20, 2014 12:37 
Water subcooled boiling  Attesz  CFX  7  January 5, 2013 04:32 
wall functions and boundary conditions in OF 1.7.1  cristobal  OpenFOAM  2  May 6, 2011 04:10 
UDF for wall slipping  HFLUENT  Fluent UDF and Scheme Programming  0  April 27, 2011 12:03 
komega simulation without wall functions  vkrastev  OpenFOAM  2  January 24, 2010 18:45 