© The Author(s), published by EDP Sciences, 2025

1 Introduction

Global warming is one of the major concerns since the beginning of this century. Climate change is primarily caused by human activity, and in particular by the transport sector, which emits large amounts of greenhouse gases. Electric cars are currently designed for emission-free modes of transport. However, electricity cannot be stored in great quantities, and its direct use through battery-powered vehicles has some limitations in terms of autonomy and recharging times. Hydrogen appears to be a viable solution for producing electricity in cars by using fuel cells. Among the different available technologies, the Proton Exchange Membrane Fuel Cell (PEMFC) seems to be the best solution for automotive applications due to its excellent dynamics, high power density, and good efficiency (>0.5). PEMFC converts hydrogen into electrical energy when the hydrogen reacts with the oxygen in the air to produce water and energy (electricity, but also heat). Although promising, PEMFC faces challenges for their large-scale commercialisation. Their durability appears to be the main issue to be solved. Degradations, which can be irreversible, are the origin of the limited lifespan of PEMFC. These degradations occur notably during an imperfect management of water and heat.

Due to their importance in PEMFC models, the balance equations of mass, momentum, and energy are derived in order to compare these equations with the models currently used in PEMFC simulations. We are concerned with two-phase multicomponent flows in porous media. The liquid water is a single component, but the gas phase mixture includes at least hydrogen, oxygen, nitrogen, and water vapour. The porous media involved in a PEMFC are the gas diffusion layer (GDL), the microporous layer (MPL), and the catalyst layer (CL). Sometimes, the bipolar plates (BPs) containing the gas distribution and water evacuation channels are also considered as porous media since these BPs contain a large number of small channels. After a brief literature survey (Sect. 2), Sections 3–5 give a complete derivation of the equations governing two-phase multicomponent flows in porous media. Section 3 presents the mathematical tools used in the paper as well as the local instantaneous balance equations. In Section 4, a volume averaging operator is applied to the balance equations derived in Section 3 in order to obtain their averaged counterparts. In Section 5, we give more usable forms of these equations, and Section 6 presents simplifications currently used in PEMFC models.

2 Brief literature survey on the models of flows in general porous media

When studying two-phase flows in porous media, a volume averaging is generally done on a representative elementary volume (REV) containing the three phases (gas, liquid, and the solid medium) (Fig. 1). In Figure 1, Agl denotes the gas-liquid interface, Ags denotes the gas-solid interface, and Als denotes the liquid-solid interface, respectively.

Figure 1 Scheme of a typical REV containing the three phases.

