Phase field modeling of convective and morphological instability during directional solidification of an alloy

(1)

Journal of Crystal Growth 295 (2006) 202–208

Phase ﬁeld modeling of convective and morphological instability during

directional solidiﬁcation of an alloy

C.W. Lan

, M.H. Lee, M.H. Chuang, C.J. Shih

Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan, ROC

Received 16 March 2006; received in revised form 22 July 2006; accepted 26 July 2006 Communicated by G.B. McFadden

Abstract

Phase field modeling is carried out to investigate the convective and morphological instability during directional solidification of a succinonitrile/acetone alloy. Considering the presence of gravity, we have found that the planar interface could become wrinkled even beyond the Mullins–Sekerka instability; this is originated from the lateral solute segregation induced by the flow. For the cases slightly above the onset of instability, morphologies of shallow cells are affected by the convection as well. The cells with different wavelengths and depths can coexist due to the flow-induced segregation. The coupling of long- (convective mode) and short wavelengths (morphological mode) is illustrated for the first time, which cannot be predicted by the linear stability theory. As the growth rate is further increased, the effect of the buoyancy decreases.

PACS: 68.70.+w; 81.30.Fb

Keywords: A1. Adaptive phase ﬁeld simulation; A1. Convection; A1. Morphological instability; A1. Solidiﬁcation

1. Introduction

During alloy directional solidification, as the constitu-tional supercooling occurs to a certain degree, cellular or dendritic patterns develop as a result of the interplay of solute, temperature, and flow fields with the solidification interface. Without convection, the planar interface be-comes unstable at a certain critical value of the control parameter, which is either the pulling speed (V) or the temperature gradient (G), as described by the classical Mullins–Sekerka (MS) theory [1]. The morphological instability may interact with the convection, and this was firstly realized by Coriell et al.[2]by using a linear stability analysis in the presence of thermosolutal convection in the melt. Two modes of instability were identified, which were distinguishable by the spatial wavelength of the interface. The morphological mode was found consistent with the prediction of the MS theory, while the long-wavelength

branch, i.e., the convective mode, was not affected by the interface. In between, some interaction region leading to an oscillatory state was found. Coriell and McFadden [3]

further performed linear stability analyses for a tin containing lead alloy, for both stabilized thermal and solutal fields. They also showed that the minima of the neutral stability curves occurred at much longer wave-lengths, sometimes on the order of ampoule diameter, than the minima of the purely morphological stability. In an interesting directional solidification experiment by Schaefer and Coriell[4]using succinonitrile containing ethanol in an ampoule, it was illustrated that the interface deformation, especially the pit formation at the center of the interface, was caused by thermal convection. The pit formation was further revealed through computer simulation by Lan and Tu[5], and the onset of the morphological instability was significantly earlier than that predicted by MS theory. Fully nonlinear analyses by using a finite element method with symmetry boundaries were also investigated by Mehrabi in Brown’s group [6]. Given the width (half of the wavelength) of the computational domain, the www.elsevier.com/locate/jcrysgro

Corresponding author. Tel./fax: +886 2 2363 3917. E-mail address:[email protected] (C.W. Lan).

(2)

interaction of the ﬂow and the morphological evolution was illustrated for the convective mode. Unfortunately, their simulation was not able to simulate the appearance of the both modes due to the limited mesh resolution and deformation. In fact, to cover a large domain accommo-dating the long (convective mode) and short (morphology mode) wavelengths and to allow complex morphological evolutions for a solidiﬁcation simulation is not trivial, especially for the front-tracking approach.

The phase field model has emerged as a powerful tool for simulating solidification problems with a complex mor-phological evolution (e.g., Refs. [7–9]). However, quanti-tative simulation of alloy solidification is not trivial. An extremely small interface thickness (on the order of 0.1 mm) is required for thermodynamic consistency, while a large domain (41000 mm) is needed to accommodate the solute boundary layer and characteristic wavelengths. In addi-tion, solute trapping[10]can be induced due to the diffuse interface, especially at high solidification speed. Fortu-nately, this problem was first resolved by Karma[11]using the anti-trapping model and by Shih et al. [12] using a simple interface model. Furthermore, with the use of the adaptive grid [13], a quantitative simulation for alloy solidification has become feasible.

In this report, the adaptive phase field method[13]using the Karma’s anti-trapping model is adopted to simulate the directional solidification of a succinonitrile (SCN)/acetone alloy under an unstable solutal field near the onset of the MS instability. The interaction of the convection and the cell development is illustrated for the first time; in fact, they are indistinguishable as pointed out by Mehrabi[6]. In the next section, the phase field model used in the simulation is briefly described. Section 3 is devoted for results and discussion, followed by a short conclusion.

