Simulation of liquid–vapour phase change process inside porous media using modified enthalpy formulation

International Journal of Thermal Sciences 105 (2016) 123e136

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Simulation of liquidevapour phase change process inside porous media using modified enthalpy formulation Subhashis Ray a, *, Omar Rafae Alomar a, b, * a €t Bergakademie Freiberg, Gustav-Zeuner-Straße 7, D-09596, Institute of Thermal Engineering, Chair of Gas and Heat Technology, Technische Universita Freiberg, Germany b Technical College of Mosul, Cultural Group Street, Mosul, Iraq

a r t i c l e i n f o

a b s t r a c t

Article history: Received 1 December 2015 Received in revised form 24 February 2016 Accepted 26 February 2016

In the present article, after critically analysing the drawbacks of the existing h-formulation based on Two Phase Mixture Model (TPMM), a modified formulation has been developed that can easily accommodate substantial density variations in the single phase regions. The results for steady one-dimensional complete liquidevapour phase change problems of water inside a porous evaporator, obtained from the modified h- and the existing H-formulations, have been compared and excellent agreements have been observed for all tested cases. It has been also observed that the modified h-formulation requires significantly less computation time, although all variants of TPMM requires smoothing of effective diffusion coefficient in order to avoid “jumps” in the predicted properties. Since the proposed formulation does not require the definition of any artificial variable and in view of the identified advantages, the method is strongly recommended for the future use. Nevertheless, it should now be extended in order to accommodate the local thermal non-equilibrium condition and multi-dimensional phase change problems in presence of substantial density variations in the single phase regions. © 2016 Elsevier Masson SAS. All rights reserved.

Keywords: Two-phase mixture model Modified enthalpy formulation Porous evaporator Complete phase change Variable density Smoothing of diffusion coefficient

1. Introduction Liquidevapour phase change processes inside porous media are encountered in diversified scientific and engineering applications [1e3]. Numerical simulations of such problems essentially rely upon the traditional Separate Flow Model (SFM) [4,5], where the liquid and the vapour phases are considered as two distinct fluids with different properties, that generally assume different velocities, satisfying separate sets of conservation equations. In this respect, the complete phase change process within porous media could be characterised by the presence of three distinct regions: sub-cooled, two-phase mixture and superheated vapour regions. Since SFM is most often extremely complex and perhaps inconvenient for the direct use in numerical simulations, based on such classification, Ramesh and Torrance [6e8] proposed the Separate Phase Model (SPM) in order to solve problems involving boiling and natural

* Corresponding authors. Institute of Thermal Engineering, Chair of Gas and Heat Technology, Technische Universit€ at Bergakademie Freiberg, Gustav-Zeuner-Straße 7, D-09596, Freiberg, Germany. E-mail addresses: [email protected], [email protected] (S. Ray), [email protected] (O.R. Alomar). http://dx.doi.org/10.1016/j.ijthermalsci.2016.02.014 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.

convection inside a fluid-saturated porous medium for each of these regions while continuously tracking the interface between them using a moving boundary approach. For practical applications, however, SPM still remains inconvenient for numerical implementation owing to the requirement of interface tracking and the presence of large number of coupled nonlinear equations. Recognising the complexities and the inherent problems of SFM, Wang and Beckermann [9] proposed an enthalpy-based Two-Phase Mixture Model (TPMM) where the different phases are viewed as the constituents of a binary mixture. As explained later in section 2, Wang et al. [10] first identified that the original enthalpy (h-) formulation [9] is also not readily suitable for numerical implementation owing to several reasons. Therefore, they proposed the first modification to the original model and presented a study on boiling and natural convection in a porous layer heated from below. Later, Peterson and Chang [11] employed this modification [10], along with the assumption of Local Thermal Non-Equilibrium (LTNE) condition in order to investigate the phase change process from the sub-cooled liquid state to the saturated mixture inside a highly conductive porous channel. However, a closer look into the first modification of TPMM [10], reveals that since some of the mixture variables and properties

124

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Nomenclature Ac al2,av2 b ~ b Cp D f Fr g h hco hfg H j J k krl, krv K l L n Nuo p P Pe Pr 00 q_ Q* R Re Reg s

cross-sectional area of the pipe per unit radian ¼ R2/2, m2 second order coefficients for liquid and vapour densities, respectively body force vector per unit mass, m/s2 normalised body force vector per unit mass ¼ b/g specific heat at constant pressure, J/kgK capillary diffusion coefficient for enthalpy, kg/ms hindrance function pffiffiffiffiffiffi Froude number ¼ uin/ug¼ uin = gR acceleration due to gravity, m/s2 specific enthalpy, J/kg convective heat transfer coefficient at the outer surface of the pipe, W/m2K latent heat of vaporization ¼ hv,sat  hl,sat, J/kg modified volumetric enthalpy ¼ r(h  2hv,sat), J/m3 diffusive mass flux vector, kg/m2s capillary pressure function thermal conductivity, W/mK relative permeabilities for liquid and vapour, respectively permeability of porous media, m2 length of individual segments, m length of the porous pipe, m exponent of saturation in the expressions for relative permeabilities Nusselt number based on liquid properties ¼ hco(2R)/kl effective pressure, Pa perimeter of the pipe per unit radian ¼ R, m Peclet number ¼ uinR/a ¼ RePr Prandtl number ¼ n/a constant heat flux, W/m2 00 normalised heat flux, ¼ q_ R=mhfg pipe radius, m Reynolds number ¼ uinR/n gravitational Reynolds number ¼ ugR/nl liquid saturation

remain undefined for the superheated vapour phase, it cannot be employed for the prediction of complete phase change process.1 In order to eliminate this drawback, Wang [12] proposed the second modification of TPMM [9], by introducing the modified volumetric enthalpy H ¼ r(h  2hv,sat) as the dependent variable for the energy conservation equation.2 Nevertheless, all versions of TPMM are characterised by the coexistence of a two-phase zone, separated by the irregular as well as moving interfaces (for transient problems) of single-phase (either liquid or vapour) regions. Since these variants are obtained for fixed numerical grid, unlike SFM (or, SPM), there is no requirement for the complex as well as time consuming interface tracking. Possibly owing to the generality and the ease of numerical implementation, the H-formulation [12] is by far the most widely used method for the simulation of phase change process inside porous media [13e20]. Although this formulation [12] is applicable only under the assumption of Local Thermal Equilibrium (LTE) condition, it has been extended further in order

1 Where phase change takes place from the sub-cooled liquid phase to the superheated vapour phase. 2 Henceforth shall be referred to as the H-formulation.

t T u ug x

time, s temperature, K velocity vector, m/s pffiffiffiffiffiffi gravitational velocity ¼ gR, m/s axial coordinates, m

Greek symbols a thermal diffusivity ¼ k/rCp, m2/s b isobaric expansion coefficient, K1 DT temperature difference over which Ghis relaxed in the single-phase regions, K Dr density difference between phases ¼ (rl  rv), kg/m3 g advection correction coefficient G diffusion coefficient for enthalpy, kg/m s (for h) or m2/s (for H) ε porosity of the medium l relative mobility m dynamic viscosity, kg/ms n kinematic viscosity, m2/s r density, kg/m3 s surface tension coefficient, N/m ~ s normalised surface tension coefficient, rl Ro s=m2l U equivalent heat capacity Subscripts eff effective h,H pertaining to hand H, respectively in inlet k kinetic l liquid i, e, h inlet, exit and heated, respectively s solid sat saturation x pertaining to axial directions v vapour w wall Superscripts * dimensionless

to accommodate the more general Local Thermal Non-Equilibrium (LTNE) condition [11,21,22].3 An additional problem that is encountered while simulating the complete liquidevapour phase change process inside porous media is the occurrence of rapid, non-physical change in the predicted properties (e.g., temperature) over an extremely short distance4 that results primarily due to the presence of discontinuities in the effective diffusion coefficient across the interfaces between the single and the two phase regions. Recognising this difficulty, different researchers addressed this issue either by modifying the manner in which the cell face diffusivities are interpolated from their adjacent nodal values [18] or by assuming the saturation enthalpies (hl,sat, hv,sat) and hence the latent heat (hfg) to be functions of the saturation pressure [23e26]. Although