Gray et al. [1] give the proof of general theorems relating spatial averages of derivatives to the derivatives of spatial averages for a two-phase medium (e.g., a porous structure embedded in a fluid). Slattery [2] considers simplified averaged mass and momentum equations for flows of viscoelastic fluids in porous media. Whitaker [3] studied the diffusion of a passive species concentration c in a porous medium. He shows that the average of the species concentration balance equation:$$\frac{\partial c}{\partial t}+\nabla ·\left(\underline{v}c\right)=\nabla ·\left(D\nabla c\right)$$(1)where $\underline{v}$ is the fluid velocity and D is the molecular diffusion coefficient, is given by the following equation:$$\frac{\partial \bar{c}}{\partial t}+\nabla ·\left(\underline{\bar{v}}\bar{c}+\underline{\psi }\right)=\nabla ·\left(D\left(\nabla \bar{c}+R\underline{\tau }\right)\right)$$(2)where the overbar denotes the spatial averaging operator, and the new terms $\underline{\psi }$ and $R\underline{\tau }$ are defined through the following equations. The vector $\underline{\psi }$ is due to the covariance between the velocity and concentration fluctuations:$$\underline{\psi }\stackrel{\scriptscriptstyle\mathrm{def}}{=}\underline{\overline{v^{\text{'}}}}\overline{c^{\text{'}}}$$(3)where the primed quantities represent fluctuations. R denotes the volumetric interfacial area of the porous medium and the vector $\underline{\tau }$ is linked to the medium tortuosity. The last term in equation (2) is a reduction of the concentration diffusion if the porous medium is tortuous (i.e., the links between the pores are not rectilinear). They are defined by:$$R\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{A_i}{V} \underline{\tau }\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\int }_{{\tilde{A}}_i}^{}\left[c\right]\underline{n} \mathrm{d}a$$(4)where the non-dimensional area element da = dA/A_i and [c] is the jump of c between the fluid and solid phases at the fluid-solid interface. Other equations are written by the author for other quantities like the fluctuations c′ and ${\underline{v}}^{\text{'}}$. The remainder of his paper is a tentative model for the modelling of these unknown terms.

Abriola et al. [4] write several mass balance equations for a four-phase two-component medium (insertion of a pollutant in a soil containing also water and air). Except for the mass balance equations, these authors use a very simplified form of the momentum equations, which is a generalized Darcy’s law for multiphase flows:$${\underline{v}}^{\alpha }=-\frac{K{K}_{\mathrm{r\alpha }}}{\epsilon s_{\alpha }{\mu }_{\alpha }}\left(\nabla P^{\alpha }-{\rho }^{\alpha }\underline{g}\right)$$(5)which gives the velocity of the α-phase ${\underline{v}}^{\alpha }$ as a function of the medium permeability K, the relative permeability K_rα, the porosity ε, the saturation s_α (which is the volumetric fraction of this phase), the viscosity μ_α, the pressure gradient ∇P ^α , and the gravitational acceleration vector $\underline{g}$. The relative permeability is often modelled as a function of the phase saturation.

In their book about convection in porous media, Nield and Bejan [5] synthesize the following mass and momentum equations for a single-phase flow in a porous medium:$$\epsilon \frac{\partial {\rho }_f}{\partial t}+\nabla ·\left({\rho }_f\underline{v}\right)=0$$(6)for the mass balance equation, where ε is the porosity (ratio of the fluid volume divided by the total volume), ρ_f is the fluid density and $\underline{v}$ is the seepage velocity, also known as the superficial velocity. It is not the real velocity but the product of it by the porosity:$$\underline{v}=\epsilon {\underline{v}}_f.$$(7)

The general momentum equation reads:$${\rho }_f\left[\frac{1}{\epsilon }\frac{\partial \underline{v}}{\partial t}+\frac{1}{\epsilon }\nabla ·\left(\frac{\underline{v}\underline{v}}{\epsilon }\right)\right]=-\nabla p-\frac{\mu \underline{v}}{K}-c_F\frac{{\rho }_f}{\sqrt{K}}\left|\underline{v}\right|\underline{v}+\frac{\mu }{\epsilon }{\nabla }^2\underline{v}.$$(8)

The LHS (Left Hand Side) of equation (8) is the inertial terms, and the terms on the RHS (Right Hand Side) are the pressure gradient, the friction term due to Darcy, an additional quadratic friction term due to Forchheimer, and a diffusion term due to Brinkman. equation (8) is rarely used entirely, but parts of it are supposed to balance, the other terms being neglected. We can therefore retrieve the original Darcy equation, or the Brinkman or Forchheimer equations from (8) by neglecting the appropriate terms.

Nield and Bejan [5] also introduce two thermal equations for the fluid and solid phases of the porous medium. This model is sometimes called the Local Thermal Non-Equilibrium (LTNE) since each phase is characterized by its own temperature equation, therefore taking the thermal imbalance between phases into account. The equations given by [5] read:$$\left(1-\epsilon \right){\left(\rho c\right)}_s\frac{\partial T_s}{\partial t}=\left(1-\epsilon \right)\nabla ·\left(k_s\nabla T_s\right)+\left(1-\epsilon \right){q^{\mathrm{\text{'}\text{'}\text{'}}}}_s+ha\left(T_f-T_s\right)$$(9) $$\epsilon {\left(\rho c_p\right)}_f\frac{\partial T_f}{\partial t}+{\left(\rho c_p\right)}_f\underline{v}·\nabla T_f=\epsilon \nabla ·\left(k_f\nabla T_f\right)+\epsilon {q^{\mathrm{\text{'}\text{'}\text{'}}}}_f+ha\left(T_s-T_f\right).$$(10)

In these equations, T_s and T_f denote the average solid and fluid temperatures, respectively. The quantities q′′′_s and q′′′_f denote volumetric source terms in each phase and h and a denote a heat exchange coefficient and the volumetric contact area between phases, respectively. If we assume that these two temperatures are equal:$$T_s= T_f=T \mathrm{thermal} \mathrm{equilibrium} \mathrm{assumption}$$(11)we can deduce a single equation for the mean temperature T by adding together the two equations (9) and (10). We thus obtain:$${\left(\rho c\right)}_m\frac{\partial T}{\partial t}+{\left(\rho c_p\right)}_f\underline{v}·\nabla T=\nabla ·\left(k_m\nabla T\right)+{q^{\mathrm{\text{'}\text{'}\text{'}}}}_m$$(12)where (ρc)_m, k_m, and q′′′_m are the physical properties of the medium and the volumetric source term in the medium, which are defined by the following equations:$${\left(\rho c\right)}_m=\left(1-\epsilon \right){\left(\rho c\right)}_s+\epsilon {\left(\rho c_p\right)}_f$$(13) $$k_m=\left(1-\epsilon \right)k_s+\epsilon k_f$$(14) $${q^{\mathrm{\text{'}\text{'}\text{'}}}}_m=\left(1-\epsilon \right){q^{\mathrm{\text{'}\text{'}\text{'}}}}_s+\epsilon {q^{\mathrm{\text{'}\text{'}\text{'}}}}_f.$$(15)

We insist on the fact that the equations (12)–(15) are only valid under the thermal equilibrium assumption (11) and that the porosity ε should be constant in space, otherwise it cannot be excluded from the divergences as in equations (9) and (10) as we will show later.

References [6–9] derived complete set of equations for two-phase flow modelling. Their equations are very detailed, unfortunately they do not deal with porous media.

Sha et al. [10] started from a result demonstrated by [1]:$$\nabla ·\langle {\underline{\phi }}_k\rangle =\frac{1}{V}{\int }_{A_e}^{}{\underline{\phi }}_k·{\underline{n}}_k\mathrm{d}a$$(16)where ${\underline{\phi }}_k$ is an arbitrary vector characterizing the k phase, $\langle {\underline{\phi }}_k\rangle$ is its average over the REV (Fig. 1) denoted by V and A_e is the part of A (the boundary surface of V) available to the fluid passage (i.e., the portion of A not contained in the solid parts of the medium). The relation demonstrated by [10] reads:$$\nabla ·\langle {\underline{\phi }}_k\rangle =\frac{\partial }{\partial x}\left({\gamma }_{\mathrm{Ax}}{\alpha }_k{\langle {\phi }_{\mathrm{kx}}\rangle }_{\mathrm{kx}}\right)+\frac{\partial }{\partial y}\left({\gamma }_{\mathrm{Ay}}{\alpha }_k{\langle {\phi }_{\mathrm{ky}}\rangle }_{\mathrm{ky}}\right)+\frac{\partial }{\partial z}\left({\gamma }_{\mathrm{Az}}{\alpha }_k{\langle {\phi }_{\mathrm{kz}}\rangle }_{\mathrm{kz}}\right)$$(17)where α_k = V_k/V_m, V_k being the instantaneous volume of phase k contained in the volume V and V_m is the fluid mixture volume contained in the volume V. The three quantities γ_Ax, γ_Ay, and γ_Az are the surface porosities contained in each of the planes yz, xz, and xy, respectively, the volume V assumed by [10] being a rectangular parallelepiped. The quantities γ_Ax, γ_Ay, and γ_Az being directional porosities, the quantities 〈φ_kx〉_kx, 〈φ_ky〉_ky, and 〈φ_kz〉_kz represent the phase-averaged vector components through the surfaces normal to x, y, and z occupied by phase k. The demonstration of the relation (17) assumes that the size of the REV is not arbitrary: it should be large in comparison to the pores, size and small in comparison to the global scale of the flow [11]. According to [11], if conditions of both stability and regularity are satisfied, then the volume-averaged and surface-averaged values coincide at every point, hence there is no need to make the distinction between the surface porosities γ_Ax, γ_Ay, and γ_Az and the volumetric one given by ε. Similarly, there is no need to make the distinction between 〈φ_kx〉_kx and 〈φ_kx〉_k which is the phase-average of the vector component φ_kx. According to these assertions, the relation (17) simplifies to:$$\nabla ·\langle {\underline{\phi }}_k\rangle =\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{\phi }}_k\rangle }_k\right)$$(18)which is simply a consequence of the chosen notations: $\langle {\underline{\phi }}_k\rangle \stackrel{\scriptscriptstyle\mathrm{def}}{=}\epsilon {\alpha }_k{\langle {\underline{\phi }}_k\rangle }_k$ (see later in the article). The proof of what we just said is also given by [12].

