[noslip]

Hi everyone

I've got a question (actually two ...) about the implementation of the wall function for the LRR model in OpenFoam. Hope you can help me there.

The way I understand the implementation is the following:
The epsilonWallfunction fixes the value of epsilon in the first cell and modifies the production term "G". Then in the R equation, the production term "P" in the near wall cells is modified using the "G" production from the wall function. However, I don't quite understand the implementation below (from LRR.C). Does anyone have any further information or literature on on this procedure?

Additionally, in Wilcox' stress omega (https://turbmodels.larc.nasa.gov/rsm-stressomega.html) model the rapid pressure strain term contains an additional part labelled "D", which is similar to "P" but uses the transpose of the velocity gradients. Would this modification of the production also affect this "D" term?

Wall cell modification from LRR.C:

// Correct the trace of the tensorial production to be consistent
// with the near-wall generation from the wall-functions

const fvPatchList& patches = this->mesh_.boundary();

forAll(patches, patchi)
{
 const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
 forAll(curPatch, facei)
{
 label celli = curPatch.faceCells()[facei];
P[celli] *= min
( 
G[celli]/(0.5*mag(tr(P[celli])) + small) , 
1.0
);
 }
}
}

[edit] formatting ...

[C-L]

I can answer part of this question ..

If you look at the source code for the LRR model you will see that the pressure-strain correlation is modelled using only two terms rather than the more complicated version found in the Wilcox model reference. If you follow the LRR.H references you will see that in the LRR paper the authors suggest the pressure-strain can be modelled using either the more complicated version you have found, or by omitting the 'D' and 'S' terms and adjusting the model constants which is what is done in LRR.C.

In terms of the implementation it looks like the 'P' tensor needs to be corrected in the near wall cells so that its trace gives the right production of TKE, 'G', but it can't simply using the value of 'G' since 'G' is a scalar. 'P' is corrected by multiplying it with the ratio of G/0.5*tr(P) which guarantees that it will have the correct trace based on the wall function.

[noslip]

Thanks for your answer! That does make sense, funny way of implementing it though. I wonder why it was done this way instead of directly replacing parts of the Reynolds production term P with correct values based on G ... Maybe to avoid having to perform transformation to wall-local coordinates?

Does anyone have any experience with this procedure for the more complicated models, where D and S are retained?

[C-L]

It can't really be done by just using 'G' since 'G' is a scalar and we need a tensor field for 'P'. Since we already have an initial guess for 'P' we can just adjust each component by using the correct 'G' from the wall function while still preserving the distribution of separate components estimated from the transport equations. I.e if it was implemented in the way you are suggesting we would get equal distribution of 'G' across u'u', v'v', and w'w'.

If you want a more complicated model I suggest just using SSG instead unless you want to add in the extra terms yourself!

[C-L]

I just noticed that the source you gave provides Pi_ij in two forms: one using 'P' and 'S' and the other in the same form as the SSG model. It wouldn't take too long to code up the full LRR model using the SSG as a baseline if you had time!

I have done a similar thing when coding up an omega based RSM model so if you wanted some guidance then I can copy my form of Pi_ij. However, the simplest method for you would probably be to just use the SSG model!

[noslip]

That's actually a very good hint, I completely missed the alternative formulation. Thanks!

On the above point: what I meant was, you could create the production term in the following way (according to Apsley): set P.xx to 2*G and P.xy to a*G and the rest to zero, where a is the ratio of normal to shear stress (which is model dependent).

I would be very interested to see your implementation! I've got my own model which works very well as a low-Re model but trying to include the wall function is a pain
I will definitely try the alternate formulation.

[C-L]

My version is essentially the same as the SSG implementation but uses b_ij instead of a_ij in the source you gave. This changes three of the coefficients of Pi_ij by a factor of 2.

    volSymmTensorField b(dev(R)/(2*k_));
    volSymmTensorField S(symm(gradU));
    volTensorField Omega(skew(gradU));

You can split Pi_ij into an implicit part:

      + fvm::Sp(betaStar()*C1()*eps_/k_, R)

and explicit:

     + alpha*rho*(
            k_*
            (
                 (betaStar()*2.0*C1()*eps_/k_)*I/3.0
              + ((4.0/3.0)*(alphaHat()+betaHat())-gammaHat())*dev(S)
              + (2.0*(alphaHat()+betaHat()))*dev(twoSymm(b&S))
              + (2.0*(alphaHat()-betaHat()))*twoSymm(b&Omega)
            )

My coeffs are damping functions but the principle is the same!

Regarding your alternative for the correction of 'P' I'm not familiar with that other method so I'm not sure which would give better results! Although it does looking like your method would force R.xy as well as just the trace - all that the openFoam version is doing is making sure that 0.5 tr(R) matches G in the first cell.

[noslip]

Thank you very much! Just one last question for my understanding: when using your formulation, the pressure strain Pi_ij is completely unaffected by the wall function; any modifications would only apply to the production term Pij, correct?

[C-L]

As far as I can tell, yes. However, wall functions will still come into play since b and S are constructed from R and U which will have wall functions applied to them (unless you integrate through the full boundary layer).

Also, I haven't included a wall-reflection term but this can be added in the same way as in the LRR.C I think.

[noslip]

This has been very helpful, thank you.

[lpz456]

Recently，I am doing some research with RSM model. Now, I start to simulate backwardFacingStep tutorial in OFv2016 version. But the results always diverged, even only simulating by k-epsilon model. So I hope to find some help on the forum, I think you have a very deep understanding of RSM model. Would you like to put forward some opinions on convergence?
Thank you much!