Modeling and Numerical Investigation of Transient Two-Phase Flow with Liquid Phase Change in Porous Media

Two-phase flow with phase change in microstructure or nanostructure is an important issue in many fronts and critical applications nowadays, but with a lack of comprehensive understanding of the mechanism. This paper numerically investigates the transient behavior of two-phase flow with liquid phase change in the porous media, which consists of a series of connected pores at micro and nanoscale with the transient form of the semi-mixed model and self-compiled programs. Transient variation and spatial distribution of structure temperature, thermal non-equilibrium characteristic, phase change location and fluid-driven pressure are obtained and analyzed, and effects of initial system temperature, structure parameter and material property on the transient behaviors of two-phase flow and fluid-structure coupling heat transfer are discussed. The numerical simulations indicate that the two-phase flow with phase change in porous media is complex and ever-changing before reaching a steady state and affected by the above-mentioned three kinds of parameters significantly. Particularly, distinct phenomena of transient heat transfer deterioration and vapor block are discovered, and it is revealed that the transient heat transfer deterioration and vapor block are more serious in a porous matrix with smaller porosity and made of materials with higher heat capacity and density.


Introduction
Two-phase flow and heat transfer at the micro and nanoscale have attracted much attention from researchers in the fields of geothermal energy recovery [1,2], fuel cell technology [3,4], building energy conservation [5,6], mobile Internet devices [7][8][9], and spacecraft thermal protection [10,11]. In transpiration cooling, which is one of the most promising heat dissipation approaches of mobile Internet devices and spacecraft [9], by virtue of porous matrix with pores at micro and nanoscale, phase change coolant exchanges heat fully with structure by flowing through, keeping the whole structure at moderate temperatures. In proton exchange membrane fuel cell, one environment-friendly power source with high energy density, back-diffusion of water from cathode to anode always appears due to the differences in pressure and concentration during the operation of fuel cells, resulting in water flooding in the micropores of anode [4]. Therefore, the two-phase flow with phase change in microstructure or nanostructure is an important issue in many fronts and critical applications nowadays.
There is much research work on the steady problems of two-phase flow with liquid phase change in microstructure or nanostructure. A great number of experiments have been conducted to explore the wide range of two-phase flow regimes in different microstructures [12][13][14]. Several widely recognized mathematical models [15][16][17][18][19][20][21][22][23][24] and a series of numerical approaches have been developed to probe into the mechanism. Based on these models and numerical approaches, the thermal non-equilibrium characteristic between fluid and solid [25], and the influences of fluid mass flow rate of fluid [20,24], heat flux distribution [26,27], and geometry configuration [19] on the two-phase flow and heat transfer characteristics have been studied.
However, there is a lack of comprehensive analysis and knowledge about the transient behavior of two-phase flow with phase change in the porous media so far. Different from the general immiscible gas-liquid two-phase flow, for example, oil-gas flow in rock layers and cracks, when phase change of liquid fluid occurs in pores and channels at micro or nanoscale, there are significant mass, momentum and energy transfers between liquid phase and vapor phase, huge amount of latent heat is released, and capillary force is generated. The resultant dramatic variations of fluid density and specific volume affect the transport characteristic of the two-phase flow and the fluid-solid coupling heat transfer greatly. Therefore, affected by the micro or nanoscales of the connected pores and the distributions of heat and pressure at structure boundaries, some special thermal responses and fluid transport characteristics would arise during the transient process. The experimental data from Huang et al. [28] indicated that there was a temporal hysteresis in the response of structure temperature due to the mismatching of the sensitivity to the variations of external heat flux and internal fluid mass flow rate. This phenomenon was confirmed by Luan et al. [29], and they further discovered that under a certain thermal environment and a range of constant fluid injection ratios, both the structure temperature and the fluid pressure fluctuated dramatically. Overall, the two-phase flow with phase change in porous media goes through an initial stage with significant variation before reaching steady state, even when the thermal environment is stable. Therefore, it is very necessary and important to develop the transient mathematical model and investigate the dynamic evolution process of phase change and the transient behavior of fluid-structure coupling heat transfer in porous media.
In this paper, the transient form of the semi-mixed model under local thermal nonequilibrium condition is proposed to describe the transient two-phase flow with phase change in porous media, and a self-compiled program is developed as the numerical approach. Through the numerical simulation, the time-varying spatial distribution of important physical parameters, such as structure temperature, thermal non-equilibrium characteristic, phase change location and fluid-driven pressure are obtained and quantitatively analyzed, and the effects of system initial temperature, microstructure parameter and material property on the transient behaviors of two-phase flow and fluid-structure coupling heat transfer are systematically discussed. Figure 1 shows the physical model for the transient problem of two-phase flow with phase change in porous media considered in this work. The homogeneous porous plate is located horizontally and has a thickness of L. Heat flux q(t) is applied on the upper surface of the porous plate, while liquid fluid at temperature T c flows into the porous plate with a mass flow rate of m c (t) from the lower surface. The liquid fluid absorbs heat from the porous matrix up to phase change, and then vapor generates and flows out of the porous plate. According to fluid states, there are three regions in the porous matrix from bottom to top: the liquid region, the mixed region of liquid and vapor phases, and the vapor region. Some mathematical models have been developed specifically to describe the characteristics of fluid flow and heat transfer in the three regions. Our previous research [16] has compared these models symmetrically, and the results indicated semi-mixed model (SMM) [19] was more accurate and flexible. Therefore, the transient form of SMM under local thermal non-equilibrium condition is proposed in this work to describe the transient two-phase flow with phase change in porous media.

