- Email: [email protected]

Quasi-steady-state modeling of dendritic growth Pavel Levin 1 Materials Engineering Department, Ben Gurion University of the Negev, PO Box 653, 84105 Beer-Sheva, Israel Received 13 February 2002; received in revised form 13 February 2002; accepted 26 February 2003 Communicated by A.R. Bishop

Abstract A new method of the moving-boundary problem analysis is developed. After proposed coordinate transformation, the quasisteady-state differential equation was reduced to equivalent Laplace equation. Transformed boundary conditions of a rescaled field distribution reflect the presence of rate-dependent sources on interface. Taking into account the local rate-dependence of the Neumann’s boundary conditions, non-equilibrium pattern formation was considered. The problem of dendrite fractal growth was reduced to one of interaction of conformal to tips charged particles. Conditions of the quasi-steady growth and the recurrence formula for k-order dendrite spacing were derived. Theoretically obtained scaling laws for interface shape, dendrite spacing, critical sidebranch distance are confirmed by available experimental data. 2003 Elsevier Science B.V. All rights reserved. PACS: 05.70.Fh; 68.70.+w; 81.10.Aj; 81.30.Fb Keywords: Solid-liquid interface; Crystal growth; Dendrite; Solute; Quasi-steady model

1. Introduction The dynamics of growing interfaces with spontaneous formation of non-equilibrium patterns has been an active field of research in the last decade. This moving-boundary problem is relevant to a broad class of applications in crystal growth, fluid displacement in porous media, dielectric breakdown, flame propagation, etc. [1,2]. As a rule, the analysis is based on well-developed methods of 2D Laplace equation solution, not considering non-uniformity of the ratedependent boundary conditions along interface and simply omitting the gradient-containing term. AccordE-mail address: [email protected] (P. Levin). 1 On leave.

ing to derived scenario of pattern growth [2], the surface tension selects out of continuum of solutions the trajectory, leading as a result of pattern competition in the weak-noise linear regime, to a single-pattern fixed point. At that the equal-pattern solution is only a saddle point on the trajectory. Such scenario simply follows from potential energy minimization principle applied to considered Laplacian system, but obviously contradicts with the real pattern behavior. Generally saying, the Laplacian approach is not applicable in such interpretation. Neglecting the local ratedependence of the Neumann’s boundary conditions, one actually considers the problem of non-equilibrium pattern formation as steady one. At that the multipattern stable attractor solution is inconsistent with the potential energy minimization principle.

0375-9601/03/$ – see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(03)00371-2

384

P. Levin / Physics Letters A 310 (2003) 383–388

The classical stability analysis of dendritic growth during solidification for initially smooth interface was carried out by Mullins and Sekerka [3,4]. In fact, this analysis was made on the base of perturbed one-dimensional solution for quasi-steady temperature distribution (zero coordinate separation constant). For considering the stable dendrite growth after initial perturbations, the 3D problem should be solved. A mathematical analysis for the single dendrite tip growth was developed by Ivantsov [5] and extended for the case of non-Fickian diffusion [6]. This analysis is based on approximation of the shape of a needle and of isotherm surfaces by paraboloids, neglecting the Gibbs–Thomson effect. Most theoretical studies use this approach, but experiments show that the tip shape deviates considerably from paraboloid [7]. The mentioned partial solution is not applicable for the study of the mutual influence of growing tips. Langer and Muller-Krumbhaar analyzed the stability of parabolic tips and found that the continuum of Ivantsov’s solutions is divided into a stable and an unstable region [8]. The marginal stability parameter, determining relationship between tip radius and tip velocity, was introduced. The universal character of this parameter was shown in works cited below. Recently, Brener with collaborators developed a microscopic solvability theory, which reduces the diffusion problem to 2D Laplacian, neglecting the local rate-dependence of the boundary conditions and introducing the anisotropy of the surface free energy as a prerequisite for the stable attractor solution [9–12]. For a pure melt, Bisang and Bilgram [7] found the experimental power-law fit to dendrite contour, close to derived by Brener shape corrections to Ivantsov’s paraboloid. However, these results are good just for neglectfully small anisotropy parameter and disagree with phase field model numerical simulations for considerable anisotropy values [13,14]. For binary systems, according to Trivedy and Mason [15], a number of well-characterized experiments on the measurements of tip radius and corresponding stability parameter show their weak dependence on anisotropy factor, that contradicts to predictions of the solvability theory. Particularly, for small characteristic radius scale, the dendrite contours with spherical tip region on photographs do not seem to be properly fitted by power-law line at all. In the Brener’s non-axisymmetric case, amplitude of sidebranches in-