2. Phase ﬁeld model

The computational domain investigated in this study is sketched inFig. 1, and a two-dimensional case is adopted for simplicity. The thermal proﬁle is assumed to be linear having a gradient G as a working approximation. As the solidiﬁcation starts, the light solute (acetone) is rejected from the solid and accumulates near the melt/solid interface to reduce the density there. As discussed in Ref.[2], for thermosolutal convection the stability criterion is more complicated than simply examining the sign of the net density gradient. The mechanism is instead based on the different rates of diffusion of heat and solute (‘‘double-diffusion’’).

Therefore, in the presence of gravity, this destabilizes the system, as illustrated inFig. 2c. Away from the interface, the system is stabilized by the thermal effect;Figs. 2a and b

describe schematically the corresponding destabilizing concentration and stabilizing thermal proﬁles, respectively. Therefore, to a certain degree of acetone accumulation, convection is induced due to hydrodynamic instability. The convection leads to lateral segregation and interacts with

the development of the local interface morphology. In order to simulate the physical problem, an adaptive phase ﬁeld model is used.

Similar to the volume-averaged derivation by Becker-mann et al. [14], the two-dimensional phase ﬁeld equation for an alloy can be obtained as the following:

Z2qf qt ¼GmA r Z 2_rf q qx ZZ 0qf qy þ q qy ZZ 0qf qx fð1 fÞð1 2fÞ d2 m_AðTA_mT mj jcLÞ 5f2ð1 fÞ2 d , ð1Þ

where the anisotropic function Z is defined for the four-fold symmetry, i.e., Z ¼ 1 þ cos 4Y while e is the intensity of the anisotropy and Y ¼ tan1_{½ðqf=qyÞ=ðqf=qxÞ is the} angle between the direction normal of the interface and the x (horizontal) axis. The phase field variable f(x, y, t) changes rapidly from 0 (solid) to 1 (liquid) at the interface. G is the Gibbs–Thomson coefficient, mA the kinetic

coefﬁcient, and m the slope of liquidus temperature line (Ti); TAm is the melting point of pure A (SCN) and cL the

concentration of the solute (acetone) at the liquid phase side. Furthermore, d is the thickness of the diffuse interface. The corresponding energy (for temperature T) and species (for the acetone concentration C) equations are

ARTICLE IN PRESS

Fig. 1. Schematic of the computational domain for directional solidiﬁca-tion under gravity.

(3)

also obtained as follows: qT qt þ r ðpðfÞVT Þ ¼ ar 2_TLA Cp qf qt, (2) qC qt þ r pðfÞC pðfÞ þ kð1 pðfÞÞV ¼ r ~D rC ð1 kÞp 0_ðfÞC pðfÞ þ k½1 pðfÞrf þ 1 ~ Djatc , ð3Þ where ~ D ¼ DSþ ðDLDSÞpðfÞ pðfÞ þ kð1 pðfÞÞ. Also, jatc¼dð1 kÞ C pðfÞ þ kð1 pðfÞÞ _qf qt rf rf

is the anti-trapping current proposed by Karma [11] to suppress the effect of solute-trapping. DS and DL are the

solute diffusivities in the solid and liquid phases, respec-tively and p(f) is an interpolation function to describe solid/melt mixture properties in the diffuse interface; p(f) ¼ f is used here. In addition, k is the equilibrium segregation coefﬁcient, a the thermal diffusivity, LA the

latent heat, and Cp the speciﬁc heat. The continuity and

momentum equations for the pressure P and the velocity vector V can further be derived as:

r ðpðfÞr_LVÞ ¼ 0, (4)

qðpðfÞr_LVÞ

qt ¼ r ðpðfÞrLVVÞ þ r ðmLrpðfÞVÞ rP þ r_LB hmLð1 pðfÞÞfð1 fÞ

d2 V, ð5Þ

where rL is the melt density, mL the viscosity, and h the

friction coefﬁcient for the solid/melt interaction in the diffuse interface; h ¼ 2.757 [14]. Since the Boussinesq approximation is adopted, the linear form of the body

force B can be approximated by

B ¼ ðbCðC CfÞ þbTðT TfÞÞg, (6)

where bC and bT are the solutal and thermal expansion

coefﬁcients, respectively, and g the gravitational accelera-tion. In this study, the gravitational direction is pointing downward, i.e., g ¼ ey. The subscript f corresponds to the

