# Numerical simulation of incompressible two-phase flows with phase change process in porous media

Aroua Ghedira, Zied Lataoui, Adel Benselama, Yves Bertin, Abdelmajid

Jemni

<figure><img src="image_1.png" /></figure>

## To cite this version:

Aroua Ghedira, Zied Lataoui, Adel Benselama, Yves Bertin, Abdelmajid Jemni. Numerical simulation of incompressible two-phase flows with phase change process in porous media. Results in engineering, 2025, 25, pp.103706. 10.1016/j.rineng.2024.103706. hal-04963740

HAL Id: hal-04963740

https://isae-ensma.hal.science/hal-04963740v1

Submitted on 24 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulation of incompressible two-phase flows with phase change process in porous media

Aroua Ghedira<sup>a,*</sup> <sup>□</sup>, Zied Lataoui<sup>a</sup>, Adel M. Benselama<sup>b</sup>, Yves Bertin<sup>b</sup>, Abdelmajid Jemni<sup>a</sup>

## ARTICLE INFO

### Keywords:

- Phase change

- Heat and mass transfer

- Porous media

- OpenFOAM

- Incompressible flow

- Transpiration cooling

## ABSTRACT

A new OpenFOAM solver is developed to simulate fluid flow in porous media with phase-change. The solver is based on equations describing the two-phase flow. The effects of porous media are expressly handled within the framework of the Representative Elementary Volume analysis. For the moving interface experiencing phase-change, the volume-averaged approach is used mainly for reasons of computational efficiency. Among enhancements, the development of the volume-averaged porous flow model adds more porous resistance factors to the momentum equation. The new solver is validated with several test cases. A great agreement is observed between the current study results and the analytical solutions and/or results of previous well documented studies. Numerical simulation of transpiration cooling in porous media proves the efficiency of the new solver in more complex applications on heat and mass transfer dynamics in porous media. The obtained results showed good agreement with previous experimental investigation data (relative errors < 6.5 %). The effect of mass flow rate on cooling efficiency is also discussed.

The present approach can be applied thoroughly, for design guidance purposes, to other problems involving more physical phenomena (combustion, chemical reactors, etc.) and complex forms of force (electrical, magnetic, inertial, etc.).

## 1. Introduction

The phenomenon of phase change, particularly in porous media, has significant importance and has garnered considerable attention from researchers across several fields [1,2]. This includes heat pipes, geothermal systems, and aircraft thermal removal applications [3-7]. Furthermore, the transpiration cooling process including phase change in porous media is one of the most effective ways to dissipate heat in porous material. This cutting-edge cooling technique can be applied to combustors [8], gas turbine blades [9,10], rocket nozzles [11-13], and other engine parts to eliminate excessive heat flux. Despite this fact, it remains a challenging task to predict exactly the distribution of the phase change and heat transfer performance in porous media because of the complexity of this process. Although many numerical approaches exist to solve gas-liquid configurations [14-16], there is still a need for efficient simulation methods. The primary reason for this situation is the limited understanding and analysis of the dynamic characteristics of two-phase flow with phase change in porous media. This point becomes particularly significant due to the fluid flow within a porous channel,

which facilitates the transfer of mass, momentum, and energy between the solid, liquid, and vapor phases. A substantial amount of latent heat accompanies this process. From pore scale and over, capillary force effects are expected to be of paramount importance in particular near wall and contact line regions. Boiling phenomena, which are not limited to industrial scale devices, are gaining further interest. It is involved in significant enhancement of the heat transfer coefficients [17].

In this topic, many experimental and numerical studies have been carried out. To evaluate the impact of the coolant injection and the particle diameter of the sintered porous media on the capacity of transpiration cooling, Huang et al. [7] conducted an experimental study of the phase change in a wind tunnel. With a particular focus on upscaling, Das and Hassanizadeh [18] have investigated mathematical and experimental methodologies to analyze multiphase fluid flow in porous media. Mottet et al. [19] conducted a heat and mass transfer study with phase change in a capillary evaporator using pore-network model and compared the obtained data with experimental results. In their experimental transpiration cooling flow through a porous medium, Ma et al. concluded that the cooling efficiency can be improved by increasing the rate of the cooling water injection [20]. They also discussed the effect of

<footnote><sup>a</sup> Laboratory of Thermal and Energetic Systems Studies (LESTE), National Engineering School of Monastir, University of Monastir, LR99ES31, Monastir 5000, Tunisia</footnote>

<footnote><sup>b</sup> Fluids-Thermal-Combustion Department, ENSMA, Institute PPRIME (UPR CNRS 3346), 1 av. Clement Ader-BP40109, Futuroscope, Chasseneuil 86961, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: aroua.ghedira@gmail.com (A. Ghedira).</footnote>
Nomenclature

C<sub>p</sub> specific heat capacity at constant pressure [kJ/kg K]

C<sub>f</sub> geometric function

d<sub>p</sub> typical pore diameter.

g gravity [m/s²]

K permeability [m²]

k thermal conductivity [W/m K]

L<sub>hg</sub> latent heat of phase change [J/kg]

ṁ mass flow rate [kg/m²s]

ṁ'' mass flow rate per unit volume [kg/m³ s]

p pressure [Pa]

R<sub>p</sub> radius of curvature [m]

T temperature [K]

V velocity [m/s]

U<sub>c</sub> compression velocity [m/s]

z axial position (m)

Greek symbols

α liquid volume fraction

ε porosity

η cooling efficiency

ρ density [kg/m³]

κ interface curvature [m⁻¹]

λ thermal diffusivity [m²/s]

μ dynamic viscosity [kg/m s]

ψ porous medium identifier

σ surface tension [N/m]

Subscripts

i ith phase in multiphase system

n number of phases

l liquid phase

f fluid

eff effective

sat saturation

v vapor phase

w wall

sm solid matrix

pore plugging on the transpiration cooling performance of perforated flat plates. At low injection ratios, the overall cooling performance was improved with the plugging ratio increase. The opposite trend was obtained for high injection ratios.

When studying phase distribution within porous media under various wetting conditions, Luis et al. [21] described a diffuse interface model. Yortsos [22] examined several elements of phase transition in porous media and placed a special emphasis on the impact of confinement in micropores and the microstructure of the porous medium on heat transfer. Whitaker [23] presented a general theory to provide a starting point for simultaneous mass, heat and mass transfer in porous media. They made some restrictive assumptions on the gas-liquid topology in the drying process. In order to investigate coupled heat, air, and moisture transmission across a multilayered porous medium, Tariku [24] created a transient model. More recently, He et al. [25] developed a numerical model with the local thermal non-equilibrium assumption to simulate two-phase flow through a double-layer porous plate. The results revealed that it could reduce the injection pressure and decrease the surface temperature improving the transpiration cooling efficiency. Otherwise, many advances have been recently made to investigate and enhance the heat transfer through porous media, for specific applications, taking advantage of the transpiration cooling process [26,27].

Wilson et al. [28] presented a numerical model and applied it to two-phase flows with different phase change rates under an interfacial superheat. The results of their analysis are in good agreement compared with the porous two-phase Stefan problem. Li et al. [29] studied the phase change in a 3D porous media of a fluid flow as a continuation of a previous study [30] in a 2D domain with a modified Kirchhoff method.

Besides, OpenFOAM, also known as Open Field Operation And Manipulation, is an open-source software package. It is designed to provide numerical solutions for problems related to Continuum Mechanics. The approach relies on the Finite Volume Method (FVM) and is particularly well-suited for computational fluid dynamics (CFD) [31]. Despite multiple uses of OpenFOAM in CFD research, its application in phase-change in the porous domain still lacks some of the most significant features. In addition, the Volume Of Fluid (VOF) approach, based on the Eulerian method, is suitable for two-phase flow simulation due to its relative simplicity and inherent ability to handle many interfaces and topology changes [32].

In their work, C. Ierardi et al. [33] introduced a modified version of the FSF (Filtered Surface Force) method. It is based on a well-consolidated VOF (Volume Of Fluid) CSF (Continuum Surface

Force) solver in OpenFOAM. As a simplification of the initial implementation, a new formulation of the filter applied to the capillary forces was developed and validated with a number of test cases.

Duong A. Hoang et al. [34] conducted a detailed investigation on the effectiveness of the Volume Of Fluid (VOF) approach, as it is implemented in OpenFOAM, in modeling the flow of confined bubbles and droplets, or “segmented flows,” in microfluidics. They examined and suggested the best computational parameters to achieve modest parasitic currents and a crisp fluid interface.

An alternative model for simulating gas-liquid systems was validated by Tronci et al. [35], which was implemented as a diffusion mixture model in OpenFOAM. This approach treats the gas and liquid components as a single phase. Samkhaniani and Ansari [36] described an implementation of the Hardt evaporation model in OpenFOAM. It was shown to be appropriate for microscale phase transition and when the interface temperature deviates from saturation temperature.

For modeling isothermal interfacial flows, the combination of the interface capturing method with single-fluid formalism has been widely applied [37,38]. The method of solving the interface can be split into two major categories: diffuse interface method and the sharp interface method [37]. For both there exist different interface resolving systems such as (FT) front tracking [39], level set (LS) [40] and Volume Of Fluid (VOF) [41]. The VOF method conserves the mass much better than LS or FT, but the interface curvature or normal vector are evaluated less accurately than the other methods [37].