creases exponentially as function of position coordinate and proportionally to thermal fluctuation strength, at that the coordinate of first sidebranch is a power function of anisotropy-dependent stability parameter. However, such anisotropy-dependence was not confirmed by mentioned experiments. Besides, sidebranch amplitude does not grow exponentially, as predicted. The secondary tip radius is of the same order as primary one, and considerable sidebranch formation starts at some critical distance from the primary tip. The microscopic solvability theory approach also cannot explain the existence of stable multi-dendrite solution. Probably, the surface free energy anisotropy and the thermal noise just serve for initial disturbing the equilibrium at a smooth surface, but do not determine the stable attractor solution for dendritic scale, including tip radius, dendrite spacing and sidebranch critical distance. The interface mobility is influenced by Gibbs free energy, determined by undercooling and solute redistribution [1,4,12,15,16]. It was found experimentally that the tip curvature radius Rt , dendrite primary and secondary spacing λ(1) and λ(2) are uniquely determined by interface rate and undercooling, according to certain scaling laws [15]. It was shown that the physics of processes which control dendritic structural scales in the binary system is analogous to the growth of pure thermal dendrites. The scaling laws, found experimentally, are considered theoretically in the present Letter.

2. Fundamental solution of the quasi-steady-state problem The temperature or concentration distribution is considered correspondingly for a pure melt or for the case of diffusion-limited solidification at a given temperature gradient ∇T . A mathematical statement for the temperature or concentration field in cylindrical coordinate system (x, r = (y − y0 )2 + (z − z0 )2 ), bound with solid/liquid interface, moving with constant rate vint along axe x, takes a form: ∂ 2 u 1 ∂u ∂ 2 u vint + 2+ + r ∂r D ∂r 2 ∂x ∂u + cos(ψint )δ(x − xint) = 0. × ∂x

(1)

P. Levin / Physics Letters A 310 (2003) 383–388

Here, u = (T − T∞ )c/L or u = (X − X∞ )/XLS is the rescaled temperature or concentration field, measured from the infinity values T∞ or X∞ . D is the diffusivity constant, c—specific heat, L—latent heat, XLS = XL − XS —solute drag at interface during solidification for given conditions (L and S refer to liquidus and solidus); xint and ψint are the interface coordinate and normal angle with axe x; δ(x) is Dirac’s function. In order to consider the rate-dependent terms in quasi-steady-state differential equation (1), the coordinate transformation was proposed [17]: vint (x − x0 ) , ξ − ξ0 = 1 − exp − D ρ − ρ0 ≡ (ψ − ψ0 )2 + (ζ − ζ0 )2 vint vint r exp − (x − x (2) = 0) . D D

385

where ξ0 is considered as a current value of the variable ξ . Υ(1) refers here to fundamental solution of the Laplace’s equation. The fundamental solution (5) does not give exactly zero for “transformed coordinate of infinity” ξ − ξ0 = 1. The Green’s function for the system (3), (4) can be obtained by the method of images, at that Υ(1) =

(2π)−1 (ξ − ξint )2 + (ρ − ρint )2

−

(2π)−1 (ξ − ξint − 2)2 + (ρ − ρint )2

.

(6)

However, solution (5) gives a good approximation for negative and small positive ξ , characterizing dendritic scale (Fig. 1). The partial solution, developed on