Mass Conservation Equation
In the single-phase region, the mass conservation equation has uniform expression, using a subscript "i" to represent liquid phase "l" and vapor phase "v": In the two-phase region, the mass conservation equations are separately described for liquid and vapor phases by liquid saturation s:

Momentum Conservation Equation
Similarly, the momentum conservation equation has uniform expression in the singlephase region, while it is separately described for liquid and vapor phases in the twophase region.
In the single-phase region: In the two-phase region:

Energy Conservation Equation
In the single-phase region, the conservation equation of fluid energy has uniform expression: In the two-phase region, the conservation equation of fluid energy is unitedly described for the mixture of liquid and vapor phases: where T f is the temperature of the fluid, and T l = T v = T f . H f is the specific enthalpy of the mixture of liquid and vapor phases, and its definition is: It can be found that when liquid saturation s = 1, fluid is in a liquid state and H f = H l ; when 0 < s < 1, fluid is in a mixed state of liquid and vapor phases; when s = 0, fluid is in vapor state and H f = H v . Therefore H f is valid in all the single-phase and twophase regions.
In all regions, energy conservation of a solid matrix is described by:

Constitutive Relations
After solving the mixed specific enthalpy, one can obtain the saturation and temperature of the fluid by the following algebraic relations: In the single-phase region, fluid in a vapor state is regarded as an ideal gas, while in the liquid state is considered to be incompressible. In the two-phase region, fluid is the mixture of liquid and vapor phases, of which the vapor phase maintains at saturation state, and the pressure of the vapor phase depends on the temperature: Then the pressure of the liquid phase can be obtained by p l = p v − p c , where p c is the capillary pressure.
Other constitutive relations used to close the solution are summarized in Table 1.

Convective heat transfer of fluid-to-solid in pores
In the single-phase regions: In two-phase regions:

Initial and Boundary Conditions
The transient problem that a porous matrix full of liquid fluid begins to be exposed to a high heat flux is considered in this paper. Hence, the initial conditions for the whole physical domain are set as: According to the discussion results in previous research [16], the most reasonable and feasible boundary conditions at the hot side and cold side are selected and employed in this work.
At cold side: At hot side:

Solution Procedure
The transient SMM under local thermal non-equilibrium condition is solved by a selfcompiled program at University of Science and Technology of China. The finite difference method is employed to discretize the governing equations: explicit scheme [30,31] is used for the time term; second-order upwind difference, central difference and linearized schemes are, respectively, used for the convection term, the diffusion term and the source term in the energy equations. After the discretization, a set of linear algebraic equations are obtained and then solved by an iterative algorithm. The discretized equations of energy, mass and momentum are solved sequentially at each time step, after achieving the convergence, the solving procedure enters the next time step. Convergence is judged by the same criteria: the relative changes of T s , T f and p between successive iterations are less than 10 −6 . Liquid water is considered as the fluid, and its thermal properties are listed in Table 2. Table 3 presents the reference values of the structure parameters and the thermal properties of the porous matrix used in the simulations.  Table 3. Reference value of parameters used in simulations.
Heat flux at the hot side (W m −2 ) q = 1.5 × 10 6 Liquid water mass flow rate (kg m −2 s −1 ) m = 0.3 Thickness of porous media (m) h c = 31.4 Liquid water temperature at the cold side (K) T c = 300