3 Two-phase flow equations in porous media

3.1 Definitions and mathematical tools

We consider a two-phase flow in a porous medium. In this porous medium, we isolate a REV containing the three phases: gas, liquid, and solid (Fig. 1). This REV, denoted by V, will be the control volume for averaging. It must be large enough in comparison to the size of the pores and links between them, and small enough in comparison to the global scales of the flow. This double condition on the size of the REV is necessary for the averaged quantities to be stable and regular [11]. It should be remarked that it is not always possible in particular applications, and this condition must be kept in mind when applying the average model to a particular situation, like PEMFC modelling. The porous media in PEMFC (gas diffusion layers, microporous layers, and catalytic layers) contain two fluid phases (water and gas) and several chemical components (oxygen, hydrogen, nitrogen, and water vapour at least), so we develop a multi-constituent two-phase model for porous media. Indexing by k the general phase, by s the solid phase, and by m the liquid-gas mixture in the pores, the porosity is defined by the following relation:$$\epsilon \stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\sum _{k\ne s}^{}{V}_k}{V}=\frac{V_m}{V}$$(19)where V_k denotes the instantaneous volume occupied by phase k (inside the REV) and V_m denotes the fluid volume in the REV, which does not depend on time since the REV is assumed to be fixed and the porous medium is assumed to be undeformable. The phase k volume fraction, sometimes called the saturation in the porous literature, is defined by the following relation:$${\alpha }_k=\frac{V_k}{\sum _{k\ne s}^{}{V}_k}=\frac{V_k}{\epsilon V}.$$(20)

From equation (20), we can deduce:$$\sum _{k\ne s}^{}{\alpha }_k=1.$$(21)

Let us define the Phase Indicator Function (PIF) by the following relation:$$X_k\left(\underline{x},t\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}\left\{\begin{array}{c}1 \mathrm{if} \mathrm{the} \mathrm{point} \underline{x} \mathrm{pertains} \mathrm{to} \mathrm{phase} k \mathrm{at} \mathrm{time} t\\ 0 \mathrm{otherwise}\end{array}\right\}.$$(22)

It can be shown by means of the generalized function theory (e.g., [1, 9]):$$\nabla X_k=-{\underline{n}}_k{\delta }_I$$(23)where ${\underline{n}}_k$ is the unit vector normal to the interfaces with the other phases and pointing outwardly from phase k. The quantity δ_I is a Dirac delta function having different interfaces as its support. δ_I is also called the local instantaneous interfacial area concentration by [13].

The volume average of an arbitrary quantity ψ is defined by the following relation:$$\langle \psi \rangle \stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V}{\int }_V^{}\psi \mathrm{d}v.$$(24)