After coordinate transformation (2), one can obtain an equivalent Laplace’s equation: ∂ 2 u 1 ∂u ∂ 2 u cos(ψint )δ(ξ − ξint ) + = 0. + + ∂ρ 2 ρ ∂ρ ∂ξ 2 (1 − (ξ − ξ0 ))

(a)

(3)

The trivial boundary conditions (zero field at infinity) after transformation (2) take a form: u|(ξ −ξ0 )2 +(ρ−ρ0 )2 →∞, ξ −ξ0 <1 = u|ξ −ξ0 =1 = 0.

(4)

Generally saying, the “surface charge density” is proportional to the local interface rate. In rescaled coordinates, the problem is analogous to the field of unit-density surface charge. The equivalent Laplacian problem (3), (4) differs from traditional one by nontrivial boundary conditions: as by “transformed coordinate of infinity” ξ −ξ0 = 1 as well by denominator in the correspondent rate-dependent term. At that an effective charge depends on considered pole coordinate (ξ0 , ρ0 ). For the unit point source with coordinates (ξint , ρint ) in half-space, the fundamental solution of an elliptic differential equation of class (1) can be represented as the generalization of a Laplace’s equation fundamental solution [18]: u(1) =

(2π)−1

(ξ − ξint + (ρ − ρint Υ(1) , ≡ 1 − (ξint − ξ0 ) )2

)2

1 − (ξint − ξ0 )

−1

(5)

(b)

(c) Fig. 1. The schema of dendrite growth: in rescaled coordinates (a), (b), and in absolute coordinates (c).

386

P. Levin / Physics Letters A 310 (2003) 383–388

the base of Green’s functions (5), (6), contains also a gradient term [18]. For a binary system, the partial solution for reduced field, which consists of derived concentration term and given temperature gradient term, should be considered.

In absolute coordinates, the well-known instability conditions, characterized by stability parameter σ ∗ = 0.019 . . .0.027 [15], can be obtained. According to (10), the marginal stability value of σ ∗ corresponds to some critical retarding of a plane solidification front, characterized by partition coefficient, k 1.

3. Initial instability analysis 4. Scaling laws of dendritic growth For the plane isothermal interface, ξint = 0, the problem after transformation (2) comes to 1D Neumann’s problem: ∂ 2 u cos(ψint )δ(ξ − ξint ) = 0. + 1 − (ξ − ξ0 ) ∂ξ 2

(7)

The interface field value, uint = uLv < 0, corresponds to a liquidus temperature at infinity concentration and given solidification rate. Analogously to (5), the partial solution of problem (7), (4) with a moving plane solidifying front can be expressed as: u ≈ uLv + uv ξ − kξ/(1 + ξ ),

(8)

where, uv is reduced field gradient, k—partition coefficient at solidifying [4] (numerically equal to rescaled temperature undercooling relatively liquidus value, uLv ). Let us consider a local thermodynamical equilibrium in a liquid at a planar interface with fluctuations in the form of spherical tips with curvature radiuses Rt (1) = Rt (ρt (1) = ρt ), Fig. 1(a). The resultant reduced field gradient at small tip height ξt can be considered as (uv − k). The planar interface becomes to be unstable, if partition coefficient k exceeds numerically reduced field gradient uv : (k − uv )ρt − 2γ /ρt 0,

(9)

where γ = vint Γ c/(DL) or γ = vint Γ /(DTLS ) = vint Γ /(DmL XLS ) is capillary length, Γ —Gibbs– Thomson coefficient, mL —liquidus temperature-concentration slope. Denoting according to [8] the expression σ ∗ = (1 − uv /k)/2

The quasi-steady state is characterized by zero resultant field flux, or, for system (8), by equity of undercooling at plane interface and overheating at infinity relatively liquidus point: k = uv ,

(12)

