- Email: [email protected]
Smoothed iterative enthalpy approach for solid-liquid phase change
International Journal of Thermal Sciences xxx (xxxx) xxx
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts
Smoothed iterative enthalpy approach for solid-liquid phase change R.R. Kasibhatla *, D. Brüggemann Chair of Engineering Thermodynamics and Transport Processes (LTTT), Center of Energy Technology (ZET), University of Bayreuth, Universit€ atsstraße 30, 95440, Bayreuth, Germany
A B S T R A C T
Numerical modelling of solid-liquid phase change problems has a wide range of applications. Here, the fictitious source based method is one of the most implemented iterative approaches to account for non-linear temperature-enthalpy coupling despite its slow convergence. In the present work a method based on the optimum approach is implemented to model the solid-liquid phase change with a fast convergence. The complete numerical model involving the heat transfer and fluid flow is implemented in the open source computational fluid dynamics software OpenFOAM. The numerical results have fairly agreed with the benchmark experiments from Gau and Viskanta. Parameters like dimensional simplification, porosity constants and the influence of temperature correction step are studied to understand the consistent deviations in the numerical modelling. The optimum approach has also shown a fast convergence compared to the fictitious source-based method. For thin melting ranges, the temperature correction step ensures energy conservation.
1. Introduction In recent years, the number of studies on phase change materials (PCM) has drastically grown due to a large choice of materials and their wide range of usage in thermal storage applications. Alongside studying such PCM applications, it is also necessary to understand the melting behaviour and patterns of the PCM. The melting behaviour of a PCM can be studied either experimentally or through numerical simulations. Here we focus on numerical modelling and simulation, as it is a flexible and cost-effective way to study parameter variation. Meanwhile, the nu merical modelling of solid-liquid phase change has also emerged as a broad field due to applications in metallurgical processes, aerospace and thermal storage units [1]. Different approaches exist to model melting and solidification. One way to classify these approaches is to distinguish between moving and fixed grid methods. Fixed grid methods are more flexible and easier to implement but the phase change boundary can be smeared over a few cells. Apart from the fixed grid approaches, the lattice-Boltzmann method is also a very good alternative to model the phase change, but not very popular [2]. In all the approaches, energy equation becomes non-linear and multiple methods have been devel oped to account for this non-linearity. The most common source based enthalpy method [3] demands a large number of iterations to converge the energy equation. To overcome this disadvantage, Swaminathan and Voller [4] came up with an optimum approach to model the solid-liquid phase change more efficiently. This method was described in a gener alized form to support its adaptability into the well known source-based method or an apparent heat capacity method [5]. Additionally, a
comparative study on the computational effort between the different enthalpy methods was also presented there including the comparison between different enthalpy-temperature curves [5]. This method how ever has drawn very less attention due to the complexity in imple menting the algorithm. However, some researchers [6,7] have compared this approach with other methods. Nevertheless, the further continuity of this research was not considerably substantial. Challenges of the proposed optimum approach by Swaminathan and Voller [5] lie in computing the slope of the enthalpy-temperature curve dH dT for every iteration step. Pham [6] treated this by computing a time-averaged and space-averaged slope of the enthalpy-temperature curve dH dT . But, this is still limits the strong coupling of the enthalpy with temperature for a given phase change problem. For computations involving complex ge ometries and sharp phase change curves with thin melting ranges, such slope formulations dH dT may demand smaller time steps and conformal computational grid. Therefore, the present article introduces a smooth slope function for an iterative optimum approach to increase stability and reduce computational effort. Apart from the new adaptation of the optimum approach, this article also studies the influence of different simulation parameters like corrector steps, Darcy term constants and choice of smooth functions on numerical results. Suppression of velocity in the solid phase is also necessary to compute a solid-liquid phase change problem. Morgan [8] tried to implement this by overwriting the velocity in solid phase to zero. However, this can lead to instabilities in advanced applications. Variable viscosity method (VVM) [9–11] is also an alternative to model the ve locity suppression in solid phase of a melting problem. In the present
* Corresponding author. E-mail addresses: [email protected] (R.R. Kasibhatla), [email protected] (D. Brüggemann). https://doi.org/10.1016/j.ijthermalsci.2019.106187 Received 20 August 2019; Received in revised form 20 October 2019; Accepted 19 November 2019 1290-0729/© 2019 Published by Elsevier Masson SAS.
Please cite this article as: R.R. Kasibhatla, https://doi.org/10.1016/j.ijthermalsci.2019.106187
D.
Brüggemann,
International
Journal
of
Thermal
Sciences,
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
Volume of the cell, m3
Nomenclature
V
a b c capp d D Dc Ds erf g H L p q_ S sech T t tanh Δt u
Greek symbols α Volume fraction of the liquid β Thermal expansion coefficient, 1/K γ Porosity in the Kozeny-Carman equation ε Tolerance criterion for convergence κ Thermal conductivity, W/(mK) λ Enthalpy jump in a melting range, m2/s2 μ Dynamic viscosity of the fluid, kg/(ms) ρ Density, kg/m3 σ Square root of variance in the peak function, K ω Under relaxation factor
Discretized energy equation coefficient Discretized energy equation coefficient Mass specific heat capacity, J/(kg K) Mass specific apparent heat capacity, J/(kg K) Discretized energy equation coefficient Darcy source term, kg/(m3s) Morphological constant in Darcy term, kg/(m3s) Zero non-divisional constant in Darcy term Error function Gravitational acceleration, m/s2 Mass specific enthalpy, m2/s2 Latent heat, m2/s2 Pressure, kg/(m s2) Heat flux, W/(m2) Source term Hyperbolic secant function Temperature, K Time, s Hyperbolic tan function Time step, s Velocity, m/s
Subscripts l m n nb nþ 1 o p s
work, the most used porosity scheme with a Darcy source term [12] is employed to account the solid velocity suppression. In recent times, Faden et al. [13] implemented such a porosity approach to model the settling of solid phase in a macro-capsule with an implicit approach for settling. Whereas, Gudibande and Iyer [14] proposed another cell-splitting method to model a similar contact melting. Another nu merical approach was proposed by Kozak and Ziskind [15] to model the melting and settling of solid, which employs a peak function to restrict the solid bulk motion.
D¼
Physical phenomena in a thermo-fluid process can be represented by the conservative equations for mass, momentum and energy [16]. During heat transfer, small density change in the liquid phase influenced by a minute temperature difference can be assumed linear. In the present work, this small density change in a natural convective melting process is accounted by the Bousinessq approximation source term in the mo mentum conservation equation [17]. The equation set solved to repre sent the thermo-fluid process in a solid-liquid phase change is Continuity equation: Momentum equation: � ∂u ρ þ u ⋅ ru ¼ rp þ μΔu þ D ⋅ u þ ρg⋅ð1 ∂t
c ¼ α ⋅ cl þ ð1
�
Energy equation: � ∂H ρ þ r ⋅ ðuHÞ ¼ r⋅ðκrTÞ: ∂t
Tm ÞÞ
(4)
αÞ⋅cs :
(5)
The specific heat capacity c in this work contributes to the sensible enthalpy and the apparent heat capacity ceff is a change of total enthalpy with respect to the temperature change dH . In this work, it is represented dT as � � dH L T Tm ceff ¼ ¼ c þ ⋅sech2 : (6) dT 2σ σ
(1)
β ⋅ ðT
αÞ2 : α3 þ Ds
ð1
where, α is the volume fraction of liquid PCM, Dc is the Darcy constant and Ds is a small constant to avoid the division with zero. In Eq. (3), H is the enthalpy of a material and can be segregated into sensible and latent parts. Approaches like the source based method [3] use the segregated representation of enthalpy. However, the enthalpy can also be repre sented as a function of temperature H ¼ fðTÞ [21]. Theoretically, melting is an isothermal process for pure metals, whereas for binary alloys it occurs over a temperature range. In this work, melting over a very small temperature range is realized in the numerical simulations. The apparent heat capacity of the melting pro cess shows a large change over a small temperature range ΔTm [22]. In this article, the phase change is mathematically represented by a peak shape function. Here, the apparent heat capacity ceff is different from c in Eq. (5), which is
2. Numerical model
r⋅u ¼ 0
Dc ⋅
and superscripts Liquid Melting point Old iteration level Neighbouring cells Actual iteration level Old time step Node of the computing cell Solid
(2)
Here, σ is the square root of the variance of the peak shape function. Where, σ ¼ ΔT4m and ΔTm is the difference between solidus and liquidus temperature Ts and Tl respectively. As shown in Fig. 1, Eq. (6) is rep resented linear in the sensible parts and experiences a peak shape at the melting temperature range. The volume fraction of the melted portion α is formulated as a function of temperature. The volume fraction of the melted portion α is obtained from the integral function of apparent heat capacity written as Z T Z T dH dT ¼ dH ¼ ceff dT: (7) Tref dT Tref
�
(3)
Here, ρ is the density and μ the dynamic viscosity of the fluid, κ the thermal conductivity, u the velocity, p the pressure, β the thermal expansion coefficient of the PCM in liquid phase, D the Darcy source term, T the temperature and Tm the melting temperature. The Darcy source term in the momentum equation suppresses the velocity of solid phase in the computational domain, whereas the interfacial mushy zone between solid and liquid phase is treated as a porous medium [18,19]. This source term in Eq. (2) reads [12,20]. 2
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
This equation can be rewritten in an iterative form as [4]. � � �� � � � � � 1� dH n o n H H T nþ1 ¼ d H ap þ d p p p p þ S: dT n
(10)
The source of fluxes from neighbouring cells S is S¼
X anb T nþ1 nb þ b: In Eq. (10), ðHnp Þ
(11) 1
is the temperature which is the inverse function of
enthalpy at the iteration step n i.e. ðHnp Þ ðρVÞp . Δt
1
¼ T np and d is a discretization
coefficient with d ¼ In Eqs. (10) and (11), n þ 1 is the current iteration step, n is the old iteration step, o is the old time step, Δt is the time step size, nb represent the neighbouring cells of the computing cell p, ðVÞp is the volume of the computing cell p and a as well as b are the coefficients of the discretized energy equation. The numerical model in the present work employs a smooth slope function instead of a slope computation in the each iterative step. Here, instead of computing the slope of the H-T curve as [21]. � � dH Hn Hn 1 ¼ n : (12) dT n T Tn 1
Fig. 1. Apparent heat capacity ceff vs. temperature.
The present numerical model computes the slope of the H-T curve or the apparent heat capacity as � � � n � dH L T Tm ¼ c þ ⋅sech2 : (13) dT n 2σ σ
Which is shown in Fig. 2. By rearranging Eq. (7) and Eq. (6) in the enthalpy H, the liquid volume fraction α is obtained by � � T Tm α ¼ 0:5 þ 0:5⋅tanh : (8)
The new proposed step as in Eq. (13) reduces the numerical in stabilities caused by the slope computation as in Eq. (12). When the temperature change between the temperature at the old iterative step Tn and the pre old iterative step Tn 1 is zero, a division by zero occurs for Eq. (12). However, this can be treated by using a small non-zero additive value in the denominator of Eq. (12). But, this small non-zero additive can affect the maximum range of the apparent heat capacity alongside creating robustness and accuracy issues. Thereby, the proposed smooth function in this numerical model treats these robustness and accuracy issues by computing a definitive slope value for every temperature point as written in Eq. (13). In the numerical model, Eq. (10) with the slope values from Eq. (13) is iterated until the consistency of enthalpytemperature coupling is achieved in every time-step [4]. Inside an every time-step,
σ
As shown in Fig. 2, the increase in liquid volume fraction α occurs over a small melting temperature range ΔTm . 2.1. Enthalpy-temperature model The proposed optimum approach is a computationally inexpensive alternative to model the solid-liquid phase change [4]. This numerical model converges faster than the standard source based fictitious method and the apparent heat capacity method [4]. Within the approach, the discretized form of the energy equation Eq. (3) is described by i h X ap Tp ¼ (9) anb Tnb þ d H op Hp þ b:
1. Initially, the old values of enthalpy and temperature, Hop and T op respectively are updated to Hnp ¼ Hop and T np ¼ T op . 2. Eq. (10) is solved for the updated enthalpy and temperature values. By solving Eq. (10), the new temperature Tnþ1 for the n þ 1 iteration p and the enthalpy Hnp for the iteration level n is obtained. 3. In an iteration step, the explicit form of the enthalpy Hnp is updated � � using the new temperature Tnþ1 and the slope of the line dH in p dT n
equation H nþ1 ¼ H np þ p
� � h dH ⋅ T nþ1 dT n p
i T np :
(14)
Inside a time step, the steps 2 and 3 are repeated until convergence, i. e. the consistency of temperature with enthalpy. The convergence of the algorithm is characterised by a convergence tolerance ε: T nþ1 p
T np � ε:
(15)
The enthalpy update step in Eq. (14) is proposed by Swaminathan and Voller [4]. In the present article, an under relaxation factor w is used to relax the sharp enthalpy update from the new temperature: � � h i dH H nþ1 (16) ¼ H np þ w⋅ ⋅ T nþ1 T np : p p dT n
Fig. 2. Volume fraction of liquid, α vs. temperature of PCM over a phase change domain. 3
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
For optimum approach, Voller [5,21] proposed a temperature correction after the enthalpy update in Eq. (14). The correction step ensures the enthalpy relation with the temperature, which reads 8 H nþ1 p > ; H nþ1 < cs Ts > p > > cs > > > > " # > > > H nþ1 cs Ts > p > < Tm þ σtanh 1 cs Ts � H nþ1 p c T þ L cs Ts nþ1 l l Tp ¼ > > > > > H nþ1 � cl Tl þ L; > p > > > > > > > H nþ1 L : p ; H nþ1 > cl Tl þ L: p cl
The rectangular cavity had dimensions of 88.9 mm in height, 63.5 mm in width and 38.1 mm in depth. The copper plates arranged at the left and right boundaries ensure the implementation of a constant homogenous surface temperature. The top and bottom of the rectangular cavity were insulated [24]. An identical geometrical configuration is developed for the numerical setup. The material properties and boundary conditions used in the numerical set-up are summarized in Table 1. In numerical simulations, the nature of melting also depends on the resistance offered against velocity in solid phase, see Kasibhatla et al. [11]. Thereby, the Darcy term parameters similar to Brent et al. [3] have been implemented, where Dc ¼ 1:6⋅106 kg/m3s and Ds ¼ 10 3 . In nu merical simulations, for the 2D case front and back walls of the cavity are assumed symmetric and as walls for the 3D case. Grid dependence is studied to depict the influence of numerical model on the grid refine ment. A thin boundary layer of 2 mm is allocated to resolve the large heat fluxes. Apart, the rest of geometry is homogeneously refined into structured hexahedral cells. Cell refinement combinations of 87�62, 131�93 and 174�123 cells in X and Y direction respectively are studied. Between the difference in melt fronts obtained, the computational grid with 87�62 cells is employed in this work. For the 2D assumption only one cell is contributed in the Z direction. Whereas, for a 3D geometrical configuration a numerical grid with 441,870 cells contributing 87�62�55 cells in X, Y and Z directions is employed. Also, a wall boundary is imposed on the front and back boundaries of the rectangular cavity for the 3D case. Convergence criterion ε ¼ 10 6 is set for all the numerical simulations in this work.
(17)
The nature of enthalpy transport with the convective term is dependent on the solid-liquid phase change problem. For pure melting with no mushy zone, sensible enthalpy is modelled in the convective term [21]. This can be derived by assuming the continuity of liquid phase after segregating the total enthalpy into sensible and latent parts [23]. For an assumption, where the velocity of solid phase can be defined to zero. The convective term formulated with the apparent heat � � capacity r⋅ u dH T is derived into as r⋅ðucl TÞ [21]. In the present work, dT to mimic the pure melting of gallium, a very thin negligible mushy zone is assumed, alongside no movement of solid particles in the thin mushy zone. Therefore, in this work only sensible enthalpy is transported in the convective term of energy equation written as r⋅ðucTÞ.
4. Results and validation of the numerical model
3. Benchmark experiments and numerical set-up
A melting range of ΔTm ¼ 0.002 K is employed to mimic a pure isothermal phase change. The assumption of a 2D geometry results in computing high magnitude convective vortices. These vortices split up into two or more vortices leading to a deeply curved phase change boundary, also see Giangi and Stella [28]. This behaviour raised ques tions if the flow was mono-cellular or multi-cellular [26]. This influence of the velocity field may lead to an enhanced melting of the substance leading to such melting patterns. It could be noted from the literature that this occurs due to the influence of 2D geometrical configuration, see Nitrityuk [29]. Therefore, a 3D case with the front and back side of the cavity considered walls is computed [25]. As seen in Figs. 4 and 3D numerical results represent the melting pattern in reality very well and similar to the experimental findings of Gau and Viskanta [24]. The small deviation of numerical results from experiments can be observed in time t ¼ 10, 17 min in Fig. 4. The phase front from the 3D geometry experi ences a smooth continuous profile unlike the curvatures due to several vortices as seen in Fig. 4. The difference in nature of vortices in a 2D and 3D geometry is discussed by Nikrityuk [29]. The numerical error be tween the energy transported into the computational domain and energy used for phase change is investigated. This error made by the present optimized approach is approximately 0.01 %. In general, the 3D simulations are computationally very expensive than 2D simulations, which also increases the efforts to study phase change numerically. But with the present smoothed optimum approach,
Investigating the melting of 99.6 % pure gallium inside a rectangular cavity by Gau and Viskanta [24] are one of the initial experimental benchmarks in the field of solid-liquid phase change. Different numer ical models developed in the recent years [7,25] were also validated with these experimental results. Despite some strong criticism [26], the issue related to numerical convergence is resolved [27]. Validation of a numerical model with these benchmark results does not just provide an absolute accuracy of the numerical model. But, also provides a relativity of other numerical approaches validated against these numerical results. A simple graphical representation of the experimental set-up is shown in Fig. 3. The rectangular cavity in the experimental set-up is filled with a very pure gallium. The convective melting of gallium inside in the cavity is measured at different instances of time.
Table 1 Parametric tabular of the numerical set-up. Material properties Density (ρ) Viscosity (μl )
4
kg/m Pa⋅ s
values 3
6093 0.00181
Specific heat capacity (c) Latent heat (L) Thermal conductivity (κ) Melting temperature (Tm )
J/kgK J/kg W/mK K K
311.15
Cold boundary temperature (Tc )
K
301.45
Heated boundary temperature (Th )
Fig. 3. Schematic diagram depicting the experimental set-up of rectangular cavity from Gau and Viskanta [24].
units
381.5 80160 32.0 302.93
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
5.1. Temperature corrector-step The numerical results with and without the temperature correction after the enthalpy update in the algorithm is studied and the comparison are shown in Fig. 6. Despite, the similar melt fronts in Fig. 6, for thin melting ranges the numerical model without a temperature corrector step violates the energy equation to a large extent. The numerical error denoted as the difference between heat trans ferred into the computational domain and the enthalpy change is compared for uncorrected and corrected models in Table 2. It can be noticed in Table 2 that the error made by the numerical model without a corrector step gradually increases with a decrease in the melting range. For thin pure melting ranges, the error made is very large. On the other hand, for melting ranges greater than 1 K the error made by the nu merical model without a temperature corrector step is relatively negli gible. This can also be noticed in the volume fraction of liquid PCM α estimated by the numerical model as in Table 3. Various factors like the size of finite volume, temperature gradient and the expansion coefficient also play a role in the numerical error. Therefore, a generalization of such a melting range value cannot be made to assert where the tem perature corrector step really is influential. To analyse this, different parameters need to studied to understand the sensitivity of temperature corrector step.
Fig. 4. Validation of the numerical model using a 3D numerical grid with the benchmark experiments [24].
the numerical effort to simulate such 3D cases including computing phase change in complex geometries is treated well. Different commercial software employ the standard source based method, due to its stability [30]. But, optimum approach with the stable slope function dH dT is a decent alternative to reduce computational effort alongside providing stability. The computational expense of the present smoothed optimum approach is compared with the fictitious source based method [3,11]. In each time step the iterations required by the present model with a temperature correction step in comparison with the source based method is shown in Fig. 5.
5.2. Choice of peak function In this section, two different slope functions are compared to un derstand the influence of numerical model on its dH dT slope function. However, at convergence the product of apparent heat capacity and � � temperature difference dH and ½Tnþ1 ðHnp Þ 1 � becomes smaller and p dT n
non-influential. This is proved by varying the choice of apparent heat capacity. For this, a Gaussian normal distribution function is employed to represent the peak behaviour of the apparent heat capacity as shown in Fig. 1. For this comparison, the value of σ ¼ ΔT4 m is set same for both the functions. Mathematically, the total integral of different peak functions remain same independent of the variance thickness σ . But, in numerical simulations a higher value in the denominator ensures a decent compression of the solid-liquid interface. Whereas, a lower value in the denominator results in a broader solid-liquid interface. Accordingly, the Gaussian normal distribution function is written in the form
5. Investigating the parametric influence Different parameters involved in the numerical model are investi gated to understand their influence on the simulation results. Numerical model parameters like the efficiency of temperature corrector step from Eq. (17), Darcy term constants and choice of the peak function are studied in this case. Due to a high computational effort by the 3D case, a 2D case with 87�62 cells is employed to study the parametric influence of the numerical results.
Fig. 6. Comparison of the melting fronts at time t ¼ 2 min and t ¼ 17 min with and without temperature correction after enthalpy update for different melting ranges.
Fig. 5. Iterations required by the optimum approach and source based method to converge shown on a logarithmic scale Eq. (15). 5
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
pressure drop. However, the empirical Kozeny-Carman equation is
Table 2 Numerical error for uncorrected and corrected temperature steps. ΔTm ðKÞ
Error T-uncorrected %
Error T-corrected %
0.01 0.1 0.25 1 2
15 2 0.26 0.191 0.159
0.0002 0.00015 0.002 0.17 0.2
rp ¼
α with T-uncorrected
α with T-corrected
0.01 0.1 0.25 1 2
0.8367 0.6817 0.6761 0.6620 0.6506
0.6704 0.6735 0.6717 0.6639 0.6501
� ceff ¼
dH ¼cþe dT
�2 ðT Tm Þ 2⋅σ
L ⋅ pffiffiffiffiffi: σ 2π
(20)
where γ is the porosity of bed, μ is the viscosity of fluid entering the packed bed, ϕ is the sphericity of packed bed particles and Db their diameter of the particles in bed. Accordingly, this equation is rewritten as Eq. (4) in the solid-liquid phase change problems. In Eq. (4), Dc is the morphological Darcy constant and Ds is defined a small value restricting the division by zero [3]. But, the selection of these values is still inconclusive. Therefore, in this article the constants Dc and Ds are varied such that the total resistance value in the solid phase remains same as in Brent et al. [3]. By substituting, selected values as Dc ¼ 1:6⋅106 kg/m3s and Ds ¼ 10 3 , the pressure gradient in the solid phase becomes 1:6⋅109 times the velocity and 0 in the liquid phase. But, this value of pressure gradient can be achieved by different choice of Dc and Ds values. These combinations of constants Dc and Ds are presented to achieve these absolute values in both the phases as shown in Fig. 8. Hereby, the pressure gradient in the mushy zone vary between them. For this, all the other simulation parameters are maintained constant as in Table 1 except the Darcy constants Dc and Ds . Fig. 9 explains the importance of total resistance offered by the Darcy term i.e. the product of Dc and D1s . Various literatures in the field of solid-liquid phase change
Table 3 Volume fraction of liquid phase for the uncorrected and corrected temperature steps at time t ¼ 20 min. ΔTm ðKÞ
150μ ð1 γÞ2 ⋅u: γ3 ϕ2 D2b
(18)
report Ds as a constant restricting the division by zero. Additionally, another numerical study is executed to depict the in fluence of total resistance on the melting patterns. For this, the numer ator Dc has been varied from Dc ¼ 1:6⋅106 to Dc ¼ 1:6⋅103 and Dc ¼ 1:6⋅1010 . Here in this study, the small non-division constant Ds is left constant at 10 3 . The results are shown in Fig. 10. The results of numerical benchmark as seen in Figs. 9 and 10 sup ports the argument that total resistance offered in the solid-phase of momentum equation is important to represent a melting pattern rather than the independent constants Dc or Ds . But, the same can play a pivotal role in modelling binary alloys and impure materials. Thereby, a need to study the nature of resistances offered in the fluid-solid combined mushy zone exists.
Accordingly, the integral of this function i.e. the liquid volume fraction is represented by an error function written as � � T Tm pffiffi : (19) α ¼ 0:5 þ 0:5⋅erf σ 2 Thereby, numerical results between Gaussian normal distribution and secant function formulation are compared in Fig. 7. As seen in Fig. 7, difference between the numerical results is very small and can be neglected. Therefore, the selection of apparent heat capacity has rather smaller influence on the solution accuracy but only on the stability. 5.3. Darcy term Darcy term implemented to suppress the velocity in solid phase is inspired from the Kozeny-Carman equation used to model the pressure drop of fluid in the packed beds [12,20]. In solid-liquid phase change, the Kozeny-Carman equation assists the suppression of velocity advec tion in solid phase. This resistance offered is implemented in the mo mentum equation with constants like Dc and Ds enabling a typical
6. Discussion and conclusion
Fig. 7. Comparison between different numerical models for the choice of the apparent heat capacity approximate.
Fig. 8. Plot showing different progressions of Darcy term with different values of Darcy constants Dc and Ds .
An optimized method with a smooth continuous phase change curve is presented to compute solid-liquid phase change. A small error with very low computational time is the topic of interest here. The numerical results employing the 3D geometry have shown a very good agreement with the experiments. Whereas, the 2D geometry case with the
6
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
specific heat capacity for phase change over temperature range. For this, as seen in Fig. 7 the choice of peak function has shown no influence on the numerical results. This is because the optimum algorithm employs the specific heat capacity function to achieve the enthalpy change till the temperature difference between the old and new iteration step becomes small, as described in Eq. (14) and Eq. (15). Another numerical parameter namely the porosity resistance to model solid-liquid phase change has been investigated in this work. Here, initially for the choice of constants Dc and Ds are chosen in such a way to achieve the same momentum resistance in the solid-phase. Re sults for numerical simulations with these different choice of constants are compared in Fig. 9. The results show a similar melt front for different choice of Dc and Ds values when the total resistance offered in solid phase is the same. This argument is further supported by an additional numerical study depicted in Fig. 10, where the value of Darcy constant Dc is varied by keeping the small non-division constant Ds the same. Different melt fronts shown in Fig. 10 depict the influence of total resistance on melting fronts. Hereby, when modelling the solid-liquid phase change of binary alloys or impure materials the particle distri bution in the interface region needs to be understood. For the same solid phase resistance as seen in Fig. 8 the source term curves takes different values for various mixture volume fractions. The amount of resistance necessary to restrict the motion of solid particle between 0 and 100% volume fraction moving in the molten liquid at the interface is still unknown. Because, for different combinations of Dc and Ds the mo mentum resistance in solid-liquid mixture at the interface takes distant values. Therefore, a detailed study on the choice of these values is necessary to study the dendrite growth for binary alloys in moving boundary problems.
Fig. 9. Comparison of numerical results with different Dc and Ds combinations but same total resistance in the full solid phase.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.ijthermalsci.2019.106187. References [1] Z. Liu, Y. Yao, H. Wu, “Numerical modeling for solid–liquid phase change phenomena in porous media: shell-and-tube type latent heat thermal energy storage, Appl. Energy 112 (2013) 1222–1232. [2] S. Chakraborty, D. Chatterjee, “An enthalpy-based hybrid lattice-Boltzmann method for modelling solid–liquid phase transition in the presence of convective transport, J. Fluid Mech. 592 (2007) 155–175. [3] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numer. Heat Tran. 13 (3) (1988) 297–318. [4] C. Swaminathan, V. Voller, A general enthalpy method for modeling solidification processes, Metall. Trans. B 23 (5) (1992) 651–664. [5] C. Swaminathan, V. Voller, On the enthalpy method, Int. J. Numer. Methods Heat Fluid Flow 3 (3) (1993) 233–244. [6] Q. Pham, Comparison of general-purpose finite-element methods for the stefan problem, Numer. Heat Tran. 27 (4) (1995) 417–435. [7] A. K€ onig-Haagen, E. Franquet, E. Pernot, D. Brüggemann, A comprehensive benchmark of fixed-grid methods for the modeling of melting, Int. J. Therm. Sci. 118 (2017) 69–103. [8] K. Morgan, A numerical analysis of freezing and melting with convection, Comput. Methods Appl. Mech. Eng. 28 (3) (1981) 275–284. [9] Y. Asako, M. Faghri, M. Charmchi, P.A. Bahrami, Numerical solution for melting of unfixed rectangular phase-change material under low-gravity environment, Numer. Heat Tran., Part A: Applications 25 (2) (1994) 191–208. [10] D.K. Gartling, Finite element analysis of convective heat transfer problems with change of phase, Comput. Method Fluids (1980) 257–284. [11] R.R. Kasibhatla, A. K€ onig-Haagen, F. R€ osler, D. Brüggemann, Numerical modelling of melting and settling of an encapsulated pcm using variable viscosity, Heat Mass Tran. (2016) 1–10. [12] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, Int. J. Heat Mass Tran. 30 (8) (1987) 1709–1719. [13] M. Faden, A. K€ onig-Haagen, S. H€ ohlein, D. Brüggemann, An implicit algorithm for melting and settling of phase change material inside macrocapsules, Int. J. Heat Mass Tran. 117 (2018) 757–767. [14] N. Gudibande, K. Iyer, Numerical simulation of contact melting using the cellsplitting modified enthalpy method, Numer. Heat Tran., Part B: Fundamentals 71 (1) (2017) 84–107.
Fig. 10. Depicting the numerical results obtained by varying Darcy constants Dc and letting the non-division constant Ds in the denominator remain constant at Ds ¼ 0.001.
symmetric boundary conditions on front and back walls is employed to investigate the numerical model parameters in this work. Contrast to the source based method, the present numerical model is a fast converging iterative approach which needs five to ten iterations in each time step to converge. This benefit assists with effortless solving of phase change problems for dense computational grids. Alongside, different simulation parameters were also studied in the present work. Here, the temperature correction after the enthalpyupdate ensures energy conservation in every time-step. Here, the very small numerical error in the energy equation results in the distorted phase front in the numerical simulations. This numerical error is large for very thin melting ranges rather than the broad melting ranges as seen in Table 2 and Table 3. As seen in Fig. 6 for a similar melting range ΔTm ¼ 1 K, the numerical results with and without temperature correc tion show a very different phase fronts. Here, it can also be seen that the numerical error increases with the decrease in melting range i.e. tem perature re-correction after an enthalpy update is very necessary when modelling pure melting problems. Whereas, for binary alloy melting problems the same can be tolerated with a small numerical error, if a very low computational expense need to be achieved. Furthermore, different mathematical functions were also assumed to model the 7
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx [22] M. Reading, D. Elliott, V. Hill, A new approach to the calorimetric investigation of physical and chemical transitions, J. Therm. Anal. Calorim. 40 (3) (1993) 949–955. [23] J. Ni, C. Beckermann, A volume-averaged two-phase model for transport phenomena during solidification, Metall. Trans. B 22 (3) (1991) 349. [24] Gau and Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Tran. 108 (1) (1986) 174–181. [25] O. Ben-David, A. Levy, B. Mikhailovich, A. Azulay, 3d numerical and experimental study of gallium melting in a rectangular container, Int. J. Heat Mass Tran. 67 (2013) 260–271. [26] J.A. Dantzig, “Modelling liquid–solid phase changes with melt convection, Int. J. Numer. Methods Eng. 28 (8) (1989) 1769–1785. [27] N. Hannoun, V. Alexiades, T.Z. Mai, Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side, Numer. Heat Tran. Part B: Fundamentals 44 (3) (2003) 253–276. [28] M. Giangi, Fulvio Stella, “Melting of a pure metal on a vertical wall: numerical simulation, Numer. Heat Tran. Part A: Applications 38 (2) (2000) 193–208. [29] P.A. Nikrityuk, Computational Thermo-Fluid Dynamics: in Materials Science and Engineering, John Wiley & Sons, 2011. [30] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, Int. J. Heat Mass Tran. 50 (9–10) (2007) 1790–1804.
[15] Y. Kozak, G. Ziskind, Novel enthalpy method for modeling of pcm melting accompanied by sinking of the solid phase, Int. J. Heat Mass Tran. 112 (2017) 568–586. [16] S.V. Patankar, Numerical Heat Transfer and Fluid Flow”, Hemisphere Publishing Corp., Washington, DC, 1980. [17] J. Boussinesq, Th�eorie de l’�ecoulement tourbillonnant et tumultueux des liquides dans les lits rectilignes a grande section, Gauthier-Villars et Fils, Paris, 1897. [18] V. Voller, N. Markatos, M. Cross, Solidification in convection-diffusion, in: N. Markatos, M. Cross, D. Tatchell, N. Rhodes (Eds.), Numerical Simulation of Fluid Flow and Heat/Mass Transfer Processes, vol. 18, Springer Berlin Heidelberg, 1986, pp. 425–432, of Lecture Notes in Engineering. [19] F. R€ osler, Modellierung und Simulation der Phasenwechselvorg€ ange in makroverkapselten latenten thermischen Speichern”. PhD thesis, University of Bayreuth, Logos Verlag Berlin, 2014. [20] P. Xu, B. Yu, “Developing a new form of permeability and kozeny–carman constant for homogeneous porous media by means of fractal geometry, Adv. Water Resour. 31 (1) (2008) 74–81. [21] W. Minkowycz, E. Sparrow, J. Murthy, Handbook of Numerical Heat Transfer, Wiley, 2006.
8
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts
Smoothed iterative enthalpy approach for solid-liquid phase change R.R. Kasibhatla *, D. Brüggemann Chair of Engineering Thermodynamics and Transport Processes (LTTT), Center of Energy Technology (ZET), University of Bayreuth, Universit€ atsstraße 30, 95440, Bayreuth, Germany
A B S T R A C T
Numerical modelling of solid-liquid phase change problems has a wide range of applications. Here, the fictitious source based method is one of the most implemented iterative approaches to account for non-linear temperature-enthalpy coupling despite its slow convergence. In the present work a method based on the optimum approach is implemented to model the solid-liquid phase change with a fast convergence. The complete numerical model involving the heat transfer and fluid flow is implemented in the open source computational fluid dynamics software OpenFOAM. The numerical results have fairly agreed with the benchmark experiments from Gau and Viskanta. Parameters like dimensional simplification, porosity constants and the influence of temperature correction step are studied to understand the consistent deviations in the numerical modelling. The optimum approach has also shown a fast convergence compared to the fictitious source-based method. For thin melting ranges, the temperature correction step ensures energy conservation.
1. Introduction In recent years, the number of studies on phase change materials (PCM) has drastically grown due to a large choice of materials and their wide range of usage in thermal storage applications. Alongside studying such PCM applications, it is also necessary to understand the melting behaviour and patterns of the PCM. The melting behaviour of a PCM can be studied either experimentally or through numerical simulations. Here we focus on numerical modelling and simulation, as it is a flexible and cost-effective way to study parameter variation. Meanwhile, the nu merical modelling of solid-liquid phase change has also emerged as a broad field due to applications in metallurgical processes, aerospace and thermal storage units [1]. Different approaches exist to model melting and solidification. One way to classify these approaches is to distinguish between moving and fixed grid methods. Fixed grid methods are more flexible and easier to implement but the phase change boundary can be smeared over a few cells. Apart from the fixed grid approaches, the lattice-Boltzmann method is also a very good alternative to model the phase change, but not very popular [2]. In all the approaches, energy equation becomes non-linear and multiple methods have been devel oped to account for this non-linearity. The most common source based enthalpy method [3] demands a large number of iterations to converge the energy equation. To overcome this disadvantage, Swaminathan and Voller [4] came up with an optimum approach to model the solid-liquid phase change more efficiently. This method was described in a gener alized form to support its adaptability into the well known source-based method or an apparent heat capacity method [5]. Additionally, a
comparative study on the computational effort between the different enthalpy methods was also presented there including the comparison between different enthalpy-temperature curves [5]. This method how ever has drawn very less attention due to the complexity in imple menting the algorithm. However, some researchers [6,7] have compared this approach with other methods. Nevertheless, the further continuity of this research was not considerably substantial. Challenges of the proposed optimum approach by Swaminathan and Voller [5] lie in computing the slope of the enthalpy-temperature curve dH dT for every iteration step. Pham [6] treated this by computing a time-averaged and space-averaged slope of the enthalpy-temperature curve dH dT . But, this is still limits the strong coupling of the enthalpy with temperature for a given phase change problem. For computations involving complex ge ometries and sharp phase change curves with thin melting ranges, such slope formulations dH dT may demand smaller time steps and conformal computational grid. Therefore, the present article introduces a smooth slope function for an iterative optimum approach to increase stability and reduce computational effort. Apart from the new adaptation of the optimum approach, this article also studies the influence of different simulation parameters like corrector steps, Darcy term constants and choice of smooth functions on numerical results. Suppression of velocity in the solid phase is also necessary to compute a solid-liquid phase change problem. Morgan [8] tried to implement this by overwriting the velocity in solid phase to zero. However, this can lead to instabilities in advanced applications. Variable viscosity method (VVM) [9–11] is also an alternative to model the ve locity suppression in solid phase of a melting problem. In the present
* Corresponding author. E-mail addresses: [email protected] (R.R. Kasibhatla), [email protected] (D. Brüggemann). https://doi.org/10.1016/j.ijthermalsci.2019.106187 Received 20 August 2019; Received in revised form 20 October 2019; Accepted 19 November 2019 1290-0729/© 2019 Published by Elsevier Masson SAS.
Please cite this article as: R.R. Kasibhatla, https://doi.org/10.1016/j.ijthermalsci.2019.106187
D.
Brüggemann,
International
Journal
of
Thermal
Sciences,
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
Volume of the cell, m3
Nomenclature
V
a b c capp d D Dc Ds erf g H L p q_ S sech T t tanh Δt u
Greek symbols α Volume fraction of the liquid β Thermal expansion coefficient, 1/K γ Porosity in the Kozeny-Carman equation ε Tolerance criterion for convergence κ Thermal conductivity, W/(mK) λ Enthalpy jump in a melting range, m2/s2 μ Dynamic viscosity of the fluid, kg/(ms) ρ Density, kg/m3 σ Square root of variance in the peak function, K ω Under relaxation factor
Discretized energy equation coefficient Discretized energy equation coefficient Mass specific heat capacity, J/(kg K) Mass specific apparent heat capacity, J/(kg K) Discretized energy equation coefficient Darcy source term, kg/(m3s) Morphological constant in Darcy term, kg/(m3s) Zero non-divisional constant in Darcy term Error function Gravitational acceleration, m/s2 Mass specific enthalpy, m2/s2 Latent heat, m2/s2 Pressure, kg/(m s2) Heat flux, W/(m2) Source term Hyperbolic secant function Temperature, K Time, s Hyperbolic tan function Time step, s Velocity, m/s
Subscripts l m n nb nþ 1 o p s
work, the most used porosity scheme with a Darcy source term [12] is employed to account the solid velocity suppression. In recent times, Faden et al. [13] implemented such a porosity approach to model the settling of solid phase in a macro-capsule with an implicit approach for settling. Whereas, Gudibande and Iyer [14] proposed another cell-splitting method to model a similar contact melting. Another nu merical approach was proposed by Kozak and Ziskind [15] to model the melting and settling of solid, which employs a peak function to restrict the solid bulk motion.
D¼
Physical phenomena in a thermo-fluid process can be represented by the conservative equations for mass, momentum and energy [16]. During heat transfer, small density change in the liquid phase influenced by a minute temperature difference can be assumed linear. In the present work, this small density change in a natural convective melting process is accounted by the Bousinessq approximation source term in the mo mentum conservation equation [17]. The equation set solved to repre sent the thermo-fluid process in a solid-liquid phase change is Continuity equation: Momentum equation: � ∂u ρ þ u ⋅ ru ¼ rp þ μΔu þ D ⋅ u þ ρg⋅ð1 ∂t
c ¼ α ⋅ cl þ ð1
�
Energy equation: � ∂H ρ þ r ⋅ ðuHÞ ¼ r⋅ðκrTÞ: ∂t
Tm ÞÞ
(4)
αÞ⋅cs :
(5)
The specific heat capacity c in this work contributes to the sensible enthalpy and the apparent heat capacity ceff is a change of total enthalpy with respect to the temperature change dH . In this work, it is represented dT as � � dH L T Tm ceff ¼ ¼ c þ ⋅sech2 : (6) dT 2σ σ
(1)
β ⋅ ðT
αÞ2 : α3 þ Ds
ð1
where, α is the volume fraction of liquid PCM, Dc is the Darcy constant and Ds is a small constant to avoid the division with zero. In Eq. (3), H is the enthalpy of a material and can be segregated into sensible and latent parts. Approaches like the source based method [3] use the segregated representation of enthalpy. However, the enthalpy can also be repre sented as a function of temperature H ¼ fðTÞ [21]. Theoretically, melting is an isothermal process for pure metals, whereas for binary alloys it occurs over a temperature range. In this work, melting over a very small temperature range is realized in the numerical simulations. The apparent heat capacity of the melting pro cess shows a large change over a small temperature range ΔTm [22]. In this article, the phase change is mathematically represented by a peak shape function. Here, the apparent heat capacity ceff is different from c in Eq. (5), which is
2. Numerical model
r⋅u ¼ 0
Dc ⋅
and superscripts Liquid Melting point Old iteration level Neighbouring cells Actual iteration level Old time step Node of the computing cell Solid
(2)
Here, σ is the square root of the variance of the peak shape function. Where, σ ¼ ΔT4m and ΔTm is the difference between solidus and liquidus temperature Ts and Tl respectively. As shown in Fig. 1, Eq. (6) is rep resented linear in the sensible parts and experiences a peak shape at the melting temperature range. The volume fraction of the melted portion α is formulated as a function of temperature. The volume fraction of the melted portion α is obtained from the integral function of apparent heat capacity written as Z T Z T dH dT ¼ dH ¼ ceff dT: (7) Tref dT Tref
�
(3)
Here, ρ is the density and μ the dynamic viscosity of the fluid, κ the thermal conductivity, u the velocity, p the pressure, β the thermal expansion coefficient of the PCM in liquid phase, D the Darcy source term, T the temperature and Tm the melting temperature. The Darcy source term in the momentum equation suppresses the velocity of solid phase in the computational domain, whereas the interfacial mushy zone between solid and liquid phase is treated as a porous medium [18,19]. This source term in Eq. (2) reads [12,20]. 2
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
This equation can be rewritten in an iterative form as [4]. � � �� � � � � � 1� dH n o n H H T nþ1 ¼ d H ap þ d p p p p þ S: dT n
(10)
The source of fluxes from neighbouring cells S is S¼
X anb T nþ1 nb þ b: In Eq. (10), ðHnp Þ
(11) 1
is the temperature which is the inverse function of
enthalpy at the iteration step n i.e. ðHnp Þ ðρVÞp . Δt
1
¼ T np and d is a discretization
coefficient with d ¼ In Eqs. (10) and (11), n þ 1 is the current iteration step, n is the old iteration step, o is the old time step, Δt is the time step size, nb represent the neighbouring cells of the computing cell p, ðVÞp is the volume of the computing cell p and a as well as b are the coefficients of the discretized energy equation. The numerical model in the present work employs a smooth slope function instead of a slope computation in the each iterative step. Here, instead of computing the slope of the H-T curve as [21]. � � dH Hn Hn 1 ¼ n : (12) dT n T Tn 1
Fig. 1. Apparent heat capacity ceff vs. temperature.
The present numerical model computes the slope of the H-T curve or the apparent heat capacity as � � � n � dH L T Tm ¼ c þ ⋅sech2 : (13) dT n 2σ σ
Which is shown in Fig. 2. By rearranging Eq. (7) and Eq. (6) in the enthalpy H, the liquid volume fraction α is obtained by � � T Tm α ¼ 0:5 þ 0:5⋅tanh : (8)
The new proposed step as in Eq. (13) reduces the numerical in stabilities caused by the slope computation as in Eq. (12). When the temperature change between the temperature at the old iterative step Tn and the pre old iterative step Tn 1 is zero, a division by zero occurs for Eq. (12). However, this can be treated by using a small non-zero additive value in the denominator of Eq. (12). But, this small non-zero additive can affect the maximum range of the apparent heat capacity alongside creating robustness and accuracy issues. Thereby, the proposed smooth function in this numerical model treats these robustness and accuracy issues by computing a definitive slope value for every temperature point as written in Eq. (13). In the numerical model, Eq. (10) with the slope values from Eq. (13) is iterated until the consistency of enthalpytemperature coupling is achieved in every time-step [4]. Inside an every time-step,
σ
As shown in Fig. 2, the increase in liquid volume fraction α occurs over a small melting temperature range ΔTm . 2.1. Enthalpy-temperature model The proposed optimum approach is a computationally inexpensive alternative to model the solid-liquid phase change [4]. This numerical model converges faster than the standard source based fictitious method and the apparent heat capacity method [4]. Within the approach, the discretized form of the energy equation Eq. (3) is described by i h X ap Tp ¼ (9) anb Tnb þ d H op Hp þ b:
1. Initially, the old values of enthalpy and temperature, Hop and T op respectively are updated to Hnp ¼ Hop and T np ¼ T op . 2. Eq. (10) is solved for the updated enthalpy and temperature values. By solving Eq. (10), the new temperature Tnþ1 for the n þ 1 iteration p and the enthalpy Hnp for the iteration level n is obtained. 3. In an iteration step, the explicit form of the enthalpy Hnp is updated � � using the new temperature Tnþ1 and the slope of the line dH in p dT n
equation H nþ1 ¼ H np þ p
� � h dH ⋅ T nþ1 dT n p
i T np :
(14)
Inside a time step, the steps 2 and 3 are repeated until convergence, i. e. the consistency of temperature with enthalpy. The convergence of the algorithm is characterised by a convergence tolerance ε: T nþ1 p
T np � ε:
(15)
The enthalpy update step in Eq. (14) is proposed by Swaminathan and Voller [4]. In the present article, an under relaxation factor w is used to relax the sharp enthalpy update from the new temperature: � � h i dH H nþ1 (16) ¼ H np þ w⋅ ⋅ T nþ1 T np : p p dT n
Fig. 2. Volume fraction of liquid, α vs. temperature of PCM over a phase change domain. 3
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
For optimum approach, Voller [5,21] proposed a temperature correction after the enthalpy update in Eq. (14). The correction step ensures the enthalpy relation with the temperature, which reads 8 H nþ1 p > ; H nþ1 < cs Ts > p > > cs > > > > " # > > > H nþ1 cs Ts > p > < Tm þ σtanh 1 cs Ts � H nþ1 p c T þ L cs Ts nþ1 l l Tp ¼ > > > > > H nþ1 � cl Tl þ L; > p > > > > > > > H nþ1 L : p ; H nþ1 > cl Tl þ L: p cl
The rectangular cavity had dimensions of 88.9 mm in height, 63.5 mm in width and 38.1 mm in depth. The copper plates arranged at the left and right boundaries ensure the implementation of a constant homogenous surface temperature. The top and bottom of the rectangular cavity were insulated [24]. An identical geometrical configuration is developed for the numerical setup. The material properties and boundary conditions used in the numerical set-up are summarized in Table 1. In numerical simulations, the nature of melting also depends on the resistance offered against velocity in solid phase, see Kasibhatla et al. [11]. Thereby, the Darcy term parameters similar to Brent et al. [3] have been implemented, where Dc ¼ 1:6⋅106 kg/m3s and Ds ¼ 10 3 . In nu merical simulations, for the 2D case front and back walls of the cavity are assumed symmetric and as walls for the 3D case. Grid dependence is studied to depict the influence of numerical model on the grid refine ment. A thin boundary layer of 2 mm is allocated to resolve the large heat fluxes. Apart, the rest of geometry is homogeneously refined into structured hexahedral cells. Cell refinement combinations of 87�62, 131�93 and 174�123 cells in X and Y direction respectively are studied. Between the difference in melt fronts obtained, the computational grid with 87�62 cells is employed in this work. For the 2D assumption only one cell is contributed in the Z direction. Whereas, for a 3D geometrical configuration a numerical grid with 441,870 cells contributing 87�62�55 cells in X, Y and Z directions is employed. Also, a wall boundary is imposed on the front and back boundaries of the rectangular cavity for the 3D case. Convergence criterion ε ¼ 10 6 is set for all the numerical simulations in this work.
(17)
The nature of enthalpy transport with the convective term is dependent on the solid-liquid phase change problem. For pure melting with no mushy zone, sensible enthalpy is modelled in the convective term [21]. This can be derived by assuming the continuity of liquid phase after segregating the total enthalpy into sensible and latent parts [23]. For an assumption, where the velocity of solid phase can be defined to zero. The convective term formulated with the apparent heat � � capacity r⋅ u dH T is derived into as r⋅ðucl TÞ [21]. In the present work, dT to mimic the pure melting of gallium, a very thin negligible mushy zone is assumed, alongside no movement of solid particles in the thin mushy zone. Therefore, in this work only sensible enthalpy is transported in the convective term of energy equation written as r⋅ðucTÞ.
4. Results and validation of the numerical model
3. Benchmark experiments and numerical set-up
A melting range of ΔTm ¼ 0.002 K is employed to mimic a pure isothermal phase change. The assumption of a 2D geometry results in computing high magnitude convective vortices. These vortices split up into two or more vortices leading to a deeply curved phase change boundary, also see Giangi and Stella [28]. This behaviour raised ques tions if the flow was mono-cellular or multi-cellular [26]. This influence of the velocity field may lead to an enhanced melting of the substance leading to such melting patterns. It could be noted from the literature that this occurs due to the influence of 2D geometrical configuration, see Nitrityuk [29]. Therefore, a 3D case with the front and back side of the cavity considered walls is computed [25]. As seen in Figs. 4 and 3D numerical results represent the melting pattern in reality very well and similar to the experimental findings of Gau and Viskanta [24]. The small deviation of numerical results from experiments can be observed in time t ¼ 10, 17 min in Fig. 4. The phase front from the 3D geometry experi ences a smooth continuous profile unlike the curvatures due to several vortices as seen in Fig. 4. The difference in nature of vortices in a 2D and 3D geometry is discussed by Nikrityuk [29]. The numerical error be tween the energy transported into the computational domain and energy used for phase change is investigated. This error made by the present optimized approach is approximately 0.01 %. In general, the 3D simulations are computationally very expensive than 2D simulations, which also increases the efforts to study phase change numerically. But with the present smoothed optimum approach,
Investigating the melting of 99.6 % pure gallium inside a rectangular cavity by Gau and Viskanta [24] are one of the initial experimental benchmarks in the field of solid-liquid phase change. Different numer ical models developed in the recent years [7,25] were also validated with these experimental results. Despite some strong criticism [26], the issue related to numerical convergence is resolved [27]. Validation of a numerical model with these benchmark results does not just provide an absolute accuracy of the numerical model. But, also provides a relativity of other numerical approaches validated against these numerical results. A simple graphical representation of the experimental set-up is shown in Fig. 3. The rectangular cavity in the experimental set-up is filled with a very pure gallium. The convective melting of gallium inside in the cavity is measured at different instances of time.
Table 1 Parametric tabular of the numerical set-up. Material properties Density (ρ) Viscosity (μl )
4
kg/m Pa⋅ s
values 3
6093 0.00181
Specific heat capacity (c) Latent heat (L) Thermal conductivity (κ) Melting temperature (Tm )
J/kgK J/kg W/mK K K
311.15
Cold boundary temperature (Tc )
K
301.45
Heated boundary temperature (Th )
Fig. 3. Schematic diagram depicting the experimental set-up of rectangular cavity from Gau and Viskanta [24].
units
381.5 80160 32.0 302.93
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
5.1. Temperature corrector-step The numerical results with and without the temperature correction after the enthalpy update in the algorithm is studied and the comparison are shown in Fig. 6. Despite, the similar melt fronts in Fig. 6, for thin melting ranges the numerical model without a temperature corrector step violates the energy equation to a large extent. The numerical error denoted as the difference between heat trans ferred into the computational domain and the enthalpy change is compared for uncorrected and corrected models in Table 2. It can be noticed in Table 2 that the error made by the numerical model without a corrector step gradually increases with a decrease in the melting range. For thin pure melting ranges, the error made is very large. On the other hand, for melting ranges greater than 1 K the error made by the nu merical model without a temperature corrector step is relatively negli gible. This can also be noticed in the volume fraction of liquid PCM α estimated by the numerical model as in Table 3. Various factors like the size of finite volume, temperature gradient and the expansion coefficient also play a role in the numerical error. Therefore, a generalization of such a melting range value cannot be made to assert where the tem perature corrector step really is influential. To analyse this, different parameters need to studied to understand the sensitivity of temperature corrector step.
Fig. 4. Validation of the numerical model using a 3D numerical grid with the benchmark experiments [24].
the numerical effort to simulate such 3D cases including computing phase change in complex geometries is treated well. Different commercial software employ the standard source based method, due to its stability [30]. But, optimum approach with the stable slope function dH dT is a decent alternative to reduce computational effort alongside providing stability. The computational expense of the present smoothed optimum approach is compared with the fictitious source based method [3,11]. In each time step the iterations required by the present model with a temperature correction step in comparison with the source based method is shown in Fig. 5.
5.2. Choice of peak function In this section, two different slope functions are compared to un derstand the influence of numerical model on its dH dT slope function. However, at convergence the product of apparent heat capacity and � � temperature difference dH and ½Tnþ1 ðHnp Þ 1 � becomes smaller and p dT n
non-influential. This is proved by varying the choice of apparent heat capacity. For this, a Gaussian normal distribution function is employed to represent the peak behaviour of the apparent heat capacity as shown in Fig. 1. For this comparison, the value of σ ¼ ΔT4 m is set same for both the functions. Mathematically, the total integral of different peak functions remain same independent of the variance thickness σ . But, in numerical simulations a higher value in the denominator ensures a decent compression of the solid-liquid interface. Whereas, a lower value in the denominator results in a broader solid-liquid interface. Accordingly, the Gaussian normal distribution function is written in the form
5. Investigating the parametric influence Different parameters involved in the numerical model are investi gated to understand their influence on the simulation results. Numerical model parameters like the efficiency of temperature corrector step from Eq. (17), Darcy term constants and choice of the peak function are studied in this case. Due to a high computational effort by the 3D case, a 2D case with 87�62 cells is employed to study the parametric influence of the numerical results.
Fig. 6. Comparison of the melting fronts at time t ¼ 2 min and t ¼ 17 min with and without temperature correction after enthalpy update for different melting ranges.
Fig. 5. Iterations required by the optimum approach and source based method to converge shown on a logarithmic scale Eq. (15). 5
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
pressure drop. However, the empirical Kozeny-Carman equation is
Table 2 Numerical error for uncorrected and corrected temperature steps. ΔTm ðKÞ
Error T-uncorrected %
Error T-corrected %
0.01 0.1 0.25 1 2
15 2 0.26 0.191 0.159
0.0002 0.00015 0.002 0.17 0.2
rp ¼
α with T-uncorrected
α with T-corrected
0.01 0.1 0.25 1 2
0.8367 0.6817 0.6761 0.6620 0.6506
0.6704 0.6735 0.6717 0.6639 0.6501
� ceff ¼
dH ¼cþe dT
�2 ðT Tm Þ 2⋅σ
L ⋅ pffiffiffiffiffi: σ 2π
(20)
where γ is the porosity of bed, μ is the viscosity of fluid entering the packed bed, ϕ is the sphericity of packed bed particles and Db their diameter of the particles in bed. Accordingly, this equation is rewritten as Eq. (4) in the solid-liquid phase change problems. In Eq. (4), Dc is the morphological Darcy constant and Ds is defined a small value restricting the division by zero [3]. But, the selection of these values is still inconclusive. Therefore, in this article the constants Dc and Ds are varied such that the total resistance value in the solid phase remains same as in Brent et al. [3]. By substituting, selected values as Dc ¼ 1:6⋅106 kg/m3s and Ds ¼ 10 3 , the pressure gradient in the solid phase becomes 1:6⋅109 times the velocity and 0 in the liquid phase. But, this value of pressure gradient can be achieved by different choice of Dc and Ds values. These combinations of constants Dc and Ds are presented to achieve these absolute values in both the phases as shown in Fig. 8. Hereby, the pressure gradient in the mushy zone vary between them. For this, all the other simulation parameters are maintained constant as in Table 1 except the Darcy constants Dc and Ds . Fig. 9 explains the importance of total resistance offered by the Darcy term i.e. the product of Dc and D1s . Various literatures in the field of solid-liquid phase change
Table 3 Volume fraction of liquid phase for the uncorrected and corrected temperature steps at time t ¼ 20 min. ΔTm ðKÞ
150μ ð1 γÞ2 ⋅u: γ3 ϕ2 D2b
(18)
report Ds as a constant restricting the division by zero. Additionally, another numerical study is executed to depict the in fluence of total resistance on the melting patterns. For this, the numer ator Dc has been varied from Dc ¼ 1:6⋅106 to Dc ¼ 1:6⋅103 and Dc ¼ 1:6⋅1010 . Here in this study, the small non-division constant Ds is left constant at 10 3 . The results are shown in Fig. 10. The results of numerical benchmark as seen in Figs. 9 and 10 sup ports the argument that total resistance offered in the solid-phase of momentum equation is important to represent a melting pattern rather than the independent constants Dc or Ds . But, the same can play a pivotal role in modelling binary alloys and impure materials. Thereby, a need to study the nature of resistances offered in the fluid-solid combined mushy zone exists.
Accordingly, the integral of this function i.e. the liquid volume fraction is represented by an error function written as � � T Tm pffiffi : (19) α ¼ 0:5 þ 0:5⋅erf σ 2 Thereby, numerical results between Gaussian normal distribution and secant function formulation are compared in Fig. 7. As seen in Fig. 7, difference between the numerical results is very small and can be neglected. Therefore, the selection of apparent heat capacity has rather smaller influence on the solution accuracy but only on the stability. 5.3. Darcy term Darcy term implemented to suppress the velocity in solid phase is inspired from the Kozeny-Carman equation used to model the pressure drop of fluid in the packed beds [12,20]. In solid-liquid phase change, the Kozeny-Carman equation assists the suppression of velocity advec tion in solid phase. This resistance offered is implemented in the mo mentum equation with constants like Dc and Ds enabling a typical
6. Discussion and conclusion
Fig. 7. Comparison between different numerical models for the choice of the apparent heat capacity approximate.
Fig. 8. Plot showing different progressions of Darcy term with different values of Darcy constants Dc and Ds .
An optimized method with a smooth continuous phase change curve is presented to compute solid-liquid phase change. A small error with very low computational time is the topic of interest here. The numerical results employing the 3D geometry have shown a very good agreement with the experiments. Whereas, the 2D geometry case with the
6
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx
specific heat capacity for phase change over temperature range. For this, as seen in Fig. 7 the choice of peak function has shown no influence on the numerical results. This is because the optimum algorithm employs the specific heat capacity function to achieve the enthalpy change till the temperature difference between the old and new iteration step becomes small, as described in Eq. (14) and Eq. (15). Another numerical parameter namely the porosity resistance to model solid-liquid phase change has been investigated in this work. Here, initially for the choice of constants Dc and Ds are chosen in such a way to achieve the same momentum resistance in the solid-phase. Re sults for numerical simulations with these different choice of constants are compared in Fig. 9. The results show a similar melt front for different choice of Dc and Ds values when the total resistance offered in solid phase is the same. This argument is further supported by an additional numerical study depicted in Fig. 10, where the value of Darcy constant Dc is varied by keeping the small non-division constant Ds the same. Different melt fronts shown in Fig. 10 depict the influence of total resistance on melting fronts. Hereby, when modelling the solid-liquid phase change of binary alloys or impure materials the particle distri bution in the interface region needs to be understood. For the same solid phase resistance as seen in Fig. 8 the source term curves takes different values for various mixture volume fractions. The amount of resistance necessary to restrict the motion of solid particle between 0 and 100% volume fraction moving in the molten liquid at the interface is still unknown. Because, for different combinations of Dc and Ds the mo mentum resistance in solid-liquid mixture at the interface takes distant values. Therefore, a detailed study on the choice of these values is necessary to study the dendrite growth for binary alloys in moving boundary problems.
Fig. 9. Comparison of numerical results with different Dc and Ds combinations but same total resistance in the full solid phase.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.ijthermalsci.2019.106187. References [1] Z. Liu, Y. Yao, H. Wu, “Numerical modeling for solid–liquid phase change phenomena in porous media: shell-and-tube type latent heat thermal energy storage, Appl. Energy 112 (2013) 1222–1232. [2] S. Chakraborty, D. Chatterjee, “An enthalpy-based hybrid lattice-Boltzmann method for modelling solid–liquid phase transition in the presence of convective transport, J. Fluid Mech. 592 (2007) 155–175. [3] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal, Numer. Heat Tran. 13 (3) (1988) 297–318. [4] C. Swaminathan, V. Voller, A general enthalpy method for modeling solidification processes, Metall. Trans. B 23 (5) (1992) 651–664. [5] C. Swaminathan, V. Voller, On the enthalpy method, Int. J. Numer. Methods Heat Fluid Flow 3 (3) (1993) 233–244. [6] Q. Pham, Comparison of general-purpose finite-element methods for the stefan problem, Numer. Heat Tran. 27 (4) (1995) 417–435. [7] A. K€ onig-Haagen, E. Franquet, E. Pernot, D. Brüggemann, A comprehensive benchmark of fixed-grid methods for the modeling of melting, Int. J. Therm. Sci. 118 (2017) 69–103. [8] K. Morgan, A numerical analysis of freezing and melting with convection, Comput. Methods Appl. Mech. Eng. 28 (3) (1981) 275–284. [9] Y. Asako, M. Faghri, M. Charmchi, P.A. Bahrami, Numerical solution for melting of unfixed rectangular phase-change material under low-gravity environment, Numer. Heat Tran., Part A: Applications 25 (2) (1994) 191–208. [10] D.K. Gartling, Finite element analysis of convective heat transfer problems with change of phase, Comput. Method Fluids (1980) 257–284. [11] R.R. Kasibhatla, A. K€ onig-Haagen, F. R€ osler, D. Brüggemann, Numerical modelling of melting and settling of an encapsulated pcm using variable viscosity, Heat Mass Tran. (2016) 1–10. [12] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, Int. J. Heat Mass Tran. 30 (8) (1987) 1709–1719. [13] M. Faden, A. K€ onig-Haagen, S. H€ ohlein, D. Brüggemann, An implicit algorithm for melting and settling of phase change material inside macrocapsules, Int. J. Heat Mass Tran. 117 (2018) 757–767. [14] N. Gudibande, K. Iyer, Numerical simulation of contact melting using the cellsplitting modified enthalpy method, Numer. Heat Tran., Part B: Fundamentals 71 (1) (2017) 84–107.
Fig. 10. Depicting the numerical results obtained by varying Darcy constants Dc and letting the non-division constant Ds in the denominator remain constant at Ds ¼ 0.001.
symmetric boundary conditions on front and back walls is employed to investigate the numerical model parameters in this work. Contrast to the source based method, the present numerical model is a fast converging iterative approach which needs five to ten iterations in each time step to converge. This benefit assists with effortless solving of phase change problems for dense computational grids. Alongside, different simulation parameters were also studied in the present work. Here, the temperature correction after the enthalpyupdate ensures energy conservation in every time-step. Here, the very small numerical error in the energy equation results in the distorted phase front in the numerical simulations. This numerical error is large for very thin melting ranges rather than the broad melting ranges as seen in Table 2 and Table 3. As seen in Fig. 6 for a similar melting range ΔTm ¼ 1 K, the numerical results with and without temperature correc tion show a very different phase fronts. Here, it can also be seen that the numerical error increases with the decrease in melting range i.e. tem perature re-correction after an enthalpy update is very necessary when modelling pure melting problems. Whereas, for binary alloy melting problems the same can be tolerated with a small numerical error, if a very low computational expense need to be achieved. Furthermore, different mathematical functions were also assumed to model the 7
R.R. Kasibhatla and D. Brüggemann
International Journal of Thermal Sciences xxx (xxxx) xxx [22] M. Reading, D. Elliott, V. Hill, A new approach to the calorimetric investigation of physical and chemical transitions, J. Therm. Anal. Calorim. 40 (3) (1993) 949–955. [23] J. Ni, C. Beckermann, A volume-averaged two-phase model for transport phenomena during solidification, Metall. Trans. B 22 (3) (1991) 349. [24] Gau and Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Tran. 108 (1) (1986) 174–181. [25] O. Ben-David, A. Levy, B. Mikhailovich, A. Azulay, 3d numerical and experimental study of gallium melting in a rectangular container, Int. J. Heat Mass Tran. 67 (2013) 260–271. [26] J.A. Dantzig, “Modelling liquid–solid phase changes with melt convection, Int. J. Numer. Methods Eng. 28 (8) (1989) 1769–1785. [27] N. Hannoun, V. Alexiades, T.Z. Mai, Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side, Numer. Heat Tran. Part B: Fundamentals 44 (3) (2003) 253–276. [28] M. Giangi, Fulvio Stella, “Melting of a pure metal on a vertical wall: numerical simulation, Numer. Heat Tran. Part A: Applications 38 (2) (2000) 193–208. [29] P.A. Nikrityuk, Computational Thermo-Fluid Dynamics: in Materials Science and Engineering, John Wiley & Sons, 2011. [30] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, Int. J. Heat Mass Tran. 50 (9–10) (2007) 1790–1804.
[15] Y. Kozak, G. Ziskind, Novel enthalpy method for modeling of pcm melting accompanied by sinking of the solid phase, Int. J. Heat Mass Tran. 112 (2017) 568–586. [16] S.V. Patankar, Numerical Heat Transfer and Fluid Flow”, Hemisphere Publishing Corp., Washington, DC, 1980. [17] J. Boussinesq, Th�eorie de l’�ecoulement tourbillonnant et tumultueux des liquides dans les lits rectilignes a grande section, Gauthier-Villars et Fils, Paris, 1897. [18] V. Voller, N. Markatos, M. Cross, Solidification in convection-diffusion, in: N. Markatos, M. Cross, D. Tatchell, N. Rhodes (Eds.), Numerical Simulation of Fluid Flow and Heat/Mass Transfer Processes, vol. 18, Springer Berlin Heidelberg, 1986, pp. 425–432, of Lecture Notes in Engineering. [19] F. R€ osler, Modellierung und Simulation der Phasenwechselvorg€ ange in makroverkapselten latenten thermischen Speichern”. PhD thesis, University of Bayreuth, Logos Verlag Berlin, 2014. [20] P. Xu, B. Yu, “Developing a new form of permeability and kozeny–carman constant for homogeneous porous media by means of fractal geometry, Adv. Water Resour. 31 (1) (2008) 74–81. [21] W. Minkowycz, E. Sparrow, J. Murthy, Handbook of Numerical Heat Transfer, Wiley, 2006.
8