Experiment and Verification
Due to the lack of proper transient experimental data in open literature, the mathematical model and numerical approach proposed in this work are validated with the steady-state experimental data from Hu et al. [32]. In the experiment, the porous cylindrical test section was consisted of sintered bronze with a porosity of 0.318 and a particle diameter of 0.2 mm and had a diameter of 40 mm and a height of 100 mm. An average heat flux 0.21 MW/m 2 was applied to the upper surface of the porous test section, and liquid water with an initial temperature of 300 K was injected from the lower surface. The wall temperature along the flow direction was measured. Simulations under the completely identical operating condition with the experiments were carried out with our mathematical model and numerical approach, and the comparison results are shown in Figure 2. The maximum relative deviation between the experimental and simulation results is less than 4%, which demonstrates the mathematical model and numerical approach suggested are valid to a certain extent.

Transient Behavior of Two-Phase Flow with Phase Change in Porous Media
The two-phase flow with phase change in porous media will eventually reach a steady state when the thermal environment and fluid inject rate are unchanged. However, before that, the temperature of the porous structure, the location of phase change and the driven pressure of fluid change acutely with time, and this stage are known as initial stage. Figure 3 shows the transient behavior of the porous system with an initial temperature of 300 K when subjected to a constant heat flux of 1.5 × 10 6 W/m 2 and cooled by liquid water with a mass flow rate of 0.3 kg/m 2 s. It can be found that the fluid temperature at the hot side and the location of the interface between two-phase region and vapor region approach certain values at approximately 300 s and then maintain unchanged, which means the whole system reach the steady state by this time. From the beginning t = 0 s to t = 300 s, the system goes through an initial stage, in which physical parameters vary with time significantly. To study the two-phase fluid flow and fluid-solid heat transfer characteristics in the initial stage, Figure 4 exhibits the distributions of the temperature, pressure and saturation of fluid within porous matrix at t = 1, 5, 25 and 100 s. As shown in Figure 4, the spatial distributions of fluid temperature and pressure have same variation trends at different times: the fluid temperature rises slowly in the liquid region, keeps almost constant in the two-phase region, and increases sharply in the vapor region; the pressure drops gradually in the liquid region, rises rapidly in the two-phase region, and decreases substantially in vapor region. Conforming to the same flow and heat transfer mechanisms, the spatial distribution laws of physical parameters obtained at different times in the initial stage are reasonable, and identical with those in the steady state obtained by steady SMM [16].
Combining Figure 3; Figure 4, it could be found that: the porous structure is full of liquid water at first, and the temperature and pressure of fluid are low; as time goes by, phase change occurs, the location of phase change moves down and vapor generates, and then the temperature and pressure of fluid increase monotonically; finally, there are three regions, i.e., vapor region, two-phase region and liquid region, sequentially formed within the porous structure. We find that the maximum temperature and pressure always appear at either the hot side or cold side. Moreover, the maximum temperature and fluid-driven pressure of the two-phase flow with phase change in porous media are most concerned in engineering applications. Therefore, we focus on the transient variation of physical parameters on the boundary in the following study. Figure 5 exhibits the transient variation of the physical parameters on the boundary obtained under different system initial temperatures when other operating parameters remain to be the same with the reference value in Table 3. As shown in Figure 5, the solid temperature at the hot side always tends to the same value when the whole system reaches steady state, regardless of the initial temperature. Similar numerical results can be obtained for the temperature difference between solid and fluid at the hot side, phase change location and fluid pressure at inlet. This is reasonable, because the steady state values of the physical parameters depend on the operation condition rather than the initial condition. For the porous matrix with fixed microstructure, when the heat flux on the boundary and the fluid mass flow rate both keep unchanged, the heat transfer between the solid matrix and the fluid will eventually achieve same equilibrium, although may spend different length of time. As the results, the spatial distributions of all the physical parameters display same values at the equilibrium or steady state.  Combining Figure 3; Figure 4, it could be found that: the porous structure is full of liquid water at first, and the temperature and pressure of fluid are low; as time goes by, phase change occurs, the location of phase change moves down and vapor generates, and then the temperature and pressure of fluid increase monotonically; finally, there are three regions, i.e., vapor region, two-phase region and liquid region, sequentially formed within the porous structure. We find that the maximum temperature and pressure always appear at either the hot side or cold side. Moreover, the maximum temperature and fluid-driven pressure of the two-phase flow with phase change in porous media are most concerned in engineering applications. Therefore, we focus on the transient variation of physical parameters on the boundary in the following study. Figure 5 exhibits the transient variation of the physical parameters on the boundary obtained under different system initial temperatures when other operating parameters remain to be the same with the reference value in Table 3. As shown in Figure 5, the solid temperature at the hot side always tends to the same value when the whole system reaches steady state, regardless of the initial temperature. Similar numerical results can be obtained for the temperature difference between solid and fluid at the hot side, phase change location and fluid pressure at inlet. This is reasonable, because the steady state values of the physical parameters depend on the operation condition rather than the initial condition. For the porous matrix with fixed microstructure, when the heat flux on the boundary and the fluid mass flow rate both keep unchanged, the heat transfer between the solid matrix and the fluid will eventually achieve same equilibrium, although may spend different length of time. As the results, the spatial distributions of all the physical However, the initial temperature of the porous system has a great influence on the transient variation of the physical parameters in the initial stage. As shown in Figure 5a, the solid temperature presents quite different variation trends in the initial stage when the system initial temperature is different. The solid temperature increases to the steady value with a low initial temperature of 300 K while briefly exceeds the steady value in the initial stage with higher initial temperatures of 400 K and 500 K. These results are very important because the temperature of the solid matrix is usually concerned most by the engineers and designers in the industry. Supposing that the melting point of the porous material is 1400 • C the steady solid temperatures obtained under the three initial conditions are the same and lower than this melting temperature. This indicates the fluid amount applied could ensure no ablation in the steady state under the operation condition considered here. However, when the system initial temperature is 400 K and 500 K, the temperature of the solid matrix has exceeded the melting point before reaching the steady state, i.e., ablation has already occurred; only when the system initial temperature is 300 K, the temperature of the solid matrix is always lower than the melting point. This is a distinct phenomenon of transient heat transfer deterioration, which alerts us that: it is not enough to only ensure the steady temperature of the whole system lower than the critical value, because ablation may occur before the system reaches steady state. Therefore, the coolant injection amount should be increased in the initial stage for the case that the system initial temperature is relatively high.  Figure 5b,c depict the effect of system initial temperature on the transient variation of phase interface location within porous structure and thermal non-equilibrium characteristic at the hot side. The thermal non-equilibrium characteristic is greatly affected by the magnitude of phase change and the phase distribution of fluid. According to the discussion in steady state [16], the temperature difference between solid and fluid increases rapidly at the beginning of the two-phase region and reaches the maximum at the interface between the two-phase and vapor regions. Therefore, the thermal non-equilibrium characteristic at the hot side is particularly significant at the beginning of the transient process of the two-phase flow with liquid phase change in porous media with an initial temperature of 300 K. When the porous structure is subjected to high heat flux, phase change briefly appears on the boundary. Solid temperature at the hot side increases rapidly, while fluid temperature there almost keeps constant due to latent heat, resulting in a locally remarkable temperature difference. When the system initial temperature is 400 K, the time of fluid phase change on the boundary shortens, and consequentially the thermal non-equilibrium characteristic there is not so strong. As the system, initial temperature increases to 500 K, fluid has no time to flow to the boundary and phase change occurs midway, so the thermal non-equilibrium characteristic at the hot side is weak.

