[maka]

Hi!

location: /OpenFOAM/OpenFOAM-1.3/src/LESmodels/incompressible/GenEddyVisc ; should also be checked in similar kind of models.

(1) In ChannelOodles Ueq we absorb laminar viscous diffusion into: sgsModel->divB(U). As a result,
B = R - 2 nu avr(D); avr is average which in the code implementation corresponds to D.
where,
R = (avr(U U) - avr(U) avr(U));
D = symm(grad(U)) = (1/2)(grad(U)+transpose(grad(U)));

(2) R = dev(R) + hyd(R);
where,
dev(R) = R - (1/3) trace(R) I; I is unit tensor;
hyd(R) = (1/3) trace(R) I = (2/3) k I;

(3) dev(R) is modeled as,
dev(R) = - 2 nuSgs dev(avr(D));

(4) using 3 in 2:
R = (2/3) k I - 2 nuSgs dev(avr(D));

(5) using 4 in 1:
B = (2/3) k I - 2 nuSgs dev(avr(D)) - 2 nu ( dev(avr(D)) + hyd(avr(D));
if we assume hyd(avr(D)) << dev(avr(D)) due to incompressibility (continuity) we have:
B = (2/3) k I - 2 nuEff dev(avr(D));
where,
nuEff = nuSgs + nu;
this is not the case in the implementation of GenEddyVisc::B() where nuSgs is used instead of nuEff [QUESTION 1].

(6) in GenEddyVisc::divB(volVectorField& U)
div(B) = - fvc::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(fvc::grad(U)().T()))

if we leave a side the implicit (fvm) and explicit (fvc) implementation of terms for a moment the above form lead to:
(I will drop the avr() for convenience and consistency with the code symbols)

div(B) = - ( div( nuEff (grad(U)) + div (nuEff dev(transpose(grad(U))) )

B = - nuEff (dev(grad(U) + hyd(grad(U))) + dev(transpose(grad(U))) )
if we assume hyd(D) << dev(D) due to incompressibility (continuity) we have:

B = - 2 nuEff ( dev(D)); this is not consistent with (5) ((2/3) k I) is missing; [QUESTION 2]


Best regards,
Maka

[henry]

> [QUESTION 1]

B currently returns the sub-grid stress generated by the sub-grid turbulence and does not include the laminar stress.

divB returns the divergence of sub-grid stress including the laminar stress for convenience and efficiency but as you say not consistency.

The options to make these more consistent would be to change the name of either B or divB (any suggestions?) or include the laminar stress in B, I am happy with either option.

> [QUESTION 2]:

For numerical reasons the ((2/3) k I) is subsumed into the pressure gradient as is common practice also in RANS modelling so again divB is deliberately inconsistent with B so perhaps divB should be renamed to make this clear.

[maka]

If B() name was changed to R(). and a comment is added to divB() to point out that it is:

divB() = div(B) - div((2/3) k I)] where the pressure is redefined

, then both inconsistencies are removed. Also this will be consistent with B definition in description part of model header files (B = 2/3*k*I - 2*nuEff*dev(D)).

Thanks for your explanations.

Best regards,
Maka.

[maka]

sorry, an after thought. divB can be renamed to somethings like divDevB (Dev for deviatoric) since this is what it is. Have a nice day!

Best regards,
Maka.

[henry]

Currently we use B to refer to the sub-grid stress and R the Reynolds stress which in LES runs is obtained by averaging the velocity fluctuations.

divDevB is as intersting idea but is this clarification and hence additional naming complexity necessary/useful? Perhaps it could be added as a comment.

[maka]

in that case no need to rename B() but rather change its definition to use nuEff instead of nuSgs.

Best regards,
Maka.

[henry]

In sub-grid stress models B is the stress corresponding to sub-grid turbulence and does not include the laminar stress so defining B from eddy viscosity models to include it would be inconsistent. I think we need a naming convention for the sub-grid stress including the laminar contribution.

[grtabor]

Total sub-grid stress vs. turbulent subgrid stress? In a manner of speaking the molecular viscosity is sub-grid-scale. Maybe B() for the latter and Beff() for the former (or something similar)

Gavin

[henry]

Yes I agree, the laminar stress is sub-grid and so there is no reason why it should not be included in the definition of sub-grid stress used by the momentum equation. I think Beff is a good idea, at least it is consistent with nuSgs and nuEff. Thanks for the thought.

[maka]

k() formulation in Smagorinsky (I did not check the rest of the models). The comment follows equation (23) in Fureby and Tabor 1997 (Theor. Comp. FD.) but the implementation uses the dev(D) instead of D. In incompressibel flow one would expect the D ~ dev(D) since hyd(D)<<dev(D) (continuity) BUT is there a specific reason for insisting on dev(D) or it is just a typing mistake? Thanks.

The comment in .H file:
//k = (2*ck/ce)*delta^2*||D||^2

The implementation.
//- Return SGS kinetic energy
// calculated from the given velocity gradient
tmp<volscalarfield> k(const tmp<voltensorfield>& gradU) const
{
return (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU)));
}

[henry]

D should be deviatoric for incompressible but numerics is not exact and div(U) != 0 exactly although div(phi) is zero to within solver tolerance. Using dev(D) ensures it deviatoric irrespective of errors in div(U).

[maka]

oneEqEddy SGS model V1.5

If we look at the form of production term in the k eq.

(1) G = -B:D

,where : is double inner product. following the above suggestion that D should always be dev(D) for incompressible flow leads to,

(2) G=-B:dev(D)

(3) B = (2/3) k I - 2 nuSgs dev(D).

putting 3 into 2,

(4) G = - (2/3) k I:dev(D) + 2 nuSgs dev(D):dev(D)

The implementation of G reads

volScalarField G = 2.0*nuSgs_*magSqr(symm(gradU));

The first part of G in (4) is missing and the second part is implemented using D instead of dev(D).

Is such deviation meant intensionally to save computational effort? Thanks.

Best regards,
Maka.

[henry]

For incompressible flow I:dev(D) = tr(D) = 0 and so dev(D) = D. However numerically on a collocated mesh this is not exactly true and there will be some difference between using D and dev(D) but it is likely to be small. Use whichever you prefer.

H

[taranov]

Hallo dear FOAMers!

Originally Posted by henry

>
> [QUESTION 2]:

For numerical reasons the ((2/3) k I) is subsumed into the pressure gradient as is common practice also in RANS modelling so again ...

If that is the case, is modified pressure (or pseudo-pressure) written to the data files "p" or is 2/3k removed and actual pressure written to file?

Best regards,
Andrey

[andrea]

Hi,
just to make it clear:
2/3 k I is an isotropic tensor just as p, so it's fine to include it in a modified pressure (Pope pag. 581):

P = p+ 2/3 k I

being k he residual kinetic energy.

I'm still not gettin why in the GenEddyVisc function divDevBeff() is added the term
- fvc::div ( nuEff * dev(div^T(U) )
maybe is it a correction term due to numerical reasons? But if it is so, since is explicitly computed, i.e. computed using the previous field, that U should be properly resolved and divergence free, there are comments on this or reference?
Thanks