[bmikuz]

Hello,

I noticed that the new version of OF 2.2.0 contains implemented v2f turbulent model. Now, I'd like to use that model for flow in a channel or for flow through a tube. The problem is that there is no prepared case using v2f turbulent model in $FOAM_TUTORIALS directory, where one can see how to set the boundary conditions (BCs).

I realized that you need 7 fields in 0 directory: U, p, nut, k, epsilon, f and v2.
Firstly, what are correct dimensions for f and v2? Is dimension (1/s) for f and (m^2/s^2) for v2 correct?
Secondly, what should be proper BC for f and v2 on the outlet? I followed the instructions for BC written in v2f.H file, but found them little incomplete...

Thanks in advance for your help or any comment.
Blaz

[cnsidero]

bmikuz,

The - model was one of the turbulence models I investigated for my doctoral research. From the literature the following is what I learned to use for and boundary conditions.

Inlet:




Outlet:
zero gradient for both

Walls:





If you want the justification and full discussion, I'll refer you to the original - paper by Paul Durbin and a PhD thesis by Andreas Sveningsson:

Durbin, P.A., 1991 “Near-Wall Turbulence Closure Modeling Without ‘Damping Functions”’ Theoretical and Computational Fluid Dynamics Vol. 3, pp. 1-13

Sveningsson, A. 2003 “Analysis of the Performance of Different - Turbulence Models in a Stator Vane Passage Flow”, Ph.D Thesis, Chalmers University of Technology ISSN 1101-9972
http://publications.lib.chalmers.se/...e-passage-flow

In particular, I found Sveningsson's thesis very helpful.

And by the way, the - requires full resolution of the boundary layer, i.e. y+ <= 1.

Hope that helps, Chris

Originally Posted by bmikuz

Hello,

I noticed that the new version of OF 2.2.0 contains implemented v2f turbulent model. Now, I'd like to use that model for flow in a channel or for flow through a tube. The problem is that there is no prepared case using v2f turbulent model in $FOAM_TUTORIALS directory, where one can see how to set the boundary conditions (BCs).

I realized that you need 7 fields in 0 directory: U, p, nut, k, epsilon, f and v2.
Firstly, what are correct dimensions for f and v2? Is dimension (1/s) for f and (m^2/s^2) for v2 correct?
Secondly, what should be proper BC for f and v2 on the outlet? I followed the instructions for BC written in v2f.H file, but found them little incomplete...

Thanks in advance for your help or any comment.
Blaz

[bscphil]

@bmikuz: I have do some simulations with the new v2f turbulence model (the incompressible version) in OpenFOAM 2.2.0 and give you my boundary conditions. I think it can help you

U:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (1 0 0);