The phase-averaged of an arbitrary quantity ψ_k pertaining to phase k is defined on the phase volume:$${\langle {\psi }_k\rangle }_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V_k}{\int }_{V_k}^{}{\psi }_k\mathrm{d}v.$$(25)

Because of the PIF definition, it is easily shown that:$$\langle X_k\rangle \stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V}{\int }_V^{}{X}_k\mathrm{d}v=\frac{V_k}{V}=\epsilon {\alpha }_k.$$(26)

The mathematical link between the volume average (24) and the phase average (25) is:$$\langle X_k{\psi }_k\rangle \stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V}{\int }_{V_k}^{}{\psi }_k\mathrm{d}v=\frac{V_k}{V}\frac{1}{V_k}{\int }_{V_k}^{}{\psi }_k\mathrm{d}v=\epsilon {\alpha }_k{\langle {\psi }_k\rangle }_k.$$(27)

Two theorems must be introduced at this stage, which allow the averaging operator to permute with the time and space derivatives. These theorems are sometimes called the limiting forms of the Leibniz and Gauss theorem [6–8]. Here we give the version of these two theorems as derived by Gray and Lee [1]. Multiplying the gradient of an arbitrary quantity by the PIF and volume averaging, we can write:$$\langle X_k\nabla {\psi }_k\rangle =\nabla \langle X_k{\psi }_k\rangle -\langle {\psi }_k\nabla X_k\rangle =\nabla \langle X_k{\psi }_k\rangle +\langle {\psi }_k{\underline{n}}_k{\delta }_I\rangle$$(28)where the last term $\langle {\psi }_k{\underline{n}}_k{\delta }_I\rangle$ represents the contributions of the interfaces between phase k and the other phases. For example, if phase k is the liquid phase, the Dirac distribution δ_I peaks to the interface with the gas on one hand, and to the interface with the solid part on the other hand. The equation (28) relates the average of a gradient to the gradient of the average. Adding a dot (scalar product) in (28), the corresponding equation for the divergence can be easily obtained. Another writing for the last term in (28) is [14]:$$\langle {\psi }_k{\underline{n}}_k{\delta }_I\rangle =\frac{1}{V}{\int }_{A_{I,k}}^{}{\psi }_k{\underline{n}}_k\mathrm{d}a=\frac{1}{V}{\int }_{A_{\mathrm{kf}}}^{}{\psi }_k{\underline{n}}_k\mathrm{d}a+\frac{1}{V}{\int }_{A_{\mathrm{ks}}}^{}{\psi }_k{\underline{n}}_k\mathrm{d}a$$(29)where the integral on the surface area A_I,k has been decomposed into the sum of two integrals: one on the interface between phase k and the other fluid phase f A_kf and another on the interface between phase k and the solid parts A_ks. A_I,k = A_kf ∪ A_ks is the total interfacial area in the volume V seen by the phase k. Putting ψ_k = 1 into the equation (28), we obtain a relation giving the gradient of the product between the porosity and the phase saturation:$$\nabla \left(\epsilon {\alpha }_k\right)=-\frac{1}{V}{\int }_{A_{I,k}}^{}{\underline{n}}_k\mathrm{d}a=-\frac{1}{V}{\int }_{A_{\mathrm{kf}}}^{}{\underline{n}}_k\mathrm{d}a-\frac{1}{V}{\int }_{A_{\mathrm{ks}}}^{}{\underline{n}}_k\mathrm{d}a.$$(30)

The PIF verifies the following transport equation [9]:$$\frac{\partial X_k}{\partial t}+{\underline{v}}_I·\nabla X_k=0$$(31)where ${\underline{v}}_I$ is the local instantaneous interface velocity field. Using equations (23) and (31), we can write:$$\langle X_k\frac{\partial {\psi }_k}{\partial t}\rangle =\frac{\partial }{\partial t}\langle X_k{\psi }_k\rangle -\langle {\psi }_k\frac{\partial X_k}{\partial t}\rangle =\frac{\partial }{\partial t}\langle X_k{\psi }_k\rangle -\langle {\psi }_k{v}_I·n_k{\delta }_I\rangle$$(32)which is the second theorem giving the link between the average of a time derivative and the time derivative of the average. The equation (32) can equivalently be written as:$$\langle X_k\frac{\partial {\psi }_k}{\partial t}\rangle =\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\psi }_k\rangle }_k\right)-\frac{1}{V}{\int }_{A_{\mathrm{kf}}}^{}{\psi }_k{v}_I·n_k\mathrm{d}a.$$(33)

Here, there is no contribution from the interface with the solid because the solid parts are assumed to be immobile. Putting ψ_k = 1 into (33), the time derivative of the phase volumetric fraction is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{kf}}}^{}{v}_I·n_k\mathrm{d}a.$$(34)

3.2 Local instantaneous balance equations for two-phase flows

