Numerical investigation of boiling and forced convection heat transfer in inclined porous enclosure using modified enthalpy formulation

Two-Phase flow in an inclined rectangular porous media, under unsteady-state condition, has been numerically investigated in this article, based on the modified h-formulation of Two-Phase Mixture Model (TPMM). The governing equations have been discretised using Finite Volume Method (FVM) and solved iteratively in a SIMPLE-like manner. Effects of various parameters on the flow and temperature fields have been investigated, which clearly demonstrated that the inclination angle strongly affect the boiling initiation. Recirculating flow has been observed when the inclination angle  0   . The results clearly indicated that the operating conditions and the porous medium properties have significant effects on the initiation and termination of phase change process. A closer inspection of the results reveals that the presence of a critical inclination angle  depending on the value of * K and l Re , which correspond to a maximum values of vapour volume. The modified h-formulation requires significantly less computational time as compared with the existing H-formulation of TPMM.


Nomenclature
Body force vector per unit mass, Two-phase flow and analysis of boiling in porous media is an important topic, which spans a broad spectrum of engineering disciplines [1][2][3].Based on Two-Phase Mixture Model (TPMM) proposed by Wang and coauthors [4 -6], which consists of mass, momentum and energy conservation equations in the modified form, three distinct regions can be observed, if a liquid is injected into a porous evaporator.They are: sub-cooled liquid region, two-phase region and superheated gaseous region.Owing to the generality and the simplicity of the numerical formulation, TPMM of Wang [6] has been most often used for the simulation of two-phase flow within porous media [6][7][8][9][10][11][12][13][14][15][16][17][18][19].
However, a major problem that is encountered when applied TPMM [6] during the simulation of two-phase flow within porous media is the occurrence of drastic, nonphysical change in the predicted temperature distribution over an extremely short distance that results primarily due to the presence of discontinuity in the effective diffusion coefficient across the interfaces.In order to avoid the aforementioned problem, a new smoothing algorithm has been developed by Alomar et al [20].The proposed algorithm has been also employed by Alomar et al. [21][22][23][24] for the successful prediction of complete phase change process of liquid water inside a divergent [21], an annular [22] and constant cross-sectional area [23] porous evaporator as well as porous channel [24] using both Local Thermal Equilibrium (LTE) and Non-Equilibrium (LTNE) conditions.
After critically analysed the drawbacks of the original enthalpy (h) formulation [3], a modified formulation has been developed by Ray and Alomar [25] that can easily accommodate substantial density variations in the single phase regions as compared to the popular modified Hformulation [6].Most importantly, the modified hformulation [25] does not require any artificial definition of the modified volumetric (H) enthalpy as the dependent variable for the energy conservation equation [5,6] while successfully eliminating the identified problems.It has been demonstrated that the modified h-formulation requires significantly less computational time as compared to the Hformulation [6] and hence it is recommended for the simulation of two-phase flow inside porous media.Later, the modified h-formulation [25] has been employed by Alomar et al. [26] in order to test the performance of nonstaggered and staggered grid layouts for the simulation of complete evaporation process inside porous evaporator.The results indicated that the solutions obtained by using non-staggered grid layout are identical with that obtained using staggered grid layout.After that Alomar et al. [27] performed a numerical study on the complete evaporation process inside anisotropic porous media using the same formulation.The results indicated that the anisotropy of permeability and solid thermal conductivity have significant effect on the complete evaporation process.
Another interesting feature of boiling and force convection effects is the influence of porous evaporator inclination.As a continuation of the previous research efforts toward a complete understanding of boiling and force convection in porous media, the two-phase flow through a porous enclosure with a finite wall heat sources has been numerically investigated in the present article to include the effect of inclination angle of the evaporator.The mathematical formulation is based on the modified hformulation [25] along with LTE model.

Problem Description
Schematic representation of the problem geometry is shown in Fig. 1 along with its dimensions.The problem under investigation is incomplete phase change process of liquid water in an inclined porous channel, under unsteadystate condition, which has been considered as a demonstrative example.The model has been applied for a two-parallel plate channel ( W L ), filled with a homogeneous as well as isotropic porous medium of known properties.The channel has unheated inlet and exit lengths of i l and e l located on the bottom wall, respectively.The length of heated section h l , where an asymmetric heat source * w q    has been applied on the bottom wall of the channel, whereas the inlet and exit lengths are considered well insulated.On the other hand, the upper plate of the channel is considered at constant temperature.An external pressure gradient drives the sub-cooled water at the inlet, with a low temperature and uniform velocity in u , to flow through the channel.The flowing water is heated as it passes through the heated surface of the channel.When the heat flux is sufficiently increased, boiling occurs at the heated section and thus a two-phase zone is formed.For numerical simulations, the modified h-formulation [25] along with LTE model has been used.Further, all the properties have been assumed to be constant for both liquid and vapour phases.