3 For LTNE condition, the solid matrix and the fluid medium can locally coexist at different temperatures and the respective energy conservation equations are coupled to each other through the volumetric heat exchange term. 4 Henceforth shall be referred to as the “jump”, which has no relation to the physical mass and heat “jump” conditions across the liquidevapour interfaces, used in Volume of Fluid (VOF)-like methods.

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

neither of these approaches could eliminate the aforementioned discontinuities, the latter approach allows explicit, temperature gradient-driven, heat diffusion even in the two-phase region.5 Such a remedy, however, does not perform equally well for cases where the pressure drop inside the porous medium is negligible owing to the extremely low mass flow rate6 and hence the saturation temperature remains nearly the same in the entire two-phase region. Nevertheless, after critically analysing the reasons for the occurrence of “jumps” in the predicted properties, Alomar et al. [22,27] suggested remedies using the smoothing of effective diffusion coefficient for both LTE [27] and LTNE [22] models and applied it successfully for the simulation of complete phase change process inside divergent porous evaporators [28]. In view of this brief literature review, it is now apparent that the most popular TPMM [12] uses the modified volumetric enthalpy as the dependent variable in the energy conservation equation, which may be regarded as a synthetically defined variable, rather than a physical quantity. It would, therefore, be worthwhile to identify the problems associated with the original enthalpy formulation [9] and propose an alternative that would perform equally well or even better than the existing H-formulation [12]. In the rest of this article, the newly proposed modification shall be referred to as the “modified h-formulation”. In order to demonstrate its usefulness, the complete phase change process of water inside a porous evaporator of known properties has been considered for two different heating conditions7 at the evaporator wall. The crosssection of the duct has been assumed to remain invariant in the axial (flow) direction and for simplicity, the flow has been considered to be one-dimensional. Nevertheless, the formulation could be easily extended for multi-dimensional problems and even involving complex geometries that require the employment of curvilinear coordinates. In order to eliminate the occurrence of non-physical “jumps” in the predicted properties, however, the smoothing of effective diffusion coefficient has been used by suitably adopting the recommendations of Alomar et al. [22,27] for Hformulation. The results obtained with the modified h-formulation have been compared with that predicted by the existing Hformulation [12,22,27] for different parametric variations. The present article is organised in six sections. The original enthalpy formulation and its proposed modification are discussed in section 2, while section 3 deals with the reduced, onedimensional formulation of the considered problem in addition to the physical description and the boundary conditions. The numerical method is fairly straightforward and hence is briefly discussed in section 4 for the sake of completeness along with the smoothing algorithm adopted for the present modification. The necessity for smoothing of effective diffusion coefficient and the results predicted by the modified h-formulation are presented and compared with that obtained with the H-formulation [12,22,27] in section 5. Finally, the conclusions and the final remarks are reported in section 6.

the conventional SFM as [9]:

vr þ V$ðruÞ ¼ 0 vt

(1)

K ru ¼  ðVp  rk bÞ n

(2)

ε

where r ¼ rls þ rv(1  s) is the mixture density and rk ¼ rlll þ rvlv is the kinetic density with ll and lv ¼ 1  ll being the relative mobilities of the liquid and the vapour phases, respectively [9]. Since the non-linear Forchheimer term as well as the contributions of inertia and viscous forces are neglected, Eq. (2) is applicable only for low mass flow rate applications. For moderate temperature change and hence density variation in the single phase regions, the Boussinesq approximation may be invoked and hence rk in the bodyeforce term of Eq. (2) may be expressed as a linear function of the temperature difference [10]:

rk ¼ rl ½1  bl ðT  Tsat Þll þ rv ½1  bv ðT  Tsat Þlv

2.2. Critical appreciation of the existing enthalpy formulation The energy conservation equation under the assumption of LTE may be written as [9]:

vT v ð1  εÞrs Cps þ ε ðrhÞ þ V$ðgh ruhÞ vt vt     dT ¼ V$ keff Vh þ V$ jhfg dh

(4)

where rh ¼ rlhls þ rvhv(1  s) is the volumetric enthalpy of the mixture. Although the form of Eq. (4) is quite similar to that finally adopted by Wang and Beckermann [9], it contains two distinct transient terms, where the first and the second terms on the left hand side signify the rate of change of stored energy within the solid medium and in the fluid volume, respectively. These two terms may be combined together and could be expressed as:

  vT v v vr þ ε ðrhÞ ¼ Uh ðrhÞ  ðUh  εÞ h vt vt vt vt

2. Original enthalpy formulation and proposed modification 2.1. Mass and momentum conservation equations

where Uh is defined as:

The governing mass and momentum conservation equations that remain identical for all variants of TPMM may be derived from

Uh ¼ ε þ

Since the saturation temperature is assumed to depend on the local pressure. 6 The mass flow rate could be as low as 100 g/hr. 7 Constant temperature [27] and prescribed heat flux [22] conditions are referred here.

(3)

where rl and rv are the densities of the liquid and the vapour phases, respectively, at T ¼ Tsat, whereas bl and bv are the isobaric expansion coefficients for the liquid and the vapour phases, respectively, at the saturation conditions. For problems, associated with large temperature differences in the single phase regions, however, more accurate models should be employed for the density variations. The permeability of porous medium K can be either evaluated from a separate detailed study [29] or may be estimated using empirical correlations [30]. In the rest of the article, however, owing to the simplicity, the permeability of the porous medium is assumed to be a known quantity, by disregarding the manner in which it is obtained or prescribed.

ð1  εÞrs Cps

5

125

  ð1  εÞrs Cps dT dh r

(5)

(6)

At this point, Wang and Beckermann [9] dropped the second term on the right hand side of Eq. (5), which is justifiable only if the densities in the single phase regions are assumed to be independent of the local temperature and hence vr/vt ¼ 0, although vT/vt may be finite. In the two-phase region, however, this term is exactly

126

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

zero since the phase change process takes place without any appreciable change in pressure, which results in dT/dh ¼ 0 and hence from Eq. (6), one obtains Uh ¼ ε. In Eq. (4), gh ¼ rhk/rh is the advection correction coefficient [9], where hk ¼ hlll þ hvlv is the kinetic enthalpy. In the single phase regions gh ¼ 1, whereas for the two-phase mixture, gh may be expressed as:

gh ¼

  ½rl s þ rv ð1  sÞ hv;sat  ll hfg rl hl;sat s þ rv hv;sat ð1  sÞ

(7)

Further in Eq. (4), j represents the total mass flux (velocity) from one phase to the other, which is a combination of the capillaryeinduced diffusive flux and the gravityeinduced migrating flux [9]:

j ¼ DVs þ f

KDr b nv

(8)

where D is the capillary diffusion coefficient and f is the hindrance coefficient [9]:

D¼

  ðεKÞ1=2 s dJ ll lv  n ds

(9a)

f ¼ krv ll

(9b)

It may be easily verified that in the single phase regions, both D and f are exactly zero. Nevertheless, if required, particularly for multi-dimensional problems and for LTNE model, Eq. (8) may be employed in order to retrieve the massevelocities for the liquid and the vapour phases as [9,22]:

rl ul ¼ ll ru þ j ¼ ð1  lv Þru þ j

(10a)

rv uv ¼ lv ru  j ¼ ð1  ll Þru  j

(10b)

Substituting j from Eq. (8), the last term on the right hand side of Eq. (4) is obtained as:

      KDrhfg V$ jhfg ¼ V$ Dhfg Vs þ V$ f b nv

(11)

In order to express Vs in terms of Vh, Wang and Beckermann [9] obtained a relationship for s as a function of r and h as follows:

s¼

  r hv;sat  h rl hfg