A general balance equation for two-phase flows before the averaging operation can be written as [6–8]:$$\frac{\partial {\rho }_k{\psi }_k}{\partial t}+\nabla ·\left({\rho }_k{\psi }_k{\underline{v}}_k\right)=-\nabla ·{\underline{J}}_k+{\varphi }_k$$(35)where ρ_k is the phase density, ${\underline{v}}_k$ is the phase velocity, ψ_k is an arbitrary quantity per unit mass of phase k, ${\underline{J}}_k$ is a diffusional flux and ϕ_k is a general source term. The equation (35) is valid inside the phase k only. The corresponding equation valid through the interfaces between two phases (sometimes called the jump condition) reads:$$\sum _k^{}\left({\dot{m}}_k{\psi }_k-{\underline{J}}_k.{\underline{n}}_k\right)={\varphi }_I$$(36)

The quantity ${\dot{m}}_k$ is the mass gain for phase k, per unit surface and unit time, due to phase change (vaporisation and condensation). It is defined by the following equation:$${\dot{m}}_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\rho }_k\left({\underline{v}}_I-{\underline{v}}_k\right)·{\underline{n}}_k.$$(37)

The quantity ϕ_I is a general source term defined on the interfaces. Table 1 gives the values of ψ_k, J_k, ϕ_k and ϕ_I for the different balance equations of interest: mass, species mass, momentum, and the different forms of energy.

Substituting these values into the balance equations (35) and (36), we thus obtain the detailed balance equations for phase k. The mass balances read:$$\frac{\partial {\rho }_k}{\partial t}+\nabla ·\left({\rho }_k{\underline{v}}_k\right)=r_k \mathrm{and} \sum _k^{}{\dot{m}}_k=0$$(38)where ρ_k is the density of phase k and the symbol Σ_k means the summation on the two phases separated by the interface (the second Eq. (38)), as all the jump conditions, is valid on an interface only). The symbol r_k is used for the source term due to chemical reactions. The species balance equations for species i in phase k read:$$\frac{\partial {\rho }_k{Y}_{\mathrm{ik}}}{\partial t}+\nabla ·\left({\rho }_k{Y}_{\mathrm{ik}}{\underline{v}}_k\right)=-\nabla ·{\underline{j}}_{\mathrm{ik}}+r_{\mathrm{ik}} \mathrm{and} \sum _k^{}\left({\dot{m}}_k{Y}_{\mathrm{ik}}-{\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k\right)=r_{\mathrm{is}}$$(39)where Y_ik is the mass fraction of species i in phase k, j _ik is the diffusive flux and r_ik is a source term due to homogeneous chemical reactions. The term r_is is the rate at which species i is produced by heterogeneous chemical reactions.

The momentum balance equation in phase k reads:$$\frac{\partial {\rho }_k{\underline{v}}_k}{\partial t}+\nabla ·\left({\rho }_k{\underline{v}}_k{\underline{v}}_k\right)=-\nabla p_k+\nabla ·{\underline{\underline{\tau }}}_k+{\rho }_k\underline{g} +r_k{v}_k^r.$$(40)

The term $r_k{v}_k^r$ was introduced by [14] and corresponds to a supply of momentum due to chemical reactions, p_k is the pressure for phase k and ${\underline{\underline{\tau }}}_k$ is the viscous stress tensor. The momentum jump condition reads [9, 15]:$$\sum _k^{}\left({\dot{m}}_k{\underline{v}}_k+{\underline{\underline{\sigma }}}_k·{\underline{n}}_k\right)={\underline{F}}_I\cong {\nabla }_s\sigma +2H_k\sigma {\underline{n}}_k.$$(41)

The vector $\underline{g}$ is the gravity acceleration, ${\underline{\underline{\sigma }}}_k$ is the total stress tensor (grouping the pressure and the viscous ones) and ${\underline{F}}_I$ is the surface tension force per unit area. In (41) σ is the surface tension coefficient and $H_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}-\left(\nabla ·{\underline{n}}_k\right)/2$ is the local mean curvature of the interface seen by phase k. The first term ∇_sσ is the surface gradient of the surface tension, which is called the Marangoni effect in the literature.