reference state and is deﬁned at the inlet condition, i.e., at y ¼ H. Also, to neglect kinetic effects, we choose the kinetic coefﬁcient mAbased on thin-interface analysis[15]:

1 mk

¼Z½5a2=4½d=DL½ðDLLA=aCpÞ þ mj jð1 kÞC, (7) where a2¼47/75.

The above equations with a given initial condition can then be solved by an adaptive ﬁnite volume method [16], where the minimum grid size (Dx) is smaller than the interface thickness d; d ¼ 0.25–1 mm is used in this study. In this study, a linear temperature proﬁle having a gradient G is assumed. Therefore, the thermal equation, Eq. (2), is not solved.

Before the results and discussion, the phase field simulation without flow is first validated through the comparison with the one by Echebarria et al.[15] for the directional solidification of a SCN alloy; all the physical and simulation parameters can be found in Ref. [15] as well. The calculated steady-state cell shapes at two different interface thicknesses, and the comparison with those obtained by Echebarria et al. [15] are shown in Fig. 3a. The morphology predicted by the well-known Saffman– Taylor shape [17] is also included for comparison. As shown, good agreement is found; the Saffaman–Taylor shape is more accurate near the tip. Further convergence tests on the steady-state tip concentration and tip radius as a function of the interface thickness are shown inFigs. 3b and c, respectively. InFig. 3b, the solidus concentration is along the cell center; the effective equilibrium coefficient keffis also shown. InFig. 3c, the tip radii are determined by

Fig. 2. Schematic illustration of the axial distribution of (a) concentration, (b) temperature, and (c) buoyancy driving forces causing the convective instability during directional solidiﬁcation; the destabilizing region exists only near the interface.

(4)

using fourth-order polynomial ﬁtting within 1 mm domain around the cell tip. Again, these properties converge nicely as d is reduced to 0.5 mm. The results obtained by using the

simple interface model[12]are also in good agreement with the results in Fig. 3. However, since Karma’s model has been widely used, for comparison purpose we shall report the simulated results based on their formulation.

3. Results and discussion

To investigate the effect of convection, we take the directional solidiﬁcation of SCN containing 0.1 mol% acetone as an example (k ¼ 0.1, m ¼ 222 K/mol frac.), while the physical properties are taken from SCN directly. The parameters used in simulation are summarized in

Table 1. A linear temperature proﬁle having a gradient G of 100 K/cm is used for simulation. In other words, the temperature is assumed to be linear and ﬁxed, i.e., the frozen temperature approximation. In reality, this assump-tion is reasonable because thermal diffusivity is much larger than the solutal diffusivity. At this concentration, at steady state, based on the MS theory, the critical pulling velocity Vc and wavelength lc are about 5.7 mm/s and

160 mm, respectively. This critical velocity and wavelength have been used to examine our simulation, and the agreement is satisfactory (Vcless than 1% and lcless than

5%). To investigate the effect of convection induced by the gravity, we have chosen the pulling speeds lying around the critical velocity for study.

Without gravity, when the pulling velocity (Vp) is set at

5 mm/s, which is lower than critical one, a planar interface is obtained and the axial concentration proﬁle is in excellent agreement with the one-dimensional analytic solution. The calculated concentration ﬁeld is shown in

Fig. 4a. The straight streamlines are due to the material movement with respect to the ﬁxed frame. Once the gravity is introduced, because the acetone is lighter than SCN (bC40), a ﬂow cell is induced (convective instability), as

shown in Fig. 4b. Because of the domain size, the

ARTICLE IN PRESS

0.7 0.6 0.5 0.4 0.3 0.2 0 0.5 1.0 1.5 δ (μm) δ (μm) Theoretical value k G-T solution keff (b) (c) (a) Cs / C0 Cs / C 0 and keff 5.0 4.4 3.8 3.2 0 0.5 1.0 1.5 Saffman-Taylor shape Tip r adius ( μ m) 0 -0.5 -1 -0.4 -0.2 0 0.2 0.4 (x-xtip) / λlength (y-y tip ) / λlength δ = 0.25 μm (present) δ = 0.50 μm (present) δ = 0.17 μm (Echebarria et. al. [15]) δ = 0.34 μm (Echebarria et. al. [15]) Saffman-Taylor shape

Fig. 3. Comparison of the simulated results for the alloy directional solidiﬁcation investigated in Ref.[13]: (a) cell shape; (b) the convergence behavior of calculated tip concentration and the effective segregation coefﬁcient with the interface thickness; and (c) the convergence behavior of the calculated tip radii with the interface thickness.

Table 1

Parameters used in simulation