(12)

This expression may be derived by substituting rv from the definition of r and using the relationship for rh while replacing hl and hv by their respective saturation values. Nevertheless, from Eq. (12), Vs may be expressed as:

r h  hv;sat Vs ¼  Vh þ Vr rl hfg rl hfg

! (13)

Substituting Vs in Eq. (11) and using it further, the final form of the energy conservation equation (4) may be obtained as follows:

vT v þ ε ðrhÞ þ V$ðgh ruhÞ vt vt 

  KDrhfg D ¼ V$ðGho VhÞ þ V$ b h  hv;sat Vr þ V$ f rl nv

ð1  εÞrs Cps

(14)

where Gho is the effective diffusion coefficient for h and is defined

as [9]:

Gho ¼

rD dT þ keff rl dh

(15)

It is obvious that the first term in Eq. (15) is absent for the single phase regions as D ¼ 0, whereas the second term does not play any role in the two-phase region since dT/dh ¼ 0. It is worth noting that other than replacing the first two terms on the left hand side of Eq. (14) by the first term on the right hand side of Eq. (5), Wang and Beckermann [9] also represented the second term on the right hand side of Eq. (14) as:

V$

    D h  hv;sat Vr ¼ V$ Gho h  hv;sat Vðln rÞ rl

(16)

Such a substitution, however, is valid only if the densities for the liquid and the vapour phases are considered to be constants, i.e., only when rl ¼ rl,sat and rv ¼ rv,sat are assumed. Therefore, Eq. (14) may be considered as the final form of the energy conservation equation, without imposing any restrictive assumption with respect to the densities in the single phase regions. In view of the original enthalpy formulation, the following comments are now in order: 1. The relationship between s and h in Eq. (12) is implicit since it contains r on the right hand side, which is an explicit function of s. During iteration, h is first obtained by solving Eq. (14) using semi-implicit treatment. Using the new h in the two-phase region, s is first evaluated from Eq. (12), since all other mixture properties explicitly depend on it. Therefore, the process would require additional computation time. 2. The second term in Eq. (14) on the right hand side represents the diffusion of r in the two-phase region, which has to be accommodated as a source in the discretised equation. Since its evaluation would lag the current calculation by one iteration, additional computation time would be required in order to achieve convergence. 3. If Eq. (14) is solved for the volumetric enthalpy rh (instead of h) as the independent variable, additional problems may arise since rh is not a monotonically increasing function of the thermodynamic state, particularly in the two-phase region [12]. This apprehension, however, would be meaningless if the mixture enthalpy h is solved from Eq. (14), as has been adopted for the modified h-formulation. 4. In the two-phase region, Gho ¼ 0 is obtained from Eq. (15) for s ¼ 1 and s ¼ 0 since both D and dT/dh are exactly equal to zero. Therefore, the formulation is liable to produce jumps in the predicted properties at the phaseechange interfaces, unless special care is taken. This problem, however, persists for all variants of TPMM [22,27,28] and may be avoided using the smoothing of effective diffusion coefficient [27]. It would, therefore, be worthwhile to explore an alternative formulation that eliminates the aforementioned drawbacks, while considering h as the dependent variable of the energy conservation equation. This formulation requires certain modifications, that are described next in this article. It is expected that the modified hformulation would still retain its benefits even for LTNE condition since the governing equations for the fluid phase remain nearly the same, except for the presence of an additional volumetric heat exchange term [22]. 2.3. Proposed modification An explicit relation between s and h can be easily obtained by

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

substituting r as function of s in Eq. (12) as:

  rv hv;sat  h    s¼  rl h  hl;sat þ rv hv;sat  h

(17)

Therefore, the first problem associated with the original enthalpy formulation is eliminated. Since r ¼ rv þ (Dr)s is also an alternative expression, one may substitute Vr ¼ (Dr)Vs in Eq. (13) and after some simple algebraic manipulations, Vs could be expressed as:

rVh    Vs ¼   rl h  hl;sat þ rv hv;sat  h

(18)

H þ rv hv;sat s¼ rl hfg þ Drhv;sat

(24a)

VH Vs ¼  rl hfg þ Drhv;sat

(24b)

Substituting Vs in Eq. (11) and following straightforward mathematical manipulation, the applicable energy conservation equation (4) is represented as:

UH

  KDrhfg vH þ V$ðgH uHÞ ¼ V$ðGH VHÞ þ V$ f b vt nv

This relation could also be derived by differentiating Eq. (17) in space while using the following relationship8:

where UH, gH and GH are defined as:

     rl rv hfg ¼ r rl h  hl;sat þ rv hv;sat  h

UH ¼ ε þ ð1  εÞrs Cps

(19)



Substituting Vs from Eq. (18), the alternative final form of Eq. (4) may be obtained using Eq. (11) as follows:

vT v þ ε ðrhÞ þ V$ðgh ruhÞ vt vt   KDrhfg ¼ V$ðGh VhÞ þ V$ f b nv

gH ¼

ð1  εÞrs Cps

(20)

where Gh is the modified effective diffusion coefficient for enthalpy transport, defined as:

Gh ¼

rDhfg dT  þ keff   dh rl h  hl;sat þ rv hv;sat  h 

GH ¼

The local temperature may be retrieved from h using Eq. (22) as:

T ¼ Tsat þ

h  hl;sat Cpl

for liquid phase

(23a)

T ¼ Tsat þ

h  hv;sat Cpv

for vapour phase

(23b)

Therefore, dT/dh ¼ 1/Cpl and 1/Cpv are obtained for the liquid and the vapour phases, respectively. For the two-phase region, however, dT/dh ¼ 0, owing to the thermodynamic constraint, as mentioned before.

2.4. Modified volumetric enthalpy formulation of Wang [12] According to the widely used H-formulation [12], Eq. (4) is rewritten with the synthetically defined H ¼ r(h  2hv,sat) as the dependent variable. Therefore, s and Vs are obtained as:

8 Since Eqs. (12) and (17) are two alternative expressions for s, their right hand sides can be equated.

(26c)

According to this formulation, explicit relations between T and H may be obtained as [12]:

¼ Tsat þ

(22b)

(26a)

(26b)

Dhfg dT þ keff dH rl hfg þ Drhv;sat

Since Eq. (20) contains only the standard diffusion term on the right hand side, the second problem, associated with the existing hformulation, is also successfully removed. Under the assumption of constant specific heats, enthalpies of the liquid and the vapour phases could be written as functions of temperature as:

hv ¼ hv;sat þ Cpv ðT  Tsat Þ

(25)



rhv;sat þ rl hfg s

T ¼ Tsat þ

(22a)

dT dH

  r hv;sat þ ll hfg

(21)

hl ¼ hl;sat  Cpl ðTsat  TÞ ¼ hl;sat þ Cpl ðT  Tsat Þ

127

T ¼ Tsat þ

  H þ rl 2hv;sat  hl;sat rl Cpl H  Hl;sat rl Cpl

for liquid phase

H þ rv hv;sat H  Hv;sat ¼ Tsat þ rv Cpv rv Cpv

(27a)

for vapour phase (27b)

where the second parts of Eq. (27) apply only when rl and rv are assumed to be constants. In this situation, dT/dH ¼ 1/rlCpl and 1/ rvCpv are obtained for the liquid and the vapour phases, respectively, and dT/dH ¼ 0 still holds for the two-phase region, owing to the assumed thermodynamic constraint. Although the expressions for T for both h- and H-formulations in Eqs. (23) and (27), respectively, are remarkably similar to each other, it will be shortly apparent that such similarity exists only when the densities in the single phase regions are assumed to be constants.

