[sam1364]

Thanks for your answer. It was so helpful.
To me, it seems that using this kind of wall function is strongly dependent on mesh selection. I guess we should select a mesh for which we can apply log law for the first node. Is that correct?!
Besides, I would be happy if you could clarify me on this part of code

00063 {
00064 scalar ypl = 11.0;
00065
00066 for (int i=0; i<10; i++)
00067 {
00068 ypl = log(E*ypl)/kappa;
00069 }
00070
00071 return ypl;

Thanks,

Originally Posted by romant

the boundary condition nutkWallFunction sets the turbulent viscosity in the first node point based on the logarithmic law of the wall http://en.wikipedia.org/wiki/Law_of_the_wall . the difference between the nutUWallFunction and the nutkWallFunction lies in the way yPlus is calculated (in the first case based on the velocity close to the wall and in the latter case based on the turbulent kinetic energy close to the wall)

[idrama]

As you can probably see is the this formula just a fixed point formula which is being solved in 10 iteration. Take a pen and paper and compute a couple of steps by using k=0.4 and E=9.8 then you will get the result: 11.89658423819896. That is the value of the laminar sublayer. Further readings can be found in Versteeg, section: wall boundary conditions.

Cheers

[Anne Lincke]

Originally Posted by sam1364

To me, it seems that using this kind of wall function is strongly dependent on mesh selection. I guess we should select a mesh for which we can apply log law for the first node. Is that correct?!

Yes, this is correct.
Like always, when using wall functions, the first grid node is supposed to be in the logarithmic layer, i.e. y+ ~> 30

[romant]

You can use the yPlusRAS utility to check your y+ in your simulation. In order to use wall functions, it should be between 30 and 150. if you have less at certain points, make the cells bigger and if you have a higher value, refine cells in that region.

[sam1364]

Thanks.

I guess y+ should be greater than 11.89 (not 30) if we want to use wall function.

[sam1364]

Hey guys,

Another thing is that even for the cavity simulation in Open Foam Tutorial, I see that yPlus for fixed walls is always less than 11. Do you have any idea how important it is to obtain yPlus greater than 11?! Why even in tutorial, yPlus values for fixed walls are less than 11.89?
As you said, OpenFoam uses wall function for yplus>11.89, but it does not perform any wall treatment for yPlus<11.89. This is also another issue.

[Anne Lincke]

If y+ is greater than 11.xx you should use wall functions, that means the law of the wall is used for k,omega,epsilon (whatever) and nut.
If y+ is less than 11, the mesh is fine enough to resolve the equations without any wall functions.
For simulations without wall functions you need to set the right boundary conditions for k, omega and nut and choose a proper, lowRe turbulence model.

[romant]

Originally Posted by sam1364

Hey guys,

Another thing is that even for the cavity simulation in Open Foam Tutorial, I see that yPlus for fixed walls is always less than 11. Do you have any idea how important it is to obtain yPlus greater than 11?! Why even in tutorial, yPlus values for fixed walls are less than 11.89?
As you said, OpenFoam uses wall function for yplus>11.89, but it does not perform any wall treatment for yPlus<11.89. This is also another issue.

Tutorials are there to give an idea of how everything should be set up, they are not real life cases and results are not validated against experimental data.

There are different considerations regarding wall functions. Pope (Turbulent Flow, 2000) states that wall functions resolve everything below 30 and that the first grid point should be around 50, whereas Wilcox (Turbulence Modeling for CFD, 2006) states that 10<y+<100 give good results.

Since wall functions are there to simplify computations and reduce the computational cost, I would say keeping it at around 50 is not such a bad idea.

[romant]

Originally Posted by Anne Lincke

If y+ is greater than 11.xx you should use wall functions, that means the law of the wall is used for k,omega,epsilon (whatever) and nut.
If y+ is less than 11, the mesh is fine enough to resolve the equations without any wall functions.
For simulations without wall functions you need to set the right boundary conditions for k, omega and nut and choose a proper, lowRe turbulence model.

In the case y+<11 it is usually adopted to the refine the mesh so that y+<1 since then the mesh is fine enough to resolve most of the phenomena in the boundary layer (Pope, Turbulent Flow, 2000). Some turbulence models won't even function without having y+<1 and will crash (this is my own experience with OF)

[Anne Lincke]

Originally Posted by romant

In the case y+<11 it is usually adopted to the refine the mesh so that y+<1 since then the mesh is fine enough to resolve most of the phenomena in the boundary layer (Pope, Turbulent Flow, 2000). Some turbulence models won't even function without having y+<1 and will crash (this is my own experience with OF)

Yes, you are right. y+<1 is , of course, better than having a value around 10 or 11.

[sam1364]

thanks for sharing your comments and experiences.
Guys, I posted a question with a different topic two weeks ago but I have not heard any response for that. I do not know if anyone noticed that or not, so I am posting that here

I am trying to change the coefficients of k-epsilon model for incompressible cavity flow. I can just change the below values in 'RASProperties':

Cmu
c1
c2
sigmaEps

it seems that sigmak is not defined for incompressible k-epsilon model. I know that sigmak for standard k-epsilon model is 1. But, I need to edit this number. Does anyone know how I can change sigmak?!

Thanks

[romant]

Originally Posted by sam1364

thanks for sharing your comments and experiences.
Guys, I posted a question with a different topic two weeks ago but I have not heard any response for that. I do not know if anyone noticed that or not, so I am posting that here

I am trying to change the coefficients of k-epsilon model for incompressible cavity flow. I can just change the below values in 'RASProperties':

Cmu
c1
c2
sigmaEps

it seems that sigmak is not defined for incompressible k-epsilon model. I know that sigmak for standard k-epsilon model is 1. But, I need to edit this number. Does anyone know how I can change sigmak?!

Thanks

you could create your own kEpsilon model, where you basically copy this whole thing and edit the part where sigmaK is missing, just follow the code and you should be fine.

[sam1364]

I added sigmak in the following format. However, my results for sigamk=0.7 is no different than the one for sigmak=1. I guess pisoFoam solver takes sigmak=1 by default. Also, when I enter some letters as the value for sigmak (sigmak aaaa, I do not see any error in the solution which further indicates sigmak is not defined in pisofoam solver. I also tried sigmaK instead of sigmak, but results were the same.


kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.776;
sigmak 0.7;
sigmaEps 1.3;
}

[Bernhard]

As Roman said, you can easily check this in the source code:
src/turbulenceModels/incompressible/RAS/kEpsilon

In kEpsilon.C you find how the equations are set up and in kEpsilon.H how the diffusivity for k is defined (DkEff). You can easily update this model to ask for a sigmak:

Step 1. Let the model ask for a sigmak, for this you can peek at how this is done with sigmaEps, it will work the same way.
Step 2. Update the definition of DkEff, which is straightforward if sigmak is defined.
Step 3. Compile this into a library.

Good luck!

[sam1364]

Thanks guys,

I do not have any exprience in making a new library, but I will try to do that. Meanwhile, I have another serious problem with my lid driven cavity Problem. I am using Ras and k-epsilon to solve the problem.
To calculoate Reynolds stress, I do not if UPrime2mean is reynolds stress or R?! there is a utility called "R" which calculates Reynolds stress at specified time step. if R components are reynolds stress tensor, then what is UPrime2Mean?!
Besides, if R is Reynolds stress, do you have any idea how I can get the field averaging of R during code running?
FieldAverage function is just used to get field averaging for U and P, How can I edit that for R?!

[sam1364]

Thanks Bernhard,

I could add that part in a new library. Do you have any idea about my question for Reynolds Stress?

Thanks,
Pooyan

Originally Posted by Bernhard

As Roman said, you can easily check this in the source code:
src/turbulenceModels/incompressible/RAS/kEpsilon

In kEpsilon.C you find how the equations are set up and in kEpsilon.H how the diffusivity for k is defined (DkEff). You can easily update this model to ask for a sigmak:

Step 1. Let the model ask for a sigmak, for this you can peek at how this is done with sigmaEps, it will work the same way.
Step 2. Update the definition of DkEff, which is straightforward if sigmak is defined.
Step 3. Compile this into a library.

Good luck!

[sam1364]

Hi

I am now aware of how nut, epsilon and production are modified for the first grid point when using nutkWallFunction in k-epsilon model. But, I also need to know where in the codes, U is modified at the first point. As you know, U at the first grid point should follow the log law equation. Do you know where I can find the equation for modifying U at the first grid point?!


Thanks,
Pooyan

[Anne Lincke]

Hey Pooyan,
I asked myself the same question a while ago.
The answer is (in my opinion): There is no extra equation or boundary condition for the velocity at the first grid point which models the log law.
But: There is an equation for computing the wall shear stress at the boundary, via the modified nut at the wall.
This leads to an equation which contributes to the coefficients of the system of equations which has to be solved.
nut * up/ yp= nut/yp * up

So the coefficient nut/yp which includes the nut (modeling the log law) contributes to the equation which is solved for getting the values for the velocity at the first grid point.

I hope this helps.

[sam1364]

Thanks Anne

I have also looked through all subroutines and I could not find a place to impose wall behavior for U and I guessed that OpenFoam solves the momentum equation to obtain U at the first grid point which you have also mentioned. But, do we have nut in momentum equation which then can influence the solution for U?!!

Also, in epsilonwallfunction we have an equation for modified production which is

G[faceCellI] =
(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);

and Uw is given as

const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];

