[ngj]

Hi all

There is a missing term in the RAS-model implementation in OpenFOAM 1.5.x, 1.5-dev and 1.6.x. This term is at least missing in the

kOmega (only present in 1.6.x)
kOmegaSST
kEpsilon

models. The problem is with the function

tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
    );
}

which should take the form

tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
      + 2.0 / 3.0 * fvc::grad(k_)
    );
}

This additional term, [ 2.0 / 3.0 * fvc::grad(k_) ], cause a driving force due to gradients in the normal part of the Reynold's stress tensor.

The other turbulence models might be affected as well, however as I am unfamiliar with those, I have not looked into the implementation.

Best regards,

Niels

P.S. My colleagues has done some preliminary tests and it seems that having this contribution co-aligned with the pressure gradients does not affect the results, however if they are perpendicular, the gradient in the normal Reynold's stresses is the only driving force.

[ngj]

Bump.

I can see that this term has not been included in the new 1.7.0 release. Any comments on post #1?

Cheers,

Niels

[mcjicpm2]

It is correct since the pressure in the simple scheme is relative so it has no affect on the integration on both side of the control volume this source term will disapear, It is correct

[ngj]

Hi Parham

I do not quite follow you. If there is a gradient in the turbulence, then there must be a driving force and this force is not included into the present releases of OpenFOAM, hence the term is important for specific purposes and needs to be included.

It has nothing to do with the simple-algorithm, as it is a more general bug inside the turbulence models.

Best regards,

Niels

[mcjicpm2]

Hi, the pressure in the code is relative isnt it? so if you add or remove -2/3k from or to it, it realy doesnt matter. I have seen this in other codes. since in the simple algorhm the pressure is relative pressure so p- or + any thing lets say 2/3 k should not matter.

[ngj]

Hi

You assume that k is constant, however as it is not, fvc::grad(k) results in a driving force.

We have verified that without this term experimental data cannot be reproduced under certain conditions, however including it yield good correlation between data and model.

The term comes from the Boussinesq approximation of the Reynolds stress tensor, hence the RANS formulation is incomplete due to its omission.

Best regards,

Niels

[l_r_mcglashan]

I agree with you. The 2/3 * k IS included in the UEqns.H file for bubbleFoam/twoPhaseEulerFoam which are the codes I use the most.

[ngj]

Hi Laurence

I see, I had not thought of looking into the individual solvers. However both simpleFoam and turbFoam uses divDevReff, however those could off course be changed instead.

Have a nice weekend,

Niels

[l_r_mcglashan]

Also, shouldn't 'dev' actually be 'dev2' in divDevReff.

This is confusing me a lot now. . I will try deriving the implementation from the theory over the weekend if I have time.

[l_r_mcglashan]

Do you have a copy of Turbulent Flows by Pope? Page 93 explains how the pressure, p, is modified to p + 2/3 k.

[ngj]

Hi Laurence

Unfortunately not, however would it be possible for you to email me the relevant page?

Best regards,

Niels

[alberto]

Hi,

you can simply re-write grad(p + 2/3k) = grad(p*). This will lead to a solution algorithm which is identical to what is implemented in OpenFOAM incompressible solvers.

@Laurence: in bubbleFoam/twoPhaseEulerFoam k is included in the fluid phase independently in each phase because its effect is different on each phase. In phase a there is the Ct coefficient. Notice also that you find the pressure gradient from the total continuity, and each phase has it scaled only by the corresponding phase fraction.

Best,

[l_r_mcglashan]

Thanks Alberto, we continued by email, apologies for leaving the thread dangling!

When would you (or anyone) alter the C_t coefficient, and do you know of any published material that has ever done that? I vaguely recall reading a paper where they did but I can't find it (if I could give 2 pieces of advice to anyone starting a cfd project, it would be BACKUP and organise your literature!).

[alberto]

Hi Laurence,

Ct should be related to the time scale of turbulence interaction between particles/bubbles and fluid turbulent structure. I do not remember where I saw this model tested changing Ct. First thing that comes to my mind is H. Rusche thesis however.

Alberto

[Andreas.Herwig]

Originally Posted by alberto

Hi,

you can simply re-write grad(p + 2/3k) = grad(p*). This will lead to a solution algorithm which is identical to what is implemented in OpenFOAM incompressible solvers.

@Laurence: in bubbleFoam/twoPhaseEulerFoam k is included in the fluid phase independently in each phase because its effect is different on each phase. In phase a there is the Ct coefficient. Notice also that you find the pressure gradient from the total continuity, and each phase has it scaled only by the corresponding phase fraction.

Best,

Hi,

I'm already searching a long time where this modified pressure is calculated grad(p + 2/3k) = grad(p*) in pisoFoam or the kEpsilon Code but can't find it. So if anyone knows ...

Thanks a lot.

Andreas

[alberto]

Originally Posted by Andreas.Herwig

Hi,

I'm already searching a long time where this modified pressure is calculated grad(p + 2/3k) = grad(p*) in pisoFoam or the kEpsilon Code but can't find it. So if anyone knows ...

Thanks a lot.

Andreas

If you write your momentum equation as a function of p*, and derive the pressure equation directly, you will find out that you do not have to do that calculation in the case of incompressible flows ;-)

[AleDR]

Dear Foamers,

Sorry for reopening this thread, but I still can't understand the point...

Ok for rewriting grad(p/rho + 2/3*k) as grad(p_*) in the momentum equation, but what about boundary conditions?

As a consequence of the pseudo-pressure formulation in /0 it will be specified one b.c. for p_* (not p) and one for k (which is part of p_*). Is this generally acceptable?

Moreover from what I see at first sight in /compressible/RAS/kEpsilon this formulation still holds for the term divDevRhoReff... will this mean that p used in perfectGas is p*=p + 2/3*k ?

I am really confused about the matter...


Thanks for any help!

.A.

[AleDR]

Hints? Suggestions? References?

...or maybe it's a silly question?

.A.