2.5. Issues related to variable properties For extremely low mass flow rate applications, the static pressure of the fluid inside the evaporator may be taken equal to the atmospheric pressure. Therefore, the variations in r and h for water may now be examined at p ¼ 101.325 kPa corresponding to Tsat ¼ 100  C [31], as presented in Figs. 1 and 2, respectively, for the single phase regions. It is evident from these figures that although h is a linear function of T in the single phase regions, which justifies the use of Eq. (22) along with the assumption of constant Cp for all variants of TPMM, the variations in r in both liquid and vapour phases are non-linear functions of T. As a result, H also turns out to be a non-linear function of temperature. The raw data for rl and rv may be expressed in the form of quadratic functions of temperature:

128

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Fig. 1. Density variations in the single phase regions (a) liquid phase, (b) vapour phase.

Fig. 2. Enthalpy variations in the single phase regions (a) liquid phase, (b) vapour phase.

h i rl ¼ rl;sat 1  bl ðT  Tsat Þ  al2 ðT  Tsat Þ2

(28a)

h i rv ¼ rv;sat 1  bv ðT  Tsat Þ þ av2 ðT  Tsat Þ2

(28b)

where bl and bv are the isobaric expansion coefficients at saturation conditions for the liquid and the vapour phases, respectively. They are obtained by fitting the raw data as9:

1 bl ¼  rl;sat 1

bv ¼  rv;sat



vr vT





vr vT

p;sat



1 ¼ rl;sat 1

p;sat



¼ rv;sat

drl dT





drv dT

¼ 8  104

(29a)

sat



significantly affect the predicted solution. Nevertheless, if property variations in the single phase regions are taken into account, rl and rv, required for different mixture properties in the two-phase region, are replaced by rl,sat and rv,sat, respectively, although the rest of the h-formulation still remains the same. For the H-formulation, however, the single phase temperatures are given only by the first part of Eq. (27). Since rl and rv are approximated by the quadratic functions of temperature, determination of T from the known values of H requires solution of cubic equations resulting from Eqs. (27) and (28) and hence demands additional computation time. The expressions for dT/dH, required in Eq. (26), are obtained as:

 3

¼ 2:47  10

(29b)

sat

The second coefficients in Eq. (28) are obtained as al2 ¼ 3:4  106 and av2 ¼ 3:25  106 . Fig. 1 also shows the variations in rl and rv according to the Boussinesq approximation. It is evident that since the density variation for the sub-cooled liquid phase is not substantial, the Boussinesq approximation or even the assumption of constant density for the liquid phase may not

9 The thermal expansion coefficients, particularly bl, are different from that used by other researchers [12e14], although these values should have been used at the saturation condition.

dT dH



dT dH

 ¼

rl r2l Cpl þ Hðdrl =dTÞ

(30a)

¼

rv r2v Cpv þ Hðdrv =dTÞ

(30b)

l

 v

where drl/dT and drv/dT are given according to Eq. (28) as:

drl ¼ rl;sat ½bl þ 2al2 ðT  Tsat Þ dT

(31a)

drv ¼ rv;sat ½bv  2av2 ðT  Tsat Þ dT

(31b)

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

129

It is now evident that although for the modified h-formulation, the expression for Gh still remains unaltered, GH for the H-formulation turns out to be a strong function of the local properties in the single phase regions. It is, therefore, obvious that the H-formulation of Wang [12] requires appropriate adaptation in order to accommodate substantial variations in r in the single phase regions. 2.6. General remarks on variants of TPMM Since Eqs. (1) and (2) remain invariant for all versions of TPMM and the energy conservation equations (20) and (25) for the h- and the H-formulations, respectively, are remarkably similar to each other, they are solved in a similar manner as elaborately described by Alomar et al. [27]. However, if r in the single phase regions are functions of T, Eq. (27) turns out to be a cubic equation. Although the roots of a cubic equation may be explicitly obtained [32e37], the expressions for dT/dH in Eq. (30) contain additional non-linear terms. Therefore, for variable properties, the H-formulation is more complex than it was originally presented for constant densities in the single phase regions. No such additional requirement, however, could be detected for the modified h-formulation and hence it should be preferred over the established H-formulation. For all variants of TPMM, keff of the composite medium, that consists of a solid matrix and the working fluid, is required. Although keff could be determined from a separate detailed analysis [38e41] of the porous structure, in the present investigation, it is obtained from the simple parallel arrangement model as:

keff ¼ ð1  εÞks þ ε½kl s þ kv ð1  sÞ

(32)

Although it has been theoretically established that at least for substantial density variations in the single phase regions, the modified h-formulation should perform better than the existing Hformulation, it would be worthwhile to demonstrate that these formulations produce identical results for a given problem with constant rl and rv, owing to their consistency. Therefore, the complete phase change problem of water inside a constant crosssectional area porous evaporator for the reduced onedimensional model has been taken up next as a demonstrative example. Although the pressure drop through the evaporator has been considered to be negligible in the present analysis and hence the phase change process has been assumed to occur at constant saturation temperature of water corresponding to the atmospheric pressure p ¼ 101.325 kPa, for multi-component evaporation or for flows with substantial pressure drop, the phase change temperature is expected to vary. These issues, along with that for the variable specific heat, may also be accommodated quite easily in the present formulation by including appropriate constitutive relations, although they are left beyond the scope of the present study. 3. Formulation of the reduced one-dimensional problem

Fig. 3. Geometry of the porous evaporator.

evaporator due to the application of an external pressure gradient, which is then heated to the superheated vapour state. 3.2. Dimensionless governing equations and boundary conditions The equations governing the respective one-dimensional formulation have been obtained by integrating the multidimensional (axi-symmetric) conservation equations (1), (2), (20) and (25) over the evaporator cross-section10. For the H-formulation [12], the procedure has been elaborately presented by Alomar et al. [27]. Since the integration process is similar also for the h-formulation, it is not repeated here for the sake of brevity. All applicable equations for constant rl and rv have been non-dimensionlised using the dimensionless variables listed in Table 1. They retain all governing equations, constitutive relations and the expressions for mixture properties (refer to Table 2 for their definitions) in the form that remain identical to their dimensional counterparts. For variable density formulation, however, rl in Table 1 has been replaced by rl,sat, whereas rl and rv appearing in Table 2 have been replaced by their saturation values. For simplicity, only the steadyestate phase change process has been considered for both h- and H-formulations, with LTE condition. Integrating Eqs. (1) and (2), over the evaporator crosssection, one obtains [27]:

d  * * r ux ¼ 0 dx*

(33)

  K * dp* * *  r b r* u*x ¼  * k x n dx*

(34)

For constant wall heat flux applications, integration of Eqs. (20) and (25) yields the following one-dimensional forms of the energy conservation equations [22,27]:

 *  d  d d * * * * vh G r u h g ¼ x h vx* þ dx* dx* h dx*

f

K * Dr* h*fg n*v

! b*x

þ

P * 00 * q_ A*c w (35)

 *  d  d d * * * vH G þ * u H g ¼ H dx* H x dx* dx vx*

f

K * Dr* h*fg n*v

! b*x

þ

P * 00 * q_ A*c w (36)

3.1. Description of the physical problem The schematic representation of the problem geometry is shown in Fig. 3. It has unheated starting and exit lengths of li and le, respectively and the length of the heated section is lh. The evaporator wall in the heated section has been either subjected to a 00 prescribed q_ w or maintained at a constant temperature Tw>Tsat. In the latter case, however, heat is assumed to be transferred from the evaporator wall to the combined medium due to a specified heat transfer coefficient hco. The sub-cooled liquid water at the inlet with a velocity uin and at a temperature Tin < Tsat flows through the

00 *

00

* =Re with Q * ¼ q_ R=mh is the prescribed dimenwhere q_ w ¼ Qw l fg w w sionless heat flux at the duct surface. During integration, the heat flux condition at the evaporator wall has been incorporated into an appropriate source term. Since all lengths have been made dimensionless with respect to the radius of the evaporator, P* ¼ 1 and A*c ¼ 1=2 are obtained per unit radian and hence P * =A*c ¼ 2 may be substituted.

10

Variations remain only in the axial direction.

130

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Table 1 Definitions of dimensionless variables [27].

Table 2 Definitions of (dimensionless) properties for two-phase mixture [9,12,27].

