Note: This documentation is a short, incomplete introduction to the modelling of a homogeneous equilibrium model in TrioCFD; being a work in progress, it will be revised.

The Homogeneous Equilibrium Model (HEM) is an approximation of a multiphase flow where the phases are mechanically and thermally coupled: they move at the same velocity and have the same temperature (which may not be their saturation temperature if the fluids are not of the same type). It is based on the Pb_Multiphase problem and uses instantaneous relaxation of a fully unequilibrated model rather than a direct approach as in [69].

Homogeneous evanescence for mixture modelling

The Pb_Multiphase framework can model a mixture by taking advantage of the evanescence operator. By choosing alpha_res = 1 and alpha_res_min = 0.5, the system reduces to the three partial differential equations of a mixture. Without any further model, this is the Homogeneous Equilibrium Model (the two phases are dynamically locked and at the same temperature); it can be the starting point of a homogeneous relaxation model. These values simplify the evanescence momentum relation to:

\[\begin{pmatrix}\mathcal{Q}_{\text{pred}}\\ \mathcal{Q}_{\text{mino}}\end{pmatrix}\Longrightarrow\begin{pmatrix}\mathcal{Q}_{\text{pred}}+\mathcal{Q}_{\text{mino}}\\ v_{\text{mino}}=v_{\text{pred}}+v_{\text{drift}}\end{pmatrix} \]

The non-viscous system of equations

A non-viscous two-phase flow can be described by Euler's equations, where every thermodynamic variable describes the mixture as the sum of each phase's contributions:

\[\begin{aligned} &\partial_t \rho +\partial_x (\rho u) = 0,\\ &\partial_t (\rho u) +\partial_x (\rho u^2 + P) = 0,\\ &\partial_t (\rho E) +\partial_x (\rho uH) = 0, \end{aligned} \]

with \(E=e+\frac{u^2}{2}\) the total energy and \(H=h+\frac{u^2}{2}\) the total enthalpy. The internal energy and enthalpy of the mixture depend on the volume fraction and the properties of each phase: \(\rho e=\alpha_G\rho_G e_G+(1-\alpha_G)\rho_L e_L\), \(\rho h=\alpha_G\rho_G h_G+(1-\alpha_G)\rho_L h_L\), and the mixture density is \(\rho=\alpha_G\rho_G+(1-\alpha_G)\rho_L\), with \(\alpha_G\in[0,1]\) the gas void fraction. The simplifications \(P_G=P_L\), \(T_G=T_L\), \(v_G=v_L\) reduce the number of equations; to close the system, an equation of state (EOS) for each phase is needed plus a thermodynamic relationship for the mixture.

The equation of state

The following EOS choices are available in Pb_HEM:

Incompressible: both phase densities are constant.
Quasi-compressible: densities depend only on temperature. For water between 50 °C and 263 °C [147] : \(\rho(T)=-0.0023864182\,T^2-0.22507878\,T+1007.1165\).
Compressible (stiffened gas): guarantees the hyperbolicity of the system (the squared sound velocity is always positive): \(P(\rho,e)=(\gamma-1)\rho(e-q)-\gamma P_{\infty}\), \(P(\rho,T)=\rho(\gamma-1)C_vT-\gamma P_{\infty}\), \(T(\rho,h)=\frac{h-q}{C_p}\), \(c^2=\gamma\frac{P+P_{\infty}}{\rho}\), with \(\gamma=C_p/C_v\).
More complex EOS will be available via an external library (e.g. RefProp) through an in-house code.

The Navier-Stokes system of equations

Source terms are added to the Euler system, \(\partial_t W+\nabla(F_c-F_v)=S\), with conservative variables \(W=(\rho,\ \rho U,\ \rho E,\ \varepsilon_{turb})^T\), convective flux \(F_c=(\rho U,\ \rho V\otimes U,\ (\rho E+P)U,\ \varepsilon_{turb}U)^T\), and viscous flux \(F_v=(0,\ \bar{\bar{\tau^v}}+\bar{\bar{\tau^t}},\ (\bar{\bar{\tau^v}}+\bar{\bar{\tau^t}})\cdot U-Q^v-Q^t,\ S_{turb})^T\), \(S\) the source terms and \(\varepsilon_{turb}\) the turbulence-model variables. The viscous stress tensor follows the Boussinesq hypothesis, and several turbulence models can be used.

Drift-flux and drift-velocity models