For the pure melt, the scale length is much larger, at that possible undercooling may be sufficient only for small local temperature gradient. Practically, uniform undercooling is required [4,7]. For the binary system, according to experimental data of Trivedy and Mason [15], the tip radius, ρt , and, therefore, diffusion coefficient, D, do not depend on interface temperature. For explanation of this phenomena, the concentration field, generated by growing tips, should be considered. Condition of Laplacian field potential energy minimization determines the semi-spherical tip shape in rescaled coordinates (Fig. 1(b)). As follows from derived above stability condition, the quasi-steady growth requires zero resultant reduced field flux at an interface. A growth rate-dependent “charge” is numerically equal to area of a tip projection, perpendicular to a growth direction. After equating generated, πρt2 , and 2 u , field fluxes, the corresponding scalexternal, ψλ(1) v ing law for rescaled primary dendrite spacing, ψλ(1) , can be derived: ψλ(1) /ρt = π/uv . (13) Obtained scaling parameter for primary dendrite spacing is inversely proportional to square root of reduced field gradient, which, according to (2), (12), can be expressed as:

(10)

as stability parameter, one immediately obtains the instability conditions for planar interface: γ k > uv , ρt = (11) . kσ ∗

k = −uLv .

uv ≈

D∇T . 2vt TLS

(14)

These results are in a good agreement with experimental scaling law [15]. For experimental conditions D = 1.05 × 10−4 mm2 /s, TLS = 3.5 . . .4 K, ∇T =

P. Levin / Physics Letters A 310 (2003) 383–388

1.19 . . .1.57 K/mm and growth rates, vt , mm/s, 3 × 10−4 , 1 × 10−2 , 4 × 10−2 , one can obtain scaling parameters, ψλ(1) /ρt , correspondingly, 6.2, 44.8, 88.2. Those are very close to correspondent experimental values of scaling parameter, λ(1)/Rt , 10.8, 54.4, 71.4. Some lower power of vt in the experimental dependence can be explained by narrowing the crystallization range, TLS , with the rate increasing. Let us consider the self-similar k-order dendrite growth with a constant tip velocity vint = vt (k) , to which the coordinate system is ascribed (Fig. 1(c)). The coordinate ξ(k) is measured from the center of the tip radius, ρt (k) ≈ Rt (k) vt (k) /D 1. The reduced field, generated by a semi-spherical tip, is determined by Eqs. (5), (6) for 3D space, and by projected tip area. Thus, the problem can be reduced to interaction of corresponding to tips charged particles, situated in their centers. One can consider effect of n = 9 tips, surrounding first tip (i = 1) with spacing ψλ(k) and corresponding distance between opposite higher-order tips, ψλ(k) = ψλ(k−1) − 2ρt (k) . The generated reduced field at the point (ξ(k) = ρt (k) , ρ(k) = 0) can be written as the sum of tip fields ut (k) : n

u(k)i

i=1

= ρt (k)

√ ρt (k) 2 1 + 1+ 2 ψλ(k) 2 2 −0.5 λ(k−1) ± 1+ λ(k)

λ(k−1) 2 −0.5 + 2+ . λ(k)

ing should be realized, that follows from condition of maximally uniform “surface charge” distribution. Therefore, the scaling law for secondary tip spacing can be represented in a full correspondence to experimental values [15]: ψλ(k) /ρt (k) = 2.4 . . .4.8.

For k = 1, “minus” is substituted, and ψλ(0) = 2. Whereas the term 1/2 in the sum in Eq. (15) reveals to the self-influence of dendrite i = 1, the effect of surrounding dendrites for case λ(k) λ(k−1) is approximately 1.71(ρt (k)/ψλ(k) ). While the selfinfluence of the tip growth is greater than total effect of others, the equal-tips solution remains to be stable. Therefore, the condition of growth stability can be ascribed as ρt (k) /ψλ(k) 0.29. Taking into account sequential retarding as well as angular deviation during secondary dendrite growth, it can be accepted for k > 1: ρt (k) /ψλ(k) 0.21 . . .0.41. Practically, during the tip arising, the minimal permitted instability spac-

(16)