Multiphase flow in porous media models can be found in a large variety of computational fluid dynamics (CFD) packages. These models are commonly based on the representative volume element approach where porosity and permeability are used to mimic the porous structure instead of the very detailed geometry.

Horgue et al. [42] developed a toolbox based on the Implicit Pressure Explicit Saturation (IMPES) algorithm, which is a segregated resolution of the coupled equations, proposed originally by Sheldon et al. [43], to simulate multiphase fluid flows in porous medium using OpenFOAM.

Various open-source simulations for porous medium flows have been created during the last decade, including Dumux [44], PFlotran [45], OpenGeoSys [46], and MRST [47]. The open-source software utilized in this study, OpenFOAM [48,49], was not included in this list because it was not created as a specialized simulator but as a general library used for solving partial differential equations.

However, OpenFOAM maturity and robustness are continuously increasing, so that its use in simulating flow through porous media can
be performed through much less additional development efforts [50, 51]. Since the porous media problem is quite general, it may be applied to facilitate the creation of novel OpenFOAM solvers. Such applications include the analysis of compressible reacting flows [52], the investigation of mass transfer in fuel cells made with solid oxide [53], and the study of the interaction between waves and coastal porous structures [54,55]. Therefore, the main purpose of this research is to develop a model using the library OpenFOAM in order to simulate two-phase fluid flows involving phase change and free interface motion capturing in porous media. It is a basic step to simulate heat and mass transfer in more complicated devices such as heat pipes [56–58].

Despite OpenFoam's many uses in CFD research, its application to phase-change in a porous domain still lacks some of the most complicated considerations. Rare modeling studies that describe the two-phase flow in the porous structure have been published in the literature using OpenFOAM.

Hence, in the present study, a new numerical solver for simulating the two-phase flow with phase change in porous media is developed. The governing equations are presented in the first part. Then, the novel OpenFOAM solver, named porousPhaseChangeFoam, is built in order to simulate heat transmission and phase change in the porous region. We propose a number of numerical tests to validate the new solver. In addition, an altogether unique phase change model is valid and uniformly used within the two-phase flow. Typical case tests are considered to validate the developed solver, for various conditions involving two-phase flow and phase change. In particular, it is used for the simulation of the cooling transpiration process. The impact of the mass rate of the coolant flow on the temperature distribution and the liquid fraction is then investigated.

## 2. Mathematical model

### 2.1. Governing equations

We present the governing equations used for the simulation of the two-phase flow and phase change heat transfer in porous media. These equations are then solved using the finite volume method (FVM) afterwards [59,60].

As mentioned earlier, Eulerian, methods are more suitable to describe two-phase flows experiencing complex and changing topology. The Level-Set (LS) [40] and volume of fluid (VOF) [32] are the two primary illustrative examples of the Eulerian approach. In the current study, only the Fluid method's volume is taken into account.

The main assumptions [60,61] are:

- both of liquid and vapor phases are considered as incompressible, Newtonian and immiscible.

- the porous media is considered homogeneous and isotropic with a uniform porosity and uniform permeability magnitude,

- all thermo-physical properties are assumed as constant

- the liquid and vapor phases are linked at the liquid-vapor interface,

- local thermal equilibrium exists between fluid and porous material.

The volume fraction of the ith phase or also named saturation, $\alpha_i$ is the ratio of the element volume of ith phase ($\Delta v_i$) to the total elemental volume of all phases ($\Delta v$):

$$ \alpha_i = \frac{\Delta v_i}{\Delta v} \qquad (1) $$

The volume fraction of all the phases (n is the phases number) must be equal to unity:

$$ \sum_{i=1}^{n} \alpha_i = 1 \qquad (2) $$

