Source Term due to evaporation in energy transport equation

styleworker · February 5, 2014, 13:05

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?

Mahdi2010 · March 24, 2014, 04:16

Hi
Would you tell me how you implement the transport and continuity source terms?
I think you added the source term into pEqu and then use the implicit MULE algorithm for
VOF transport equation. If so, I think it is important to check the order solver reads.
I mean the heat transfer source term is in Tequ while it might be called before pEqu.

jameswilson620 · September 3, 2014, 22:24

Would you please post your code for your pEqn and alphaEqn? I am working on this problem and am having similar difficulties.

Thanks, James

saicharan662000@gmail.com · September 7, 2022, 04:09

Hi style worker ,
How did you get d_int?
Thanks in Advance

February 5, 2014, 13:05	Source Term due to evaporation in energy transport equation	#1
styleworker Member Join Date: May 2012 Posts: 55 Rep Power: 15	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? tom_flint2012 and Kummi like this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
add source term to energy equation by udf	tomas	Fluent UDF and Scheme Programming	7	September 10, 2013 03:08
[swak4Foam] build problem swak4Foam OF 2.2.0	mcathela	OpenFOAM Community Contributions	14	April 23, 2013 14:59
[swak4Foam] funkySetFields compilation error	tayo	OpenFOAM Community Contributions	39	December 3, 2012 06:18
Is it possible to use divergence in the source term of an scalar transport equation?	jannnesss	CFX	0	January 8, 2010 20:53
UDFs for Scalar Eqn - Fluid/Solid HT	Greg Perkins	FLUENT	0	October 11, 2000 04:43

March 24, 2014, 04:16		#2
Mahdi2010 Member Mahdi Join Date: Jul 2012 Posts: 53 Rep Power: 14	Hi Would you tell me how you implement the transport and continuity source terms? I think you added the source term into pEqu and then use the implicit MULE algorithm for VOF transport equation. If so, I think it is important to check the order solver reads. I mean the heat transfer source term is in Tequ while it might be called before pEqu.

September 3, 2014, 22:24		#3
jameswilson620 Member james wilson Join Date: Aug 2014 Location: Orlando, Fl Posts: 39 Rep Power: 12	Would you please post your code for your pEqn and alphaEqn? I am working on this problem and am having similar difficulties. Thanks, James

September 7, 2022, 04:09		#4
saicharan662000@gmail.com Member hari charan Join Date: Sep 2021 Location: India,hyderabad Posts: 97 Rep Power: 5	Hi style worker , How did you get d_int? Thanks in Advance