- Email: [email protected]
Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points
Accepted Manuscript
Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points ˇ Edin Berberovi´c, Markus Schremb, Zeljko Tukovi´c, Suad Jakirli´c, Cameron Tropea PII: DOI: Reference:
S0301-9322(17)30646-8 10.1016/j.ijmultiphaseflow.2017.11.005 IJMF 2676
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
29 August 2017 16 October 2017 2 November 2017
ˇ Please cite this article as: Edin Berberovi´c, Markus Schremb, Zeljko Tukovi´c, Suad Jakirli´c, Cameron Tropea, Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.11.005
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
ACCEPTED MANUSCRIPT
Highlights • Development of a new computational framework for modeling the solidification of supercooled liquid water. • Utilization of the combination of the phase-field front propagation approach with the thermal diffusion equations in the phases.
CR IP T
• Decoupling at the moving phase-interface by using immersed points to impose the melting/freezing temperature at the interface. • Robust and efficient implementation within the object-oriented framework in foam-extend, the extenR sion of OpenFOAM • Extension of the numerical model to include conjugate heat transfer with the neighboring, thermally coupled solid wall.
AN US
• The validation of the numerical model demonstrates that it has a good capability to compute the solidification in the absence of walls, or far from the wall.
• At the wall, the problem formulation imposes a singularity as a fundamental difficulty at the triple contact line, formed during solidification of liquid over a solid substrate, which is also confirmed.
AC
CE
PT
ED
M
• The numerical model can be used as a framework for investigations of the interaction of the solidifying supercooled water with the thermally coupled solid substrate, by inclusion of and combining with suitably formulated plausible models, which would eventually remove the singularity.
1
ACCEPTED MANUSCRIPT
Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points ˇ Edin Berberovi´ca,∗, Markus Schrembb , Zeljko Tukovi´cc , Suad Jakirli´cb,∗, Cameron Tropeab a Polytechnic faculty, University of Zenica, Zenica, Bosnia and Herzegovina of Fluid Mechanics and Aerodynamics, Technische Universit¨ at Darmstadt, Darmstadt, Germany c Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Zagreb, Croatia
CR IP T
b Institute
Abstract
AN US
Computational modeling of phase change due to solidification, besides correctly capturing the heat transfer, requires accurate evaluation of mass transfer at the phase-interface. In the present study a framework is developed for modeling the solidification of supercooled water using a phase-field approach for the propagation of the interface between ice and supercooled water. Energy equations in the solid and liquid phases are decoupled at the interface by using immersed points to impose the melting temperature as a moving boundary condition. The propagation of the interface is determined by the local energy balance across the interface, which is accounted for in the reconstructed interface points. The model is validated using the known theoretical solutions for the two-phase Stefan problem. An extension of the model to incorporate conjugate heat transfer with a neighboring solid wall is presented, establishing a framework that can be used for further modeling of the interaction of solidifying supercooled water on a solid substrate. Keywords: supercooled water, freezing process, phase-field method, thermal diffusion model, three-phase contact point
M
1. Introduction
15
20
25
ED
PT
CE
10
AC
5
Pure water freezes at a temperature of 273.15 K, but under certain conditions it can be cooled down to almost 40 K lower temperatures, while retaining its liquid ther- 30 modynamic state [1]. The process of subsequent solidification of supercooled water is of essential interest for diverse engineering applications, where solid surfaces are exposed to different environmental conditions, common examples being wind turbines, or electric power line sys- 35 tems. A familiar practical application where the freezing of supercooled droplets poses a serious problem is in aviation, where droplets of supercooled water levitate in clouds. Aircraft passing through such clouds collides with the droplets, triggering their abrupt solidification, which 40 can lead to the formation of ice layers on aircraft surfaces, within the engine, or on instruments. In the freezing of supercooled water several stages are identified [2], and the main influencing parameters include the degree of supercooling, the environmental conditions, 45 and the solid surface temperature and structure [3, 4, 5, 6]. Understanding the processes in solidifying supercooled water and their interaction with a solid surface is of practical significance for defining methodologies for their prevention and deicing. However, even in the absence of phase change, 50 the non-isothermal two-phase flow in a supercooled water ∗ Corresponding
authors Email addresses: [email protected] (Edin Berberovi´ c ), [email protected] (Suad Jakirli´ c) Preprint submitted to International Journal of Multiphase Flow
55
drop impacting a surface represents a challenge, due to its high unsteadiness and the small characteristic length and time scales [7], which may span several orders of magnitude. The situation is additionally complicated if one of the phases undergoes changes in its thermodynamic state. The small scales render direct and detailed experimental access practically impossible, which is why research relies extensively on numerical tools. Therefore, there is considerable motivation for devising computational methods for solving the solidification of supercooled water, which are thermally coupled with the solid surfaces and enable examination of the governing physical mechanisms of the interaction. In the commonly applied numerical approaches, the ice front is treated as a liquid-solid phase-interface which evolves in time and space driven by the phase change [8, 9, 10]. The propagation of the interface is represented using either Lagrangian-based interface tracking or Eulerian-based volume capturing approach [11], coupled with heat transfer and a model for the mass transfer due to the phase change. The interface tracking provides a sharp representation of the interface by definition, whereas the volume capturing relies on the homogeneous mixture model and yields a smeared interface. The widely used volume-of-fluid approach has become an attractive methodology, utilizing an algebraic scheme for the propagation of the interface [12, 13], which can be combined with the enthalpy-porosity technique, or some if its variants, to model mass transfer [14, 15]. In spite of the work conducted in recent years using different modeling and numerical techniques, there is still a need to improve November 4, 2017
ACCEPTED MANUSCRIPT
85
90
95
100
105
110
CR IP T
iso-surface by virtually connecting the tips of dendrites in a propagating dendritic cloud. The main idea was to combine the numerical solution of the heat diffusion equation in the liquid around this cloud at a larger, mesoscopic scale, with a local analytical solution for the growth of a single dendrite tip, relevant at smaller, microscopic scales. The analytical solution is matched with the computationally obtained temperature field at a predefined distance from the iso-surface. In a modified version, the model was also successfully used in [22] for the solidification of an undercooled binary alloy. The analytical solution, applied for the velocity of the iso-surface, is the so-called stagnant film solution and represents a modification of the Ivantsov solution for the growth of a dendrite in the form of an isothermal paraboloid [23, 24]. Matching the analytical with the numerical solution is achieved by supplying the analytical solution with the value for the supercooling ∆T , obtained from the numerical solution for the heat diffusion equation, at a small normal distance from the isosurface δf , called the stagnant film thickness. Although the mesoscopic model, as will be shown later, is restricted to only small supercoolings and small interface velocities, the phase-field approach is considered suitable to be used for the interface propagation. In the present study, a computational model for freezing of supercooled water is presented within the standard cellcentered finite-volume CFD framework [25]. The model utilizes a phase-field diffuse interface model to track the liquid-solid (ice-water) interface. Energy equations in liquid and solid phases are expressed in terms of temperature. The interface points are reconstructed, and the solution for the temperature is decoupled at the interface, by using immersed points, and imposing the melting temperature as a moving boundary condition. The velocity of propagation of the interface is determined by the local interfacial energy balance, given by the sum of the normal-to-theinterface heat fluxes at both sides. With the calculated velocity, the diffused interface is explicitly evolved in time and space. The capability of the model is confirmed by calculating two-phase solidification with known analytical solutions. The model is further extended to couple the conjugate heat transfer in the solid substrate, comprising a framework which can be used to incorporate further models for the interaction of the solidifying supercooled water with the underlying solid substrate. The model enables computation of the solidification of supercooled water as an isolated phase-change process. The proper coupling with conjugate heat transfer enables computation of the interaction with the solid substrate and quantification of its thermal influence on the icing outcome. The previous findings related to the difficulties in the singularity at the contact line can be tested within the finite-volume framework for continuum mechanics in conjunction with the phase-field model, and the model can be used as a tool for further modeling attempts to overcome this difficulty. It should be noted that the model is aimed at, but not restricted to, computing the icing of
AN US
M
80
ED
75
PT
70
CE
65
AC
60
capabilities, especially if the conjugate heat transfer problem associated with the substrate is of importance. Investigations of cases of supercooled water freezing in115 contact with a solid substrate, revealed an interesting observation of a very thin ice film developing and traveling over the substrate surface, prior to the dendritic solidification [3, 4, 6]. The velocity of the tip of this ice layer was determined to be higher than that of dendritic solidi-120 fication in free liquid, depending strongly on the degree of supercooling and the substrate material. A similar phenomenon was studied in [16] by investigating the development of a molten contact line, which forms when a droplet of microcrystalline wax spreads and solid-125 ifies on a cold solid substrate, made of the same material. It was found that the advancing molten contact line approaches an apparent dynamic contact angle, which was modeled as a function of the properties of the substrate material and thermal conditions. This was used in a re-130 lated study in [17], where an attempt was made to compute the spreading and solidification of the molten material over a cold, solid substrate made of the same material. It was shown however, that the conventional continuum formulation for the problem has no meaningful solution. The135 contact angle of the solidified material film at the substrate is determined by the heat flux at the contact line, where the heat flux in the mathematical model has a singularity. The singularity is attributed to a discontinuity of the boundary condition for the energy equation at the140 contact line, which is also present in the case of pure heat transfer in a wedge geometry with sharp corners [18]. An attempt was made in [17] to remove this singularity, by estimating an ad-hoc cut-off length scale for the unknown physical mechanism, which would limit the heat flux and145 remove the singularity at the contact line, in order to yield agreement with the experimental data for a specific material. However, the length scale was obtained artificially, it was not physically justified and not universal, but rather case-specific, thus not providing insight into the underly-150 ing physical mechanism. In [19] an extended systematic computational study was conducted in order to overcome the singularity. A finite element method was used with several independent approaches: different forms of the predefined interface shape;155 different contact angles; imposing an overall heat transfer coefficient as a boundary condition at the interface; allowing temperature change along the interface; and using an artificial resistive layer at the substrate surface. However, again, none of these approaches solves the problem. It160 should be noted that the problem in [16, 17, 19] is simpler, since the solidifying material is the same as the substrate, compared to freezing of supercooled water, which is in contact with the substrate of different kind. In the present study, the phase-field approach is used165 for tracking the ice-water interface. The idea for using the phase-field is borrowed from [20, 21], where this approach was used to simulate the solidification of a supercooled melt. The interface was defined as the phase-field
3
ACCEPTED MANUSCRIPT
175
180
185
supercooled water, other fluids may be used as well. The205 local interfacial energy balance is applied by utilizing a quadratic interpolation algorithm for the temperature in probe points located at small distances from both sides from the reconstructed interface, thus allowing the effects from the surrounding temperature field in the evaluation210 of the local heat fluxes normal to the interface to be included. With regard to the methodology and physics, to the authors’ knowledge such a computational model does not exist in the present literature. Although it might be relevant only in regions of extremely small scales and ex-215 tremely large interface curvatures, the possibility to include the capillary undercooling is also readily available in the model, by simply replacing the melting temperature with the curvature-dependent one in the moving boundary condition. 2. Description of the model
100,00
ED
M
195
Ivantsov -5 stgf, δf=10 m -4 -3
ρl cp,l
-2
0,1
1 ∆T, K
∂T = ∇ · (kl ∇T ), ∂t
225
230
10 235
Figure 1: The dependence of the predicted front propagation velocity on the supercooling at various values for the stagnant-film thickness.
For small supercoolings, say ∆T < 1 K, the theoretical velocity can be predicted by using different values240 of the stagnant-film thickness, in which the temperature drop is a different portion of the global supercooling, as shown in Fig. 1. For higher values of the supercooling, the stagnant-film thickness must be lowered under the (larger) 4
(1)
∂T = ∇ · (ks ∇T ), (2) ∂t where the subscripts l and s refer to liquid (water) and the solid (ice), respectively. The symbols ρ, cp and k denote the density, specific heat and conductivity, respectively, and T is the temperature. The temperature at the interface between the solid and the liquid phase is the melting temperature Tm , and the local energy balance in the normal direction to the interface, for the two-phase moving boundary problem with ρs < ρl , is given by [26] ρs cp,s
ρs Lvn = ks ∇Ts − kl ∇Tl ,
stgf, δf=10 m
AC
0,01
stgf, δf=10 m
CE
v, cm/s
1,00
PT
stgf, δf=10 m
200
2.1. The basic phase-field model for solidification In the solidification of supercooled water, in the absence of bulk motion of the liquid water, the energy diffusion equations in the solid and liquid read
AN US
190
The aforementioned mesoscopic model uses the local temperature drop near the interface to determine its propagation velocity, and is applicable to small values of the global supercooling and small velocities. To illustrate this, the dependence of the predicted front propagation velocity on the supercooling and the stagnant film thickness220 is shown in Fig. 1. The values for the velocity were obtained by numerically inverting the analytical expression from the stagnant-film model coupled with the selection criterion, as described in [21], and using different values for the stagnant-film thickness, and constant thermophysical properties of water at 273.15 K.
mesoscopic scale down to microscopic scales, and the computationally obtained temperature field would have to be resolved on much smaller scales, which is what the mesoscopic model actually tries to avoid. In a recent study [3], utilizing a Hele-Shaw cell to investigate two-dimensional freezing of supercooled water in contact with the underlying solid substrate, the velocity of the ice-layer tip at the wall was found to be much higher than that of the dendritic front far from the wall, thus rendering the analytical stagnant-film solution for the solidification unsuitable. In spite of that, the interface propagation utilizing the phase-field approach provides a basis, that can be coupled to other phase change models and used to calculate the freezing.
CR IP T
170
(3)
where L is the latent heat of melting and the temperature gradients are evaluated on both sides of the interface. The capillary undercooling for curved interfaces can be accounted for with the Gibbs-Thomson relation, which has effects only for nanoparticles with very large curvatures [27] and is neglected here. For most materials the density of the liquid phase is usually smaller than that of the solid one, water being an exception in which the density of the solid phase is lower than that of the liquid one. In such a case, the motion of the phase-change interface is also affected by the volume expansion upon solidification. However, its relative influence on the freezing front velocity is determined by (ρwater /ρice − 1) [26], which has a relatively low amount of about 9% at T = 0◦ C, and decreases with decreasing temperature, since the liquid density decreases, while the ice density increases with decreasing temperature. Therefore, the volume expansion is of minor importance and will be neglected in order to simplify the analysis. In order to propagate the interface with the obtained normal velocity, the phase-field approach from [28] is used
ACCEPTED MANUSCRIPT
in the present study. The sharp interface between the solid and the liquid is represented by a thin transitional region, across which the indicator φ smoothly changes from the value of 1 to 0. The distribution of φ is obtained by solving275 the equation h = b ∇2 φ −
∂φ + vn · ∇φ = ∂t φ(1−φ)(1−2φ) (w/2)2
i ∇φ − |∇φ|∇ · − |∇φ| . (4)
The change of the phase-field variable φ from 1 to 0 in the diffused interface region is initially defined as the hyperbolic tangent profile n i 1h dist 1 − tanh , (5) φn = 2 w
obtained, the temperature gradients are calculated in the normal direction on both sides of the interface. The propagation velocity of the interface is determined from Eq. (3), and the new distribution of the phase-field is calculated explicitly from Eq. (4). The stabilization factor b in Eq. (4) influences the accuracy of the phase-field evolution and must be set using the stability limit blim that depends on the time step size and the mesh resolution [22, 28, 29]. In the present model, the time step size is self-adjusted during the simulation and the parameter b is set as
CR IP T
245
b = 0.05blim ,
blim =
M (6)
The term on the right-hand-side of Eq. (4) represents285 a stabilization term which acts towards ensuring the hyperbolic tangent profile of the phase-field, as discussed in [28, 29]. It is derived by decomposing the normal interface velocity into parts which are independent of and proportional to the interface curvature. This allows the curvature290 to be expressed as a function of the phase-field, which is given by the first two terms in the square brackets, while the last term in the square brackets is the counter-term that cancels out curvature-driven motion at leading order for the interface of finite thickness. In the case of a flat295 interface, the third term in the square brackets will vanish and the first two terms will yield the hyperbolic tangent profile in the stationary solution. During interface motion in absence of external advection, the mutual effects of the left- and right-hand sides of Eq. (4) relax the phase field300 distribution to a hyperbolic tangent profile across the interface and act towards sustaining the profile during the interface motion. The basic model consists of the equations (1)–(4). The solution algorithm starts by solving the energy equations305 (1)–(2). After the solution for the temperature has been
260
265
270
AC
255
CE
PT
250
∇φ . |∇φ|
ED
n=−
5
(8)
where the index t refers to the value in the current time step, index t − ∆t refers to the value from the previous time step, and ndim takes into account the dimensionality of the problem, being equal to 2 for two-dimensional and to 3 for three-dimensional cases, respectively. The maximum Courant number is based on the interface propagation velocity and is calculated as |vn |∆t . (9) Comax = max |d|
AN US
with ndist being the normal distance from the middle of the profile at which φ = 0.5. The characteristic width w, required to achieve sufficient accuracy in representing the hyperbolic transition, depends on the resolution of the finite-volume numerical mesh, and is calculated as w ≥ √ 2|d|, where d is the spacing between two adjacent cell centres in the numerical mesh. It corresponds to the cellsize ∆x in the case of a uniform hexagonal mesh. Thus, the diffused interface is spread across at least three cells at both sides of the iso-surface φ = 0.5. The solid phase is identified by 0.5 ≤ φ ≤ 1, whereas the liquid phase is determined by 0 ≤ φ < 0.5. It should be noted that the material properties in all equations are set according to this distinction. Using the phase-field indicator, the direction of the unit normal to the interface is evaluated280 from the expression
|d|2 1 [1 − 2Comax,t−∆t ] 2ndim ∆tt−∆t
(7)
For the solution of Eq. (4), the velocity should be calculated at all cell-centres. However in practice, it is enough to calculate the velocity only in cells within the transitional band between φmin = 0.01 and φmax = 0.99 [22]. The value 0.05 used in the formulation of the parameter b is smaller than the one recommended in [29] for similar reasons explained in [22], since this recommendation was obtained for a uniform normal interface velocity. In the present work, the problem is much more complex, since the interface deforms due to spatially and temporally nonuniform velocity. The influence of b on the solution and the interface shape was quantified in [22] by varying the ratio b/blim in a wide range. The variations of the width of the diffused interface over the range 0.05 < b/blim < 0.4 were found to be around 2.5%, but it was observed that a value for b higher than 0.05 constrains the interface shape too much. The present model is implemented in the open-source software foam-extend [30], an extended version of the R widely used software OpenFOAM [31]. The divergence and gradient terms in the transport equations are discretized using the Gaussian theorem, and temporal integration is Euler implicit. The new distribution for the phase-field is calculated explicitly, using old-time values for φ. The present model is similar to that in [22]; however, there are some major differences: instead of the analytical solution, the explicit energy balance at the interface is used to model mass transfer; transport equations for heat diffusion in both phases are used, instead of a conservation equation for the solute in the liquid phase; the interface
ACCEPTED MANUSCRIPT
310
points are calculated for the cells at both sides of the interface; the heat diffusion equations are discretized and solved as a single set of algebraic equations, with a proper decoupling at the interface, provided by adding the contributions to the coefficients in the solution matrix. These details are explained in the following subsections.
where the indices 1 and 2 indicate the edge start and end points, respectively. Thus, a set of the edges is detected that is connected and crossed by the interface, as shown in Fig. 2. The position vectors of the points where the interface crosses each mesh edge are calculated from xint,i = x1 +
2.2. Reconstruction of the interface points
320
xp = xP − δint ntri ,
with the same meaning of indices 1 and 2 as above. In order to accelerate the searching process, the width of the local transition region in which the search is confined is estimated as (1 − φint )φmin , (13) rmax = 1.25w ln φint (1 − φmin )
where the factor 1.25 is used to slightly extend the analytical value obtained by using the inverse of the hyperbolic tangent profile, Eq. (5). Based on that, a set of the nearest iso-surface points xint,i to each cell-centre xP within the band is found by calculating the distances between the cell-centre and the interface points lying within rmax
AN US
(10)
where xP are the position vectors of the cell-centres within the [φmin , φmax ] transition band, and ntri is the local unit normal to the interface. To calculate xp , the appropriate distances δint = xP − xp from all the cell-centres within the transition band to the corresponding points at the interface in the normal direction must be determined.
M
P
d int
fmin
PT
p
ED
ntri
interface
dP i = |xint,i − xP | < rmax .
CE
x int,1
x int,2
fmax
fint 325
Figure 2: Illustration of the algorithm for finding the points at the interface. The identified mesh edges crossed by the iso-surface are denoted by dashed lines.
The algorithm for finding the points at the interface,330 defined by φint = 0.5, is based on detecting the intersections of the mesh-edges with the interface. First the values for the phase-field variable φ stored at cell-centres are interpolated to the nodes of the cells in the mesh by using the least squares volume-to-point interpolation described335 in [32]. The edges of the mesh intersecting the interface, are extracted using the condition (φ1 − φint ) · (φ2 − φint ) ≤ 0,
(11) 6
(14)
From the set of the nearest interface points, the three nearest interface points are extracted with the smallest distances dP i from the cell-centre of interest, xP . Using the three interface points, a local unit normal to the interface is determined from the vector products xint,3 − xint,1 xint,2 − xint,1 × , (15) ntri = |xint,2 − xint,1 | |xint,3 − xint,1 |
where the indices 1, 2 and 3 indicate the three iso-surface points with the smallest distances dP i . Finally, the distance between the cell-centre of interest to the corresponding interface point is calculated by projecting the nearest distance dP i to the direction of the local unit normal, mapped to the centre of the cell P δint = |(xint,1 − xP ) · ntri |.
AC
315
(12)
CR IP T
The energy balance at the interface must be applied to calculate the velocity in Eq. (3) for all cells within the band 0.01 ≤ φ ≤ 0.99. Therefore, the corresponding interface points for all cells must be found in order to calculate the heat fluxes via the temperature gradients at both sides of the interface. During calculation, the interface points and the temperature gradients are mapped to all cells within the interfacial band. Referring to Fig. 2, the locations of the points p at the interface are determined from
φint − φ1 (x2 − x1 ), φ2 − φ1
(16)
If the cell P is in the solid, φP ≥ 0.5, the value δint is multiplied by -1, thus representing the signed normal distance from each cell-centre P within the band to the interface, being negative for the cells in the solid phase and positive for the cells in the liquid phase. The described procedure is performed for each cell within the predefined band φmin ≤ φ ≤ φmax on both sides of φint and the corresponding interface point is calculated from Eq. (10) for each of those cells. Currently, the described model is two-dimensional, and therefore only three points are necessary to determine the local iso-surface segment within the cell, according to R and in foam-extend, the mesh Eq. (15). In OpenFoam , is always three-dimensional, having at least one layer of cells. Thus, the mesh has a ’thickness’ of one cell, and simulations in one or two dimensions are enabled by suppressing the boundary conditions in the other dimensions. Therefore, the third required point xint,3 always lies at the edge at the opposite side of the cell.
ACCEPTED MANUSCRIPT
(17)
M
P
N
365
d
immersed points
ip N
x ip
x
fint
d
xP
Figure 3: Immersed points used in the discretization of the heat fluxes in the energy equations for setting the melting temperature Tm . The gray colored cells are identified as ’solid’ and the white colored cell as a ’liquid’ cell.
For a one-dimensional geometry, the change of φ around the ghost cell-face between the adjacent cell-centres P and N can be approximated as linear and, considering for example the cell P, the distance x can be calculated as φP − φint x= 1− |d|. (18) φP − φN However, for a two-dimensional geometry, the change of φ around the ghost cell-faces between adjacent cells can be arbitrary and therefore, x is approximated using a cubic interpolation a0 + a1 x + a2 x2 + a3 x3 = 0,
(20)
d · ∇φP , |d|
(21)
a2 =
3(φN − φP ) − d · (∇φN + 2∇φP ) , |d|2
(22)
a3 =
d · (∇φN + ∇φP ) − 2(φN − φP ) . |d|3
(23)
CR IP T
a1 =
The distance between the cell of interest P and the immersed point ip is then |xP − xip |, and similarly for the neighboring solid cell N |xN − xip |. In order to account for the melting temperature in the immersed points, referring to the upper two cells P and N in Fig. 3, the heat flux through the immersed cell-faces is discretized in the energy equations as kP
TP − Tm Sif , |xP − xip |
(24)
for the flux entering the cell P from the left, and
immersed faces Sif ip
a0 = φP − φint ,
kN
ghost faces
interface
with the following coefficients
AN US
xip
d = xN + x . |d|
ED
355
PT
350
CE
345
2.3. Decoupling of the energy equations For the discretization of the energy equations (1)–(2) in the liquid and the solid phase, an approach utilizing immersed points ip and immersed cell-faces is used. The ’liquid’ and ’solid’ cells are identified as those for which φ < φint and φ ≥ φint , respectively. The melting temperature Tm is explicitly set in the immersed points, which are mapped to the intersected cell-faces. The intersected (referred to as ghost) cell-faces are those, for which one adjacent cell sharing that face is identified as liquid, and the other as solid (faces with dashed edges in Fig. 2). The term ghost face originates from the fact that the discretized heat flux through such a cell-face is eliminated in the dis-360 cretization of the energy equation and replaced by the corresponding heat flux through the immersed faces. This enables to discretize the heat fluxes entering the interfacecrossed cells from the liquid-solid interface; hence, the immersed points have the purpose of acting as an immersed moving boundary at which Tm is set. The position vectors of the immersed points, shown in Fig. 3, are determined from
AC
340
(19) 7
Tm − TN Sif , |xip − xN |
(25)
for the flux entering the cell N from the right, with Sif being the area of the immersed cell-face equal to that of the ghost cell-face between the cells P and N. Thus, the corresponding contribution to the diagonal coefficient in the solution matrix for the temperature for the ’liquid’ cell P is kP Sif /|xP − xip |, and the source-term contribution is kP Tm Sif /|xP − xip |, and similar for the ’solid’ cell N. In the calculation of the velocity, it is important to evaluate correctly the temperature gradients at the interface, for which a more accurate approximation is required than a simple interpolation of the gradients from the surrounding cells. For a one-dimensional geometry, the temperature gradients at the interface could be calculated directly, by using the calculated temperatures from the cellcentres adjacent to the interface and their distances to the interface points. However, if the geometry is not onedimensional, this will not be accurate and the gradients will not be smooth enough along the interface, as shown in [33]. Therefore, for calculation of the temperature gradients at the interface in the present study, an additional set of points, referred to as probe points pp, is placed in the normal direction to the interface on both sides and at the same distance δpp . With reference to Fig. 4, the locations of the probe points for a cell P within the interfacial band are calculated xpp,l = xp + δpp nP ,
(26)
xpp,s = xp − δpp nP ,
(27)
ACCEPTED MANUSCRIPT
for the liquid and solid side respectively, with nP being the unit normal calculated using Eq. (6). The probe points are used to interpolate the temperature from the surrounding cell-centres and to calculate explicitly the normal-to-the-interface temperature gradients with the following expressions ∇Tl =
Tpp,l − Tm n, δpp
(28)
∇Ts =
Tm − Tpp,s n, δpp
(29)
weighted least squares method TP 0 − Tm c0 c1 ... c2 = [(XT ·W·X)−1 ·XT ·W]· TP j − Tm c3 ... TP n − Tm c4
L
interface
S
S
L
L
x pp,l
S
S
S
L
S
S
S
ppl p pps
xp
x pp ,s
dpp
390
395
M
dpp
ED
p – interface point, ppl/s – probe points in liquid/solid L/S – liquid/solid interpolation cells stencils
400
Figure 4: Probe points near the interface used for the calculation of the temperature gradients at the interface. The gray colored cells comprise the interpolation stencil for the probe point at the solid side, and the white colored cells are for the liquid side.
385
PT
If the temperature gradients at the interface would be calculated directly using the temperatures in the adjacent cell-centres, the region around the interface should be fully thermally resolved in order to obtain a good approximation for the temperature gradients. The thickness of the thermal diffusion layer in the normal direction to the interface is estimated as δt ∼ αl /vn , with αl being the thermal diffusivity of the liquid, and vn the magnitude of the interface velocity. In order to estimate the required spatial resolution, the velocity vn must be known. For demonstration, the velocity can be related to the thermal conditions from the steady state energy balance using a coordinate system positioned at, and moving with the interface with the constant velocity equal to vn . For the case of a onedimensional one-phase solidification of a supercooled liquid, shown schematically in Fig. 5, the energy equation in the liquid has the form
The distance δpp , as the input parameter, should be chosen such that the probe points are close enough to the interface to obtain a good approximation of the temperature gradients, but not too close in order to avoid numerical instabilities. The distance can be set to a portion of the cell size √ δpp = C∆x, where C ≤ 0.5 for one dimension, or C ≤ 0.5 2 for two-dimensions. The interpolation of the temperature in the probe points is performed using a second-order quadratic spatial function, defining the distribution of the temperature T over a region around the interface point p, which has the form
CE
380
3. Results and discussion
AC
375
which is first multiplied by XT · W from the left side, and afterwards again by (XT · W · X)−1 from the left side. The weights in the diagonal weighting matrix W are calculated based on the distance R of each cell-centre in the interpolation stencil to the interface point p, from the expression 0.5{1 + cos[(πR)/(1.1Rmax )]}/lref , where Rmax is the maximum distance from the cell-centres in the interpolation stencil, and the reference length is taken as lref = 2∆x. The interpolation stencils at the liquid and solid sides consist of two layers of cells closest to the interface point, although additional layers may be added. Upon obtaining the temperatures in the probe points Tpp from Eq. (30), the gradients at the interface are explicitly calculated using equations (28)–(29) and substituted to Eq. (3) to calculate the velocity, mapped to each cell within the diffused interface band.
AN US
L
, (31)
where X is the n × 5 matrix whose row j consists of the following entries: (xP j − xp ), (yP j − yp ), (xP j − xp )(yP j − yp ), (xP j − xp )2 , (yP j − yp )2 , and n is the number of the neighboring cells used in the interpolation stencil. Equation (31) follows from TP 0 − Tm c0 c1 ... T − T c (32) = X· m , 2 Pj c3 ... TP n − Tm c4
at the liquid and solid sides respectively. The distance to the probe points δpp = |xp − xpp,l/s | is equal for the liquid and the solid side of the interface.
L
CR IP T
370
T (x, y) = Tm + c0 (x − xp ) + c1 (y − yp )
+c2 (x − xp )(y − yp ) + c3 (x − xp )2 + c4 (y − yp )2 (30) where Tm is the melting temperature at the interface point p, and (xp , yp ) are the coordinates of that point. The unknown coefficients ci can be expressed in terms of the temperature TP j in the surrounding cell-centres using
−vn 8
∂T ∂2T = αl 2 . ∂x ∂x
(33)
ACCEPTED MANUSCRIPT
274 theory
273
Tm
-3 ∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
-vn
T, K
Tf vn
x
271 270 269 268
Figure 5: Sketch for the one-dimensional energy balance.
267 0
Integrating Eq. (33) once, the distribution of the temperature gradient in the normal direction to the interface is obtained ∂T vn ∂T = exp − x , (34) ∂x ∂x 0 αl
274 273
271
0,02
AN US
T, K
272
CR IP T
ice
272 water
where (∂T /∂x)0 is the temperature gradient at the interface x = 0. Integrating again, the temperature distribution in the normal direction is obtained ∂T αl vn T = T∞ − exp − x , (35) ∂x 0 vn αl
ED
M
420
425
0,08
0,1
theory -3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
270
267 0
In order to calculate the interface temperature gradient, a one-dimensional numerical setup was used, with one boundary declared as the interface and having a fixed temperature equal to Tm = 273.15 K, and at the opposite boundary the temperature T∞ = 267.15 K is set, yielding430 a global supercooling of 6 K. The velocity in the domain is constant and equal to 0.01 m/s, which corresponds to the experimentally obtained value in [34] for a supercooling of 6 K. The analytical value for the temperature gradient at 435 5 the interface equals −( ∂T ∂x )0 = 4.65 · 10 K/m, at x = 0. Equation (33) is solved using the standard finite-volume discretization in the software foam-extend, with the following values for the properties of water at ∆T = 6 K: kl = 0.5487, W/(mK), ρl = 999.08, kg/m3 and cpl = 4246 J/(kgK). With these values, the thermal diffusivity of water is αl = 1.29 · 10−7 m2 /s, and the estimated thickness of the thermal diffusion layer is δt = 12.9 · 10−5 m. The computed results for the temperature and the temperature gradient in the normal direction to the interface are shown in Fig. 6 and Fig. 7, where the results at the bottom of the figures are the zoomed profiles near the interface. The obtained profiles for T (x) and ∂T ∂x (x) approach the analytical solutions, Eq. (35) and Eq. (34), respectively,
CE
PT
415
0,06
268
AC
410
x, mm
269
where T∞ is the far-field temperature at x → ∞. Using ∆T = T − T∞ , an analytical relation is obtained between the global supercooling and the temperature gradient and the front propagation velocity at the interface (x = 0) ∂T ∆T vn − = . (36) ∂x 0 αl 405
0,04
9
0,01 x, mm
0,02
Figure 6: The computed temperature in the normal direction to the interface: temperature profile (top) is zoomed near the interface (bottom).
as the thermal diffusion layer in the normal direction to the interface becomes better resolved. From the results in Fig. 7 it follows that the thermal diffusion layer has to be resolved with more than 20 computational cells in the normal direction to the interface. At higher velocities, or higher degrees of supercooling, the thickness of the thermal diffusion layer will decrease, and the cell size required to resolve this layer would have to be much smaller. For a geometry which is not one-dimensional and at higher velocities, this leads to a significant increase of the computational effort. However, when the temperature gradients are calculated with the previously described procedure using the interpolated temperatures to the probe points, the required mesh resolution is expected to be lower. In order to test the capabilities of the model using the least square interpolation of the temperature in the probe points, the onedimensional two-phase Stefan problem is used. The liquid is initially supercooled and, at the initial time, the temperature at one boundary suddenly changes to Tm , triggering solidification, upon which the interface starts moving. The problem is described by equations (1)–(3), and the analyt-
∇T, K/m
ACCEPTED MANUSCRIPT
√ α vn = λ √ s . t
0 theory -3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
-1e+05 -2e+05
440
-3e+05
(41)
In the numerical setup, the length of the mesh in the horizontal direction, normal to the interface, is 10 mm, with an initial supercooling of 20 K obtained by specifying fixed temperatures Ts0 = Tl0 = 253.15 K at the opposite boundaries. The properties of water at the given supercooling, and that of ice, are listed in Table 1, and the latent heat of melting is L = 333000 J/kg.
CR IP T
-4e+05 Table 1: Thermophysical properties of water at ∆T = 20 K and ice.
0,02
0,04
x, mm
0,06
0,08
0,1
water ice
0 theory
-1e+05 -2e+05
445
-3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
450
-3e+05 -4e+05
455
-5e+05 0
0,02
M
0,01 x, mm
PT
ED
Figure 7: The computed temperature gradient in the normal direction to the interface: profile of the gradient (top) is zoomed near the460 interface (bottom).
ical solution is given by
CE
Tl = Tl0 + (Tm − Tl,0 )
erfc ( 2√xα t )
Ts = Ts0 + (Tm − Ts,0 )
k, W/(mK) 0.5263 2.22
ρ, kg/m3 996.32 916.2
cp , J/(kgK) 4318 2050
If the temperature at the left boundary suddenly changes to Tm , solidification starts theoretically from the left boundary and the interface moves from left to right. Since the phase-field is initialized using Eq. (5), the initial position of the interface xinit cannot be at the left boundary (in that case the velocity would be infinite). Instead of that, the interface φ = 0.5 is placed at a small distance equal to 6.5w from the left boundary, in order to obtain a suitable initial φ–profile. The theoretical initial temperature distribution according to equations (37)–(38) is prescribed in the liquid and the solid, respectively. Thus, the initial position of the interface has a time shift equal to tinit = x2init /(4λ2 αs ), which is accounted for in the output values from the computation. The computed interface position and interface velocity are shown in Fig. 8, and the computed temperature distribution normal to the interface is shown in Fig. 9, all compared to the corresponding analytical solutions. The position of the interface is given by the maximum position vector of the cross-section of the interface with the mesh edges, and the velocity is taken directly from the cells within the interfacial band, since they are all equal in a time step for the one-dimensional case. As expected, using the interpolation of the temperature in the probe points allows a coarser mesh resolution. The thermal diffusivity of water is α = 1.22 · 10−7 m2 /s, and the thermal diffusion layer increases from a lower value at the beginning to δt ∼ 10−3 m at later times, as the interface velocity decreases. The coarsest mesh has about 4 times less cells in the thermal diffusion layer, which is sufficient for a high accuracy of the calculated interface temperature gradients. The computed interface position approaches the analytical solution with increasing mesh resolution, as the upper diagram in Fig. 8 shows. The relative error in the computed interface position is evaluated using the expression
AN US
∇T, K/m
-5e+05 0
l
erfc (νλ)
erf ( 2√xαs t ) erf (λ)
465
,
(37)
,
(38)
AC
where ν = αs /αl , and Tl0 and Ts0 are the temperatures in the liquid phase at x → ∞ and in the solid phase at x = 0, respectively [26]. The characteristic root λ is obtained by numerically inverting the following expression Stl Sts + eλ2 erf (λ) νe(ν 2 λ2 ) erfc (νλ)
√ = λ π,
470
475
(39)
where Sts = cp,s (Tm − Ts,0 )/L and Stl = cp,l (Tm − Tl,0 )/L are the characteristic Stefan numbers for the solid and the liquid, respectively. The position and the velocity of the interface are given by the expressions √ x = 2λ αs t, (40)
ε=
|xi − xth | · 100%, xth
(42)
where xi is the computed position of the interface obtained with different cell sizes and xth is the analytical solution. 10
ACCEPTED MANUSCRIPT
4,0 Table 2: The computed interface position and the relative error at t = 30 s.
i 1 2 3
2,0
Analytical -3 ∆x = 200.10 mm -3 ∆x = 100.10 mm -3 . ∆x = 50 10 mm
1,0
0
10
0,2
30
ε=
Analytical -3 ∆x = 50.10 mm δpp = 0.1∆x 485
0,1
0,1
Table 3: The computed interface position and the relative error in the temperature profile at t = 5 s.
5,0
10,0
15,0 t, s
20,0
25,0
i 1 2
30,0
Figure 8: The computed interface position (top) and magnitude of the interface velocity (bottom) vs. time for different cell sizes.
Analytical δpp = ∆x
275
PT
δpp = 0.1∆x
AC
255
250 0,0
495
simulation time = 5 s -3 ∆x = 50.10 mm
CE
T, K
270
260
490
ED
280
265
1,0
2,0
x, mm
δpp ∆x 0.1∆x
xi , mm 1.4488 1.4589
ε, % 1.64 0.96
The performance of the model in two dimensions is demonstrated using the setup with a geometry having dimensions of 1 mm × 1 mm, discretized with an equal number of cells in the horizontal and vertical directions, and having adiabatic boundaries on all sides. The phase-field φ is initialized according to Eq. (5) as a circle with the centre located at the centre of the mesh, the initial supercooling is ∆T = 6 K everywhere, the properties of water and ice are listed in Table 4, and the latent heat of melting is L = 333000 J/kg. Table 4: Thermophysical properties for water at ∆T = 6 K and ice.
3,0
4,0
water ice
5,0
Figure 9: The computed temperature distribution normal to the500 interface for different distances of the probe points.
480
(43)
where xδpp ,i is the computed interface position in the temperature profile obtained with different distances from the probe points and xth,T = 1.473 mm is the analytical solution at time t = 5 s. The values of the relative error are given in Table 3, which shows that the result is slightly better if the probe points are placed closer to the interface.
M
0,0 0,0
ε, % 1.08 0.94 0.68
|xδpp ,i − xth,T | · 100%, xth,T
AN US
vn, cm/s
0,2
20
t, s
xi , mm 3.3908 3.3958 3.4049
The results in Fig. 9 show that choosing a larger distance of the probe points from the interface, comparable to the cell size, yields slightly less accurate result. The relative error for the computed interface position in the temperature profile at time t = 5 s is evaluated using the expression
δpp = 0.1∆x 0,0
∆x, mm 50·10−3 100·10−3 200·10−3
CR IP T
x, mm
3,0
The values of the computed interface position and the corresponding relative error obtained with the analytical so-505 lution xth = 3.428 mm at t = 30 s are given in Table 2, which shows that the relative error decreases with increasing mesh resolution. 11
k, W/(mK) 0.5487 2.22
ρ, kg/m3 999.08 916.2
cp , J/(kgK) 4246 2050
The event of freezing along a solid substrate is very fast and finishes within only several milliseconds at a distance of a few millimeters [3]. In the present study, the simulation time is set to 1 s and the computed final temperature field is shown in Fig. 10 for three mesh resolutions, with the indicated φ = 0.5 iso-contour. With the coarsest mesh, the interface starts developing as a circular shape that tends to evolve to a rectangular shape at later times, as can be seen in the top snapshot in Fig. 10. This is attributed to
ACCEPTED MANUSCRIPT
ED
M
AN US
CR IP T
the fact that the thickness of the thermal diffusion layer therefore, there is a greater total temperature difference quickly increases, as the interface velocity decreases. 520 between Tm at the interface and the temperature at the corners. Consequently, this leads to higher temperature gradients at the interface and higher interface velocity, according to Eq. (3), leading to the faster progression of the interface along the diagonals towards the corners. By in525 creasing the mesh resolution, the interface shape is closer to being circular, as can be seen in the middle and the bottom snapshots in Fig. 10. This can be seen more clearly in Fig. 11, in which separate final iso-surface contours are plotted with different 530 colors for different mesh resolutions used in the simulations.
535
AC
CE
PT
Figure 11: The final iso-contour shapes for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
Figure 10: The computed final temperature field for different mesh540 resolutions: 50×50 cells (top), 100×100 cells (middle) and 200×200 cells (bottom). 510
515
At later times, the thermal diffusion layer reaches the domain boundaries, and the iso-contour shape starts de-545 forming. The deformation occurs in the diagonal directions, since the temperature distribution is steeper along the diagonal, compared to the horizontal and vertical directions. Since the distance from the interface to the corners of the mesh is greater in the diagonal directions, it550 takes more time for heat to diffuse in that directions, and 12
The computed temporal evolution of the interface and the temperature field at early times is shown in Fig. 12 for the finest mesh resolution with 200×200 cells. As expected, due to the higher conductivity and thermal diffusivity of ice, the temperature within the ice reaches the melting temperature Tm at a shorter time, whereas the temperature in the liquid water requires a longer time to increase. The reconstruction of the interface points is demonstrated in Fig. 13, which shows the computed cross-section points of the iso-surface with the mesh edges, according to Eq. (12), and the calculated interface points, according to Eq. (10). It can be seen that the interface points are located between the cross-section points of the interface with the mesh edges. The zoomed detail in the bottom snapshot in Fig. 13 shows this more clearly. The diffused interface band contains several cells in the direction normal to the interface, and the corresponding interface point needs to be found for all cells within the band in both phases (in order to apply the condition in Eq. (3)). This is why sev-
ACCEPTED MANUSCRIPT
CR IP T AN US M ED CE
PT
Figure 13: The computed cross-section points of the iso-surface with the mesh edges (Eq. (12), black symbols) and the interface points (Eq. (10), red symbols) for the mesh with 200×200 cells. The bottom snapshot is a zoomed detail from the top snapshot.
AC
555
eral interface points may be located in each segment of the interface between two consecutive cross-section points of the iso-surface with the mesh edges. They are mapped to the corresponding cells within the band, depending on the orientation of the interface according to Eq. (10).
560
Figure 12: The computed temporal evolution of the interface and the temperature field for the mesh with 200×200 cells.
13
In order to quantify the model performance, the interface velocity is evaluated twofold, as the velocity of a single interface point, and the mean velocity within the diffused interface band. First, the velocity of the interface point moving along the horizontal direction is calculated from the expression xt − xt−∆t , (44) vn,t = ∆t where xt and xt−∆t are the instantaneous distances of the interface point tracked along the horizontal direction at the time t and at the previous time t − ∆t, respectively, from its initial position. The position of that point is readily available from the model, as the maximum position vector of the cross-section of the iso-surface with the mesh edges. Fig. 14 shows the magnitude of the computed velocity during the simulation time, and it is seen that the veloc-
ACCEPTED MANUSCRIPT
575
580
0,04 0,02
400
t, ms
600
800
1000
600
PT
Figure 14: The magnitude of the computed velocity of the interface point along the horizontal direction for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
800
1000
where P/π is a measure for the mean diameter D of the iso-contour. In the case of an ideal circular iso-contour, D would be equal to its real diameter. Additionally, the same quantity is calculated in another way, by summing all the volumes of the cells (m in total) identified as ’solid’, for which φ ≥ 0.5, and dividing the sum by the mesh thickness, in which case the amount of ice is determined by Pmthe two-dimensional area covered by the solid cells, A = i=1 ∆Ai . The computed temporal evolution of the isocontour perimeter is shown in Fig. 16, and the evolution of the computed ratio AD /A is plotted in Fig. 17.
CE
1,0 P, mm
590
600
1,2
Next, the interface velocity is represented in terms of the mean velocity in the cells, which are located within the diffused interface band in a single time step. The arithmetic mean of the velocity is calculated by summing the velocity magnitudes in all cells lying within the [φmin , φmax ] interfacial band and dividing the sum by the number of these interfacial cells. The computed mean interface velocity magnitude within the band is plotted in Fig. 15, in which it is seen that the mean velocity does not show oscillations. In order to further quantify the predictive capabilities of the model, in each time step the amount of the new formed ice (solidified water) is calculated in two different ways. First, in each time step the iso-contour perimeter is computed n X P = πD ≈ ∆li , (45)
AC
585
t, ms
Figure 15: The computed mean interface velocity magnitude within the diffused interface band for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
595
200
400
crossing the mesh-edges. It is equal to the area of the isosurface divided by the thickness of the mesh. The amount of the solidified water is then approximated as the area in two dimensions 2 π P π 2 , (46) AD = D = 4 4 π
ED
0,00 0
200
AN US
0,06
0,04
0,00 0
M
vn , cm/s
0,08
0,06
0,02
0,10 50x50 100x100 200x200
50x50 100x100 200x200
0,08
CR IP T
570
0,10
ity oscillates slightly around a mean value. These oscillations are attributed to the time instants when the interface passes the cell-centres, i.e. the cell that was previously identified as liquid now becomes solid. For the transient term in the energy equation, the temperature Tm is set as the old-time value in that cell when this occurs. The exact time instant when the iso-surface is located exactly at the cell-centre cannot be determined within a time step; otherwise some kind of temporal extrapolation for the temperature would have to be applied. In addition, the stencils for the interpolation of the temperature in the probe points change, all leading to somewhat non-gradual temporal change of the temperature gradients at the interface and the oscillating velocity. In spite of that, it should be noted that the oscillations are not noticeable in the onedimensional simulation, and the amplitudes of the oscillations are relatively small, also in the two-dimensional case.
vn,mean , cm/s
565
0,8 50x50 100x100 200x200 0,6 0
200
400
t, ms
600
800
1000
Figure 16: The computed temporal evolution of the iso-contour perimeter for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
i=1
by summing all the interface line segments ∆li (n in total)
14
ACCEPTED MANUSCRIPT
1,3
AD / A, -
1,1 1,0
4Ld
symmetry
1,2
adiabatic
adiabatic 50x50 100x100 200x200
0,8 0,7 0
adiabatic
200
400
t, ms
600
800
1000
2Ld
Figure 17: The computed ratio AD /A for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
Figure 18: Initial configuration for the simulation of the ice dendrite growth.
635
615
630
phase-field φ is initialized according to Eq. (5), where ndist is calculated as the minimum distance from the cell-centres within the phase-field band to the profile of the parabolic dendrite with φ = 0.5. At a supercooling of 10 K, the dendrite tip velocity predicted by the MST is 4.9 cm/s with the corresponding theoretical tip radius of R = 2.45 · 10−7 m. The initial parabolic φ–profile is y = ax2 +bx+c, where b = 0, c = Ld , and the coefficient a = −1/(2R) is obtained from the condition that the radius of curvature is R in the vertex of the parabola located at (0, Ld ). The simulation is run until the temperature at the adiabatic boundaries (right and top) starts changing, indicating that the thermal diffusion layer has reached the boundary. Hence, the assumption of a semi-infinite space in which the dendrite grows is not violated and the steady-state dendrite tip velocity is reached. The computed interface position and velocity of the interface point at the symmetry plane for different mesh resolutions is plotted in Fig. 19. The velocity at a time t is calculated using Eq. (44) by replacing x with the maximum coordinate y of the iso-surface point at the symmetry plane. At the beginning of the simulation, the temperature gradients at the interface are highest; hence, the dendrite tip moves very fast. Shortly after that, the temperature field within both phases develops and the dendrite tip velocity decreases. The computed dendrite tip velocity quickly converges to the theoretical value, as the lower diagram in Fig. 19 shows. The computed interface position converges with increasing mesh resolution, as can be seen in the upper diagram in Fig. 19. The relative error is quantified in terms of the value of the computed interface position obtained with the finest mesh, from the expression
AN US
AC
625
CE
PT
620
M
610
As can be seen in Fig. 16, the plot for the temporal evolution of the perimeter of the iso-contour has the characteristic shape, similar to that obtained in the one-dimensional case in Fig. 8. More important for justification is the ratio AD /A, which converges to unity with mesh refinement, as640 Fig. 17 shows. Finally, the described model is used to compute the growth of an ice dendrite into supercooled water. The transition from planar stable freezing to the unstable dendritic growth is characterized in the marginal stability the-645 ory (MST) [35, 36, 37]. If the solid-liquid interfacial energy is neglected, the product of the velocity of the dendrite tip and its radius at a specific supercooling is determined by the solution of Ivantsov [23] for the growth of a parabolicshaped dendrite in semi-infinite space. When the interfa-650 cial energy is accounted for, the individual radius of the dendrite tip, which is dependent on the given supercooling, is selected by using MST. Basically, the dendrite tip velocity is a single-valued function of the tip radius at a given supercooling, and every selected value of the dendrite tip655 radius determines a corresponding exact steady-state tip velocity. The setup for the simulation is sketched in Fig. 18. The geometry is 3.5 µm ×7 µm, discretized with cells having equal dimensions in horizontal and vertical directions, with660 all adiabatic boundaries, except the left boundary, which is a symmetry plane. The initial supercooling is ∆T = 10 K everywhere, the properties of water and ice are listed in Table 5, and the latent heat of melting is L = 333000 J/kg.
ED
605
CR IP T
Ld
0,9
Table 5: Thermophysical properties for water at ∆T = 10 K and ice.
water ice
k, W/(mK) 0.5392 2.22
ρ, kg/m3 998.13 916.2
cp , J/(kgK) 4272 2050
ε=
|yi − y4 | · 100%, y4
(47)
where yi is the position of the interface point at the symmetry plane obtained in simulations with coarser meshes
The length of the initial dendrite is Ld = 1.75 µm. The 15
ACCEPTED MANUSCRIPT
0,15
expression 50x100 100x200 200x400 400x800
675
0,05
0,00 0e+00
2e-04
4e-04
100
t, ms
6e-04
1e-03
680
60
AN US
theory (MST) 50x100 100x200 200x400 400x800
80
vn , cm/s
8e-04
40 20 0 0e+00
4e-04
t, ms
6e-04
8e-04
1e-03
M
2e-04
CR IP T
y, µm
0,10
y2 − y3 y3 − y4 p= = 1.307, (48) ln 2 where the factor 2 in the denominator is given by the mesh refinement ratio ∆xi /∆xi−1 . The spatial accuracy is not second order, which is expected, since the discretization of the spatial terms in the transport equations is not second order. Finally, the computed interface shape with the developed temperature field and the velocity field at t = 0.001 ms is plotted in Fig. 20 for the mesh with 200×400 cells. Shortly after the beginning of the computation, the temperature has reached the value Tm within the dendrite, and heat is diffused into the liquid, producing parabolic isotherms around the φ = 0.5 iso-contour. The highest computed velocity is in the region at the tip of the dendrite. ln
PT
CE
670
(50×100, 100×200 and 200×400), and y4 is the one obtained in the simulation with the finest mesh (400×800) at the same simulation time. The values of the computed position of the interface point at the symmetry plane at t = 0.001 s, and the corresponding relative error calculated with Eq. (47), are given in Table 6, which shows that the relative error decreases with increasing mesh resolution.
Table 6: The computed position of the interface point at the symmetry plane and the relative error at t = 0.001 s.
AC
665
ED
Figure 19: The computed interface position (top) and magnitude of the interface velocity (bottom) of the interface point at the symmetry plane vs. time for different mesh resolutions.
i 1 2 3 4
mesh 50×100, ∆x 100×200, ∆x/2 200×400, ∆x/4 400×800, ∆x/8
yi , µm 0.083 0.0911 0.0958 0.0977
ε, % 15.05 6.76 1.94 0
Figure 20: The computed temperature field and the velocity field within the diffused interface band for the mesh with 200×400 cells. 685
The values of the relative error in Table 6 indicate a690 convergence rate higher than second order, however this should not be confused with the accuracy of the results. The order of accuracy is estimated using the results for y for the three finest meshes (2, 3 and 4 in Table 6) from the 16
It can be concluded from the results and discussion that the model is capable of predicting the process of solidification of supercooled water with the accompanying thermal and phase-change effects on the moving ice-water interface. This is demonstrated in particular by the correct prediction of the existing analytical results for the Stefan moving boundary problem and the theoretical results for the ice-dendrite growth into supercooled water as the commonly used representative analytical solutions. The
ACCEPTED MANUSCRIPT
(50)
which is internally conserved without the need to explicitly transfer the fluxes at the shared boundary during the computation. This functionality was successfully used previously for the simulation of non-isothermal drop impact [7] and is also utilized here. The effective solidliquid conductivity at the wall cell-faces within the band φmin < φ < φmax , is calculated from the expression kls = φks + (1 − φ)kl , which, combined with the conductivity of the wall material, is then used for the harmonic interpolation of the conductivity at the wall cellfaces. The interpolation stencils in Eq. (31) are extended to include the values for the temperature at the cell-faces at the coupled boundary as well. At the coupled boundary the phase-field is extrapolated with zeroth order from the cells adjacent to the substrate surface, thus imposing a zero-gradient boundary condition. The thickness of the thin ice film developing along the cold substrate was not reported in [3], whereas it was estimated at 20 µm in [6]. In the present study, the propagation of the thin ice film developing along the cold underlying substrate is computed for the initial supercooling of ∆T = 6.25 K. The two-dimensional computational domain shown in Fig. 22 has two regions: the upper for the water and ice and the lower for the substrate. Both have dimensions in the vertical plane of 50 µm × 50 µm, discretized with equal number of cells in the horizontal and vertical directions.
AN US
Modeling the evolution of the thin ice film due to solidification of supercooled water in a thin region over a cold substrate requires solving energy equations in all three phases, i.e. within water and ice, and in the substrate as720 well, in order to account for the conjugate heat transfer at the wall surface. When the underlying substrate, shown in the sketch in Fig. 21, has a much greater thermal diffusivity compared to supercooled water and ice, and its thermal influence on the dynamics of the spreading of the725 thin ice film is assumed to be purely as a heat absorber.
730
vn
substrate (wall)
M
735
Figure 21: Sketch of the thin ice layer due to solidification of supercooled water over a cold substrate.
ED
The additional transport equation to be solved is the energy diffusion equation in the substrate (wall), which reads ∂T ρw cpw = ∇ · (kw ∇T ), (49) ∂t
PT
T =T¥
where the subscript w denotes the wall. The energy equation (49) is coupled with the basic model defined by equations (1)–(4). Considering the discretization of the energy equations, the software package foam-extend provides a functionality to solve a single energy equation for the entire domain, consisting of the mesh for the fluid, where the phase change takes place, and the mesh in the substrate. The coupling between the two meshes makes use of the software’s functionality of connecting the two meshes for a single solution of the temperature field. The discretization of the energy equation follows the standard finite-volume procedure. Special care is given to the cell-face interpolation of the conductivity in the diffusion term at the cell-faces shared by both meshes. According to [38], the cell-face value of the conductivity740 is evaluated using harmonic interpolation. Thus, the energy equation is solved for both meshes simultaneously, by extending the solution matrix for the energy equation to include the mesh for the substrate as well. The coupling between the ice-water and the substrate domains is745 17
T =T¥
water (liquid) ice (solid)
coupled
T =T¥
interface, Tm
CE
710
715
AC
705
−kl/s (∇T |⊥ )l/s = −kw (∇T |⊥ )w ,
CR IP T
4. Extension to conjugate heat transfer
achieved by fulfilling the requirement that the temperature be continuous at the wall surface, and the heat fluxes be conserved in the normal direction to the wall surface
symmetry
700
same problems were also recently used for the assessment of the volume-of-fluid and the level-set methods, and by comparison with these results (c.f. for example the results in [10]), it can be seen that the present model performs well in terms of accuracy of the predicted position, shape and velocity of the interface.
symmetry
695
T =T¥ Figure 22: Initial configuration for the simulation of the thin ice film propagation over the cold underlying substrate.
At the initial time the ice nucleus is initialized using Eq. (5), having the shape of a quarter of a circle, with the centre placed in the lower left corner. Melting temperature Tm is set within the ice at the beginning, and the uniform temperature T∞ = 266.9 K corresponding to the initial supercooling is prescribed outside. At all outer
ACCEPTED MANUSCRIPT
k, W/(mK) ρ, kg/m3 cp , J/(kgK)
water, ∆T = 6.25 K 0.5487 999.08 4246
ice 2.22 916.2 2050
copper 398 8954 384
760
765
The computed interface position and temperature distribution around the ice film at different times are shown in Fig. 23 for the mesh with 100 cells in both directions.
CE
PT
ED
M
AN US
770
During the evolution of the ice film, heat is being transported from the ice-water region into the substrate. Hence, the temperature within the ice film does not stay at the initial value of Tm , but it decreases from Tm at the interface to a lower value at the wall surface. The temperature change normal to the interface in Fig. 23 is much less sharp far from the tip of the ice film at the wall. Around the tip, the thermal diffusion layer in the normal direction to the interface remains much thinner, yielding much higher temperature gradients in the normal direction to the interface close to the wall. This is expected, since the underlying substrate has much higher thermal diffusivity and absorbs effectively the heat released during the solidification of water. Therefore, the horizontal propagation velocity of the ice film tip at the wall is much higher compared to the velocity of the interface in the vertical direction, leading to the obtained elongated shape of the developed ice film. The computed velocity distribution around the tip of the ice film at t = 0.5 ms is shown in Fig. 24, for the same mesh. In accordance with the above explanation, the velocities in Fig. 24 within the band [φmin , φmax ] are highest in the computational cells closest to the wall surface.
CR IP T
Table 7: Thermophysical properties for water, ice and substrate.
AC
750
boundaries the temperature is kept constant and equal to T∞ , except at the left boundary, which is the plane of symmetry. Thermophysical properties for the supercooled755 water, ice and substrate made of copper used in the calculations are listed in Table 7, and the latent heat of melting is L = 333000 J/kg.
Figure 24: The computed velocity distribution around the tip of the ice film at t = 0.5 ms, for the mesh with 100×100 cells, the supercooling ∆T = 6.25 K and copper substrate. 775
780
Figure 23: The computed interface position and temperature distribution around the ice film at different times for the mesh with 100×100 cells, the supercooling ∆T = 6.25 K and copper substrate. The time sequence is from top to bottom, left to right.
785
18
The computed temporal evolution of the interface is shown in Fig. 25, in which the evolution of the iso-contour point moving along the wall, and of the point moving along the symmetry plane is plotted for different mesh resolutions. The magnitude of the computed corresponding propagation velocity for the same two points is plotted in Fig. 26. The computed distribution of the temperature along the wall and the symmetry plane is shown in Fig. 27 at different time instants for different mesh resolutions. From Figs. 25–27, several conclusions can be drawn. The velocity at a time t is calculated from Eq. (44) and the computed interface propagation velocity at the wall again exhibits oscillations, but the amplitude is much higher in this
ACCEPTED MANUSCRIPT
20
0,05
0 0
0,02
3 0,3 0,4 t, ms
0,5
0,6
0,7
vn , cm/s
0,2
0,2
AN US
0,1
0,007
0,3 0,4 t, ms
0,5
0,6
0,7
100×100 200×200 300×300
2
1
0,006 y, mm
0,1
4
100×100 200×200 300×300
0,008
0,005
M
0,004
100×100 200×200 300×300 0,1
0,2
ED
0,003 0,002 0
CR IP T
x, mm
5
0,03
0,01
0,3 0,4 t, ms
0,5
0,6
0,7
CE
PT
Figure 25: The computed temporal evolution of the iso-surface points at the wall (top) and at the symmetry plane (bottom) for different800 mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate.
AC
However, the mean velocity, which was determined in experiments to be constant at a specific supercooling and a specific substrate [3], also stays constant in the simula-805 tion (except at the beginning, before the temperature field develops). Therefore, the magnitude of the mean velocity at a time t, plotted with solid lines in Fig. 26, is calculated in the simulation from the expression vn,t
795
15 10
0,04
0 0
exp. at T = 266.9 K 100×100 200×200 300×300
25
vn , cm/s
790
30
case, as seen in the dashed-line plot in the upper diagram in Fig. 26. As previously explained, the oscillations are attributed to the time instants when the interface passes the cell-centres. Similar behavior is detected for the velocity at the symmetry plane, but with much lower amplitudes.
810
t
1 X = Pt (xp,t − xp,t−∆t ). 0 ∆t 0
(51)
Another observation is the difference in the computed values for the interface velocity at the wall and at the815 symmetry plane. While the computed values at the symmetry plane converge with increasing mesh resolution, the values at the wall are not mesh independent. With increasing mesh resolution, the interface propagation veloc19
0 0
0,1
0,2
0,3 0,4 t, ms
0,5
0,6
0,7
Figure 26: The magnitude of the computed interface propagation velocity at the wall (top) and at the symmetry plane (bottom) for different mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate. Solid lines are plots of the mean velocity, and the dashed line is the instantaneous velocity.
ity at the wall also increases and does not converge towards a finite value (the experimentally determined value about 15 cm/s for this case [3]). This is further illustrated in Fig. 28, which shows that the interface velocity at the wall linearly increases with a linear increase in the mesh resolution. Similar behavior is detected also in simulations with a higher supercooling of ∆T = 12.35 K, the results of which are plotted in Fig. 29. On the other hand, the mesh independence of the solution at the symmetry plane is confirmed in all results for the interface position, its propagation velocity and the temperature distribution along the symmetry plane. As stated above, the mean velocities, shown in Fig. 26 and Fig. 29, are constant in time, and the small increase at later times is solely because the interface approaches the boundary of the computational domain, where the temperature is kept fixed. Inspection of the results in Fig. 27 reveals that the temperature distribution at the wall has a similar shape for all mesh resolutions at different times. With increasing the mesh resolution, the computed temperatures in cells clos-
ACCEPTED MANUSCRIPT
30
276
25
T, K
vn , cm/s
272 270
0,02
0,01
274
x, mm
0,03
0,04
10 0
0,05
red: 100×100 blue: 200×200 green: 300×300
272
825
268
0,02
0,03
0,04
0,05 830
M
y, mm
PT
exp. at T = 266.9 K comp.
840
CE
20 15
AC
vn , cm/s
25
10
ED
Figure 27: The computed temperature distribution along the wall (top) and the symmetry plane (bottom) at different time instants for different mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate. 835
30
100
200 300 number of cells at the wall
845
400 850
Figure 28: The increase in magnitude of the computed interface propagation velocity at the wall with increasing mesh resolution, for the supercooling ∆T = 6.25 K and copper substrate. 855
820
0,4 t, ms
0,6
0,8
expressions (28)–(29). At first sight, the problem could be overcome by applying a constant distance δpp for each mesh. However, in that case the thermal diffusion layer could be skipped for a fine enough mesh, and a further mesh refinement would then even lead to a decrease in the calculated temperature gradient. The mesh dependence of the results at the wall is not attributed to the fixed-value temperature condition at the boundaries. The computations with the same initial conditions were performed using adiabatic boundaries instead of T = T∞ in Fig. 22, except at the bottom of the substrate. The results (not shown here) were similar to the ones obtained with the fixed temperature at the boundaries. This is explained as the consequence that the thermal diffusion layer around the interface is very thin and does not reach the right boundary for the most of the simulation time. The reason for the mesh dependence of the solution at the wall is also not the small dimensions of the computational domain. The thermal diffusion layer in the water for the specific conditions is estimated as δt ∼ αl /vn = 8.3 · 10−7 m, and even with the small computational domain, having dimensions of 50 µm × 50 µm, the solution should not be influenced by the conditions at the boundaries, at least at the beginning of the simulation. This is confirmed by the results at the symmetry plane, that do not show such mesh dependence. It is further demonstrated in Fig. 30, which shows the results for the computed interface velocity, obtained in simulations with a computational domain that has 20 times longer horizontal and vertical dimensions. The interface propagation velocity is only slightly smaller with the much larger computational domain, indicating that the outer boundaries, toward which the interface moves, have negligible influence on the solution for the evolution of the ice-film tip at the wall. The freezing of water directly at the substrate surface is determined by the local thermal conditions around the tip, thus confirming the assumption for the substrate acting purely as an
AN US
T, K
270
0,01
0,2
Figure 29: The magnitude of the computed interface propagation velocity at the wall for different mesh resolutions, the supercooling ∆T = 12.35 K and copper substrate.
solid line: t = 0.1 ms dashed line: t = 0.2 ms
266 0
20
15
268 266 0
100×100 200×200 300×300
solid line: t = 0.1 ms dashed line: t = 0.2 ms
CR IP T
274
red: 100×100 blue: 200×200 green: 300×300
est to the interface do not change much, but the calculated temperature gradients increase with mesh refinement due do smaller distances to the probe points, according to the 20
ACCEPTED MANUSCRIPT
895
25 exp. at T = 266.9 K -5 domain (5×5).10 m
vn , cm/s
20
-3
domain (1×1).10 m 900
15 10
905
5 0 0
0,1
0,2
t, ms
0,3
0,4
0,5 910
AN US
Figure 30: The comparison of the magnitude of the computed interface propagation velocity at the wall for the supercooling ∆T = 6.25 K and copper substrate, with that computed using a computational domain with 20 times longer dimensions. 860
880
885
890
M
ED
PT
875
CE
870
915
The conclusion is that the described model exhibits a good capability to compute the two-phase solidification far from the wall. However, the solidification at the wall, in particular the contact line (triple point water-ice-substrate in a two-dimensional case), imposes a singularity, which is920 indicated by the continuously increasing interface velocity when the triple point is approached by refining the mesh. Similar behavior was detected also in the aforementioned previous related studies, where the interface was even not perpendicular to the wall surface. In a recent925 study [39], where numerical simulations of a drop solidifying on a cold substrate were conducted, a three-phase point which might appear similar to the present study was encountered. However, it is related to a different underlying physics, since the studied cases include slow solidification930 with several orders of magnitude larger time scales, and the triple point, which connects the same liquid and solid material, moves away from the substrate and not along it. Thus, the triple point is of completely different nature compared to the problem in the present study. Obviously,935 in order to model this phenomenon for three materials of different kind, as a prerequisite a better understanding is required of the physics in the local area around the icing front at the substrate.
AC
865
balance, accounted for in the discretization by using reconstructed interface points. The model has been tested by computing two-phase solidification problems with existing analytical solutions. An extension of the computational model to incorporate conjugate heat transfer has been presented that can be used as a framework for combining additional models to account for the interaction of the freezing of the supercooled water with the cold, solid substrate. It has been demonstrated that the model has a good capability to compute the solidification in the absence of walls, or far from the wall. At the wall, the problem formulation imposes a singularity, indicated by diverging velocities at the contact line, which is formed during solidification of liquid over a solid substrate. The heat flux singularity at the contact line presents a fundamental difficulty for modeling the freezing within the finite-volume continuum mechanics. In accordance to the previous studies, this intrinsic singularity is also reproduced with the present model. New understanding of the underlying physics must be acquired, in order to formulate plausible models which would remove the singularity and provide means for using them within the conventional continuum framework. Experimental data at extremely small length and time scales around the contact line would be a very useful basis, though hardly accessible in measurements, in particular in the case of a contact line moving at relatively high speeds. If such new models would require including the capillary undercooling in regions of very large curvatures of the interface, the present model can be easily extended by calculating the curvature-dependent temperature and replacing the melting temperature in the points at the interface. The model can be used for computing the solidification of supercooled water, in which the phase-change process can be isolated and analyzed separately without the bulk liquid motion. As opposed to the standard interfacetracking algorithms, such as volume-of-fluid, where the Navier-Stokes flow solver is coupled with an algorithm for the proper advection of the tracked interface by the underlying velocity field, the interface in the present model propagates with the normal velocity obtained from the interfacial energy balance calculated only within the diffused interface band. Thus, the interface propagation is determined here using an algorithm based on a different underlying physics compared to the standard volume-of-fluid interface tracking. In terms of the computational costs, the model will by definition perform faster than the level-set method, since the re-initialization of the phase-field variable and the computation of the extended velocity field, required for the re-initialization in the entire computational domain, is not necessary in the present model.
CR IP T
effective heat absorber.
940
5. Conclusions
A computational model for the freezing of supercooled water within the standard cell-centered finite-volume CFD framework has been presented. The model uses the phasefield approach for tracking the interface between the solid945 Acknowledgement and the liquid phase. Energy equations in both phases are decoupled at the interface by imposing the melting temWe gratefully acknowledge financial support from the perature at the interface as a boundary condition. The Deutsche Forschungsgemeinschaft within the collaborative velocity of the interface is determined by the local energy research project SFB-TRR 75 (TP-C3). 21
ACCEPTED MANUSCRIPT
980
985
990
995
1000
1005
1010
1015
[23]
[24] [25]
[26]
[27] [28]
[29]
[30] [31]
equiaxed dendritic solidification of a binary alloy, Computational Materials Science 112 (2016) 304–317. G. Ivantsov, Temperature field around spherical, cylindrical, and needle-shaped crystals which grow in supercooled melt, Doklady Akademii Nauk SSSR 58 (1947) 567–569. B. Cantor, A. Vogel, Dendritic solidification and fluid flow, Journal of Crystal Growth 41 (1977) 109–123. F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, An Advanced InR and Matlab , R Springer Intertroduction with OpenFOAM national Publishing Switzerland, 2016. V. Alexiades, A. Solomon, Mathematical modeling of melting and freezing processes, Hemisphere publishing corporation, Washington, 1993. F. Martinez, Beyond the classical stefan problem, Ph.D. thesis, Universitat Polit` ecnica de Catalunya (2014). Y. Sun, C. Beckermann, A two-phase diffuse-interface model for hele-shaw flows with large property contrasts, Physica D 237 (2008) 3089–3098. Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, Journal of Computational Physics 220 (2007) 626–653. R The OpenFOAM extend project, http://www.extendproject.de/. R OpenFOAM , The open source CFD toolbox, http://www.openfoam.com/. ˇ Tukovi´ P. Cardiff, Z. c, H. Jasak, A. Ivankovi´ c, A block-coupled finite volume methodology for linear elasticity and unstructured meshes, Computers and Structures 175 (2016) 100–122. A. Criscione, Influence of ice formation on drop dynamics, Ph.D. thesis, Technische Universit¨ at Darmstadt (2014). A. Shibkov, Y. Golovin, M. Zheltov, A. Korolev, A. Leonov, Morphology diagram of nonequilibrium patterns of ice crystals growing in supercooled water, Physica A 319 (2003) 65–79. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– I. elements of a stability analysis, Acta Metallurgica 26 (1978) 1681–1687. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– II. instabilities in the limit of vanishing surface tension, Acta Metallurgica 26 (1978) 1689–1695. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– III. effects of surface tension, Acta Metallurgica 26 (1978) 1697– 1708. S. Patankar, Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, Washington, DC, 1980. T. Vu, G. Tryggvason, S. Homma, J. Wells, Numerical investigations of drop solidification on a cold plate in the presence of volume change, International Journal of Multiphase Flow 76 (2015) 73–85.
AN US
975
M
970
ED
965
PT
960
CE
955
[1] E. Moore, V. Molinero, Structural transformation in supercooled water controls the crystallization rate of ice, Nature 479 (2011) 506–508. [2] F. Tavakoli, S. Davis, H. Kavehpour, Freezing of supercooled1025 water drops on cold solid substrates: initiation and mechanism, Journal of Coatings Technology and Research 12 (2015) 869– 875. [3] M. Schremb, C. Tropea, Solidification of supercooled water in the vicinity of a solid wall, Physical Review E 94 (2016) 052804.1030 [4] M. Schremb, J. Campbell, H. Christenson, C. Tropea, Ice layer spreading along a solid substrate during solidification of supercooled water: Experiments and modeling, Langmuir 33 (2017) 4870–4877. [5] S. Jung, M. Tiwari, N. Vuong Doan, D. Poulikakos, Mechanism1035 of supercooled droplet freezing on surfaces, Nature communications 3 (2012) 615. [6] W. Kong, H. Liu, A theory on the icing evolution of supercooled water neat solid substrate, International Journal of Heat and Mass transfer 91 (2015) 1217–1236. 1040 [7] M. Schremb, S. Borchert, E. Berberovi´ c, S. Jakirli´ c, I. Roisman, C. Tropea, Computational modelling of flow and conjugate heat transfer of a drop impacting onto a cold wall, International Journal of Heat and Mass Transfer 109 (2017) 971–980. [8] J. Blake, D. Thompson, T. Strobl, Simulating the freezing1045 of supercooled water droplets impacting a cooled substrate, The American Institute of Aeronautics and Astronautics AIAA Journal 53 (2015) 1725–1739. [9] P. Rauschenberger, B. Weigand, A volume-of-fluid method with interface reconstruction for ice growth in supercooled water,1050 Journal of Computational Physics 282 (2015) 98–112. [10] P. Rauschenberger, A. Criscione, K. Eisenschmidt, D. Kintea, ˇ Tukovi´ S. Jakirli´ c, Z. c, I. Roisman, B. Weigand, C. Tropea, Comparative assessment of volume-of-fluid and level-set methods by relevance to dendritic ice growth in supercooled water,1055 Computers & Fluids 79 (2013) 44–52. [11] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, Cambridge, UK, 2011. [12] S. Deshpande, L. Anumolu, M. Trujillo, Evaluating the perfor-1060 mance of the two-phase flow solver interfoam, Computational Science & Discovery 5 (2012) 014016. [13] T. Waclawczyk, T. Koronowicz, Comparison of CICSAM and HRIC high-resolution schemes for interface capturing, Journal of Theoretical and Applied Mechanics 46 (2008) 325–345. 1065 [14] A. Bourdillon, P. Verdin, C. Thompson, Numerical simulations of water freezing processes in cavities and cylindrical enclosures, Applied Thermal Engineering 75 (2015) 839–855. [15] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, International Journal of Heat and Mass Transfer 30 (1987) 1709–1719. [16] S. Schiaffino, A. Sonin, Motion and arrest of a molten contact line on a cold surface: An experimental study, Physics of Fluids 9 (1997) 2217–2226. [17] S. Schiaffino, A. Sonin, On the theory for the arrest of an advancing molten contact line on a cold solid of the same material, Physics of Fluids 9 (1997) 2227–2233. [18] D. Anderson, S. Davis, Local fluid and heat flow near contact lines, Journal of Fluid Mechanics 268 (1994) 231–265. [19] M. Quintella, Numerical investigation of the heat flux singularity at an advancing molten contact line, Master’s thesis, Massachusetts Institute of Technology (1999). [20] I. Steinbach, C. Beckermann, B. Kauerauf, Q. Li, J. Guo, Threediensional modeling of equiaxed dendritic growth on a mesoscopis scale, Acta materialia 47 (1999) 971–982. [21] I. Steinbach, H. Diepers, C. Beckermann, Transient growth and interaction of equiaxed dendrites, Journal of Crystal Growth 275 (2005) 624–638. [22] Y. Souhar, V. De Felice, C. Beckermann, H. Combeau, M. Zaloˇ znik, Three-dimensional mesoscopic modeling of
AC
950
1020
CR IP T
References
22
[32]
[33] [34]
[35]
[36]
[37]
[38] [39]
Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points ˇ Edin Berberovi´c, Markus Schremb, Zeljko Tukovi´c, Suad Jakirli´c, Cameron Tropea PII: DOI: Reference:
S0301-9322(17)30646-8 10.1016/j.ijmultiphaseflow.2017.11.005 IJMF 2676
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
29 August 2017 16 October 2017 2 November 2017
ˇ Please cite this article as: Edin Berberovi´c, Markus Schremb, Zeljko Tukovi´c, Suad Jakirli´c, Cameron Tropea, Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.11.005
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
ACCEPTED MANUSCRIPT
Highlights • Development of a new computational framework for modeling the solidification of supercooled liquid water. • Utilization of the combination of the phase-field front propagation approach with the thermal diffusion equations in the phases.
CR IP T
• Decoupling at the moving phase-interface by using immersed points to impose the melting/freezing temperature at the interface. • Robust and efficient implementation within the object-oriented framework in foam-extend, the extenR sion of OpenFOAM • Extension of the numerical model to include conjugate heat transfer with the neighboring, thermally coupled solid wall.
AN US
• The validation of the numerical model demonstrates that it has a good capability to compute the solidification in the absence of walls, or far from the wall.
• At the wall, the problem formulation imposes a singularity as a fundamental difficulty at the triple contact line, formed during solidification of liquid over a solid substrate, which is also confirmed.
AC
CE
PT
ED
M
• The numerical model can be used as a framework for investigations of the interaction of the solidifying supercooled water with the thermally coupled solid substrate, by inclusion of and combining with suitably formulated plausible models, which would eventually remove the singularity.
1
ACCEPTED MANUSCRIPT
Computational modeling of freezing of supercooled water using phase-field front propagation with immersed points ˇ Edin Berberovi´ca,∗, Markus Schrembb , Zeljko Tukovi´cc , Suad Jakirli´cb,∗, Cameron Tropeab a Polytechnic faculty, University of Zenica, Zenica, Bosnia and Herzegovina of Fluid Mechanics and Aerodynamics, Technische Universit¨ at Darmstadt, Darmstadt, Germany c Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Zagreb, Croatia
CR IP T
b Institute
Abstract
AN US
Computational modeling of phase change due to solidification, besides correctly capturing the heat transfer, requires accurate evaluation of mass transfer at the phase-interface. In the present study a framework is developed for modeling the solidification of supercooled water using a phase-field approach for the propagation of the interface between ice and supercooled water. Energy equations in the solid and liquid phases are decoupled at the interface by using immersed points to impose the melting temperature as a moving boundary condition. The propagation of the interface is determined by the local energy balance across the interface, which is accounted for in the reconstructed interface points. The model is validated using the known theoretical solutions for the two-phase Stefan problem. An extension of the model to incorporate conjugate heat transfer with a neighboring solid wall is presented, establishing a framework that can be used for further modeling of the interaction of solidifying supercooled water on a solid substrate. Keywords: supercooled water, freezing process, phase-field method, thermal diffusion model, three-phase contact point
M
1. Introduction
15
20
25
ED
PT
CE
10
AC
5
Pure water freezes at a temperature of 273.15 K, but under certain conditions it can be cooled down to almost 40 K lower temperatures, while retaining its liquid ther- 30 modynamic state [1]. The process of subsequent solidification of supercooled water is of essential interest for diverse engineering applications, where solid surfaces are exposed to different environmental conditions, common examples being wind turbines, or electric power line sys- 35 tems. A familiar practical application where the freezing of supercooled droplets poses a serious problem is in aviation, where droplets of supercooled water levitate in clouds. Aircraft passing through such clouds collides with the droplets, triggering their abrupt solidification, which 40 can lead to the formation of ice layers on aircraft surfaces, within the engine, or on instruments. In the freezing of supercooled water several stages are identified [2], and the main influencing parameters include the degree of supercooling, the environmental conditions, 45 and the solid surface temperature and structure [3, 4, 5, 6]. Understanding the processes in solidifying supercooled water and their interaction with a solid surface is of practical significance for defining methodologies for their prevention and deicing. However, even in the absence of phase change, 50 the non-isothermal two-phase flow in a supercooled water ∗ Corresponding
authors Email addresses: [email protected] (Edin Berberovi´ c ), [email protected] (Suad Jakirli´ c) Preprint submitted to International Journal of Multiphase Flow
55
drop impacting a surface represents a challenge, due to its high unsteadiness and the small characteristic length and time scales [7], which may span several orders of magnitude. The situation is additionally complicated if one of the phases undergoes changes in its thermodynamic state. The small scales render direct and detailed experimental access practically impossible, which is why research relies extensively on numerical tools. Therefore, there is considerable motivation for devising computational methods for solving the solidification of supercooled water, which are thermally coupled with the solid surfaces and enable examination of the governing physical mechanisms of the interaction. In the commonly applied numerical approaches, the ice front is treated as a liquid-solid phase-interface which evolves in time and space driven by the phase change [8, 9, 10]. The propagation of the interface is represented using either Lagrangian-based interface tracking or Eulerian-based volume capturing approach [11], coupled with heat transfer and a model for the mass transfer due to the phase change. The interface tracking provides a sharp representation of the interface by definition, whereas the volume capturing relies on the homogeneous mixture model and yields a smeared interface. The widely used volume-of-fluid approach has become an attractive methodology, utilizing an algebraic scheme for the propagation of the interface [12, 13], which can be combined with the enthalpy-porosity technique, or some if its variants, to model mass transfer [14, 15]. In spite of the work conducted in recent years using different modeling and numerical techniques, there is still a need to improve November 4, 2017
ACCEPTED MANUSCRIPT
85
90
95
100
105
110
CR IP T
iso-surface by virtually connecting the tips of dendrites in a propagating dendritic cloud. The main idea was to combine the numerical solution of the heat diffusion equation in the liquid around this cloud at a larger, mesoscopic scale, with a local analytical solution for the growth of a single dendrite tip, relevant at smaller, microscopic scales. The analytical solution is matched with the computationally obtained temperature field at a predefined distance from the iso-surface. In a modified version, the model was also successfully used in [22] for the solidification of an undercooled binary alloy. The analytical solution, applied for the velocity of the iso-surface, is the so-called stagnant film solution and represents a modification of the Ivantsov solution for the growth of a dendrite in the form of an isothermal paraboloid [23, 24]. Matching the analytical with the numerical solution is achieved by supplying the analytical solution with the value for the supercooling ∆T , obtained from the numerical solution for the heat diffusion equation, at a small normal distance from the isosurface δf , called the stagnant film thickness. Although the mesoscopic model, as will be shown later, is restricted to only small supercoolings and small interface velocities, the phase-field approach is considered suitable to be used for the interface propagation. In the present study, a computational model for freezing of supercooled water is presented within the standard cellcentered finite-volume CFD framework [25]. The model utilizes a phase-field diffuse interface model to track the liquid-solid (ice-water) interface. Energy equations in liquid and solid phases are expressed in terms of temperature. The interface points are reconstructed, and the solution for the temperature is decoupled at the interface, by using immersed points, and imposing the melting temperature as a moving boundary condition. The velocity of propagation of the interface is determined by the local interfacial energy balance, given by the sum of the normal-to-theinterface heat fluxes at both sides. With the calculated velocity, the diffused interface is explicitly evolved in time and space. The capability of the model is confirmed by calculating two-phase solidification with known analytical solutions. The model is further extended to couple the conjugate heat transfer in the solid substrate, comprising a framework which can be used to incorporate further models for the interaction of the solidifying supercooled water with the underlying solid substrate. The model enables computation of the solidification of supercooled water as an isolated phase-change process. The proper coupling with conjugate heat transfer enables computation of the interaction with the solid substrate and quantification of its thermal influence on the icing outcome. The previous findings related to the difficulties in the singularity at the contact line can be tested within the finite-volume framework for continuum mechanics in conjunction with the phase-field model, and the model can be used as a tool for further modeling attempts to overcome this difficulty. It should be noted that the model is aimed at, but not restricted to, computing the icing of
AN US
M
80
ED
75
PT
70
CE
65
AC
60
capabilities, especially if the conjugate heat transfer problem associated with the substrate is of importance. Investigations of cases of supercooled water freezing in115 contact with a solid substrate, revealed an interesting observation of a very thin ice film developing and traveling over the substrate surface, prior to the dendritic solidification [3, 4, 6]. The velocity of the tip of this ice layer was determined to be higher than that of dendritic solidi-120 fication in free liquid, depending strongly on the degree of supercooling and the substrate material. A similar phenomenon was studied in [16] by investigating the development of a molten contact line, which forms when a droplet of microcrystalline wax spreads and solid-125 ifies on a cold solid substrate, made of the same material. It was found that the advancing molten contact line approaches an apparent dynamic contact angle, which was modeled as a function of the properties of the substrate material and thermal conditions. This was used in a re-130 lated study in [17], where an attempt was made to compute the spreading and solidification of the molten material over a cold, solid substrate made of the same material. It was shown however, that the conventional continuum formulation for the problem has no meaningful solution. The135 contact angle of the solidified material film at the substrate is determined by the heat flux at the contact line, where the heat flux in the mathematical model has a singularity. The singularity is attributed to a discontinuity of the boundary condition for the energy equation at the140 contact line, which is also present in the case of pure heat transfer in a wedge geometry with sharp corners [18]. An attempt was made in [17] to remove this singularity, by estimating an ad-hoc cut-off length scale for the unknown physical mechanism, which would limit the heat flux and145 remove the singularity at the contact line, in order to yield agreement with the experimental data for a specific material. However, the length scale was obtained artificially, it was not physically justified and not universal, but rather case-specific, thus not providing insight into the underly-150 ing physical mechanism. In [19] an extended systematic computational study was conducted in order to overcome the singularity. A finite element method was used with several independent approaches: different forms of the predefined interface shape;155 different contact angles; imposing an overall heat transfer coefficient as a boundary condition at the interface; allowing temperature change along the interface; and using an artificial resistive layer at the substrate surface. However, again, none of these approaches solves the problem. It160 should be noted that the problem in [16, 17, 19] is simpler, since the solidifying material is the same as the substrate, compared to freezing of supercooled water, which is in contact with the substrate of different kind. In the present study, the phase-field approach is used165 for tracking the ice-water interface. The idea for using the phase-field is borrowed from [20, 21], where this approach was used to simulate the solidification of a supercooled melt. The interface was defined as the phase-field
3
ACCEPTED MANUSCRIPT
175
180
185
supercooled water, other fluids may be used as well. The205 local interfacial energy balance is applied by utilizing a quadratic interpolation algorithm for the temperature in probe points located at small distances from both sides from the reconstructed interface, thus allowing the effects from the surrounding temperature field in the evaluation210 of the local heat fluxes normal to the interface to be included. With regard to the methodology and physics, to the authors’ knowledge such a computational model does not exist in the present literature. Although it might be relevant only in regions of extremely small scales and ex-215 tremely large interface curvatures, the possibility to include the capillary undercooling is also readily available in the model, by simply replacing the melting temperature with the curvature-dependent one in the moving boundary condition. 2. Description of the model
100,00
ED
M
195
Ivantsov -5 stgf, δf=10 m -4 -3
ρl cp,l
-2
0,1
1 ∆T, K
∂T = ∇ · (kl ∇T ), ∂t
225
230
10 235
Figure 1: The dependence of the predicted front propagation velocity on the supercooling at various values for the stagnant-film thickness.
For small supercoolings, say ∆T < 1 K, the theoretical velocity can be predicted by using different values240 of the stagnant-film thickness, in which the temperature drop is a different portion of the global supercooling, as shown in Fig. 1. For higher values of the supercooling, the stagnant-film thickness must be lowered under the (larger) 4
(1)
∂T = ∇ · (ks ∇T ), (2) ∂t where the subscripts l and s refer to liquid (water) and the solid (ice), respectively. The symbols ρ, cp and k denote the density, specific heat and conductivity, respectively, and T is the temperature. The temperature at the interface between the solid and the liquid phase is the melting temperature Tm , and the local energy balance in the normal direction to the interface, for the two-phase moving boundary problem with ρs < ρl , is given by [26] ρs cp,s
ρs Lvn = ks ∇Ts − kl ∇Tl ,
stgf, δf=10 m
AC
0,01
stgf, δf=10 m
CE
v, cm/s
1,00
PT
stgf, δf=10 m
200
2.1. The basic phase-field model for solidification In the solidification of supercooled water, in the absence of bulk motion of the liquid water, the energy diffusion equations in the solid and liquid read
AN US
190
The aforementioned mesoscopic model uses the local temperature drop near the interface to determine its propagation velocity, and is applicable to small values of the global supercooling and small velocities. To illustrate this, the dependence of the predicted front propagation velocity on the supercooling and the stagnant film thickness220 is shown in Fig. 1. The values for the velocity were obtained by numerically inverting the analytical expression from the stagnant-film model coupled with the selection criterion, as described in [21], and using different values for the stagnant-film thickness, and constant thermophysical properties of water at 273.15 K.
mesoscopic scale down to microscopic scales, and the computationally obtained temperature field would have to be resolved on much smaller scales, which is what the mesoscopic model actually tries to avoid. In a recent study [3], utilizing a Hele-Shaw cell to investigate two-dimensional freezing of supercooled water in contact with the underlying solid substrate, the velocity of the ice-layer tip at the wall was found to be much higher than that of the dendritic front far from the wall, thus rendering the analytical stagnant-film solution for the solidification unsuitable. In spite of that, the interface propagation utilizing the phase-field approach provides a basis, that can be coupled to other phase change models and used to calculate the freezing.
CR IP T
170
(3)
where L is the latent heat of melting and the temperature gradients are evaluated on both sides of the interface. The capillary undercooling for curved interfaces can be accounted for with the Gibbs-Thomson relation, which has effects only for nanoparticles with very large curvatures [27] and is neglected here. For most materials the density of the liquid phase is usually smaller than that of the solid one, water being an exception in which the density of the solid phase is lower than that of the liquid one. In such a case, the motion of the phase-change interface is also affected by the volume expansion upon solidification. However, its relative influence on the freezing front velocity is determined by (ρwater /ρice − 1) [26], which has a relatively low amount of about 9% at T = 0◦ C, and decreases with decreasing temperature, since the liquid density decreases, while the ice density increases with decreasing temperature. Therefore, the volume expansion is of minor importance and will be neglected in order to simplify the analysis. In order to propagate the interface with the obtained normal velocity, the phase-field approach from [28] is used
ACCEPTED MANUSCRIPT
in the present study. The sharp interface between the solid and the liquid is represented by a thin transitional region, across which the indicator φ smoothly changes from the value of 1 to 0. The distribution of φ is obtained by solving275 the equation h = b ∇2 φ −
∂φ + vn · ∇φ = ∂t φ(1−φ)(1−2φ) (w/2)2
i ∇φ − |∇φ|∇ · − |∇φ| . (4)
The change of the phase-field variable φ from 1 to 0 in the diffused interface region is initially defined as the hyperbolic tangent profile n i 1h dist 1 − tanh , (5) φn = 2 w
obtained, the temperature gradients are calculated in the normal direction on both sides of the interface. The propagation velocity of the interface is determined from Eq. (3), and the new distribution of the phase-field is calculated explicitly from Eq. (4). The stabilization factor b in Eq. (4) influences the accuracy of the phase-field evolution and must be set using the stability limit blim that depends on the time step size and the mesh resolution [22, 28, 29]. In the present model, the time step size is self-adjusted during the simulation and the parameter b is set as
CR IP T
245
b = 0.05blim ,
blim =
M (6)
The term on the right-hand-side of Eq. (4) represents285 a stabilization term which acts towards ensuring the hyperbolic tangent profile of the phase-field, as discussed in [28, 29]. It is derived by decomposing the normal interface velocity into parts which are independent of and proportional to the interface curvature. This allows the curvature290 to be expressed as a function of the phase-field, which is given by the first two terms in the square brackets, while the last term in the square brackets is the counter-term that cancels out curvature-driven motion at leading order for the interface of finite thickness. In the case of a flat295 interface, the third term in the square brackets will vanish and the first two terms will yield the hyperbolic tangent profile in the stationary solution. During interface motion in absence of external advection, the mutual effects of the left- and right-hand sides of Eq. (4) relax the phase field300 distribution to a hyperbolic tangent profile across the interface and act towards sustaining the profile during the interface motion. The basic model consists of the equations (1)–(4). The solution algorithm starts by solving the energy equations305 (1)–(2). After the solution for the temperature has been
260
265
270
AC
255
CE
PT
250
∇φ . |∇φ|
ED
n=−
5
(8)
where the index t refers to the value in the current time step, index t − ∆t refers to the value from the previous time step, and ndim takes into account the dimensionality of the problem, being equal to 2 for two-dimensional and to 3 for three-dimensional cases, respectively. The maximum Courant number is based on the interface propagation velocity and is calculated as |vn |∆t . (9) Comax = max |d|
AN US
with ndist being the normal distance from the middle of the profile at which φ = 0.5. The characteristic width w, required to achieve sufficient accuracy in representing the hyperbolic transition, depends on the resolution of the finite-volume numerical mesh, and is calculated as w ≥ √ 2|d|, where d is the spacing between two adjacent cell centres in the numerical mesh. It corresponds to the cellsize ∆x in the case of a uniform hexagonal mesh. Thus, the diffused interface is spread across at least three cells at both sides of the iso-surface φ = 0.5. The solid phase is identified by 0.5 ≤ φ ≤ 1, whereas the liquid phase is determined by 0 ≤ φ < 0.5. It should be noted that the material properties in all equations are set according to this distinction. Using the phase-field indicator, the direction of the unit normal to the interface is evaluated280 from the expression
|d|2 1 [1 − 2Comax,t−∆t ] 2ndim ∆tt−∆t
(7)
For the solution of Eq. (4), the velocity should be calculated at all cell-centres. However in practice, it is enough to calculate the velocity only in cells within the transitional band between φmin = 0.01 and φmax = 0.99 [22]. The value 0.05 used in the formulation of the parameter b is smaller than the one recommended in [29] for similar reasons explained in [22], since this recommendation was obtained for a uniform normal interface velocity. In the present work, the problem is much more complex, since the interface deforms due to spatially and temporally nonuniform velocity. The influence of b on the solution and the interface shape was quantified in [22] by varying the ratio b/blim in a wide range. The variations of the width of the diffused interface over the range 0.05 < b/blim < 0.4 were found to be around 2.5%, but it was observed that a value for b higher than 0.05 constrains the interface shape too much. The present model is implemented in the open-source software foam-extend [30], an extended version of the R widely used software OpenFOAM [31]. The divergence and gradient terms in the transport equations are discretized using the Gaussian theorem, and temporal integration is Euler implicit. The new distribution for the phase-field is calculated explicitly, using old-time values for φ. The present model is similar to that in [22]; however, there are some major differences: instead of the analytical solution, the explicit energy balance at the interface is used to model mass transfer; transport equations for heat diffusion in both phases are used, instead of a conservation equation for the solute in the liquid phase; the interface
ACCEPTED MANUSCRIPT
310
points are calculated for the cells at both sides of the interface; the heat diffusion equations are discretized and solved as a single set of algebraic equations, with a proper decoupling at the interface, provided by adding the contributions to the coefficients in the solution matrix. These details are explained in the following subsections.
where the indices 1 and 2 indicate the edge start and end points, respectively. Thus, a set of the edges is detected that is connected and crossed by the interface, as shown in Fig. 2. The position vectors of the points where the interface crosses each mesh edge are calculated from xint,i = x1 +
2.2. Reconstruction of the interface points
320
xp = xP − δint ntri ,
with the same meaning of indices 1 and 2 as above. In order to accelerate the searching process, the width of the local transition region in which the search is confined is estimated as (1 − φint )φmin , (13) rmax = 1.25w ln φint (1 − φmin )
where the factor 1.25 is used to slightly extend the analytical value obtained by using the inverse of the hyperbolic tangent profile, Eq. (5). Based on that, a set of the nearest iso-surface points xint,i to each cell-centre xP within the band is found by calculating the distances between the cell-centre and the interface points lying within rmax
AN US
(10)
where xP are the position vectors of the cell-centres within the [φmin , φmax ] transition band, and ntri is the local unit normal to the interface. To calculate xp , the appropriate distances δint = xP − xp from all the cell-centres within the transition band to the corresponding points at the interface in the normal direction must be determined.
M
P
d int
fmin
PT
p
ED
ntri
interface
dP i = |xint,i − xP | < rmax .
CE
x int,1
x int,2
fmax
fint 325
Figure 2: Illustration of the algorithm for finding the points at the interface. The identified mesh edges crossed by the iso-surface are denoted by dashed lines.
The algorithm for finding the points at the interface,330 defined by φint = 0.5, is based on detecting the intersections of the mesh-edges with the interface. First the values for the phase-field variable φ stored at cell-centres are interpolated to the nodes of the cells in the mesh by using the least squares volume-to-point interpolation described335 in [32]. The edges of the mesh intersecting the interface, are extracted using the condition (φ1 − φint ) · (φ2 − φint ) ≤ 0,
(11) 6
(14)
From the set of the nearest interface points, the three nearest interface points are extracted with the smallest distances dP i from the cell-centre of interest, xP . Using the three interface points, a local unit normal to the interface is determined from the vector products xint,3 − xint,1 xint,2 − xint,1 × , (15) ntri = |xint,2 − xint,1 | |xint,3 − xint,1 |
where the indices 1, 2 and 3 indicate the three iso-surface points with the smallest distances dP i . Finally, the distance between the cell-centre of interest to the corresponding interface point is calculated by projecting the nearest distance dP i to the direction of the local unit normal, mapped to the centre of the cell P δint = |(xint,1 − xP ) · ntri |.
AC
315
(12)
CR IP T
The energy balance at the interface must be applied to calculate the velocity in Eq. (3) for all cells within the band 0.01 ≤ φ ≤ 0.99. Therefore, the corresponding interface points for all cells must be found in order to calculate the heat fluxes via the temperature gradients at both sides of the interface. During calculation, the interface points and the temperature gradients are mapped to all cells within the interfacial band. Referring to Fig. 2, the locations of the points p at the interface are determined from
φint − φ1 (x2 − x1 ), φ2 − φ1
(16)
If the cell P is in the solid, φP ≥ 0.5, the value δint is multiplied by -1, thus representing the signed normal distance from each cell-centre P within the band to the interface, being negative for the cells in the solid phase and positive for the cells in the liquid phase. The described procedure is performed for each cell within the predefined band φmin ≤ φ ≤ φmax on both sides of φint and the corresponding interface point is calculated from Eq. (10) for each of those cells. Currently, the described model is two-dimensional, and therefore only three points are necessary to determine the local iso-surface segment within the cell, according to R and in foam-extend, the mesh Eq. (15). In OpenFoam , is always three-dimensional, having at least one layer of cells. Thus, the mesh has a ’thickness’ of one cell, and simulations in one or two dimensions are enabled by suppressing the boundary conditions in the other dimensions. Therefore, the third required point xint,3 always lies at the edge at the opposite side of the cell.
ACCEPTED MANUSCRIPT
(17)
M
P
N
365
d
immersed points
ip N
x ip
x
fint
d
xP
Figure 3: Immersed points used in the discretization of the heat fluxes in the energy equations for setting the melting temperature Tm . The gray colored cells are identified as ’solid’ and the white colored cell as a ’liquid’ cell.
For a one-dimensional geometry, the change of φ around the ghost cell-face between the adjacent cell-centres P and N can be approximated as linear and, considering for example the cell P, the distance x can be calculated as φP − φint x= 1− |d|. (18) φP − φN However, for a two-dimensional geometry, the change of φ around the ghost cell-faces between adjacent cells can be arbitrary and therefore, x is approximated using a cubic interpolation a0 + a1 x + a2 x2 + a3 x3 = 0,
(20)
d · ∇φP , |d|
(21)
a2 =
3(φN − φP ) − d · (∇φN + 2∇φP ) , |d|2
(22)
a3 =
d · (∇φN + ∇φP ) − 2(φN − φP ) . |d|3
(23)
CR IP T
a1 =
The distance between the cell of interest P and the immersed point ip is then |xP − xip |, and similarly for the neighboring solid cell N |xN − xip |. In order to account for the melting temperature in the immersed points, referring to the upper two cells P and N in Fig. 3, the heat flux through the immersed cell-faces is discretized in the energy equations as kP
TP − Tm Sif , |xP − xip |
(24)
for the flux entering the cell P from the left, and
immersed faces Sif ip
a0 = φP − φint ,
kN
ghost faces
interface
with the following coefficients
AN US
xip
d = xN + x . |d|
ED
355
PT
350
CE
345
2.3. Decoupling of the energy equations For the discretization of the energy equations (1)–(2) in the liquid and the solid phase, an approach utilizing immersed points ip and immersed cell-faces is used. The ’liquid’ and ’solid’ cells are identified as those for which φ < φint and φ ≥ φint , respectively. The melting temperature Tm is explicitly set in the immersed points, which are mapped to the intersected cell-faces. The intersected (referred to as ghost) cell-faces are those, for which one adjacent cell sharing that face is identified as liquid, and the other as solid (faces with dashed edges in Fig. 2). The term ghost face originates from the fact that the discretized heat flux through such a cell-face is eliminated in the dis-360 cretization of the energy equation and replaced by the corresponding heat flux through the immersed faces. This enables to discretize the heat fluxes entering the interfacecrossed cells from the liquid-solid interface; hence, the immersed points have the purpose of acting as an immersed moving boundary at which Tm is set. The position vectors of the immersed points, shown in Fig. 3, are determined from
AC
340
(19) 7
Tm − TN Sif , |xip − xN |
(25)
for the flux entering the cell N from the right, with Sif being the area of the immersed cell-face equal to that of the ghost cell-face between the cells P and N. Thus, the corresponding contribution to the diagonal coefficient in the solution matrix for the temperature for the ’liquid’ cell P is kP Sif /|xP − xip |, and the source-term contribution is kP Tm Sif /|xP − xip |, and similar for the ’solid’ cell N. In the calculation of the velocity, it is important to evaluate correctly the temperature gradients at the interface, for which a more accurate approximation is required than a simple interpolation of the gradients from the surrounding cells. For a one-dimensional geometry, the temperature gradients at the interface could be calculated directly, by using the calculated temperatures from the cellcentres adjacent to the interface and their distances to the interface points. However, if the geometry is not onedimensional, this will not be accurate and the gradients will not be smooth enough along the interface, as shown in [33]. Therefore, for calculation of the temperature gradients at the interface in the present study, an additional set of points, referred to as probe points pp, is placed in the normal direction to the interface on both sides and at the same distance δpp . With reference to Fig. 4, the locations of the probe points for a cell P within the interfacial band are calculated xpp,l = xp + δpp nP ,
(26)
xpp,s = xp − δpp nP ,
(27)
ACCEPTED MANUSCRIPT
for the liquid and solid side respectively, with nP being the unit normal calculated using Eq. (6). The probe points are used to interpolate the temperature from the surrounding cell-centres and to calculate explicitly the normal-to-the-interface temperature gradients with the following expressions ∇Tl =
Tpp,l − Tm n, δpp
(28)
∇Ts =
Tm − Tpp,s n, δpp
(29)
weighted least squares method TP 0 − Tm c0 c1 ... c2 = [(XT ·W·X)−1 ·XT ·W]· TP j − Tm c3 ... TP n − Tm c4
L
interface
S
S
L
L
x pp,l
S
S
S
L
S
S
S
ppl p pps
xp
x pp ,s
dpp
390
395
M
dpp
ED
p – interface point, ppl/s – probe points in liquid/solid L/S – liquid/solid interpolation cells stencils
400
Figure 4: Probe points near the interface used for the calculation of the temperature gradients at the interface. The gray colored cells comprise the interpolation stencil for the probe point at the solid side, and the white colored cells are for the liquid side.
385
PT
If the temperature gradients at the interface would be calculated directly using the temperatures in the adjacent cell-centres, the region around the interface should be fully thermally resolved in order to obtain a good approximation for the temperature gradients. The thickness of the thermal diffusion layer in the normal direction to the interface is estimated as δt ∼ αl /vn , with αl being the thermal diffusivity of the liquid, and vn the magnitude of the interface velocity. In order to estimate the required spatial resolution, the velocity vn must be known. For demonstration, the velocity can be related to the thermal conditions from the steady state energy balance using a coordinate system positioned at, and moving with the interface with the constant velocity equal to vn . For the case of a onedimensional one-phase solidification of a supercooled liquid, shown schematically in Fig. 5, the energy equation in the liquid has the form
The distance δpp , as the input parameter, should be chosen such that the probe points are close enough to the interface to obtain a good approximation of the temperature gradients, but not too close in order to avoid numerical instabilities. The distance can be set to a portion of the cell size √ δpp = C∆x, where C ≤ 0.5 for one dimension, or C ≤ 0.5 2 for two-dimensions. The interpolation of the temperature in the probe points is performed using a second-order quadratic spatial function, defining the distribution of the temperature T over a region around the interface point p, which has the form
CE
380
3. Results and discussion
AC
375
which is first multiplied by XT · W from the left side, and afterwards again by (XT · W · X)−1 from the left side. The weights in the diagonal weighting matrix W are calculated based on the distance R of each cell-centre in the interpolation stencil to the interface point p, from the expression 0.5{1 + cos[(πR)/(1.1Rmax )]}/lref , where Rmax is the maximum distance from the cell-centres in the interpolation stencil, and the reference length is taken as lref = 2∆x. The interpolation stencils at the liquid and solid sides consist of two layers of cells closest to the interface point, although additional layers may be added. Upon obtaining the temperatures in the probe points Tpp from Eq. (30), the gradients at the interface are explicitly calculated using equations (28)–(29) and substituted to Eq. (3) to calculate the velocity, mapped to each cell within the diffused interface band.
AN US
L
, (31)
where X is the n × 5 matrix whose row j consists of the following entries: (xP j − xp ), (yP j − yp ), (xP j − xp )(yP j − yp ), (xP j − xp )2 , (yP j − yp )2 , and n is the number of the neighboring cells used in the interpolation stencil. Equation (31) follows from TP 0 − Tm c0 c1 ... T − T c (32) = X· m , 2 Pj c3 ... TP n − Tm c4
at the liquid and solid sides respectively. The distance to the probe points δpp = |xp − xpp,l/s | is equal for the liquid and the solid side of the interface.
L
CR IP T
370
T (x, y) = Tm + c0 (x − xp ) + c1 (y − yp )
+c2 (x − xp )(y − yp ) + c3 (x − xp )2 + c4 (y − yp )2 (30) where Tm is the melting temperature at the interface point p, and (xp , yp ) are the coordinates of that point. The unknown coefficients ci can be expressed in terms of the temperature TP j in the surrounding cell-centres using
−vn 8
∂T ∂2T = αl 2 . ∂x ∂x
(33)
ACCEPTED MANUSCRIPT
274 theory
273
Tm
-3 ∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
-vn
T, K
Tf vn
x
271 270 269 268
Figure 5: Sketch for the one-dimensional energy balance.
267 0
Integrating Eq. (33) once, the distribution of the temperature gradient in the normal direction to the interface is obtained ∂T vn ∂T = exp − x , (34) ∂x ∂x 0 αl
274 273
271
0,02
AN US
T, K
272
CR IP T
ice
272 water
where (∂T /∂x)0 is the temperature gradient at the interface x = 0. Integrating again, the temperature distribution in the normal direction is obtained ∂T αl vn T = T∞ − exp − x , (35) ∂x 0 vn αl
ED
M
420
425
0,08
0,1
theory -3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
270
267 0
In order to calculate the interface temperature gradient, a one-dimensional numerical setup was used, with one boundary declared as the interface and having a fixed temperature equal to Tm = 273.15 K, and at the opposite boundary the temperature T∞ = 267.15 K is set, yielding430 a global supercooling of 6 K. The velocity in the domain is constant and equal to 0.01 m/s, which corresponds to the experimentally obtained value in [34] for a supercooling of 6 K. The analytical value for the temperature gradient at 435 5 the interface equals −( ∂T ∂x )0 = 4.65 · 10 K/m, at x = 0. Equation (33) is solved using the standard finite-volume discretization in the software foam-extend, with the following values for the properties of water at ∆T = 6 K: kl = 0.5487, W/(mK), ρl = 999.08, kg/m3 and cpl = 4246 J/(kgK). With these values, the thermal diffusivity of water is αl = 1.29 · 10−7 m2 /s, and the estimated thickness of the thermal diffusion layer is δt = 12.9 · 10−5 m. The computed results for the temperature and the temperature gradient in the normal direction to the interface are shown in Fig. 6 and Fig. 7, where the results at the bottom of the figures are the zoomed profiles near the interface. The obtained profiles for T (x) and ∂T ∂x (x) approach the analytical solutions, Eq. (35) and Eq. (34), respectively,
CE
PT
415
0,06
268
AC
410
x, mm
269
where T∞ is the far-field temperature at x → ∞. Using ∆T = T − T∞ , an analytical relation is obtained between the global supercooling and the temperature gradient and the front propagation velocity at the interface (x = 0) ∂T ∆T vn − = . (36) ∂x 0 αl 405
0,04
9
0,01 x, mm
0,02
Figure 6: The computed temperature in the normal direction to the interface: temperature profile (top) is zoomed near the interface (bottom).
as the thermal diffusion layer in the normal direction to the interface becomes better resolved. From the results in Fig. 7 it follows that the thermal diffusion layer has to be resolved with more than 20 computational cells in the normal direction to the interface. At higher velocities, or higher degrees of supercooling, the thickness of the thermal diffusion layer will decrease, and the cell size required to resolve this layer would have to be much smaller. For a geometry which is not one-dimensional and at higher velocities, this leads to a significant increase of the computational effort. However, when the temperature gradients are calculated with the previously described procedure using the interpolated temperatures to the probe points, the required mesh resolution is expected to be lower. In order to test the capabilities of the model using the least square interpolation of the temperature in the probe points, the onedimensional two-phase Stefan problem is used. The liquid is initially supercooled and, at the initial time, the temperature at one boundary suddenly changes to Tm , triggering solidification, upon which the interface starts moving. The problem is described by equations (1)–(3), and the analyt-
∇T, K/m
ACCEPTED MANUSCRIPT
√ α vn = λ √ s . t
0 theory -3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
-1e+05 -2e+05
440
-3e+05
(41)
In the numerical setup, the length of the mesh in the horizontal direction, normal to the interface, is 10 mm, with an initial supercooling of 20 K obtained by specifying fixed temperatures Ts0 = Tl0 = 253.15 K at the opposite boundaries. The properties of water at the given supercooling, and that of ice, are listed in Table 1, and the latent heat of melting is L = 333000 J/kg.
CR IP T
-4e+05 Table 1: Thermophysical properties of water at ∆T = 20 K and ice.
0,02
0,04
x, mm
0,06
0,08
0,1
water ice
0 theory
-1e+05 -2e+05
445
-3
∆x = 2.10 mm -3 ∆x = 1.10 mm -3 ∆x = 0.5.10 mm
450
-3e+05 -4e+05
455
-5e+05 0
0,02
M
0,01 x, mm
PT
ED
Figure 7: The computed temperature gradient in the normal direction to the interface: profile of the gradient (top) is zoomed near the460 interface (bottom).
ical solution is given by
CE
Tl = Tl0 + (Tm − Tl,0 )
erfc ( 2√xα t )
Ts = Ts0 + (Tm − Ts,0 )
k, W/(mK) 0.5263 2.22
ρ, kg/m3 996.32 916.2
cp , J/(kgK) 4318 2050
If the temperature at the left boundary suddenly changes to Tm , solidification starts theoretically from the left boundary and the interface moves from left to right. Since the phase-field is initialized using Eq. (5), the initial position of the interface xinit cannot be at the left boundary (in that case the velocity would be infinite). Instead of that, the interface φ = 0.5 is placed at a small distance equal to 6.5w from the left boundary, in order to obtain a suitable initial φ–profile. The theoretical initial temperature distribution according to equations (37)–(38) is prescribed in the liquid and the solid, respectively. Thus, the initial position of the interface has a time shift equal to tinit = x2init /(4λ2 αs ), which is accounted for in the output values from the computation. The computed interface position and interface velocity are shown in Fig. 8, and the computed temperature distribution normal to the interface is shown in Fig. 9, all compared to the corresponding analytical solutions. The position of the interface is given by the maximum position vector of the cross-section of the interface with the mesh edges, and the velocity is taken directly from the cells within the interfacial band, since they are all equal in a time step for the one-dimensional case. As expected, using the interpolation of the temperature in the probe points allows a coarser mesh resolution. The thermal diffusivity of water is α = 1.22 · 10−7 m2 /s, and the thermal diffusion layer increases from a lower value at the beginning to δt ∼ 10−3 m at later times, as the interface velocity decreases. The coarsest mesh has about 4 times less cells in the thermal diffusion layer, which is sufficient for a high accuracy of the calculated interface temperature gradients. The computed interface position approaches the analytical solution with increasing mesh resolution, as the upper diagram in Fig. 8 shows. The relative error in the computed interface position is evaluated using the expression
AN US
∇T, K/m
-5e+05 0
l
erfc (νλ)
erf ( 2√xαs t ) erf (λ)
465
,
(37)
,
(38)
AC
where ν = αs /αl , and Tl0 and Ts0 are the temperatures in the liquid phase at x → ∞ and in the solid phase at x = 0, respectively [26]. The characteristic root λ is obtained by numerically inverting the following expression Stl Sts + eλ2 erf (λ) νe(ν 2 λ2 ) erfc (νλ)
√ = λ π,
470
475
(39)
where Sts = cp,s (Tm − Ts,0 )/L and Stl = cp,l (Tm − Tl,0 )/L are the characteristic Stefan numbers for the solid and the liquid, respectively. The position and the velocity of the interface are given by the expressions √ x = 2λ αs t, (40)
ε=
|xi − xth | · 100%, xth
(42)
where xi is the computed position of the interface obtained with different cell sizes and xth is the analytical solution. 10
ACCEPTED MANUSCRIPT
4,0 Table 2: The computed interface position and the relative error at t = 30 s.
i 1 2 3
2,0
Analytical -3 ∆x = 200.10 mm -3 ∆x = 100.10 mm -3 . ∆x = 50 10 mm
1,0
0
10
0,2
30
ε=
Analytical -3 ∆x = 50.10 mm δpp = 0.1∆x 485
0,1
0,1
Table 3: The computed interface position and the relative error in the temperature profile at t = 5 s.
5,0
10,0
15,0 t, s
20,0
25,0
i 1 2
30,0
Figure 8: The computed interface position (top) and magnitude of the interface velocity (bottom) vs. time for different cell sizes.
Analytical δpp = ∆x
275
PT
δpp = 0.1∆x
AC
255
250 0,0
495
simulation time = 5 s -3 ∆x = 50.10 mm
CE
T, K
270
260
490
ED
280
265
1,0
2,0
x, mm
δpp ∆x 0.1∆x
xi , mm 1.4488 1.4589
ε, % 1.64 0.96
The performance of the model in two dimensions is demonstrated using the setup with a geometry having dimensions of 1 mm × 1 mm, discretized with an equal number of cells in the horizontal and vertical directions, and having adiabatic boundaries on all sides. The phase-field φ is initialized according to Eq. (5) as a circle with the centre located at the centre of the mesh, the initial supercooling is ∆T = 6 K everywhere, the properties of water and ice are listed in Table 4, and the latent heat of melting is L = 333000 J/kg. Table 4: Thermophysical properties for water at ∆T = 6 K and ice.
3,0
4,0
water ice
5,0
Figure 9: The computed temperature distribution normal to the500 interface for different distances of the probe points.
480
(43)
where xδpp ,i is the computed interface position in the temperature profile obtained with different distances from the probe points and xth,T = 1.473 mm is the analytical solution at time t = 5 s. The values of the relative error are given in Table 3, which shows that the result is slightly better if the probe points are placed closer to the interface.
M
0,0 0,0
ε, % 1.08 0.94 0.68
|xδpp ,i − xth,T | · 100%, xth,T
AN US
vn, cm/s
0,2
20
t, s
xi , mm 3.3908 3.3958 3.4049
The results in Fig. 9 show that choosing a larger distance of the probe points from the interface, comparable to the cell size, yields slightly less accurate result. The relative error for the computed interface position in the temperature profile at time t = 5 s is evaluated using the expression
δpp = 0.1∆x 0,0
∆x, mm 50·10−3 100·10−3 200·10−3
CR IP T
x, mm
3,0
The values of the computed interface position and the corresponding relative error obtained with the analytical so-505 lution xth = 3.428 mm at t = 30 s are given in Table 2, which shows that the relative error decreases with increasing mesh resolution. 11
k, W/(mK) 0.5487 2.22
ρ, kg/m3 999.08 916.2
cp , J/(kgK) 4246 2050
The event of freezing along a solid substrate is very fast and finishes within only several milliseconds at a distance of a few millimeters [3]. In the present study, the simulation time is set to 1 s and the computed final temperature field is shown in Fig. 10 for three mesh resolutions, with the indicated φ = 0.5 iso-contour. With the coarsest mesh, the interface starts developing as a circular shape that tends to evolve to a rectangular shape at later times, as can be seen in the top snapshot in Fig. 10. This is attributed to
ACCEPTED MANUSCRIPT
ED
M
AN US
CR IP T
the fact that the thickness of the thermal diffusion layer therefore, there is a greater total temperature difference quickly increases, as the interface velocity decreases. 520 between Tm at the interface and the temperature at the corners. Consequently, this leads to higher temperature gradients at the interface and higher interface velocity, according to Eq. (3), leading to the faster progression of the interface along the diagonals towards the corners. By in525 creasing the mesh resolution, the interface shape is closer to being circular, as can be seen in the middle and the bottom snapshots in Fig. 10. This can be seen more clearly in Fig. 11, in which separate final iso-surface contours are plotted with different 530 colors for different mesh resolutions used in the simulations.
535
AC
CE
PT
Figure 11: The final iso-contour shapes for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
Figure 10: The computed final temperature field for different mesh540 resolutions: 50×50 cells (top), 100×100 cells (middle) and 200×200 cells (bottom). 510
515
At later times, the thermal diffusion layer reaches the domain boundaries, and the iso-contour shape starts de-545 forming. The deformation occurs in the diagonal directions, since the temperature distribution is steeper along the diagonal, compared to the horizontal and vertical directions. Since the distance from the interface to the corners of the mesh is greater in the diagonal directions, it550 takes more time for heat to diffuse in that directions, and 12
The computed temporal evolution of the interface and the temperature field at early times is shown in Fig. 12 for the finest mesh resolution with 200×200 cells. As expected, due to the higher conductivity and thermal diffusivity of ice, the temperature within the ice reaches the melting temperature Tm at a shorter time, whereas the temperature in the liquid water requires a longer time to increase. The reconstruction of the interface points is demonstrated in Fig. 13, which shows the computed cross-section points of the iso-surface with the mesh edges, according to Eq. (12), and the calculated interface points, according to Eq. (10). It can be seen that the interface points are located between the cross-section points of the interface with the mesh edges. The zoomed detail in the bottom snapshot in Fig. 13 shows this more clearly. The diffused interface band contains several cells in the direction normal to the interface, and the corresponding interface point needs to be found for all cells within the band in both phases (in order to apply the condition in Eq. (3)). This is why sev-
ACCEPTED MANUSCRIPT
CR IP T AN US M ED CE
PT
Figure 13: The computed cross-section points of the iso-surface with the mesh edges (Eq. (12), black symbols) and the interface points (Eq. (10), red symbols) for the mesh with 200×200 cells. The bottom snapshot is a zoomed detail from the top snapshot.
AC
555
eral interface points may be located in each segment of the interface between two consecutive cross-section points of the iso-surface with the mesh edges. They are mapped to the corresponding cells within the band, depending on the orientation of the interface according to Eq. (10).
560
Figure 12: The computed temporal evolution of the interface and the temperature field for the mesh with 200×200 cells.
13
In order to quantify the model performance, the interface velocity is evaluated twofold, as the velocity of a single interface point, and the mean velocity within the diffused interface band. First, the velocity of the interface point moving along the horizontal direction is calculated from the expression xt − xt−∆t , (44) vn,t = ∆t where xt and xt−∆t are the instantaneous distances of the interface point tracked along the horizontal direction at the time t and at the previous time t − ∆t, respectively, from its initial position. The position of that point is readily available from the model, as the maximum position vector of the cross-section of the iso-surface with the mesh edges. Fig. 14 shows the magnitude of the computed velocity during the simulation time, and it is seen that the veloc-
ACCEPTED MANUSCRIPT
575
580
0,04 0,02
400
t, ms
600
800
1000
600
PT
Figure 14: The magnitude of the computed velocity of the interface point along the horizontal direction for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
800
1000
where P/π is a measure for the mean diameter D of the iso-contour. In the case of an ideal circular iso-contour, D would be equal to its real diameter. Additionally, the same quantity is calculated in another way, by summing all the volumes of the cells (m in total) identified as ’solid’, for which φ ≥ 0.5, and dividing the sum by the mesh thickness, in which case the amount of ice is determined by Pmthe two-dimensional area covered by the solid cells, A = i=1 ∆Ai . The computed temporal evolution of the isocontour perimeter is shown in Fig. 16, and the evolution of the computed ratio AD /A is plotted in Fig. 17.
CE
1,0 P, mm
590
600
1,2
Next, the interface velocity is represented in terms of the mean velocity in the cells, which are located within the diffused interface band in a single time step. The arithmetic mean of the velocity is calculated by summing the velocity magnitudes in all cells lying within the [φmin , φmax ] interfacial band and dividing the sum by the number of these interfacial cells. The computed mean interface velocity magnitude within the band is plotted in Fig. 15, in which it is seen that the mean velocity does not show oscillations. In order to further quantify the predictive capabilities of the model, in each time step the amount of the new formed ice (solidified water) is calculated in two different ways. First, in each time step the iso-contour perimeter is computed n X P = πD ≈ ∆li , (45)
AC
585
t, ms
Figure 15: The computed mean interface velocity magnitude within the diffused interface band for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
595
200
400
crossing the mesh-edges. It is equal to the area of the isosurface divided by the thickness of the mesh. The amount of the solidified water is then approximated as the area in two dimensions 2 π P π 2 , (46) AD = D = 4 4 π
ED
0,00 0
200
AN US
0,06
0,04
0,00 0
M
vn , cm/s
0,08
0,06
0,02
0,10 50x50 100x100 200x200
50x50 100x100 200x200
0,08
CR IP T
570
0,10
ity oscillates slightly around a mean value. These oscillations are attributed to the time instants when the interface passes the cell-centres, i.e. the cell that was previously identified as liquid now becomes solid. For the transient term in the energy equation, the temperature Tm is set as the old-time value in that cell when this occurs. The exact time instant when the iso-surface is located exactly at the cell-centre cannot be determined within a time step; otherwise some kind of temporal extrapolation for the temperature would have to be applied. In addition, the stencils for the interpolation of the temperature in the probe points change, all leading to somewhat non-gradual temporal change of the temperature gradients at the interface and the oscillating velocity. In spite of that, it should be noted that the oscillations are not noticeable in the onedimensional simulation, and the amplitudes of the oscillations are relatively small, also in the two-dimensional case.
vn,mean , cm/s
565
0,8 50x50 100x100 200x200 0,6 0
200
400
t, ms
600
800
1000
Figure 16: The computed temporal evolution of the iso-contour perimeter for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
i=1
by summing all the interface line segments ∆li (n in total)
14
ACCEPTED MANUSCRIPT
1,3
AD / A, -
1,1 1,0
4Ld
symmetry
1,2
adiabatic
adiabatic 50x50 100x100 200x200
0,8 0,7 0
adiabatic
200
400
t, ms
600
800
1000
2Ld
Figure 17: The computed ratio AD /A for different mesh resolutions: 50×50 cells (black line), 100×100 cells (red line) and 200×200 cells (blue line).
Figure 18: Initial configuration for the simulation of the ice dendrite growth.
635
615
630
phase-field φ is initialized according to Eq. (5), where ndist is calculated as the minimum distance from the cell-centres within the phase-field band to the profile of the parabolic dendrite with φ = 0.5. At a supercooling of 10 K, the dendrite tip velocity predicted by the MST is 4.9 cm/s with the corresponding theoretical tip radius of R = 2.45 · 10−7 m. The initial parabolic φ–profile is y = ax2 +bx+c, where b = 0, c = Ld , and the coefficient a = −1/(2R) is obtained from the condition that the radius of curvature is R in the vertex of the parabola located at (0, Ld ). The simulation is run until the temperature at the adiabatic boundaries (right and top) starts changing, indicating that the thermal diffusion layer has reached the boundary. Hence, the assumption of a semi-infinite space in which the dendrite grows is not violated and the steady-state dendrite tip velocity is reached. The computed interface position and velocity of the interface point at the symmetry plane for different mesh resolutions is plotted in Fig. 19. The velocity at a time t is calculated using Eq. (44) by replacing x with the maximum coordinate y of the iso-surface point at the symmetry plane. At the beginning of the simulation, the temperature gradients at the interface are highest; hence, the dendrite tip moves very fast. Shortly after that, the temperature field within both phases develops and the dendrite tip velocity decreases. The computed dendrite tip velocity quickly converges to the theoretical value, as the lower diagram in Fig. 19 shows. The computed interface position converges with increasing mesh resolution, as can be seen in the upper diagram in Fig. 19. The relative error is quantified in terms of the value of the computed interface position obtained with the finest mesh, from the expression
AN US
AC
625
CE
PT
620
M
610
As can be seen in Fig. 16, the plot for the temporal evolution of the perimeter of the iso-contour has the characteristic shape, similar to that obtained in the one-dimensional case in Fig. 8. More important for justification is the ratio AD /A, which converges to unity with mesh refinement, as640 Fig. 17 shows. Finally, the described model is used to compute the growth of an ice dendrite into supercooled water. The transition from planar stable freezing to the unstable dendritic growth is characterized in the marginal stability the-645 ory (MST) [35, 36, 37]. If the solid-liquid interfacial energy is neglected, the product of the velocity of the dendrite tip and its radius at a specific supercooling is determined by the solution of Ivantsov [23] for the growth of a parabolicshaped dendrite in semi-infinite space. When the interfa-650 cial energy is accounted for, the individual radius of the dendrite tip, which is dependent on the given supercooling, is selected by using MST. Basically, the dendrite tip velocity is a single-valued function of the tip radius at a given supercooling, and every selected value of the dendrite tip655 radius determines a corresponding exact steady-state tip velocity. The setup for the simulation is sketched in Fig. 18. The geometry is 3.5 µm ×7 µm, discretized with cells having equal dimensions in horizontal and vertical directions, with660 all adiabatic boundaries, except the left boundary, which is a symmetry plane. The initial supercooling is ∆T = 10 K everywhere, the properties of water and ice are listed in Table 5, and the latent heat of melting is L = 333000 J/kg.
ED
605
CR IP T
Ld
0,9
Table 5: Thermophysical properties for water at ∆T = 10 K and ice.
water ice
k, W/(mK) 0.5392 2.22
ρ, kg/m3 998.13 916.2
cp , J/(kgK) 4272 2050
ε=
|yi − y4 | · 100%, y4
(47)
where yi is the position of the interface point at the symmetry plane obtained in simulations with coarser meshes
The length of the initial dendrite is Ld = 1.75 µm. The 15
ACCEPTED MANUSCRIPT
0,15
expression 50x100 100x200 200x400 400x800
675
0,05
0,00 0e+00
2e-04
4e-04
100
t, ms
6e-04
1e-03
680
60
AN US
theory (MST) 50x100 100x200 200x400 400x800
80
vn , cm/s
8e-04
40 20 0 0e+00
4e-04
t, ms
6e-04
8e-04
1e-03
M
2e-04
CR IP T
y, µm
0,10
y2 − y3 y3 − y4 p= = 1.307, (48) ln 2 where the factor 2 in the denominator is given by the mesh refinement ratio ∆xi /∆xi−1 . The spatial accuracy is not second order, which is expected, since the discretization of the spatial terms in the transport equations is not second order. Finally, the computed interface shape with the developed temperature field and the velocity field at t = 0.001 ms is plotted in Fig. 20 for the mesh with 200×400 cells. Shortly after the beginning of the computation, the temperature has reached the value Tm within the dendrite, and heat is diffused into the liquid, producing parabolic isotherms around the φ = 0.5 iso-contour. The highest computed velocity is in the region at the tip of the dendrite. ln
PT
CE
670
(50×100, 100×200 and 200×400), and y4 is the one obtained in the simulation with the finest mesh (400×800) at the same simulation time. The values of the computed position of the interface point at the symmetry plane at t = 0.001 s, and the corresponding relative error calculated with Eq. (47), are given in Table 6, which shows that the relative error decreases with increasing mesh resolution.
Table 6: The computed position of the interface point at the symmetry plane and the relative error at t = 0.001 s.
AC
665
ED
Figure 19: The computed interface position (top) and magnitude of the interface velocity (bottom) of the interface point at the symmetry plane vs. time for different mesh resolutions.
i 1 2 3 4
mesh 50×100, ∆x 100×200, ∆x/2 200×400, ∆x/4 400×800, ∆x/8
yi , µm 0.083 0.0911 0.0958 0.0977
ε, % 15.05 6.76 1.94 0
Figure 20: The computed temperature field and the velocity field within the diffused interface band for the mesh with 200×400 cells. 685
The values of the relative error in Table 6 indicate a690 convergence rate higher than second order, however this should not be confused with the accuracy of the results. The order of accuracy is estimated using the results for y for the three finest meshes (2, 3 and 4 in Table 6) from the 16
It can be concluded from the results and discussion that the model is capable of predicting the process of solidification of supercooled water with the accompanying thermal and phase-change effects on the moving ice-water interface. This is demonstrated in particular by the correct prediction of the existing analytical results for the Stefan moving boundary problem and the theoretical results for the ice-dendrite growth into supercooled water as the commonly used representative analytical solutions. The
ACCEPTED MANUSCRIPT
(50)
which is internally conserved without the need to explicitly transfer the fluxes at the shared boundary during the computation. This functionality was successfully used previously for the simulation of non-isothermal drop impact [7] and is also utilized here. The effective solidliquid conductivity at the wall cell-faces within the band φmin < φ < φmax , is calculated from the expression kls = φks + (1 − φ)kl , which, combined with the conductivity of the wall material, is then used for the harmonic interpolation of the conductivity at the wall cellfaces. The interpolation stencils in Eq. (31) are extended to include the values for the temperature at the cell-faces at the coupled boundary as well. At the coupled boundary the phase-field is extrapolated with zeroth order from the cells adjacent to the substrate surface, thus imposing a zero-gradient boundary condition. The thickness of the thin ice film developing along the cold substrate was not reported in [3], whereas it was estimated at 20 µm in [6]. In the present study, the propagation of the thin ice film developing along the cold underlying substrate is computed for the initial supercooling of ∆T = 6.25 K. The two-dimensional computational domain shown in Fig. 22 has two regions: the upper for the water and ice and the lower for the substrate. Both have dimensions in the vertical plane of 50 µm × 50 µm, discretized with equal number of cells in the horizontal and vertical directions.
AN US
Modeling the evolution of the thin ice film due to solidification of supercooled water in a thin region over a cold substrate requires solving energy equations in all three phases, i.e. within water and ice, and in the substrate as720 well, in order to account for the conjugate heat transfer at the wall surface. When the underlying substrate, shown in the sketch in Fig. 21, has a much greater thermal diffusivity compared to supercooled water and ice, and its thermal influence on the dynamics of the spreading of the725 thin ice film is assumed to be purely as a heat absorber.
730
vn
substrate (wall)
M
735
Figure 21: Sketch of the thin ice layer due to solidification of supercooled water over a cold substrate.
ED
The additional transport equation to be solved is the energy diffusion equation in the substrate (wall), which reads ∂T ρw cpw = ∇ · (kw ∇T ), (49) ∂t
PT
T =T¥
where the subscript w denotes the wall. The energy equation (49) is coupled with the basic model defined by equations (1)–(4). Considering the discretization of the energy equations, the software package foam-extend provides a functionality to solve a single energy equation for the entire domain, consisting of the mesh for the fluid, where the phase change takes place, and the mesh in the substrate. The coupling between the two meshes makes use of the software’s functionality of connecting the two meshes for a single solution of the temperature field. The discretization of the energy equation follows the standard finite-volume procedure. Special care is given to the cell-face interpolation of the conductivity in the diffusion term at the cell-faces shared by both meshes. According to [38], the cell-face value of the conductivity740 is evaluated using harmonic interpolation. Thus, the energy equation is solved for both meshes simultaneously, by extending the solution matrix for the energy equation to include the mesh for the substrate as well. The coupling between the ice-water and the substrate domains is745 17
T =T¥
water (liquid) ice (solid)
coupled
T =T¥
interface, Tm
CE
710
715
AC
705
−kl/s (∇T |⊥ )l/s = −kw (∇T |⊥ )w ,
CR IP T
4. Extension to conjugate heat transfer
achieved by fulfilling the requirement that the temperature be continuous at the wall surface, and the heat fluxes be conserved in the normal direction to the wall surface
symmetry
700
same problems were also recently used for the assessment of the volume-of-fluid and the level-set methods, and by comparison with these results (c.f. for example the results in [10]), it can be seen that the present model performs well in terms of accuracy of the predicted position, shape and velocity of the interface.
symmetry
695
T =T¥ Figure 22: Initial configuration for the simulation of the thin ice film propagation over the cold underlying substrate.
At the initial time the ice nucleus is initialized using Eq. (5), having the shape of a quarter of a circle, with the centre placed in the lower left corner. Melting temperature Tm is set within the ice at the beginning, and the uniform temperature T∞ = 266.9 K corresponding to the initial supercooling is prescribed outside. At all outer
ACCEPTED MANUSCRIPT
k, W/(mK) ρ, kg/m3 cp , J/(kgK)
water, ∆T = 6.25 K 0.5487 999.08 4246
ice 2.22 916.2 2050
copper 398 8954 384
760
765
The computed interface position and temperature distribution around the ice film at different times are shown in Fig. 23 for the mesh with 100 cells in both directions.
CE
PT
ED
M
AN US
770
During the evolution of the ice film, heat is being transported from the ice-water region into the substrate. Hence, the temperature within the ice film does not stay at the initial value of Tm , but it decreases from Tm at the interface to a lower value at the wall surface. The temperature change normal to the interface in Fig. 23 is much less sharp far from the tip of the ice film at the wall. Around the tip, the thermal diffusion layer in the normal direction to the interface remains much thinner, yielding much higher temperature gradients in the normal direction to the interface close to the wall. This is expected, since the underlying substrate has much higher thermal diffusivity and absorbs effectively the heat released during the solidification of water. Therefore, the horizontal propagation velocity of the ice film tip at the wall is much higher compared to the velocity of the interface in the vertical direction, leading to the obtained elongated shape of the developed ice film. The computed velocity distribution around the tip of the ice film at t = 0.5 ms is shown in Fig. 24, for the same mesh. In accordance with the above explanation, the velocities in Fig. 24 within the band [φmin , φmax ] are highest in the computational cells closest to the wall surface.
CR IP T
Table 7: Thermophysical properties for water, ice and substrate.
AC
750
boundaries the temperature is kept constant and equal to T∞ , except at the left boundary, which is the plane of symmetry. Thermophysical properties for the supercooled755 water, ice and substrate made of copper used in the calculations are listed in Table 7, and the latent heat of melting is L = 333000 J/kg.
Figure 24: The computed velocity distribution around the tip of the ice film at t = 0.5 ms, for the mesh with 100×100 cells, the supercooling ∆T = 6.25 K and copper substrate. 775
780
Figure 23: The computed interface position and temperature distribution around the ice film at different times for the mesh with 100×100 cells, the supercooling ∆T = 6.25 K and copper substrate. The time sequence is from top to bottom, left to right.
785
18
The computed temporal evolution of the interface is shown in Fig. 25, in which the evolution of the iso-contour point moving along the wall, and of the point moving along the symmetry plane is plotted for different mesh resolutions. The magnitude of the computed corresponding propagation velocity for the same two points is plotted in Fig. 26. The computed distribution of the temperature along the wall and the symmetry plane is shown in Fig. 27 at different time instants for different mesh resolutions. From Figs. 25–27, several conclusions can be drawn. The velocity at a time t is calculated from Eq. (44) and the computed interface propagation velocity at the wall again exhibits oscillations, but the amplitude is much higher in this
ACCEPTED MANUSCRIPT
20
0,05
0 0
0,02
3 0,3 0,4 t, ms
0,5
0,6
0,7
vn , cm/s
0,2
0,2
AN US
0,1
0,007
0,3 0,4 t, ms
0,5
0,6
0,7
100×100 200×200 300×300
2
1
0,006 y, mm
0,1
4
100×100 200×200 300×300
0,008
0,005
M
0,004
100×100 200×200 300×300 0,1
0,2
ED
0,003 0,002 0
CR IP T
x, mm
5
0,03
0,01
0,3 0,4 t, ms
0,5
0,6
0,7
CE
PT
Figure 25: The computed temporal evolution of the iso-surface points at the wall (top) and at the symmetry plane (bottom) for different800 mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate.
AC
However, the mean velocity, which was determined in experiments to be constant at a specific supercooling and a specific substrate [3], also stays constant in the simula-805 tion (except at the beginning, before the temperature field develops). Therefore, the magnitude of the mean velocity at a time t, plotted with solid lines in Fig. 26, is calculated in the simulation from the expression vn,t
795
15 10
0,04
0 0
exp. at T = 266.9 K 100×100 200×200 300×300
25
vn , cm/s
790
30
case, as seen in the dashed-line plot in the upper diagram in Fig. 26. As previously explained, the oscillations are attributed to the time instants when the interface passes the cell-centres. Similar behavior is detected for the velocity at the symmetry plane, but with much lower amplitudes.
810
t
1 X = Pt (xp,t − xp,t−∆t ). 0 ∆t 0
(51)
Another observation is the difference in the computed values for the interface velocity at the wall and at the815 symmetry plane. While the computed values at the symmetry plane converge with increasing mesh resolution, the values at the wall are not mesh independent. With increasing mesh resolution, the interface propagation veloc19
0 0
0,1
0,2
0,3 0,4 t, ms
0,5
0,6
0,7
Figure 26: The magnitude of the computed interface propagation velocity at the wall (top) and at the symmetry plane (bottom) for different mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate. Solid lines are plots of the mean velocity, and the dashed line is the instantaneous velocity.
ity at the wall also increases and does not converge towards a finite value (the experimentally determined value about 15 cm/s for this case [3]). This is further illustrated in Fig. 28, which shows that the interface velocity at the wall linearly increases with a linear increase in the mesh resolution. Similar behavior is detected also in simulations with a higher supercooling of ∆T = 12.35 K, the results of which are plotted in Fig. 29. On the other hand, the mesh independence of the solution at the symmetry plane is confirmed in all results for the interface position, its propagation velocity and the temperature distribution along the symmetry plane. As stated above, the mean velocities, shown in Fig. 26 and Fig. 29, are constant in time, and the small increase at later times is solely because the interface approaches the boundary of the computational domain, where the temperature is kept fixed. Inspection of the results in Fig. 27 reveals that the temperature distribution at the wall has a similar shape for all mesh resolutions at different times. With increasing the mesh resolution, the computed temperatures in cells clos-
ACCEPTED MANUSCRIPT
30
276
25
T, K
vn , cm/s
272 270
0,02
0,01
274
x, mm
0,03
0,04
10 0
0,05
red: 100×100 blue: 200×200 green: 300×300
272
825
268
0,02
0,03
0,04
0,05 830
M
y, mm
PT
exp. at T = 266.9 K comp.
840
CE
20 15
AC
vn , cm/s
25
10
ED
Figure 27: The computed temperature distribution along the wall (top) and the symmetry plane (bottom) at different time instants for different mesh resolutions, the supercooling ∆T = 6.25 K and copper substrate. 835
30
100
200 300 number of cells at the wall
845
400 850
Figure 28: The increase in magnitude of the computed interface propagation velocity at the wall with increasing mesh resolution, for the supercooling ∆T = 6.25 K and copper substrate. 855
820
0,4 t, ms
0,6
0,8
expressions (28)–(29). At first sight, the problem could be overcome by applying a constant distance δpp for each mesh. However, in that case the thermal diffusion layer could be skipped for a fine enough mesh, and a further mesh refinement would then even lead to a decrease in the calculated temperature gradient. The mesh dependence of the results at the wall is not attributed to the fixed-value temperature condition at the boundaries. The computations with the same initial conditions were performed using adiabatic boundaries instead of T = T∞ in Fig. 22, except at the bottom of the substrate. The results (not shown here) were similar to the ones obtained with the fixed temperature at the boundaries. This is explained as the consequence that the thermal diffusion layer around the interface is very thin and does not reach the right boundary for the most of the simulation time. The reason for the mesh dependence of the solution at the wall is also not the small dimensions of the computational domain. The thermal diffusion layer in the water for the specific conditions is estimated as δt ∼ αl /vn = 8.3 · 10−7 m, and even with the small computational domain, having dimensions of 50 µm × 50 µm, the solution should not be influenced by the conditions at the boundaries, at least at the beginning of the simulation. This is confirmed by the results at the symmetry plane, that do not show such mesh dependence. It is further demonstrated in Fig. 30, which shows the results for the computed interface velocity, obtained in simulations with a computational domain that has 20 times longer horizontal and vertical dimensions. The interface propagation velocity is only slightly smaller with the much larger computational domain, indicating that the outer boundaries, toward which the interface moves, have negligible influence on the solution for the evolution of the ice-film tip at the wall. The freezing of water directly at the substrate surface is determined by the local thermal conditions around the tip, thus confirming the assumption for the substrate acting purely as an
AN US
T, K
270
0,01
0,2
Figure 29: The magnitude of the computed interface propagation velocity at the wall for different mesh resolutions, the supercooling ∆T = 12.35 K and copper substrate.
solid line: t = 0.1 ms dashed line: t = 0.2 ms
266 0
20
15
268 266 0
100×100 200×200 300×300
solid line: t = 0.1 ms dashed line: t = 0.2 ms
CR IP T
274
red: 100×100 blue: 200×200 green: 300×300
est to the interface do not change much, but the calculated temperature gradients increase with mesh refinement due do smaller distances to the probe points, according to the 20
ACCEPTED MANUSCRIPT
895
25 exp. at T = 266.9 K -5 domain (5×5).10 m
vn , cm/s
20
-3
domain (1×1).10 m 900
15 10
905
5 0 0
0,1
0,2
t, ms
0,3
0,4
0,5 910
AN US
Figure 30: The comparison of the magnitude of the computed interface propagation velocity at the wall for the supercooling ∆T = 6.25 K and copper substrate, with that computed using a computational domain with 20 times longer dimensions. 860
880
885
890
M
ED
PT
875
CE
870
915
The conclusion is that the described model exhibits a good capability to compute the two-phase solidification far from the wall. However, the solidification at the wall, in particular the contact line (triple point water-ice-substrate in a two-dimensional case), imposes a singularity, which is920 indicated by the continuously increasing interface velocity when the triple point is approached by refining the mesh. Similar behavior was detected also in the aforementioned previous related studies, where the interface was even not perpendicular to the wall surface. In a recent925 study [39], where numerical simulations of a drop solidifying on a cold substrate were conducted, a three-phase point which might appear similar to the present study was encountered. However, it is related to a different underlying physics, since the studied cases include slow solidification930 with several orders of magnitude larger time scales, and the triple point, which connects the same liquid and solid material, moves away from the substrate and not along it. Thus, the triple point is of completely different nature compared to the problem in the present study. Obviously,935 in order to model this phenomenon for three materials of different kind, as a prerequisite a better understanding is required of the physics in the local area around the icing front at the substrate.
AC
865
balance, accounted for in the discretization by using reconstructed interface points. The model has been tested by computing two-phase solidification problems with existing analytical solutions. An extension of the computational model to incorporate conjugate heat transfer has been presented that can be used as a framework for combining additional models to account for the interaction of the freezing of the supercooled water with the cold, solid substrate. It has been demonstrated that the model has a good capability to compute the solidification in the absence of walls, or far from the wall. At the wall, the problem formulation imposes a singularity, indicated by diverging velocities at the contact line, which is formed during solidification of liquid over a solid substrate. The heat flux singularity at the contact line presents a fundamental difficulty for modeling the freezing within the finite-volume continuum mechanics. In accordance to the previous studies, this intrinsic singularity is also reproduced with the present model. New understanding of the underlying physics must be acquired, in order to formulate plausible models which would remove the singularity and provide means for using them within the conventional continuum framework. Experimental data at extremely small length and time scales around the contact line would be a very useful basis, though hardly accessible in measurements, in particular in the case of a contact line moving at relatively high speeds. If such new models would require including the capillary undercooling in regions of very large curvatures of the interface, the present model can be easily extended by calculating the curvature-dependent temperature and replacing the melting temperature in the points at the interface. The model can be used for computing the solidification of supercooled water, in which the phase-change process can be isolated and analyzed separately without the bulk liquid motion. As opposed to the standard interfacetracking algorithms, such as volume-of-fluid, where the Navier-Stokes flow solver is coupled with an algorithm for the proper advection of the tracked interface by the underlying velocity field, the interface in the present model propagates with the normal velocity obtained from the interfacial energy balance calculated only within the diffused interface band. Thus, the interface propagation is determined here using an algorithm based on a different underlying physics compared to the standard volume-of-fluid interface tracking. In terms of the computational costs, the model will by definition perform faster than the level-set method, since the re-initialization of the phase-field variable and the computation of the extended velocity field, required for the re-initialization in the entire computational domain, is not necessary in the present model.
CR IP T
effective heat absorber.
940
5. Conclusions
A computational model for the freezing of supercooled water within the standard cell-centered finite-volume CFD framework has been presented. The model uses the phasefield approach for tracking the interface between the solid945 Acknowledgement and the liquid phase. Energy equations in both phases are decoupled at the interface by imposing the melting temWe gratefully acknowledge financial support from the perature at the interface as a boundary condition. The Deutsche Forschungsgemeinschaft within the collaborative velocity of the interface is determined by the local energy research project SFB-TRR 75 (TP-C3). 21
ACCEPTED MANUSCRIPT
980
985
990
995
1000
1005
1010
1015
[23]
[24] [25]
[26]
[27] [28]
[29]
[30] [31]
equiaxed dendritic solidification of a binary alloy, Computational Materials Science 112 (2016) 304–317. G. Ivantsov, Temperature field around spherical, cylindrical, and needle-shaped crystals which grow in supercooled melt, Doklady Akademii Nauk SSSR 58 (1947) 567–569. B. Cantor, A. Vogel, Dendritic solidification and fluid flow, Journal of Crystal Growth 41 (1977) 109–123. F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, An Advanced InR and Matlab , R Springer Intertroduction with OpenFOAM national Publishing Switzerland, 2016. V. Alexiades, A. Solomon, Mathematical modeling of melting and freezing processes, Hemisphere publishing corporation, Washington, 1993. F. Martinez, Beyond the classical stefan problem, Ph.D. thesis, Universitat Polit` ecnica de Catalunya (2014). Y. Sun, C. Beckermann, A two-phase diffuse-interface model for hele-shaw flows with large property contrasts, Physica D 237 (2008) 3089–3098. Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, Journal of Computational Physics 220 (2007) 626–653. R The OpenFOAM extend project, http://www.extendproject.de/. R OpenFOAM , The open source CFD toolbox, http://www.openfoam.com/. ˇ Tukovi´ P. Cardiff, Z. c, H. Jasak, A. Ivankovi´ c, A block-coupled finite volume methodology for linear elasticity and unstructured meshes, Computers and Structures 175 (2016) 100–122. A. Criscione, Influence of ice formation on drop dynamics, Ph.D. thesis, Technische Universit¨ at Darmstadt (2014). A. Shibkov, Y. Golovin, M. Zheltov, A. Korolev, A. Leonov, Morphology diagram of nonequilibrium patterns of ice crystals growing in supercooled water, Physica A 319 (2003) 65–79. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– I. elements of a stability analysis, Acta Metallurgica 26 (1978) 1681–1687. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– II. instabilities in the limit of vanishing surface tension, Acta Metallurgica 26 (1978) 1689–1695. J. Langer, R. M¨ uller-Krumbhaar, Theory of dendritic growth– III. effects of surface tension, Acta Metallurgica 26 (1978) 1697– 1708. S. Patankar, Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, Washington, DC, 1980. T. Vu, G. Tryggvason, S. Homma, J. Wells, Numerical investigations of drop solidification on a cold plate in the presence of volume change, International Journal of Multiphase Flow 76 (2015) 73–85.
AN US
975
M
970
ED
965
PT
960
CE
955
[1] E. Moore, V. Molinero, Structural transformation in supercooled water controls the crystallization rate of ice, Nature 479 (2011) 506–508. [2] F. Tavakoli, S. Davis, H. Kavehpour, Freezing of supercooled1025 water drops on cold solid substrates: initiation and mechanism, Journal of Coatings Technology and Research 12 (2015) 869– 875. [3] M. Schremb, C. Tropea, Solidification of supercooled water in the vicinity of a solid wall, Physical Review E 94 (2016) 052804.1030 [4] M. Schremb, J. Campbell, H. Christenson, C. Tropea, Ice layer spreading along a solid substrate during solidification of supercooled water: Experiments and modeling, Langmuir 33 (2017) 4870–4877. [5] S. Jung, M. Tiwari, N. Vuong Doan, D. Poulikakos, Mechanism1035 of supercooled droplet freezing on surfaces, Nature communications 3 (2012) 615. [6] W. Kong, H. Liu, A theory on the icing evolution of supercooled water neat solid substrate, International Journal of Heat and Mass transfer 91 (2015) 1217–1236. 1040 [7] M. Schremb, S. Borchert, E. Berberovi´ c, S. Jakirli´ c, I. Roisman, C. Tropea, Computational modelling of flow and conjugate heat transfer of a drop impacting onto a cold wall, International Journal of Heat and Mass Transfer 109 (2017) 971–980. [8] J. Blake, D. Thompson, T. Strobl, Simulating the freezing1045 of supercooled water droplets impacting a cooled substrate, The American Institute of Aeronautics and Astronautics AIAA Journal 53 (2015) 1725–1739. [9] P. Rauschenberger, B. Weigand, A volume-of-fluid method with interface reconstruction for ice growth in supercooled water,1050 Journal of Computational Physics 282 (2015) 98–112. [10] P. Rauschenberger, A. Criscione, K. Eisenschmidt, D. Kintea, ˇ Tukovi´ S. Jakirli´ c, Z. c, I. Roisman, B. Weigand, C. Tropea, Comparative assessment of volume-of-fluid and level-set methods by relevance to dendritic ice growth in supercooled water,1055 Computers & Fluids 79 (2013) 44–52. [11] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, Cambridge, UK, 2011. [12] S. Deshpande, L. Anumolu, M. Trujillo, Evaluating the perfor-1060 mance of the two-phase flow solver interfoam, Computational Science & Discovery 5 (2012) 014016. [13] T. Waclawczyk, T. Koronowicz, Comparison of CICSAM and HRIC high-resolution schemes for interface capturing, Journal of Theoretical and Applied Mechanics 46 (2008) 325–345. 1065 [14] A. Bourdillon, P. Verdin, C. Thompson, Numerical simulations of water freezing processes in cavities and cylindrical enclosures, Applied Thermal Engineering 75 (2015) 839–855. [15] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, International Journal of Heat and Mass Transfer 30 (1987) 1709–1719. [16] S. Schiaffino, A. Sonin, Motion and arrest of a molten contact line on a cold surface: An experimental study, Physics of Fluids 9 (1997) 2217–2226. [17] S. Schiaffino, A. Sonin, On the theory for the arrest of an advancing molten contact line on a cold solid of the same material, Physics of Fluids 9 (1997) 2227–2233. [18] D. Anderson, S. Davis, Local fluid and heat flow near contact lines, Journal of Fluid Mechanics 268 (1994) 231–265. [19] M. Quintella, Numerical investigation of the heat flux singularity at an advancing molten contact line, Master’s thesis, Massachusetts Institute of Technology (1999). [20] I. Steinbach, C. Beckermann, B. Kauerauf, Q. Li, J. Guo, Threediensional modeling of equiaxed dendritic growth on a mesoscopis scale, Acta materialia 47 (1999) 971–982. [21] I. Steinbach, H. Diepers, C. Beckermann, Transient growth and interaction of equiaxed dendrites, Journal of Crystal Growth 275 (2005) 624–638. [22] Y. Souhar, V. De Felice, C. Beckermann, H. Combeau, M. Zaloˇ znik, Three-dimensional mesoscopic modeling of
AC
950
1020
CR IP T
References
22
[32]
[33] [34]
[35]
[36]
[37]
[38] [39]