interCondensatingEvaporatingFoam: pressure equation mass transfer source term

Summary

I was going through interCondensatingEvaporatingFoam equations derivation including constant phase change model. Considering the mass balance equation of the liquid phase (e.g. \alpha_{1}) and positive mass transfer defined as being from the vapor to the liquid:

\frac{\partial \left ( \alpha_{1}\rho_{1} \right )}{\partial t} + \bigtriangledown\cdot \left ( \alpha_{1}\rho_{1}\vec{U} \right ) = -\dot{m_{12}} + \dot{m_{21}}

with \dot{m_{12}} , \dot{m_{21}}> 0.

By summing over both phases one can derive the statement for total mass balance equation:

\bigtriangledown\cdot \left ( \vec{U} \right ) = \left [ -\dot{m_{12}} + \dot{m_{21}} \right ] \left ( \frac{1}{\rho_{1}} - \frac{1}{\rho_{2}} \right )

Converting to the solver implementation notations, we have:

\dot{m_{v}} = -\dot{m_{12}}

\dot{m_{c}} = \dot{m_{21}}

Keeping in mind that the current implementation returns the mass transfer calculated for each phase with the correct sign, we can re-write the total mass balance equation with OpenFOAM notations as follow:

\bigtriangledown\cdot \left ( \vec{U} \right ) = \left [ -\dot{m_{12}} + \dot{m_{21}} \right ] \left ( \frac{1}{\rho_{1}} - \frac{1}{\rho_{2}} \right ) = \dot{v_{21}} - \dot{v_{12}} = \dot{v_{c}} + \dot{v_{v}}

However, instead of \dot{v_{c}} + \dot{v_{v}} the current pEqn implementation reads:

        fvScalarMatrix p_rghEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, p_rgh)
          ==
            vDotv - vDotc
        );

Can you explain why the solver is using \dot{v_{v}} - \dot{v_{c}}?

It is worth to mention that the mass transfer source terms for alphaEqn and TEqn are following the correct analogy as the sign of the terms are consistent. However, there seems to be a another problem regarding how energy sources are calculated for TEqn which I try to explain below.

Steps to reproduce

1. Run condensatingVessel tutorial with the current configuration
2. Run condensatingVessel tutorial with the coeffE = 0

The result is completely different.

Example case

condensatingVessel

What is the current bug behaviour?

When simulating a condensation only problem, the evaporation coefficient affects the solution.

What is the expected correct behavior?

The evaporation coefficient should not affect the condensation only problems and vice-versa.

Environment information

• OpenFOAM version : develop
• Operating system : Ubuntu
• Hardware info : irrelevant
• Compiler : gcc-7

Possible fixes

The reason for this behaviour is that the constant model defines the evaporation and condensation sources simultaneously. Based on the implemented mass transfer mechanisms, a cell can only have one of the following conditions:

1. no phase change
2. evaporation
3. condensation

However, the current implementation includes phase change energy sources without considering the cell condition:

    const volScalarField Vcoeff
    (
        coeffE_*mixture_.rho1()*limitedAlpha1*L
    );
    const volScalarField Ccoeff
    (
        coeffC_*mixture_.rho2()*limitedAlpha2*L
    );

    TSource =
        fvm::Sp(Vcoeff, T) - Vcoeff*TSat
      - fvm::Sp(Ccoeff, T) + Ccoeff*TSat;

I think the correct implementation should consider the cell condition as follow, e.g. we should not have evaporation energy source when the cell is actually condensing.

    const volScalarField Vcoeff
    (
        coeffE_*mixture_.rho1()*limitedAlpha1*L*pos(T - TSat)
    );
    const volScalarField Ccoeff
    (
        coeffC_*mixture_.rho2()*limitedAlpha2*L*pos(TSat - T)
    );

    TSource =
        fvm::Sp(Vcoeff, T) - Vcoeff*TSat
      - fvm::Sp(Ccoeff, T) + Ccoeff*TSat;