$$ \alpha_i = \begin{cases} 0, & \text{in a gas region} \\ 0; 1[, & \text{unsaturated, i.e. partly liquid - filled region} \\ 1, & \text{in the liquid region} \end{cases} $$

A solid matrix having several linked voids, or pores, is referred to as a porous media. Because the voids are interconnected, they contain one or more fluids that can travel across the medium. The term “porosity” $\epsilon'$ is most usually used to describe the solid matrix’s void fraction.

$$ \epsilon = \begin{cases} 0; 1[, & \text{in the porous region} \\ 1, & \text{in the solid - free region} \end{cases} $$

Along with the solid matrix, the governing equations must be stated for liquid and vapor. Applications involving heat and mass transfer heavily rely on the interactions between solids, liquids, and vapors. Therefore, the governing equation for the liquid and the vapor phases should take into account any potential phase changes between liquid and vapor. The single component's transport equations that coexist in both the liquid, l, and vapor, v, phases are provided [60].

### Continuity equation

In the porous media, both phases (named here vapor/liquid) are considered. The continuity equation is given by:

$$ \bullet \text{ Vapor} \quad \frac{\partial}{\partial t}(\epsilon \rho_v) + \nabla(\epsilon \rho_v V_v) = \dot{m}_{lv}'' \qquad (3) $$

$$ \bullet \text{ Liquid} \quad \frac{\partial}{\partial t}(\epsilon \rho_l) + \nabla(\epsilon \rho_l V_l) = -\dot{m}_{lv}''' \qquad (4) $$

For incompressible flow, these equations reduce to:

$$ \nabla \cdot V = \frac{\dot{m}_{lv}''}{\epsilon} \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right) \qquad (5) $$

Where $\dot{m}_{lv}''$ is the volumetric mass transfer term (kg/m³s). When it is positive and when it is negative, it shows the rate of evaporation and condensation, respectively and V is the pondered velocity:

$$ V = \epsilon \alpha_l V_l + \epsilon(1 - \alpha_l) V_v $$

### Momentum equation

At the macroscopic scale, the momentum equation in porous media takes into consideration the dispersion effect of the amount of movement on the flow (Brinkman term) and the inertial drag component (Forchheimer term). It is given by the following equation: (see Appendix A)

$$ \rho_f \left[ \frac{1}{\epsilon} \frac{\partial V}{\partial t} + \frac{V}{\epsilon^2} \nabla V \right] - \frac{\mu_f}{\epsilon} \nabla^2 V = - \nabla p + \rho_f g - \psi \left[ \frac{\mu_f}{K} V + \frac{\rho_f C_f}{\sqrt{K}} |V| V \right] \\ - \psi \frac{2\sigma}{R_p} \nabla \alpha_l + \frac{1}{\epsilon} \sigma \kappa \nabla \alpha_l \qquad (6) $$

With:

$$ \psi = \begin{cases} 1, & \text{in porous region} \\ 0, & \text{in non-porous region} \end{cases} \qquad (7) $$

$$ \rho_f = \alpha \rho_l + (1 - \alpha) \rho_v \qquad (8) $$

$$ \mu_f = \alpha \mu_l + (1 - \alpha) \mu_v \qquad (9) $$

the geometric function ($C_f$) and the permeability of porous media (K) are given by, respectively [62]:

$$ C_f = \frac{1.75}{\sqrt{150\epsilon^3}} \qquad (10) $$

$$ K = \frac{\epsilon^3 d_p^2}{150(1 - \epsilon)^2} \qquad (11) $$

Using the continuous surface force model (CSF) [63], the last term in Eq. (6) is surface tension.
The interface curvature $κ$ is:

$$κ = -\nabla \left( \frac{\nabla \tilde{\alpha}_L}{|\nabla \tilde{\alpha}_L|} \right) \qquad (12)$$

$\tilde{\alpha}_L$ is the smoothed VOF function determined by filtering $\alpha_L$ across a limited region at the interface. Over a thin region, the fluid interface in the VOF method changes sharply.

### Energy equation

Assuming that the flow is incompressible, ignoring the viscosity dissipation term, and taking into account that fluid and solid are in local thermal equilibrium, the energy equation is provided by:

$$(\rho C_p)_{\text{eff}} \frac{\partial T}{\partial t} + (\rho C_p)_{\text{eff}} V \nabla T = \nabla k_{\text{eff}} \nabla T + \dot{m}''_{L\text{hg}} \qquad (13)$$

Where the effective heat capacity is:

$$(\rho C_p)_{\text{eff}} = (1 - \varepsilon)(\rho C_p)_{\text{sm}} + \varepsilon[\alpha_L(\rho C_p)_1 + (1 - \alpha_L)(\rho C_p)_v] \qquad (14)$$

and the thermal conductivity is:

$$k_{\text{eff}} = (1 - \varepsilon)k_{\text{sm}} + \varepsilon[\alpha_L k_L + (1 - \alpha_L)k_v] \qquad (15)$$

## 2.2. Pressure velocity coupling

With OpenFOAM solvers we have different pressure-velocity coupled algorithms already available. Pressure implicit split operator PISO [64] for transient problems, semi-implicit method for pressure linked equation SIMPLE [65] for steady state cases or Pimple [66] algorithm which is a combination of these two algorithms are used. Such a combination allows the user to set larger time steps. These algorithms are inherently iterative techniques for linking momentum and mass equations and therefore may experience some computing-cost shortcoming. To ensure mass conservation and satisfy momentum conservation, they solve a pressure equation while explicitly correcting the velocity. The main difference between them is how they loop over the equation, even if they all answer the identical governing equation.

In our case, we used the Pimple algorithm, so we need to develop the pressure-Poisson equation for the domain.

By applying the divergence operator to the momentum equation Eq. (6), the pressure-Poisson expression can be obtained:

$$\nabla \left( \frac{\nabla p}{a_p} \right) = \nabla \left( \frac{1}{a_p} H(U) \right) - \nabla \left( \frac{\rho g x}{a_p} \right) - \nabla \left( \frac{\psi}{a_p} \left( \frac{\mu}{K} V + \frac{\rho C_f}{\sqrt{K}} |V| \right) \right) + \nabla \left( \frac{1}{\varepsilon a_p} \sigma \kappa \nabla \alpha_L \right) - \nabla \left( \frac{\psi}{a_p} \frac{2\sigma}{R} \nabla \alpha_L \right) - \dot{m}'' \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right) \qquad (16)$$

Where $a_p$ is the diagonal coefficient of the velocity matrix and $H(U)$ incorporates source terms and off-diagonal coefficients. These terms are evaluated from the previous time step or the previous iteration of the velocity field.

## 2.3. Phase transport equation for volume fraction

The volume fraction is expected to capture the interface because it is moving with the flow.

In order to develop the transport equation of the volume fraction, we explore the global continuity equation.

$$\frac{\partial \rho_f}{\partial t} + \nabla \rho_f V = 0 \qquad (17)$$

Where the fluid density ($\rho_f$) is given by:

$$\rho_f = \alpha_L \rho_l + (1 - \alpha_L) \rho_v \qquad (18)$$

By replacing $\rho_f$ by its expression (Eq. (18)) into Eq. (17) and rear-ranging, we obtain the following equation: (see Appendix B)

$$\frac{\partial \alpha_L}{\partial t} + V \nabla \alpha_L + \alpha_L \nabla V = - \frac{\rho_v \nabla V}{(\rho_L - \rho_v)} \qquad (19)$$

In our case with phase change process in porous media, the transport equation is determined from the continuity equation (Eq. (5)). We derive the following transport equation for $\alpha_L$:

$$\frac{\partial \alpha_L}{\partial t} + V \nabla \alpha_L = - \frac{\dot{m}''_{lv}}{\varepsilon} \left( \frac{1}{\rho_l} - \alpha_L \left( \frac{1}{\rho_l} - \frac{1}{\rho_v} \right) \right) \qquad (20)$$

The additional divergence term is included in the transport equation in OpenFOAM. This term only makes a contribution in the interface area ($0 < \alpha_L < 1$) [67]:

$$\frac{\partial \alpha_L}{\partial t} + V \nabla \alpha_L + \nabla (\alpha_L(1 - \alpha_L)U_c) = - \frac{\dot{m}''_{lv}}{\varepsilon} \left( \frac{1}{\rho_l} - \alpha_L \left( \frac{1}{\rho_l} - \frac{1}{\rho_v} \right) \right) \qquad (21)$$

Where $U_c$ is the artificial compressive velocity. It is determined, to avoid any excessive interface thickness dispersion, in the normal direction to that interface [68].

## 2.4. Phase change model

The main purpose is to study the mass and energy transfer for liquid-gas phase change under specific conditions. Thus, a new phase-change model can be rapidly incorporated. Here, we have proposed the change phase model 'Lee model' [69].

This model was developed under the assumption of a quasi-thermo-equilibrium condition and with the assumption that mass is transported in phase change flow at a constant pressure. The volumetric transferred mass rate is calculated as follows:

$$\dot{m}'' = r_c(1 - \alpha_L)\rho_v \frac{T - T_{\text{sat}}}{T_{\text{sat}}} \quad \text{for condensation } T < T_{\text{sat}} \qquad (22)$$

$$\dot{m}'' = r_c \alpha_L \rho_L \frac{T - T_{\text{sat}}}{T_{\text{sat}}} \quad \text{for evaporation } T > T_{\text{sat}} \qquad (23)$$

Where the empirical coefficients $r_c$ and $r_e$ (in s⁻¹) are known as the mass transfer intensity factors. According to the literature, their value should be set between 0.1 and 10⁷ s⁻¹ [70,71]. It is also suitable to choose a ratio $r_e/r_c$ almost equal to the fluid phase's density ratio $\rho_v/\rho_L$ to insure numerical stability of the solution [70].

## 3. Numerical details

We present in this section details for numerical implementation in the open-source software OpenFOAM. It is a flexible and powerful C++ toolbox that includes software (solvers and utilities) and operating functions for scalar, vector, and tensor fields [49]. It depends mainly on the Finite Volume Method and is highly customizable and extensible for new applications involving any set of partial differential equations.

The current solver porousPhaseChangeFoam is developed from an existing solver, namely interPhaseChangeFoam (OpenFOAM, version 6). The native solver can be used to tackle fluid flow problems involving cavitation. It is able to handle the flow of two incompressible and isothermal immiscible fluids with possible mass transfer. The VOF (Volume Of Fluid) phase-fraction-based method is used to capture the moving interface.

Although the set of phase-change models is offered to simulate cavitation, the same solver framework can also be used to implement additional phase-change models. In our case, we implemented the Lee model already described in Section 2.4.

The porousPhaseChangeFoam solver is implemented in OpenFOAM to simulate two-phase flow in porous media in non-isothermal conditions by solving Eqs. (6), (13), (16) and (21). In fact, this equation set is implicitly discretized using the finite volume method, and solved at each time step following the algorithm outlined in Appendix C.
The main features of the new solver are summarized as follows:

1. The transfer equation for thermal energy (file TEqn.H) is added to the solver to estimate the temperature distribution at each time step.

2. New library for the phase change model of the Lee model (Lee.C) is considered to solve evaporation and condensation phenomena and may be chosen by the user during operation.

3. Additionally, the governing equations and volume fraction transport equation (file alphaEqn.H) are adjusted with the proper source terms to take into consideration porous terms, by adapting the original velocity equation (file UEqn.H) and the pressure equation (file pEqn. H).

It is also worth mentioning that the liquid phase is considered the primary-phase. A transient solution with a self-adjusted time step is carried out in the model. A maximum global Courant number of 0.2 was set to ensure the stability of the simulation solutions. The pseudo-steady state regime is reached when the main parameters (temperature, pressure, and velocity) are almost stable (less than 3 % of time-step fluctuations). A (16-core in CPU, 12Go in RAM) computer is used for the development, the test, and the simulation of the investigated numerical cases.

## 4. Validation

In order to validate the relevance of the proposed solver, different thermal flow problems are proposed. For the phase change process, the Stefan problem and the Horizontal film condensation test cases are investigated. For the fluid flow through porous media, we suggest these different test cases:

- The first case: the incompressible isothermal fluid flow in a blank pipe driven by external constant forces

- The second case: the incompressible isothermal fluid flow in a porous conduit driven by external constant forces

- The third case: the incompressible isothermal fluid flow through a porous pipe at a constant velocity against the pipe wall

- The fourth case: the incompressible isothermal fluid flow in a conduit filled partially with porous medium.

### 4.1. Phase change validation

#### 4.1.1. Stefan problem

The 1D Stefan problem is commonly used [72,73] as a standard test case for phase change flows. The saturated liquid is separated from a superheated wall by a vapor film. Initially, liquid and vapor are in equilibrium, but then, evaporation forces are expected to remove the liquid from the hot wall.

In Fig. 1, a sketch of the Stefan problem is presented. The superheated wall's temperature is 10 K higher than the saturation

<figure><img src="image_3.png" /><figcaption>Fig. 1. Sketch of Stefan problem.</figcaption></figure>

temperature.

The analytical solutions for the position of interface $\theta_{an}(t)$ and temperature $T(x, t)$ on the vapor side are as following [72]:

$$ \theta_{an}(t) = 2\eta\sqrt{\lambda_v}t \qquad (24) $$

$$ T(x, t) = T_w + \frac{T_{sat} - T_w}{\text{erf}(\eta)} \text{erf}\left(\frac{x}{2\sqrt{\lambda_v}t}\right) \qquad (25) $$

$\eta$ is determined as following:

$$ \eta \text{erf}(\eta^2) \text{erf}(\eta) = \frac{C_{pG}(T_w - T_{sat})}{\sqrt{\pi H_{LG}}} \qquad (26) $$

The thermophysical properties of liquid and vapor are those of water at the saturation pressure ($P_{sat} = 1.0$ MPa). To maintain a constant coefficient of mass flow rate during the phase transition process in the CFD model, the specific heat values of the liquid and vapor phases are considered equal.

Fig. 2 shows the great agreement between the CFD results and the analytical solution of the described Stephan test case.

The maximum relative error is less than 0.83 %.

The temperature distribution along the domain, at various time values, is depicted in Fig. 3. As it can be seen, the temperature on the vapor side predicts the analytical solution in an almost perfectly quiescent fluid, and it is equivalent, as expected, to the saturation temperature on the liquid side. These test results confirm the validity of the applied method to simulate simple, one dimensional, heat transfer through a liquid-vapor domain.

#### 4.1.2. Horizontal film condensation

In this test case, at a temperature $T_w$ on the top surface of an isothermal plate, the vapor condenses into a liquid film, as shown in Fig. 4. The temperature profile is assumed to be linear from $T_w$ at the bottom wall to $T_{sat}$ at the interface; an analytical solution for the film thickness is obtained from a control volume study [74]:

$$ \delta_{an}(t) = \left[ 2t \frac{k_L}{\rho_L C_{pL}} \left( \frac{1}{2} + \frac{H_{LG}}{C_{pL} \Delta T} \right)^{-1/2} \right]^{1/2} \qquad (27) $$

The thermophysical properties chosen are for the saturated water at $P_{sa} = 1$ MPa.

In the directions tangential to the subcooled wall, a quasi-one-dimensional simulation domain with a single grid cell is taken into consideration. $T_{sat} = 453.03$ K and $T_w = 423.03$ K are taken into

<figure><img src="image_4.png" /><figcaption>Fig. 2. Comparison of the interface position in Stefan problem between the CFD model and the exact solution.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Temperature distribution at different time values comparison between CFD results (solid lines) and analytical results (symbols).</figcaption></figure>

<figure><img src="image_6.png" /></figure>

$$T_W = T_{sat} - \Delta T_{sub}$$

<caption>Fig. 4. Sketch of the horizontal film condensation.</caption>

consideration for a 30 K subcooled wall case. Every liquid that enters through the top side has the same characteristics as vapor. A very thin liquid film is taken into account at the beginning of the simulation in the computational domain to ensure the numerical stability of the solution procedure.

Fig. 5 shows the results of the liquid film thickness development. They imply that the CFD and the analytical solutions are in very good

<figure><img src="image_7.png" /><figcaption>Fig. 5. Film thickness results for the horizontal film condensation case.</figcaption></figure>

agreement.

The maximum relative error is less than 0.6 %.

## 4.2. Validation of fluid flow simulation through porous media

### 4.2.1. Incompressible and isothermal fluid flow in a blank pipe driven by external constant forces

In this test case, an external force $g=(g_x, 0)$ is constant and drives an incompressible isothermal fluid through a conduit. Fig. 6 depicts the proposed problem in illustration form.

The porosity $\epsilon$ is set to 1 and the permeability $K \to \infty$ (i.e., a very large number) so that the porous pipe degenerates into a blank pipe. The fluid flow problem is assumed to be axisymmetric, and the simulation domain is $0 < x < 1$ and $0 < r < R$.

The boundary conditions are:

$$\frac{\partial u_x}{\partial r} = u_r = \frac{\partial T}{\partial r} = 0 \text{ at } r = 0 \qquad (28)$$

$$u_x = u_r = 0, \ T = T_w = \text{constant, at } r = R \qquad (29)$$

where $u_x$ and $u_r$ are the velocity components.

The analytical answer to the Hagen Poiseuille flow's axial velocity is provided by:

$$u_x(r) = u_0 \left( 1 - \frac{r^2}{R^2} \right) \qquad (30)$$

Where

$$u_0 = g_x \frac{R^2}{4\nu} \qquad (31)$$

$u_0$ is the maximum velocity, and $g_x$ is an external force that is constant and exerted on the flow. Different test cases are simulated using the developed code and under the same working kinetic viscosity of $\nu=0.1 \text{ m}^2/\text{s}$.

The velocity profile is shown in Fig. 7 for various external force values $g_x$, and it can be seen that there is good consistency between simulation and analytical results, demonstrating the effectiveness of the numerical model. The external force term directly affects velocity, as can be seen from the analytical solution: the maximum velocity rises with the force term.

The maximum relative error is less than 1.2 %.

### 4.2.2. Incompressible isothermal fluid flow in a porous conduit driven by external constant forces

The next test case, as illustrated in Fig. 8, is an incompressible, isothermal fluid flow in a porous conduit that is driven by a constant external force $g_x$.

The same previous test case's initial and boundary conditions are applied here. The kinetic viscosity and the porosity are set to $0.01 \text{ m}^2/\text{s}$ and $0.1$, respectively.

The permeability $K$ ranges from $10^{-4}$ to $10^{-2} \text{ m}^2$, and the constant force $g_x$ is equal to $10 \text{ m s}^{-2}$.

By removing the Forchheimer term, which captures the nonlinear impact on flows across the porous medium, the problem can be analytically resolved [75]:

<figure><img src="image_8.png" /><figcaption>Fig. 6. Sketch of the first case.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Comparison of the axial velocity profiles resulted from the CFD simulation (solid lines) for the first test case for different external force values with the analytical solution (symbols).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Sketch of the second test case.</figcaption></figure>

$$u_x(r) = \text{Real} \left\{ \frac{g_x K}{\nu} - \frac{g_x K J_0((-\epsilon/K)^{1/2}r)}{J_0((-\epsilon/K)^{1/2}R)\nu} \right\} \quad (32)$$

Where 'Real' refers to the real portion of a complex variable, while $J_0$ is the first type's zeroth-order Bessel function.

In Fig. 9, the analytical solution and numerical results of the CFD model are compared in terms of the axial velocity profile for various permeability values. As noticeable, a very good agreement between the analytical and numerical data is obtained. In addition, the permeability has an important impact on the velocity profile and magnitude: increasing the permeability value causes an increase in the maximum

<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison of the axial velocity obtained from the CFD simulation (solid lines) for the second test case for different permeability with the analytical solution (symbols) without Forchheimer term.</figcaption></figure>

velocity, which is explained by a reduction in the Darcy force effect itself as a consequence of the increased fluid's passage through the porous medium.

The maximum relative error is less than 0.9 %.

The flow through porous media can also be significantly affected by inertial forces. Therefore, a deeper examination of the Forchheimer term's impact is made. We evaluated our numerical results against those obtained using finite difference solutions [76] for various permeability values. By comparing the findings from the present model with those provided by the central difference scheme [76,77], we obtained a good agreement (Fig. 10).

The maximum relative error is less than 2 %.

We can notice from Figs. 9 and 10 that using the Forchheimer term affects significantly the evolution of the axial velocity: when increasing the permeability, the axial velocity decreases. In addition, the velocity profile becomes more flattened when the Forchheimer term is considered. For this case such a term should not be neglected.

### 4.2.3. Incompressible isothermal fluid flow through a porous pipe at a constant velocity against the pipe wall

This test case consists in an isothermal fluid flow that is incompressible and propelled by a constant velocity on a porous pipe wall. Fig. 11 shows the problem's configuration, with the flow, in this example, being pushed by a constant velocity against the wall, $u_0 = u(x, R) = 1$ m/s. The remaining boundary conditions are unchanged compared to the first and second cases.

First, the issue is examined without the Forchheimer term, and the analytical solution is provided as follows [75]:

$$u_x(r) = \text{Real} \left\{ \frac{u_0 J_0((-\epsilon/K)^{1/2}r)}{J_0((-\epsilon/K)^{1/2}R)\nu} \right\} \quad (33)$$

A striking agreement is obtained in this test case as well for various K values when the numerical findings are contrasted to the analytical solution (Fig. 12). Along the radial axis, the axial velocity rises until it reaches its maximum value at the pipe wall. Additionally, like in the previous case, we can see that lowering the permeability results in a more flattened velocity profile.

The maximum relative error is less than 2 %.

The impact of the Forchheimer factor in the momentum equation, is then appreciable. There is indeed a good agreement found when comparing the computational results from the current model with the finite difference approach data [75,77] (Fig. 13).

<figure><img src="image_12.png" /><figcaption>Fig. 10. The axial velocity profiles resulted from the CFD simulation (solid lines) for the second test case for different permeability compared with analytical solution (symbols) with Forchheimer term.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11. Sketch of the third case.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Comparison of the axial velocity obtained from the CFD simulation (solid lines) for the third test case for different permeability values with analytical solution (symbols) without Forchheimer term.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Comparison of the axial velocity resulted from the CFD simulation (solid lines) for the third test case for different permeability values with finite difference solution (symbols) with Forchheimer term [75].</figcaption></figure>

The maximum relative error is less than 3 %.

When the Forchheimer term is taken into account, the velocity profile is considerably flatter in Figs. 12 and 13.

### 4.2.4. Incompressible isothermal fluid flow in a conduit partially filled with porous medium

The next test scenario, as shown in Fig. 14, involves the flow of an incompressible heated fluid in a pipe that is partially filled with a porous medium. It was fully described by Rong et al. [78]. R<sub>p</sub> is the radius of the porous media, and R<sub>0</sub> is the entire pipe radius. By assuming the local

<figure><img src="image_16.png" /><figcaption>Fig. 14. Sketch of the fourth case.</figcaption></figure>

equilibrium, the porous media is homogeneous and saturated with a single-phase fluid.

The numerical solutions are compared to those reported by Rong et al. [78] for the following set values: ε = 0.4, Da=10<sup>-3</sup>, R<sub>p</sub> = 0.5 cm, Re = 100, and Pr = 0.7.

A good agreement is found between the CFD results and those suggested by Rong et al [78], as shown in Figs. 15 and 16.

Fig. 15 shows that the axial velocity variations are lower in the porous zone than in the clear zone. The existence of the porous structure generates forces that make the fluid circulation along the pipe radius more uniform. However, the equivalent viscosity is lower. Thus, the axial velocity is higher in the clear zone (the region free of a porous structure). The maximum relative error is less than 1.5 %.

Fig. 16 shows the temperature distribution along the domain, which appears to be constant in the clear zone and rises in the porous region to reach the wall temperature value. The smaller velocity gradient in the porous zone and the higher effective thermal conductivity in the porous material can both be used to explain the strong link between the velocity and temperature profiles. The maximum relative error is less than 0.002 %.

## 5. The application of the new solver to liquid coolant phase change during transpiration cooling through a porous medium

In the previous sections, we have shown the effectiveness of the developed numerical tool in the porous medium and independently validated the phase change model. It is then interesting to test the entire model implemented in the solver to assess its efficiency to simulate a more complicated phase change process in the porous zone.

Transpiration cooling with phase change is a highly effective method of heat dissipation, in which a coolant penetrates the porous structure, absorbs heat, and produces a protective film as it flows out, maintaining a moderate temperature across the entire structure.

The transpiration-cooling model shown in Fig. 17 is taken as an example. In this process, a porous matrix is injected with a liquid coolant (temperature T<sub>c</sub>, mass flow rate m<sub>c</sub>). The heat applied to the upper wall

<figure><img src="image_17.png" /><figcaption>Fig. 15. The axial velocity of the fourth test case.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 16. Comparison of the temperature profile of the fourth test case.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Physical model for liquid coolant phase change during transpiration cooling in the porous medium.</figcaption></figure>

causes a phase change.

The boundary conditions are:

At the upper side:

$$P_v|_{y=0} = P_0$$

$$k_{\text{eff}} \left. \frac{\partial T}{\partial t} \right|_{y=0} = Q$$

left and right vertical sides:

$$k_{\text{eff}} \left. \frac{\partial T}{\partial t} \right|_{y=0} = 0 \quad (\text{Adiabatic})$$

Coolant inlet: at y=-L:

$$\rho_1 \varepsilon V_1 = \dot{m}(x) \quad (\text{Coolant mass flow rate})$$

## 5.1. Mesh independence test

Numerical tests using three different grid numbers were carried out. Various mesh stats points are considered: Mesh 1: 204, Mesh 2: 304, Mesh 3: 404. The obtained numerical results are used to examine the mesh independence.

The main features of the simulation are set: the porosity (0.35), the heat flux (1.5 MW/m²) and the mass flux of coolant (0.4 kg/m²s). Fig. 18 shows the temperature distribution along the porous matrix. It is obvious that the simulation results provided by the mesh configurations

<figure><img src="image_20.png" /><figcaption>Fig. 18. Temperature distribution under three different meshes.</figcaption></figure>

2 and 3 are similar. Hence, in the following numerical simulations, the mesh configuration 2 is preferred in order to save computational resources and ensure numerical accuracy.

## 5.2. Verification of numerical model through experimentation

The experimental set up [79] is illustrated in Fig. 19. The experimental data in Ref. [79] are utilized to validate the proposed numerical model. Made of sintered metal copper with a porosity of 0.318 and a particle diameter of 0.20 mm, the porous matrix was a cylinder with a radius of 40 mm and a length of 100 mm. Above the cylinder surface, a xenon lamp is mounted to produce a 0.21 MW/m² heat flux. At a constant temperature of 300 K, two water mass flow rate values were used: 5.67 g/min and 7.45 g/min. To measure the temperature of the porous matrix, seven thermocouples were arranged along the direction of the coolant flow.

As shown in Fig. 20, the simulation results and experimental data are quite similar and show a good agreement. The maximum relative error is less than 6.5 %. It is worth noticing that the accordance of the results is better when the coolant mass flow rate is increased. In such a case, the temperature within the porous matrix is sensibly lowered, mainly in the top region (z > 0.05 m). The numerical approach and its mathematical model are considered valid.

## 5.3. Results and discussions

In this section, a porous section with a thickness of 0.1 m, a pore diameter of 5.10⁻⁴ m, a porosity of 0.35, and a thermal conductivity of 30 W/mK is considered. The fluid flows opposite to gravity (Fig. 17). Indeed, liquid water was injected from the lower surface into the porous matrix with different mass flow rates (0.3, 0.4, and 0.5 kg/m²s). A heat flux of 1.5 MW/m² is applied to the upper surface. The thermophysical properties of the coolant fluid are summarized in Table 1.

The longitudinal variation of the temperature in the porous media is illustrated in Fig. 21 at various mass flow rates. Figs. 21 and 22 show the influence of this parameter (coolant mass flow rate) on the main porous matrix characteristics (the temperature and the saturation distribution, respectively). It has been demonstrated that increasing the coolant injection mass delays the occurrence of phase changes within the porous medium. It reduces therefore the area of the vapor phase and the two-phase sections. The surface temperature exposed to the heat flow thus decreases. Also, we can notice that from the flow rate value (ṁ = 0.5 kg/m²s), the temperature remains unchanged despite the variation in coolant injection rate.
<figure><img src="image_21.png" /><figcaption>Fig. 19. schematic model of experimental setup [79].</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Comparison of numerical results and experimental data.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 21. Temperature distribution with various coolant mass flow rate.</figcaption></figure>

<table><caption>**Table 1**

Thermophysical properties of saturated water at the reference temperature (373.15 K).</caption><thead><tr><th>Properties</th><th>Liquid</th><th>Vapor</th><th>Units</th></tr></thead><tbody><tr><td>Density</td><td>1000</td><td>0.6</td><td>kg/m<sup>3</sup></td></tr><tr><td>Specific heat</td><td>4217</td><td>2029</td><td>J/kg K</td></tr><tr><td>Thermal conductivity</td><td>0.68</td><td>0.02</td><td>W/m K</td></tr><tr><td>Dynamic viscosity</td><td>279 × 10<sup>-6</sup></td><td>12.02 × 10<sup>-6</sup></td><td>kg/m s</td></tr><tr><td>Latent heat of evaporation</td><td>2,675,710</td><td></td><td>J/kg</td></tr></tbody></table>

Fig. 22 represents the saturation distribution using the sharp interface method, which justifies the very thin liquid-vapor interface. Let us just remind that the VOF method used here describes only sharp interfaces of numerically zero thickness. Therefore, the different physical properties, like viscosity and density, are necessarily discontinuous through that interface. It can be somewhat considered as a native limitation of the considered numerical approach.

Fig. 23 shows the CFD results of the saturation distribution for the case of $\dot{m} = 0.3$ kg/m²s as a function of time: the different phases are clearly represented. This figure demonstrates the condensation phenomena: the phase change process occurs. It is emphasized by the evolution of the liquid phase and consequently the interface position over the time progress ($t = 0$ s and $t = 50$ s).

Meanwhile, the transpiration cooling efficiency can be defined as [80]:

<figure><img src="image_24.png" /><figcaption>Fig. 22. Saturation distribution with various coolant mass flow rate.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 23. Saturation distribution at different time steps.</figcaption></figure>

<table><caption>**Table 2**
The average cooling efficiency variations according to the coolant mass flow.</caption><tr><td>ṁ (kg/m²s)</td><td>0.3</td><td>0.4</td><td>0.5</td><td>0.7</td></tr><tr><td>T<sub>w</sub> (K)</td><td>1600</td><td>933</td><td>567</td><td>569.56</td></tr><tr><td>η</td><td>0.17</td><td>0.62</td><td>0.86</td><td>0.86</td></tr></table>

$$ \eta = \frac{T_{\text{adia}} - T_{\text{w}}}{T_{\text{adia}} - T_{\text{c}}} \qquad (34) $$

Where T_c is the coolant inlet temperature. T_adia and T_w are respectively the top surface temperatures in the adiabatic case (without cooling) and the transpiration cooling case.

For the case (Q = 1.5 MW/m² and e = 0.35), the obtained transpiration cooling efficiency values are listed in Table 2.

It illustrates the cooling performance evolution when the coolant mass flow rate is increased. It decreases with the increase of the top wall temperature, resulting from heat accumulation and relatively bad heat transfer. It can be seen that the transpiration cooling is more efficient when the coolant flow rate increases until reaching the relative maximum (η = 0.86) with a flow rate equal to 0.5 kg/m²s; over that no variation is noticed. The same observation was noticed in previous experimental investigations on transpiration cooling [7]. The obtained values of the cooling efficiency η are comparable with the results of previous studies [7,20]. Values over 0.9 had indeed been reached when the coolant flow rate was increased [7].

For an optimization purpose of transpiration cooling, a mass flow rate of 0.5 kg/m²s should be set for the investigated case (Q = 1.5 MW/m² and e = 0.35). However, the improvement of the cooling performance that can be obtained by simply increasing the mass flow of the input coolant is limited. Other linked variables, like the cooling fluid inlet temperature, the porosity and the dimensions of the porous matrix, are expected to have deep effects on the mean efficiency of such a particular cooling process. Thus, a more in-depth study is necessary to perform an optimization investigation of the phase change transpiration cooling in porous media.

## 6. Conclusions

Developing numerical simulation has proven to be an interesting means to understand detailed mechanisms that occur within porous media, especially when two-phase fluid undergoes phase change flows. In this study, we successfully developed and validated a numerical solver (porousPhaseChangeFoam) for two-phase flow in porous media based on the OpenFOAM toolbox. This solver was successfully verified by comparing its results with some analytical solutions and/or the data from previous CFD-based studies. It involved: the two-phase Stefan problem, the Horizontal film condensation case, the incompressible and isothermal single-phase fluid flow in a blank pipe, and then in a porous conduit. We have also successfully tested the cases of the incompressible isothermal single-phase fluid flow through a porous pipe at a constant velocity against the pipe wall, and the incompressible heated single-phase fluid flow in a pipe partially filled with a porous medium. In particular, the impact of the Forchheimer factor, in the momentum equation, is shown to be appreciable. There is indeed a good agreement found when comparing the computational results from the current model with the finite difference approach data of other numerical test investigations [75,77]. Furthermore, the velocity profile is considerably flattened when the Forchheimer term is taken into account.

Finally, the phase change transpiration cooling in porous media was investigated. Thus, the obtained results demonstrated the efficiency of the new solver in studying accurately the phase change heat transfer problem in porous medium. The effect of mass flow rate on cooling efficiency was then examined. The results show that increasing the coolant flow rate can improve the efficiency of the process until a limit value at which point the variation effect becomes slight. However, more in-depth optimization study is needed since other linked parameters (porosity, matrix thickness, and coolant temperature) have a deep impact on the cooling efficiency of such a process.

The developed numerical tool can be easily implemented in simulation models of many other complex engineering devices such as heat pipes and their variants. In fact, the key element of such heat transfer devices is usually the porous medium where the phenomena of two-phase flow, phase change, and moving interface must be taken into account.

## CRediT authorship contribution statement

Aroua Ghedira: Writing – original draft, Validation, Supervision, Methodology, Investigation, Conceptualization. Zied Lataoui: Writing – original draft, Investigation. Adel M. Benselama: Visualization, Methodology. Yves Bertin: Visualization. Abdelmajid Jemni: Writing – original draft, Supervision, Investigation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This research was carried out as part of the PHC-Utique program supported both by the French and Tunisian governments under contract CMCU no. 23G1104.
## Appendix A. The Momentum equation derivation

The differential balance equation is given by [59]:

$$ \frac{\partial \rho_f \Phi_f}{\partial t} + \nabla (V_f \rho_f \Phi_f) = -\nabla J_f + \rho_f X_f \quad (A.1) $$

The first term of the above equation is the time rate of change of the quantity per unit volume, whereas the second term is the rate of convection per unit volume. The right-hand side terms represent the surface flux and the volume source.

The conservation of momentum can be obtained from Eq. (A.1) by introducing the surface stress tensor $T_k$ and the body force $g_k$, thus we set:

$$ \Phi_f = V_f $$

$$ J_f = -T_k = p_k I - \tilde{C}_k $$

$$ X_f = g_f $$

where $I$ is the unit tensor. Here we have split the stress tensor into the pressure term and the viscous stress $\tilde{C}_k$. In view of Eq. (A.1) we have:

$$ \frac{\partial \rho_f V_f}{\partial t} + \nabla (\rho_f V_f V_f) = -\nabla p_f + \tilde{C}_k + \rho_f g_f \quad (A.2) $$

In general, two-phase flow systems with transport of mass, momentum and energy are characterized by the existence of two different densities and velocities.

The VOF model used in this study can be applied to model two immiscible fluids with a clearly defined interface between phases and is used for surface tracking applied to a fixed mesh. Also, in VOF technique, one set of Navier Stokes equations are solved through the computational domain and used to track the motion of different phases by defining the volume fraction of each phase.

The volume fraction of the ith phase or also named saturation, $\alpha_i$ in the VOF model is the ratio of the element volume of ith phase ($\Delta v_i$) to the total elemental volume of all phases ($\Delta v$):

$$ \alpha_i = \frac{\Delta v_i}{\Delta v} \quad (A.3) $$

The volume fraction of all the phases (n is the phases number) must be equal to unity:

$$ \sum_{i=1}^{n} \alpha_i = 1 \quad (A.4) $$

$$ \alpha_i = \begin{cases} 0, & \text{in a gas region} \\ 0; 1[, & \text{unsaturated, i.e. partly liquid - filled region} \\ 1, & \text{in the liquid region} \end{cases} $$

So, we can replace terms in Eq. (A.2) by the following [59,61]:

$$ \rho_f = \alpha_L \rho_l + (1 - \alpha_L) \rho_v $$

$$ V_f = \alpha_L V_l + (1 - \alpha_L) V_v $$

$$ p_f = \alpha_L p_l + (1 - \alpha_L) p_v $$

$$ \frac{\partial}{\partial t}(\alpha_i \rho_l V_l) + \nabla (\alpha_i \rho_l V_l V_l) = -\nabla \alpha_i p_l + \alpha_i \tilde{C}_i + \alpha_i \rho_l g_l $$

i is the phase (liquid or vapor)

The forces acting in the fluid are the gravity, pressure and surface tension [61].

Considering that the fluid is incompressible, we can obtain the following equation:

$$ \rho_f \left[ \frac{\partial V_f}{\partial t} + V_f \nabla V_f \right] = -\nabla p + \rho_f g + \sigma \kappa \nabla \alpha_i \quad (A.5) $$

The porous region is assumed to be assumed homogenous, isotropic and saturated with a single-phase fluid that is in equilibrium with the solid.

The following details presents a development of the momentum equation [60].

The term “porosity” $\epsilon$ is most usually used to describe the solid matrix’s void fraction:

$$ \epsilon = \begin{cases} 0; 1[, & \text{in the porous region} \\ 1, & \text{in the solid - free region} \end{cases} $$

Henry Darcy developed what is now known as Darcy's law, which describes the link between the pressure gradients/resistance and the dynamic viscosity, $\mu_f$, the extrinsic phase-averaged velocity, and the permeability, $K$:

$$ \rho_f g - \frac{1}{\epsilon} \nabla (\epsilon P) = \frac{\mu_f}{K} \langle V_f \rangle \quad (A.6) $$

To take into account for non-linearity in porous zone, Ward proposed a quadratic drag term dependent of $\rho_f K^{1/2}$
$$ \rho_f g - \frac{1}{\varepsilon} \nabla (eP) = \frac{\mu_f}{K} \langle V_f \rangle + \frac{\rho_f C_f}{K^{1/2}} |\langle V_f \rangle| \langle V_f \rangle \quad (A.7) $$

where C_f is a dimensionless drag constant.

As in the case of packed beds composed of spheres, the effects of the nearby spheres are included in the Stokes drag force on a sphere in an infinite domain. The drag exerted by a single sphere in an infinite domain applies to a flow with negligible inertia is known as the Stokes drag force. This was achieved by overlaying the Darcy flow with the Stokes flow, known as the Brinkman equation:

$$ \rho_f g - \frac{1}{\varepsilon} \nabla (eP) = \frac{\mu_f}{K} \langle V_f \rangle - \frac{\mu_f}{\varepsilon} \nabla^2 \langle V_f \rangle \quad (A.8) $$

Combining the inertial drag component, Eq. (A.7), Brinkman's equation, Eq. (A.8), and Darcy's law, Eq. (A.6), can heuristically yield an empirical momentum equation:

$$ \rho_f g - \frac{1}{\varepsilon} \nabla (eP) = \frac{\mu_f}{K} \langle V_f \rangle + \frac{\rho_f C_f}{K^{1/2}} |\langle V_f \rangle| \langle V_f \rangle - \frac{\mu_f}{\varepsilon} \nabla^2 \langle V_f \rangle \quad (A.9) $$

Using the continuity equation, and assuming constant porosity, the momentum equation becomes

$$ \rho_f \left[ \frac{1}{\varepsilon} \frac{\partial \langle V_f \rangle}{\partial t} + \frac{\langle V_f \rangle}{\varepsilon^2} \nabla \langle V_f \rangle \right] = - \nabla P + \frac{\mu_f}{\varepsilon} \nabla^2 \langle V_f \rangle + \rho_f g - \frac{\mu_f}{K} \langle V_f \rangle - \frac{\rho_f C_f}{K^{1/2}} |\langle V_f \rangle| \langle V_f \rangle \quad (A.10) $$

The governing equations must be specified for the liquid and vapor as well as for the solid matrix when the porous medium is saturated with both liquid and vapor. Solid-liquid-vapor interactions play an important role in heat and mass transfer applications. Therefore, we must account for a possible phase change between the liquid and vapor in the governing equation for the liquid and vapor phases. The liquid and vapor are considered to be in local thermal equilibrium. We should also add the pressure gradient due to the change in capillary pressure, p_c, as a function of volume fraction of liquid [60]:

$$ pc = \frac{2\sigma}{R_p} \nabla \alpha_l $$

in this case, considering porous media, the velocity and the pressure become:

$$ V_f = \varepsilon \alpha_L V_L + \varepsilon (1 - \alpha_L) V_V $$

$$ p_f = \varepsilon \alpha_L p_L + \varepsilon (1 - \alpha_L) p_V $$

Finally, we can obtain the momentum equation in porous media with phase change. As mentioned above, one set of Navier Stokes equations are solved.

Our objective is to build a solver that can be adapted to any problem, i.e., a solver that can study different flows where there is a or no phase change phenomena in or not a porous medium, such as the example of mesh heat pipes. So, we add the term to the different porous terms in the equation

$$ \psi = \begin{cases} 1, & \text{in porous region} \\ 0, & \text{in non-porous region} \end{cases} $$

Finally, we can obtain the following equation:

$$ \rho_f \left[ \frac{1}{\varepsilon} \frac{\partial V}{\partial t} + \frac{V}{\varepsilon^2} \nabla V \right] - \frac{\mu_f}{\varepsilon} \nabla^2 V = - \nabla p + \rho_f g - \psi \left[ \frac{\mu_f}{K} V + \frac{\rho_f C_f}{\sqrt{K}} |V| V \right] - \psi \frac{2\sigma}{R_p} \nabla \alpha_l + \frac{1}{\varepsilon} \sigma \kappa \nabla \alpha_l \quad (A.11) $$

## Appendix B. Details for Eq. (37)

The continuity equation in the fluid phase is given by Eq. (35):

$$ \frac{\partial \rho_f}{\partial t} + \nabla (\rho_f V) = 0 \quad (35) $$

Considering the volume fractions of liquid and vapor contained in the medium, the density of the fluid is written as follows:

$$ \rho_f = \alpha_L \rho_L + (1 - \alpha_L) \rho_V \quad (36) $$

The development of this equation leads to:

$$ \frac{\partial (\alpha_L \rho_L + (1 - \alpha_L) \rho_V)}{\partial t} + \nabla (\alpha_L \rho_L + (1 - \alpha_L) \rho_V) V = 0 $$

$$ \frac{\partial (\alpha_L \rho_L)}{\partial t} + \frac{\partial \rho_V}{\partial t} - \frac{\partial (\alpha_L \rho_V)}{\partial t} + (\alpha_L \rho_L + (1 - \alpha_L) \rho_V) \nabla \vec{V} + \vec{V} \nabla (\alpha_L \rho_L + (1 - \alpha_L) \rho_V) = 0 $$

For incompressible flow densities of the liquid (ρ_L) and of the vapor (ρ_V) are supposed to be constant, then:

$$ \Rightarrow (\rho_L - \rho_V) \frac{\partial \alpha_L}{\partial t} + \alpha_L (\rho_L - \rho_V) \nabla \vec{V} + \rho_V \nabla \vec{V} + \vec{V} \rho_L \nabla \alpha_L + \vec{V} \nabla \rho_V - \vec{V} \rho_V \nabla \alpha_L = 0 $$

ρ_V is constant then (∇ρ_V = 0):
$$ \Rightarrow (\rho_L - \rho_V) \frac{\partial \alpha_L}{\partial t} + (\rho_L - \rho_V) \alpha_L \nabla \vec{V} + \rho_V \nabla \vec{V} + (\rho_L - \rho_V) \vec{V} \nabla \alpha_L = 0 $$

Then

$$ \frac{\partial \alpha_L}{\partial t} + \alpha_L \nabla \vec{V} + \vec{V} \nabla \alpha_L = \frac{-\rho_V \nabla \vec{V}}{(\rho_L - \rho_V)} \qquad (37) $$

## Appendix C. The diagram of the implemented solver porousPhaseChangeFoam

<figure><img src="image_26.png" /></figure>

## Data availability

The authors do not have permission to share data.

## References

[1] S. Bekri, O. Vizika, J.F. Thovert, et al., Binary two-phase flow with phase change in porous media, s.l. Int. J. Multiph. Flow 27 (3) (2001) 477-526.

[2] S. Weise, S. Meinicke, T. Wetzel, et al., Modelling the pressure drop of two-phase flow through solid porous media, Int. J. Multiph. Flow 112 (2019) 13-26.

[3] et T. Brahim, A. Jemni, Effect of the heat pipe adiabatic region, J. Heat Transf. 136 (4) (2014) 042901.

[4] P. Cheng, Heat transfer in geothermal systems. Advances in Heat Transfer, Elsevier, 1979, pp. 1-105.

[5] G. Gökçe, C. Kurt, G. Odabasi, et al., Comprehensive three-dimensional hydrodynamic and thermal modeling of steady-state operation of a flat grooved heat pipe, Int. J. Multiph. Flow 160 (2023) 104370.

[6] L.L. Vasiliev, Micro and miniature heat pipes-electronic component coolers, Appl. Therm. Eng. 28 (4) (2008) 266-273.

[7] G. Huang, Y. Zhu, Z. Liao, et al., Experimental investigation of transpiration cooling with phase change for sintered porous plates, Int. J. Heat Mass Transf. 114 (2017) 1201-1213.

[8] K.D. Song, S.H. Choi, S.J. Scotti, Transpiration cooling experiment for scramjet engine combustion chamber by high heat fluxes, J. Propuls. Power 22 (1) (2006) 96-102.

[9] Y.H. Andoh, B. Lips, Prediction of porous walls thermal protection by effusion or transpiration cooling. An analytical approach, Appl. Therm. Eng. 23 (15) (2003) 1947-1958.

[10] C. Trevin, A. Medina, Analysis of the transpiration cooling of a thin porous plate in a hot laminar convective flow, Eur. J. Mech. B Fluids 18 (2) (1999) 245-260.

[11] D. Keener, J. Lenertz, R. Bowersox, J. Bowman, Transpiration cooling effects on nozzle heat transfer and performance, J. Spacecr. Rockets 32 (6) (1995) 981-985.

[12] Transpiration cooling applied to C/C liners of cryogenic liquid rocket engines. Greuel, D., Herbertz, A., Haidn, O., Ortelt, M., & Hald, H. p. 2004.

[13] J. Landis, W. Bowman, Numerical study of a transpiration cooled rocket nozzle, in: Proceedings of the 32nd Joint Propulsion Conference and Exhibit, 1996.

[14] et S. Becker, A. Sokolichin, G. Eigenberger, Gas-liquid flow in bubble columns and loop reactors: Part II. Comparison of detailed experiments and flow simulations, Chem. Eng. Sci. 49 (24) (1994) 5747-5762.

[15] et E. Delnoij, J.A.M. Kuipers, W.P.M. Van Swaaij, Dynamic simulation of gas-liquid two-phase flow: effect of column aspect ratio on the flow structure, Chem. Eng. Sci. 52 (21-22) (1997) 3759-3772.

[16] C.F. Tai, J.N. Chung, A direct numerical simulation of axisymmetric cryogenic two-phase flows in a pipe with phase change, Comput. Fluids 48 (1) (2011) 163-182.

[17] et W. Wibel, U. Schygulla, J.J. Brandner, Micro device for liquid cooling by evaporation of R134a, Chem. Eng. J. 167 (2-3) (2011) 705-712.

[18] et D.B. Das, S.M. Hassanizadeh, Upscaling Multiphase Flow in Porous Media, Springer, Berlin, Germany, 2005. et.

[19] et L. Mottet, T. Coquard, M. Prat, Three dimensional liquid and vapour distribution in the wick of capillary evaporators, Int. J. Heat Mass Transf. 83 (2015) 636-651.

[20] H. Ma, L. Yang, Y. Luan, et al., Effect of pore plugging on transpiration cooling performance implemented with perforated flat plates, Int. Commun. Heat Mass Transf. 159 (Part C) (2024) 108248.

[21] et L. Cueto-Felgueroso, X. Fu, R. Juanes, Pore-scale modeling of phase change in porous media, Phys. Rev. Fluids 3 (8) (2018) 084302.

[22] et Y.C. Yortsos, A.K. Stubos, Phase change in porous media, Curr. Opin. Colloid Interface Sci. 6 (3) (2001) 208-216.

[23] S. Whitaker, Simultaneous heat, mass, and momentum transfer in porous media: a theory of drying. Advances in Heat Transfer, Elsevier, 1977, pp. 119-203.

[24] et F. Tariku, K. Kumaran, P. Fazio, Transient model for coupled heat, air and moisture transfer through multilayered porous media, Int. J. Heat Mass Transf. 53 (15-16) (2010) 3035-3044.

[25] F. He, T. Liu, R. Ding, J. Wang, Numerical investigation on the double-layer porous plate of transpiration cooling with phase change under different heat flux, Results Eng. 20 (2023) 101549.

[26] X. Wang, X. Fan, W. Li, et al., Exploration of local transpiration cooling method and numerical study of cooling performance, Int. Commun. Heat Mass Transf. 159 (Part C) (2024) 108229.

[27] F. Wang, X. Xia, W. Song, X. Chen, Experimental investigation on the cooling performance of a solid hybrid gel with self-transpiration cooling at high temperature, Appl. Therm. Eng. (2024) 124913.

[28] Wilson, A. James, Majid Haghshenas, and Ranganathan Kumar. Phase-change mechanism for evaporation in porous media using volume of fluid: Implicit formulation of interfacial temperature, Int. Commun. Heat Mass Transf. 103 (2019) 90-99.

[29] H.Y. Li, K.C. Leong, L.W. Jin, et al., Three-dimensional numerical simulation of fluid flow with phase change heat transfer in an asymmetrically heated porous channel, Int. J. Therm. Sci. 49 (12) (2010) 2363-2375.

[30] H.Y. Li, K.C. Leong, L.W. Jin, et al., Transient two-phase flow and heat transfer with localized heating in porous media, Int. J. Therm. Sci. 49 (7) (2010) 1115-1127.

[31] Description and utilization of interFoam multiphase solver. Damian, S.M. 2012, International Center for Computational Methods in Engineering, pp. 1-64.

[32] et C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201-225.

[33] C. Ierardi, A. Della Torre, G. Montenegro, et al., Modeling of two-phase flows at low Capillary number with VoF method, Comput. Fluids 252 (2023) 105772.

[34] DA. Hoang, V. Van Steijn, L.M. Portela, et al., Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method, Comput. Fluids 86 (2013) 28-36.
[35] G. Tronci, A. Buffo, M. Vanni, et al., Validation of the diffusion mixture model for the simulation of bubbly flows and implementation in OpenFOAM, Comput. Fluids 227 (2021) 105026.

[36] et N. Samkhaniani, M. Ansari, A VOF method to phase change modeling, in: Proceedings of the ICHMT, 2012.

[37] M. Wörner, Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications, Microfluid. Nanofluid. 12 (6) (2012) 841–886.

[38] et P. Angeli, A. Gavriliidis, Hydrodynamics of Taylor flow in small channels: a review, Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 222 (5) (2008) 737–751.

[39] et S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1) (1992) 25–37.

[40] et M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159.

[41] et J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354.

[42] P. Horgue, C. Soulaine, J. Franc, et al., An open-source toolbox for multiphase flow in porous media, Comput. Phys. Commun. 187 (2015) 217–226.

[43] et J.W. Sheldon, W.T. Cardwell Jr., One-dimensional, incompressible, noncapillary, two-phase fluid flow in a porous medium, Trans. AIME 216 (1) (1959) 290–296.

[44] B. Flemisch, M. Darcis, K. Erbertseder, et al., DuMex: DUNE for multi-(phase, component, scale, physics,...) flow and transport in porous media, Adv. Water Resour. 34 (9) (2011) 1102–1112.

[45] PFLOTRAN user manual. Lichtner, P.C., Hammond, G.E., Lu, C., et al. In: Technical Report : s.n., 2013.

[46] O. Kolditz, S. Bauer, L. Bilke, et al., OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media, Environ. Earth Sci. 67 (2012) 589–599.

[47] K.A. Lie, S. Krogstad, I.S. Ligaarden, et al., Open-source MATLAB implementation of consistent discretisations on complex grids, Comput. Geosci. 16 (2012) 297–322.

[48] J. Hrvoje. Error analysis and estimation in the Finite Volume method with applications to fluid flows, 1996.

[49] H.G. Weller, G. Tabor, H. Jasak, et al., A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys. 12 (6) (1998) 620–631.

[50] et R. Aziz, V. Joekar-Niasar, P. Martinez-Ferrer, Pore-scale insights into transport and mixing in steady-state two-phase flow in porous media, Int. J. Multiph. Flow 109 (2018) 51–62.

[51] Y. Yang, S. Cai, J. Yao, et al., Pore-scale simulation of remaining oil distribution in 3D porous media affected by wettability and capillarity based on volume of fluid method, Int. J. Multiph. Flow 143 (2021) 1037.

[52] et F. Piscaglia, A. Montorfano, A. Onorati, Multi-dimensional computation of compressible reacting flows through porous media to apply to internal combustion engine simulation, Math. Comput. Model. 52 (7-8) (2010).

[53] V. Novaresio, M. García-Camprubí, S. Izquierdo, et al., An open-source library for the numerical modeling of mass-transfer in solid oxide fuel cells, Comput. Phys. Commun. 183 (1) (2012) 125–146.

[54] et P. Higuera, J.L. Lara, I.J. Losada, Three-dimensional interaction of waves and porous coastal structures using OpenFOAM®. Part I: formulation and validation, Coast. Eng. 83 (2014) 243–258.

[55] et P. Higuera, J.L. Lara, I.J. Losada, Three-dimensional interaction of waves and porous coastal structures using OpenFOAM®. Part II: application, Coast. Eng. 83 (2014) 259–270.

[56] W. Sanhan, K. Vafai, N. Kammuang-Lue, et al., Numerical simulation of flattened heat pipe with double heat sources for CPU and GPU cooling application in laptop computers, J. Comput. Des. Eng. 8 (2) (2021).

[57] K. Grissa, A.M. Benselama, Z. Lataoui, et al., Investigations of the thermal performance of a cylindrical wicked heat pipe, Int. J. Energy Res. 42 (9) (2018) 3048–3058.

[58] et T. Brahim, A. Jenmi, CFD analysis of hotspots copper metal foam flat heat pipe for electronic cooling applications, Int. J. Therm. Sci. 159 (2021) 106583.

[59] et M. Ishii, T. Hibiki, Thermo-Fluid Dynamics of Two-Phase Flow, Springer Science & Business Media, 2010, et.

[60] et A. Faghri, Y. Zhang, Transport Phenomena in Multiphase Systems, Elsevier, 2006, et.

[61] F.J. Carrillo, I.C. Bourg, C. Soulaine, Multiphase flow modeling in multiscale porous media: an open-source micro-continuum approach, J. Comput. Phys. X 8 (2020) 100073.

[62] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (2) (1952) 89.

[63] et R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1) (1999) 567–603.

[64] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1) (1986) 40–65.

[65] L.S. Caretto, A.D. Gosman, S.V. Patankar, et al., Two calculation procedures for steady, three-dimensional flows with recirculation. Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics: Vol. II Problems of Fluid Mechanics, Springer Berlin Heidelberg, 1973, pp. 60–68.

[66] T. Holzmann, Mathematics, Numerics, Derivations and OpenFOAM®, Holzmann CFD, Loeben, Germany, 2016.

[67] et J. Klostermann, K. Schaake, R. Schwarze, Numerical simulation of a single rising bubble by VOF with surface compression, Int. J. Numer. Methods Fluids 71 (8) (2013) 960–982.

[68] A. Albadawi, D.B. Donoghue, A.J. Robinson, et al., Influence of surface tension implementation in volume of fluid and coupled volume of fluid with level set methods for bubble growth and detachment, Int. J. Multiph. Flow 53 (2013) 11–28.

[69] L. Wen Ho, A pressure iteration scheme for two-phase flow modeling, Multiphase transport fundamentals, reactor safety, applications 1 (1980) 407–431.

[70] et S.C.K. De Schepper, G.J. Heynderickx, G.B. Marin, Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker, Comput. Chem. Eng. 33 (1) (2009) 122–132.

[71] et A. Alizadehdakhel, M. Rahimi, A.A. Alsaairafi, CFD modeling of flow and heat transfer in a thermosyphon, Int. Commun. Heat Mass Transf. 37 (3) (2010) 312–318.

[72] et S.W.J. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, J. Comput. Phys. 160 (2) (2000) 662–682.

[73] et S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source terms, J. Comput. Phys. 227 (11) (2008) 5871–5895.

[74] et A.S. Rattner, S. Garimella, Simple mechanistically consistent formulation for volume-of-fluid based computations of condensing flows, J. Heat Transf. 136 (7) (2014) 071501.

[75] F. Rong, Z. Guo, Z. Chai, et al., A lattice Boltzmann model for axisymmetric thermal flows through porous media, Int. J. Heat Mass Transf. 53 (23-24) (2010) 5519–5527.

[76] et F. Rong, B. Shi, Incompressible lattice Boltzmann model for axisymmetric flows through porous media, Int. J. Mod. Phys. C 26 (04) (2015) 1550036.

[77] K. Grissa, R. Chaabane, Z. Lataoui, et al., Lattice Boltzmann model for incompressible axisymmetric thermal flows through porous media, Phys. Rev. E 94 (4) (2016) 043306.

[78] F. Rong, W. Zhang, B. Shi, et al., Numerical study of heat transfer enhancement in a pipe filled with porous media by axisymmetric TLB model based on GPU, Int. J. Heat Mass Transf. 70 (2014) 1040–1049.

[79] H. Hu, P. Jiang, X. Ouyang, C. Zhao, R. Xu, A modified energy equation model for flow boiling in porous media and its application to transpiration cooling at low pressures with transient effect, Int. J. Heat Mass Transf. 158 (2020) 119745.

[80] L. Shen, J. Wang, W. Dong, J. Pu, J. Peng, D. Qu, L. Chen, An experimental investigation on transpiration cooling with phase change under supersonic condition, Appl. Therm. Eng. 105 (2016) 549–556.