CFD Online Discussion Forums - The density cross the interface is incorrect in the KEpsilon model in interFoam?

Recently I have found a paper (OMAE2015-41812) discussed about the KEpsilon model in the interFoam (2.3.0)， which claimed that "The existing formulations treat the density as a constant without considering its variations crossing the air and water free surface which is not correct. "

See the attached pic for more details.

I then checked KEpsilon.C, for the Dissipation equation

Code:

 

    tmp<fvScalarMatrix> epsEqn

    (

        fvm::ddt(epsilon_)

      + fvm::div(phi_, epsilon_)

      - fvm::laplacian(DepsilonEff(), epsilon_)

     ==

        C1_*G*epsilon_/k_

      - fvm::Sp(C2_*epsilon_/k_, epsilon_)

    );

and the turbulent kinetic energy equation

Code:

    tmp<fvScalarMatrix> kEqn

    (

        fvm::ddt(k_)

      + fvm::div(phi_, k_)

      - fvm::laplacian(DkEff(), k_)

     ==

        G

      - fvm::Sp(epsilon_/k_, k_)

    );

and G , G = $P_k$

Code:

volScalarField G(GName(), nut_*2*magSqr(symm(fvc::grad(U_))));

It seems the equation described in this paper is correct. However assuming the fluid is incompressible, for each phase (water or air), the density could be treated as a constant and remove from the Dissipation equation and energy equation. It does not sounds like a mistake. But of course we need to check the different density of water and air is accounted in the momentum equation.

Then I checked the UEqn in interFoam (2.3.0), which reads,

Code:

fvVectorMatrix UEqn

    (

        fvm::ddt(rho, U)

      + fvm::div(rhoPhi, U)

      + turbulence->divDevRhoReff(rho, U)

    );

and the divDevRhoReff(rho,U) also in KEpsilon.C,

Code:

tmp<fvVectorMatrix> kEpsilon::divDevRhoReff

(

    const volScalarField& rho,

    volVectorField& U

) const

{

    volScalarField muEff("muEff", rho*nuEff());



    return

    (

      - fvm::laplacian(muEff, U)

      - fvc::div(muEff*dev(T(fvc::grad(U))))

    );

}

It seems in the UEqn, the same rho is passed to the turbulence model.

Now I have the questions, which is this "const volScalarField& rho"?
It is rho_water, rho_air, or rho = alpha.water*rho_water + alpha.air*rho_air, with alpha.air = 1 - alpha.water?

Do you guys think the treatment in interFoam is correct or not?

The paper is linked here:
https://www.dropbox.com/s/53nvxo79m8...41812.pdf?dl=0