# Wall functions - questions about implementation

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Display Modes
 October 29, 2011, 08:08 Wall functions - questions about implementation #1 Senior Member   Robert Sawko Join Date: Mar 2009 Posts: 117 Rep Power: 13 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]); Are these comments right and if they are what is the point of having laminar viscosity in shear stress calculation and then calculating the gradient from log-law anyway? 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()); Q3) What is db object? It appears from for example here: Code:  volScalarField& G = const_cast(db().lookupObject(GName_)); Q4) Am I right that k wall function in OF is just zeroGradient? Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re 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: 154
Rep Power: 7
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:
 Originally Posted by AlmostSurelyRob 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?
I do not know nutUWallFunction but I use nutWallFunction.
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

nutWallFunction

Quote:
 Originally Posted by AlmostSurelyRob 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]); Are these comments right and if they are what is the point of having laminar viscosity in shear stress calculation and then calculating the gradient from log-law anyway?
The comments are right. The production term is \bar uv dU/dy.
\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:
 Originally Posted by AlmostSurelyRob Q4) Am I right that k wall function in OF is just zeroGradient?
Yes you are right, the k wall function implements a zero Gradient boundary condition.

Quote:
 Originally Posted by AlmostSurelyRob Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re 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 would use "calculated" for nut when using a low-Re model.

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: Stockholm, Sweden
Posts: 364
Rep Power: 12
Quote:
 Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re 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.
The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum)
__________________
~roman

May 15, 2012, 06:20
#4
Senior Member

Robert Sawko
Join Date: Mar 2009
Posts: 117
Rep Power: 13
Thanks for the comments Anne.

Quote:
 Originally Posted by Anne Lincke I do not know nutUWallFunction but I use nutWallFunction. This wall function is compatible with both, kqRWallFunction and omegaWallFunction/epsilonWallFunction.
Which version of OF are you using? I am on OF2.1.x and I do not have nutWallFunction at all. Only:

nutUWallFunction
nutkWallFunction

and their rough-wall 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 Launder-Spalding 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:
 Originally Posted by Anne Lincke 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?
No, I don't think it's valid. In fact when you look at some tutorials that use low-Re models you will see that that omega/epsilon and k have different BC. BC for k is fixedValue with value=0 and epsilon depending on the model is zero gradient or a value. Wilcox book has a good chapter about this.

Wilcox, Turbulence Modeling for CFD, Chapter 4 section 9 pp 138 Low Reynolds-Number 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 CFD-wiki about this (feel free to correct it you find anything!)
http://www.cfd-online.com/Wiki/Near-...k-omega_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:
 Originally Posted by Anne Lincke I would use "calculated" for nut when using a low-Re model. Anne
I was doing the same but romant comment... made me see the light!
Quote:
 Originally Posted by romant The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum) Anne
Many thanks romant.

Quote:
 Originally Posted by Anne Lincke 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?
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: 154
Rep Power: 7
Hey Robert,

Quote:
 Originally Posted by AlmostSurelyRob Which version of OF are you using? I am on OF2.1.x and I do not have nutWallFunction at all. Only: nutUWallFunction nutkWallFunction and their rough-wall equivalents.
I am using OpenFOAM-1.7.0. The nutWallFunction implementd there equals nutkWallFunction.

Quote:
 Originally Posted by AlmostSurelyRob 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 Launder-Spalding 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.
I think you are right. The basic difference between nutUWallFunction and nutkWallfunction is the computation of y+. In nutkWallFunction it is computed from the turbulent kinetic energy and the constant c_mu, whereas in nutUWallFunction y+ is simply read in.

For the implementation of the production term have a look at this thread

difference between kOmegaSST (OF-1.7.x) and kOmegaSST_lowRe (OF-1.6-ext)

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:
 Originally Posted by AlmostSurelyRob 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 CFD-wiki about this (feel free to correct it you find anything!) http://www.cfd-online.com/Wiki/Near-...k-omega_models 2) Read the description to kOmegaSST in OF. For some reason OF developers do not believe this sort of tricks.
There is a blending procedure in OpenFOAM, which is equal to the technique in the thread you posted. It switches between lowRe and highRe, so the omegaWallFunction can be used for both, lowRe and highRe.
Look at the source code of OF-2.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));
This is exactly the same as in the thread which was mentioned before.