boundaryField
{
    inlet
    {
      
        type            fixedValue;
        value           nonuniform List<vector> 
64
(
(0.111708 6.41525e-06 0)
(0.322171 -1.26309e-05 0)
(0.488089 5.5452e-06 0)
(0.601945 -1.16204e-07 0)
(0.681378 1.39714e-06 0)
(0.74007 3.13196e-06 0)
(0.785756 4.7772e-06 0)
(0.822917 6.06065e-06 0)
(0.854232 6.85654e-06 0)
(0.881348 7.05289e-06 0)
(0.905324 6.60611e-06 0)
(0.926875 5.52305e-06 0)
(0.946502 3.84961e-06 0)
(0.964567 1.65758e-06 0)
(0.98134 -9.65316e-07 0)
(0.997026 -3.92283e-06 0)
(1.01178 -7.11429e-06 0)
(1.02574 -1.04334e-05 0)
(1.03898 -1.37726e-05 0)
(1.0516 -1.70166e-05 0)
(1.0636 -2.00434e-05 0)
(1.075 -2.27193e-05 0)
(1.08579 -2.48985e-05 0)
(1.0959 -2.64265e-05 0)
(1.10529 -2.71472e-05 0)
(1.11386 -2.69128e-05 0)
(1.12153 -2.55954e-05 0)
(1.12819 -2.3057e-05 0)
(1.13373 -1.94304e-05 0)
(1.13809 -1.47821e-05 0)
(1.14115 -9.3311e-06 0)
(1.14277 -3.1689e-06 0)
(1.14277 3.30393e-06 0)
(1.14115 9.46547e-06 0)
(1.13809 1.49146e-05 0)
(1.13373 1.9562e-05 0)
(1.12819 2.31885e-05 0)
(1.12153 2.57287e-05 0)
(1.11386 2.70509e-05 0)
(1.10529 2.72948e-05 0)
(1.0959 2.65897e-05 0)
(1.08579 2.50852e-05 0)
(1.075 2.29391e-05 0)
(1.0636 2.03078e-05 0)
(1.0516 1.73391e-05 0)
(1.03899 1.41684e-05 0)
(1.02574 1.09186e-05 0)
(1.01178 7.70432e-06 0)
(0.997028 4.63039e-06 0)
(0.981344 1.79609e-06 0)
(0.964572 -7.09707e-07 0)
(0.946509 -2.80824e-06 0)
(0.926884 -4.43494e-06 0)
(0.905335 -5.54543e-06 0)
(0.88136 -6.12197e-06 0)
(0.854244 -6.18138e-06 0)
(0.822926 -5.77808e-06 0)
(0.785759 -5.0163e-06 0)
(0.740064 -4.00941e-06 0)
(0.68136 -2.93789e-06 0)
(0.601916 -1.92725e-06 0)
(0.488061 -7.47706e-06 0)
(0.322159 1.1563e-05 0)
(0.111704 -6.52362e-06 0)
)
;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    frontAndBack
    {
        type            empty;
    }
}

p:

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

epsilon:

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 3.825E-5;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 3.825E-5;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

k:

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 3.75E-3;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 3.75E-3;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            fixedValue;
        value           uniform 1.0E-15;
    }
    frontAndBack
    {
        type            empty;
    }
}

nut:

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
    walls
        {
//        type            calculated;
//        value           uniform 0;

      type            zeroGradient;

//      type            nutLowReWallFunction;
//      value           uniform 0;
    }

    frontAndBack
    {
        type            empty;
    }
}

v2: v2(inlet) = 2/3*k(inlet)

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 2.5e-3;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 2.5e-3;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            fixedValue;
        value           uniform 1e-10;
    }
    frontAndBack
    {
        type            empty;
    }
}

f:

dimensions      [0 0 -1 0 0 0 0];

internalField   uniform 1e-10;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            fixedValue;
        value           uniform 1e-10;
    }
    frontAndBack   
    {
        type            empty;
    }
}

I hope this helps you. Other assuptions for my case are:

1. the case is a "low-Re" turbulence modeling case without any wallFunctions
2. the walls have patch-type "wall", to use the "nutLowReWallFunction" for calculating the yPlus values
3. the yPlus values are between 0.3 and 7 (calculated with the "nutLowReWallFunction" --> yPlus should be approx 1 !! but results are ok...
4. for the velocity inlet condition i use a fully developed velocity profile. This is calculated befor e.g. with the "LaunderSharmaKE" turbulence model.
5. i have compared my results with experimental results from the ERCOFTAC Case 8.2 (Flow through an Asymetric diffuser) http://www.ercoftac.org/fileadmin/us...8.2/index.html
6. the results are good for the mean-velocity and the distributions of skin friction coefficients are plausible (i see transition and re-attechment points)

@cnsidero: thanks for the references

[bmikuz]

Thanks a lot @cnsidero and @bscphil. You were of great help!

I used @bscphil's BCs and a "tutorial" simulation works now. Since my real case is quite big and contains a lot of vanes inside a tube, I cannot afford very dense mesh in the boundary layer. So my next step is using wall functions, which are also implemented in that v2f turbulent model...

[cnsidero]

Good to hear.

Seems kind of dubious to me to use wall functions with the v2-f. The whole purpose of adding the v2 and f equations to the k-e model was to provide better scaling of the turbulent kinetic energy near the wall. The Durbin paper I reference discusses this at length. My point is, if you're going to use a wall-function with the v2-f, which will effectively override the benefit of the v2-f model, it's seems pointless to use it.

Just my $0.02.

Originally Posted by bmikuz

Thanks a lot @cnsidero and @bscphil. You were of great help!

I used @bscphil's BCs and a "tutorial" simulation works now. Since my real case is quite big and contains a lot of vanes inside a tube, I cannot afford very dense mesh in the boundary layer. So my next step is using wall functions, which are also implemented in that v2f turbulent model...

[bmikuz]

Ok, good to know that. I don't know much about the v2f turbulent model yet, but I checked the Description in the v2f.H file (attached to the post) and there are somehow recommended wall boundary conditions as follows:

k=kLowReWallFunction
epsilon=epsilonLowReWallFunction
v2=v2WallFunction
f=fWallFunction

There must be a good reason for implementing these wall functions in v2f model or isn't it?

[cnsidero]

Originally Posted by bmikuz

There must be a good reason for implementing these wall functions in v2f model or isn't it?

My guess is to make it more robust for practical applications. Specifically I mean for complex applications it may be difficult to generate mesh that meets y+ requirement everywhere.

But in general, the formulation of the v2-f was develop to handle resolved flows better than other 2 equation formulations of the k-e. I did notice in the v2f.H header they reference some newer papers, which I haven't read, but it's possible that they discuss near wall handling of v2 and f for y+ > 1.

So while it can handle y+ > 1, my understanding of its benefit is that you should try to keep y+ <= 1. Of course, some validation against some know applications that demonstrate the v2-f benefits might be a good exercise.

FYI, the near wall v2 and f behaviors are defined in v2WallFunctionFvPatchScalarField.C and fWallFunctionFvPatchScalarField.C where it checks for the yPlusLam value (edge of viscous sub-layer) and switches accordingly.

[bscphil]

Originally Posted by cnsidero

But in general, the formulation of the v2-f was develop to handle resolved flows better than other 2 equation formulations of the k-e. I did notice in the v2f.H header they reference some newer papers, which I haven't read, but it's possible that they discuss near wall handling of v2 and f for y+ > 1.

Ok, my understanding is the same, the v2f turbulence model is a socalled "low"-Re turbulence model (also for high Re-Numbers in the core flow) with no use of any wall functions, because of a resolved boundary layer. The differences between "high"- and "low"-Re models is in the attached figure "itw_vs_wtf.jpg".

For me it make no sense to use wall functions in the v2f-model (don't know what the OpenFOAM developer want to do with these ?? ). I have read all newer and older papers about the v2f model and i dosn't find any near wall handling of v2 and f for y approx 30 like in an "high"-Re model.

I read the v2f.H-file too and spend some time to simulate my diffuser testcase with the use of all wall functions ...

the results with the use of wall functions are very bad:

1. f and v2 becomes infinity field values
2. the residuals swing a lot and there was no convergence
3. the velocity and pressure field looks normal, but i don't have the near wall seperation bubble like in my simulation without wall functions
4. the mean velocity and skin friction distributions in comparison with the experimental data are very bad (see atteched plots "Ux_..." and "cf_..." --> the red lines are without wall functions and the green lines are with wall functions)

I have two possibilites for you:

1. refine your mesh at the no-slip walls to have y+<1 (or firstly y+<5-7)...use some boundary layers at your walls (5-10 layers maybe)
2. use another turbulence model, for example any of the "high"-Re turbulence models and then u can use your wall functions

@all: when i don't right with my assumptions then let me know, please !!!

[bmikuz]

Thanks for instructions to both of you. Since I already did a simulation with k-omega SST turbulence model (with wall functions), I'll have to follow your first suggestion @bscphil and add more layers at the walls. I hope that the mesh will not become significantly bigger than it is now (~40 millions cells).

About the @bscphil's testing of wallfunctions on diffusor: were the meshes of the same size in both cases (with and without wallfunctions)?

[bscphil]

Originally Posted by bmikuz

About the @bscphil's testing of wallfunctions on diffusor: were the meshes of the same size in both cases (with and without wallfunctions)?

I have good news the velocity field with the use of wall functions have good comparison to the experimental data

To your question, no the meshes don't have the same size, because of wall functions i use a coarse mesh at the walls. I think about my case and i find my mistake, i used for the simulation with wall functions the mapped field of my first post for the velocity, this field was calculated with the "low-Re" LaunderSharmaKE model and represent a fully developed channel flow with near wall resolution --> i cancel this and use for velocity at the inlet:

inlet
    {     
        type            fixedValue;
        value           uniform (1 0 0);
    }

Now i get better results with the use of wall functions, but i have some "building sites":

1. my residuals for v2, f and p are not so good, they swing a lot and are high
2. the fields for v2 and f have some differences to the same fields with the calculation without wall functions (see atteched pictures)
3. the skin friction distributions are ok for the lower wall (the seperation point is exactly the point of experiment, a wonder? ) and not so good for the upper wall in comparison with the experiments
4. the seperation bubble is at the right position
5. i don't have plot the pressure coefficients, maybe i can do this?

You said that you have a channel or tube flow with some vanes inside...i have some questions to you:

1. have you a full 3D tube? when yes, why don't you use a little part like a 10-45 degree piece of the tube for your calculation (with cyclic/periodic boundary conditions at the cutting planes)...then you can have a finer mesh
2. What's the Re-number for your case?
3. Which solver you use? --> just describe us your case, please

I hope i can help you for the first time?

(--> the red lines are without wall functions and the green lines are with wall functions !!)

[bscphil]

ok the limit for upload stopped me, here is the rest....

[bmikuz]

First, my apology for late reply.
Thanks for sharing your results, which are more promising for me.

My "tutorial" case was flow in a circular tube just to initiate/try v2f model. Now I want to use v2f for my "real" case, which is flow through fuel bundle of nuclear reactor. To be more specific, in my case fuel bundle is an array of 5×5 fuel rods which are held together with spacer grid (in the middle of sketch geo.jpg). At the outlet of spacer grid there are mixing vanes (spacer1.jpg), which serve for increasing the turbulence of the flow... So, the coolant (water) is flowing in a square channel along the rods. Reynolds number is around 50000, I want to run transient simulation with pisoFoam. Benchmark will be done with experimental measurements of velocity fields behind mixing vanes, so I'm not interested in skin frictions...

Since I'm interested in velocity fields, the simulation with wall functions may be good enough? I'll probably prepare both simulations (with & without wall functions) but it will take some more time.

[Maff]

Since I've been working with the v2f model for a while, I wanted to share with you some thoughts about the wall conditions.
For the v2f model implemented in openfoam, the boundary conditions at the wall for low-re cases are:
k=1e-10
v2=1e-10
f=1e-10 (because of the code-friendly modification by Lien)
nut=nutLowReWallFunction(i.e. =0)
epsilon is not zero nor zeroGradient.
epsilon = 2*nu*k1/y1^2, where k1 and y1 are the cell-center values at the wall.
The wall function implemented in openfoam is supposed to take care of it (epsilonLowReWallFunction).
But I had some problems using it, it makes the simulation diverge.
I had a look at the code:

if (yPlus > yPlusLam_)
        {
            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
        }
        else
        {
            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
        }

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

I think it should be:

if (yPlus > yPlusLam_)
        {
            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
			
            G[cellI] =
                w
               *(nutw[faceI] + nuw[faceI])
               *magGradUw[faceI]
               *Cmu25*sqrt(k[cellI])
               /(kappa_*y[faceI]);
        }
        else
        {
            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
        }

With this modification, the wall function is working for me.
If you want it, I'm attaching it to this post.
it is called myEpsilonLowReWallFunction.
just compile it and then add this to controldict:

libs ( "myEpsilonLowReWallFunction.so" );

If you want to avoid this, you should get the same results with groovyBC:

        type            groovyBC;
        valueExpression "2*nu*internalField(k)/sqr(mag(delta()))";
        value           uniform 1e-10;

[canopus]

Hi, I tried to use the v2f model in OF and for that I ran boundaryFoam and compared to Channel395 DNS.
For High Reynolds Number (HRN) model (y+ > 30) I got quite good results with OpenFOAM recommended settings i.e. as attached

Wall settings

k       = kLowReWallFunction                                        
 epsilon = epsilonLowReWallFunction 
v2      = v2WallFunction                                         
f       = fWallFunction
nut     = nutUSpaldingWallFunction;

But on trying to use in Low Reynolds Number (LRN) version (y+ <1 ) it gives weird result.

After googling and trying all possible combinations I got the best (yet not good) from following combination

Wall settings

k       = kLowReWallFunction                                         
epsilon = fixedValue 1e-10 
v2      = v2WallFunction                                         
f       = fWallFunction
nut     = zeroGradient

My observations are that most important factors are
setting nut = zeroGradient and epsilon = fixedValue (1e-10) at walls.

If you use any other as nutLowReWallFunction or epsilonLowReWallFunction results are way off.
For k,v2,f you can also use fixedValue (1e-10) and results are same.

Anyone else has same observations?
Or wants to test can ask for case files.

[Maff]

In my experience, results like that in your second figure are usually caused by an early divergence of the turbulence variables (you might see huge values for k or v2 near the walls).
It may still give you a result, but it will not be a physical one.
I wouldn't play too much with the boundary conditions, there are physical reasons to use nutLowReWallFunction and epsilonLowReWallFunction.
The initialization is usually the key: try and use a low-re k-epsilon model first, like LaunderSharmaKE (remember to use epsilon=0 at the wall with that model) or AbeKondohNaganoKE (http://www.cfd-online.com/Forums/ope...lon-model.html).

[Kire]

Hi

I have a question about the v2 boundary condition in v2WallFunction.

   // Set v2 wall values
    forAll(v2, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar uTau = Cmu25*sqrt(k[faceCellI]);

        scalar yPlus = uTau*y[faceI]/nuw[faceI];

        if (yPlus > yPlusLam_)
        {
            scalar Cv2 = 0.193;
            scalar Bv2 = -0.94;
            v2[faceI] = Cv2/kappa_*log(yPlus) + Bv2;
        }
        else
        {
            scalar Cv2 = 0.193;
            v2[faceI] = Cv2*pow4(yPlus);
        }

        v2[faceI] *= sqr(uTau);
    }

From the references Lien and Kalitzin (2001), it seems v2 approaches O(y^4) as y ---> 0, but how is v2[faceI] = Cv2*pow4(yPlus) derived? Where does this Cv = 0.193 come from?

Also, due to my poor c++ skills, I wonder why v2 gets multiplied by sqr(uTau) and assigned the value after the if loop? What is the purpose of that?

Thanks.

[Maff]

you will find everything in this paper:
Kalitzin et al. 2004. Near-wall behavior of RANS turbulence models and implications for wall functions.
In the if loop v2 is actually v2Plus, that's why you multiply it by sqr(uTau).

[karlli]

Dear all,
I'd just like to add my small piece of insight in this issue, after also struggling to understand the OF wall functions.

It turns out that there was a bug in the epsilonLowReWallFunction up to and including version 2.4.x, causing erratic behavior of this wall function. Just as canopus, I got some very wierd profiles for velocity and turbulent kinetic energy for a channel flow test case despite good initialization.

Using the wall function implementation in v.3.0.0 gives much better results. The alternative approach (setting v2,f,k to machine zero, epsilon to zeroGradient as explained here) also give reasonable results (green curve).

So to conclude, it seems like the new wall functions reduce to the low-Re approach (more or less) for small y+ values.

[akidess]

For sake of completeness, I'm guessing you are referring to bug ID 0001852, resolved by 9bf69ecd0bb19deb1bdfa5d2bb7034b6559242ea:

http://www.openfoam.org/mantisbt/view.php?id=1852

Originally Posted by karlli

Dear all,
I'd just like to add my small piece of insight in this issue, after also struggling to understand the OF wall functions.

It turns out that there was a bug in the epsilonLowReWallFunction up to and including version 2.4.x, causing erratic behavior of this wall function. Just as canopus, I got some very wierd profiles for velocity and turbulent kinetic energy for a channel flow test case despite good initialization.

Using the wall function implementation in v.3.0.0 gives much better results. The alternative approach (setting v2,f,k to machine zero, epsilon to zeroGradient as explained here) also give reasonable results (green curve).

So to conclude, it seems like the new wall functions reduce to the low-Re approach (more or less) for small y+ values.

[canopus]

"epsilon to zeroGradient" for Low Re Turbulence models?

Are you sure ?
In the tutorial folder

tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0

for walls epsilon is set to

lowerWall
    {
        type            fixedValue;
        value           uniform 1e-08;
    }