Dimensionless variables

Expressions

Variables

Length Density Velocity Pressure

x* ¼ x/R; L* ¼ L/R; l* ¼ l/R r* ¼ r=rl ; r*l ¼ 1 u* ¼ u=uin

Density

Temperature

T * ¼ TCpl =hfg h* ¼ h=hfg ; h*fg ¼ 1 H* ¼ H=rl hfg ¼ r* ðh*  2h*v;sat Þ

Enthalpy Kinetic density

p* ¼ p=rl u2in

Enthalpy Modified volumetric enthalpy Specific heat Dynamic viscosity Thermal conductivity Body force per unit mass

Cp* m*

* ¼1 Cp =Cpl ; Cpl m=ml Rel ; m*l ¼ 1=Rel

Isobaric expansion coefficient

b ¼ bhfg =Cpl

Wall heat flux

* =Re ; Q * ¼ q_ R=mh q_ w ¼ q_ w =rl uin hfg ¼ Qw l fg w w K* ¼ K/R2

¼ ¼ k* ¼ k=rl Cpl Ruin ¼ k=kl Pel ; k*l ¼ 1=Pel ~x gR=u2 ¼ b ~x =Fr2 b* ¼ b x *

00 *

Permeability

in

00

00

For constant wall temperature, the wall heat flux may be 00 modelled as q_ w ¼ hco ðTw  TÞ [27] and hence in the dimensionless form, Eqs. (35) and (36) may be obtained as:

   d  d dh* d gh r* u*x h* ¼ * G*h * þ * * dx dx dx dx  Nuo  * Tw  T * þ Pel  *  d  d d * * * dH G þ * u H g ¼ H x H dx* dx* dx* dx  Nuo  * Tw  T * þ Pel

f

K * Dr* h*fg n*v

! b*x

r* ¼ r*l s þ r*v ð1  sÞ

r* h* ¼ r*l h*l s þ r*v h*v ð1  sÞ

* Þl þ r* ½1  b* ðT *  T * Þl r*k ¼ r*l ½1  b*l ðT *  Tin v l v v sat

Kinetic enthalpy

h*k ¼ ll h*l þ ð1  ll Þh*v

Kinematic viscosity

n* ¼ ½krl =n*l þ krv =n*v 1

Effective thermal conductivity Advection correction coefficient Advection correction coefficient Effective diffusion coefficient Effective diffusion coefficient Capillary diffusion coefficient Surface tension coefficient Relative mobilities

k*eff ¼ ð1  εÞk*s þ ε½k*l s þ k*v ð1  sÞ h*

gh ¼ hk* ¼

r* ½h*l ll þh*v ð1ll Þ r*l h*l sþr*v h*v ð1sÞ

h* 2h*

gH ¼ hk* 2h*v;sat ¼ v;sat

G*h ¼ r* ðh* h* l

r* ðh*v;sat þll h*fg Þ r* h*v;sat þr*l h*fg s

r* D* h*fg þ Þþr*v ðh*v;sat h* Þ

l;sat

D* h*

*

k*eff dT dh*

*

G*H ¼ r* h* þDrfg* h*

dT þ k*eff dH *   ðεK * Þ1=2 s* dJ * D ¼ ll ð1  ll Þ ds n* l

v;sat

fg

~=Re2l s* ¼ s=rl u2in R ¼ s ll ¼ k

krl =n*l * * rl =nl þkrv =nv

krv =n*v * * rl =nl þkrv =nv

; lv ¼ k

¼ 1  ll

krl krv =n*l krl =n*l þkrv =n*v

Hindrance coefficient

f ¼ krv ll ¼

Relative permeabilities Capillary pressure function

krl ¼ sn ; krv ¼ ð1  sÞn ; n ¼ 1 J ¼ 1.417(1  s)2.120(1  s)2 þ 1.263(1  s)3

(37)

f

K * Dr* h*fg n*v

! b*x (38)

where Nuo ¼ hco(2R)/kl is the characteristic Nusselt number. In unheated sections of the evaporator, i.e., for 0  x*  l*i and l*i þ l*h  x*  L* , however, the third term on the right hand sides of Eqs. (35)e(38) has been set to zero. The relations for evaluating T* and s from the solution of either h or H are already summarised in section 2 in their dimensional form and hence they are not presented here for the sake of brevity. In order to solve Eqs. (33)e(38), boundary conditions for both hand H-formulations are required only in the axial direction at the inlet and the exit. At x* ¼ 0, u*x ¼ 1 has been set. From the known * < T * , h* for the h-formulation and H * for the H-formulation Tin sat in in have been specified according to Eqs. (22a) and (27a), respectively, * from Eq. (28a). At x* ¼ L*, the where r*in has been obtained at Tin second derivatives of all variables are set to zero11 and hence the outlet h* and H* have been obtained by linear extrapolation [27]. Integrating Eq. (33), the local mass flux may be obtained as:

r* u*x ¼ r*in u*in

(39)

Therefore, if required, ¼ may be used, where r* is the local density. According to the modified h-formulation, since r* u*x explicitly appears only in the convective terms of Eqs. (35) and (37) and remains constant in the axial direction, irrespective of whether r*l and r*v are constants or strong functions of T*, identical solutions are obtained. Such a simple conclusion, however, cannot be drawn for the H-formulation since i) instead of the mass flux, only u*x appears u*x

Expressions

r*in u*in =r*

11 This condition is less stringent as compared to the fully-developed condition for shorter exit lengths.

in the convective terms of Eqs. (36) and (38), which strongly depends on r*, ii) evaluation of T* from H* in the single phase regions according to Eq. (27) explicitly depends on r* and iii) variable r*l and r*v also significantly modifies G*H in the single phase regions (see Eqs. (26), (30) and (31)). Therefore, the only way to reach an identical conclusion from the H-formulation would be to simultaneously simulate variable as well as constant density (for the single phase regions) cases, while considering identical problems and obtain identical temperature12 profiles from both solutions. This could be regarded as one of the drawbacks of the H-formulation [12]. For multi-dimensional problems, however, the assumption of variable r*l and r*v may lead to different solutions as compared to that obtained by considering them to be constant. Similarly, for LTNE condition, since the volumetric heat exchange coefficient between the solid and the fluid phases depends strongly on the magnitude of local velocity, and hence the local density as already demonstrated in Eq. (39), the solutions obtained for constant and variable r*l and r*v are expected to differ from each other, even though the problem may still be treated as one-dimensional. Nevertheless, since the major objective of the present article is to propose a modified h-formulation based on LTE assumption that performs equally well or even better than the existing H-formulation, which requires the definition of an artificial modified volumetric enthalpy, these issues have been left beyond the scope of the present investigation and should be taken up later in subsequent articles. Nevertheless, since the variations in rl and rv do not affect the temperature solution for steady one-dimensional problems, they have been assumed to be constants (equal to their saturation values) for all demonstrative examples, presented in this article. Therefore, r*in u*in has been set to unity and hence the local velocity, which is required for the H-formulation, has been obtained as

12 These solutions would not produce identical H*, but would yield identical h* and hence T*.

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

