CompressibleInterFoam: theory, implementation and examples

This page is under development

CompressibleInterFoam is a default multiphase solver in OpenFOAM. It is capable of solving two compressible, non-isothermal immiscible fluids, with/without cavitation. The phase interface is captured using volume of fluid (VoF) approach.

Derivation from conservation law

The mass conservation for individual species are

\begin{equation} \label{eq:mass_alpha1} \frac{\partial\left(\rho_1 \alpha_1\right)}{\partial t}+\nabla \cdot\left(\rho_1 \boldsymbol{U} \alpha_1\right)=\dot{m} \end{equation}

\begin{equation} \label{eq:mass_alpha2} \frac{\partial\left(\rho_2 \alpha_2\right)}{\partial t}+\nabla \cdot\left(\rho_2 \boldsymbol{U} \alpha_2\right)=-\dot{m} \end{equation}

where $\alpha_{1}+\alpha_{2} = 1$; $\dot{m}$ is the source term from phase change. If there is no phase change, $\dot{m}$ is zero.

Consider only the volume fraction in the l.h.s., and move other parts to the r.h.s. to be the source term. Eq. \ref{eq:mass_alpha1} and \ref{eq:mass_alpha2} were re-written as

\begin{equation} \label{eq:mass_alpha1_comp} \frac{\partial \alpha_1}{\partial t}+\nabla \cdot\left(\boldsymbol{U} \alpha_1\right)=-\frac{\alpha_1}{\rho_1} \frac{D \rho_1}{D t} + \frac{\dot{m}}{\rho_1} \end{equation}

\begin{equation} \label{eq:mass_alpha2_comp} \frac{\partial \alpha_2}{\partial t}+\nabla \cdot\left(\boldsymbol{U} \alpha_2\right)=-\frac{\alpha_2}{\rho_2} \frac{D \rho_2}{D t} - \frac{\dot{m}}{\rho_2} \end{equation}

The sum of these two equations lead to

\begin{equation} \label{eq:divU} \nabla \cdot\boldsymbol{U}=-\left(\frac{\alpha_1}{\rho_1} \frac{D \rho_1}{D t}+\frac{\alpha_2}{\rho_2} \frac{D \rho_2}{D t}\right)+\dot{m}\left(\frac{1}{\rho_1}-\frac{1}{\rho_2}\right) \end{equation}

The 1st and 2nd term in the r.h.s. are the compressibility and phase change effects.

Insert it to Eq. \ref{eq:mass_alpha1_comp} by adding and substracting $\alpha_{1}\nabla \cdot\boldsymbol{U}$

\begin{equation} \label{eq:mass_alpha1_final} \frac{\partial \alpha_1}{\partial t}+\nabla \cdot\left(\boldsymbol{U} \alpha_1\right)=\alpha_1 \alpha_2 \left(\frac{1}{\rho_1} \frac{D \rho_1}{D t}-\frac{1}{\rho_2} \frac{D \rho_2}{D t}\right) +\alpha_1 \nabla \cdot \boldsymbol{U}+\dot{m}\left(\frac{\alpha_2}{\rho_1}+\frac{\alpha_1}{\rho_2}\right) \end{equation}

The compressibility, $\alpha_1 \nabla \cdot \boldsymbol{U}$ and phase change are source terms for the alpha equation. They are separated into $S_{p}$ and $S_{u}$ in the MULES algorithm. $S_{p}$ is implicit source term, $S_{u}$ is explicit source term.

MULES method

In Eq. \ref{eq:mass_alpha1_final}, it was defined

\begin{equation} \mathrm{dgdt}=-\alpha_1 \alpha_2 \left(\frac{1}{\rho_1} \frac{D \rho_1}{D t}-\frac{1}{\rho_2} \frac{D \rho_2}{D t}\right) \end{equation}

dgdt was calculated in the pEqn, i.e.

dgdt =
(
    alpha1*(p_rghEqnComp2 & p_rgh)
    - alpha2*(p_rghEqnComp1 & p_rgh)
);

In solving Eq. \ref{eq:mass_alpha1_final}, dgdt was separated according to its sign in each cell

if (dgdt[celli] > 0.0)
{
    SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
    SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
}
else if (dgdt[celli] < 0.0)
{
    SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
}

The code means, if dgdt > 0

\[S_p = \text{phaseChangeS}[1] -\frac{\mathrm{dgdt}}{\alpha_{2}}, S_u = \text{phaseChangeS}[0] + \frac{\mathrm{dgdt}}{\alpha_{2}}\]

if dgdt < 0

\[S_p = \text{phaseChangeS}[1] + \frac{dgdt}{\alpha_{1}}, S_u = \text{phaseChangeS}[0]\]

$S_p$, $S_u$ together with other terms were sending to

// In src/twoPhaseModels/twoPhaseMixture/VoF/alphaEqn.H
if (divU.valid())
{
    MULES::explicitSolve
    (
        geometricOneField(), // 1
        alpha1,
        phiCN,  // phi for Euler ddt method, i.e. U
        alphaPhi10, // alpha*U
        Sp(),
        (Su() + divU()*min(alpha1(), scalar(1)))(),
        oneField(),
        zeroField()
    );
}