Governing Conservation Equations
In the present study, the LTE model has been employed for the present simulations 2 .The existing model considers unsteady, two-dimensional incomplete evaporation process within an isotropic porous medium due to localized asymmetric heating in an inclined porous channel.In order to construct a modified h-formulation [25] applicable under very low flow velocity, Darcy's law is introduced.All applicable equations have been made dimensionless using the definitions in Table 1 in order to reduce the number of parameters those affect final results 3 .By adopting the modified h-formulation [25], the conservation equations of mass, momentum and energy by including the transient term and the effects of inclination angle are given as: (2) 2 This model assumes that the solid and the fluid phases are in thermal equilibrium with each other and they coexist at the same temperature. 3The dimensionless parameters have been carefully chosen such that they retain all governing equations along with expressions for the mixture variables in the same form of modified h-formulation [25].
where symbols with superscript '*' denote the dimensionless values and h is the enthalpy.Further definitions of dimensionless mixture properties in the twophase region, appearing in Eq. ( 1), (2), and (3) are provided in Tables 1 and 2. Table 1.Definitions of Dimensionless Variables [20,25].
It may be noted here that according to TPMM, the momentum balance given in Eq. ( 2) considers only the body force term and Darcy's law [4][5][6] for the additional pressure drop due to the presence of a porous matrix.Since Darcy's law is used to model the flow in porous media, Eq. ( 2) is valid only for low mass flow rate applications 4 , where the fluid Reynolds number based on the permeability of the porous medium, , is less than one [17].In Eqs. ( 2) and (3), the permeability of the porous medium ) can be evaluated either from a separate detailed study or by using empirical correlation where K is expressed as a function of the porosity  and average pore diameter of porous matrix [28,29].In Eq. 4 Since the convective and the diffusive effects have been neglected.
(2), * k  is the dimensionless kinetic density, which according to the Boussinesq approximation and it could be expressed using the linear expansions of single phase densities as a function of local temperature around their saturation temperature [6]: Additionally,
The expression of * h  in Eq. ( 3) is presented in Table 2, which contains the effective thermal conductivity of the combined medium * eff k in its dimensionless form.A simple estimation for * eff k can be obtained from a separate numerical and/or experimental study [30 -33].In the present investigation, it has been evaluated according to the parallel arrangement model as [6,20]: where * l k , * v k and * s k are the thermal conductive of liquid, vapour and solid, respectively.The velocities of individual phases can be evaluated using the following relations [6]: where * j is the total diffusive flux, which can be written as: The combined temperature * T in the single phase region and the liquid saturation s in the two-phase zone can be calculated from the solution of * h as shown in Table 3.It is evident from Table 1 that the dimensional temperatures can be retrieved from their dimensionless values using only the properties of working fluid [20].Table 3. Supplementary relations between enthalpy, temperature, and liquid saturation [25].

Boundary Conditions
For the present investigation, the following conditions have been employed at different boundaries of the computational domain in order to solve Eqs. ( 1 Therefore, the prescribed boundary conditions at the inlet may be summarized as: iii) At the evaporator outlet, i.e., at , the axial diffusion in Eq. ( 3) has been neglected and hence the second derivatives of * h in the flow direction has been set zero 6 .Therefore, * h at the outlet has been obtained by linearly extrapolating the values at the neighbouring interior nodes that ensure 5 Since sub-cooled liquid enters the channel, 1 *  in  has been set. 6For shorter unheated exit length, or even in its absence, as long as the flow is assumed to be thermally fully-developed for constant wall heat flux, the present condition would always hold, indicating linear rise in temperature, although the latter condition may never be met.