December 27, 2012, 12:55
#6
Member

Malik
Join Date: Dec 2012
Location: Austin, USA
Posts: 52
Rep Power: 4
Quote:
 Originally Posted by romant The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum)
Hi,
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 k-epsilon 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: 154 Rep Power: 7 Hey, the nutWallFunction implements the wall function via the boundary condition for nut. Look at this thread nutWallFunction Kind Regards Anne

 January 8, 2013, 11:16 #8 Member   Malik Join Date: Dec 2012 Location: Austin, USA Posts: 52 Rep Power: 4 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: 52
Rep Power: 4
Quote:
 Originally Posted by Anne Lincke Yes you are right, the k wall function implements a zero Gradient boundary condition.
Hi,
I tried to understand why this wall-boundary-condition 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: 13 This near-wall 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 log-region, 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: 52 Rep Power: 4 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: 52 Rep Power: 6 Hi everybody, Is there any way to see the functions implementation in Openfoam? I want to see the definition of some turbulence-related 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: 41 Rep Power: 8 Hello Ali, they are hidden at /opt/OpenFOAM-2.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: 52 Rep Power: 6 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: 41 Rep Power: 8 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,274 Rep Power: 18 Hi all, Could anyone explain where the "k" equation (let's say in the k-epsilon model) does the equivalent to Code: epsilon_.boundaryField().updateCoeffs(); ... epsEqn().boundaryManipulate(epsilon_.boundaryField()); or (in k-omega) Code: omega_.boundaryField().updateCoeffs(); ... omegaEqn().boundaryManipulate(omega_.boundaryField()); , i.e. where in the code is the selected k-wall function actually called? __________________ The skeleton ran out of shampoo in the shower.

September 20, 2013, 07:56
#17
Senior Member

Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 364
Rep Power: 12
Quote:
 Originally Posted by RodriguezFatz Hi all, Could anyone explain where the "k" equation (let's say in the k-epsilon model) does the equivalent to Code: epsilon_.boundaryField().updateCoeffs(); ... epsEqn().boundaryManipulate(epsilon_.boundaryField()); or (in k-omega) Code: omega_.boundaryField().updateCoeffs(); ... omegaEqn().boundaryManipulate(omega_.boundaryField()); , i.e. where in the code is the selected k-wall function actually called?
If you take a look at the implementation of the kqRWallFunction, which is the corresponding wall function for k, you will notice, that k does not have a function updateCoeffs. The wall function also describes what it does, which is to set a zeroGradient boundary condition.

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_);
The zeroGradient boundary condition only has a member function evaluate
Code:
template<class Type>
void kqRWallFunctionFvPatchField<Type>::evaluate
(
const Pstream::commsTypes commsType
)
{
zeroGradientFvPatchField<Type>::evaluate(commsType);
}
The zeroGradient condition can be applied to many different solution variables, like U, p and T, which all do not use an updateCoeff() function before or after there solution.
__________________
~roman

 September 20, 2013, 07:59 #18 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,274 Rep Power: 18 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: Stockholm, Sweden
Posts: 364
Rep Power: 12
Quote:
 Originally Posted by RodriguezFatz 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 omega or epsilon wall function do not just set a boundary condition for the differential equations in the common sense, like dk/dy=0 or k=const. But they manipulate the value at the first point or on the wall depending on the properties or solution variables close to the wall, which is different at every point of a wall. For example, take a look at omegaWallFunction
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]);
}
where the indices faceI and faceCellI refer to the location of the center of the face on the wall and the locate of the center of the cell next to the wall, respectively.

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,274 Rep Power: 18 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 Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post syler3321 Fluent UDF and Scheme Programming 2 September 20, 2014 12:37 Attesz CFX 7 January 5, 2013 04:32 cristobal OpenFOAM 2 May 6, 2011 04:10 HFLUENT Fluent UDF and Scheme Programming 0 April 27, 2011 12:03 vkrastev OpenFOAM 2 January 24, 2010 18:45

All times are GMT -4. The time now is 17:14.

 Contact Us - CFD Online - Top