Considering the total enthalpy of the phase k

$H_k=e_k+\frac{1}{2}u_k^2+\frac{P}{\rho_k}$   (13)

The energy equation is written

$\frac{\partial}{\partial t}(\alpha_k \rho_k H_{k})+\frac{\partial}{\partial x_j}(\alpha_k \rho_k H_{k}U_{k,j}) = \frac{\partial}{\partial x_j} (\alpha_k \tau_{k,ij} U_{k,i})-\frac{\partial}{\partial x_j}\alpha_kQ'_{k,j}-\alpha_k\frac{\partial P}{\partial x_i}+\alpha_k \rho_k U_{k,i} g_{i}+\Pi_k+Q_{wall→k}+I_{p→k,i}$   (14)

Where $Q'_k=\lambda_kT_k$ and $\lambda_k$ represents the thermal conductivity.

$Q{wall →k}$ denotes the heat exchange with boundaries and is described by the nucleate boiling model. It takes into account bubble creation and satisfies:

 $\sum_{k=1}^{n phase}Q{wall→k}=Q_{wall}$   (15)

where $Q_{wall}$ is the total imposed heat flux.

Furthermore, $\Pi_k$ represents the bulk interface heat transfer, sum of the interface transfer between phase p and phase k, which complies with conservation relation

$\Pi_k= \sum_{p≠k} \Pi_{p→k}$   (16)

 where

$\Pi_{p→k}+\Pi_{k→p}=0$   (17)

It should be noted that the code does not consider the  terms and $\alpha_k \tau_{k,ij} U_{k,i}$ and $I_{p→k,i}$