Dendrite growth takes place in direction of external field gradient. At that a channel width between needles in transformed coordinates, ψλ(k−1) , remains constant. In the absolute coordinates, that corresponds to a needle inclination of order (ψλ(k)/2 − ρ(ξ )). As far as can be judged by photographs for considered above experiments [15], the real angles are in a good correlation with theoretical ones for different vt (scale lengths D/v, µm, 140, 22.8, 10.5, 2.6). The field at the needle surface is determined by 3D field of tip “charge”, ut (1) = ρt2 /(4ξ ), and by external reduced gradient, uv . For a quasi-steady problem, the mentioned fields are fixed relatively tip (see Fig. 1(c)). According to (11), the secondary tips should arise with the same radius order as primary ones: Rt (2) ≈ Rt . At that undercooling, sufficient for secondary tip formation, can be reached just at distance ξs from the primary tip, which corresponds to resultant rescaled field value uLv = −uv . Therefore, the first sidebranch coordinate is determined by ratio: ξs /ρt = (4uv )−0.5 .

(15)

387

(17)

Let us consider the case of single dendrite growth from a pure melt [7], with undercooling 25 . . . 142 mK (uLv = (5 . . . 27) × 10−4 ). The experimental values xs /Rt = 23 . . . 12 are in a good correspondence with theoretical range ξs /ρt = 22.7 . . .9.5. As can be seen from photographs for considered above binary system [15] (uv = (4 . . . 314) × 10−4 ), the experimental values are xs /Rt = 17 . . .3.5. That is in a reasonable correspondence with theoretical range ξs /ρt = 24.8 . . .2.8.

5. Conclusions The present solution of the 3D moving-boundary problem is consistent with the deterministic character of parameter selection in a fractal dendritic growth, at that the surface free energy anisotropy parameter should not be introduced. By considering the local

388

P. Levin / Physics Letters A 310 (2003) 383–388

rate-dependence of the Neumann’s boundary conditions in the equivalent Laplacian system and zero resultant field flux condition, one can derive the main relationships for the quasi-steady fractal growth. Obtained scaling laws for interface shape, dendrite spacing, critical sidebranch distance are in a good correspondence with available experimental data for the pure melt and for the binary systems.

Acknowledgement The author thanks Dr. A. Bershadskii for his help. References [1] W. van Saarlos, Phys. Rep. 301 (1998) 9. [2] J. Casademunt, F.X. Magdaleno, Phys. Rep. 337 (2000) 1.

[3] W.W. Mullins, R.F. Sekerka, J. Appl. Phys. 35 (1964) 444. [4] W. Kurz, D.J. Fisher, Fundamentals of Solidification, Trans. Tech. Publications, Switzerland, 1986. [5] G.P. Ivantsov, Dokl. Akad. Nauk SSSR 58 (1947) 567. [6] P.K. Galenko, D.A. Danilov, Phys. Lett. A 272 (2000) 207. [7] U. Bisang, J.H. Bilgram, Phys. Rev. E 54 (1996) 5309. [8] J.S. Langer, H. Muller-Krumbhaar, Acta Metall. 26 (1978) 1681. [9] E. Brener, Phys. Rev. Lett. 71 (1993) 3653. [10] E. Brener, D. Temkin, Phys. Rev. E 51 (1995) 351. [11] M. Ben Amar, E. Brener, Physica D 98 (1996) 128. [12] E. Brener, H. Muller-Krumbhaar, D. Temkin, Phys. Rev. E 54 (1996) 2714. [13] J.A. Warren, W.J. Boetinger, Acta Metall. Mater. 43 (1995) 689. [14] A. Karma, W.-J. Rappel, Phys. Rev. E 57 (1998) 4323. [15] R. Trivedi, J.T. Mason, Metall. Trans. A 22 (1991) 235. [16] M.J. Aziz, W.J. Boetinger, Acta Metall. Mater. 42 (1994) 527. [17] P. Levin, Heat Mass Transfer, submitted for publication. [18] S.L. Sobolev, Partial Differential Equations of Mathematical Physics, Pergamon Press, London, 1964, 430 p.