In MULES::explicitSolve, the time step was calculated, then limiter was applied, last another explicitSolve was called.

// In src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false);
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); 
// rho - geometricOneField()
// psi - alpha1
// phiPsi - alphaPhi10
// Sp - SpRef
// Su - SuRef + alpha1*divU()

In the latter explicitSolve called, $\alpha_1$ for next time step was obtained via

// In src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
psiIf =
(
    rho.oldTime().field()*psi0*rDeltaT
    + Su.field()
    - psiIf
)/(rho.field()*rDeltaT - Sp.field());

This means:

\[\alpha_1(t+\Delta t)=\frac{\alpha_1(t)/\Delta t + \mathrm{Su.field} - \frac{1}{V_P} \sum_f U\alpha_1 S_{f}}{1/\Delta t-\mathrm{Sp.field} }\]

The corresponding equation was

\[\frac{\partial \alpha_1}{\partial t}+\nabla \cdot(\boldsymbol{U} \alpha_1)=\alpha_1 S_{p} + S_{u} + \alpha_{1}\nabla\cdot\boldsymbol{U}\]

(Note: Su.field() = $S_u+\alpha_{1}\nabla\cdot\boldsymbol{U}$, $\alpha_1 S_{p} + S_{u}$ is implemented in alphaSuSp.H in the solver)

Pressure-velocity coupling

The velocity field is firstly estimated from the momentum equation, then the velocity field is inserted to Eq. \ref{eq:divU}.

\[\begin{equation} \nabla \cdot (\mathbf{H}_{\boldsymbol{U}} - \frac{1}{a_p} \nabla p_{\mathrm{rgh}}) = -\left(\frac{\alpha_1}{\rho_1} \frac{D \rho_1}{D t} + \frac{\alpha_2}{\rho_2} \frac{D \rho_2}{D t}\right) + \dot{m}\left(\frac{1}{\rho_1} - \frac{1}{\rho_2}\right) \end{equation}\]

The source code for pressure-possion equation is

fvScalarMatrix p_rghEqnIncomp
(
    fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
    == Sp_rgh
);

solve
(
    p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp
);

where p_rghEqnIncomp means the incompressible part while p_rghEqnComp1 and p_rghEqnComp2 are compressible parts corresponding to the two phases.

fvc::div(phiHbyA): $\nabla \cdot \mathbf{H}_{\boldsymbol{U}}$
fvm::laplacian(rAUf, p_rgh): $\nabla \cdot ( \frac{1}{a_p} \nabla p_{\mathrm{rgh}})$
Sp_rgh: source term from phase change, i.e. $\dot{m}\left(\frac{1}{\rho_1}-\frac{1}{\rho_2}\right)$
p_rghEqnComp1(): $\frac{\alpha_1}{\rho_1} \frac{D \rho_1}{D t}$
p_rghEqnComp2(): $\frac{\alpha_2}{\rho_2} \frac{D \rho_2}{D t}$

Consideration of the compressibility

coming soon

Phase change model

The phase change model was included in the $$ \ alpha$ equation via

Pair<tmp<volScalarField::Internal>> phaseChangeS
(
    phaseChange.Salpha(alpha1)
);

if (phaseChangeS[0].valid())
{
    Su = phaseChangeS[0];
    Sp = phaseChangeS[1];
}

and pressure equation via

fvScalarMatrix Sp_rgh(phaseChange.Sp_rgh(rho, gh, p_rgh));

The Schnerr-Sauer model

Source term of alphaEqn

Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
Foam::twoPhaseChangeModels::SchnerrSauer::mDotAlphal() const
{
    ...
    const volScalarField::Internal pCoeff(this->pCoeff(p));
    ...
    return Pair<tmp<volScalarField::Internal>>
    (
        Cc_*limitedAlpha1*pCoeff*max(p - pSat(), p0_),
        Cv_*(1.0 + alphaNuc() - limitedAlpha1)*pCoeff*min(p - pSat(), p0_)
    );
}
    // Note: 
    // pCoeff = (3*rho1()*rho2())*sqrt(2/(3*rho1()))
    //   *rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));

$\dot{m}=\underbrace{\alpha_{v}\dot{m}_{c}}_{\text{condensation}}+\underbrace{\alpha_{l}\dot{m}_{v}}_{\text{vaporization}}$