The total energy balance equation in phase k reads:$$\frac{\partial {\rho }_k\left(e_k+\frac{{v_k}^2}{2}\right)}{\partial t}+\nabla ·\left({\rho }_k\left(e_k+\frac{{v_k}^2}{2}\right){\underline{v}}_k\right)=-\nabla ·\left({\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k\right)+{\rho }_k\underline{g}·{\underline{v}}_k$$(42)and the total energy jump condition reads:$$\sum _k^{}\left({\dot{m}}_k\left(e_k+\frac{{v_k}^2}{2}\right)-\left({\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k\right)·{\underline{n}}_k\right)=E_I$$(43)where e_k is the internal energy per unit mass, ${\underline{q}}_k$ is the heat flux density and E_I is a possible interfacial energy source. References [16, 17] give the exact expression for the interfacial source term E_I:$$E_I={\rho }_I\frac{D_I\left(e_I+{v_I}^2/2\right)}{\mathrm{Dt}}+{\dot{m}}_I\left(e_I+{v_I}^2/2\right)-{\rho }_I{\underline{v}}_I·\underline{g}+{\nabla }_s·{\underline{q}}_I-{\nabla }_s·\left({\sigma \underline{v}}_{\mathrm{It}}\right).$$(44)

In this equation, ρ_I is the interface density per unit surface, e_I is the interfacial internal energy per unit mass and ${\underline{v}}_I$ is the interface velocity, ${\sigma \underline{v}}_{\mathrm{It}}$ is the product of the surface tension coefficient and the tangential interfacial velocity vector. The interfacial mass source term ${\dot{m}}_I$ is defined by:$${\dot{m}}_I\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{D_I{\rho }_I}{\mathrm{Dt}}+{\rho }_I{\nabla }_s·{\underline{v}}_I.$$(45)

In this study, we neglect the interface properties except the surface tension and latent heat, therefore the approximate form of equation (44) is:$${\rho }_I={\underline{q}}_I=0\to E_I\cong -{\nabla }_s·\left({\sigma \underline{v}}_{\mathrm{It}}\right).$$(46)

The advantage of the generic balance equation (35) is that we can apply the averaging operator and theorems once to it to obtain the corresponding generic averaged balance equation. After this is done, the detailed averaged balance equations will be deduced from the generic one by employing the values given in Table 1.

The enthalpy balance equation reads:$$\frac{\partial {\rho }_k{h}_k}{\partial t}+\nabla ·\left({\rho }_k{h}_k{\underline{v}}_k\right)=-\nabla ·{\underline{q}}_k+\frac{D_k{p}_k}{\mathrm{Dt}}+{\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k +r_k{h}_k^r$$(47)where $h_k^r$ is a reaction enthalpy and $\frac{D_k{p}_k}{\mathrm{Dt}}$ is the material derivative of the pressure:$$\frac{D_k{p}_k}{\mathrm{Dt}}\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{\partial p_k}{\partial t}+{\underline{v}}_k·\nabla p_k.$$(48)

The term ${\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k$ is the dissipation of mechanical energy into heat due to internal friction. The interfacial enthalpy balance reads:$$\sum _k^{}\left({\dot{m}}_k{h}_k-{\underline{q}}_k·{\underline{n}}_k\right)=H_I.$$(49)

References [16, 17] derive the following expression for the interfacial balance of enthalpy:$$\sum _k^{}\left({\dot{m}}_k\left(h_k+\frac{{v_k}^2}{2}-{\underline{v}}_k·{\underline{v}}_I\right)-{\underline{q}}_k·{\underline{n}}_k+{\underline{\underline{\tau }}}_k·{\underline{n}}_k·\left({\underline{v}}_k-{\underline{v}}_I\right)\right)=-{\rho }_I\frac{D_I{h}_I}{\mathrm{Dt}}-\frac{D_I\sigma }{\mathrm{Dt}}-{\dot{m}}_I\left(h_I-\frac{{v_I}^2}{2}\right)-{\nabla }_s·{\underline{q}}_I.$$(50)

The comparison of (49) and (50) allows us to find the exact expression of the term H_I. Neglecting the interfacial properties as we have done in (46) and neglecting the mechanical effects in comparison to the thermal effects in the LHS of (50), we obtain the following approximate expression for H_I:$$H_I\cong -\frac{D_I\sigma }{\mathrm{Dt}}=-\frac{\partial \sigma }{\partial t}-{\underline{v}}_I·{\nabla }_s\sigma .$$(51)

The internal energy balance equation reads:$$\frac{\partial {\rho }_k{e}_k}{\partial t}+\nabla ·\left({\rho }_k{e}_k{\underline{v}}_k\right)=-\nabla ·{\underline{q}}_k{-p}_k\nabla ·{\underline{v}}_k+{\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k+r_k{h}_k^r$$(52)where the source term due to homogeneous chemical reactions is proportional to the reaction enthalpy $h_k^r$ [14]. This can be easily verified by comparing (52) and (47). The interfacial balance equation (36) reads:$$\sum _k^{}\left({\dot{m}}_k{e}_k-{\underline{q}}_k·{\underline{n}}_k\right)={ϵ}_I.$$(53)

References [16, 17] derive the following expression for the interfacial balance of internal energy:$$\sum _k^{}\left({\dot{m}}_k\left(e_k+\frac{{v_k}^2}{2}-{\underline{v}}_k·{\underline{v}}_I\right)-{\underline{q}}_k·{\underline{n}}_k+{\underline{\underline{\sigma }}}_k·{\underline{n}}_k·\left({\underline{v}}_k-{\underline{v}}_I\right)\right)=-{\rho }_I\frac{D_I{e}_I}{\mathrm{Dt}}-{\dot{m}}_I\left(e_I-\frac{{v_I}^2}{2}\right)-{\nabla }_s·{\underline{q}}_I+\sigma {\nabla }_s·{\underline{v}}_{\mathrm{It}}.$$(54)

The comparison of (53) and (54) allows to find the exact expression of the term ϵ_I. Neglecting the interfacial properties as we have done before and neglecting the mechanical effects in comparison to the thermal effects in the LHS of (54), we obtain the following approximate expression for ϵ_I:$${ϵ}_I\cong \sigma {\nabla }_s·{\underline{v}}_{\mathrm{It}}.$$(55)

The entropy balance equation reads:$$\frac{\partial {\rho }_k{s}_k}{\partial t}+\nabla ·\left({\rho }_k{s}_k{\underline{v}}_k\right)=-\nabla ·\left(\frac{{\underline{q}}_k}{T_k}\right)+{∆}_k+r_k{s}_k^r.$$(56)where ∆_k is the entropy source term, which is given by:$${∆}_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{{\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k}{T_k}+{\underline{q}}_k·\nabla \left(\frac{1}{T_k}\right)\ge 0.$$(57)

The fact that ∆_k ≥ 0 is the expression of the second law of thermodynamics for phase k.

References [16, 17] give the following expression of the second law of thermodynamics at the interface:$$\sum _k^{}\left({\dot{m}}_k{s}_k-\frac{{\underline{q}}_k}{T_k}·{\underline{n}}_k\right)+{\rho }_I\frac{D_I{s}_I}{\mathrm{Dt}}+{\dot{m}}_I{s}_I+{\nabla }_s·\left(\frac{{\underline{q}}_I}{T_I}\right)\ge 0.$$(58)

Equation (36) gives for the particular case of entropy (see Table 1):$$\sum _k^{}\left({\dot{m}}_k{s}_k-\frac{{\underline{q}}_k}{T_k}·{\underline{n}}_k\right)={∆}_I.$$(59)

Comparing these two expressions and neglecting the interfacial properties as we have done before, the interfacial entropy source is equal to zero:$${∆}_I\cong 0.$$(60)

One important consequence of (60) is that the local temperatures are equal at the interface [7, 8]:$$T_L=T_G=T_I.$$(61)

Equation (47), (52), and (56) are equivalent to the primary total energy equation (42). The choice of the energy equation is just a question of convenience.

4 Two-phase volume-averaged equations

Multiplying the phase general equation (35) by X_k, taking the volume average of the resulting equation and using the theorems and definitions of the previous section, the following averaged general balance equation is easily obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k{\psi }_k\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k{\psi }_k{\underline{v}}_k\rangle }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k{\psi }_k-{\underline{J}}_k·{\underline{n}}_k\right)\mathrm{d}a-\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{J}}_k\rangle }_k\right)+\epsilon {\alpha }_k{\langle {\varphi }_k\rangle }_k.$$(62)

In the same manner, the interfacial average of equation (36) gives:$$\sum _k^{}\frac{1}{V}{\int }_{A_I}^{}\left({\dot{m}}_k{\psi }_k-{\underline{J}}_k·{\underline{n}}_k\right)\mathrm{d}a=\frac{1}{V}{\int }_{A_I}^{}{\varphi }_I\mathrm{d}a.$$(63)

4.1 Mass balance equation

Putting ${\psi }_k=1, {\underline{J}}_k=0, {\varphi }_k=r_k$ in the equation (62), the following phase k mass balance equation is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k{\underline{v}}_k\rangle }_k\right)=\epsilon {\alpha }_k{\langle r_k\rangle }_k+\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}{\dot{m}}_k\mathrm{d}a.$$(64)