u*x ¼ 1=r* . It is also evident that since for the present steadyestate, oneedimensional problem, the velocity solution does not depend on the pressure field, Eq. (34) may be dropped unless the solution for the pressure distribution is required [27]. Therefore, the numerical solutions are required only for Eqs. (35)e(38), which is discussed briefly in the next section. 4. Numerical solution In the present investigation, the Finite Volume Method (FVM) along with the staggered grid layout has been used for the discretisation of Eqs. (35)e(38), where only the velocities are stored at the cell faces and all other variables are defined at the cell centres. For the H-formulation, the densities are required at the cell faces, which have been obtained by linear interpolation, although there is no such requirement for the h-formulation since r* u*x remains constant in the axial direction. Further, the comparison of Eqs. (35) and (36), as well as Eqs. (37) and (38) clearly shows that the sets of energy conservation equations for h- and H-formulations, with identical heating condition at the wall, remain nearly the same, except for the expressions of g and G*. The discretisation procedure for the Hformulation has been discussed in details by Alomar et al. [27] and since the procedure for the modified h-formulation is also similar, only a brief description is presented here for the sake of brevity. The first order accurate Upwind Differencing Scheme (UDS) and the second order accurate Central Differencing Scheme (CDS) have been employed in order to express the convective and the diffusive terms, respectively [42]. In addition, G*h and G*H at the cell faces have been obtained using the harmonic mean approximation, which ensures the balance of diffusive energy fluxes at the cell faces [42]. The advection correction coefficients gh and gH, given in Eqs. (7) and (26b), as well as the hindrance function f at the cell faces have been obtained by the linear interpolation. The discretised equations along with the boundary conditions have been solved employing the Thomas algorithm [42] for tridiagonal matrix. However, owing to the presence of strong nonlinearity in Eqs. (35)e(38), the coefficients in the discretised equations have been obtained in a semi-implicit manner and hence the discretised equations have been solved iteratively. During iterations, however, the discretised equations have been further under-relaxed and typical under-relaxation factors of 0.4 and 0.1 have been used for the h- and the H-formulations, respectively, for most of the simulations reported in this article. For all tested cases, the convergence criterion has been set to 106, which has been judiciously chosen such that the choice does not significantly affect the final solution. The special treatment for G*H , which is required for avoiding jumps in the predicted properties, has been elaborately explained by Alomar et al. [22,27,28] for the H-formulations resulting from both LTE and LTNE assumptions. In the present investigation, a similar strategy has been adopted also for the modified h-formulation and hence it is not repeated here for the sake of brevity. The application of smoothing algorithm restricts both G*h and G*H to finite values by preventing them from going to zero, at s ¼ 0 and s ¼ 1 in the two-phase region. 5. Results and discussion 5.1. Ranges of parameters and general remarks In order to demonstrate the applicability of the proposed enthalpy formulation, water has been considered as the working fluid and its properties are listed in Table 3, where Tsat has been taken as 100+C [27,31,43] at the atmospheric pressure of 101.325 kPa. For all tested cases reported in this article, a circular

131

pipe of radius R ¼ 25 mm (R* ¼ 1) and with an overall length L ¼ 500 mm (L* ¼ 20) has been considered, where the length of the heated section depends on li and le. The inlet length has been varied between 0 and 0.2L, whereas the exit length has been kept fixed at 0.1L. For the reference case, however, both li and le have been taken as 0.1L and hence lh ¼ 0.8L. The properties of porous media, i.e., ε, K and ks for the compressed wire-mesh structures prepared from low-carbon steel [43] have been varied from 0.2 to 0.4, 1013 m2 to 1010 m2 and 10 W/mK to 30 W/mK, respectively. Sub-cooled liquid water at Tin ¼ 20 oC enters the evaporator, with an extremely low mass flow rate that varies between 100 g/hr and 500 g/hr. Therefore, uin has been varied from 0.015 mm/s to 0.075 mm/s and hence Rel ¼ uinR/nl varies between 1.25 and 4. At 00 the heated section of the evaporator, q_ w for the constant wall heat flux applications has been varied from 2 kW/m2 to 3 kW/m2, whereas for the constant wall temperature cases, Tw has been varied between 140 oC and 200 oC. Further for the latter cases, Nuo ¼ 4 has been assumed, which is required for solving Eqs. (37) and (38) [27]. In addition, since the body force vector per unit mass acts in the ~x ¼ 1 has been specified. The gravitational negative x*-direction, b Reynolds number Reg, required for the calculation of Fr and the ~ have been prescribed as normalised surface tension coefficient s 4.25104 and 1.81107, respectively, based on g ¼ 9.81 m/s2, s ¼ 0.0589 N/m and the properties of liquid water listed in Table 3. Therefore, in the present investigation, the following ranges of dimensionless parameters have been considered:

Starting and exit lengths : l*i ¼ 0  4ðvariableÞ; l*e ¼ 2ðfixedÞ Overall length : L* ¼ 20ðfixedÞ Inlet Reynolds number : Rel ¼ 1:25  4ðvariableÞ Porosity : ε ¼ 0:2  0:4ðvariableÞ Darcy number : K * ¼ 109  107 ðvariableÞ Thermal conductivity of the solid phase : k*s;ref ¼ 4:25  12:5 ðvariableÞ * ¼ 8:5  102  10:5  102 ðvariableÞ Heat Flux : Qw

Other than the fixed values, the reference case has been chosen with the following parameters: Rel ¼ 2 (the inlet condition, other than Tin ¼ 20  C), ε ¼ 0.3, K* ¼ 108, k*s ¼ 8:25 (properties of the * ¼ 9:5  102 (for constant heat flux), porous medium), Qw Tw ¼ 160 oC (for constant wall temperature, other than Nuo ¼ 4), L* ¼ 20, l*i ¼ l*e ¼ 2 and hence l*h ¼ 16 (the problem geometry). These parameters have been kept unaltered for all cases, unless

Table 3 Thermo-physical properties of water [27,31 and 43]. Property

Liquid

Density, r (kg/m3) Specific heat, Cp (J/kgK) Dynamic viscosity, m (kg/ms) Thermal conductivity, k (W/mK) Saturation enthalpy, h (J/kg) Prandtl number Pr ¼ mCp/k Latent heat of evaporation, hfg (J/kg) Surface tension coefficient, s (N/m)

957.85 0.5978 4190.2 2029.0 4 2.79  10 1.202 105 0.68 0.0248 419.02  103 2676.05  103 1.72 0.983 2257.03  103 0.0589

Vapour

132

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Fig. 4. Effects of heating and flow conditions on the axial temperature distribution for constant wall heat flux cases with ε ¼ 0.3, K* ¼ 108, k*s ¼ 8:25 and l*i ¼ 2: (a) effect of applied *. heat flux for Rel ¼ 2 and (b) effect of inlet Reynolds number for different Qw

* ¼ 9:5  102 and l* ¼ 2: (a) effect of porosity (b) Fig. 5. Effects of porous media properties on the axial temperature distribution for constant wall heat flux cases with Rel ¼ 2, Qw i effect of Darcy number and (c) effect of thermal conductivity of the solid phase.

their effects on the phase change process have been analysed. Prior to obtaining the results, a detailed grid independence study has been carried out and it has been observed that at least 400 uniform control volumes are required in order to obtain “nearly” gridindependent results and further grid refinement does not significantly modify the presented solutions. The reference thermal conductivity of the solid phase k*s;ref requires a special mention. The definition of k*s (or, for that matter, k*) in Table 1 depends on Rel. Therefore, its numerical value would

differ for different Rel even if the thermal conductivities of the solid and the fluid phases are kept unaltered. In order to compare the data obtained for different Rel for identical properties of the solid and the fluid phases, k*s;ref has been introduced, which corresponds to k*s for Rel,ref ¼ 2. Therefore, for a prescribed Rel which is different from Rel,ref, the true value of k*s can be easily obtained as k*s ¼ k*s;ref Rel;ref =Rel [22]. For both modified h- and existing H-formulations, it has been observed that unless the smoothing algorithm [27], is

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Fig. 6. Effect of inlet length on the axial temperature distribution for constant wall * ¼ 9:5  102 , ε ¼ 0.3, K* ¼ 108 and.k* ¼ 8:25. heat flux cases with Rel ¼ 2, Qw s