Physical properties of SCN–acetone alloy[6]

Mass diffusivity (liquid side), DL 109(m2/s)

Mass diffusivity (solid side), DS 1014(m2/s)

Kinetic viscosity, n 106_(m2_/s)

Gravity acceleration, g 10 (m/s2₎

Solutal expansion coefﬁcient, bC 1.8 103(mol %1)

Thermal expansion coefﬁcient, bT 8 10 4

(1/K) Melting temperature of pure melt, Tm

A

331.24 (K) Gibbs–Thompson coefﬁcient, G 6.48 108(K/m)

Anisotropy, e 0.007

Equilibrium partition coefﬁcient, k 0.1

Slope, |m| 2.22 (K/mol%)

Simulation parameters

Domain size 500 2500 (mm2)

Initial interface position, yi0 500 (mm)

Temperature gradient, G 10000 (K/m) Pulling velocity, Vp 5.0, 6.5, 12 (mm/s)

Reference length, l 106_(m)

Bulk concentration, C0 0.10 (mol %)

(5)

wavelength of the flow is two times of the domain width, i.e., 1000 mm. The lateral acetone segregation induced by the flow also leads to the large deformation of the interface; again, the wavelength of the deformation is 1000 mm. Interestingly, if we take a close-up look of the interface near the edge in Fig. 4b, a shallow cell (morphological instability) is observed. In other words, the convection due to the gravity, not only causes the long-wavelength deformation (convective mode) as a result of the convective instability, but also the short-wavelength one (morpholo-gical mode) due to the morpholo(morpholo-gical instability. The later is caused by the lateral acetone segregation. The higher acetone concentration makes the morphological instability happen earlier, even though the pulling speed is less than the critical one based on the MS theory. The liquid concentration profile at the interface (f ¼ 0.5) along the x-axis is further plotted inFig. 4c. It is clear that the lateral concentration gradients are caused by the convection. The dashed-line indicating C0/k is the average concentration,

which is also the same as the one without gravity because a steady state is reached. In addition, the cell formation also results in a local acetone accumulation. Furthermore, the local interface concentration near the morphological instability seen in Fig. 4c is in good agreement with the constitutional supercooling criterion[18],

C0 k Vp4 GDL jmjð1 kÞ. (8)

In Fig. 4c, our CLVp for obtaining the morphological

instability is found to be in good agreement with the above equation. In fact, the stability function [19] due to the interfacial energy for SCN is near unity. Hence, Eq. (8) gives a good criterion for the morphological breakdown.

If the pulling speed is increased to 6.5 mm/s (4Vc), the

calculated steady-state interface is no longer planar even without gravity. A shallow-cell structure appears as shown inFig. 5ain the absence of gravity. The wavelength of the shallow cells is about 2lc/3, i.e., about 100 mm, which is

Fig. 4. Effect of gravity for Vp¼ 5 mm/s (oVc) on the ﬂow and concentration ﬁelds and the interface shape: (a) without gravity; (b) with gravity; (c) the

(6)

consistent with a bifurcated branch near the onset in Ref. [20]. With buoyancy convection, the streamlines are slightly distorted; if the pulling is subtracted from the streamline values, a counterclockwise flow cell appears. Similar to the one inFig. 4b, the flow induces the lateral concentration segregation. As a result, the shallow cells disappear on the left of the domain due to the reduced acetone concentration. Meanwhile, two cells remain and become much deeper in the right corner, again as a result of higher acetone concentration there. Some pitch-off solute droplets are also observed in the bottom of the cell glooves. As a result of the physical solute trapping, the lateral solute gradients for the buoyancy flow are slightly reduced. And this could be a reason that the overall flow intensity in the present case seems to be smaller than that inFig. 4. On the other hand, due to the straightened effect of the higher pulling speed, the flow cell is not obvious.

Fig. 5cpresents the acetone concentration obtained from

Figs. 5a and balong the x-axis. Because the pulling speed is larger, the critical concentration is smaller indicating that

the region having morphological instability is wider, as compared to that in Fig. 4c. Again, the critical CLVp for

having the morphological instability is found in good agreement with Eq. (8).

Furthermore, if the pulling speed is further increased, the solutal boundary layer is thinner, on the order of DL/Vp,

and the convection in the lateral direction could be suppressed due to the smaller destabilized regime, as shown inFig. 6at Vp¼12 mm/s. The destabilizing driving

force, the Rayleigh number, is proportional to (DL/Vp)3

for solutal convection. Therefore, as the pulling speed increases, the convective speed due to the gravity decreases. Beside the increase in the cell depth, the effect of gravity becomes insigniﬁcant. The wavelength is about 77 mm, which is about lc/2. In the present simulation, the buoyant

melt speed is on the order of the crystal growth speed (pulling speed).

Finally, as we plot the MS boundary in Fig. 7, the physical feature of the previous simulations becomes clear. Since the length in the x-direction is ﬁxed at 500 mm, the

ARTICLE IN PRESS

Fig. 5. Effect of gravity for Vp¼ 6.5 mm/s (4Vc) on the ﬂow and concentration ﬁelds and the interface shape: (a) without gravity; (b) with gravity; (c) the

concentration distributions along the x-axis obtained from (a) and (b) at the interface.

(7)

simulated ﬂow cell has a wavelength of 1000 mm, and the induced interface deformation has the same wavelength as well. This wavelength caused by the buoyancy effect is much longer than that caused by the constitutional supercooling (about lc/2). In other words, the coexistence

of the long- and short-wavelength modes is found in the

present investigation, which is not seen from the previous linear stability prediction and the non-linear ﬁnite element simulations.

4. Conclusions

We have successfully used the phase field model to simulate the thermal-solutal effect on the interface mor-phology during directional solidification of an SCN/ acetone alloy. By using the phase field simulation, the deformation of the interface and the species/momentum transports can be treated easily. In the presence of gravity, the solidification interface can become unstable even at a stable regime predicted by Mullins–Sekerka theory due to the lateral solutal segregation induced by the hydrody-namic instability. The shallow cells near the onset of the instability are affected by the flows as well. The cells with different wavelengths and depths can coexist due to the flow-induced segregation. The coupling and co-existence of the long- (convective) and short- (morphological) wave-length modes are found and are illustrated for the first time by the phase field simulation.

Acknowledgment

This research is sponsored by the National Science Council of Taiwan.

References

[1] W.W. Mullins, R.F. Sekerka, J. Appl. Phys. 35 (1964) 444. [2] S.R. Coriell, M.R. Cordes, W.J. Boettinger, R.F. Sekerka, J. Crystal

Growth 49 (1980) 13.

[3] S.R. Coriell, G.B. McFadden, J. Crystal Growth 94 (1989) 513. [4] R.J. Schaefer, S.R. Coriell, Metall. Trans. A 15 (1984) 2109. [5] C.W. Lan, C.Y. Tu, J. Crystal Growth 220 (2000) 619.

[6] M. R. Mehrabi, Ph.D. Thesis, Massachusetts Institute of Technol-ogy, 1994.

[7] J.A. Warren, C. Beckermann, A. Karma, Annu. Rev. Mater. Res. 32 (2002) 163.

[8] A. Karma, W.J. Rappel, Phys. Rev. E 53 (1996) 3017. [9] C.W. Lan, Y.C. Chang, J. Crystal Growth 250 (2003) 525. [10] N.A. Ahmad, A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys.

Rev. E 58 (1998) 3436.

[11] A. Karma, Phys. Rev. Lett. 87 (2001) 115701.

[12] C.J. Shih, M.H. Lee, C.W. Lan, J. Crystal Growth 282 (2005) 515. [13] C.W. Lan, C.J. Shih, M.H. Lee, Acta Mater. 53 (8) (2005) 2285. [14] C. Beckermann, H.-J. Diepers, I. Steinbach, A. Karma, X. Tong,

J. Comput. Phys. 154 (1999) 468.

[15] B. Echebarria, R. Folch, A. Karma, M. Plapp, Phys. Rev. E 70 (2004) 061604.

[16] C.W. Lan, C.C. Liu, C.M. Hsu, J. Comput. Phys. 178 (2002) 464. [17] P.G. Saffman, G.I. Taylor, Proc. Roy. Soc. London A 254 (1958)

312.

[18] W.A. Tiller, K.A. Jackson, J.W. Rutter, B. Chalmers, Acta Matall. 1 (1953) 428.

[19] R.F. Sekerka, in: W. Bardsley, D.T.J. Hurle, J.B. Mullin (Eds.), Crystal Growth: An Introduction, Elsevier, North Holland (Chapter 15), 1973.

[20] K. Tsiveriotis, R.A. Brown, Phys. Rev. B 48 (1993) 13495. Fig. 6. Calculated ﬂow and solute ﬁelds, as well as the interface shape, for

Vp¼ 12 mm/s (bVc).

Fig. 7. Calculated wavelengths and the comparison with the MS boundary; the cross symbols indicate the convective mode, circle symbols the morphological mode, and the square symbols the morphological mode under buoyancy convection.