Numerical Solution
The Finite Volume Method (FVM) on staggered grid layout has been used in order to discretise the conservation equations ( 1), (2), and (3).According to FVM, the entire computational domain is first divided into a number of nonoverlapping control volumes.In the present study, only uniform control volumes have been used for all considered cases.The required pressure gradient and the kinetic density at the cell faces have been obtained by employing the second order Central Differencing Scheme (CDS) and 7 This condition imply that the gradient of * h in the axial direction is assumed to remain constant at the exit, which, would allow finite heat transfer at the exit.As will be shortly apparent from the predicted results, owing to the presence of unheated exit length, heat transfer at the exit of the evaporator turns out to be negligible and hence by linearly interpolating the adjacent nodal values, respectively.For the energy conservation, the first order Upwind Differencing Scheme (UDS) and second order CDS have been used for expressing the convective terms and the diffusive terms in Eq. (3), respectively.Discretised equations for all the nodes, along with boundary conditions, exhibit a tri-diagonal matrix structure and hence they have been solved by the Thomas algorithm [34].In addition, diffusivities for the enthalpy * h  at cell faces have been obtained using the harmonic mean approximation, which ensures the balance of diffusive energy flux.On the other hand, the advection correction coefficient h  and the hindrance coefficient f at the cell faces have been obtained using linear interpolation from the adjacent nodal values.
In the present study, the smoothing algorithm proposed by Alomar et al. [20,35] along with the recommended smoothing parameters has been employed for * h  to avoid the non-physical jump in the predicted temperature distributions 10 and to obtain the converged solutions for all considered cases.A grid independence study has been carried out and the results show that a mesh of 300 180   CV N control volumes produce a grid-independent solutions.In this regards, all subsequent have been performed using this mesh size.Since the mixture properties, listed in .Owing to the strong non-linearity in the governing conservation equations and the related constitutive relations arising out of the interdependence of dependent variable and since during discretisation, the discretised equations have been further under-relaxed and a relaxation factor of 1 .0 has been generally used for most of the simulations in order to achieve convergence with the convergence criterion set to 5  10  .

Results and Discussion
For all studied cases, the length and height of the evaporator/channel have been kept fixed at ) 11 .Therefore, the length of the heated section on the bottom plate has been taken as ). Sub-cooled liquid water has been used as the working fluid and its properties are given in Table 4, where the saturation temperature sat T of water 10 The non-physical jump occurs primarily due to the discontinuity in the effective diffusion coefficient across the saturated liquid condition. 11The exit length has been carefully chosen in order to avoid the influence of exit boundary condition on the internal solution.
has been considered to be C 100 o [20,36,37].Water enters the domain with uniform velocity and temperature For the porous medium, the porosity has been assumed to be around 0.35.For the present demonstrative examples, the permeability of the porous medium has been assumed to follow the Carman-Kozeny type relation [28,29].Based on this, the permeability has been varied from 879 g , respectively, assuming a porous matrix made of glass bead system [36,37].At the heated section, it has been assumed that heat is added from an external source, supplying a constant heat flux ).In addition, the acceleration due to gravity vector has been assumed to act in both negative xand y-directions and hence   .Wang [6] denoted the interface that separates the sub-cooled liquid region from the two-phase zone as a condensation front.The condensation front is appreciably inclined toward the outlet due to incoming flow of sub-cooled liquid.In order to identify the interface separating the sub-cooled liquid and the two-phase regions, isotherm is presented for . This isotherm is, however, marked as that corresponding to C 100 o  sat T in all relevant figures, where the negative sign refer to the saturation isotherm from the sub-cooled liquid zone.It is evident from the Figs. 2 and 3 that the two-phase zone is extended in both axial and transverse directions as compared with the horizontal case due to more energy added to the domain associated with increasing the buoyancy force.The buoyancy force in the transverse direction increases and reaches its maximum value at  0   and then the condensation front is moved considerably towards the cooled wall.Most importantly, the axial and transverse diffusions remain present even in the two-phase region as there exists a finite gradient of * h , although the temperature remains unchanged at saturation condition.Since * h increases upward and in the downstream direction due to heating and buoyancy effect, energy is transferred in the backward direction owing to the axial diffusion [20].The liquid velocity indicates parallel flow of the subcooled liquid at the inlet but deflected movement towards the upper surface of the channel when the liquid approaches the condensation front, so as to bypass the vapour region surrounding the bottom heated surface, as shown in Fig. 4 for , recirculating flow has been observed as shown in Fig. 4. For inclined and vertical positions, ascending forced flow assists buoyancy-induced flow near the heated wall and the heated zone is more extended and is stretched toward the outlet.When  0   , the buoyancy force in the longitudinal direction becomes higher and evacuated heat from the heat source is increased.The vapour velocity vectors presented in Fig. 5 show a primarily upward movement upon the vapour generated on the heated surface due to local buoyancy.The max , K Re for liquid and vapour phases is evaluated according to the maximum velocity in the respective phase in order to ensure that the assumption of Darcy flow is still valid.
In order to further demonstrate the effect of inclination angle of the channel, results have been generated considering a variation in l Re and * K .These results are summarized in Figs. 6 -9.The effect of * K on the temperature, liquid saturation, velocity vectors for both liquid and vapour phases is presented in Fig. 6 when   In order to demonstrate the usefulness of the modified h-formulation [25] for multi-dimensional problems, under unsteady-state condition, the solutions of two-phase flow obtained using the current formulation has been compared with those obtained using the H-formulation [9] under the same conditions.The comparison is presented in terms of the evaporated volume fraction ev as a function of the inclination angle, as shown in Fig. 10.The ev has been used to describe the importance of boiling and it has been defined in dimensional form as [9] 12 : 12 Najjari and Nasrallah [9] pointed out that the evaporated volume fraction is a direct measure of the boiling strength and the extent of the resulting two-phase and vapour regions.The average absolute difference in the evaporated volume fraction has been found to be 0.05%.Particularly for these cases, the modified h-formulation requires approximately 45% less computational time as compared to the existing H-formulation [6] and hence it is recommended for the simulation of the two-phase flow inside porous media.

