Numerical simulation of solidification around a circular cylinder with natural convection

Vietnam Journal of Mechanics, VAST, Vol. 38, No. 4 (2016), pp. 295 – 306 DOI:10.15625/0866-7136/7442 NUMERICAL SIMULATION OF SOLIDIFICATION AROUND A CIRCULAR CYLINDERWITH NATURAL CONVECTION Vu Van Truong∗, Truong Viet Anh, Hoang Thi Bich Ngoc Hanoi University of Science and Technology, Vietnam ∗E-mail: vuvantruong.pfae@gmail.com Received November 21, 2015 Abstract. An analysis is carried out for solidification around a cold cylinder in a rectangu- lar cavity using numerical simulations. The transient influences of solidification accompa- nied by natural convection are investigated in detail. The governing equations, in terms of one-fluid formulation, including the Navier-Stokes and energy equations for incompress- ible Newtonian fluids are written for the whole domain. The no-slip and Dirichlet-type isothermal temperature boundary conditions are both implemented using a linear interpo- lation technique (immersed boundary method). The solidification interface is represented by connected elements that move on the fixed background grid. Code validations are car- ried out through various problems. Finally, the temporal dependence of the solid area ratio of the solid phase to the cylinder upon various dimensionless parameters, such as Rayleigh number, Prandtl number, Stefan number, thermal property ratios as well as a parameter indicating the effect of superheat is studied. Keywords: Solidification, cylinder, natural convection, front tracking, immersed boundary method. 1. INTRODUCTION Solid-liquid phase change finds application in many systems and in nature in- cluding latent heat energy storage, metallurgy, and food and pharmaceutical process- ing. Solidification around cooled cylinders appears in thermal energy storage systems, freezing of soil and others. Therefore, understanding solidification heat transfer around cooled cylinders plays an important role in designing and operating such systems. Ac- cordingly, there have been many works concerned with solidification around cylinders. Sasaguchi and co-workers [1] presented a numerical method for the solid-liquid phase change around a single cylinder and two horizontal cylinders in a rectangular cavity. They used a generalized coordinate system and their results were restricted to the freez- ing of water. A more complicated system in which there exist three or more cylinders has been investigated by Sugawara and Beer [2]. The authors used the PHOENICS code of a commercial software package. In another work [3], an enthalpy formulation-based c© 2016 Vietnam Academy of Science and Technology 296 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc fixed-grid approach was employed to investigate the solidification around a number of staggered cylinders. However, in the above-mentioned papers, the results were restricted to water, and effects of some dimensionless parameters have not been considered. In the present study, we present a numerical method for simulating the solidifi- cation around a cooled cylinder with natural convection. The method utilizes the front tracking technique [4] to represent the solid-liquid interface and an interpolation tech- nique, i.e. an immersed boundary method, to deal with the no-slip and constant isother- mal temperature boundary conditions. Effects of some dimensionless parameters are also investigated. 2. MATHEMATICAL FORMULATION Fig. 1 shows the investigated problem, a solidification layer formed around a cooled circular cylinder in a rectangular cavity. The fluid and thermal properties of each phase are assumed constant, and the densities of the solid and liquid phases are identical (i.e. no volume change). In the liquid region, the flow assumed incompressible is driven by buoyancy-induced natural convection. We treat all phases as one fluid with variable properties such as density ρ, viscosity µ, thermal conductivity k and heat capacity Cp. In terms of the single-field representation, the momentum and thermal energy equations are ∂ ∂t (ρ0u) +∇ · (ρ0uu) = −∇p +∇ · µ(∇u+∇uT) + ρ0f+ ρ0g [1− β (T − T0)] , (1) ρ0Cp ∂T ∂t + ρ0Cp∇ · (Tu) = ∇ · (k∇T) + ∫ f q˙δ ( x− x f ) dS + ρ0Cph, (2) ∇ · u = 0. (3) (a) (b) Fig. 1. Solidification around a cold cylinder in a rectangular cavity: (a) computational domain and (b) front-tracking representation for the solidification interface Numerical simulation of solidification around a circular cylinder with natural convection 297 Here, u is the velocity vector, p is the pressure, g is the gravitational acceleration. T and the superscriptTdenote the temperature and the transpose. D/Dt is the material de- rivative. f is the momentum forcing used to impose no-slip condition on the solid-liquid interface [5]. The last term in Eq. (1) is the Boussinesq approximation for density changes due to thermal gradients [6]. h is the energy forcing term used to impose a constant tem- perature on the cylinder boundary. ρ0 is the density at the reference temperature T0, e.g. fusion temperature. The Dirac delta function δ(x− x f ) is zero everywhere except for a unit impulse at the solidification interface x f . β is the thermal expansion coefficient of the melt. q˙ is the heat source at the solidification interface, given as q˙ = ks ∂T ∂n ∣∣∣∣ s − kl ∂T∂n ∣∣∣∣ l = −ρsVnLh , (4) where the subscripts s and l represent solid and liquid, respectively. Vn is the velocity normal to the solidification front (Fig. 1) and Lh is the latent heat. These above-mentioned equations are solved by the front-tracking method com- bined with interpolation techniques on a staggered grid with second order accuracy in time and space. 3. NUMERICAL METHOD AND VALIDATION The numerical technique used in the present study is based on the front-tracking method [4] and an interpolation technique similar to the immersed boundary proposed in [5]. The momentum and mass conservation equations are discretized using an explicit predictor-corrector time-integration method and a second-order centered difference ap- proximation for the spatial derivatives. The discretized equations are solved on a fixed, staggered grid using the MAC method [7]. The interface between the solid and liquid phases are represented by finite discrete points on a stationary grid, as shown in Fig. 1(b). These interface (or front) points propagate by xn+1f = x n f + n f Vn∆t, (5) where Vn is simply given by Vn = −q˙ f / (ρsLh), (6) and the heat source q˙ is calculated using a normal probe technique [8, 9]. Positions of the interface points are used to determine an indicator function I which is zero in the solid phase and one the liquid phase ∇I = ∫ f δ ( x− x f ) n f dS, (7) where n f is the unit normal to the interface. Accordingly, I = 0.5 corresponds to the interface. The indicator function I is reconstructed from the new position of the interface points. The values of the material properties are then found: φ = φs (1− I) + Iφl, where φ stands for ρ, µ, Cp, or k. To impose no-slip condition on the solid-liquid interface as well as zero velocities within the solid phase, we use a method similar to the direct volume 298 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc forcing method of Liao et al. [5]. We use a similar technique to determine the energy forcing inside, on and near the boundary to satisfy a constant temperature Tc on the cylinder surface. Detailed solution procedures, which based on the projection method, are as follows. Suppose n time steps have been completed, to calculate the solution at time level n + 1 carry out the following steps: 1. Calculate the heat source, and update the solidification interface points (Eq. (5)). 2. Reconstruct the indicator function and update the material properties. 3. Calculate a provisional velocity field u∗∗ = un − ∆t { ∇ · (uu)n + µ∇ · [( ∇un +∇ ( uT )n)]/ ρ0 + g [1− β (Tn − T0)] } . (8) 4. Calculate a provisional temperature field T∗= ρ0Cnp Tn + ∆t −ρ0∇ · (Cnp Tnun)+∇ · (kn∇Tn) + ∫ f q˙δ ( x− x f ) dS  /(ρ0Cn+1p ). (9) 5. Establish locations of the momentum forcing points near the solidification in- terface and locations of the energy forcing points near the cylinder interface (see [5]for details). 6. Obtain the forcing term f f = ( u f − u∗∗ )/ ∆t, (10) with u f is the velocity, at the forcing points, interpolated from u∗∗, or zero, at points inside the solid and cylinder [5]. 7. Obtain the energy forcing term h h = ( Tf − T∗ )/ ∆t, (11) with Tf is the temperature, at the energy forcing points near the cylinder interface, inter- polated from T∗, or Tc, at the points inside the cylinder. 8. Update the temperature field at time n + 1 Tn+1 = T∗ + h∆t. (12) 9. Calculate an intermediate velocity field u∗ = u∗∗ + f∆t. (13) 10. Find the pressure field by solving the Poisson equation( ∇ · un+1 −∇2 p )/ ρ0 = −(∇ · u∗) / ∆t, (14) with ∇ · un+1 = 0. 11. Compute divergence-free time level n+1 fluid velocity un+1 = u∗ − ∆t∇h p / ρ0. (15) Numerical simulation of solidification around a circular cylinder with natural convection 299 This solution procedure for time integration is first order, to produce a second- order scheme, we use a technique described in Esmaeeli and Tryggvason [10] and in [11]. To validate the method, we first compare the computational results with exact so- lutions for a two-dimensional Stefan problem [12] in which a line heat source Q causes a circular solid seed at the center to evolve in the direction of increasing the radius of the seed. This comparison is to validate the front-tracking method for solidification. Fig. 2 shows that the computational results agree well with the exact solutions at Stefan number = 0.1, heat source Q = 40, Cpsl = 0.5 and ksl = 2. Detailed description of this comparison and others can be found in [4]. (a) (b) Fig. 2. Comparisons with exact solutions of Carslaw and Jaeger [12] for a 2D Stefan solidification problem: (a) evolution of the solidification front and (b) temperature along y = 0.5. The data is plotted at time t = 0.03, 0.11, 0.19 and 0.27 Next, to validate the technique of the momentum and energy forcing, we per- formed computations of natural convection in an annulus [13, 14]. We use the same con- figuration and parameters as described in Wang and co-workers [14] (see Fig. 3). The inner cylinder has a radius Ri = 0.625 and temperature Th, while the outer cylinder has a radius R0 = 1.625 and temperature Tc. The two cylinders are placed at the cen- ter of the 5 × 5 domain. The Prandtl number is Pr = 0.7 and the Rayleigh number is Ra = gβ(Th − Tc)(R0 − Ri)3 / (να) = 5× 104 where ν and α are the kinematic viscosity and the thermal diffusivity, respectively. The distribution of the temperature field in the annulus at the angle ϕ = 600 is shown in Fig. 3(b). Good agreement between our calcula- tion and the result from Kuehn and Goldstein [13], as well as the result of Wang et al. [14] (see Fig. 3(b)) supports the accuracy of the present method. 300 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc (a) (b) Fig. 3. Natural convection in an annulus: (a) domain for computation and (b) temperature variation in the annulus at ϕ = 600 in comparison with data reported in Kuehn and Goldstein [13] and in Wang et al. [14] 4. RESULTS AND DISCUSSION We turn to our main problem, i.e. solidification around a circular cylinder (Fig. 1). We choose the diameter d of the cylinder as a scaling length, and τc = ρlCpld 2 / kl as the characteristic time scale. The characteristic velocity scale is thus taken to be Uc = d / τc. With these above choices, it is possible to show that the dynamics of the problem is governed by the Prandtl number Pr, the Rayleigh number Ra, the dimensionless initial temperature in the liquid θ0, and the ratios of the thermal properties of the solid to liquid St = Cpl (Tm − Tc) Lh , Pr = Cplµl kl , Ra = gβ(Ti − Tm)d3 νlαl , θ0 = Ti − Tc Tm − Tc , Cpsl = Cps Cpl , ksl = ks kl . (16) The temperature is non-dimensionalized as θ = (T − Tc) / (Tm − Tc). The dimen- sionless time is τ = t/τc. The center of the cylinder is at x/d = (0, 3) in a domain of 2 × 6. For parametric investigations of phase change problems, the values of these parameters are varied in various ranges. For example, for liquid metals or semiconduc- tors Pr is varied in the range of 0.01 to 0.1, Ra is varied in the range of 103 to 105, St is varied in the range of 0.01 to 1.0, and θ0 is varied in the range of 0.1 to 2.0 [15–17]. Accordingly, in the present study, the typical values of the parameters are chosen as St = 0.1, Pr = 0.1, Ra = 1 × 105, Cpsl = 1.0, ksl = 1.0, and θ0 = 2. These values cor- respond to a liquid metal [15, 16] or a semiconductor material [17]. We first perform grid refinement study for our solidification problem (Fig. 1) to choose the best grid. To do so, we perform simulations with the typical values of the parameters, i.e. St = 0.1, Pr = 0.1, Ra = 1 × 105, Cpsl = 1.0, ksl = 1.0, and θ0 = 2, Numerical simulation of solidification around a circular cylinder with natural convection 301 for three grid resolutions 64× 192, 128× 384 and 256× 768 as shown in Fig. 4. Fig. 4(a) shows the variation with respect to time of the ratio As/Ac of the solid phase area to the cylinder area. Comparison of the front locations obtained with two finer grids is shown in Fig. 4(b). As shown in this figure, the grid 64× 192 yields the results different from those obtained by the two finer grids. Contrarily, the results obtained from 128× 384 and from 256× 768 are almost identical. Accordingly, we use 128× 384 for all computations presented below. (a) (b) Fig. 4. Grid refinement study with St = 0.1, Pr = 0.1, Ra = 1× 105, Cpsl = 1.0, ksl = 1.0, and θ0 = 2: (a) temporal variation of the area ratio of the solid phase to the cylinder (As/Ac); and (b) location of the solidification front at time τ = 4.2, 8.4 and 12.6 (a) (b) (c) Fig. 5. Solidification around a cylinder with St = 0.1, Pr = 0.1, Ra = 1× 105, Cpsl = 1.0, ksl = 1.0, and θ0 = 2. In each frame, the left shows the isotherms (plotted every ∆θ = 0.4) and the right shows the velocity field normalized by Uc (the reference vector has a magnitude of 50). The dash line is the solidification front 302 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc Fig. 5 shows the temporal evolution of the solidification front with temperature and velocity fields. At early times of the solidification process, downward flow arises along the solid–liquid interface (the dash line in Fig. 5), at which the temperature is θm = 1.0, because the density of the liquid increases with a decrease in the temperature. As shown in Fig. 5(a) and 5(b), the downward flow is strong with the cooled liquid accumulating at the bottom part of the cavity. As a result, a thermally stratified region is gradually formed. As time progresses, the temperature of the liquid phase decreases to near the fusion temperature θm. Consequently, the downward flow is suppressed. At τ = 6.6 (Fig. 5(c)), the liquid temperature is uniform in the entire cavity, and almost no flows are evident at this time, and the solidification process is merely controlled by conduction. As shown in Fig. 5(a) and 5(b) that the solidification growth rate is slightly faster around the lower part of the cylinder due to the thermal stratification. Next we consider the effects of some parameters on the solidification process. 4.1. Effects of the initial temperature Fig. 6(a) compares the velocity field in the lower part of the cavity for two cases: θ0 = 1.25 (left) and 2.25 (right) at an early time of the solidification process, i.e. τ = 1.68. It is evident that the solidification rate for θ0 = 2.25 is smaller than that for θ0 = 1.25. This is understandable since the warm liquid flows more strongly towards the solidification interface for θ0 = 2.25 than it does for θ0 = 1.25 (Fig. 6(a)). As a result, smaller initial liquid temperatures lead to larger solid area ratios As/Ac as shown in Fig. 6(b) . (a) (b) Fig. 6. Effects of the initial temperature on the solidification: (a) velocity field at τ = 1.68 for θ0 = 1.25 (left) and θ0 = 2.25 (right); (b) evolution of the solid area ratio. St = 0.1, Pr = 0.1, Ra = 1× 105, Cpsl = 1.0, and ksl = 1.0 4.2. Effects of Ra and Pr Fig. 7(a) shows the effect of Ra on the solidification process. At early times, as Ra increases the warm fluid flows more strongly towards the solidification interface, result- ing in a decrease in the solidification rate (see the lower set in Fig. 7(a)). However, strong Numerical simulation of solidification around a circular cylinder with natural convection 303 circulation corresponding to high values of Ra causes liquid temperature to fast reach the fusion temperature of the solidification interface. Accordingly, as time progresses the solidification rate becomes larger for higher Ra (the upper set in Fig. 7(a)). Consequently, the solid area ratio (As/Ac) is larger for higher Ra as shown in Fig. 7(a). (a) (b) Fig. 7. Effects of (a) the Rayleigh number with Pr = 0.1, and (b) the Prandtl number with Ra = 1× 105. Other parameters: St = 0.1, Cpsl = 1.0, ksl = 1.0, and θ0 = 2 The influence of the Prandtl number Pr on the area ratio is shown in Fig. 7(b). As Pr increases the solid area ratio increases. This is understandable since increasing Pr corresponds to increasing the viscous effect. Accordingly, the warm liquid flows more strongly towards the solidification interface at smaller Pr, resulting in a decrease in the solidification rate (Fig. 7(b)). 4.3. Effects of St and ksl The effect of St on the solidification process is shown in Fig. 8(a). According to Eq. (6), increasing St corresponding to decreasing the latent heat released (see Eq. (16)) results in an increase in the growth rate of the solidification front, and thus increases the solid area ratio. We expect that the solidification process is considerably affected by the solid-to- liquid thermal conductivity ratio ksl . Fig. 8(b) confirms our expectation. According to Eq. (4), a higher solid-to-liquid thermal conductivity ratio produces higher heat flux at the solidification interface. Thereby, increasing this ratio results in an increase in the solidification rate. As a result, the solid area ratio is larger at higher ksl , as clearly shown in Fig. 8(b). 4.4. Effects of Cpsl A comparison of velocity vectors in the lower part of the cavity (at τ = 1.8) for Cpsl = 0.5 (left) with those for Cpsl = 1.5 (right) is shown in Fig. 9(a). We see two cir- culations on the left (Cpsl = 0.5) while stronger ones appear on the right (Cpsl = 1.5). In other words, the warm fluid flows more strongly towards the solidification interface for 304 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc (a) (b) Fig. 8. Effects of (a) the Stefan number with ksl = 1.0, and (b) the solid-to-liquid thermal conductivity ratio ksl with St = 0.1. Other parameters: Ra = 105, Pr = 0.1, Cpsl = 1.0 and θ0 = 2 higher Cpsl , resulting in a decrease in the solidification rate. Accordingly, the solid-liquid interface proceeds faster as Cpsl decreases. The solid area ratio is thus smaller at higher Cpsl , as clearly shown in Fig. 9(b). (a) (b) Fig. 9. Effects of the solid-to-liquid heat capacity ratio Cpsl with St = 0.1, Ra = 105, Pr = 0.1, ksl = 1.0, and θ0 = 2: (a) velocity field normalized by Uc at τ = 1.8 presented for Cpsl = 0.5 (left) and for Cpsl = 1.5 (right); and (b) temporal variation of the solid area ratio 5. CONCLUTION An investigation has been carried out for the solidification, under the influence of natural convection, around a cooled cylinder in a rectangular enclosure. The numerical Numerical simulation of solidification around a circular cylinder with natural convection 305 method including the front-tracking technique representing for the solidification inter- face and an immersed boundary technique for the no-slip and Dirichlet-type isothermal temperature boundary conditions has been used. The Navier-Stokes and enery equations have been solved on the whole domain with the fixed and rectangular grid. The method has been validated by comparisons with the exact solutions for a 2D Stefan solidification problem, and with the natural convection solutions in an annulus reported in [13, 14]. For the problem of solidification around a cylinder, the influences of various di- mensionless parameters on the solidification process have been investigated. Numerical results show that increasing Ra, Pr, St or ksl increases the solidification rate. Contrarily, the solidification interface is found to proceed faster as the initial liquid temperature θ0 or Cpsl decreases. ACKNOWLEDGMENT This research is funded by Vietnam National Foundation for Science and Technol- ogy Development (NAFOSTED) under grant number 107.03-2014.21. We are grateful to Prof. John C. Wells at Ritsumeikan University (Japan) for facilitating our computing re- sources. REFERENCES [1] K. Sasaguchi, K. Kusano, and R. Viskanta. A numerical analysis of solid-liquid phase change heat transfer around a single and two horizontal, vertically spaced cylinders in a rectan- gular cavity. International Journal of Heat and Mass Transfer, 40, (6), (1997), pp. 1343–1354. doi:10.1016/s0017-9310(96)00210-4. [2] M. Sugawara and H. Beer. Numerical analysis for freezing/melting around verti- cally arranged four cylinders. Heat and Mass Transfer, 45, (9), (2009), pp. 1223–1231. doi:10.1007/s00231-009-0492-y. [3] Y.-C. Shih and H. Chou. Numerical study of solidification around staggered cylinders in a fixed space. Numerical Heat Transfer, Part A: Applications, 48, (3), (2005), pp. 239–260. doi:10.1080/10407780590945579. [4] T. V. Vu, G. Tryggvason, S. Homma, J. C. Wells, and H. Takakura. A front-tracking method for three-phase computations of solidification with volume change. Journal of Chemical Engi- neering of Japan, 46, (11), (2013), pp. 726–731. doi:10.1252/jcej.13we169. [5] C.-C. Liao, Y.-W. Chang, C.-A. Lin, and J. McDonough. Simulating flows with moving rigid boundary using immersed-boundary method. Computers & Fluids, 39, (1), (2010), pp. 152– 167. doi:10.1016/j.compfluid.2009.07.011. [6] H. Gan, J. Chang, J. J. Feng, and H. H. Hu. Direct numerical simulation of the sedimentation of solid particles with thermal convection. Journal of Fluid Mechanics, 481, (2003), pp. 385–411. doi:10.1017/s0022112003003938. [7] F. H. Harlow, J. E. Welch, et al. Numerical calculation of time-dependent viscous incom- pressible flow of fluid with free surface. Physics of fluids, 8, (12), (1965), pp. 2182–2189. doi:10.1063/1.1761178. [8] N. Al-Rawahi and G. Tryggvason. Numerical simulation of dendritic solidification with con- vection: two-dimensional geometry. Journal of Computational Physics, 180, (2), (2002), pp. 471– 496. doi:10.1006/jcph.2002.7092. 306 Vu Van Truong, Truong Viet Anh, Hoang Thi Bich Ngoc [9] A. Esmaeeli and G. Tryggvason. Direct numerical simulations of bubbly flows. part 1. low reynolds number arrays. Journal of Fluid Mechanics, 377, (1998), pp. 313–345. doi:10.1017/s0022112098003176. [10] A. Esmaeeli and G. Tryggvason. Computations of film boiling. part i: numerical method. International Journal of Heat and Mass Transfer, 47, (25), (2004), pp. 5451–5461. doi:10.1016/j.ijheatmasstransfer.2004.07.027. [11] G. Tryggvason, R. Scardovelli, and S. Zaleski. Direct numerical simulations of gas-liquid multi- phase flows. Cambridge University Press, (2011). doi:10.1017/cbo9780511975264. [12] H. S. Carslaw and J. C. Jaeger. Conduction of heat in solids. Oxford University Press, USA, (1986). [13] T. H. Kuehn and R. J. Goldstein. An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders. Journal of Fluid Mechanics, 74, (04), (1976), pp. 695–719. doi:10.1017/s0022112076002012. [14] Z. Wang, J. Fan, K. Luo, and K. Cen. Immersed boundary method for the simulation of flows with heat transfer. International Journal of Heat and Mass Transfer, 52, (19), (2009), pp. 4510– 4518. doi:10.1016/j.ijheatmasstransfer.2009.03.048. [15] N. Ramachandran, J. Gupta, and Y. Jaluria. Thermal and fluid flow effects during solidifica- tion in a rectangular enclosure. International Journal of Heat and Mass Transfer, 25, (2), (1982), pp. 187–194. doi:10.1016/0017-9310(82)90004-7. [16] D. Johnson and R. Narayanan. A note on the stability of a solidifying material in an under- cooled melt in the presence of convection. In Interactive Dynamics of Convection and Solidifica- tion, pp. 25–32. Springer, (2001). doi:10.1007/978-94-015-9807-1 4. [17] B. T. Murray, S. R. Coriell, G. B. McFadden, A. A. Wheeler, and B. V. Saunders. Gravitational modulation of thermosolutal convection during directional solidification. Journal of Crystal Growth, 129, (1), (1993), pp. 70–80. doi:10.1016/0022-0248(93)90435-y.