- Email: [email protected]
An enhanced numerical approach for convection phase change problems: A solution of tin melting problem
Journal Pre-proof An enhanced approach for convection phase change problems: A solution of tin melting problem Salaheddine Kaba, Khalid Achoubir, Abdelkhalek Cheddadi PII:
S2214-157X(19)30429-0
DOI:
https://doi.org/10.1016/j.csite.2020.100585
Reference:
CSITE 100585
To appear in:
Case Studies in Thermal Engineering
Received Date: 26 October 2019 Revised Date:
15 December 2019
Accepted Date: 4 January 2020
Please cite this article as: S. Kaba, K. Achoubir, A. Cheddadi, An enhanced approach for convection phase change problems: A solution of tin melting problem, Case Studies in Thermal Engineering (2020), doi: https://doi.org/10.1016/j.csite.2020.100585. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
AN ENHANCED APPROACH FOR CONVECTION PHASE CHANGE PROBLEMS: A SOLUTION OF TIN MELTING PROBLEM Salaheddine Kaba1, Khalid Achoubir1 and Abdelkhalek Cheddadi1,* 1
Mohammed V University in Rabat, Mohammadia School of Engineers, Ibn Sina Str., P.O. Box 765 Agdal, Rabat, Morocco * Corresponding author, Email address : [email protected] Abstract
A comprehensive and efficient numerical phase change model compatible with the stream function-vorticity formulation is developed. The model is based on the source term of the enthalpy-porosity method derived from Darcy's law. The results of the model are compared with those of tin melting at = 0.02, = 2.5 × 10 , = 0.01, in this case the method demonstrates its capability to deal with this controversial case by using an optimal grid. The comparison with results obtained by using this approach was satisfactory. Keywords: melting, phase change materials (PCM), stream function-vorticity, enthalpy-porosity.
Nomenclature Alphabetically listed terms Darcy source term Constant of Darcy source term Thermal capacity (J/kg K) Liquid fraction g Gravity acceleration (m/s²) Thermal conductivity (W/m K) Latent heat (J/kg) Reference length (m) Pr Prandtl number Rayleigh number Ste Stefan number S Source term t Dimensionless time T Dimensionless temperature u Horizontal velocity v Vertical velocity x,y Cartesian coordinates
Greek symbols Thermal diffusivity (m2/s) Thermal expansion coefficient (K-1) υ Kinematic viscosity (m2/s) Vorticity Density (kg/m3) Stream function ∆" Temperature difference between walls δ" Dimensionless half width of mushy zone temperature range Subscripts and superscripts 0 Reference value # Liquidus $ Solidus % Hot Cold ′ Dimensional value
1. Introduction The numerical modeling of phase change remains since several decades an open field of investigation, due to the presence of melting and solidification techniques in diverse technologic domains such as: the fabrication of metallic alloys of high purity, thermal energy storage in term of latent heat, etc.
2
The major challenge of the mathematical modeling of phase change is the moving of the solid-liquid interface. Indeed, the phase change is characterized by a mobile front accompanied by the absorption or desorption of the latent heat with thermophysical properties variation, including the thermal conductivity, the specific heat, the viscosity and the density. Furthermore, most of phase change applications take place in a multidimensional configuration. In addition, the natural convection in liquid phase plays an important role in the heat transfer control which further complicates the problem. Usually, there are two groups of numerical model for studying this kind of problem: the fixed grid schemes and moving grid schemes [1]. The last group uses a transformed system of coordinates and sophisticated mathematics for the interface tracking, which may cause a lot of problems in programming and ineluctably results in a higher computational cost. On the other hand, in the fixed grid schemes, the position of the interface is deduced implicitly from the enthalpy/temperature distribution, leading to relatively simple mathematical models. Three principal methods were proposed for the formulation of the conservation equations of energy: the enthalpy method [2], the equivalent thermal capacity method [3], and the temperature transforming model [4]. Before the application of each of these models, one must express the velocity in the two phases, solid and liquid. The velocity must be of course zero in the solid phase and calculated from the equations of motion in the liquid phase and in the mushy zone in the case of non isothermal phase change. Several methods were proposed for the treatment of this problem. We cite the variable viscosity method [5], the source term method [6-11] and the switch-off method [12]. All of these methods were developed for the equations of motion in terms of velocity-pressure formulation. The using of stream function-vorticity formulation in the phase change problems (melting or solidification) is limited. Giangi et al. [13] used this formulation to study the phase change of a pure element. For stopping the velocity they used a source term without presenting explicitly the method used to derive it. Timchenko et al. [14] used this formulation to simulate the solidification of an alloy at eutectic composition under microgravity conditions. Those authors [13,14] have used an adequate approach for the cases where the interface is very thin and boxed in one calculation cell; the stream function , the vorticity and the velocities were considered nil in the solid phase and calculated from the definition of the stream function in the liquid phase. It should be noted that these works [13,14], didn’t present a method for the treatment of non isothermal phase change, where the existence of mushy zone is indispensable and request an additional relation to identify it. The objective of the present work is to apply an approach for the correction of velocity, which will be adequate with the stream function-vorticity formulation. This method is based on the frequently used Darcy source term approach, especially in the numerical studies of the phase change materials (PCM) dedicated to the thermal energy storage [15,20]. The capability of this approach is demonstrated in the comparison with the controversial case of tin melting. 2. Mathematical model and approximation 2.1 Numerical approximation The fluid considered is supposed Newtonian and incompressible. The flow is laminar and two dimensional with gravity parallel to the ascendant y axis. The thermophysical properties are considered constant except the density which is approached in the gravity term with the Boussinesq hypothesis: = ' (1 − (" ′ − "'′ ), (1) ′ where "' is the reference temperature, ' is the reference density, is the thermal expansion coefficient.
3
The dimensionless variables are defined by: ′ -′ .′ (-, .) = / , 0, = 1 5′
(2, 3) = 4
6
,
7′ 6
8, " =
9 ′ :9;′ ∆9 ′
is the reference length, is taken as
6
>?@A
.
(2)
is the thermal diffusivity, ∆" ′ = "<′ − "=′ . The reference velocity
2.2. Treatment of the equation of motion The stream function-vorticity formulation ( , ) proposed in this paper is derived from the equations of motion (in dimensionless form) written in primitive variables that contain the source term as proposed in [2]: B5 BC B7 BC
B5
B5
B7
B7
+ 2 BE + 3 BF =
+2
BE
+ 3
Where L
=6
BF
=
BH
G 1 2 − BE + 2
G13 −
BH BF
+
(3) (" − "' ) + 3
is a parameter used to correct the velocity,
is the Prandtl number.
(4) =
IJ∆9 ′ >K?@A L6
is the Rayleigh number and
In the context of this work, the final shape of the parameter depends on the approximation used for the flow in the mushy zone (where liquid and solid phases are mixed). In [2, 21], the expression of is obtained by considering the cavity as a porous medium. However in the case of isothermal phase changes, the same model can be used with a dummy mushy zone with an infinitesimal thickness and since the depends on the liquid fraction , a suitable relationship between expression of the source term temperature and liquid fraction must be established as will be seen later to treat this type of transition. Voller and Prakash [2] proposed the following expression which claims to prevent fluid flow in porous media: =−
=(M: N )O ( NK PQ)
(5)
Where is an arbitrary constant (large), and R is also a constant (very small), introduced to avoid the possible division by zero. When the liquid fraction is zero = 0, the control volume is totally solid, therefore the source term tends to the infinity forcing the velocity to switch-off. On the other side when = 1, the source term vanishes and the velocity is calculated through the equations of motion in fluid medium. The expression (5) is used primarily to treat non-isothermal phase changes (wide mushy zone). According to [2], this expression forces the fluid to follow the Carman-Kozny law, which is experimentally validated for flow in dendritic regions. One of the numerical disadvantages of this formulation is the high computational cost (presence of two powers with a division), hence the use of a less complex expression, Eq. (6), is required as was proposed by [11]. = (1 − ) (6) Note that when = 0, one obtains = . If the constant is not well chosen, it is possible that the flow will not be sufficiently turned off. In the present work, we first inject Eq. (6) in Eqs. (3) and (4). By differentiating Eq. (3) with respect to y, Eq. (4) with respect to x and subtracting one from the other (3) - (4). we obtain the stream functionvorticity formulation( , ) with an adequate source term T as follows: BT BC
BT
BT
+ 2 BE + 3 BF =
G1 −
B9 BE
+
T
(7)
4
= − G 1
2=−
BU
With:
T
BF
; 3 =
BU BE
(8)
(9)
= [(1 − ) + 3
B N BE
−2
B N ] BF
(10)
2.3. Treatment of energy equation: According to [22], the dimensional energy conservation equation can be expressed directly in terms of temperature. B9 ′ BC ′
+
(2′
B9 ′ BE ′
+ 3′
B9 ′
BF ′
) = G 1 " ′ + ′ 9
(11)
Where ′ 9 is the source term of energy equation: ′9 = −
> B N
=YN BC ′
−
>
=YN
B
B
(2′ BEN′ + 3 ′ BFN′ )
(12)
Introducing the dimensionless temperature defined in Eq.(2) leads to the following conservation equation in dimensionless form: B9 BC
B9
B9
M B N BC
+ 2 BE + 3 BF = G 1 " − ZC
Where
=
=Y ∆9 ′ >
M
− ZC (2
B N BE
+ 3
B N ) BF
(13)
is the Stefan number.
The convective part of the source term 9 (second term in the right member of Eq. 12) will be neglected in the case of isothermal phase change (e.g melting of pure elements), and will be introduced in the case of non-isothermal phase changes (e.g alloys solidification) as operated in [21]. 2.4. Method of numerical resolution In the present work the equations of motion are discretized with a centered finite difference method, with an adequate implementation of source term given by Eq. (10) in each compute node. On the other side the energy conservation Eq. (13) is discretized with a finite volume method by adopting a staggered grid for " compared to that used for and (Figure 1). The convective terms of Eq. (7,8) are discretized according to a centered scheme while those of Eq. (13) are discretized according to an Upwind scheme. The temporal part of source term of Eq. (13) is treated according to an Euler explicit scheme while the convective part is dicretized based on the method proposed by Patankar [23]. For the resolution we adopt the ADI method (Alterning Direction Implicit) with internal iteration, which generates tridiagonal systems of linear equations. These systems of equations are solved with Thomas algorithm (TDMA). Thus the resolution is performed with a considerable reduction in calculation time; Figure 2 shows a general algorithm of the major steps of the used code. The convergence is established at each time step when the maximum relative deviations are less than -5 10 for the dynamic and thermal fields. The dimensionless time step was between (10-5 - 10-6), the mesh was between (61×61 and 101×101). The refinement of the time and space step is carried out at the beginning of the convective melting regime and at the zone of transition (cell merging). The grid used is given for each simulation (See figure captions)
5
Figure 1. Presentation of mesh with staggered grid. Iterative updating of liquid fraction: To ensure the convergence of Eq. (13), special treatment of the source term is performed at each time step to adapt the liquid fraction in each cell with temperature calculated through the equation of energy. An important property which strongly affects the convergence of the internal iteration process is the procedure used for updating the liquid fraction from one iteration to another. Based on the method presented in [22], using an algorithm considered adequate for the update of the liquid fraction in the case of binary alloy solidification, we propose the following method. The temperature-liquid fraction relationship may be found by the experiment, or drawn from one of the metallurgical relationships available in literature, e.g. Scheil equation [24]. In the present work, the liquid fraction is taken under a linear form as a function of temperature: =
9 [ : 9\[ 9][ :9\[
(14)
In the system which is undergoing a phase change, the total enthalpy may be written as follows: % = "′ + (15) Where at the nth iteration we get: ^
%^ = " ′ + ^ (16) At each time step, and during the updating process by internal iterations, the total enthalpy of each node Eq. (15), is conserved between iterations (n) and (n+1). Only the distribution of the latent and the sensible enthalpy (respectively and " ′ ) will be different. Taking guidance from Eq. (14) and (16), we obtain: ^
"′ + With:
^
=
"′
^PM
^PM
+
^PM
"′ = ^PM (" ′ − " ′ _ ) + " ′ _ From which one gets: ^PM
=
` =Y 49 ′ :9 ′ a 8P N` >
>P1b9 ′ =Y
(17) (18)
(19)
Where c" ′ is a numerical value representing the half width of mushy zone temperature range. In the case of melting/solidification of a mixture of multiconstituents, this parameter is defined from the composition of each constituent. In the case of pure substances, a fictive value of c" ′ is introduced, which does not represent the physical reality, although it is used as a parameter for bypassing the discontinuity associated with the transition from one phase to another and avoid the numerical complexity associated with this transition. Practically c" ′ is selected so as to have " − "_ ∿ 0.
6
Finally, Eq. (19) can be written in the dimensionless form as follows: ^PM
=
ZC (9 ` :b9)P N` MP1 b9 ZC
Where c" is the dimensionless half width of mushy zone temperature range.
(20)
3. Results and discussion The simulations are undertaken in a square cavity (Figure 2) where the reference length is taken equal to the side of cavity. The Tin melting in a square cavity differentially heated, Bertrand et al. [25], is characterized by the presence of instabilities at the start of the convective regime. These instabilities have been identified in the work of Dantzig on the Gallium melting [26], who links this behavior to the instabilities encountered at the end of the conductive regime.
Figure 2. Geometric configuration of melting problem. In the case of melting, the experiment starts when one of the walls is heated to the melting point of Tin, then the heat is conveyed through a liquid layer [21]. Therefore the physical properties of the liquid will govern the melting process. The Tin physical properties are shown in Table 1. Table 1. Tin thermal parameters for melting problem Parameter Cp K ρ υ L TH TC Lref g Ra Pr Ste
Value 200.0 60.0 7.5 × 103 4.0 × 10-5 8.0 × 10-7 2.67 × 10-4 6 × 104 508.0 505.0 0.1 10 2.5 × 105 0.02 0.01
Unit J/kg K W/m K Kg/m3 m2/s m2/s K-1 J/kg K K m m/s2
The thermal parameters used as cited in the benchmark of Bertrand et al. [25] are: = 0.02, = 2.5 × 10 , = 0.01. The grid used in our study varies between 61×61 for the majority of cases and 101×101 for the detection of instabilities at the start of the convective regime.
7
At t = 1, Figure 3-a shows a good concordance regarding the shape and the location of the interface of fusion, between the participants in the benchmark on the one hand and the present work on the other. Meanwhile, a clear discrepancy between the participants is noticed at evolved time, especially at t = 4 (Figure 3-b). Some of them found an interface with two bumps which correspond to two co-rotative convective cells, while the others obtained an interface with one bump resulting from one single convective cell. The obtained results within the framework of the present paper are in agreement with the first group of participants (Figure 3-b). The interface shape and the convective cells number (Figure 4) confirms this observation. This contradiction was discussed by Gobin and Lequere [27] in a synthesis that followed this comparison exercise. The authors concluded that the simulation at low number tends to confirm the result obtained by Dantzig [26]. A more recent numerical investigation [28] was dedicated to resolve the controversy. It showed that by refining the grid (from 100×100 up to 400×400) and by using various numerical schemes (Figure 4), it becomes clear that the solution giving two cells is numerically more stable and that the one cell regime appears only when a coarse grid is used. These comparisons show the capacity of the approach used in this work to detect the instabilities that are detectable only at very fine grids in previous works [28, 29] (Figure 4). These instabilities, manifested by the shape of the interface and the number of convective cells which appear at successive times, are solved in the present work with a relatively coarse grid compared with the one used in [28]. We note as well from Figure 4 that the number of co-rotating cells is high at the beginning of the convective regime, and starts to reduce (recombining of convective vortices) to give rise to one cell dominated structure toward t = 8 (Figures 4 and 3-c). In Figure 5, we represent the evolution of the flow intensity efE as a function of time, showing the flow structure at each transition. According to the shape of the curve representing efE evolution, we note that the transitions, which correspond to a merging of convective cells, are associated with more or less important fluctuation in the value of efE . Table 2 presents the transition times and the resulting regime. Table 2. Regimes and instants of cells merging. Grid 100×100 Time
Flow regime
0.3
Beginning of 4 cells regime
0.6-0.8
Merging : 4 cells into 3 cells
1.8-2
Merging : 3 cells into 2 cells
4.5-5.5
Merging : 2 cells into 1 cell
8
(a) t = 1
(b) t = 4
(c) t = 8
Figure 3. Evolution of melting interface: (a,b) Comparison with Bertrand et al. [25] c) Comparison with Hannoun et al. [28]. = 0.02, = 2.5 × 10 , = 0.01. Grid: 61×61.
(a) t = 0.8
(b) t = 4
(c) t = 8
Figure 4. Stream lines and melting interface (Part of cavity) at different times. Pr = 0.02 , Ra = 2.5 × 10 , Ste = 0.01. Hannoun et al [28] (left) Grid: (100×100 up to 400×400). Present work (right) Grid: 101×101 (t = 0.8). Grid: 61×61 (t = 4 and t = 8).
Figure 5. Evolution of maximum value of stream function efE as function of time. = 2.5 × 10 , = 0.01. Grid: 101×101
= 0.02 ,
9
4. Conclusion In this work, a numerical method for the correction of the velocity adequate with the stream functionvorticity formulation was adopted, based on the source term approach. This method is valid for modeling both the isothermal and non-isothermal phase changes. The comparison with the case of the melting (low Prandlt number) is investigated, = 0.01. The results were in agreement with the reference solutions available in literature [28,29]. According to the results of this work, we can conclude that the coupling of the enthalpy porosity techniques with the stream function-vorticity formulation gives an acceptable solution for the controversial case of Tin melting problem with a relatively simple implementation of the code and a very optimal computational time.
References [1] A. A. Samarskii, P. N. Vabischchevich, O. P. Iliev, A. G. Churbanov.”Numerical Simulation of Convection/Diffusion Problems: A Review”. Int. J. Heat Mass Transfer 36 , no. 17 (1993): 40954106, [2] V. R. Voller, C. Prakash.”A fixed grid numerical modeling methodology for convection– diffusion mushy region phase-change problems”. Int. J. Heat Mass Transfer 30, no. 12 (1987): 1079–1718. [3] K. Morgan.” A numerical analysis of freezing and melting from a vertical wall”. Comput. Methods Appl. Mech. Eng 28, no. 3 (1981): 275–284,. [4] Y. Cao, A. Faghri.”A numerical analysis of phase change problem including natural convection”. ASME J. Heat Transfer 112, no. 3 (1990): 812–815. [5] D. K. Gartling. Finite element analysis of convective heat transfer problems with change of phase. Computer Methods in Fluids. Pentech, 1980. [6] B. Binet, L. Lacroix. ”Melting from heat sources flush mounted on a conducting vertical wall”. Int. J. Num. Meth. Heat & Fluid Flow 10, no. 3 (2000): 286-306. [7] V. R. Voller, M. Cross, N. C. Markatos. ”An enthalpy method for convection/diffusion phase change”. Int. J. Num. Meth. Eng 24, no. 1 (1987): 271-84. [8] P. A. Galione, O. Lehmkuhl, J. Rigola, A. Oliva. “Fixed grid modeling of solid-liquid phase change in unstructured meshes using explicit times schemes”. Numerical heat transfer part B: Fundamentals 65, no. 1 (2014): 27-52. [9] M. Ghalambaz, A. Doostanidezfuli, H. Zargartalebi, A. J. Chamkha. “MHD phase change heat transfer in an inclined enclosure: effect of a magnetic field and a cavity inclination”. Numerical heat transfer part A: Applications 71, no. 1 (2017): 91-109. [10] Y. Belhamadia, A. S. Kane, A. Fortin. “An enhanced mathematical model for phase change problems with natural convection”. Int. J. of numerical analysis and modeling, series B 3, no. 3 (2012): 192-206. [11] A. A. Al-abidi, S. B. Mat, K. Sopian, M. Y. Sulaiman, A. Th. Mohammed. “CFD applications for latent heat thermal energy storage: a review”. Renewable and Sustainable Energy Reviews 20, (2013): 353-363. [12] M. Yang, W. Q. Tao. ”Numerical study of natural convection heat transfer in a cylindrical envelope with internal concentric slotted hollow cylinder, Numerical Heat Transfer, Part A 22, no. 3 (1992): 289-305.
10
[13] M. Giangi, F. Stella, T. A. Kowalewski. “Melting Phase change problems with free convection: fixed grid numerical simulation”. Computing and Visualization in Science 2, (1999): 123-130. [14] V. Timchenko, P. Y. P. Chen, E. Leonardi, G. De Vahl Davis, R. Abbaschian. “A computational study of transient plane front solidification of alloys in a Bridgman apparatus under microgravity conditions”. Int. J. Heat Mass Transfer 43, no. 6 (2000): 963-980. [15] M. Mahdaoui, T. Kousksou, S. Blancher, A. Ait Msaad, T. El Rhafiki, M. Mouqallid. “A numerical analysis of solid–liquid phase change heat around a horizontal cylinder”. Applied Mathematical Modelling 38, no. 3 (2014): 1101-1110. [16] M. Abdollahzadeh, M. Esmaeilpour. “Enhancement of phase change material (PCM) based latent heat system with nano fluid and wavy surface”. Int. J. Heat Mass Transfer 80, 2015 ): 376-385. [17] S. Tiari, S. Qiu, M. Mahdavi. “Numerical study of finned heat pipe-assisted thermal energy storage system with high temperature phase change material”. Energy Conversion and Management 89, (2015): 833-842. [18] M. Auriemma, A. Lazzetta. “Numerical analysis of melting paraffin wax with Al2O3 ZnO and CuO nanoparticule in rectangular enclosure”. Indian journal of science and technology 9, no 3, (2016). [19] F. Fornarelli, S. M. Camporeale, B. Fortunato, M. Torresi, T. P. Oresta, L. Magliocchetti, A. Miliozzi, G. Santo. “CFD analysis of melting process in a shell-and-tube latent heat storage for concentrated solar power plants”. Applied Energy 164, (2016): 711-722. [20] M. Parsazadeh, X. Duan. “Numerical and Statistical Study on Melting of Nanoparticle Enhanced Phase Change Material in a Shell-and-Tube Thermal Energy Storage System”. Applied thermal engineering 111, (2017): 950-960. [21] A. D. Brent, V. R. Voller, K. J. Reid. “Enthalpy-porosity technique for modeling convectiondiffusion phase change: Application to the melting of pure Metal”. Numerical Heat transfer 13, no. 3 (1988): 297-318. [22] V. R. Voller, C. Prakash. “On the numerical solution of the mixture model equations describing binary solid-liquid phase change”. Numerical Heat transfer 15, no. 2 (1989): 171-1989. [23] S. V. Patankar. Numerical Heat Transfer and Fluid Flow. McGraw-Hill, 1980. [24] W. Kurz, D. J. Fisher. Fundamentals of solidification. Trans Tech, 1986. [25] O. Bertrand, B. Binet, H. Combeau, S. Couturier, Y. Delannoy, D. Gobin, M. Lacroix, P. Le Quere, M. Medale, J. Mencinger, H. Sadat, G. Vieira. “Melting driven by natural convection. A comparison exercise: First result”. Int. J. of thermal science 38, no. 1 (1999): 5-26. [26] J. A. Dantzig. “Modelling Liquid-Solid Phase Changes with Melt Convection”. Int. J. Numer. Meth. Eng 28, no. 8 (1989): 1769–1785. [27] D. Gobin, P. Lequere. “Melting from an isothermal vertical wall. Synthesis of a numerical comparison exercise”. Computer Assisted Mechanics and Engineering Sciences 7, no. 3 (2000): 289-306. [28] N. Hannoun, V. Alexiades, T.Z.Mai. “Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side”. Numerical heat transfer part B: Applications 44, no. 3 (2003): 253–276. [29] N. Hannoun, V. Alexiades, T. Z. Mai. “A reference solution for phase change with convection”. Int. J. Numer. Meth. in fluid 48, no. 11 (2005): 1283–1308.
Paper ID: CSITE_2019_424 The three authors have all participated in all stages of the preparation of this article.
Prof A. Cheddadi Corresponding author
Declaration of interests X The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
On behalf of the authors Prof Abdelkhalek CHEDDADI
S2214-157X(19)30429-0
DOI:
https://doi.org/10.1016/j.csite.2020.100585
Reference:
CSITE 100585
To appear in:
Case Studies in Thermal Engineering
Received Date: 26 October 2019 Revised Date:
15 December 2019
Accepted Date: 4 January 2020
Please cite this article as: S. Kaba, K. Achoubir, A. Cheddadi, An enhanced approach for convection phase change problems: A solution of tin melting problem, Case Studies in Thermal Engineering (2020), doi: https://doi.org/10.1016/j.csite.2020.100585. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
AN ENHANCED APPROACH FOR CONVECTION PHASE CHANGE PROBLEMS: A SOLUTION OF TIN MELTING PROBLEM Salaheddine Kaba1, Khalid Achoubir1 and Abdelkhalek Cheddadi1,* 1
Mohammed V University in Rabat, Mohammadia School of Engineers, Ibn Sina Str., P.O. Box 765 Agdal, Rabat, Morocco * Corresponding author, Email address : [email protected] Abstract
A comprehensive and efficient numerical phase change model compatible with the stream function-vorticity formulation is developed. The model is based on the source term of the enthalpy-porosity method derived from Darcy's law. The results of the model are compared with those of tin melting at = 0.02, = 2.5 × 10 , = 0.01, in this case the method demonstrates its capability to deal with this controversial case by using an optimal grid. The comparison with results obtained by using this approach was satisfactory. Keywords: melting, phase change materials (PCM), stream function-vorticity, enthalpy-porosity.
Nomenclature Alphabetically listed terms Darcy source term Constant of Darcy source term Thermal capacity (J/kg K) Liquid fraction g Gravity acceleration (m/s²) Thermal conductivity (W/m K) Latent heat (J/kg) Reference length (m) Pr Prandtl number Rayleigh number Ste Stefan number S Source term t Dimensionless time T Dimensionless temperature u Horizontal velocity v Vertical velocity x,y Cartesian coordinates
Greek symbols Thermal diffusivity (m2/s) Thermal expansion coefficient (K-1) υ Kinematic viscosity (m2/s) Vorticity Density (kg/m3) Stream function ∆" Temperature difference between walls δ" Dimensionless half width of mushy zone temperature range Subscripts and superscripts 0 Reference value # Liquidus $ Solidus % Hot Cold ′ Dimensional value
1. Introduction The numerical modeling of phase change remains since several decades an open field of investigation, due to the presence of melting and solidification techniques in diverse technologic domains such as: the fabrication of metallic alloys of high purity, thermal energy storage in term of latent heat, etc.
2
The major challenge of the mathematical modeling of phase change is the moving of the solid-liquid interface. Indeed, the phase change is characterized by a mobile front accompanied by the absorption or desorption of the latent heat with thermophysical properties variation, including the thermal conductivity, the specific heat, the viscosity and the density. Furthermore, most of phase change applications take place in a multidimensional configuration. In addition, the natural convection in liquid phase plays an important role in the heat transfer control which further complicates the problem. Usually, there are two groups of numerical model for studying this kind of problem: the fixed grid schemes and moving grid schemes [1]. The last group uses a transformed system of coordinates and sophisticated mathematics for the interface tracking, which may cause a lot of problems in programming and ineluctably results in a higher computational cost. On the other hand, in the fixed grid schemes, the position of the interface is deduced implicitly from the enthalpy/temperature distribution, leading to relatively simple mathematical models. Three principal methods were proposed for the formulation of the conservation equations of energy: the enthalpy method [2], the equivalent thermal capacity method [3], and the temperature transforming model [4]. Before the application of each of these models, one must express the velocity in the two phases, solid and liquid. The velocity must be of course zero in the solid phase and calculated from the equations of motion in the liquid phase and in the mushy zone in the case of non isothermal phase change. Several methods were proposed for the treatment of this problem. We cite the variable viscosity method [5], the source term method [6-11] and the switch-off method [12]. All of these methods were developed for the equations of motion in terms of velocity-pressure formulation. The using of stream function-vorticity formulation in the phase change problems (melting or solidification) is limited. Giangi et al. [13] used this formulation to study the phase change of a pure element. For stopping the velocity they used a source term without presenting explicitly the method used to derive it. Timchenko et al. [14] used this formulation to simulate the solidification of an alloy at eutectic composition under microgravity conditions. Those authors [13,14] have used an adequate approach for the cases where the interface is very thin and boxed in one calculation cell; the stream function , the vorticity and the velocities were considered nil in the solid phase and calculated from the definition of the stream function in the liquid phase. It should be noted that these works [13,14], didn’t present a method for the treatment of non isothermal phase change, where the existence of mushy zone is indispensable and request an additional relation to identify it. The objective of the present work is to apply an approach for the correction of velocity, which will be adequate with the stream function-vorticity formulation. This method is based on the frequently used Darcy source term approach, especially in the numerical studies of the phase change materials (PCM) dedicated to the thermal energy storage [15,20]. The capability of this approach is demonstrated in the comparison with the controversial case of tin melting. 2. Mathematical model and approximation 2.1 Numerical approximation The fluid considered is supposed Newtonian and incompressible. The flow is laminar and two dimensional with gravity parallel to the ascendant y axis. The thermophysical properties are considered constant except the density which is approached in the gravity term with the Boussinesq hypothesis: = ' (1 − (" ′ − "'′ ), (1) ′ where "' is the reference temperature, ' is the reference density, is the thermal expansion coefficient.
3
The dimensionless variables are defined by: ′ -′ .′ (-, .) = / , 0, = 1 5′
(2, 3) = 4
6
,
7′ 6
8, " =
9 ′ :9;′ ∆9 ′
is the reference length, is taken as
6
>?@A
.
(2)
is the thermal diffusivity, ∆" ′ = "<′ − "=′ . The reference velocity
2.2. Treatment of the equation of motion The stream function-vorticity formulation ( , ) proposed in this paper is derived from the equations of motion (in dimensionless form) written in primitive variables that contain the source term as proposed in [2]: B5 BC B7 BC
B5
B5
B7
B7
+ 2 BE + 3 BF =
+2
BE
+ 3
Where L
=6
BF
=
BH
G 1 2 − BE + 2
G13 −
BH BF
+
(3) (" − "' ) + 3
is a parameter used to correct the velocity,
is the Prandtl number.
(4) =
IJ∆9 ′ >K?@A L6
is the Rayleigh number and
In the context of this work, the final shape of the parameter depends on the approximation used for the flow in the mushy zone (where liquid and solid phases are mixed). In [2, 21], the expression of is obtained by considering the cavity as a porous medium. However in the case of isothermal phase changes, the same model can be used with a dummy mushy zone with an infinitesimal thickness and since the depends on the liquid fraction , a suitable relationship between expression of the source term temperature and liquid fraction must be established as will be seen later to treat this type of transition. Voller and Prakash [2] proposed the following expression which claims to prevent fluid flow in porous media: =−
=(M: N )O ( NK PQ)
(5)
Where is an arbitrary constant (large), and R is also a constant (very small), introduced to avoid the possible division by zero. When the liquid fraction is zero = 0, the control volume is totally solid, therefore the source term tends to the infinity forcing the velocity to switch-off. On the other side when = 1, the source term vanishes and the velocity is calculated through the equations of motion in fluid medium. The expression (5) is used primarily to treat non-isothermal phase changes (wide mushy zone). According to [2], this expression forces the fluid to follow the Carman-Kozny law, which is experimentally validated for flow in dendritic regions. One of the numerical disadvantages of this formulation is the high computational cost (presence of two powers with a division), hence the use of a less complex expression, Eq. (6), is required as was proposed by [11]. = (1 − ) (6) Note that when = 0, one obtains = . If the constant is not well chosen, it is possible that the flow will not be sufficiently turned off. In the present work, we first inject Eq. (6) in Eqs. (3) and (4). By differentiating Eq. (3) with respect to y, Eq. (4) with respect to x and subtracting one from the other (3) - (4). we obtain the stream functionvorticity formulation( , ) with an adequate source term T as follows: BT BC
BT
BT
+ 2 BE + 3 BF =
G1 −
B9 BE
+
T
(7)
4
= − G 1
2=−
BU
With:
T
BF
; 3 =
BU BE
(8)
(9)
= [(1 − ) + 3
B N BE
−2
B N ] BF
(10)
2.3. Treatment of energy equation: According to [22], the dimensional energy conservation equation can be expressed directly in terms of temperature. B9 ′ BC ′
+
(2′
B9 ′ BE ′
+ 3′
B9 ′
BF ′
) = G 1 " ′ + ′ 9
(11)
Where ′ 9 is the source term of energy equation: ′9 = −
> B N
=YN BC ′
−
>
=YN
B
B
(2′ BEN′ + 3 ′ BFN′ )
(12)
Introducing the dimensionless temperature defined in Eq.(2) leads to the following conservation equation in dimensionless form: B9 BC
B9
B9
M B N BC
+ 2 BE + 3 BF = G 1 " − ZC
Where
=
=Y ∆9 ′ >
M
− ZC (2
B N BE
+ 3
B N ) BF
(13)
is the Stefan number.
The convective part of the source term 9 (second term in the right member of Eq. 12) will be neglected in the case of isothermal phase change (e.g melting of pure elements), and will be introduced in the case of non-isothermal phase changes (e.g alloys solidification) as operated in [21]. 2.4. Method of numerical resolution In the present work the equations of motion are discretized with a centered finite difference method, with an adequate implementation of source term given by Eq. (10) in each compute node. On the other side the energy conservation Eq. (13) is discretized with a finite volume method by adopting a staggered grid for " compared to that used for and (Figure 1). The convective terms of Eq. (7,8) are discretized according to a centered scheme while those of Eq. (13) are discretized according to an Upwind scheme. The temporal part of source term of Eq. (13) is treated according to an Euler explicit scheme while the convective part is dicretized based on the method proposed by Patankar [23]. For the resolution we adopt the ADI method (Alterning Direction Implicit) with internal iteration, which generates tridiagonal systems of linear equations. These systems of equations are solved with Thomas algorithm (TDMA). Thus the resolution is performed with a considerable reduction in calculation time; Figure 2 shows a general algorithm of the major steps of the used code. The convergence is established at each time step when the maximum relative deviations are less than -5 10 for the dynamic and thermal fields. The dimensionless time step was between (10-5 - 10-6), the mesh was between (61×61 and 101×101). The refinement of the time and space step is carried out at the beginning of the convective melting regime and at the zone of transition (cell merging). The grid used is given for each simulation (See figure captions)
5
Figure 1. Presentation of mesh with staggered grid. Iterative updating of liquid fraction: To ensure the convergence of Eq. (13), special treatment of the source term is performed at each time step to adapt the liquid fraction in each cell with temperature calculated through the equation of energy. An important property which strongly affects the convergence of the internal iteration process is the procedure used for updating the liquid fraction from one iteration to another. Based on the method presented in [22], using an algorithm considered adequate for the update of the liquid fraction in the case of binary alloy solidification, we propose the following method. The temperature-liquid fraction relationship may be found by the experiment, or drawn from one of the metallurgical relationships available in literature, e.g. Scheil equation [24]. In the present work, the liquid fraction is taken under a linear form as a function of temperature: =
9 [ : 9\[ 9][ :9\[
(14)
In the system which is undergoing a phase change, the total enthalpy may be written as follows: % = "′ + (15) Where at the nth iteration we get: ^
%^ = " ′ + ^ (16) At each time step, and during the updating process by internal iterations, the total enthalpy of each node Eq. (15), is conserved between iterations (n) and (n+1). Only the distribution of the latent and the sensible enthalpy (respectively and " ′ ) will be different. Taking guidance from Eq. (14) and (16), we obtain: ^
"′ + With:
^
=
"′
^PM
^PM
+
^PM
"′ = ^PM (" ′ − " ′ _ ) + " ′ _ From which one gets: ^PM
=
` =Y 49 ′ :9 ′ a 8P N` >
>P1b9 ′ =Y
(17) (18)
(19)
Where c" ′ is a numerical value representing the half width of mushy zone temperature range. In the case of melting/solidification of a mixture of multiconstituents, this parameter is defined from the composition of each constituent. In the case of pure substances, a fictive value of c" ′ is introduced, which does not represent the physical reality, although it is used as a parameter for bypassing the discontinuity associated with the transition from one phase to another and avoid the numerical complexity associated with this transition. Practically c" ′ is selected so as to have " − "_ ∿ 0.
6
Finally, Eq. (19) can be written in the dimensionless form as follows: ^PM
=
ZC (9 ` :b9)P N` MP1 b9 ZC
Where c" is the dimensionless half width of mushy zone temperature range.
(20)
3. Results and discussion The simulations are undertaken in a square cavity (Figure 2) where the reference length is taken equal to the side of cavity. The Tin melting in a square cavity differentially heated, Bertrand et al. [25], is characterized by the presence of instabilities at the start of the convective regime. These instabilities have been identified in the work of Dantzig on the Gallium melting [26], who links this behavior to the instabilities encountered at the end of the conductive regime.
Figure 2. Geometric configuration of melting problem. In the case of melting, the experiment starts when one of the walls is heated to the melting point of Tin, then the heat is conveyed through a liquid layer [21]. Therefore the physical properties of the liquid will govern the melting process. The Tin physical properties are shown in Table 1. Table 1. Tin thermal parameters for melting problem Parameter Cp K ρ υ L TH TC Lref g Ra Pr Ste
Value 200.0 60.0 7.5 × 103 4.0 × 10-5 8.0 × 10-7 2.67 × 10-4 6 × 104 508.0 505.0 0.1 10 2.5 × 105 0.02 0.01
Unit J/kg K W/m K Kg/m3 m2/s m2/s K-1 J/kg K K m m/s2
The thermal parameters used as cited in the benchmark of Bertrand et al. [25] are: = 0.02, = 2.5 × 10 , = 0.01. The grid used in our study varies between 61×61 for the majority of cases and 101×101 for the detection of instabilities at the start of the convective regime.
7
At t = 1, Figure 3-a shows a good concordance regarding the shape and the location of the interface of fusion, between the participants in the benchmark on the one hand and the present work on the other. Meanwhile, a clear discrepancy between the participants is noticed at evolved time, especially at t = 4 (Figure 3-b). Some of them found an interface with two bumps which correspond to two co-rotative convective cells, while the others obtained an interface with one bump resulting from one single convective cell. The obtained results within the framework of the present paper are in agreement with the first group of participants (Figure 3-b). The interface shape and the convective cells number (Figure 4) confirms this observation. This contradiction was discussed by Gobin and Lequere [27] in a synthesis that followed this comparison exercise. The authors concluded that the simulation at low number tends to confirm the result obtained by Dantzig [26]. A more recent numerical investigation [28] was dedicated to resolve the controversy. It showed that by refining the grid (from 100×100 up to 400×400) and by using various numerical schemes (Figure 4), it becomes clear that the solution giving two cells is numerically more stable and that the one cell regime appears only when a coarse grid is used. These comparisons show the capacity of the approach used in this work to detect the instabilities that are detectable only at very fine grids in previous works [28, 29] (Figure 4). These instabilities, manifested by the shape of the interface and the number of convective cells which appear at successive times, are solved in the present work with a relatively coarse grid compared with the one used in [28]. We note as well from Figure 4 that the number of co-rotating cells is high at the beginning of the convective regime, and starts to reduce (recombining of convective vortices) to give rise to one cell dominated structure toward t = 8 (Figures 4 and 3-c). In Figure 5, we represent the evolution of the flow intensity efE as a function of time, showing the flow structure at each transition. According to the shape of the curve representing efE evolution, we note that the transitions, which correspond to a merging of convective cells, are associated with more or less important fluctuation in the value of efE . Table 2 presents the transition times and the resulting regime. Table 2. Regimes and instants of cells merging. Grid 100×100 Time
Flow regime
0.3
Beginning of 4 cells regime
0.6-0.8
Merging : 4 cells into 3 cells
1.8-2
Merging : 3 cells into 2 cells
4.5-5.5
Merging : 2 cells into 1 cell
8
(a) t = 1
(b) t = 4
(c) t = 8
Figure 3. Evolution of melting interface: (a,b) Comparison with Bertrand et al. [25] c) Comparison with Hannoun et al. [28]. = 0.02, = 2.5 × 10 , = 0.01. Grid: 61×61.
(a) t = 0.8
(b) t = 4
(c) t = 8
Figure 4. Stream lines and melting interface (Part of cavity) at different times. Pr = 0.02 , Ra = 2.5 × 10 , Ste = 0.01. Hannoun et al [28] (left) Grid: (100×100 up to 400×400). Present work (right) Grid: 101×101 (t = 0.8). Grid: 61×61 (t = 4 and t = 8).
Figure 5. Evolution of maximum value of stream function efE as function of time. = 2.5 × 10 , = 0.01. Grid: 101×101
= 0.02 ,
9
4. Conclusion In this work, a numerical method for the correction of the velocity adequate with the stream functionvorticity formulation was adopted, based on the source term approach. This method is valid for modeling both the isothermal and non-isothermal phase changes. The comparison with the case of the melting (low Prandlt number) is investigated, = 0.01. The results were in agreement with the reference solutions available in literature [28,29]. According to the results of this work, we can conclude that the coupling of the enthalpy porosity techniques with the stream function-vorticity formulation gives an acceptable solution for the controversial case of Tin melting problem with a relatively simple implementation of the code and a very optimal computational time.
References [1] A. A. Samarskii, P. N. Vabischchevich, O. P. Iliev, A. G. Churbanov.”Numerical Simulation of Convection/Diffusion Problems: A Review”. Int. J. Heat Mass Transfer 36 , no. 17 (1993): 40954106, [2] V. R. Voller, C. Prakash.”A fixed grid numerical modeling methodology for convection– diffusion mushy region phase-change problems”. Int. J. Heat Mass Transfer 30, no. 12 (1987): 1079–1718. [3] K. Morgan.” A numerical analysis of freezing and melting from a vertical wall”. Comput. Methods Appl. Mech. Eng 28, no. 3 (1981): 275–284,. [4] Y. Cao, A. Faghri.”A numerical analysis of phase change problem including natural convection”. ASME J. Heat Transfer 112, no. 3 (1990): 812–815. [5] D. K. Gartling. Finite element analysis of convective heat transfer problems with change of phase. Computer Methods in Fluids. Pentech, 1980. [6] B. Binet, L. Lacroix. ”Melting from heat sources flush mounted on a conducting vertical wall”. Int. J. Num. Meth. Heat & Fluid Flow 10, no. 3 (2000): 286-306. [7] V. R. Voller, M. Cross, N. C. Markatos. ”An enthalpy method for convection/diffusion phase change”. Int. J. Num. Meth. Eng 24, no. 1 (1987): 271-84. [8] P. A. Galione, O. Lehmkuhl, J. Rigola, A. Oliva. “Fixed grid modeling of solid-liquid phase change in unstructured meshes using explicit times schemes”. Numerical heat transfer part B: Fundamentals 65, no. 1 (2014): 27-52. [9] M. Ghalambaz, A. Doostanidezfuli, H. Zargartalebi, A. J. Chamkha. “MHD phase change heat transfer in an inclined enclosure: effect of a magnetic field and a cavity inclination”. Numerical heat transfer part A: Applications 71, no. 1 (2017): 91-109. [10] Y. Belhamadia, A. S. Kane, A. Fortin. “An enhanced mathematical model for phase change problems with natural convection”. Int. J. of numerical analysis and modeling, series B 3, no. 3 (2012): 192-206. [11] A. A. Al-abidi, S. B. Mat, K. Sopian, M. Y. Sulaiman, A. Th. Mohammed. “CFD applications for latent heat thermal energy storage: a review”. Renewable and Sustainable Energy Reviews 20, (2013): 353-363. [12] M. Yang, W. Q. Tao. ”Numerical study of natural convection heat transfer in a cylindrical envelope with internal concentric slotted hollow cylinder, Numerical Heat Transfer, Part A 22, no. 3 (1992): 289-305.
10
[13] M. Giangi, F. Stella, T. A. Kowalewski. “Melting Phase change problems with free convection: fixed grid numerical simulation”. Computing and Visualization in Science 2, (1999): 123-130. [14] V. Timchenko, P. Y. P. Chen, E. Leonardi, G. De Vahl Davis, R. Abbaschian. “A computational study of transient plane front solidification of alloys in a Bridgman apparatus under microgravity conditions”. Int. J. Heat Mass Transfer 43, no. 6 (2000): 963-980. [15] M. Mahdaoui, T. Kousksou, S. Blancher, A. Ait Msaad, T. El Rhafiki, M. Mouqallid. “A numerical analysis of solid–liquid phase change heat around a horizontal cylinder”. Applied Mathematical Modelling 38, no. 3 (2014): 1101-1110. [16] M. Abdollahzadeh, M. Esmaeilpour. “Enhancement of phase change material (PCM) based latent heat system with nano fluid and wavy surface”. Int. J. Heat Mass Transfer 80, 2015 ): 376-385. [17] S. Tiari, S. Qiu, M. Mahdavi. “Numerical study of finned heat pipe-assisted thermal energy storage system with high temperature phase change material”. Energy Conversion and Management 89, (2015): 833-842. [18] M. Auriemma, A. Lazzetta. “Numerical analysis of melting paraffin wax with Al2O3 ZnO and CuO nanoparticule in rectangular enclosure”. Indian journal of science and technology 9, no 3, (2016). [19] F. Fornarelli, S. M. Camporeale, B. Fortunato, M. Torresi, T. P. Oresta, L. Magliocchetti, A. Miliozzi, G. Santo. “CFD analysis of melting process in a shell-and-tube latent heat storage for concentrated solar power plants”. Applied Energy 164, (2016): 711-722. [20] M. Parsazadeh, X. Duan. “Numerical and Statistical Study on Melting of Nanoparticle Enhanced Phase Change Material in a Shell-and-Tube Thermal Energy Storage System”. Applied thermal engineering 111, (2017): 950-960. [21] A. D. Brent, V. R. Voller, K. J. Reid. “Enthalpy-porosity technique for modeling convectiondiffusion phase change: Application to the melting of pure Metal”. Numerical Heat transfer 13, no. 3 (1988): 297-318. [22] V. R. Voller, C. Prakash. “On the numerical solution of the mixture model equations describing binary solid-liquid phase change”. Numerical Heat transfer 15, no. 2 (1989): 171-1989. [23] S. V. Patankar. Numerical Heat Transfer and Fluid Flow. McGraw-Hill, 1980. [24] W. Kurz, D. J. Fisher. Fundamentals of solidification. Trans Tech, 1986. [25] O. Bertrand, B. Binet, H. Combeau, S. Couturier, Y. Delannoy, D. Gobin, M. Lacroix, P. Le Quere, M. Medale, J. Mencinger, H. Sadat, G. Vieira. “Melting driven by natural convection. A comparison exercise: First result”. Int. J. of thermal science 38, no. 1 (1999): 5-26. [26] J. A. Dantzig. “Modelling Liquid-Solid Phase Changes with Melt Convection”. Int. J. Numer. Meth. Eng 28, no. 8 (1989): 1769–1785. [27] D. Gobin, P. Lequere. “Melting from an isothermal vertical wall. Synthesis of a numerical comparison exercise”. Computer Assisted Mechanics and Engineering Sciences 7, no. 3 (2000): 289-306. [28] N. Hannoun, V. Alexiades, T.Z.Mai. “Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side”. Numerical heat transfer part B: Applications 44, no. 3 (2003): 253–276. [29] N. Hannoun, V. Alexiades, T. Z. Mai. “A reference solution for phase change with convection”. Int. J. Numer. Meth. in fluid 48, no. 11 (2005): 1283–1308.
Paper ID: CSITE_2019_424 The three authors have all participated in all stages of the preparation of this article.
Prof A. Cheddadi Corresponding author
Declaration of interests X The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
On behalf of the authors Prof Abdelkhalek CHEDDADI