So, the Uw here represents the first grid point velocity and magGradUw[faceI] represents the normal gradient of Uw. this is why it was important for me to know how the first point velocity is calculated in OpenFoam. even for the normal gradient, I expect that it will be derived based on log law formulation, but I do not know how this parameter is calculated in OpenFoam.

[Anne Lincke]

Hey Sam,

Originally Posted by sam1364

Thanks Anne

I have also looked through all subroutines and I could not find a place to impose wall behavior for U and I guessed that OpenFoam solves the momentum equation to obtain U at the first grid point which you have also mentioned. But, do we have nut in momentum equation which then can influence the solution for U?!!

The diffusive term in the momentum equation contains nu_eff which is the sum of molecular and turbulent viscosity. So nut does have an influence on the solution for U.

Originally Posted by sam1364

Also, in epsilonwallfunction we have an equation for modified production which is

G[faceCellI] =
(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);

and Uw is given as

const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];

So, the Uw here represents the first grid point velocity and magGradUw[faceI] represents the normal gradient of Uw. this is why it was important for me to know how the first point velocity is calculated in OpenFoam. even for the normal gradient, I expect that it will be derived based on log law formulation, but I do not know how this parameter is calculated in OpenFoam.

Good question. I could not find a method for computing Uw, either... it looks like Uw coming from the turbulence model. It has to contain the log law but I could not find where :/