CFD Online Discussion Forums - Source Term due to evaporation in energy transport equation

Hello foamers,

I'm struggling with the implementation of an evaporation model in interFoam. So far I was able to implement a interface reconstruction algorithm, which selects interface cells and their neighboring cells. With this algorithm I'm able to calculate the local evaporating mass flux:

$j_{int} = \frac{\dot{q}_l+\dot{q}_v}{\Delta h_v} \quad \text{ with } \dot{q}_l=\kappa_l\frac{\text{d}T_l}{\text{d}x_{int,n}}\quad \text{ and } \quad\dot{q}_v=-\kappa_v\frac{\text{d}T_v}{\text{d}x_{int,n}}$

The interface reconstruction provides the geometric data needed to calculate the temperature gradient of the liquid and vapor side. The calculation of the temperature gradient is not performed directly in the cells that contain a part of the interface but in the neighboring cells that contain pure liquid or gas.

$\frac{\text{d}T_l}{\text{d}x_{int,n}} = \frac{T_{c,l}-T_{sat}}{d_{int}},$

$\frac{\text{d}T_v}{\text{d}x_{int,n}} = \frac{T_{c,v}-T_{sat}}{d_{int}},$
where d_int is the distance from the interface to the centre of the neighboring cell.

The mass flux (mDot) through the liquid-vapor interface can be transferred to a energy source term

$S = -\dot{m}^{'''}\Delta h_v \quad \text{with}\quad \dot{m}^{'''}= j_{int}\frac{A_{cell,int}}{V_{cell}},$
which is included in the energy transport equation:

$\frac{\partial\left(\rho c T\right)}{\partial t}+\nabla\cdot\left(\rho c T \vec{u}\right)-\nabla\cdot\left(\kappa\nabla T\right) = S$

The source term for the continuity equation accounts for the divergence of the velocity field in cells where mass is removed or added:

$\nabla\cdot\vec{u}=\frac{\dot{m}^{'''}}{\rho}$
Due to the fact that the velocity field is not free of divergence in cells which contain source terms, an additional source term is added to the transport equation:

$\frac{\text{d}F}{\text{d}t} \nabla \cdot \left(\vec{u}F\right)+\nabla\cdot\left(c_F\vec{u}\vec{n}\left[ F \left( F-1\right)\right]F\right) = \frac{\dot{m}^{'''}}{\rho}F$

The source terms in the continuity and transport equation are smeared over a band of cells to avoid numerical errors as Hardt&Wondra mentioned in their paper:
http://www.sciencedirect.com/science...21999108001228

The implementation of the source terms in the transport and continuity equation seems to be qualitatively correct. I've tested the code with the stefan problem test case. The liquid is evaporating, if I turn off the source term S in the energy equation. But if source term S is turned on, the temperature of the interface is oscillating, which causes an oscillating evaporation flux. Depending on the time step, the fluid is evaporating or condensating.

The corresponding code for solving the energy equation is:

Code:

    fvScalarMatrix TEqn

    (

        fvm::ddt(rhoCp,T)

          + fvm::div(rhoCpPhi,T)

          - fvm::laplacian(k12,T)

        ==

           S //sharp explicit source term (only in interface cells)





    );

I assume the source of the oscillations is the explicit source term. How can I solve the energy equation so that the interface temperature is at saturation temperature?