Effect of Initial Temperature
The magnitude of phase change and phase distribution of fluid also has great influence on the driven pressure of the fluid. From Figure 5c, it could be found that: when T 0 = 300 K, the area of vapor region in the initial stage is always smaller than that in the steady state, but when T 0 = 400 K and 500 K, the area of vapor region in the initial stage expands quickly and once exceeds the steady value. For the flow resistance of vapor is much higher than that of liquid, larger vapor region inevitably corresponds to larger fluid pressure at inlet. As shown in Figure 5d, when T 0 = 400 K and 500 K, the fluid pressure at inlet in the initial stage likewise exceeds the steady value. This is a distinct phenomenon of transient vapor block, which is worthy of attention. The fluid-driven pressure is an extremely important parameter in the design of a two-phase flow system and is usually set according to steadystate operating conditions. However, based on the results obtained here, the driven force of fluid required in the initial stage could actually be much greater than that in the steady state. If lacking driven force in the initial stage, the fluid supply must break down, so in this situation, the cooling system may already fail in the initial stage. Figure 6 exhibits the effect of porosity on the transient variations of solid temperature at the hot side and fluid pressure at inlet. The system initial temperature is set as 500 K to embody the distinct phenomena of transient heat transfer deterioration and vapor block. From Figure 6a, the maximum of solid temperature at the hot side could be observed in the initial stage, i.e., transient heat transfer deterioration, which is consistent with the result in Figure 5. With the increase of porosity, the solid temperature at the hot side increases while the transient heat transfer deterioration reduces. As listed in Table 1, the specific surface of a pore is linear function of the porosity, i.e., α = 6(1−ε)/d p , which indicates that fluid-to-solid heat transfer coefficient decreases with the porosity. Hence, more heat is stored in the solid matrix with larger porosity, resulting in a higher solid temperature. From Figure 6b, it could be found that the fluid pressure at inlet and the transient vapor block effect both decrease with the porosity. This is because solid matrix with larger porosity has smaller effective thermal conductivity k s.eff . As the results, the temperature gradient within the solid matrix becomes larger and vapor region becomes smaller. Accordingly, the resistance of fluid flowing through the vapor region reduces, i.e., the driven pressure decreases. It should be particularly noted that the porosity not only affects the transient variation of physical parameters in the initial stage but also changes their steady values. Large porosity design could greatly reduce the fluid-driven force and impair the transient heat transfer deterioration and vapor block, but at the cost of increasing the system temperature. Therefore, porosity is an important structure parameter of the two-phase flow with liquid phase change in porous media and should be given special consideration in the practice of engineering design.  Figure 7 shows the effect of material heat capacity on the transient variations of solid temperature at the hot side and fluid pressure at inlet. The system's initial temperature is 500 K, while other operating parameters remain to be the same with the reference values in Table 3. As the heat capacity of the solid matrix increases, the system spends more time reaching steady state, and the transient heat transfer deterioration in the initial stage is severer. This is reasonable because the larger the heat capacity of the solid matrix, the slower the variation of solid temperature, and accordingly, the longer the time required by the system to reach a steady state. On the other hand, solid matrix with larger heat capacity stores more heat, which is hard to release in a short time and inevitably causes severer transient heat transfer deterioration. The maximum fluid pressure in Figure 7b shows the same variation trend as that of solid temperature in Figure 7a, i.e., the transient vapor block effect become more significant with the increase of material heat capacity. After all, a higher solid temperature on the hot side corresponds to a larger vapor region and thereby greater flow resistance. As a result, the fluid pressure at inlet is higher.  Figure 8 shows the effect of material density on the transient variations of solid temperature at the hot side and fluid pressure at inlet. For the similar reason, as the density of the solid matrix increases, the system spends more time reaching steady state, and the transient heat transfer deterioration and the transient vapor block effect in the initial stage both become more significant. Overall, material properties have significant influence on the initial stage of the two-phase flow with liquid phase change in porous media, but no influence on the steady value of structure temperature and fluid-driven pressure. Particularly, the transient heat transfer deterioration and the transient vapor block effect are reduced when a material with low heat capacity and density is employed. Therefore, this provides an approach for the engineer and designer to impair the transient heat transfer deterioration and vapor block and avoid the corresponding negative effects.

