[usv001]

Hello everyone,

I am a bit confused about the energy equation implementation in rhoCentralFoam. Particularly, why does the viscous term 'sigmaDotU' appear in the (supposedly) Euler equation?

surfaceScalarField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
              + (mesh.Sf() & fvc::interpolate(tauMC))
            )
            & (a_pos*U_pos + a_neg*U_neg)
        );

        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

Here's why I think it should not:
- Kuragnov flux scheme is derived for strictly hyperbolic equations which are the Euler equations in this case
- 'a_pos' and 'a_neg' are calculated using the eigenvalues of Euler equations and hence, the flux scheme should not be applied to the viscous term 'sigmaDotU'

However, when I modified the solver by removing the viscous term from this equation and introducing it in the internal energy/enthalpy equation later, the results were not as good as the current solver (more diffusion, shocks were not sharp, unexpected flow separation). So, I am guessing that whatever the current solver is doing is correct but I can't seem to understand why.

There are two other threads that I managed to find asking the exact question but both remain open:
1. Viscous terms in rhoCentralFoam
2. energy equation in rhoCentralFoam

I am hoping that someone would point in the right direction. Thanks to all!

[mkhm]

Did you find any answer for your questions ? I have the same exact questions.

[usv001]

Hi Mary,

I am afraid that I do not have a complete answer but let me give my opinion and, hopefully, someone will confirm/correct it.

The first thing to note is that if the case is inviscid, then muEff=0. So, the term 'sigmaDotU' does not have any effect on the results for inviscid flows.

The second thing to note is that the energy equation is solved using a predictor step for rhoE:

  fvm::ddt(rhoE)
+ fvc::div(phiEp)
- fvc::div(sigmaDotU)

and a corrector step for e:

  fvm::ddt(rho, e) - fvc::ddt(rho, e)
- fvm::laplacian(turbulence->alphaEff(), e)

When you add these two steps together you'll recover the correct energy equation. So, there is nothing overtly wrong with this splitting. However, for some reason the developers of rhoCentralFoam decided to place the 'sigmaDotU' term in the predictor step. Perhaps, they were aiming to get a good estimate of rhoE in the predictor step by including all the possible terms (if so, why didn't they do the same for momentum equation with 'tauMC'?).

a_pos and a_neg are essentially interpolation factors to compute an intermediate state at each face based on the pos and neg values obtained by TVD reconstruction. So, (a_pos*U_pos + a_neg*U_neg) can be seen as some intermediate velocity Ustar. This intermediate velocity is used in the calculation of 'sigmaDotU'.

surfaceScalarField sigmaDotU
(
    "sigmaDotU",
    (
        fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
      + (mesh.Sf() & fvc::interpolate(tauMC))
    )
    & (a_pos*U_pos + a_neg*U_neg)
);

Hope that clarifies at least some of your doubts.

Cheers,
USV

[mkhm]

Thanks a lot USV. Your explanation clarifies my question very well.