CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Pre-Processing (https://www.cfd-online.com/Forums/openfoam-pre-processing/)
-   -   Boundary Conditions for v2f turbulent model (https://www.cfd-online.com/Forums/openfoam-pre-processing/115572-boundary-conditions-v2f-turbulent-model.html)

bmikuz April 2, 2013 12:05

Boundary Conditions for v2f turbulent model
 
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 April 2, 2013 12:59

bmikuz,

The v^2-f 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 v^2 and f boundary conditions.

Inlet:

v^2 = \frac{2}{3}\,k_{inlet}
\frac{\partial f}{\partial n} = 0

Outlet:
zero gradient for both

Walls:

v^2 = 0.0

f = \frac{20\,\nu^2\,v^2}{\varepsilon_{inlet}\,y^4}

If you want the justification and full discussion, I'll refer you to the original v^2-f 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 v^2-f 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 v^2-f requires full resolution of the boundary layer, i.e. y+ <= 1.

Hope that helps, Chris

Quote:

Originally Posted by bmikuz (Post 417886)
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 April 2, 2013 16:18

@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:
Code:

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

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

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

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

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)
Code:

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

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 April 3, 2013 10:42

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 April 3, 2013 10:50

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.

Quote:

Originally Posted by bmikuz (Post 418112)
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 April 3, 2013 11:29

1 Attachment(s)
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 April 3, 2013 11:53

Quote:

Originally Posted by bmikuz (Post 418131)
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 April 3, 2013 18:52

3 Attachment(s)
Quote:

Originally Posted by cnsidero (Post 418139)
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 :rolleyes:
@all: when i don't right with my assumptions then let me know, please !!!

bmikuz April 4, 2013 05:34

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 April 4, 2013 19:30

5 Attachment(s)
Quote:

Originally Posted by bmikuz (Post 418274)
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:
Code:

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:rolleyes:
  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:confused:
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 April 4, 2013 19:33

3 Attachment(s)
ok the limit for upload stopped me, here is the rest....

bmikuz April 8, 2013 04:48

2 Attachment(s)
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 March 24, 2014 11:53

1 Attachment(s)
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:
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:
Code:

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

libs ( "myEpsilonLowReWallFunction.so" );
If you want to avoid this, you should get the same results with groovyBC:
Code:

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


canopus November 20, 2015 05:22

3 Attachment(s)
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
Code:

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

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 November 21, 2015 06:24

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 November 23, 2015 08:34

v2 wall boundary condition
 
Hi

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

Code:

  // 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 November 23, 2015 09:36

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 November 24, 2015 13:39

2 Attachment(s)
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 November 25, 2015 03:01

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

Quote:

Originally Posted by karlli (Post 574810)
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 January 16, 2016 16:36

"epsilon to zeroGradient" for Low Re Turbulence models?

Are you sure ?
In the tutorial folder

Quote:

tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0
for walls epsilon is set to

Code:

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



All times are GMT -4. The time now is 06:15.