CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

A new wall function based on the velocity field----strange result

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   September 3, 2015, 22:37
Default A new wall function based on the velocity field----strange result
  #1
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Hi,all:
Recently, I was implementing a new type of wall function into OF based velocity. As done in OF, the wall function is implemented based on nut or nuSgs for eddy viscosity. While for non-linear model, it's not applicable. Instead, it is needed to implement a wall function based on velocity so that it is more general.
As in wall function theory, the tangential velocity is modified, the idea is to devide the velocity into the wall-normal velocity and tangential velocity:
U=Uver+Ut
the wall-normal velocity can be obtained by:
Code:
const vectorField wallnormal=U.snGrad()/max(mag(U.snGrad()),VSMALL);
    const vectorField Uver = (U&wallnormal)*wallnormal;
or you can use directly:
Code:
const vectorField wallnormal=patch.nf();
Then the tangential direction is:
Code:
const vectorField Ut=(U-Uver)/max(mag(U-Uver),VSMALL);
Now, we should calculate yPlusrefer to the wallShearStress utility)

Code:
 const vectorField patchnormal=-patch().Sf()/patch().magSf();
    const vectorField wallshearstress=patchnormal&devBeff;
scalar walls=mag(wallshearstress[facei]);scalar utau=sqrt(walls);
        scalar yplus=utau/(ry[facei]*nuw[facei]);//
At last, apply the judge condition based on yplus:
y>11, u+=ln(Ey+)/K
y<=11,u+=y+
Code:
        if (yplus>11)
        {
         Uw[facei]=(log(E_*yplus)/kappa_)*utau*Ut[facei]+Uver[facei];
        }
        else
        {
        Uw[facei]=yplus*utau*Ut[facei]+Uver[facei];
        }
The code is successfully compiled and run with channel flow case in turtorial(OF211), however, the result is incorrect, the pressure gradient is twice the DNS data

Does anyone have any idea?

here is the code:
USpaldingWallFunctionFvPatchVectorField.H
USpaldingWallFunctionFvPatchVectorField.C

Xianbei
huangxianbei is offline   Reply With Quote

Old   September 4, 2015, 22:35
Default
  #2
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Another problem I'm facing is the way I calculate the wallshearstress used for yplus calculation. In my code, the magnitude of wallshearstress vector is used, while if I change it to:
Ut&wallshearstress
to get the tangential component of wallshearstress, then the calculation break down
I don't know how do this occur!!

edit1: this is due to that the wallshearstress may be 0 and sqrt(0) is not allowed.However, the result is not correct. Does anyone have similar experience?

Last edited by huangxianbei; September 5, 2015 at 22:49.
huangxianbei is offline   Reply With Quote

Old   September 6, 2015, 07:48
Default
  #3
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Well, I found that it may be not possible to implement such a wall function in OF because that the boundary field of U would not correct during the calculation. Therefore, no matter how I change the formulation of Uw, the result is not affected. So it would only be possible to modify the SGS stress just like the process in nuSgsUSpaldingWallFunction.

Another problem is that the gradient operation can be used as following:
Code:
const fvPatchVectorField& U = lesModel.U().boundaryField()[patchi];
const symmTensorField S = symm(fvc::grad(U));
Is there anyone who can tell me how to use the velocity gradient in a wall function as nuSgsUSpaldingWallFunction?

Edit 1:solved, use "fvc:;grad(lesModel.U())"

Last edited by huangxianbei; September 6, 2015 at 09:03.
huangxianbei is offline   Reply With Quote

Old   September 7, 2015, 05:42
Default
  #4
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Hi, all:
After further investigation, I found that the reason why wall function on U got incorrect result was that the field never changed! Here is two steps of the calculations:

case1: without wall function:
Code:
Time = 0.2

Courant Number mean: 0.101138 max: 0.121578
DILUPBiCG:  Solving for Ux, Initial residual = 0.0152247, Final residual = 8.25801e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 3.17434e-06, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.0458469, Final residual = 1.27448e-07, No Iterations 3
DICPCG:  Solving for p, Initial residual = 1, Final residual = 0.0439203, No Iterations 4
time step continuity errors : sum local = 6.97808e-06, global = -4.07795e-19, cumulative = -4.07795e-19
DICPCG:  Solving for p, Initial residual = 0.0352028, Final residual = 5.3519e-07, No Iterations 37
time step continuity errors : sum local = 1.72154e-10, global = -1.35659e-19, cumulative = -5.43454e-19
Uncorrected Ubar = 0.133433	pressure gradient = 0.000336634
ExecutionTime = 0.04 s  ClockTime = 0 s

Time = 0.4

Courant Number mean: 0.101256 max: 0.121545
DILUPBiCG:  Solving for Ux, Initial residual = 0.0386084, Final residual = 1.19873e-08, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.333311, Final residual = 1.1725e-07, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.055974, Final residual = 4.27262e-08, No Iterations 3
DICPCG:  Solving for p, Initial residual = 0.883414, Final residual = 0.0426073, No Iterations 4
time step continuity errors : sum local = 6.88827e-06, global = -2.91983e-19, cumulative = -8.35436e-19
DICPCG:  Solving for p, Initial residual = 0.0802868, Final residual = 5.44705e-07, No Iterations 38
time step continuity errors : sum local = 6.15091e-11, global = -1.27589e-18, cumulative = -2.11133e-18
Uncorrected Ubar = 0.133557	pressure gradient = -8.86117e-05
ExecutionTime = 0.06 s  ClockTime = 0 s
case 2: with wall function on U boundary condition:
Code:
Time = 0.2

