[jaswi]

Dear OpenFOAM Users

I am trying to change interFoam for turbulent flows. I am aware of the fact that there exists rasinterFoam as well for that purpose.

My questions are based upon the Henrik Rusche's thesis, where he, at length, explains both:
1)two-fluid model -->implemented in bubbleFoam,
2)interface capturing method --> implemented in interFoam

Now I have done some reading, and taken apart the basic equations and this is what I understood. Please take a look and comment:

according to Rusche thesis (eq 3.23) the Effective Reynolds stress is defined as :

R_phi_effective = -nuEff_phi(grad(U_phi) + grad(U_phi) * T( )) - 2.0/3.0 * deltai,j * grad(U_phi) + 2.0/3.0 * deltai,j * k_phi -- (1)

In general, the viscous term of momentum equation for compressible flow reads as:

div[nu(U i,j + U j,i) - 2.0/3.0 * deltai,j * U k,k)] -- (2)

for an incompressible flow the dilatational part can be neglected and that leaves :

div[nu(U i,j + U j,i)] --(3) --> in OpenFOAM -->interFoam this in UEqn.H is represented as :

- [fvm::laplacian(muf, U) + (fvc::grad(U) & fvc::grad(muf))] --(4)

From RANS point of view, for a time averaged momentum equation for incompressible flow the viscous term reads as:

div[nu(U i,j + U j,i) - rho*(UiUj) ] --(5)

Now according to Boussinesq assumption reynolds stresses are linked to velocity gradients via turbulent viscosity thus replacing

rho*(UiUj) = { -nut(U i,j + U j,i) } + { 2.0/3.0 * deltai,j * rho * k } --(6)

where second term { 2.0/3.0 * deltai,j * rho * k } is to account for the contraction.

inserting for rho*(UiUj) into (5) we get:

div[ nu(U i,j + U j,i) - { -nut(U i,j + U j,i) } + { 2.0/3.0 * deltai,j * rho * k } ] --(7)

=>

div[ nu(U i,j + U j,i) + nut(U i,j + U j,i) - { 2.0/3.0 * deltai,j * rho * k } ]

=>

div[ { (nu + nut) (U i,j + U j,i) } - { 2.0/3.0 * deltai,j * rho * k } ]

=>

div[ (nu + nut) (U i,j + U j,i) ] - div [ 2.0/3.0 * deltai,j * rho * k ] --(8)

-------------------------------------------------------------------------------- ----------------------------
Now switching to rasInterFoam :

surfaceScalarField muf = twoPhaseProperties.muf() + fvc::interpolate(rho*turbulence->nut());

and then UEqn the viscous term is calculated as

- [fvm::laplacian(muf, U) + (fvc::grad(U) & fvc::grad(muf))]

which is equivalent to the 1st term of Eqn (8).

-------------------------------------------------------------------------------- ---------

My question are thus :

1) Is it correct to include the 2nd term div [ 2.0/3.0 * deltai,j * rho * k ] into UEqn

2) will that improve the turbulence calculation and thus the calculation of Effective Reynolds stresses

I would be grateful if someone could please provide some feedback.

With Best Regards
Jaswinder