\[\begin{aligned} & \dot{m}_{c} = C_{c}\alpha_{l}\frac{3\rho_{l}\rho_{v}}{\rho R_{b} \sqrt{|P-P_{s}|}}\sqrt{\frac{2}{3\rho_{l}}}max(P-P_{s},0) \geq 0 \\ &\dot{m}_{v} = C_{v}(1+\alpha_{Nuc}-\alpha_{l})\frac{3\rho_{l}\rho_{v}}{\rho R_{b} \sqrt{|P-P_{s}|}}\sqrt{\frac{2}{3\rho_{l}}}min(P-P_{s},0) \leq 0\\ &\text{pcoeff} = \frac{3\rho_{l}\rho_{v}}{\rho R_{b} \sqrt{|P-P_{s}|}}\sqrt{\frac{2}{3\rho_{l}}} \end{aligned}\]

Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
Foam::twoPhaseChangeModels::cavitationModel::Salpha
(
    volScalarField& alpha
) const
{
    ...
    const volScalarField::Internal vDotcAlphal(alphalCoeff*mDotAlphal[0]());
    const volScalarField::Internal vDotvAlphal(alphalCoeff*mDotAlphal[1]());
    return Pair<tmp<volScalarField::Internal>> // mdot =  \alpha_{l}S_{p} + S_{u}
    (
        1.0*vDotcAlphal, // Su
        vDotvAlphal - vDotcAlphal // Sp
    );
}

The source term in $\alpha_{l}$ equation, $S_{\alpha}$ was seperated as

\[\begin{aligned} S_{\alpha} &= \alpha_{l,coeff}\dot{m} \\ &= \alpha_{l,coeff}(\alpha_{l}\dot{m}_{v} + \alpha_{v}\dot{m}_{c}) \\ &= \alpha_{l,coeff}[\alpha_{l}(\dot{m}_{v}-\dot{m}_{c}) + \dot{m}_{c}]\\ &= \underbrace{\alpha_{l}\alpha_{l,coeff}(\dot{m}_{v}-\dot{m}_{c})}_{\alpha_{l} \text{Sp} <0} + \underbrace{\alpha_{l,coeff}\dot{m}_{c}}_{\text{Su} >0} \\ &= \alpha_{l}\text{Sp} + \text{Su} \end{aligned}\] \[\alpha_{l,coeff} = \left(\frac{1}{\rho_l}-\alpha_{l}\left(\frac{1}{\rho_l}-\frac{1}{\rho_v}\right)\right) > 0\]

In OpenFOAM, fvm::Sp handles negative source term. fvm::Sp makes the negative source term implicit so it contributes to the diagonal and increase the . This can help convergence when the source term is negative on the rhs (sink term). (ref.: OpenFOAM® Beginner training session)

Source term of pEqn

The mass transfer rate in pEqn, $\dot{m}_{p}$ was defined as

\[\dot{m}_{p} = -\frac{\dot{m}}{P-P_{s}} = \underbrace{-\frac{\alpha_{l}\dot{m}_{v}}{P-P_{s}}}_{\dot{m}_{p,v} \leq 0 } - \underbrace{\frac{(1-\alpha_{l})\dot{m}_{c}}{P-P_{s}}}_{\dot{m}_{p,c} \geq 0 } = \dot{m}_{p,v} - \dot{m}_{p,c}\]

Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
Foam::twoPhaseChangeModels::SchnerrSauer::mDotP() const
{
    const volScalarField::Internal apCoeff(limitedAlpha1*pCoeff);

    return Pair<tmp<volScalarField::Internal>>
    (
        Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff,

        (-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff
    );
}

\[\begin{aligned} & \dot{m}_{p,c} = C_{c}(1-\alpha_{l})*\underbrace{\alpha_{l}\text{pcoeff}}_{\text{apcoeff}}*\text{pos0}(P-P_{s}) =\frac{1-\alpha_{l}}{P-P_{s}}\dot{m}_{c} \geq 0 \\ &\dot{m}_{p,v} = -C_{v}(1+\alpha_{Nuc}-\alpha_{l})*\underbrace{\alpha_{l}\text{pcoeff}}_{\text{apcoeff}}*\text{neg}(P-P_{s}) = -\frac{\alpha_{l}}{P-P_{s}}\dot{m}_{v} \leq 0 \end{aligned}\]

Foam::tmp<Foam::fvScalarMatrix>
Foam::twoPhaseChangeModels::cavitationModel::Sp_rgh
(
    const volScalarField& rho,
    const volScalarField& gh,
    volScalarField& p_rgh
) const
{
    const volScalarField::Internal pCoeff(1.0/rho1() - 1.0/rho2());
    const Pair<tmp<volScalarField::Internal>> mDotP = this->mDotP();

    const volScalarField::Internal vDotcP(pCoeff*mDotP[0]);
    const volScalarField::Internal vDotvP(pCoeff*mDotP[1]);

    return
        (vDotvP - vDotcP)*(pSat() - rho()*gh())
      - fvm::Sp(vDotvP - vDotcP, p_rgh);
}

\[\begin{aligned} S_{P}&= \left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)\dot{m}\\ & = \left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)(\dot{m}_{p,v} - \dot{m}_{p,c})(P_{s}-P)\\ & = \underbrace{\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)(\dot{m}_{p,v} - \dot{m}_{p,c})P_{s}}_{\text{Su}>0}\quad \underbrace{-\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)(\dot{m}_{p,v} - \dot{m}_{p,c})P}_{P*\text{Sp}<0}\\ \end{aligned}\]

(Note: ${1}/{\rho_{l}}-{1}/{\rho_{v}} <0$)

Temperature Equation