Courant Number mean: 0.101138 max: 0.121578
DILUPBiCG:  Solving for Ux, Initial residual = 0.0152247, Final residual = 8.25801e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 3.17434e-06, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.0458469, Final residual = 1.27448e-07, No Iterations 3
DICPCG:  Solving for p, Initial residual = 1, Final residual = 0.0439203, No Iterations 4
time step continuity errors : sum local = 6.97808e-06, global = -4.07795e-19, cumulative = -4.07795e-19
DICPCG:  Solving for p, Initial residual = 0.0352028, Final residual = 5.3519e-07, No Iterations 37
time step continuity errors : sum local = 1.72154e-10, global = -1.35659e-19, cumulative = -5.43454e-19
Uncorrected Ubar = 0.133433	pressure gradient = 0.000336634
ExecutionTime = 0.11 s  ClockTime = 0 s

Time = 0.4

Courant Number mean: 0.101256 max: 0.121545
DILUPBiCG:  Solving for Ux, Initial residual = 0.0386084, Final residual = 1.19873e-08, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.333311, Final residual = 1.1725e-07, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.055974, Final residual = 4.27262e-08, No Iterations 3
DICPCG:  Solving for p, Initial residual = 0.883414, Final residual = 0.0426073, No Iterations 4
time step continuity errors : sum local = 6.88827e-06, global = -2.91983e-19, cumulative = -8.35436e-19
DICPCG:  Solving for p, Initial residual = 0.0802868, Final residual = 5.44705e-07, No Iterations 38
time step continuity errors : sum local = 6.15091e-11, global = -1.27589e-18, cumulative = -2.11133e-18
Uncorrected Ubar = 0.133557	pressure gradient = -8.86117e-05
ExecutionTime = 0.13 s  ClockTime = 0 s
As can be seen ,they are exactly the same!

That means the boundary field is not modified during the calculation when wall function is applied! How can I make it take effect during the calculation?

In the code, fixedValueFvPatchVectorField::evaluate() is called. In the description, this function is to evaluate the patch field. How can it be that the patch field is not changed???

Any idea is welcome!

Xianbei
huangxianbei is offline   Reply With Quote

Old   September 9, 2015, 22:25
Default
  #5
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
As U.correctBoundaryConditions() is called in the solver, it's strange that U is not corrected. The correctBoundaryConditions is explained here, it will call the evaluate() function in the boundary condition applied. So in my code, U.correctBoundaryConditions() should call evaluate() and correct the velocity field, however, no change is seen.
huangxianbei is offline   Reply With Quote

Old   September 30, 2015, 03:54
Default
  #6
Member
 
Jason Tan
Join Date: Sep 2014
Posts: 47
Rep Power: 4
tzqfly is on a distinguished road
Quote:
Originally Posted by huangxianbei View Post
As U.correctBoundaryConditions() is called in the solver, it's strange that U is not corrected. The correctBoundaryConditions is explained here, it will call the evaluate() function in the boundary condition applied. So in my code, U.correctBoundaryConditions() should call evaluate() and correct the velocity field, however, no change is seen.
hello,
I read your thread and I know what you research is about DNS, and so do I, Now I face a hard question is that how do calculate y+ in DNS,you know there is no any codes about calculating y+ in DNS. I saw your answer in other thread about y+ about DNS,but you didn't show how,Can you show me? my email is : tzqfly2009@163.com
tzqfly is offline   Reply With Quote

Old   September 30, 2015, 07:52
Default
  #7
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 6
huangxianbei is on a distinguished road
Quote:
Originally Posted by tzqfly View Post
hello,
I read your thread and I know what you research is about DNS, and so do I, Now I face a hard question is that how do calculate y+ in DNS,you know there is no any codes about calculating y+ in DNS. I saw your answer in other thread about y+ about DNS,but you didn't show how,Can you show me? my email is : tzqfly2009@163.com
Hi:
For a DNS calculation , a simple way to achieve this is using the velocity gradient. If the mesh is fine, you can use gradU=u1/y1, here, 1 represents the value at the first node near-wall. tau_wall=nu*gradU, u_tau=sqrt(tau_wall), so you can calculate the y+

Xianbei
huangxianbei is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Define periodic velocity profile by field function at the inlet uhlondon STAR-CCM+ 3 June 25, 2015 08:44
Difficulty in calculating angular velocity of Savonius turbine simulation alfaruk CFX 8 December 3, 2013 17:51
non-orthogonal faces and incorrect orientation? nennbs OpenFOAM Native Meshers: blockMesh 7 April 17, 2013 05:42
Compile problem ivanyao OpenFOAM Running, Solving & CFD 1 October 12, 2012 09:31
Wall function in adverse pressure gradients stephane baralon Main CFD Forum 11 September 2, 1999 04:05


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