Conclusions
Numerical simulations of two-phase flow of water inside an inclined porous evaporator, based on the modified hformulation [25], have been investigated in this study in order to assess the effect of inclination angle of the evaporator.The flow has been assumed to be unsteady as well as two-dimensional.All the governing equations have been discretised using the Finite Volume Method (FVM) and the resultant coupled set of discretised equations have been solved iteratively using SIMPLE-like algorithm [34].Numerical simulations of incomplete phase change process have been obtained for similar geometry with identical heat flux value of Najjari and Nasrallah [9] for the purpose of comparison.The effects of inclination angle, inlet Reynolds number, Darcy number and total evaporated volume fraction have been investigated.The results obtained from modified h-and H-formulation [6] have been compared.The major conclusions from the present study may be summarised as follows: 1.The results clearly demonstrated that the boiling depends strongly on inclination angle of the channel.

Fig. 1 .
Fig. 1.Schematic representation of the phase change problem.
since the phase change process in an inclined porous media has been considered.

,
) -(3): i) At the inlet 0 *  t , the porous medium has been assumed to be at constant temperature * in T and the liquid velocity is uniform: the porous medium has been assumed to be saturated with sub-cooled liquid water that enters the evaporator with uniform velocity * in u and specified constant temperature of the fluid *

7 .
On the other hand, the flow has been assumed to be hydro-dynamically fully-developed at the exit8 .Therefore, the exit boundary condition for the flow may be written as: imposed9 .As far as the energy conservation equation is concerned, the specified heat flux boundary condition has been employed at the walls, where a vanishing local heat flux has been enforced at the unheated segments, i.e.
imposed and the wall has been kept at constant temperature * in T .
all tested cases.
study.The inlet velocity has been carefully chosen in order to satisfy the criterion of the fluid Reynolds number 1  K Re for the velocity field of both liquid and vapour phases.With these inputs and the properties of water in Table4 other hand, the properties of solid phase, i.e., thermal conductivity for all the cases.The inclination angle of the channel has been considered to vary from  0 to  90 for all cases studied.The gravitational Reynolds number g Re , required for the calculation of the Froude number Fr , and the normalized surface tension  ~ have been taken as for the calculation of relative permeabilities, n=1 has been adopted without any exception.

.
It is evident from Fig. 6 that * K has significant effects on temperature distribution, particularly in the two-phase region.An increase in * K leads to a decrease in the size of two-phase region due to increase the effective diffusion coefficient in this region.As Conferences 240, 01001 (2018) https://doi.org/10.1051/matecconf/201824001001ICCHMT 2018 * K is decreased, the two-phase region initiates earlier and the size of this region is extended in the axial and transverse directions, as may be observed from Figs. 6(a) and 6(b).An early initiation of phase change processes for explained by the reduced capillary diffusion affecting the effective diffusion coefficient in the two-phase region.As * K is increased, the volume of vapour is decreased by increasing  and the recirculating flow is increased, as shown in Figs. 6 -9.
mentioned here that all reported computations have been carried out on a desktop personal computer (Intel Core i7 2.80 GHz CPU with 4 cores, each with 4 GB DDR3 1333 MHz RAM made by Hynix).The typically required computation time on this machine for the present study has been recorded to vary from 1 hour to 2 hours.It is evident from Fig.10that the results obtained from modified h-formulation provide an excellent agreement with that obtained by using H-formulation[9] for all tested cases.It is clear from the figure that the volume of vapour MATEC Web of Conferences 240, 01001 (2018) https://doi.org/10.1051/matecconf/201824001001ICCHMT 2018 decreases for a higher * K and l Re as well as for an increasing inclination angle  .It has been found that the h- formulation requires significantly less computation time.

Fig. 10 .
Fig. 10.Comparison of evaporated volume fraction determined in the present study and results of Najjari and Nasrallah [9] for various  at different * K and

2 .
Recirculating flow has been observed except when the channel in a horizontal position.3. The results also indicated that the presence of a critical inclination angle  depending on the value of * K and l Re which correspond to a maximum values of vapour volume.The results show that l Re and * K have a strong effect on the initiation and termination of phase change process.
3 m kg