The drift-flux model is appropriate when the two-phase flow is strongly coupled with only a small kinematic disequilibrium [77]. Zuber and Findlay [192] expressed the void-fraction-weighted mean gas velocity as \(\overrightarrow{V_g}=\frac{\langle\alpha_g\overrightarrow{v_g}\rangle}{\langle\alpha_g\rangle}=C_0\overrightarrow{J}+\overrightarrow{V_{g0}}\), with \(C_0\) the distribution parameter, \(\overrightarrow{V_{g0}}\) the mean drift velocity and \(\overrightarrow{J}=\langle\overrightarrow{j}\rangle\) the area-averaged mixture velocity ( \(\overrightarrow{j}=\alpha_g\overrightarrow{v_g}+\alpha_l\overrightarrow{v_l}\)). The constitutive equations, originally for vertical upward flows [192], [82], were extended to 3D [114], [113].

The drift velocity is implemented in Vitesse_derive_base (filling the relative-velocity table vr and its derivative dvr). Available models:

Model	Used	Validated	Test case
Constant bubble	Yes	Yes	TRUST: canal bouillant drift & TrioCFD: Drift flux
Ishii-Hibiki	Yes	Yes	TRUST: canal bouillant drift & TrioCFD: Drift flux
Spelt	Yes	No
Forces	Yes	Yes	TRUST: canal bouillant drift & TrioCFD: Drift flux

Constant [81] (Vitesse_derive_constante, params C0, vg0_x/y/z): user-defined \(C_0\) and drift velocity vector.
Ishii-Hibiki (bubbly flow) [76] (Vitesse_derive_Ishii; defaults \(C_\infty=1.2\), \(\theta=1.75\), \(\zeta=18\)): \(C_0=(C_\infty+(1-C_\infty)\sqrt{\frac{\rho_g}{\rho_l}})(1-\texttt{sb}\,e^{-\zeta\alpha_g})\), \(\vec{V_{g0}}=-\sqrt{2}(\frac{(\rho_l-\rho_g)g\sigma}{\rho_l^2})^{1/4}(1-\alpha_g)^{\theta}\frac{\vec{g}}{|g|}\).
Spelt-Biesheuvel [163] : \(\vec{V_{g0}}=(-(u_g-u_l)\frac{\vec{g}}{|g|}+\frac{\nu_{\text{Spelt}}\nabla\alpha_g}{\max(\alpha_g,10^{-4})})(1-C_0\alpha_g)\).
Forces (Vitesse_derive_Forces): the drift velocity is derived from a balance of the dispersion and lift forces against the drag.

At the component scale (e.g. nuclear core, steam generator), a porous-medium equivalent is used and the drift parameters \(C_0\), \(\overrightarrow{V_{g0}}\) are either user-defined constants or fitted (Hibiki-Ishii). At the local scale, the 3D drift velocity is derived (Manninen [114]) from a zero force balance on the dispersed phase (lift + buoyancy + dispersion + drag), currently for a dispersed gas phase in a continuous liquid phase in a vertical pipe (centrifugal acceleration neglected).

Two-phase frictional multiplier

The frictional pressure drop in gas-liquid flow is expressed through a two-phase friction multiplier [51], \(\Phi_{k}^2=\frac{(\partial p/\partial z)|_{\text{two-phase}}}{(\partial p/\partial z)|_{\text{single phase }k}}\), based on the pure-liquid friction \(f_l\) and pure-gas friction \(f_g\). The base class (Multiplicateur_diphasique_base) fills a coeff table (multipliers for the single-phase and mix friction factors). Available models:

Homogeneous: \(\Phi^2=1+x(\frac{\rho_l}{\rho_g}-1)\), with \(x\) the gas mass quality.
Friedel [59] (horizontal/vertical smooth tubes, \(\mu_l/\mu_g<1000\)): \(\Phi^2=E+\frac{3.24FH}{Fr^{0.0454}We^{0.035}}\), with \(E=(1-x)^2+x^2\frac{\rho_l f_g}{\rho_g f_l}\), \(F=x^{0.78}(1-x)^{0.224}\), \(H=(\frac{\rho_l}{\rho_g})^{0.91}(\frac{\mu_g}{\mu_l})^{0.19}(1-\frac{\mu_g}{\mu_l})^{0.7}\), and the mixture Froude and Weber numbers.
Lottes-Flinn [110] : sodium two-phase pressure drop.
Müller-Steinhagen [122] : air-water, water-hydrocarbons and refrigerants in pipes, using \(f_m^*=(\frac{f_l}{\rho_l}+a x^b(\frac{f_g}{\rho_g}-\frac{f_l}{\rho_l}))(1-x)^{1/c}+\frac{f_g}{\rho_g}x^c\).

(The exact coeff matrix-filling and the void-fraction blending fractions Frag_l/Frag_g for each multiplier are given in the source / implementation.)