Now, two different notations will be introduced. First, the phase average of the product between the density and velocity will be rewritten as:$${\langle {\rho }_k{\underline{v}}_k\rangle }_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}{\langle {\rho }_k\rangle }_k{\stackrel{\text{̿}}{{\underline{v}}_k}}^k$$(65)which defines the Favre averaged velocity ${\stackrel{\text{̿}}{{\underline{v}}_k}}^k$. The Favre average has been introduced in the compressible flow theories in order to avoid the appearance of double correlations between density fluctuations and velocity fluctuations. The mean velocity ${\stackrel{\text{̿}}{{\underline{v}}_k}}^k$ is the center-of-mass velocity of the k phase included in the REV. The last term in the RHS of equation (64) corresponds to the phase k volumetric production rate due to phase change, and we introduce the following special notation for it:$${\mathrm{\Gamma }}_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}{\dot{m}}_k\mathrm{d}a.$$(66)

Multiplying the last equation (38) by δ_I followed by volume averaging, we obtain:$$\sum _k^{}{\mathrm{\Gamma }}_k=0$$(67)showing that there is no mass accumulation on the interfaces.

4.2 Species mass balance equation

Putting ${\psi }_k=Y_{\mathrm{ik}}, {\underline{J}}_k={\underline{j}}_{\mathrm{ik}}, {\varphi }_k=r_{\mathrm{ik}}$ in the equation (62), the averaged species mass balance equation is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k{Y}_{\mathrm{ik}}\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k{Y}_{\mathrm{ik}}{\underline{v}}_k\rangle }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k{Y}_{\mathrm{ik}}-{\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k\right)\mathrm{d}a-\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{j}}_{\mathrm{ik}}\rangle }_k\right)+\epsilon {\alpha }_k{\langle r_{\mathrm{ik}}\rangle }_k.$$(68)