Conclusions
This work is aimed at providing a comprehensive understanding of transient twophase flow with a liquid phase change at the micro and nanoscale. A new mathematical model and corresponding numerical approach are developed, through which the dynamic evolution of phase change and transient behavior of fluid-structure coupling heat transfer are analyzed and the influences of initial system temperature, structure parameter and material property are discussed. From the simulations and analysis, the following conclusions can be drawn: (1) For the porous matrix with fixed microstructure, when the heat flux on the boundary and the fluid mass flow rate both keep unchanged, the two-phase flow with liquid phase change in porous media will eventually achieve equilibrium, which is known as the steady state. Before reaching steady state, the transient two-phase flow with liquid phase change in porous media goes through an initial stage, in which the physical parameters vary with time significantly; (2) The system initial temperature has no influence on the steady state of the two-phase flow with liquid phase change in porous media. However, when the initial system temperature is relatively high, distinct transient heat transfer deterioration and vapor block occur in the initial stage, which may cause unexpectedly premature cooling failure and structure ablation; (3) The porosity not only affects the transient variation of the physical parameter in the initial stage but also changes their steady values. Large porosity design could greatly reduce the fluid-driven force and impair the transient heat transfer deterioration and vapor block, but at the cost of increasing the system temperature; (4) As the heat capacity and the density of the solid material increases, the system spends more time reaching a steady state, and the transient heat transfer deterioration and vapor block effect in the initial stage are severer. Therefore, porous matrix made of a material with low heat capacity and density could reduce the transient heat transfer deterioration and the transient vapor block effect without affecting the steady state.