implemented, the solutions always produced “jumps” in the predicted properties, particularly close to the saturated vapour condition. The results are, however, not presented in this article for the sake of brevity. Nevertheless, all results presented in this article have been finally obtained using a code that works “without” smoothing, while providing the predictions obtained “with” smoothing as their guessed solutions. Therefore, it can be safely concluded that although the proposed smoothing algorithm [27] for the modified h-formulation provides a successful remedy for the occurrence of “jumps” in the predicted properties, it does not modify the expected solutions. It may be mentioned here that all reported computations have been carried out on a desktop personal computer that has Intel Core i7 2.80 GHz CPU with 4 cores, each consisting of 4 GB DDR3 1333 MHz RAM made by Hynix semiconductor, as its configuration, although its parallel processing features have not been utilised during the present investigation. The typically required computation time on this machine for the present onedimensional phase change problems has been recorded to vary between 1 and 12 min. 5.2. Comparisons of results The axial temperature distributions T(x*) obtained from both h* in Figs. 4e6 that and H-formulations are compared for constant Qw

133

exhibit nearly identical predictions and excellent agreements. For the reference case, the average absolute difference in the predicted T* (or T) has been found to be 0.04%, which expectedly decreases with the increase in employed grid size. Particularly for this case, the modified h-formulation requires approximately 50% less computational time as compared to the existing H-formulation and hence is a better option for the simulation of the complete phase change process inside porous media. The effects of heating and flow conditions on T(x*) in Fig. 4 show * for a fixed Re enhances the outlet temthat the increase in Qw l perature Te, whereas with the increase in Rel, Te decreases considerably. Fig. 5 shows the effects of porous media properties on Te from which it is clear that although K* (Darcy number) and k*s have stronger effects on T(x*), ε has only marginal influence on Te. The effect of l*i on T(x*) in Fig. 6 is rather expected. Nevertheless, the axial diffusion in the upstream direction for all cases plays an important role on T(x*). Since these observations have been elaborately explained by Alomar et al. [22] for the existing H-formulation, they are not repeated here for the sake of brevity. In Figs. 7e9 the comparisons of T(x*), predicted by both h- and H-formulations, are presented for prescribed Tw cases, where Fig. 7 shows the effects of Tw and Rel, whereas the influences of porous media properties and l*i are shown in Figs. 8 and 9, respectively. The figures clearly demonstrate excellent agreements in the predicted T(x*) for all tested cases with constant Tw that ascertains the consistency of the modified h-formulation. The physical explanations for such variations have been already provided by Alomar et al. [17] and hence they are not repeated here for the sake of brevity. The average absolute difference in T* has been recorded to be 0.08% for the reference case, for which the present formulation requires approximately 12.5% less computational time as compared to the existing H-formulation. For the present one-dimensional problems, the heating conditions at the wall have been incorporated as source terms in Eqs. (35)e(38) that result directly from the integration process. Therefore, the constant Tw applications require iterative calculations of the source term and hence less computational time is saved * cases, where the source term can be compared to the constant Qw directly evaluated. For multi-dimensional problems, however, it is expected that problems with the Dirichlet boundary conditions (for constant Tw) in the transverse direction would require less computational time as compared to that with the Neumann * ) and hence the modified hboundary conditions (for constant Qw formulation would be even more beneficial. As mentioned in section 3, the issues related to the multi-dimensional problems shall

Fig. 7. Effects of heating and flow conditions on the axial temperature distribution for constant wall temperature cases with ε ¼ 0.3, K* ¼ 108, k*s ¼ 8:25 and l*i ¼ 2: (a) effect of wall temperature for Rel ¼ 2 and (b) effect of inlet Reynolds number for.Tw ¼ 160  C.

134

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

Fig. 8. Effects of porous media properties on the axial temperature distribution for constant wall temperature cases with Rel ¼ 2, Tw ¼ 160  C and l*i ¼ 2: (a) effect of porosity (b) effect of Darcy number and (c) effect of thermal conductivity of the solid phase.

formulation, the one-dimensional complete liquidevapour phase change problem of water inside a circular porous evaporator with known properties has been considered by employing two different heating conditions. The results obtained using the present modification and the existing H-formulation [12] have been compared and excellent agreements have been observed for all tested cases. The major conclusions from the present research effort may be summarised as follows:

Fig. 9. Effect of inlet length on the axial temperature distribution for constant wall temperature cases with Rel ¼ 2, Tw ¼ 160  C, ε ¼ 0.3, K* ¼ 108 and.k*s ¼ 8:25.

be taken up later in separate publications. 6. Conclusions and final remarks In the present investigation, the existing h- and H-formulations based on TPMM [9,12], have been critically analysed. After recognising their drawbacks, a modified h-formulation has been proposed that eliminates the identified problems. In order to demonstrate the consistency and the accuracy of the proposed

1. The final form of the existing h-formulation [9] is restrictive since it is applicable only when rl and rv are assumed to be constants. Further, the formulation relies upon an implicit relation between s and h in the two-phase region and it contains an additional diffusion-like term, which, being unconventional, has to be treated as a source in Eq. (14). Since its evaluation would lag the current calculation by one iteration, additional computation time would be required. 2. The proposed modification eliminates the problems associated with the existing h-formulation [9]. It can be applied even in the presence of substantial variations in rl and rv, without requiring any major modification, other than replacing them by their respective saturation values in the expressions for two-phase mixture properties. 3. The theoretical analyses in section 2 clearly demonstrates that the existing H-formulation [12] cannot be applied in its present form if rl and rv are strong functions of the local temperature, particularly when the temperature difference in the vapour phase is substantial. The required modifications have been also

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136

4.

5.

6.

7.

presented that indicate the necessity for iterative solution of T from the computed H. For steady one-dimensional problems with prescribed mass flow rate, the assumption of either constant or variable rl and rv does not affect the solution of h and hence T(x), although it would strongly affect the velocity solution and the resulting pressure drop in the evaporator. Such a simple conclusion, which is evident from the modified h-formulation, is not apparent from the H-formulation. For LTNE assumption, however, the solutions obtained for constant and variable rl and rv are expected to differ from each other, even for one-dimensional idealisation. For multi-dimensional problems with large temperature variations, however, the constant density assumption or even the Boussinesq approximation would be invalid, although it has been invoked by several researchers [12,21,44]. For such cases, the modified h-formulation should be used. Nearly identical solutions predicted by both h- and H-formulations for all tested cases demonstrate the consistency and the accuracy of the proposed modification. It has been also observed that the present formulation requires considerably less computation time as compared to that for the existing Hformulation. For multi-dimensional problems, the present formulation is expected to be even more beneficial with respect to the required computation time. All variants of TPMM, including the modified h-formulation, require the smoothing of effective diffusion coefficient [27], which justifies its necessity and the usefulness.

In view of the advantages and since the proposed h-formulation does not rely upon the definition of any synthetic variable like H [12], it is strongly recommended for the future use. Nevertheless, it is also apparent that the method should be further employed in order to simulate the multi-dimensional phase change problems inside porous media without invoking any simplifying assumption with respect to the single phase densities. Simulation of multidimensional phase change problems and the effects of density variation, however, have been left beyond the scope of the present study and should be taken up in the future, owing to its importance for more accurate modelling. It is also important to note that the modified h-formulation relies upon the validity of LTE assumption, according to which, artificial insulation layers are formed across the interfaces separating the single and the two-phase regions, primarily owing to the extremely low effective diffusion coefficients in the two-phase region for sz0 and sz1 [22]. Since any realistic prediction of the multi-dimensional complete phase change process within porous media requires the employment of LTNE assumption [22], the present formulation should be extended in order to accommodate such issues. It may, however, be expected that the proposed hformulation, along with LTNE assumption, should also perform better than the existing H-formulation and should be beneficial particularly for variable rl and rv, since other than the additional volumetric heat exchange term, the main features of the energy conservation equation for the fluid phase remain nearly the same [22]. As a final remark, it may be mentioned that when the existing Hformulation has been used along with the variable density assumption for the single phase regions after adopting the required modifications, in spite of the best efforts made by the present authors, no converged solution could be obtained since severe numerical oscillations in the residues for the energy conservation equation could only be detected. Therefore, no meaningful comparison in this respect could be reported in the article. Nevertheless, instead of identifying this observation as a real drawback of the

135

existing H-formulation with variable rl and rv and reporting it in the article, this issue has been left out for the future debate, particularly for the researchers who would still prefer to use the existing H-formulation, instead of the proposed modification. Acknowledgement One of us, Mr. Omar Rafae Alomar, would like to gratefully acknowledge the financial support from the German Academic Exchange Service (DAAD) and the Ministry of Higher Education and Scientific Research of Iraq for their scholarship grant to him. The authors also gratefully acknowledge the German Research Foundation (DFG) for supporting the subproject B02 of the Collaborative Research Centre CRC 920. References [1] Cheng P. Heat transfer in geothermal systems. In: Hartnett JP, Irvine Jr JF, editors. Advances in heat transfer, vol. 14. New York: Academic Press; 1979. 1e105. [2] Marie CM. Multiphase flow in porous media. Gulf, Houston. 1981. [3] Wang CY, Cheng P. Multiphase flow and heat transfer in porous media. In: Advances in heat transfer, vol. 30. New York: Academic Press; 1997. 93e196. [4] Bear J. Dynamics of fluids in porous media. New York: Elsevier; 1972. [5] Scheidegger AE. The physics of flow through porous Media. 3rd edn. Toronto: University of Toronto Press; 1974. [6] Ramesh PS, Torrance KE. Stability of boiling in porous media. Int J Heat Mass Transf 1990a;33:1895e908. [7] Ramesh PS, Torrance KE. Numerical algorithm for problems involving boiling and natural convection in porous materials. Numer Heat Transf B 1990b;17: 1e24. [8] Ramesh PS, Torrance KE. Boiling in a porous layer heated from below: effects of natural convection and a moving liquid two-phase interface. J Fluid Mech 1993;257:289e309. [9] Wang CY, Beckermann C. A two-phase mixture model of liquid-gas flow and heat transfer in capillary porous media, I. Formulation. Int J Heat Mass Transf 1993;36(11):2747e58. [10] Wang CY, Beckermann C, Fan C. Numerical study of boiling and natural convection in capillary porous media using the two-phase mixture model. Numer Heat Transf Part A 1994;26:375e98. [11] Peterson GP, Chang CS. Numerical heat analysis and evaluation for two-phase flow in porous-channel heat sinks. Numer Heat Transf Part A 1997;31(2): 113e30. [12] Wang CY. A fixed-grid numerical algorithm for two-phase flow and heat transfer in porous media. Numer Heat Transf Part B 1997;32:85e105. [13] Zhao TS, Liao Q. Mixed convective boiling heat transfer in a vertical capillary structure heated asymmetrically. J Thermophys Heat Transf 1999;13(3): 302e7. [14] Zhao TS, Cheng P, Wang CY. Buoyancy-induced flows and phase-change heat transfer in a vertical capillary structure with symmetric heating. Chem Eng Sci 2000;55:2653e61. [15] Najjari M, Nasrallah SB. Numerical study of boiling in inclined porous layer. J Porous Media 2003;6(1):71e81. [16] Najjari M, Nasrallah SB. Effects of latent heat storage on heat transfer in a forced flow in a porous layer. Int J Therm Sci 2007;47:825e33. [17] Najjari M, Nasrallah SB. Heat transfer between a porous layer and a forced flow: influence of layer thickness. Dry Technol 2009;27:336e43. [18] Li HY, Leong KC, Jin LW, Chai JC. Transient two-phase flow and heat transfer with localized heating in graphite foams. Int J Therm Sci 2010;49:1115e27. [19] Li HY, Leong KC, Jin LW, Chai JC. Transient behavior of fluid flow and transfer with phase change in vertical porous media. Int J Heat Mass Transf 2010;53: 5209e22. [20] Yuki K, Abei J, Hashizume H, Toda S. Numerical investigation of thermo-fluid flow characteristics with phase change against high heat flux in porous media. ASME J Heat Transf 2008;130. 012602-1e012602-12. [21] Shi JX, Wang JH. Numerical investigation of transpiration cooling with liquid coolant phase change. Transp Porous Media 2011;87:703e16. [22] Alomar OR, Mendes MAA, Trimis D, Ray S. Simulation of complete liquidvapour phase change inside porous evaporator using local thermal nonequilibrium model. Int J Therm Sci 2015;94:228e41. [23] Wei K, Wang J, Mao M. Model discussion of transpiration cooling with boiling. Transp Porous Media 2012;94:303e18. [24] He F, Wang J, Xu L, Wang X. Modeling and simulation of transpiration cooling with phase change. Appl Therm Eng 2013;58:173e80. [25] Xin C, Rao Z, You X, Song Z, Han D. Numerical investigation of vaporeliquid heat and mass transfer in porous media. Energy Convers Manag 2014;78:1e7. [26] He F, Wang J. Numerical investigation on critical heat flux and coolant volume required for transpiration cooling with phase change. Energy Convers Manag 2014;80:591e7. [27] Alomar OR, Mendes MAA, Trimis D, Ray S. Numerical simulation of complete

136

[28]

[29]

[30] [31] [32] [33] [34]

[35] [36]

S. Ray, O.R. Alomar / International Journal of Thermal Sciences 105 (2016) 123e136 liquid-vapour phase change process inside porous media using smoothing of diffusion coefficient. Int J Therm Sci 2014;86:408e20. Alomar OR, Mendes MAA, Trimis D, Ray S. Simulation of complete liquidvapor phase change inside divergent porous evaporator. Int J Mater Mech Manuf 2014;2(3):223e9. Werzner E, Laurinat M, Herrmann A, Ray S, Trimis D. Performance of synthetic geometries for computationally efficient modelling of flow through highly porous open-cell foams. In: International Conference on numerical and mathematical modeling of flow and transport in porous media e NM2 porous media, Dubrovnik, Croatia; September 29 e October 3, 2014. Kaviany M. Principles of heat transfer in porous media. 2nd edn. New York: Springer-Verlag; 1995. Borgnakke C, Sonntag RE. Fundamentals of thermodynamics. John Wiley & Sons; 2009. Dickson LE. A new solution of the cubic equation. Amer Math Mon 1898;5: 38e9. Kennedy EC. A note on the roots of a cubic. Amer Math Mon 1933;40:411e2. 9th printing. Dover, New York. In: Abramowitz M, Stegun IA, editors. Handbook of mathematical functions with formulas, graphs, and mathematical Tables; 1972. Berndt BC. Ramanujan's notebooks, part IV. New York: Springer-Verlag; 1994. lyi T. Cubic equations in polynomials and polynomial inBorwein P, Erde equalities. New York: Springer-Verlag; 1995.

€user; 1996. [37] King RB. Beyond the quartic equation. Boston, MA: Birkha [38] Mendes MAA, Ray S, Trimis D. A simple and efficient method for the evaluation of effective thermal conductivity of open-cell foam-like structures. Int J Heat Mass Transf 2013;66:412e22. [39] Mendes MAA, Ray S, Trimis D. Evaluation of effective thermal conductivity of porous foams in presence of arbitrary working fluid. Int J Therm Sci 2014;79: 260e5. [40] Mendes MAA, Ray S, Trimis D. An improved model for the effective thermal conductivity of open-cell porous foams. Int J Heat Mass Transf 2014;75: 224e30. [41] Mendes MAA, Ray S, Trimis D. Determination of effective thermal conductivity of open-cell porous foams for arbitrary working fluid using single point measurement. In: Proc. 22nd National & 11th ISHMT-ASME Heat and Mass Transfer Conf. Kharagpur, India: Indian Institute of Technology; December 28 e 31, 2013. [42] Patankar SV. Numerical heat transfer and fluid flow. Hemisphere, New York. 1980. [43] Incropera FP, Dewitt DP, Bergman TL, Lavine AS. Fundamentals of heat and mass transfer. John Wiley & Sons; 2007. [44] Lindner F, Mundt Ch, Pfitzner M. Fluid flow and heat transfer with phase change and local thermal non-equilibrium in vertical porous media. Transp Porous Media 2015;106:201e20.