Multiplying the second equation (39) by δ_I followed by volume averaging, the following averaged jump condition is obtained:$$\sum _k^{}\left({\mathrm{\Gamma }}_k{\stackrel{̅}{\overline{Y_{\mathrm{ik}}}}}^I-a_{\mathrm{Ik}}{\langle {\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k\rangle }_I\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}{r}_{\mathrm{is}}\mathrm{d}a.$$(69)where the following averages are defined:$$\langle {\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k{\delta }_I\rangle \stackrel{\scriptscriptstyle\mathrm{def}}{=}{a}_{\mathrm{Ik}}{\langle {\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k\rangle }_I$$(70) $$\langle {\dot{m}}_k{Y}_{\mathrm{ik}}{\delta }_I\rangle ={\mathrm{\Gamma }}_k{\stackrel{̅}{\overline{Y_{\mathrm{ik}}}}}^I$$(71)where a_Ik = 〈δ_I〉 is the interfacial area concentration seen by phase k (including the wall interface), ${\langle {\underline{j}}_{\mathrm{ik}}·{\underline{n}}_k\rangle }_I$ is an interfacial average of the normal diffusive flux and ${\stackrel{̅}{\overline{Y_{\mathrm{ik}}}}}^I$ is the interfacial average of the mass fraction weighted by phase change.

4.3 Momentum balance equation

Putting ${\psi }_k={\underline{v}}_k, {\underline{J}}_k=-{\underline{\underline{\sigma }}}_k, {\varphi }_k={\rho }_k\underline{g}+r_k{\underline{v}}_k^r$ in the equation (62), the averaged momentum balance equation is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k{\underline{v}}_k\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k{\underline{v}}_k{\underline{v}}_k\rangle }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k{\underline{v}}_k+{\underline{\underline{\sigma }}}_k·{\underline{n}}_k\right)\mathrm{d}a+\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{\underline{\sigma }}}_k\rangle }_k\right)+\epsilon {\alpha }_k{\langle {\rho }_k\rangle }_k\underline{g}+\epsilon {\alpha }_k{\langle r_k{\underline{v}}_k^r\rangle }_k.$$(72)

We will introduce a new notation for the set of momentum interfacial terms:$${\underline{M}}_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}\langle \left({\dot{m}}_k{\underline{v}}_k+{\underline{\underline{\sigma }}}_k·{\underline{n}}_k\right){\delta }_I\rangle =\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k{\underline{v}}_k+{\underline{\underline{\sigma }}}_k·{\underline{n}}_k\right)\mathrm{d}a.$$(73)

Multiplying the jump condition (41) by δ_I and averaging, the following averaged jump condition is obtained:$$\sum _k^{}{\underline{M}}_k={\underline{M}}_I\stackrel{\scriptscriptstyle\mathrm{def}}{=}\langle {\underline{F}}_I{\delta }_I\rangle$$(74)where ${\underline{M}}_I$ is the interfacial momentum source due to the surface tension.

4.4 Total energy balance equation

Putting ${\psi }_k=e_k+\frac{{v_k}^2}{2}, {\underline{J}}_k={\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k, {\varphi }_k={\rho }_k\underline{g}·{\underline{v}}_k$ in the equation (62), the averaged total energy balance equation is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k\left(e_k+\frac{{v_k}^2}{2}\right)\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k\left(e_k+\frac{{v_k}^2}{2}\right){\underline{v}}_k\rangle }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k\left(e_k+\frac{{v_k}^2}{2}\right)-\left({\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k\right)·{\underline{n}}_k\right)\mathrm{d}a-\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k\rangle }_k\right)+\epsilon {\alpha }_k{\langle {\rho }_k{\underline{v}}_k\rangle }_k·\underline{g}.$$(75)

Now, we can define:$$Q_k\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k\left(e_k+\frac{{v_k}^2}{2}\right)-\left({\underline{q}}_k-{\underline{\underline{\sigma }}}_k·{\underline{v}}_k\right)·{\underline{n}}_k\right)\mathrm{d}a.$$(76)

Multiplying the jump condition (43) by δ_I and averaging the resulting equation, the averaged form of the total energy jump condition is obtained:$$\sum _k^{}{Q}_k=Q_I\stackrel{\scriptscriptstyle\mathrm{def}}{=}\langle E_I{\delta }_I\rangle .$$(77)

4.5 Enthalpy balance equation

Putting ${\psi }_k=h_k, {\underline{J}}_k={\underline{q}}_k, {\varphi }_k=\frac{D_k{p}_k}{\mathrm{Dt}}+{\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k+r_k{h}_k^r$ in the equation (62), the averaged enthalpy balance equation is obtained:$$\frac{\partial }{\partial t}\left(\epsilon {\alpha }_k{\langle {\rho }_k{h}_k\rangle }_k\right)+\nabla ·\left(\epsilon {\alpha }_k{\langle {\rho }_k{h}_k{\underline{v}}_k\rangle }_k\right)=\frac{1}{V}{\int }_{A_{\mathrm{Ik}}}^{}\left({\dot{m}}_k{h}_k-{\underline{q}}_k·{\underline{n}}_k\right)\mathrm{da}-\nabla ·\left(\epsilon {\alpha }_k{\langle {\underline{q}}_k\rangle }_k\right)+\epsilon {\alpha }_k{\langle \frac{D_k{p}_k}{\mathrm{Dt}}+{\underline{\underline{\tau }}}_k:\nabla {\underline{v}}_k+r_k{h}_k^r\rangle }_k.$$(78)

And the enthalpy mean jump condition reads (Eq. (49)):$$\sum _k^{}\frac{1}{V}{\int }_{A_I}^{}\left({\dot{m}}_k{h}_k-{\underline{q}}_k·{\underline{n}}_k\right)\mathrm{d}a=\frac{1}{V}{\int }_{A_I}^{}{H}_I\mathrm{d}a\cong 0$$(79)where the term involving H_I can be neglected according to [15].