International Journal of Heat and Mass Transfer 127 (2018) 302–312 Contents lists available at ScienceDirect International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt A numerical study of a liquid drop solidifying on a vertical cold wall Vinh Nguyen Duy a,b, Truong V. Vu c,⇑ a Department for Management of Science and Technology Development, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Faculty of Electrical & Electronics Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam c School of Transportation Engineering, Hanoi University of Science and Technology, 1 Dai Co Viet, Hai Ba Trung, Hanoi, Viet Nam b a r t i c l e i n f o Article history: Received 21 May 2018 Received in revised form 29 July 2018 Accepted 9 August 2018 Keywords: Front-tracking Drop Solidification Numerical simulation Vertical wall a b s t r a c t Solidification of a liquid drop on a vertical wall is a typical phase change heat transfer problem that exists widely in natural and engineering situations. In this study, we present the fully resolved two-dimensional simulations of such a problem by a front-tracking/finite difference method. Because of gravity, the liquid drop assumed stick to the cold wall shifts to the bottom during solidification. Numerical results show that the conical tip at the solidified drop top is still available with the presence of volume expansion, but the tip location is shifted downward, resulting in an asymmetric drop after complete solidification. The tip shift, height and shape of the solidified drop are investigated under the influences of various parameters such as the Prandtl number Pr, the Stefan number St, the Bond number, the Ohnesorge number Oh, and the density ratio of the solid to liquid phases qsl. We also consider the effects of the growth angle /gr (at the triple point) and the initial drop shape (in terms of the contact angle /0) on the solidification process. The most influent parameter is Bo whose increase in the range of 0.1–3.16 makes the drop more deformed with the tip shift linearly increasing with Bo. The tip shift also increases with an increase in /0. However, increasing Oh (from 0.001 to 0.316), St (from 0.01 to 1.0) or qsl (from 0.8 to 1.2) leads to a decrease in the tip shift. Concerning time for completing solidification, an increase in Bo, /gr or /0, or a decrease in St, or qsl results in an increase in the solidification time. The effects of these parameters on the drop height are also investigated. Ó 2018 Elsevier Ltd. All rights reserved. 1. Introduction Solid–liquid phase change of drops is of crucial importance because the process widely exists in nature and engineering phenomena. For instance, water drops freeze on leaves , airplane , or wind turbine blades . Metal drops solidify in deposition manufacturing  or in atomization . Molten semiconductor drops crystallized during falling or on cold substrates have been used for solar cell applications [6–9]. Accordingly, understanding the complex heat transfer and solidification phenomena is extremely important to advance the above-mentioned technologies. Concerning liquid drops solidifying on a horizontal plate, a vast number of works have been conducted. Experimentally, one can find this problem investigation in one of the pioneering works, Anderson et al. , in which the authors froze a water drop to demonstrate the accuracy of the proposed dynamic contact angle model at the tri-junction. An interesting feature observed from the experiment is the formation of an apex at the drop top after complete solidification. After then, a few experiments focusing ⇑ Corresponding author. E-mail address: firstname.lastname@example.org (T.V. Vu). https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.031 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved. on this singularity formed on frozen water drops have been done, e.g. [11–13]. Recently, this problem has been rapidly growing and getting much of attention. Jin et al.  froze a water drop sessile on a plate by lowering the plate temperature to a value below the water freezing point. The ice layer initially formed at the plate propagated to the drop top to complete the solidification. Similar works have been done in Zhang et al.  and Zhang et al. [16,17]. Focusing on not only water but also some other semiconductor materials (e.g. silicon, germanium and indium antimonide), Satunkin  reported the formation of conical solidified drops induced by volume expansion and growth angles at the trijunction. Basing on theoretical analysis supported by the experiments, Satunkin found that the growth angle at the tri-junction is almost constant except near the end of the solidification process. Another work concerning the crystallization of a molten silicon drop can be found in Itoh et al. . Theoretical studies on these problems can be found in [19–21], in which the authors mostly paid attention to the solidified drop profile and the temporal evolution of the solidification front assumed to be flat. Concerning numerical simulations of drops solidifying on a plate, a few studies have been performed by a boundary integral method , a finite element method , an enthalpy-based method , a volume of V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 fluid method [17,25], and a front-tracking method in our recent works [26–29]. However, in these works, the authors considered only drops attached to a horizontal wall. Related to drops solidifying on inclined surfaces, there have been many studies focusing on drops impacting, spreading and solidifying on the plates. For instance, Jin et al. [30,31] performed experiments for the freezing process of water on different inclined surfaces of different materials. To study the freezing process of asymmetric water drops, Ismail and Waghmare  tilted a cold plate and established the relationship between the asymmetry (in terms of contact angles) and the tip position and height of the frozen drops by scaling analysis and experiments. Wang and Matthys  experimentally investigated the heat transfer between the solidifying drops, of nikel and copper, and a metallic substrate tilted at an angle of 45°. Numerically, Zhang et al.  used a method of smoothed particle hydrodynamics to find the distribution of the temperature and movement of the solid–liquid interface of a drop impacting and solidifying on an inclined interface. Yao et al.  used the volume of fluid method implemented in an open source code OpenFOAM  to simulate the freezing process of a water drop. Some other numerical works can be found in [37,38]. Concerning drops on a vertical wall, Podgorski et al.  reported controlled experiments to show various shapes with remarkable temporal patterns of a water drop on a vertical plane. Smolka and SeGall  presented the formation of a fingering pattern, which tends to break up into drops, on the surface outside of a vertical cylinder. Li and Chen  experimentally formed frozen drops on a vertical wall and used the ultrasonic vibration to remove them from the wall. Numerically, Schwartz et al.  used a long-wave or lubrication approximation-based model to simulate three-dimensional unsteady motion of a drop on a vertical wall. Tilehboni et al.  used a lattice Boltzmann method to simulate a liquid drop moving on a vertical wall. However, the abovementioned works have not considered solidification heat transfer. Even though there have been many numerical studies on the phase change heat transfer of liquid drops, they mostly focused on the horizontal plate. The drop solidification on a vertical wall is rarely found. Accordingly, filling this gap is the main purpose of the present study because of its importance in academic and engineered applications [44–47]. We here use a two-dimensional front-tracking method, for tracking interfaces, combined with an interpolation technique, for dealing with the non-slip boundary on the solid surface, [26,48–50] to investigate the deformation with the solidification of a liquid drop attached to a vertical cold wall. Various parameters are investigated to reveal their effects on the solidification process. 2. Mathematical formulation and numerical method @ qu @ qu2 @ quv @p @ @u @ @u @ v þ þ l þ þ ¼ þ 2l @t @x @x @x @y @y @x @x @y ! þ qg x þ i ðFm þ Fr Þ Fig. 1. A two-dimensional liquid drop, with an initial center of mass C0(xC0, yC0), solidifies from a vertical cold wall. Cl(xCl, yCl) is the center of mass of the liquid part with e = yC0 yCl called ‘‘shift” of the liquid part. Has is the averaged height of the solidification front. @ qv @ quv @ qv 2 @p @ @u @ v @ @v þ þ þ ¼ þ l þ 2l @y @x @y @x @y @t @x @y @y ! þ qg y ½1 bðT T m Þ þ j ðFm þ Fr Þ @ðqC p TÞ @ðqC p TuÞ @ðqC p T v Þ þ þ @t @x @y Z @ k@T @ k@T þ þ q_ f dðx xf ÞdS ¼ @x @x @y @y f @u @u þ ¼0 @x @y ð1Þ ð2Þ ð3Þ ð4Þ Here, u = (u, v) is the velocity vector. p is the pressure, g = (gx, gy) is the acceleration due to gravity. T is the temperature. b is the ! ! Boussinesq coefficient. f denotes interface. i and j are respectively the unit vectors on the x and y axes. Fm is the momentum forcing term used for enforcing the non-slip boundary condition on the solid surface [48,49]. Fr is the interfacial tension force acting on the liquid–gas interface : Z Fr ¼ f rjdðx xf Þnf dS ð5Þ where r is the interfacial tension coefficient, j is the curvature, nf is the unit vector normal to the interface. d(x xf) is the Dirac delta function, which is zero everywhere except for a unit impulse at the interface xf. q_ is the heat flux at the solidification interface, given as q_ ¼ ks We consider a liquid drop that attaches to a vertical wall kept at constant temperature Tc (Fig. 1). The liquid has a fusion temperature Tm higher than the temperature of the wall, Tm > Tc, and thus at the beginning, a thin liquid layer at the wall changes into solid. As time progresses, this solid layer evolves to the top of the drop. Here, we consider only a two-dimensional drop. We assume that the gas and liquid phases are incompressible, immiscible and Newtonian. In addition, the thermal (thermal conductivity k and heat capacity Cp) and fluid (density q and viscosity m) properties are assumed constant in each phase. The one-fluid representation gives 303 @T @n kl s @T @n l with the subscripts s, l and g (when available) denoting solid, liquid and gas. In the present study, volume change induced by the density difference between the solid and liquid phases is also taken into account, and thus Eq. (4) is modified as  Z @u @ v 1 1 1 _ þ dðx xf ÞqdS ¼ @x @y Lh qs ql f ð6Þ where Lh is the latent heat of fusion. Initially, we assume the drop as a section of a sphere with a contact angle at the wall /0 and the apparent radius R ¼ ½3V 0 =ð4pÞ1=3 , where V0 is the initial volume of the drop. To simplify the problem, the temperature in the entire domain is set to T0 equal to the fusion temperature, T0 = Tm. In addition, a thin solid layer at Tc with a thickness of 0.02R is present at the plate at t = 0 . The boundary conditions are shown in Fig. 1. 304 V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 To solve this problem, we use the front-tracking method that has been extensively used in our previous works for the drop solidification on a plate with free and forced convection (e.g. [26,27,29,48]). Accordingly, we here briefly describe the method. The momentum and energy equations are approximated by a second-order central difference with a second-order predictorcorrector scheme for time integration. The discretized equations are solved on the fixed, rectangular staggered grid by MAC method . The interfaces separating different phases are approximated by a chain of points that moves on the fixed grid (Fig. 1). For the phase change interface, the points propagate with the velocity Vn calculated from the heat balance at this interface, i.e. _ qs Lh Þ. For the liquid–gas interface, its points are updated V n ¼ q=ð by the velocity interpolated from the velocities at the nearest fixed grid points. At the triple points, where these two interfaces meet, their positions are corrected by imposing a constant growth angle, /gr = constant [23,26]. In addition, to simplify the problem, we assume that the growth angles are the same at the upper and lower triple points. When the solidification process proceeds, the triple points form the solid–gas interface. Details of the numerical method can be found elsewhere, e.g. [26,27,29,48]. 3. Numerical parameters, grid refinement test and method validation 2 We scale the length, time and velocity by d = 2R, sc ¼ ql C pl d =kl and Uc = d/sc, respectively. Thus, the dimensionless time is s = t/sc with the dimensionless temperature defined as h ¼ ðT T c Þ= ðT m T c Þ. It is easy to prove that the following parameters govern the dynamics of the problem : Pr ¼ C pl ll ; kl St ¼ C ps ðT m T c Þ ; Lh 3 Ra ¼ g y bl ðT m T c Þd ml al ; ql g y d2 ; r l l ﬃ Oh ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ql d r h0 ¼ T0 Tc ; Tm Tc bgl ¼ bg ; bl ksl ¼ ks ; kl kg ; kl C psl ¼ kgl ¼ Bo ¼ qsl ¼ C ps ; C pl ð7Þ lg qg qs ; qgl ¼ ; lgl ¼ ql ql ll C pgl ¼ C pg C pl ð8Þ ð9Þ Here, Pr, St, Bo, Ra and Oh are respectively the Prandtl, Stefan, Bond, Rayleigh and Ohnesorge numbers. h0 is the initial dimensionless temperature. Other parameters include the ratios of the thermal expansion coefficients (bgl), densities (qsl and qgl), viscosities (mgl), thermal conductivities (ksl and kgl) and heat capacities (Cpsl and Cpgl). We have found that the variations of Ra in the range of 10–10,000 and bgl in the range of 0.0005–0.5 have almost no effect on the shape of the solidified drop as shown in Fig. 2. Thus, we keep these parameters fixed Ra = 1000, and bgl = 0.005. Similarly to our previous works , we here focus on the effects of Pr, St, Bo, Oh and the solid-to-liquid density ratio qsl with various initial liquid drop shape (i.e. various /0). Accordingly, other parameters are also kept fixed: h0 = 1.0, ksl = 1.0, kgl = 0.005, Cpsl = Cpgl = 1.0, and qgl = lgl = 0.05. The computational domain is chosen as W L = 2d 6d with the grid resolution of 256 768 based on the grid refinement study shown below. The values of the parameters investigated in this study correspond to drops with diameter in the order of a few millimeters (see the tables shown in [26,27] for reference). As demonstrated in our previous work for the case of solidification affected by a laminar crossing flow , the position of the drop has a very minor effect on the solidification process. Accordingly, in this study we skip the investigation of the drop location, and thus the initial position of the drop is also kept fixed, yC0 = 3.5d. Fig. 2. Effects of (a) the Rayleigh number and (b) the thermal expansion coefficient ratio of the gas to liquid phases on the solidified drop profile. Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, qgl = 0.9, h0 = 1.0, ksl = 1.0, kgl = 0.005, Cpsl = Cpgl = 1.0, qgl = lgl = 0.05, /gr = 0° and /0 = 105°. To verify the numerical method applied to the present problem, we perform a grid refinement study with Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. Two grid resolutions 128 384 and 256 768 are used. Fig. 3 shows a comparison of the drop profiles at different moments of solidification and the averaged distance from the solidification front to the wall (i.e. the averaged height Has defined below) computed by two grids. Fig. 3 indicates the results with complete agreement between two grids. Accordingly, we use 256 768 for the rest of the computations presented in this paper. The numerical method used in the present study has been extensively used in our previous works. The validations of the method were also carried out [26,48–50,53]. However, to ease reference, we here present the solidification result of a drop sessile on the horizontal plate and compare with the experiments of Zhang et al. . This choice of comparison is due to the lack of a detailed experiment on the drop solidification on a vertical wall. The initial volume of the water drop was 10 mL, and the plate temperature was set to 15.6 °C. Fig. 4 show the results yielded by our method in comparison with the results reported in  (see Figs. 4 and 6 in ). Concerning the conical angle at the top of the frozen drop, our computation yielded an angle of 140.62° whereas that observed in  was 140°, resulting in an error of 0.44%. Once again, reasonable agreement is found, supporting the accuracy of the method used in the present study. Other validation can be found elsewhere, e.g. [26,48–50]. 4. Results and discussion Fig. 5 illustrates the deformation of the solidifying drop against the vertical wall with Pr = 0.01, St = 0.1, Bo = 1.78, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. xn and yn in Fig. 5 and the following figures are the normalized coordinates xn ¼ x=d; yn ¼ ðy yC0 Þ=d ð10Þ This typical case with /0 = 105° corresponds to a drop attached to a hydrophobic cold wall [9,15] with purposes of, for example, reducing ice beds on walls  or the formation of ‘‘spherical” crystallized drops [6,7,9]. At s = 0.18 (Fig. 5a), the gravity starts to deform the drop whereas the solid–liquid interface parallel to the wall grows fast because of the large temperature difference between the wall and the liquid . This results in the variation of Has with a high slope at the beginning (Fig. 5f), where Has, the averaged height of the solidification front (Fig. 1), is defined as Pns Has ¼ i¼1 xfsi ns ð11Þ V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 305 Fig. 3. Grid refinement study with Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°: (a), (b), (c) drop profile at the different stages of solidification, and (d) temporal evolution of Has - the averaged height of all points on the solidification front away from the wall (i.e. the average of the x coordinates). Fig. 4. Comparisons between the present computation and the experiment of Zhang et al.  for a water drop solidifying on a horizontal plate: (a) drop profiles at different stages of solidification and (b) temporal variation of the drop volume Vd normalized by the initial volume V0. In Eq. (11), xfs and ns are the x coordinate of a point and the number of points representing the solid–liquid interface. Accordingly, the solidified drop height Hs is Has at the end of the solidification process. Because of the initially gravitational effect, the drop slightly deforms with a downward flow in the liquid and an upward flow in the surrounding gas, i.e., a counterclockwise circulation formed. Thus, the shift e of the liquid phase (i.e. the shift in the vertical position of the center of mass of the liquid phase) defined as e ¼ yC0 yCl ð12Þ gradually increases. At s = 0.48 (Fig. 5b), the drop becomes more deformed due to gravity, and the liquid is still falling whereas the solid–liquid interface continues moving to the top of the drop. This results in a continuous increase in the shift of the liquid part. At a later time (s = 0.9, Fig. 5c), the interfacial tension force tending to hold the drop in a circular shape is against the gravity force and pulls the liquid upwards. Thereby, a clockwise circulation is formed around the top of the drop. This leads to a decrease in the shift of the liquid part, as shown in Fig. 5f. The gravity then pulls down the liquid phase and forms a counterclockwise circulation again. This, falling and moving up, happens until the whole drop solidifies, resulting in a damping wave-like pattern of the liquid phase shift (Fig. 5f). An interesting feature here is that the movement of the solidification interface is almost independent of the damped oscillation of the liquid adhered to it (Has in Fig. 5f), and a cone forms at the top of the drop at the end of solidification because of volume expansion (Fig. 5e). It means that, like the drop sessile on the plate with free convection, the solidified drop attached to the vertical wall still presents a conical tip at the drop top, but the tip location shifts by distance et = 0.26d (where et is the shift of the solidified drop tip, i.e. the tip shift) because of gravity (see the definition in Fig. 5e). Next, we see how the solidified shape and some other interesting features of the drop are affected by some dimensionless parameters. 4.1. Effect of the Bond number Bo One of the most affecting parameters in the present problem is the Bond number, which is the ratio of the force induced by gravity to the interfacial tension force. Gravity tends to deform the drop during solidification. In contrast, the tension force tries to keep the drop as spherical as possible. Accordingly, as Bo < 1.0, the gravitational force plays a weak role, and thus the drop keeps almost spherical during solidification, as shown in the left frame of Fig. 6a. Increasing the value of Bo to a value higher than 1.0, i.e. 306 V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 Fig. 5. Evolution of the solidification front with the temperature and velocity fields. Has and e are defined in the text. The parameters are Pr = 0.01, St = 0.1, Bo = 1.78, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. The circles in (f) correspond to the times shown in (a)–(e). The velocity in (a) and all following figures is normalized by Uc. Fig. 6. Effect of Bo: (a) flow and temperature fields at s = 0.9 for Bo = 0.32 and 3.16, (b) solidified drop shape for various Bo at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to Bo. The dash-lines in (d) represent the corresponding fittings with ss = 0.0676Bo2 0.0445Bo + 2.4156 and et/d = 0.1723Bo 0.0149. Pr = 0.01, St = 0.1, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. Bo = 3.16 (right frame of Fig. 6a), the situation changes dramatically with large deformation of the solidifying drop. This results in the difference in the solidification front shape and in the temperature field around it. That is, increasing Bo from 0.32 to 3.16 makes the phase change front change from almost flat shape to a S-line shape. Thus, the thermal boundary layer around the solidliquid interface becomes inclined for the higher Bo. As a result, the solidified drop with large deformation is produced with a high V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 Bo, i.e. Bo = 3.16 (Fig. 6b). Fig. 6c confirms that increasing Bo results in the liquid part more deformed to the bottom, and thus increases its shift. Consequently, when Bo increases, the solidification process takes longer time to complete, and the solidified drop tip shifts more down to the bottom, as shown in Fig. 6d. The curve fittings indicate that the tip shift is linearly proportional to Bo (et/d 0.1723Bo 0.0149) whereas the solidification time increases with the square of Bo (ss 0.0676Bo2 0.0445Bo + 2.4156). The increase of the solidification time with Bo is also understandable from Fig. 6a in which the higher Bo (Bo = 3.16) corresponds to the smaller liquid-solid phase change area. As a result, it takes more time to complete the solidification process as Bo increases . However, the solidified drop height is almost independent of the variation of the Bond number in the range of 0.1–3.16. 4.2. Effect of the Prandtl number Pr For Pr = 0.01, at s = 0.36, the liquid part is going down due to gravity, no oscillation has happened before (left frame of Fig. 7a). Increasing Pr to a value 10 times higher than (i.e. Pr = 0.1), the liquid is also flowing down but with a higher velocity, resulting in larger deformation and thus in a higher shift e, as shown in Fig. 7c. Fig. 7c also indicates that the liquid is merely shifting for Pr = 0.01whereas Pr = 0.1 has made the liquid phase experience a little oscillation. However, this oscillation occurs only at the initial stage of solidification and damps soon after, whereas in the lower Pr cases, the oscillation happens later and lasts longer (Fig. 7c). This is understandable since higher Pr corresponds to higher viscous resistance that suppresses the oscillation of the drop as the solidification progresses. Consequently, the interface of the solidified drop becomes smoother (or less wavy) as the value of Pr increases from 0.01 to 1.0, as shown in Fig. 7b. However, an increase of Pr in the range of 0.01–1.0 just slightly decreases the solidified drop height with a slightly increase in et and thus makes the solidifica- 307 tion process finish a little bit sooner, as depicted in Fig. 7d. This is understandable since, as previously explained, an increase in e corresponding to an increase in Pr (Fig. 7b and c) tends to increase the solidification time (Fig. 6). However, as demonstrated in our previous work , increasing Pr also tends to decrease time for completing solidification. Accordingly, the combination of the both contrast effects results in a minor effect on the solidification time (Fig. 7d). 4.3. Effect of the Stefan number St Fig. 8 shows the effects of the Stefan number on the solidification process on a vertical wall with Pr = 0.01, Bo = 1.0, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. In comparison to St = 0.032 (left frame of Fig. 8a), the solid–liquid phase change front for St = 0.316 (right frame of Fig. 8a) advances further at s = 0.24 because of a lower latent heat released resulting a higher growth rate . Because of the fast solidification, St = 0.316 leads to a solidified drop with less deformation than St = 0.032 (Fig. 8b) and with a considerably reduction in the number of oscillation before complete solidification (Fig. 8c). As a result of less deformation, the solidified drop of St = 0.316 is more spherical than that of the lower one, thus St = 0.316 yields a drop with a lower height than St = 0.032 because of mass conservation. That is, increasing the Stefan number in the range of 0.01 to 1.0 results in a slight decrease in the height and tip shift of the solidified drop with a strongly decrease in the solidification time (ss 0.2853St0.949), as shown in Fig. 8d. 4.4. Effect of Ohnesorge number Oh Fig. 9a compares the drop profile, the temperature field and the velocity field at s = 0.36 for two Oh numbers: Oh = 0.00316 (left) and Oh = 0.0316 (right), in which the stronger flow field corresponds to the lower value of Oh. It seems to be contrast to our Fig. 7. Effect of Pr: (a) flow and temperature fields at s = 0.36 for Pr = 0.01 and 0.10, (b) solidified drop shape for various Pr at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to Pr. Bo = 1.0, St = 0.1, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. 308 V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 Fig. 8. Effect of St: (a) flow and temperature fields at s = 0.24 for St = 0.032 and 0.316, (b) solidified drop shape for various St at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to St. The dash-line in (d) represents the corresponding fitting with ss = 0.2853St0.949. Pr = 0.01, Bo = 1.0, Oh = 0.01, qgl = 0.9, /gr = 0° and /0 = 105°. Fig. 9. Effect of Oh: (a) flow and temperature fields at s = 0.36 for Oh = 0.00316 and 0.0316, (b) solidified drop shape for various Oh at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to Oh. Pr = 0.01, St = 0.1, Bo = 1.0, qgl = 0.9, /gr = 0° and /0 = 105°. V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 expectation since Oh is the ratio of the viscous force to the interfacial and inertial force. Thus, if Oh increased the drop should deform more. However, the tendency shown in Fig. 9a–c is understandable since we keep all parameters, especially Bo, unchanged except for Oh. Decreasing Oh corresponds to increasing the interfacial tension force and thus increasing the gravitational force since Bo is fixed. As a result of the gravitational effect, a decrease in Oh results in more deformation and in an increase in the number of oscillations before complete solidification (Fig. 9c). Consequently, the height and tip shift of the solidified drop slightly decrease as Oh increases. However, the variation of Oh in the range of 0.001–0.316 has a minor effect on the solidification time (Fig. 9d). 309 three-phase line arranges in such a way that the slopes of the liquid–gas and solid–gas interfaces (to the vertical line) decrease (to induce the increase in growth angle). Thereby, the solidified phase becomes smaller in the vertical direction but higher in the horizontal direction [18,26,27] (Fig. 10b). As a result of increasing the drop height (Fig. 10d), the time for completing solidification also linearly increases with the growth angle  (ss 1.0524/gr + 2.4428 where /gr is measured in radian). However, varying the growth angle has a very minor effect on the number of oscillations of the liquid phase during solidification, as shown in Fig. 10c, and on the tip shift (Fig. 10d). 4.6. Effect of the density ratio of the solid to liquid phases qsl 4.5. Effect of the growth angle /gr Our literature survey [18,23,27,28] indicates that the growth angle almost keeps constant during solidification and strongly affects the shape of the solidified drop. The value of the growth angle is dependent on the drop material (e.g., /gr = 0° for water, /gr = 12° for silicon, /gr = 14° for germanium, and /gr = 28° for indium antimonide [18,23]). For the solidifying drop against the vertical wall, we also observe such effect of the growth angle as shown in Fig. 10, in which the angle is varied in the range of 0–20°. The comparison in the velocity field for /gr = 4° and /gr = 20° indicates that the growth angle has a minor effect on the flow field, but the higher /gr produces a more conical solidified drop. As a result, the drop height of the solidified drop linearly increases with the growth angle (Hs/d 0.3556/gr + 0.8188 where the unit of /gr is radian). This can be understandable since at the triple point, three phases meet, and the interfacial tension energy balance locally appears with the solid–liquid front almost perpendicular to the liquid– gas interface. Accordingly, when the growth angle increases the It has been reported that the density difference between the liquid and solid plays an important role in the final product of the drop solidification when the drop size is in a range of a few millimeters [18,28]. After complete solidification, the solidified drop tends to possess a dimple if the solid-to-liquid density ratio is greater than unity whereas the conical top surface is present if the solid density is lower than the liquid. This is evidently shown in Fig. 11b where the solidified drop profiles for three density ratios qsl = 0.8 (expansion), 1.0 (no volume change) and 1.2 (shrinkage) are depicted. The reason is that, for small drops, the interfacial tension force can hold the drop as a section of a sphere, and other force, i.e. the gravity in this study, slightly deforms the drop. Then, in such situation, the volume change would impose its evident impact on the drop by expanding the drop for qsl < 1.0 or shrinking it for qsl > 1.0 during solidification. Such effect results in the flow away (qsl < 1.0) or towards (qsl > 1.0) the solid phase near the solid–liquid interface (Fig. 11a). During the final stages of solidification, the liquid drift away (qsl < 1.0) or towards (qsl > 1.0) the solidifying front leads to the formation of an apex Fig. 10. Effect of /gr: (a) flow and temperature fields at s = 0.9 for /gr = 4° and 16°, (b) solidified drop shape for various /gr at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to /gr (the linear fittings give Hs/d = 0.3556/gr + 0.8188, ss = 1.0524/gr + 2.4428 and et/d = 0.0566/gr + 0.14, with /gr measured in radian). Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, qgl = 0.9 and /0 = 105°. 310 V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 Fig. 11. Effect of qsl: (a) flow and temperature fields at s = 0.84 for qsl = 0.8 and 1.2, (b) solidified drop shape for various qsl at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to qgl. The dash-lines in (d) represent the corresponding fittings with Hs/d = 1.2617q2 3.4247q + 2.8788, ss = 2.3028q2 6.6865q + 6.6019 and et/d = –0.1776q + 0.3048, where q = qsl. Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, /gr = 0° and /0 = 105°. Fig. 12. Effect of /0: (a) flow and temperature fields at s = 0.3 for /0 = 60° and 120°, (b) solidified drop shape for various /0 at the end of solidification, (c) temporal variation of the shift of the liquid part and (d) variations of Hs, ss and et with respect to /0 (the fittings, i.e. dash-lines, give Hs/d = 0.3385/ + 0.1855, ss = 0.4679/3 1.3536/2 + 2.6273/ 0.7035 and et/d = 0.061/3 0.1442/2 + 0.0975/ + 0.0726, with / = /0 measured in radian). Pr = 0.01, St = 0.1, Bo = 1.0, Oh = 0.01, qsl = 0.9 and /gr = 0°. V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312 or dimple at the top of the drop [18,28]. Consequently, the drop is more deformed (Fig. 11c) with an increase in the height and in the solidification time as the density ratio is decreased as shown in Fig. 11d. If we define the tip of the solidified drop as the last solidification point or the farthest point from the wall, then we have the variation of the tip shift with respect to qsl, as shown in Fig. 11d, in which the tip shift linearly decreases with an increase in qsl (et/d = 0.1776qsl + 0.3048). This tendency is consistent with intuition that increasing qsl results in smaller solidified drops (Fig. 11b). 4.7. Effect of the contact angle /0 To conclude, we examine how the drop shape, in terms of the contact angle induced by the wall wettability, affects the process. Fig. 12a compares the drop behavior with high (left) and low (right) wettability. Concerning the solidifying interface, at this time (s = 0.3) it has a concave-up arc near the triple point for /0 = 60° whereas /0 = 120° produces a concave-down shape at this point since the solid–liquid front tends to be perpendicular to the liquid–gas front . However, when half of the drop is solidified, the solid–liquid front for /0 = 120° will behave similarly to /0 = 60°, that is, a concave-up arc near the triple point [13,28]. As a result of the high contact angle resulting in a high drop, the liquid part experiences more oscillation (Fig. 12c) with a more deformed solidified drop after complete solidification (Fig. 12b) when the contact angle increases. This leads to an increase in the solidified drop height as well as the solidification time. As a result, the location of the tip of the solidified drop shifts more to the bottom with the increase in /0 in the range of 45–135° (Fig. 12d). Fig. 12b also suggests that, large contact angles induced by drops on hydrophobic walls could cause the drops to easily fall off the wall after complete solidification because of gravity, and thus enhance the deicing process . 5. Conclusion We have presented a fully resolved two-dimensional numerical investigation of a liquid drop solidifying on a vertical cold wall by the front-tracking method to track the temporal movement of the interfaces. The non-slip velocity boundary condition on the solid surface is satisfied by a linear interpolation technique. The gravitational force results in an asymmetric drop whose conical tip induced by volume expansion is shifted down after complete solidification. We varied the Prandtl number Pr (in the range of 0.01– 1.0), the Stefan number St (in the range of 0.01–1.0), the Bond number Bo (in the range of 0.1–3.16), the Ohnesorge number Oh (in the range of 0.001–0.316), the density ratio of the solid to liquid phases qsl (in the range of 0.8–1.2), and the growth angle /gr (in the range of 0–20°) and the contact angle /0 at the wall (in the range of 45–135°) to show their influences on the shape, height and tip shift of the solidified drop and the solidification time. The numerical results show that the solidified drop shape is strongly influenced by Bo, St, Oh, qsl, /gr and /0. The location of the solidified drop tip shifts more to the bottom as we increase any one of Bo and /0 or decrease any one of St, Oh and qsl. For instance, the tip shift et varies with respect to Bo by et/d 0.1723Bo 0.0149 and with respect to qsl by et/d 0.1776qsl + 0.3048. The height Hs of the solidified drop strongly increases with an increase in /gr (Hs/d 0.3556/gr + 0.8188 with /gr measured in radian) and /0 (Hs/d 0.3385/0 + 0.1855 with /0 measured in radian) or with a decrease in qsl (Hs/d 1.2617q2 3.4247q + 2.8788). Concerning time ss for completing the solidification process, the most affecting parameters are Bo (it increases with Bo by ss 0.0676Bo2 0.0445Bo + 2.4156), St (it decreases with St by ss 0.2853St0.949), /gr (it increases with /gr in radian by ss 1.0524/gr + 2.4428), qsl 311 (it decreases with qsl by ss 2.3028q2sl 6.6865qsl + 6.6019) and /0 (it increases with /0 in radian by ss 0.4679/30 1.3536 /20 + 2.6273/0 0.7035). Even though our study is very detailed, many questions are still unresolved. For instance, the present calculations are merely twodimensional, and indeed for a ridge rather than a drop, and thus three-dimensional (3D) simulations would circumvent the present limitation. However, one of the most difficulties in 3D computations is to handle the growth angle (or the contact angles) at the tri-junction (i.e. three-phase line). In addition, the contact angles might be varied along the three-phase line because of the gravitational effect. The drop might slide on the wall before starting solidification, and thus the simulation requires the dynamic contact angle and nucleation models. In some cases, the drop might pinch off before completing solidification, and thus further investigations need to be done to reveal the parameter ranges in which the drop pinch-off happens. Acknowledgments This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.03-2017.01. The authors are grateful to Prof. John C. Wells at Ritsumeikan University, Japan for facilitating computing resources. Conflict of interest We have no conflict of interest to declare. References  L.C. Mayo, S.W. McCue, T.J. Moroney, W.A. Forster, D.M. Kempthorne, J.A. Belward, I.W. Turner, Simulating droplet motion on virtual leaf surfaces, Roy. Soc. Open Sci. 2 (2015) 140528.  Y. Cao, Z. Wu, Y. Su, Z. Xu, Aircraft flight characteristics in icing conditions, Prog. Aerosp. Sci. 74 (2015) 62–80.  N. Dalili, A. Edrisy, R. Carriveau, A review of surface engineering issues critical to wind turbine performance, Renew. Sust. Energy Rev. 13 (2009) 428–438.  C.H. Amon, K.S. Schmaltz, R. Merz, F.B. Prinz, Numerical and experimental investigation of interface bonding via substrate remelting of an impinging molten metal droplet, J. Heat Transfer 118 (1996) 164–172.  T.V. Vu, H. Takakura, J.C. Wells, T. Minemoto, Production of hollow spheres of eutectic tin-lead solder through a coaxial nozzle, J. Solid Mech. Mater. Eng. 4 (2010) 1530–1538.  T. Minemoto, H. Takakura, Fabrication of spherical silicon crystals by dropping method and their application to solar cells, Jpn. J. Appl. Phys. 46 (2007) 4016– 4020.  A.V. Hariharan, J. Ravi, Laser conversion of high purity silicon powder to densified granular forms, WO2008057483 A2, 2008. <http://www.google.com. na/patents/WO2008057483A2> (accessed June 28, 2017).  M.G. Tsoutsouva, T. Duffar, D. Chaussende, M. Kamguem, Undercooling measurement and nucleation study of silicon droplets on various substrates, J. Cryst. Growth. 451 (2016) 103–112.  H. Itoh, H. Okamura, C. Nakamura, T. Abe, M. Nakayama, R. Komatsu, Growth of spherical Si crystals on porous Si3N4 substrate that repels Si melt, J. Cryst. Growth. 401 (2014) 748–752.  D.M. Anderson, M.G. Worster, S.H. Davis, The case for a dynamic contact angle in containerless solidification, J. Cryst. Growth 163 (1996) 329–338.  O.R. Enríquez, Á.G. Marín, K.G. Winkels, J.H. Snoeijer, Freezing singularities in water drops, Phys. Fluids 24 (2012), 091102-091102-2.  J.H. Snoeijer, P. Brunet, Pointy ice-drops: How water freezes into a singular shape, Am. J. Phys. 80 (2012) 764–771.  A.G. Marin, O.R. Enriquez, P. Brunet, P. Colinet, J.H. Snoeijer, Universality of tip singularity formation in freezing water drops, Phys. Rev. Lett. 113 (2014) 054301.  Z. Jin, X. Cheng, Z. Yang, Experimental investigation of the successive freezing processes of water droplets on an ice surface, Int. J. Heat Mass Transfer 107 (2017) 906–915.  H. Zhang, Y. Zhao, R. Lv, C. Yang, Freezing of sessile water droplet for various contact angles, Int. J. Therm. Sci. 101 (2016) 59–67.  X. Zhang, X. Wu, J. Min, Freezing and melting of a sessile water droplet on a horizontal cold plate, Exp. Therm. Fluid Sci. 88 (2017) 1–7.  X. Zhang, X. Liu, X. Wu, J. Min, Simulation and experiment on supercooled sessile water droplet freezing with special attention to supercooling and volume expansion effects, Int. J. Heat Mass Transfer 127 (2018) 975–985. 312 V.N. Duy, T.V. Vu / International Journal of Heat and Mass Transfer 127 (2018) 302–312  G.A. Satunkin, Determination of growth angles, wetting angles, interfacial tensions and capillary constant values of melts, J. Cryst. Growth 255 (2003) 170–189.  A. Sanz, The crystallization of a molten sphere, J. Cryst. Growth 74 (1986) 642– 655.  M. Nauenberg, Theory and experiments on the ice–water front propagation in droplets freezing on a subzero surface, Eur. J. Phys. 37 (2016) 045102.  X. Zhang, X. Wu, J. Min, X. Liu, Modelling of sessile water droplet shape evolution during freezing with consideration of supercooling effect, Appl. Therm. Eng. 125 (2017) 644–651.  W.W. Schultz, M.G. Worster, D.M. Anderson, Solidifying sessile water droplets, in: Interactive Dynamics of Convection and Solidification, Kluwer Academic Publishers, 2001, pp. 209–226.  A. Virozub, I.G. Rasin, S. Brandon, Revisiting the constant growth angle: estimation and verification via rigorous thermal modeling, J. Cryst. Growth 310 (2008) 5416–5422.  G. Chaudhary, R. Li, Freezing of water droplets on solid surfaces: an experimental and numerical study, Exp. Therm. Fluid Sci. 57 (2014) 86–93.  X. Zhang, X. Liu, X.M. Wu, J.C. Min, Numerical simulations of freezing process of a sessile supercooled water droplet using Eulerian method, in: Proceedings of the 9th International Symposium on Heat Transfer. Beijing, Beijing, China, 2016.  T.V. Vu, G. Tryggvason, S. Homma, J.C. Wells, Numerical investigations of drop solidification on a cold plate in the presence of volume change, Int. J. Multiphase Flow 76 (2015) 73–85.  T.V. Vu, Fully resolved simulations of drop solidification under forced convection, Int. J. Heat Mass Transfer 122 (2018) 252–263.  T.V. Vu, C.T. Nguyen, D.T. Khanh, Direct numerical study of a molten metal drop solidifying on a cold plate with different wettability, Metals 8 (2018) 47.  T.V. Vu, Deformation and breakup of a pendant drop with solidification, Int. J. Heat Mass Transfer 122 (2018) 341–353.  Z. Jin, D. Sui, Z. Yang, The impact, freezing, and melting processes of a water droplet on an inclined cold surface, Int. J. Heat Mass Transfer 90 (2015) 439– 453.  Z. Jin, Z. Wang, D. Sui, Z. Yang, The impact and freezing processes of a water droplet on different inclined cold surfaces, Int. J. Heat Mass Transfer 97 (2016) 211–223.  M.F. Ismail, P.R. Waghmare, Universality in freezing of an asymmetric drop, Appl. Phys. Lett. 109 (2016) 234105.  G.-X. Wang, E.F. Matthys, Experimental determination of the interfacial heat transfer during cooling and solidification of molten metal droplets impacting on a metallic substrate: effect of roughness and superheat, Int. J. Heat Mass Transfer 45 (2002) 4967–4981.  M.Y. Zhang, H. Zhang, L.L. Zheng, Simulation of droplet spreading, splashing and solidification using smoothed particle hydrodynamics method, Int. J. Heat Mass Transfer 51 (2008) 3410–3419.  Y. Yao, C. Li, H. Zhang, R. Yang, Modelling the impact, spreading and freezing of a water droplet on horizontal and inclined superhydrophobic cooled surfaces, Appl. Surf. Sci. 419 (2017) 52–62.  C. Greenshields, OpenFOAM. <https://openfoam.org/> (accessed May 15, 2018).  H. Zhang, L.L. Zheng, V. Prasad, T.Y. Hou, A curvilinear level set formulation for highly deformable free surface problems with application to solidification, Numer. Heat Tr. B-Fund. 34 (1998) 1–30.  M. Raessi, J. Mostaghimi, Three-dimensional modelling of density variation due to phase change in complex free surface flows, Numer. Heat Tr. B-Fund. 47 (2005) 507–531.  T. Podgorski, J.-M. Flesselles, L. Limat, Corners, cusps, and pearls in running drops, Phys. Rev. Lett. 87 (2001) 036102.  L.B. Smolka, M. SeGall, Fingering instability down the outside of a vertical cylinder, Phys. Fluids 23 (2011) 092103.  D. Li, Z. Chen, Experimental study on instantaneously shedding frozen water droplets from cold vertical surface by ultrasonic vibration, Exp. Therm. Fluid Sci. 53 (2014) 17–25.  L.W. Schwartz, D. Roux, J.J. Cooper-White, On the shapes of droplets that are sliding on a vertical wall, Physica D 209 (2005) 236–244.  S.M. Tilehboni, E. Fattahi, H.H. Afrouzi, M. Farhadi, Numerical simulation of droplet detachment from solid walls under gravity force using lattice Boltzmann method, J. Mol. Liq. 212 (2015) 544–556.  Z. Liu, X. Zhang, H. Wang, S. Meng, S. Cheng, Influences of surface hydrophilicity on frost formation on a vertical cold plate under natural convection conditions, Exp. Therm. Fluid Sci. 31 (2007) 789–794.  M. Fang, S. Chandra, C.B. Park, Building three-dimensional objects by deposition of molten metal droplets, Rapid Prototyping J. 14 (2008) 44–52.  C.-H. Wang, H.-L. Tsai, Y.-C. Wu, W.-S. Hwang, Investigation of molten metal droplet deposition and solidification for 3D printing techniques, J. Micromech. Microeng. 26 (2016) 095012.  F. Wang, C. Liang, X. Zhang, Research of anti-frosting technology in refrigeration and air conditioning fields: a review, Renew. Sust. Energy Rev. 81 (2018) 707–722.  T.V. Vu, J.C. Wells, Numerical simulations of solidification around two tandemly-arranged circular cylinders under forced convection, Int. J. Multiphase Flow 89 (2017) 331–344.  T.V. Vu, A.V. Truong, N.T.B. Hoang, D.K. Tran, Numerical investigations of solidification around a circular cylinder under forced convection, J. Mech. Sci. Technol. 30 (2016) 5019–5028.  T.V. Vu, Three-phase computation of solidification in an open horizontal circular cylinder, Int. J. Heat Mass Transfer 111 (2017) 398–409.  G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759.  F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182.  T.V. Vu, G. Tryggvason, S. Homma, J.C. Wells, H. Takakura, A front-tracking method for three-phase computations of solidification with volume change, J. Chem. Eng. Jpn. 46 (2013) 726–731.