[AleDR]

Hi FOAMers!

I'm doing some programming works on k-Omega SST alternative forms in OF and I found some strange implementations in the incompressible versions of turbulence models.

In particular I'm wondering if the formula for the Reynolds stress tensor R, the deviatoric part of the effective stress tensor devReff and its divergence divDevReff are correct.

To my knowledge, the formulas for all of these quantities are the same for compressible and incompressible fluids, being the latter a particular case of the first (with the assumption that div(U_i) = tr(S_ij) = 0 , being S_ij = 1/2(d u_i/d x_j + d u_j/d x_i) the symmetric part of the grad U_i tensor, a.k.a the strain rate tensor).

Anyway the formulas in OF for the incompressible models seems a little puzzling because:

1. Reynolds stress tensor :
R_ij = 2/3 I_ij k - nut (2 S_ij - 2/3 div U_i I_ij)
in which for the incompressible case the last term drops.
The OF formulation is: R = ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)) where 2/3 div U_i I_ij has been dropped.

2. Deviatoric part of the effective stress tensor :
D_ij = nuEff (2/3 div U_i I_ij - 2 S_ij)
in which for the incompressible case the first term drops.
The OF formulation is: devReff = - nuEff()*dev(twoSymm(fvc::grad(U_))) where the divergence term has NOT been dropped, because of the dev operator.

3. Divergence of the deviatoric part of the effective stress tensor :
div(D_ij) = 2/3 d/d x_i( nuEff d U_k/d x_k ) - 2 d/d x_j( nuEff S_ij ) = 2/3 d/d x_i( nuEff d U_k/d x_k ) - d/d x_j( nuEff d U_i/d x_j + nuEff dU_j/d x_i )
in which for the incompressible case the first term drops.
The OF formulation is: divDevReff = - fvm::laplacian(nuEff, U) - fvc::div(nuEff()*dev(fvc::grad(U)().T())) where the divergence term has NOT been dropped, because of the dev operator.
Moreover it has been erroneously included as 1/3 div U_i , not as 2/3 div U_i I_ij, because operator dev is used instead of dev2 and there isn't any twoSymm in the argument, as for D_ij above!
From the point of view of strict numerics this shouldn't have any impact on the solution, because for incompressible fluids div U_i = 0. Anyway for the mathematics this is not correct.
Maybe if div U_j is not near machine zero, could this cause an error?

I'd like to know your opinions and in particular I'd like to know:
- Are there reasons for not keeping the divergence term in Reynolds stress tensor, but keeping it in the effective stress tensor?
- Can anybody explain me why there is dev in the divergence of the deviatoric part of the stress tensor instead of dev2 operator?

Thanks and let's wait for some answers!


Alex.