Ind. Eng. Chem. Res. 2008, 47, 5087–5091
5087
Growth Constants in Solidification Saurabh Agarwal, Lewis E. Johns, and Ranga Narayanan* Department of Chemical Engineering, UniVersity of Florida, GainesVille, Florida 32611
In phase-change problems, there are domain and surface contributions to the dynamics of a moving front. It is ordinarily assumed that the surface contributions are in complete control. This is likely to be true in the case of planar fronts, but it may or may not be true in the cylindrical case. We study the growth rates of disturbances imposed on the surface of a cylindrical solid in equilibrium with its melt and deduce the conditions under which time derivatives of bulk-phase variables are important. They are important in the case of a whisker and are even more important if a solute is added. 1. Introduction Time derivatives appear in many places in the models that are used to predict the growth rates of disturbances in phasechange problems where a front advances as one phase turns into another. They appear in the heat and mass balances across the surface, and they appear in the diffusion equations, accounting for solute and heat migration in the bulk phases. In the case where the growth of a perturbation imposed on a base state is the issue, it is often assumed that the time derivatives in the bulk phases are not important. This is the case, in many studies, of growth rates in electrodeposition (e.g., see Kahanda et al.1) and in solidification (e.g., see Mullins and Sekerka2). This is a fair approximation in the case of a planar front, especially if the wavenumbers of interest lie near the greatest growth rate, which is often the case. For instance, if a pure solid is in equilibrium with its melt across a planar interface and a deflection is imposed on the surface, it decays at a rate σ, where σ is given by σ ) -Lcapk + √Lcap2k2 + 1 σS and σS is the rate at which the perturbation would decay if the time derivatives on the domain were neglected. In this formula, the capillary length (Lcap) is given by Lcap )
( Lγ )( Lκ ⁄ k
TM cond
)
where γ is the surface tension, L the latent heat per unit volume, TM the melting point, κ the thermal diffusivity, and kcond the thermal conductivity; k denotes the wavenumber of the perturbation. It is the small value of Lcapstypically of the order of 10-8 m in solidificationsthat makes the time derivatives at the surface dominate over the time derivatives in the bulk phases. For ordinary lengths, the values of the wavenumber in the neighborhood of the wavenumber at the greatest growth rate are not large enough to offset the small values of Lcap. Our aim is to present a case where, in a physically realizable experiment, domain effects must be taken into account. The problem is a whisker growing from a subcooled melt; but, for simplicity, we take the case of a whisker of radius R0 in equilibrium with its melt, upon which a surface displacement is imposed. * To whom correspondence should be addressed. Tel.: +1 352 392 9103. Fax: +1 352 392 9513. E-mail address:
[email protected].
Figure 1. Sketch used to explain the instability of a whisker in equilibrium with its melt.
A cylindrical front is more interesting than a planar front. In the planar case, the curvature is stabilizing, with temperatures falling at crests and rising at troughs. This is not so in the cylindrical case. An axisymmetric disturbance is of most interest. Then, a series of crests followed by troughs in the axial direction also implies diameter changes along this direction. A sketch of this is presented in Figure 1. In the figure, the cylinder is at equilibrium with its melt at temperature TM - (γTM/L)(1/R0), whereupon a sinusoidal disturbance of the surface is imposed, changing the temperature along the surface, sometimes in a way to weaken the imposed disturbance, other times to strengthen it. The curvature in the longitudinal direction decreases the temperature at a crest, increasing it at a trough. This is stabilizing, melting back a crest, building up a trough, as heat is conducted from the bulk phases to a crest, or away from a trough into the bulk phases. The strength of this effect is dependent on the wavelength of the disturbance, being stronger the shorter the wavelength. Opposing this, and independent of the wavelength, is the transverse diameter effect, raising the temperature at a crest and reinforcing the crest, lowering it at a trough. These effects come into balance at a critical wavelength; at longer wavelengths, the disturbance is unstable and a cylindrical rod of solid breaks. Our interest is in the growth rate of an unstable disturbance and, in particular, in the conditions under which the usual assumption that the bulk dynamics is unimportant fails. The cylindrical case offers us a possibility of seeing this, no matter how small the value of Lcap, because a cylindrical solid introduces its radius as an additional controllable variable.
10.1021/ie070838w CCC: $40.75 2008 American Chemical Society Published on Web 01/03/2008
5088 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008
The equilibrium base solution is denoted by the subscript zero. It is given by T0 ) T/0 ) TM -
γTM 1 L R0
(5)
where R0 is the radius of the cylindrical whisker. To find the growth rate σ corresponding to a displacement of the surface of wavelength 2π/k in the z-direction, we must solve the eigenvalue problem ∂T1 ) κ∇2T1 ∂t
(6)
and ∂T1/ ) κ∇2T1/ (7) ∂t where T1 is bounded at r ) 0 and T1* is bounded as r f ∞, and where, at the surface r ) R0, we have T1 ) T/1 )
γTM 2H1 L
(8)
and ∂T1 ∂T/1 - kcond ) Lu1 ∂r ∂r The formulas for 2H1 and u1 are
(9)
kcond
Figure 2. Diagram of a whisker in equilibrium with its melt.
2H1 )
2. The Problem The sketch in Figure 2 illustrates a whisker in equilibrium with its melt at a uniform temperature below the equilibrium value TM, corresponding to a plane surface. Some notation is given in the figure, as well as the formulas for the surface normal, twice its mean curvature, and its speed. We assume that the displacement of the surface is axisymmetric. The domain equations are given by ∂T ) κ∇2T ∂t
(1)
R1
∂2R1
R0
∂z2
+ 2
(10)
and ∂R1 (11) ∂t where R1 denotes the displacement of the surface. These equations obtain on introducing small corrections to the base state variables into the nonlinear equations and dropping quadratic and higher-order terms. Solving these equations by writing u1 )
and
T1 ) Tˆ1(r) cos(kz) exp(σt)
(12)
∂T* ) κ∇2T* (2) ∂t where T is bounded at r ) 0 and T* is bounded as r f ∞. At the surface r ) R(z,t), phase equilibrium implies
T/1 ) Tˆ/1(r)
(13)
( )
γTM 2H (3) T ) T * ) TM + L where the term (γTM/L)2H corrects the equilibrium temperature for curvature and is at the heart of the problem. The heat balance at the surface is given by kcondn · ∇ T - kcondn · ∇ T * ) L n · u f
f
f f
cos(kz) exp(σt)
and R1 ) Rˆ1 cos(kz) exp(σt) we find that the eigenvalue σ is given, implicitly, by
(
) [
Lcap 1 - k2R02 λR0 I′0(λR0) K′0(λR0) σ ) R0 kR0 kR0 I0(λR0) K0(λR0) κk2
]
(15)
where λR0 and Lcap are given by
(4)
The assumption has been made that the thermal properties of the solid and liquid phases are nearly the same. They have been written as solid-phase properties. This assumption is not essential, and the work goes through if the thermal properties differ; however, the resulting formulas are not so simple. Thisassumption is, in fact, true in many systems (cf. Wollkind and Notestine3). We have also assumed the densities of the solid and liquid phases to be the same; otherwise, we would have to account for a small liquid side velocity at the surface. We neglect this.
(14)
λR0 ) kR0 and Lcap )
(
1+
σ κk2
(16)
)
(17)
TM γ L Lκ ⁄ kcond
We do not introduce scaled variables, which would make our formulas look a little better. This is due to the fact that we could not find a scaling of the nonlinear equations that leads to the correct conclusions about the importance of the various time
Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5089
The model we use is that of Mullins and Sekerka2 and, to account for the effect of the solute, one correction to the temperature problem must be made. A term must be added to eq 3 to account for a reduction of the equilibrium temperature due to the presence of the solute. Hence, in place of eq 3, we now have, at the surface, T ) T/ ) TM +
Figure 3. Curves of constant σ/σS for the solidification of a pure material.
derivatives. Yet, the reader might like to see the results in some scaling; hence, we do not cancel a number of R0’s whereupon the dimensionless groups σ/(κk2), kR0, Lcap/R0, etc. appear. Our main interest is in the range 0 < kR0 < 1, where there is one positive growth constant, the growth constant that would be of interest in an experiment. We first inquire as to the value of σ if the domains have no dynamic role. Denoting this value of σ by σS, we have
(
)[
]
(18)
I′0(λR0) K′0(λR0) I0(λR0) K0(λR0) I′0(kR0) K′0(kR0) I0(kR0) K0(kR0)
(19)
Lcap 1 - k2R02 I′0(kR0) K′0(kR0) ) R0 kR0 I′0(kR0) K0(kR0) κk2 σS
[
Then σ/σS is given, implicitly, by
σ ) σS
( )( ) 1+
σS
σ σS κk2
where
]
1+
2
σ σS
κk The first factor on the right-hand side of eq 19 appears to tell us that σ/σS will be ∼1 unless σS/(κk2) is of the order of 1 or greater, and, from eq 18, we can see that σS/(κk2) is a multiple of Lcap. Now the factor[(I0′(kR0)/I0(kR0)) - (K0′(kR0)/K0(kR0))] increases as kR0 decreases, with the increase becoming rapid when kR0 < 1/4, and this is due to -(K′0(kR0)/K0(kR0)), not (I′0(kR0)/ I0(kR0)), which is negligible, implying that time derivatives in the solid phase are not important. Figure 3 presents our results as curves drawn at constant values of σ/σS on a graph whose axes are kR0 and R0/Lcap. The information presented is this: analysis of experiments is simplest if σ/σS is ∼1. This figure then tells us, given a value of σ/σS, that all points on the plane below and to the left of the curve have values of σ/σS higher than the specified value. For example, given Lcap, then to any R0, the range of kR0 for which σ/σS exceeds the specified value can be read from the graph. Assuming that the growth is in the form of a single crest, we then have the lengths of the whisker that we ought to avoid at a given value of R0, viz, L ) (π/k), whence the range of πR0/L is known. 3. The Case of a Dilute Binary Alloy The values of R0 at which time derivatives of domain variables must be taken into account increase if a dilute solute is added to the whisker studied in section 2.
(20)
where c* denotes the solute concentration in the melt and m denotes the slope of the liquidus line (m < 0). The solute concentration must then satisfy ∂c/ ) D/∇2c/ ∂t
(21)
on the melt side, and, at the surface, c/ )
cS kdist
(22)
and -D/n · ∇ c/ ) - [cS - c/]n · u f
( )( )
λR0 ) kR0
σS
γTM 2H + mc/ L
f f
(23)
where cS denotes the solute concentration in the solid and kdist denotes the equilibrium distribution ratio. The base solution then corresponds to uniform temperatures and compositions, viz, c0/ )
cS kdist
(24)
γTM 1 + mc/0 L R0
(25)
and T0 ) T/0 ) TM -
where R0 is the radius of the cylindrical whisker at equilibrium. Again, we seek the growth rate σ of a disturbance of wavelength (2π/k). We do not perturb the solute concentration on the solid side, assuming solid-phase diffusion of the solute to be very slow. This might present a problem, due to the small values of R0 of interest, were it not for the fact that the values of k of interest are such that kR0 < 1, whence the concentration variation in the liquid is confined to an equally thin region near the surface. The equations satisfied by T1 and T*1 are as before except that eq 8 must be replaced by T1 ) T/1 -
γTM 2H1 + mc/1 L
(26)
and to these equations must be added ∂c/1 ) D/∇2c/1 ∂t
(27)
on the melt side (r > R0) and -D/
∂c/1 ) (1 - kdist)c/0u1 ∂r
(28)
at the surface (r ) R0). Solving these equations as before, we find that the eigenvalue σ is now given by
5090 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008
( [
) [ ] [
]
Lcap 1 - k2R02 λR0 I′0(λR0) K′0(λR0) R0 kR0 kR0 I0(λR0) K0(λR0) σ ) 2 I′0(λR0) K′0(λR0) κk kcond m(kdist - 1)c/0 λR0 I0(λR0) K0(λR0) 1+ / K′0(µR0) L µR0 D K0(µR0)
[
]
]
(29)
where µ is given by
µR0 ) kR0
1+
σ D/k2
(30)
and where the second term in the denominator on the righthand side is due to the presence of the solute. Again positive values of σ occur for 0 < kR0 < 1 and we have k < λ < µ. There are five places where σ appears in the model: there are three domain terms and two surface terms in which σ appears. Dropping the domain σ terms and denoting the result σS, we have
σS 2
( [
Lcap R0
)
κk
1+
kcond L
)[
[
] ]
] [
where σ/σS is given by
[ [
]
1 - k2R20 I′0(kR0) K′0(kR0) kR0 I0(kR0) K0(kR0) I′0(kR0) K′0(kR0) m(kdist - 1)c/0 I0(kR0) K0(kR0) / K′0(kR0) D K0(kR0)
] ]
I′0(λR0) K′0(λR0) I0(λR0) K0(λR0) σ ) × σS I′0(kR0) K′0(kR0) kR0 I0(kR0) K0(kR0) λR0
1+
1+
kcond L
[
[
[]
]
I′0(kR0) K′0(kR0) I0(kR0) K0(kR0) K′ kR D 0( 0) K0(kR0) I′0(λkR0) K′0(λR0) / m(kdist - 1)c0 λR0 I0(λR0) K0(λR0) / K′ µR µR D 0 0( 0) K0(µR0)
kcond m(kdist - 1)c/0 / L
[
] [ [
(31)
]
]
[
]
Table 1. Properties of Succinonitrile property
value
melting point, TM thermal conductivity of solid, kcond thermal diffusivity of solid, κ surface tension, γ latent heat of fusion, L
330 K 0.23 W/(m K) 1.2 × 10-7 m2/s 9 × 10-3 J/m2 0.5 × 108 J/m3
Table 2. Properties of SuccinonitrilesAcetone Alloy property
value
equilibrium distribution coefficient, kdist liquidus slope, m liquid mass diffusivity, D*
0.1 mol %/mol % -2.2 (K/mol %) 1.3 × 10-9 m2/s
than 5 µm if solute is added, σ/σS will be greater than 1.1. Third, we add the observation that the heat required to melt a curved surface differs from the heat required to melt a plane surface, being less, and, in fact, much less for a cylinder of radius 10-6 m (see Lai et al.6 for experimental results and Wollkind and Notestine3 for a theoretical explanation). In our case, the base surface is curved and Lcap might be 1 or 2 orders of magnitude higher than 10-8 m. This increases the area of the kR0 - (R0/ Lcap) plane, where time derivatives in the bulk phases become important. 4. Conclusion
]
(32)
In these formulas, all of the terms in I′0/I0 can be dropped, visa-vis, those in K′0/K0, as long as 0 < kR0 < 1/4. Curves of constant values of σ/σS, with solute added, are given in Figure 4 as R0/Lcap vs kR0. They are drawn for kcond m(kdist - 1)c/0 ) 0.7 L D/
Figure 4. Curves of constant σ/σS, where a solute has been added.
(33)
corresponding to an equilibrium solute concentration of 0.1 mol %. Physical properties are given in Tables 1 and 2 (cf. Fang et al.4 and Li and Beckermann5). First, we see that time derivatives on the domain are more important in the presence of a solute than in the earlier solutefree case. Second, to illustrate this, we give some results at σ/σS ) 1.1, Lcap ) 10-8 m, and R0 ) 10-6 m. For whisker lengths of greater than 40 µm in the case of a pure material, or greater
We have presented results that indicate the importance of time derivatives in the bulk phases, vis-a-vis, time derivatives at the surface in the growth of a perturbation imposed on a phase boundary. Ordinarily, surface time derivatives control the growth, due to the small capillary lengths that are usually encountered in phase-change problems; however, cylindrical geometries offer the prospect of increasing the importance of bulk-phase time derivatives. Indeed, at a cylinder radius of 10-6 m, we can find whisker lengths where the bulk phase has a non-negligible role in physically realizable experiments, i.e., reasonable whisker lengths. A radius of 10-6 m is of the order of the radii of the whiskers studied by Kirill et al.7 We also find that the addition of solute increases the importance of domain dynamics. Experiments on solidification in cylindrical geometries have been carried out by Hardy et al.8 However, these are not equilibrium experiments and, hence, there is no steady base state. Our work ought to lead to experiments on equilibrium systems. Nomenclature D* ) liquid mass diffusivity (m2/s) 2H ) twice the mean curvature, as defined in Figure 2 (m-1) (I0, K0) ) modified Bessel function of order zero k ) wavenumber of perturbation (m-1)
Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5091 kdist ) equilibrium distribution coefficient (mol%/mol%) kcond ) thermal conductivity (W/(m K)) Lcap ) capillary length, as defined by eq 17 (m) m ) slope of liquidus line (K/(mol%)) b n ) unit normal vector defined in Figure 2 R0 ) equilibrium radius of a whisker (m) (r, z) ) unscaled cylindrical coordinates TM ) melt temperature (K) Greek Letters γ ) surface tension (J/m2) κ ) thermal diffusivity (m2/s) (λR0) ) defined by eq 16 (µR0) ) defined by eq 30 σ ) growth constant of perturbation, as defined by eqs 12-14 (s-1) σs ) growth constant when time derivatives in domain differential equations are neglected (s-1) Other Characters L ) latent heat of fusion per unit volume (J/m3) 0 ) equilibrium state 1 ) perturbed state
Literature Cited (1) Kahanda, G. L. M. K. S.; Zou, X.-q.; Farrel, R.; Wong, P.-z. Columnar growth and kinetic roughening in electrochemical deposition. Phys. ReV. Lett. 1992, 68, 3741. (2) Mullins, W. W.; Sekerka, R. F. Stability of a planar interface during solidification of a dilute binary alloy. J. Appl. Phys. 1964, 35, 444. (3) Wollkind, D. J.; Notestine, R. D. A nonlinear stability analysis of the solidification of a pure substance. IMA J. Appl. Math. 1981, 27, 85. (4) Fang, Q. T.; Glicksman, M. E.; Coriell, S. R.; McFadden, G. B.; Boisvert, R. F. Convective influence on the stability of a cylindrical solidsliquid interface. J. Fluid Mech. 1985, 151, 121. (5) Li, Q.; Beckermann, C. Modeling of a free dendritic growth of succinonitrilesacetone alloys with thermosolutal melt convection. J. Cryst. Growth 2002, 236, 482. (6) Lai, S. L.; Guo, J. Y.; Petrova, V.; Ramanath, G.; Allen, L. H. Sizedependent melting properties of small tin particles: nanocalorimetric measurements. Phys. ReV. Lett. 1996, 77, 99. (7) Kirill, D. J.; Davis, S. H.; Miksis, M. J.; Voorhees, P. W. Morphological instability of a whisker. Proc. R. Soc. London, Ser. A 1999, 455, 3825. (8) Hardy, S. C.; Coriell, S. R. Morphological stability and the icewater interfacial free energy. J. Cryst. Growth 1968, 3, 569.
ReceiVed for reView June 19, 2007 ReVised manuscript receiVed September 18, 2007 Accepted September 18, 2007 IE070838W
