- Email: [email protected]
A numerical method for predicting crystal micro-structure
A numerical method for predicting crystal micro-structure Stuart J. Galloway and Noel F. Smyth Department of Mathematics and Statistics, University of Edinburgh, Edinburgh, Scotland, UK
The numerical solution of the Stefan problem in a two-dimensional region is considered using a boundary integral technique. Heat jlow both ahead of and behind the front is included. This numerical formulation is found to have a number of aduantages ouer finite difference and finite element techniques. For example, complex boundaries can be considered and the discretisation of the Stefan condition is not needed. The numerical method is also extended from that ofpreuious work to enable the front to be found when the ratio of the thermal diffsivities of the liquid phase to the solid phase is large. Using a simple model of crystal growth in a freezing material the numerical scheme will be used to predict crystal size in a frozen material. Predictions using this model will be compared with experimental results. 0 1997 by Elsevier Science Inc.
1. Introduction Situations in which a material undergoes a change of phase are widespread in science, engineering, and industry. Industrial applications range from metals processing to food processing to blood plasma storage. While the equations governing the phase change of a material, consisting of the heat equation in the two phases plus an energy conservation condition at the phase change front, the so-called Stefan condition, are well understood, outside of one-dimensional situations there are very few exact solutions (see for example Hill and Dewynne’). This is because the position of the phase change front is unknown and is coupled with the temperature fields in the two phases. Therefore the details of a phase change in a two- or higher dimensional material can, in general, only be calculated using numerical methods. As the motion of the phase change front must be determined along with the temperature fields in the two phases, most straightforward numerical methods suffer from unphysical oscillations in the temperature around the front. These oscillations are due to the discontinuous derivative in the temperature at the front (due to the Stefan condition). Successful numerical methods to solve phase change problems must therefore include a way to accommodate the discontinuous temperature derivative at the phase change front. A number of these numerical methods are
Address reprint requests to Dr. N. Smyth, University of Edinburgh, Department of Mathematics and Statistics, Mayfield Road, The King’s Buildings, Edinburgh EH9 352, Scotland, UK. Received 1997
30 September
1996; revised
2 May 1997; accepted
Appl. Math. Modelling 1997, 21: 569-578, September 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
23 May
given in the review articles by Fox’ and Voller et al.3 One successful method is the enthalpy method in which the discontinuous temperature derivative is by-passed by solving for the enthalpy instead of the temperature. The enthalpy method and its numerical implementation is discussed by Elliot and Ockendon,“ in the review article on finite element and finite volume methods for Stefan problems by Voller et al.3 and in the book on the numerical solution of moving boundary problems by Finlayson.’ Finite difference methods based on the enthalpy method work well for rectangular regions’ but are difficult to apply when the domain is not rectangular. Finite element methods overcome this drawback.3,5,7,s As an alternative to finite element methods, Coleman’ used a boundary integral method to solve phase change problems in which the temperature ahead of the phase change front is kept constant at the fusion temperature. This method has a great advantage in that it can be applied easily to regions of any shape and, furthermore, the Stefan condition can be fitted naturally into the numerical scheme without discretisation. The boundary integral method reduces the problem of solving for the temperature field in the two phases to solving for unknowns on the boundary of the material and on the phase change front only. Thus a two-dimensional problem is reduced to a one-dimensional problem on two curves, with a resulting great saving in calculation. Coleman’ applied his boundary integral method to solving phase change problems on squares and L-shaped regions, and it was found to give fast, accurate results in good agreement with other methods. The boundary integral method does, however, have a drawback in that it can only be applied to well-posed Stefan problems and not to supercooled situations, unlike the enthalpy method, which can be applied in both cases.3
0307-904x/97/$17.00 PII s0307-904x(97)00054-1
Crystal microstructure
prediction: S. J. Galloway
and N. F. Smyth
In the present work the boundary integral method of Coleman’ will be extended in two ways. The boundary integral formulation will be extended to include heat flow ahead of the phase change front, so that more general and realistic initial conditions can be used. Furthermore, Coleman’ used Newton’s method to find the position of the phase change front. It is found that for the case of the freezing of a liquid, when the ratio of the thermal diffusivities of the liquid and the solid phases is large, Newton’s method has difficulty converging due to the resulting small temperature derivatives. An example of this is the familiar case of water freezing to ice, as the ratio of the thermal diffusivities of water and ice is 8.9649. The numerical method presented in the present work overcomes this by using a bisection method to determine the position of the phase change front when Newton’s method has difficulty in converging. By symmetry, Newton’s method will also have difficulty in converging for a melting problem for which the ratio of the thermal diffusivities of the liquid and the solid phases is small. An application of the numerical method is made to predict the size of crystals formed when a material freezes, based on a simple model of crystal growth.“, l1 The predictions of this model are compared with experimental results.
2. Governing
equations
To be specific, let us consider the freezing of a liquid contained in an arbitrary two-dimensional region R. The fixed boundary of the region R will be denoted by Bf. Inside the region R there will be a moving phase change front, which will be denoted by BP. Furthermore let us denote the temperature of the liquid by T,(r, t) and the temperature of the solid by T,(r, t). Let the thermal diffusivity and thermal conductivity of the liquid and the solid be v, and K, and us and KS, respectively. The equation governing the heat flow in the solid is then
JT, - = vsV2Ts
(1)
at
for points r lying between the fixed boundary Bf and the phase change front B,,. The equation governing the heat flow in the liquid is
JT,
-
dt
= vlV2T,
(2)
for points r lying ahead of the phase change front BP.’ Initially the region R is totally filled with liquid with some temperature distribution, so that the initial condition is Tl(r, 0) = T,,(r) in R.
(3)
The boundary condition on the fixed boundary general, a Newton cooling condition n.VT, = p(T, - T,) on Bf,
570
Appl.
Math.
Modelling,
Bf is, in
(4)
1997,
Vol. 21, September
where n is the outward normal to the boundary and T, is the temperature of the surroundings. On the moving boundary BP the Stefan condition K,z-Ks-$=
-Lpn*V
(5)
holds, where n is the outward normal to the solid region, V is the velocity of the phase change front, p is the density of the solid, and L is the latent heat of fusion.’ The Stefan condition expresses energy conservation across the phase change front with energy being liberated due to the latent heat release as the liquid freezes. In deriving the Stefan condition it has been assumed that the density of the solid and the density of the liquid are the same for simplicity. If this assumption is not made, then the bulk motion of the solid and/or liquid must be included, which increases the complexity of the governing equations greatly. The governing equations (1) and (2), initial condition (31, and boundary conditions (4) and (5) can be nondimensionalised with respect to the parameter values for the liquid phase. Let temperature be measured with respect to the fusion temperature Tf of the liquid and scaled with respect to some typical temperature T, of the initial temperature distribution. A nondimensional temperature T’ is then (6)
If a is a typical length of the region R, then nondimensional space and time variables can be defined by r’
=
f
a
and
t’=-.
v
Using these nondimensional variables gives the nondimensional equations
JT, =
-
at
(7)
a2
and dropping primes
vV2Ts
for the solid phase between the fixed boundary phase change front BP and
JT, = V2T,
-
at
(8) Bf and the
(9)
for the liquid phase ahead of the phase change front BP. The boundary condition (4) becomes n-VT, = aPT,(T, - T,)
(10)
on the fixed boundary Bf, where T, is the nondimensional temperature of the surroundings. The Stefan condition (5) becomes
aT, --_K&& an
(11)
Crystal microstructure prediction: S. J. Galloway and N. F. Sm yth on the moving phase change front BP. Here the nondimensional diffusivity v and conductivity K are
v=%
and
GJr, r’, t, t’)
K K=&.
(12)
4
VI
Let G, be the Green’s function
= H(t - t’) The Stefan number
’ ev( - ,‘:,‘i:,) 47rv(t - t’)
a is
(16)
PV,L
(131
a=KIT,* The nondimensional equations with the Stefan condition (ll), the constant > 0, and the boundary constant < 0 on the fixed boundary formally to p = 03 in [lo]) have an dimension. This solution is
(8) and (9) together .. initial condition T,, = condition T, = Tb = Bf (which corresponds exact solution in one
so that
z
+ vV’*G,=
-S(t-t’)s(r-r’)
(171
where V’ refers to derivatives with respect to r’. Then it can be shown from equations (8) and (17) that
-$(G,cl
+ vV’.(T,‘V’G,
- G,V’T,‘)
= -T,‘i?(t-t’)6(r-r’) (14)
TI = T,, -
T,O
lc)\ erfc
(see Hill and Dewynne’). The phase change cated at x = e\ll, where 13 is the solution of
T,o
e-eZ/4 +
KTb
e-
front
where T,’ = T,(r’, t’). The fixed boundary Bf and the phase change front BP can be extended as surfaces Sr and S, in space-time, respectively, as shown in Figure I. The phase change front is then described by an equation of the form t-f(r)
is lo-
= 0
in space-time. given by
eZ/(4v)
(18)
(191
The outward unit normal to S, is therefore
(Vf, -1) = -
$k
This solution will be needed to start the numerical tion to be described in the next section.
3. Numerical
np=
(20)
d_’
(15) solu-
scheme
In this section a boundary integral scheme to numerically solve the system (8) and (9) together with the Stefan condition (11) will be described. This scheme is an extension of the method of Coleman’ to include a variable temperature ahead of the phase change front and the possibility that Newton’s method, which is used to determine the position of the phase change front, will not converge. Using the boundary integral method the system of equations (8) and (9) is converted into a pair of coupled integral equations for the temperature ahead of and behind the front. Let us consider the integral equation for the temperature behind the front first.
lime
t?L Y
X
Fixed
Figure
1.
tended
as surfaces
Sf
X-scaiOO
sf boundary
B,
S, and S,
Appl. Math. Modelling,
1997,
and
phase
change
front
B,
ex-
in space-time.
Vol. 21, September
571
Crystal microstructure
prediction:
S. J. Galloway and IV. F. Smyth
Integrating equation (18) over space-time between the surfaces S, and S,, and using the divergence theorem in space-time then gives -hT,(r,f)
vi:“‘/,;icz -G~z) ds’df’
=
(21)
The parameter A takes the value 1 at points away from the boundary and the value l/2 at regular boundary points. Projecting the integral over the surface Sp down onto the x -y plane as the surface B, finally gives -,+~,(r,t)
=
vl+‘Lf(T;z-Gsz)ds’df’ - u
/
GJr,r’,t,f(x’,y’)) B,!
xV’T;V’f
]If_f(x’,y+‘dy’.
(22)
This integral equation alone is not enough to determine the temperature behind the phase change front as the gradients VT,‘. V’f ItI=f(x’,y’) on the phase change front are not known. These unknown gradients are related to the same gradients ahead of the front by the Stefan condition (11). Hence the temperature field ahead of the front also needs to be considered. The equation governing the temperature field ahead of the phase change front is (9), and so the Green’s function appropriate to the region ahead of the front is
on noting that the outward normal to the surface Sp is now -np. The surface S, is the region inside the frxed boundary Bf. The integral equations (22) behind the phase change front and (26) ahead of the front together with the Stefan condition (11) completely determine the temperature field within the material. The integral equation (22) involves both the temperature T, and the temperature gradient ilTs/an at the fixed boundary Bs. Both of these will not be known and if, say, a fixed temperature boundary condition is imposed on Bf, then the gradient aT,,/an on Bf will have to be determmed as part of the solution. Hence this gradient, the gradient V’T/*V’flt~=f(X~,Y~)on one side of the phase change front (with this gradient on the other side of the phase change front determined via the Stefan condition [ll]) and the position of the phase change front itself have to be determined as part of the numerical solution. Once these are known the temperature at any point in the region can be determined using equations (22) and (26). The numerical method used to determine these unknowns will now be described. Let time be discretised into intervals of length At, so that ti = iht, and let the temperature at time ti be ryi(r) behind the front and Tli(r) ahead of the front. Then approximating the time integrals in (22) using the midpoint rule we have that the temperature on the boundary at time step k is
- ;
Tsk(=) k-l =
ZJA~c i=
?is4(r,r’,(k
/
(i - 1/2)At) =H(t
-t’)
4T(:_t’) exp[ -
- G,(r,r’,
:(-:;;j(23)
I
- 1/2)At,
1 &f
+
(k - 1,‘2)At,
(i - 1/2)At),,,
d Tii
G, then satisfies
$ +V”G,
+ v = - i?(t - t’)t?(r - r’)
so that from equation ;(G,T;)
(24) x GJr, r’, (k - 1/2)At,
(9)
+ V’.(T;V’G/
t’) dt’
(k-1'2)ArGs(r,r',
- G,V’T;)
= - T/6(t - t’)iS(r - r’).
(k - I/Z)
(25)
At, t’) dt' In a similar manner the front, integrating the surface S, gives -AT,(r,t)
- v
= $7
+
r’,O)G,(r,
r’, t,O) dw’ dy’
xV’T; V’f\r+~,~~)
572
Appl.
Math.
Modelling,
1997,
As in Coleman,”
h’ W
/
B,
GJr,r’,(k
- 1/2)At,f(x’,y’))
x V’T;a V’f[t~=f(r’,Y’)du’ dy’.
IB, Gr(r,r’,t,f(x’,Y’)) (26)
Vol. 21, September
ds' i
as for the temperature field behind this equation in space-time ahead of
(27)
the integral over the last time step has not been evaluated using the midpoint rule due to the singular nature of the Green’s function there. Instead,
Crystal microstructure over the last time step, the approximations
(k - 1/2)A1
/ (k-
r,r’, (k - 1/2)At,
ZZ-
27TV
T,‘(r’,O)G,( r,r’,(k-
t’) dt’
(r - r’1.n (r -
e-~r~r”2/(2~A~)
/
G)(r,r’, (k - 1/21At,f(x’,
y’))
4
~V’T~.V’fI,,=~(~,,~,)dx’dy’
yr12
and 1’2)AtGs(r,r’, (k - l/2)&, / (kUC- 1)At
t’) dt’
are used, where E, is the exponential integral.‘* Now that time has been discretised, space is also discretised. To evaluate the spatial integrals in equation (27) the boundaries Bf and BP are approximated by N-sided polygons and T, and ~?T,/an are approximated by their midpoint values on the sides. The integrals are then evaluated using the midpoint rule, except for when G, and aG,Jan’ become singular. This occurs when r lies on the midpoint of a side of Bf, in which case by expanding about this point it can be shown that the term aTs/dn’ in (27) has coefficient AL(y - 2 + log(AL2/(8vAt)))/(47rv), where y is the Euler constant and AL is the length of the side.” At these singular points the term in equation (27) involving T, on Bf is zero.’ The integral over the surface B, in equation (27) is evaluated in the following manner. The phase change front E$ is approximated by a polygon at each time step. By joinmg the vertices of these polygons for times 0 up to kAt together a set of quadrilaterals is formed that cover the surface B,. The integral over each of these polygons is then approximated by the value of the integrand at the centroid times the area of the quadrilateral, which is the midpoint rule in two dimensions. This approximation is valid for all quadrilaterals except those between the (k 11th and kth time step, due to G, becoming singular on the phase change front BP. As in Coleman,’ the integral over these quadrilaterals is approximated by
(area of quadrilateral)
1/2)At,O)dx’dy’
su +
1
and IV. F. Sm yth
(26) gives
0 = -/
a dn’GS(
1)Ar
equation
prediction: S. J. Galloway
x
V’T,‘.v’.f
/r’=(k-
1/2)Ar
4rvAt
(29)
where r’ is the quadrilateral centroid. The function f is approximated by (i - 1/2)At for any quadrilateral between the (i - 1)th and ith time step. To evaluate the unknown gradients V’T,‘-V’fIr,=fcx,,y,, the integral equation (26) is applied at the phase change front BP where T, = 0. Setting T, = 0 at time t = kAt,
(301
where r lies on the phase change front BP. The second integral in this expression is evaluated in the same manner as the similar integral in equation (27). The first integral over the initial condition is evaluated by dividing the region inside Bf up into a grid. The integral is then approximated by the midpoint rule as the sum of the area of each grid box times the value of T,‘(r’,O)G, at the centroid. The numerical procedure to solve equations (27) and (30) then works in the following manner. For definiteness, let us suppose that the boundary condition on Bf is the fixed temperature boundary condition T, = Tsb < 0. We then assume that we know the temperature gradients V’T,‘.V’f(t,=f~x~,Y,) on the phase change front for times 0 up to (k - 1)At. The position of the phase change front at time kAt is then approximated by extrapolating from its previous positions. The linear system arising from the integral equation (30) can then be solved for the gradients V’T,’. V’f]r,=,-(x,,y,) on the front at time kAt. From these gradients and the Stefan condition (11) the gradients V’T,‘v’fIr,=f(x’,Y’) behind the phase change front can then be determined. Substituting these gradients behind the phase change front into the linear system resulting from the integral equation (27) the unknown normal derivatives aT,/an on the boundary Bf can then be found. Once these normal temperature derivatives are known, Newton’s method with the temperature TY given by the discretised form of equation (22) can then be used to find where the temperature T, is zero. This then gives a revised position of the phase change front at time kbt. The whole process is repeated until the difference between the positions of successive estimates of the front is within a prescribed tolerance. To start the process off the one-dimensional solution (14) and (15) is used to give the front position and gradients VT;. V’f If8=P(~‘,~‘) at time At. This will give a good approximation to these quantities at time At if At is small. To keep the phase change front smooth, when the iterations have converged, the points on the front are passed through a three-point smoothing scheme. It was found that the position of the phase change front could be successfully determined using Newton’s method to search for when T, = 0 unless v was large. When v is large, the derivative of the temperature T, is small (see the onedimensional solution [14]), and this causes Newton’s method to have difficulty in converging. To enable solutions to be found when v is large a bisection method was used to determine the position of the phase change front (searching for when T, = 0) when Newton’s method had difficulty in converging. In applying the numerical scheme outlined in this section, Newton’s method was used to determine the position of the phase change front, unless it did not converge within a fixed number of iterations. In
Appl. Math.
Modelling,
1997,
Vol. 21, September
573
Crystal microstructure
prediction: S. J. Galloway and N. F. Smyth
this case the front finding method was then switched over to the bisection method.
the phase change front will start to move instantaneously. As a check on the numerical scheme the boundary Bf was taken to be a square of side length 2, and To was set equal to zero. This then reduces to one of the examples considered by Coleman’ as the temperature ahead of the phase change front is constant at Tl = 0. With Tb = - 1, Q = 1, or= 1, and K = 1 the time for the complete freezing of the region was found to be 0.474695, which agrees with the value 0.475 found by Coleman.’ Figure 2 shows the solution for Tb = - 1.0 and T,, = 0.65 with the parameter values v = 8.9649, CY= 3.5732, and
4. Results The boundary integral method outlined in the previous section was applied to a number of examples of the freezing of a region. For simplicity the boundary condition applied on Bf was taken to be T, = Tb < 0 and the initial condition was taken to be T, = T,, a 0, where Tb and T,, are constants. With these initial and boundary conditions
1 0.5 T
0 -0.5 -1
0
0
x
.I
-
.6
-
..
-
.2 -
.
I-
-.2
-
-.1
-.6
-.a
-
Cc) Figure 2. (a) Contour plot of the temperature for (x, y) in the fourth quadrant of a square of side length 2 at time ~0.0875 for Tb= - 1 .O, Tc=O.65, v=8.9649, o=3.5732, and K=4.1516 with T,=2O”C. (b) Surface plot of the temperature for Ix, y) in the fourth quadrant of a square of side length 2 at time t=0.0875 for the parameters of (a). (c) Phase change front positions for a square of side length 2 at time intervals of 0.0175 up until complete freezing for the same parameter values as (a).
574
Appl.
Math.
Modelling,
1997,
Vol. 21, September
Crystal microstructure K = 4.1516. These parameter values are those appropriate for water freezing to ice with T, = 20”C.‘3 For these parameter values, Newton’s method had difficulty in converging, and a bisection method was used to determine the position of the phase change front. The region inside Bf was taken to be a square of side length 2. Figures 2(a) and 2(b) show contour and surface plots for the temperature field at time t = 0.0875. The phase change front can be clearly seen, particularly in the contour plot. The time evolution of the phase change front is shown in Figure 26~) at time intervals of 0.0175. It can be seen that the front is evolving to a circular shape as time increases. The freezing time for the parameter values of Figure 2 is 0.325134, which is smaller than the freezing time for the example of Coleman.’ This is due to the larger values of the thermal conductivity K and the thermal diffusivity V, which means that heat is more easily transported to the boundary Bf. The Stefan number a is larger than for the example of Coleman,’ which means that more heat is liberated on freezing. This and the larger value of the initial temperature TO would cause the freezing time to increase over the value of Coleman.’ However the larger values of the thermal diffusivity and the thermal conductivity have the greater effect, and the freezing time decreases. The fact that the freezing time is not much reduced over the value of Coleman’ shows that the increased values of V, K, (Y, and T,, have effects that nearly balance out. Figure 3(a) shows the variation of the freezing time tf with TOfor Tb = - 1.0, u = 8.9649, CY= 1.0 and K = 4.1516. These parameter values are those for water, except for (Y. For this high value of v Newton’s method has difficulty in converging, and the bisection method is used to locate the position of the phase change front. As expected the freezing time increases with the initial temperature TO, with the increase being nearly linear. The increase in the freezing time tf as the initial temperature T,, ranges from 0 to 1 is not very large because the boundary temperature is relatively low at Tb = - 1. This low boundary temperature and the large relative thermal diffusivity Y of the solid phase then dominate the heat flow. Figure 3(b) shows the variation of the freezing time tf with the initial temperature TO for Tb = - 1.0, v = 0.5, (Y= 3.5732, and K = 4.1516. The value of (Y used for this figure is the value appropriate for water. The freezing time again increases with increasing initial temperature T,, in a nearly linear fashion. The freezing times in Figure 3(b) are larger than those in Figure 3(a) due both to the smaller value of v and to the larger value of (Y. The smaller value of v means that heat diffuses out of the region at a slower rate, and the larger value of (r means that more heat is liberated on freezing, which then has to diffuse out of the region. Figure 4 shows the variation of the freezing time tf with T,, for T,, = - 1.0 and the water parameters (v = 8.9649, (Y= 3.5732, and K = 4.1516). The freezing time again increases with T,,, as expected, but in this case the increase is highly nonlinear, with the rate of increase slowing with increasing TO. The freezing times are between those of Figures 3(a) and 3(b). The freezing times are increased over those of Figure 3(a) due to the larger value of (Y
prediction:
S. J. Galloway and N. F. Smyth
I
I
0.15 _ 0.145 0.14 0.135 0.13 tt
0.125 0.12 -
0.54;
0.2
0.4
.l-
0.6
I 1
0.8
‘0
W Figure 3. (a) Plot of freezing time t, as a function of To for a square of side length 2 for the parameter values Tb= - 1.0, v=8.9649, a= 1.0, and K=4.1516 with T,=2O”C. (b) Plot of freezing time t, as a function of To for a square of side length 2 for the parameter values Tb= - 1.0, v=O.5, (~=3.5732, and K=4.1516 with T,=2o”C.
leading to more heat liberated on freezing and are decreased over those of Figure 3fb) due to the larger value of v leading to more rapid loss of heat.
5. Crystal formation The numerical method outlined in Section 3 can be used, together with a simple model of crystal formation, to predict the sizes of crystals formed when a liquid freezes. If the growth of the crystal is radially symmetric, Frank” and Coriell and McFadden” showed that the size of the crystal formed when a liquid freezes is given in dimensional units by
(31)
Appl.
Math.
Modelling,
1997,
Vol. 21, September
575
Crystal microstructure
Figure
4.
square
of side
v=8.9649,
prediction: S. J. Galloway and N. F. Smyth
Plot of freezing length
a=3.5732,
K are those for water
2 for
time the
t,
and K=4.1516 freezing
as a function
parameter
values
(the values
to ice with
of
T, T,=
for
a
- 1 .O,
of Y. a, and
T,= 20°C).
Here K, is the thermal conductivity of the liquid, L is the latent heat, T, is the freezing temperature of the liquid, and T, is the temperature far ahead of the growing crystal. It should be emphasised that this expression is valid for the stable freezing of a liquid and is not valid for the dendritic growth that occurs when the liquid is undercooled.14 The size expression (31) holds for a single crystal that has been growing for a time t and was found to give crystal sizes in accord with experimental results.” When a region freezes, as in the numerical solutions of the previous section, a large number of crystals form, the crystal at a given point being formed when the phase change front passes that point. Hence to apply the crystal growth expression (31) to the freezing of a region, this expression needs to be applied to each crystal formed (i.e., at each point of the region), where the time t is the time taken for a given crystal to form. In the present case the crystal size expression (31) is most easily applied if the time t is replaced by the front velocity. This is in agreement with models of microstructure formation that show that the size of microstructural features is primarily determined by the speed of the phase change front.‘5,‘6 If the phase change front were planar, then the nondimensional front position is given by x = 0fi (see equation [15]), so that the nondimensional front velocity is I/ = 0/(2fi). Eliminating time in favor of front velocity in the crystal size expression (31) and using equation (7) to convert to dimensional variables gives the alternative crystal size expression
(32)
where change planar phase tional
576
cI is the specific heat of the liquid. While the phase fronts considered in the present work are not it was noted by Coleman’ that even for nonplanar change fronts the position of the front is proporto fi to a good approximation. Therefore the
Appl.
Math.
Modelling,
1997,
Vol. 21, September
crystal size expression (32) will be used in the present work on the freezing of two-dimensional regions. The crystal sizes formed when the liquid in a region freezes are then calculated in the following manner using the boundary integral method of Section 3. Using the positions of the phase change front calculated by the numerical scheme of Section 3 the normal velocity V, of a point on the front at a given time step can be calculated. This normal velocity V, is then used as the velocity V’ in the crystal size expression (32). The crystal size expression (32) is now employed to calculate the crystal size in each quadrilateral used to calculate the surface integral B,. At each time step the parameter 0 is determined from the transcendental equation (151, with Tb set equal to the boundary temperature and T,, set equal to the temperature at the centre of the square of Section 4. While this transcendental equation is strictly valid only for a planar phase change front the crystal size was found to be only weakly dependent on the value of 19. Indeed the crystal sizes were found not to change to a significant extent if 0 was set equal to its initial value. The value of 0 was calculated using equation (19, so that the effect of the changing temperature of the liquid could be included. The temperature T, is taken to be the boundary temperature T,, and, for water, the freezing temperature is T, = 0, these temperatures being made dimensional using equation (6). Hence the crystal size in each quadrilateral used to calculate the surface integral B, can then be calculated. Adding up all such local crystal sizes over all positions of the phase change front until the region has completely frozen and dividing by the area of the region then gives the average crystal size $. Figure 5(u) shows a plot of the mean crystal size ?? as a function of the freezing time tf for Tb = - 1.0 and the water parameters with T, = 20°C and a = 5.0 cm. The mean crystal size is nearly linear in tf, except near tf = 0.34 (ie. except for initial temperatures near To = 1.0: see Figure 4). The rapid increase in the mean crystal size near To = 1 is related to the slow increase of the freezing time with initial temperature To near To = 1, which can be seen in Figure 4. The increase of the mean crystal size with freezing time is expected since a longer freezing time means a smaller mean front speed, which in turn leads to a larger crystal size. The actual mean crystal sizes, while a little on the high side, are in reasonable accord with the ice crystal sizes reported by Coriell and McFadden.” One reason for the mean crystal sizes being somewhat higher than experimental values is that the fluid flow ahead of the phase change front has been neglected. This flow is expected to produce faster freezing times and hence smaller mean crystal sizes. The inclusion of fluid flow is the subject of current work. The experimental measurements of Miyawaki et al. ” of the freezing of soy protein curd show a basically linear dependence of the mean crystal size on freezing time (with a large amount of experimental scatter), which is in qualitative agreement with our simple model of crystal formation on freezing. Unfortunately the exact experimental conditions were not reported in this work, so that quantitative comparisons with their results cannot be made. Figure 5(b) shows the
Crystal microstructure
prediction: S. J. Galloway and N. F. Sm yth
change front is taken into account, extending the work of Coleman,’ which considered heat flow behind the phase change front only. It is found that this boundary integral method formulation of the Stefan problem has the advantages of (i) being a fast numerical method that yields accurate solutions and (ii) that it can easily incorporate nonrectangular boundaries. When parameter values appropriate for water were used it was found that the standard method for determining the front position, based on Newton’s method, was unstable, and so a new method for determining the front position, based on the bisection method, was developed for these parameter values. Using a simple model for the growth of a crystal from a freezing liquid the numerical method was extended to predict the sizes of the crystals in the crystalline microstructure formed when a body of liquid freezes. The predictions of this model for crystal formation were found to be in qualitative agreement with experimental results. While this model of crystal size was developed in the context of water it can easily be applied to other areas in which the prediction of crystal size is important, for example the growing of silicon crystals for the computer industry.
7.6 7.4 7.2 76.6 6.6 6.4 -
6.6 r 6.5 6.4 6.3 6.2 . 6.1 .
Acknowledgments
6-
The authors the referees presentation EPSRC for dentship.
5.9 5.6 5.7 -/ 5.6 I -1
-0.9
-0.6
0.7
-0.6
-0.5
-0.4
would like to express their appreciation to for pointing out references that improved the of this paper. S.J.G. would like to thank the funding through an earmarked Ph.D. stu-
I -0.2
-0.3
5
(W
Figure 5. the
freezing
v=8.9649, 5.0
(al
Plot of the
time
t,
a=3.5732,
cm. (b) Plot of the
boundary v=8.9649.
temperature (~=3.5732,
for
mean T,=
and mean
crystal
- 1.0 K=4.1516 crystal
Tb for T,=O.l and
K=4.1516
and
size the with
3 as a function water
size 3 as a function with
1. Hill, J. M. and Dewynne, J. N. Heat Conducfion. Blackwell, Oxford, UK, 1987 2. Fox, L. The Stefan problem: Moving boundaries in parabolic equations. A Surrvy of Numerical Methods for Partial Differential Equafions, eds. 1. Gladwell and R. Wait, Clarendon Press, Oxford, UK, 1979 3. Voller, V. R., Swaminathan, C. R., and Thomas, B. G. Fixed grid techniques for phase change problems: A review. Inr. J. Num.
parameters
T,=2o”C
and the water
References
of
and
a=
of the
parameters
T,=2O”C
and
a=
5.0 cm.
Method. Engin. 1990, 30, 875-898 4. Elliot, C. M. and Ockendon, J. R. Weak
dependence of the mean crystal size 5 on the boundary temperature Tb for fixed initial temperature T,, = 0.1 and the water parameters with T, = 20°C and a = 5.0 cm. The mean crystal size decreases with decreasing boundary temperature, which is expected since a lower boundary temperature implies a larger mean front speed, which leads to a smaller crystal size. The mean crystal size dependence shown in Figure S(b) shows a good qualitative similarity to the experimental measurements for soy protein curd obtained by Ref. 17.
5. 6.
7.
8.
9.
6. Conclusions
and variational
methods
for moving boundary problems. (Pitman Advanced Publishing Program) Res. Notes Math. 1982, 59 Finlayson, B. A. Numerical Methods for Problems with Mol+zg Fronts. Ravenna Park Publishing, Seattle, WA, 1992 Voller, V. and Cross, M. Accurate solutions of moving boundary problems using the enthalpy method. Int. J. Heat Mass Transfer 1981, 24, 545-556 Lynch, D. R. Unified approach to simulation on deforming elements with application to phase change problems. J. Camp. Phys. 1982, 47, 387-411 Chow, P. and Cross, M. An enthalpy control-volume-unstructured-mesh (CV-UM) algorithm for solidification by conduction only. ht. J. Num. Method. Eng. 1992, 35, 1849-1870 Coleman, C. J. A boundary integral formulation of the Stefan problem. Appl. Math. Mode&g 1986, 10, 445-449
10. Frank, F. C. Radially symmetric phase growth controlled diffusion. Proc. Roy. Sot. Land. A 1950, 201, 586-599
It has been shown that a boundary integral method can be used to find a numerical solution of the Stefan problem when heat flow both ahead of and behind the phase
11. Coriell,
S. R. and
McFadden,
G. B. Morphological
Handbook of Ctysral tirowth I, Fundamentals, and Stability, ed. D. T. J. Hurle, North-Holland,
Appl.
Math.
Modelling,
1997,
by
Stability.
Part B: Transport
Amsterdam,
Vol. 21, September
1993
577
Crystal microstructure
prediction:
S. J. Galloway and N. F. Smyth
12. Abramowitz, M. and Stegun, I. Handbook of Mathematical Functions. Dover Publications, New York, 1972 13. Lide, D. R., ed. CRC Handbook of Chemistry and Physics, 74th edition, Chemical Rubber Publishing Company, Boca Raton, FL, 1993 14. Langer, J. S. Instabilities and pattern formation in crystal growth. Reo. Mod. Phys. 1980, 52(l), l-28 15. Rappaz, M., Carrupt, B., Zimmermann, M., and Kurz, W. Nu-
578
Appl.
Math.
Modelling,
1997,
Vol. 21, September
merical-simulation ment of materials.
of eutectic solidification in the laser He/u. Phys. Acta 1987, 60, 924-936
16. Rappaz, M. Modelling of microstructure formation tion processes. fat. Mater. Reu. 1989, 34(j), 93-123
treat-
in solidifica-
17. Miyawaki, O., Abe, T., and Yano, T. Freezing and ice structure formed in protein gels. Biosci. Biotech. Biochem. 1992, 56(6), 953-957
The numerical solution of the Stefan problem in a two-dimensional region is considered using a boundary integral technique. Heat jlow both ahead of and behind the front is included. This numerical formulation is found to have a number of aduantages ouer finite difference and finite element techniques. For example, complex boundaries can be considered and the discretisation of the Stefan condition is not needed. The numerical method is also extended from that ofpreuious work to enable the front to be found when the ratio of the thermal diffsivities of the liquid phase to the solid phase is large. Using a simple model of crystal growth in a freezing material the numerical scheme will be used to predict crystal size in a frozen material. Predictions using this model will be compared with experimental results. 0 1997 by Elsevier Science Inc.
1. Introduction Situations in which a material undergoes a change of phase are widespread in science, engineering, and industry. Industrial applications range from metals processing to food processing to blood plasma storage. While the equations governing the phase change of a material, consisting of the heat equation in the two phases plus an energy conservation condition at the phase change front, the so-called Stefan condition, are well understood, outside of one-dimensional situations there are very few exact solutions (see for example Hill and Dewynne’). This is because the position of the phase change front is unknown and is coupled with the temperature fields in the two phases. Therefore the details of a phase change in a two- or higher dimensional material can, in general, only be calculated using numerical methods. As the motion of the phase change front must be determined along with the temperature fields in the two phases, most straightforward numerical methods suffer from unphysical oscillations in the temperature around the front. These oscillations are due to the discontinuous derivative in the temperature at the front (due to the Stefan condition). Successful numerical methods to solve phase change problems must therefore include a way to accommodate the discontinuous temperature derivative at the phase change front. A number of these numerical methods are
Address reprint requests to Dr. N. Smyth, University of Edinburgh, Department of Mathematics and Statistics, Mayfield Road, The King’s Buildings, Edinburgh EH9 352, Scotland, UK. Received 1997
30 September
1996; revised
2 May 1997; accepted
Appl. Math. Modelling 1997, 21: 569-578, September 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
23 May
given in the review articles by Fox’ and Voller et al.3 One successful method is the enthalpy method in which the discontinuous temperature derivative is by-passed by solving for the enthalpy instead of the temperature. The enthalpy method and its numerical implementation is discussed by Elliot and Ockendon,“ in the review article on finite element and finite volume methods for Stefan problems by Voller et al.3 and in the book on the numerical solution of moving boundary problems by Finlayson.’ Finite difference methods based on the enthalpy method work well for rectangular regions’ but are difficult to apply when the domain is not rectangular. Finite element methods overcome this drawback.3,5,7,s As an alternative to finite element methods, Coleman’ used a boundary integral method to solve phase change problems in which the temperature ahead of the phase change front is kept constant at the fusion temperature. This method has a great advantage in that it can be applied easily to regions of any shape and, furthermore, the Stefan condition can be fitted naturally into the numerical scheme without discretisation. The boundary integral method reduces the problem of solving for the temperature field in the two phases to solving for unknowns on the boundary of the material and on the phase change front only. Thus a two-dimensional problem is reduced to a one-dimensional problem on two curves, with a resulting great saving in calculation. Coleman’ applied his boundary integral method to solving phase change problems on squares and L-shaped regions, and it was found to give fast, accurate results in good agreement with other methods. The boundary integral method does, however, have a drawback in that it can only be applied to well-posed Stefan problems and not to supercooled situations, unlike the enthalpy method, which can be applied in both cases.3
0307-904x/97/$17.00 PII s0307-904x(97)00054-1
Crystal microstructure
prediction: S. J. Galloway
and N. F. Smyth
In the present work the boundary integral method of Coleman’ will be extended in two ways. The boundary integral formulation will be extended to include heat flow ahead of the phase change front, so that more general and realistic initial conditions can be used. Furthermore, Coleman’ used Newton’s method to find the position of the phase change front. It is found that for the case of the freezing of a liquid, when the ratio of the thermal diffusivities of the liquid and the solid phases is large, Newton’s method has difficulty converging due to the resulting small temperature derivatives. An example of this is the familiar case of water freezing to ice, as the ratio of the thermal diffusivities of water and ice is 8.9649. The numerical method presented in the present work overcomes this by using a bisection method to determine the position of the phase change front when Newton’s method has difficulty in converging. By symmetry, Newton’s method will also have difficulty in converging for a melting problem for which the ratio of the thermal diffusivities of the liquid and the solid phases is small. An application of the numerical method is made to predict the size of crystals formed when a material freezes, based on a simple model of crystal growth.“, l1 The predictions of this model are compared with experimental results.
2. Governing
equations
To be specific, let us consider the freezing of a liquid contained in an arbitrary two-dimensional region R. The fixed boundary of the region R will be denoted by Bf. Inside the region R there will be a moving phase change front, which will be denoted by BP. Furthermore let us denote the temperature of the liquid by T,(r, t) and the temperature of the solid by T,(r, t). Let the thermal diffusivity and thermal conductivity of the liquid and the solid be v, and K, and us and KS, respectively. The equation governing the heat flow in the solid is then
JT, - = vsV2Ts
(1)
at
for points r lying between the fixed boundary Bf and the phase change front B,,. The equation governing the heat flow in the liquid is
JT,
-
dt
= vlV2T,
(2)
for points r lying ahead of the phase change front BP.’ Initially the region R is totally filled with liquid with some temperature distribution, so that the initial condition is Tl(r, 0) = T,,(r) in R.
(3)
The boundary condition on the fixed boundary general, a Newton cooling condition n.VT, = p(T, - T,) on Bf,
570
Appl.
Math.
Modelling,
Bf is, in
(4)
1997,
Vol. 21, September
where n is the outward normal to the boundary and T, is the temperature of the surroundings. On the moving boundary BP the Stefan condition K,z-Ks-$=
-Lpn*V
(5)
holds, where n is the outward normal to the solid region, V is the velocity of the phase change front, p is the density of the solid, and L is the latent heat of fusion.’ The Stefan condition expresses energy conservation across the phase change front with energy being liberated due to the latent heat release as the liquid freezes. In deriving the Stefan condition it has been assumed that the density of the solid and the density of the liquid are the same for simplicity. If this assumption is not made, then the bulk motion of the solid and/or liquid must be included, which increases the complexity of the governing equations greatly. The governing equations (1) and (2), initial condition (31, and boundary conditions (4) and (5) can be nondimensionalised with respect to the parameter values for the liquid phase. Let temperature be measured with respect to the fusion temperature Tf of the liquid and scaled with respect to some typical temperature T, of the initial temperature distribution. A nondimensional temperature T’ is then (6)
If a is a typical length of the region R, then nondimensional space and time variables can be defined by r’
=
f
a
and
t’=-.
v
Using these nondimensional variables gives the nondimensional equations
JT, =
-
at
(7)
a2
and dropping primes
vV2Ts
for the solid phase between the fixed boundary phase change front BP and
JT, = V2T,
-
at
(8) Bf and the
(9)
for the liquid phase ahead of the phase change front BP. The boundary condition (4) becomes n-VT, = aPT,(T, - T,)
(10)
on the fixed boundary Bf, where T, is the nondimensional temperature of the surroundings. The Stefan condition (5) becomes
aT, --_K&& an
(11)
Crystal microstructure prediction: S. J. Galloway and N. F. Sm yth on the moving phase change front BP. Here the nondimensional diffusivity v and conductivity K are
v=%
and
GJr, r’, t, t’)
K K=&.
(12)
4
VI
Let G, be the Green’s function
= H(t - t’) The Stefan number
’ ev( - ,‘:,‘i:,) 47rv(t - t’)
a is
(16)
PV,L
(131
a=KIT,* The nondimensional equations with the Stefan condition (ll), the constant > 0, and the boundary constant < 0 on the fixed boundary formally to p = 03 in [lo]) have an dimension. This solution is
(8) and (9) together .. initial condition T,, = condition T, = Tb = Bf (which corresponds exact solution in one
so that
z
+ vV’*G,=
-S(t-t’)s(r-r’)
(171
where V’ refers to derivatives with respect to r’. Then it can be shown from equations (8) and (17) that
-$(G,cl
+ vV’.(T,‘V’G,
- G,V’T,‘)
= -T,‘i?(t-t’)6(r-r’) (14)
TI = T,, -
T,O
lc)\ erfc
(see Hill and Dewynne’). The phase change cated at x = e\ll, where 13 is the solution of
T,o
e-eZ/4 +
KTb
e-
front
where T,’ = T,(r’, t’). The fixed boundary Bf and the phase change front BP can be extended as surfaces Sr and S, in space-time, respectively, as shown in Figure I. The phase change front is then described by an equation of the form t-f(r)
is lo-
= 0
in space-time. given by
eZ/(4v)
(18)
(191
The outward unit normal to S, is therefore
(Vf, -1) = -
$k
This solution will be needed to start the numerical tion to be described in the next section.
3. Numerical
np=
(20)
d_’
(15) solu-
scheme
In this section a boundary integral scheme to numerically solve the system (8) and (9) together with the Stefan condition (11) will be described. This scheme is an extension of the method of Coleman’ to include a variable temperature ahead of the phase change front and the possibility that Newton’s method, which is used to determine the position of the phase change front, will not converge. Using the boundary integral method the system of equations (8) and (9) is converted into a pair of coupled integral equations for the temperature ahead of and behind the front. Let us consider the integral equation for the temperature behind the front first.
lime
t?L Y
X
Fixed
Figure
1.
tended
as surfaces
Sf
X-scaiOO
sf boundary
B,
S, and S,
Appl. Math. Modelling,
1997,
and
phase
change
front
B,
ex-
in space-time.
Vol. 21, September
571
Crystal microstructure
prediction:
S. J. Galloway and IV. F. Smyth
Integrating equation (18) over space-time between the surfaces S, and S,, and using the divergence theorem in space-time then gives -hT,(r,f)
vi:“‘/,;icz -G~z) ds’df’
=
(21)
The parameter A takes the value 1 at points away from the boundary and the value l/2 at regular boundary points. Projecting the integral over the surface Sp down onto the x -y plane as the surface B, finally gives -,+~,(r,t)
=
vl+‘Lf(T;z-Gsz)ds’df’ - u
/
GJr,r’,t,f(x’,y’)) B,!
xV’T;V’f
]If_f(x’,y+‘dy’.
(22)
This integral equation alone is not enough to determine the temperature behind the phase change front as the gradients VT,‘. V’f ItI=f(x’,y’) on the phase change front are not known. These unknown gradients are related to the same gradients ahead of the front by the Stefan condition (11). Hence the temperature field ahead of the front also needs to be considered. The equation governing the temperature field ahead of the phase change front is (9), and so the Green’s function appropriate to the region ahead of the front is
on noting that the outward normal to the surface Sp is now -np. The surface S, is the region inside the frxed boundary Bf. The integral equations (22) behind the phase change front and (26) ahead of the front together with the Stefan condition (11) completely determine the temperature field within the material. The integral equation (22) involves both the temperature T, and the temperature gradient ilTs/an at the fixed boundary Bs. Both of these will not be known and if, say, a fixed temperature boundary condition is imposed on Bf, then the gradient aT,,/an on Bf will have to be determmed as part of the solution. Hence this gradient, the gradient V’T/*V’flt~=f(X~,Y~)on one side of the phase change front (with this gradient on the other side of the phase change front determined via the Stefan condition [ll]) and the position of the phase change front itself have to be determined as part of the numerical solution. Once these are known the temperature at any point in the region can be determined using equations (22) and (26). The numerical method used to determine these unknowns will now be described. Let time be discretised into intervals of length At, so that ti = iht, and let the temperature at time ti be ryi(r) behind the front and Tli(r) ahead of the front. Then approximating the time integrals in (22) using the midpoint rule we have that the temperature on the boundary at time step k is
- ;
Tsk(=) k-l =
ZJA~c i=
?is4(r,r’,(k
/
(i - 1/2)At) =H(t
-t’)
4T(:_t’) exp[ -
- G,(r,r’,
:(-:;;j(23)
I
- 1/2)At,
1 &f
+
(k - 1,‘2)At,
(i - 1/2)At),,,
d Tii
G, then satisfies
$ +V”G,
+ v = - i?(t - t’)t?(r - r’)
so that from equation ;(G,T;)
(24) x GJr, r’, (k - 1/2)At,
(9)
+ V’.(T;V’G/
t’) dt’
(k-1'2)ArGs(r,r',
- G,V’T;)
= - T/6(t - t’)iS(r - r’).
(k - I/Z)
(25)
At, t’) dt' In a similar manner the front, integrating the surface S, gives -AT,(r,t)
- v
= $7
+
r’,O)G,(r,
r’, t,O) dw’ dy’
xV’T; V’f\r+~,~~)
572
Appl.
Math.
Modelling,
1997,
As in Coleman,”
h’ W
/
B,
GJr,r’,(k
- 1/2)At,f(x’,y’))
x V’T;a V’f[t~=f(r’,Y’)du’ dy’.
IB, Gr(r,r’,t,f(x’,Y’)) (26)
Vol. 21, September
ds' i
as for the temperature field behind this equation in space-time ahead of
(27)
the integral over the last time step has not been evaluated using the midpoint rule due to the singular nature of the Green’s function there. Instead,
Crystal microstructure over the last time step, the approximations
(k - 1/2)A1
/ (k-
r,r’, (k - 1/2)At,
ZZ-
27TV
T,‘(r’,O)G,( r,r’,(k-
t’) dt’
(r - r’1.n (r -
e-~r~r”2/(2~A~)
/
G)(r,r’, (k - 1/21At,f(x’,
y’))
4
~V’T~.V’fI,,=~(~,,~,)dx’dy’
yr12
and 1’2)AtGs(r,r’, (k - l/2)&, / (kUC- 1)At
t’) dt’
are used, where E, is the exponential integral.‘* Now that time has been discretised, space is also discretised. To evaluate the spatial integrals in equation (27) the boundaries Bf and BP are approximated by N-sided polygons and T, and ~?T,/an are approximated by their midpoint values on the sides. The integrals are then evaluated using the midpoint rule, except for when G, and aG,Jan’ become singular. This occurs when r lies on the midpoint of a side of Bf, in which case by expanding about this point it can be shown that the term aTs/dn’ in (27) has coefficient AL(y - 2 + log(AL2/(8vAt)))/(47rv), where y is the Euler constant and AL is the length of the side.” At these singular points the term in equation (27) involving T, on Bf is zero.’ The integral over the surface B, in equation (27) is evaluated in the following manner. The phase change front E$ is approximated by a polygon at each time step. By joinmg the vertices of these polygons for times 0 up to kAt together a set of quadrilaterals is formed that cover the surface B,. The integral over each of these polygons is then approximated by the value of the integrand at the centroid times the area of the quadrilateral, which is the midpoint rule in two dimensions. This approximation is valid for all quadrilaterals except those between the (k 11th and kth time step, due to G, becoming singular on the phase change front BP. As in Coleman,’ the integral over these quadrilaterals is approximated by
(area of quadrilateral)
1/2)At,O)dx’dy’
su +
1
and IV. F. Sm yth
(26) gives
0 = -/
a dn’GS(
1)Ar
equation
prediction: S. J. Galloway
x
V’T,‘.v’.f
/r’=(k-
1/2)Ar
4rvAt
(29)
where r’ is the quadrilateral centroid. The function f is approximated by (i - 1/2)At for any quadrilateral between the (i - 1)th and ith time step. To evaluate the unknown gradients V’T,‘-V’fIr,=fcx,,y,, the integral equation (26) is applied at the phase change front BP where T, = 0. Setting T, = 0 at time t = kAt,
(301
where r lies on the phase change front BP. The second integral in this expression is evaluated in the same manner as the similar integral in equation (27). The first integral over the initial condition is evaluated by dividing the region inside Bf up into a grid. The integral is then approximated by the midpoint rule as the sum of the area of each grid box times the value of T,‘(r’,O)G, at the centroid. The numerical procedure to solve equations (27) and (30) then works in the following manner. For definiteness, let us suppose that the boundary condition on Bf is the fixed temperature boundary condition T, = Tsb < 0. We then assume that we know the temperature gradients V’T,‘.V’f(t,=f~x~,Y,) on the phase change front for times 0 up to (k - 1)At. The position of the phase change front at time kAt is then approximated by extrapolating from its previous positions. The linear system arising from the integral equation (30) can then be solved for the gradients V’T,’. V’f]r,=,-(x,,y,) on the front at time kAt. From these gradients and the Stefan condition (11) the gradients V’T,‘v’fIr,=f(x’,Y’) behind the phase change front can then be determined. Substituting these gradients behind the phase change front into the linear system resulting from the integral equation (27) the unknown normal derivatives aT,/an on the boundary Bf can then be found. Once these normal temperature derivatives are known, Newton’s method with the temperature TY given by the discretised form of equation (22) can then be used to find where the temperature T, is zero. This then gives a revised position of the phase change front at time kbt. The whole process is repeated until the difference between the positions of successive estimates of the front is within a prescribed tolerance. To start the process off the one-dimensional solution (14) and (15) is used to give the front position and gradients VT;. V’f If8=P(~‘,~‘) at time At. This will give a good approximation to these quantities at time At if At is small. To keep the phase change front smooth, when the iterations have converged, the points on the front are passed through a three-point smoothing scheme. It was found that the position of the phase change front could be successfully determined using Newton’s method to search for when T, = 0 unless v was large. When v is large, the derivative of the temperature T, is small (see the onedimensional solution [14]), and this causes Newton’s method to have difficulty in converging. To enable solutions to be found when v is large a bisection method was used to determine the position of the phase change front (searching for when T, = 0) when Newton’s method had difficulty in converging. In applying the numerical scheme outlined in this section, Newton’s method was used to determine the position of the phase change front, unless it did not converge within a fixed number of iterations. In
Appl. Math.
Modelling,
1997,
Vol. 21, September
573
Crystal microstructure
prediction: S. J. Galloway and N. F. Smyth
this case the front finding method was then switched over to the bisection method.
the phase change front will start to move instantaneously. As a check on the numerical scheme the boundary Bf was taken to be a square of side length 2, and To was set equal to zero. This then reduces to one of the examples considered by Coleman’ as the temperature ahead of the phase change front is constant at Tl = 0. With Tb = - 1, Q = 1, or= 1, and K = 1 the time for the complete freezing of the region was found to be 0.474695, which agrees with the value 0.475 found by Coleman.’ Figure 2 shows the solution for Tb = - 1.0 and T,, = 0.65 with the parameter values v = 8.9649, CY= 3.5732, and
4. Results The boundary integral method outlined in the previous section was applied to a number of examples of the freezing of a region. For simplicity the boundary condition applied on Bf was taken to be T, = Tb < 0 and the initial condition was taken to be T, = T,, a 0, where Tb and T,, are constants. With these initial and boundary conditions
1 0.5 T
0 -0.5 -1
0
0
x
.I
-
.6
-
..
-
.2 -
.
I-
-.2
-
-.1
-.6
-.a
-
Cc) Figure 2. (a) Contour plot of the temperature for (x, y) in the fourth quadrant of a square of side length 2 at time ~0.0875 for Tb= - 1 .O, Tc=O.65, v=8.9649, o=3.5732, and K=4.1516 with T,=2O”C. (b) Surface plot of the temperature for Ix, y) in the fourth quadrant of a square of side length 2 at time t=0.0875 for the parameters of (a). (c) Phase change front positions for a square of side length 2 at time intervals of 0.0175 up until complete freezing for the same parameter values as (a).
574
Appl.
Math.
Modelling,
1997,
Vol. 21, September
Crystal microstructure K = 4.1516. These parameter values are those appropriate for water freezing to ice with T, = 20”C.‘3 For these parameter values, Newton’s method had difficulty in converging, and a bisection method was used to determine the position of the phase change front. The region inside Bf was taken to be a square of side length 2. Figures 2(a) and 2(b) show contour and surface plots for the temperature field at time t = 0.0875. The phase change front can be clearly seen, particularly in the contour plot. The time evolution of the phase change front is shown in Figure 26~) at time intervals of 0.0175. It can be seen that the front is evolving to a circular shape as time increases. The freezing time for the parameter values of Figure 2 is 0.325134, which is smaller than the freezing time for the example of Coleman.’ This is due to the larger values of the thermal conductivity K and the thermal diffusivity V, which means that heat is more easily transported to the boundary Bf. The Stefan number a is larger than for the example of Coleman,’ which means that more heat is liberated on freezing. This and the larger value of the initial temperature TO would cause the freezing time to increase over the value of Coleman.’ However the larger values of the thermal diffusivity and the thermal conductivity have the greater effect, and the freezing time decreases. The fact that the freezing time is not much reduced over the value of Coleman’ shows that the increased values of V, K, (Y, and T,, have effects that nearly balance out. Figure 3(a) shows the variation of the freezing time tf with TOfor Tb = - 1.0, u = 8.9649, CY= 1.0 and K = 4.1516. These parameter values are those for water, except for (Y. For this high value of v Newton’s method has difficulty in converging, and the bisection method is used to locate the position of the phase change front. As expected the freezing time increases with the initial temperature TO, with the increase being nearly linear. The increase in the freezing time tf as the initial temperature T,, ranges from 0 to 1 is not very large because the boundary temperature is relatively low at Tb = - 1. This low boundary temperature and the large relative thermal diffusivity Y of the solid phase then dominate the heat flow. Figure 3(b) shows the variation of the freezing time tf with the initial temperature TO for Tb = - 1.0, v = 0.5, (Y= 3.5732, and K = 4.1516. The value of (Y used for this figure is the value appropriate for water. The freezing time again increases with increasing initial temperature T,, in a nearly linear fashion. The freezing times in Figure 3(b) are larger than those in Figure 3(a) due both to the smaller value of v and to the larger value of (Y. The smaller value of v means that heat diffuses out of the region at a slower rate, and the larger value of (r means that more heat is liberated on freezing, which then has to diffuse out of the region. Figure 4 shows the variation of the freezing time tf with T,, for T,, = - 1.0 and the water parameters (v = 8.9649, (Y= 3.5732, and K = 4.1516). The freezing time again increases with T,,, as expected, but in this case the increase is highly nonlinear, with the rate of increase slowing with increasing TO. The freezing times are between those of Figures 3(a) and 3(b). The freezing times are increased over those of Figure 3(a) due to the larger value of (Y
prediction:
S. J. Galloway and N. F. Smyth
I
I
0.15 _ 0.145 0.14 0.135 0.13 tt
0.125 0.12 -
0.54;
0.2
0.4
.l-
0.6
I 1
0.8
‘0
W Figure 3. (a) Plot of freezing time t, as a function of To for a square of side length 2 for the parameter values Tb= - 1.0, v=8.9649, a= 1.0, and K=4.1516 with T,=2O”C. (b) Plot of freezing time t, as a function of To for a square of side length 2 for the parameter values Tb= - 1.0, v=O.5, (~=3.5732, and K=4.1516 with T,=2o”C.
leading to more heat liberated on freezing and are decreased over those of Figure 3fb) due to the larger value of v leading to more rapid loss of heat.
5. Crystal formation The numerical method outlined in Section 3 can be used, together with a simple model of crystal formation, to predict the sizes of crystals formed when a liquid freezes. If the growth of the crystal is radially symmetric, Frank” and Coriell and McFadden” showed that the size of the crystal formed when a liquid freezes is given in dimensional units by
(31)
Appl.
Math.
Modelling,
1997,
Vol. 21, September
575
Crystal microstructure
Figure
4.
square
of side
v=8.9649,
prediction: S. J. Galloway and N. F. Smyth
Plot of freezing length
a=3.5732,
K are those for water
2 for
time the
t,
and K=4.1516 freezing
as a function
parameter
values
(the values
to ice with
of
T, T,=
for
a
- 1 .O,
of Y. a, and
T,= 20°C).
Here K, is the thermal conductivity of the liquid, L is the latent heat, T, is the freezing temperature of the liquid, and T, is the temperature far ahead of the growing crystal. It should be emphasised that this expression is valid for the stable freezing of a liquid and is not valid for the dendritic growth that occurs when the liquid is undercooled.14 The size expression (31) holds for a single crystal that has been growing for a time t and was found to give crystal sizes in accord with experimental results.” When a region freezes, as in the numerical solutions of the previous section, a large number of crystals form, the crystal at a given point being formed when the phase change front passes that point. Hence to apply the crystal growth expression (31) to the freezing of a region, this expression needs to be applied to each crystal formed (i.e., at each point of the region), where the time t is the time taken for a given crystal to form. In the present case the crystal size expression (31) is most easily applied if the time t is replaced by the front velocity. This is in agreement with models of microstructure formation that show that the size of microstructural features is primarily determined by the speed of the phase change front.‘5,‘6 If the phase change front were planar, then the nondimensional front position is given by x = 0fi (see equation [15]), so that the nondimensional front velocity is I/ = 0/(2fi). Eliminating time in favor of front velocity in the crystal size expression (31) and using equation (7) to convert to dimensional variables gives the alternative crystal size expression
(32)
where change planar phase tional
576
cI is the specific heat of the liquid. While the phase fronts considered in the present work are not it was noted by Coleman’ that even for nonplanar change fronts the position of the front is proporto fi to a good approximation. Therefore the
Appl.
Math.
Modelling,
1997,
Vol. 21, September
crystal size expression (32) will be used in the present work on the freezing of two-dimensional regions. The crystal sizes formed when the liquid in a region freezes are then calculated in the following manner using the boundary integral method of Section 3. Using the positions of the phase change front calculated by the numerical scheme of Section 3 the normal velocity V, of a point on the front at a given time step can be calculated. This normal velocity V, is then used as the velocity V’ in the crystal size expression (32). The crystal size expression (32) is now employed to calculate the crystal size in each quadrilateral used to calculate the surface integral B,. At each time step the parameter 0 is determined from the transcendental equation (151, with Tb set equal to the boundary temperature and T,, set equal to the temperature at the centre of the square of Section 4. While this transcendental equation is strictly valid only for a planar phase change front the crystal size was found to be only weakly dependent on the value of 19. Indeed the crystal sizes were found not to change to a significant extent if 0 was set equal to its initial value. The value of 0 was calculated using equation (19, so that the effect of the changing temperature of the liquid could be included. The temperature T, is taken to be the boundary temperature T,, and, for water, the freezing temperature is T, = 0, these temperatures being made dimensional using equation (6). Hence the crystal size in each quadrilateral used to calculate the surface integral B, can then be calculated. Adding up all such local crystal sizes over all positions of the phase change front until the region has completely frozen and dividing by the area of the region then gives the average crystal size $. Figure 5(u) shows a plot of the mean crystal size ?? as a function of the freezing time tf for Tb = - 1.0 and the water parameters with T, = 20°C and a = 5.0 cm. The mean crystal size is nearly linear in tf, except near tf = 0.34 (ie. except for initial temperatures near To = 1.0: see Figure 4). The rapid increase in the mean crystal size near To = 1 is related to the slow increase of the freezing time with initial temperature To near To = 1, which can be seen in Figure 4. The increase of the mean crystal size with freezing time is expected since a longer freezing time means a smaller mean front speed, which in turn leads to a larger crystal size. The actual mean crystal sizes, while a little on the high side, are in reasonable accord with the ice crystal sizes reported by Coriell and McFadden.” One reason for the mean crystal sizes being somewhat higher than experimental values is that the fluid flow ahead of the phase change front has been neglected. This flow is expected to produce faster freezing times and hence smaller mean crystal sizes. The inclusion of fluid flow is the subject of current work. The experimental measurements of Miyawaki et al. ” of the freezing of soy protein curd show a basically linear dependence of the mean crystal size on freezing time (with a large amount of experimental scatter), which is in qualitative agreement with our simple model of crystal formation on freezing. Unfortunately the exact experimental conditions were not reported in this work, so that quantitative comparisons with their results cannot be made. Figure 5(b) shows the
Crystal microstructure
prediction: S. J. Galloway and N. F. Sm yth
change front is taken into account, extending the work of Coleman,’ which considered heat flow behind the phase change front only. It is found that this boundary integral method formulation of the Stefan problem has the advantages of (i) being a fast numerical method that yields accurate solutions and (ii) that it can easily incorporate nonrectangular boundaries. When parameter values appropriate for water were used it was found that the standard method for determining the front position, based on Newton’s method, was unstable, and so a new method for determining the front position, based on the bisection method, was developed for these parameter values. Using a simple model for the growth of a crystal from a freezing liquid the numerical method was extended to predict the sizes of the crystals in the crystalline microstructure formed when a body of liquid freezes. The predictions of this model for crystal formation were found to be in qualitative agreement with experimental results. While this model of crystal size was developed in the context of water it can easily be applied to other areas in which the prediction of crystal size is important, for example the growing of silicon crystals for the computer industry.
7.6 7.4 7.2 76.6 6.6 6.4 -
6.6 r 6.5 6.4 6.3 6.2 . 6.1 .
Acknowledgments
6-
The authors the referees presentation EPSRC for dentship.
5.9 5.6 5.7 -/ 5.6 I -1
-0.9
-0.6
0.7
-0.6
-0.5
-0.4
would like to express their appreciation to for pointing out references that improved the of this paper. S.J.G. would like to thank the funding through an earmarked Ph.D. stu-
I -0.2
-0.3
5
(W
Figure 5. the
freezing
v=8.9649, 5.0
(al
Plot of the
time
t,
a=3.5732,
cm. (b) Plot of the
boundary v=8.9649.
temperature (~=3.5732,
for
mean T,=
and mean
crystal
- 1.0 K=4.1516 crystal
Tb for T,=O.l and
K=4.1516
and
size the with
3 as a function water
size 3 as a function with
1. Hill, J. M. and Dewynne, J. N. Heat Conducfion. Blackwell, Oxford, UK, 1987 2. Fox, L. The Stefan problem: Moving boundaries in parabolic equations. A Surrvy of Numerical Methods for Partial Differential Equafions, eds. 1. Gladwell and R. Wait, Clarendon Press, Oxford, UK, 1979 3. Voller, V. R., Swaminathan, C. R., and Thomas, B. G. Fixed grid techniques for phase change problems: A review. Inr. J. Num.
parameters
T,=2o”C
and the water
References
of
and
a=
of the
parameters
T,=2O”C
and
a=
5.0 cm.
Method. Engin. 1990, 30, 875-898 4. Elliot, C. M. and Ockendon, J. R. Weak
dependence of the mean crystal size 5 on the boundary temperature Tb for fixed initial temperature T,, = 0.1 and the water parameters with T, = 20°C and a = 5.0 cm. The mean crystal size decreases with decreasing boundary temperature, which is expected since a lower boundary temperature implies a larger mean front speed, which leads to a smaller crystal size. The mean crystal size dependence shown in Figure S(b) shows a good qualitative similarity to the experimental measurements for soy protein curd obtained by Ref. 17.
5. 6.
7.
8.
9.
6. Conclusions
and variational
methods
for moving boundary problems. (Pitman Advanced Publishing Program) Res. Notes Math. 1982, 59 Finlayson, B. A. Numerical Methods for Problems with Mol+zg Fronts. Ravenna Park Publishing, Seattle, WA, 1992 Voller, V. and Cross, M. Accurate solutions of moving boundary problems using the enthalpy method. Int. J. Heat Mass Transfer 1981, 24, 545-556 Lynch, D. R. Unified approach to simulation on deforming elements with application to phase change problems. J. Camp. Phys. 1982, 47, 387-411 Chow, P. and Cross, M. An enthalpy control-volume-unstructured-mesh (CV-UM) algorithm for solidification by conduction only. ht. J. Num. Method. Eng. 1992, 35, 1849-1870 Coleman, C. J. A boundary integral formulation of the Stefan problem. Appl. Math. Mode&g 1986, 10, 445-449
10. Frank, F. C. Radially symmetric phase growth controlled diffusion. Proc. Roy. Sot. Land. A 1950, 201, 586-599
It has been shown that a boundary integral method can be used to find a numerical solution of the Stefan problem when heat flow both ahead of and behind the phase
11. Coriell,
S. R. and
McFadden,
G. B. Morphological
Handbook of Ctysral tirowth I, Fundamentals, and Stability, ed. D. T. J. Hurle, North-Holland,
Appl.
Math.
Modelling,
1997,
by
Stability.
Part B: Transport
Amsterdam,
Vol. 21, September
1993
577
Crystal microstructure
prediction:
S. J. Galloway and N. F. Smyth
12. Abramowitz, M. and Stegun, I. Handbook of Mathematical Functions. Dover Publications, New York, 1972 13. Lide, D. R., ed. CRC Handbook of Chemistry and Physics, 74th edition, Chemical Rubber Publishing Company, Boca Raton, FL, 1993 14. Langer, J. S. Instabilities and pattern formation in crystal growth. Reo. Mod. Phys. 1980, 52(l), l-28 15. Rappaz, M., Carrupt, B., Zimmermann, M., and Kurz, W. Nu-
578
Appl.
Math.
Modelling,
1997,
Vol. 21, September
merical-simulation ment of materials.
of eutectic solidification in the laser He/u. Phys. Acta 1987, 60, 924-936
16. Rappaz, M. Modelling of microstructure formation tion processes. fat. Mater. Reu. 1989, 34(j), 93-123
treat-
in solidifica-
17. Miyawaki, O., Abe, T., and Yano, T. Freezing and ice structure formed in protein gels. Biosci. Biotech. Biochem. 1992, 56(6), 953-957