CFD Online Discussion Forums - A new wall function based on the velocity field----strange result

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - A new wall function based on the velocity field----strange result (https://www.cfd-online.com/Forums/openfoam-programming-development/158806-new-wall-function-based-velocity-field-strange-result.html)

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

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 yPlus:(refer 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:
Attachment 41875
Attachment 41876

Xianbei

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:eek:
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?

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())"

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

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.

Quote:

Originally Posted by huangxianbei (Post 563228)

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

Quote:

Originally Posted by tzqfly (Post 565932)

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