Wall functions  questions about implementation
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] = 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 = 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. 
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 
Quote:

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:
Quote:

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])); 
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. 
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 
Thanks a lot, this thread answered clearly my questions and even more !

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 ! 
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! :) 
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 :D 
function definitions
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 
Hello Ali,
they are hidden at /opt/OpenFOAM2.1.0/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions (your path may differ.) Regards, David 
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 
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 
Hi all,
Could anyone explain where the "k" equation (let's say in the kepsilon model) does the equivalent to Code:
epsilon_.boundaryField().updateCoeffs(); Code:
omega_.boundaryField().updateCoeffs(); 
Quote:
I am not completely sure, but I think the zeroBoundary condition is evaluated while solving kEqn in Code:
kEqn().relax(); Code:
template<class Type> 
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? 
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? 
All times are GMT 4. The time now is 06:43. 