Steady needle growth with 3-D anisotropic surface tension

Abstract

The effect of the anisotropic interfacial energy on dendritic growth has been an important subject, and has preoccupied many researchers in the field of materials science and condensed matter physics. The present paper is dedicated to the study of the effect of full 3-D anisotropic surface tension on the steady state solution of dendritic growth. We obtain the analytical form of the first order approximation solution in the regular asymptotic expansion around the Ivantsov’s needle growth solution, which extends the steady needle growth solution of the system with isotropic surface tension obtained by Xu and Yu (J. J. Xu and D. S. Yu, J. Cryst. Growth, 1998, 187: 314; J. J. Xu, Interfacial Wave Theory of Pattern Formation: Selection of Dendrite Growth and Viscous Fingering in a Hele-Shaw Flow, Berlin: Springer-Verlag, 1997).

The solution is expanded in the general Laguerre series in any finite region around the needle-tip, and it is also expanded in a power series in the far field behind the tip. Both solutions are then numerically matched in the intermediate region. Based on this global valid solution, the dependence of Peclet number Pe and the interface’s morphology on the anisotropy parameter of surface tension as well as other physical parameters involved are determined. On the basis of this global valid solution, we explore the effect of the anisotropy parameter on the Peclet number of growth, as well as the morphology of the interface.

Graphical abstract

Keywords

dendrite growth / steady state solution / Laguerre series / asymptotic expansion / anisotropy of surface energy

1 Introduction

In studying free dendritic growth, the effect of anisotropic interfacial energy on pattern formation and selection bas been an important subject, which has preoccupied many researchers in the field of materials science and condensed matter physics for quite a long time now [1–4]. During the past decades, most theoretical investigations on this subject were focused on the system of 2-D dendritic growth with anisotropy of surface tension, or the system of 3-D axi-symmetric dendritic growth with axial anisotropy of surface tension, neglecting the three dimensionality of the anisotropy of surface tension [5–15]. With these efforts, the effect of the anisotropy of surface tension on the morphology of the steady state, interfacial stability mechanisms and selection of dendritic pattern formation for 2-D systems are now well explored analytically.

A realistic system of dendritic growth is always three dimensional. The anisotropy of surface tension at the interface is also three dimensional. Up to now in literature, there are very few analytical works on dendritic growth in systems with full 3-D anisotropic surface tension. As a consequence, the effects of three dimensionality of anisotropy on the steady state solution as well as the stability mechanisms of dendritic growth are still unclear.

In order to explore the effect of 3-D anisotropic surface tension on the stability mechanisms and the selection condition of dendritic growth, as the first step, it is necessary to study the effect of 3-D anisotropic surface tension on the steady state of dendritic growth. The goal of our effort is to attack this problem.

We obtain the analytical solution in the form of regular asymptotic expansion around Ivantsov’s needle growth solution, which extends the work on steady needle growth solution of the system with isotropic surface tension conducted by Xu and Yu [7, 14]. The solution is expended in the general Laguerre series in any finite region around the needle-tip and is also expanded in a power series in the far field behind the tip. Both solutions are then matched in the intermediate region numerically. On the basis of this global valid solution, we explore the effect of the anisotropy parameter on the Peclet number of growth, as well as the morphology of the interface.

The paper is arranged as follows. In Section 2, we formulate the problem in the moving paraboloidal coordinate system; in Section 3, we derive the asymptotic solution in the finite region behind the tip in the form of a general Laguerre series; in Section 4, we derive the solution in the region far away from the tip in the power series; then in Section 5, we match two solutions in the intermediate region numerically; in Section 6, we show the numerical results for various cases; finally in Section 7, we give the discussions and conclusions.

2 Mathematical formulation of the problem of steady needle crystal growth

Let us consider a system of free dendritic growth from a pure substance at the later stage of evolution, when the dendrite is fully developed with a sufficiently long stem. At this stage, the effects of the initial growth conditions as well as the situation at the root in the finite region around the tip on the solution will diminish to a minimum. Thus, for a mathematical model of the steady state of dendritic growth at the later stage, one may consider a needle-like crystal growing with a constant tip-velocity U into an undercooled pure melt with the undercooling temperature T_∞ < T_M0, where T_M0 is the melting temperature of a flat interface (refer to Refs.[7, 14]). The full 3-D anisotropic surface tension at the interface between liquid and solid phase is taken into account, so that the needle-like crystal under consideration is non-axisymmetrical. For simplicity, we assume that in the growth system, the mass density ρ, the specific heat c_p, and other thermal characteristic constants of the solid state are the same as that of the liquid state (the so-called symmetric model); gravity is negligible, so no convection is involved.

As usual, we use the growth speed of dendrite U, as the scale of the velocity, and the thermal length {\mathcal{l}}_{\text{T}}={\kappa }_{\text{T}}/U as the length scale, {\mathcal{l}}_{\text{T}}/U as the time scale and ΔH/(c_pρ) as the temperature scale, where ΔH is the latent heat per unit volume of the solid.

We adopt the paraboloidal coordinate system moving together with the growing need-crystal (refer to Ref.[14]), which can be defined through the Cartesian coordinate system (x, y, z) by (see Fig. 1):

\begin{array}{}\text{(2.1)}& \frac{x}{{\eta }_0^2}=\xi \eta \cos \theta \mathrm{,}\frac{y}{{\eta }_0^2}=\xi \eta \sin \theta \mathrm{,}\frac{z}{{\eta }_0^2}=\frac{1}{2}({\xi }^2-{\eta }^2)\end{array}

where the constant {\eta }_0^2 is to be chosen, so that the steady interface shape satisfies

\begin{array}{}\text{(2.2)}& {\eta }_{\text{s}}(0)=1\end{array}

It is known that this constant is just the Peclet number, Pe₀, of the system with zero surface tension, fully determined by the undercooling T_∞ via the Ivantsov solution. In general, the Peclet number, Pe, is defined as the ratio of the tip radius {\mathcal{l}}_{\text{t}} and the thermal diffusion length {\mathcal{l}}_{\text{T}}, and is a function of all physical parameters involved in the system, including the anisotropy parameter.

The governing equation for three-dimensional dendritic growth can be written in the forms:

(1) In the liquid region:

\begin{array}{}\text{(2.3)}& \begin{array}{ll}\hfill & \frac{{\partial }^{2}T}{\partial {\xi }^2}+\frac{{\partial }^{2}T}{\partial {\eta }^2}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)\frac{\partial T}{\partial \xi }\hfill \\ \hfill & +\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)\frac{\partial T}{\partial \eta }+\left(\frac{1}{{\xi }^2}+\frac{1}{{\eta }^2}\right)\frac{{\partial }^{2}T}{\partial {\theta }^2}=0\hfill \end{array}\end{array}

(2) In the solid region:

\begin{array}{}\text{(2.4)}& \begin{array}{ll}\hfill & \frac{{\partial }^2{T}_{\text{S}}}{\partial {\xi }^2}+\frac{{\partial }^2{T}_{\text{S}}}{\partial {\eta }^2}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)\frac{\partial T_{\text{S}}}{\partial \xi }\hfill \\ \hfill & +\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)\frac{\partial T_{\text{S}}}{\partial \eta }+\left(\frac{1}{{\xi }^2}+\frac{1}{{\eta }^2}\right)\frac{{\partial }^2{T}_{\text{S}}}{\partial {\theta }^2}=0\hfill \end{array}\end{array}

The boundary conditions are:

(1) The up-stream far-field condition:

\begin{array}{}\text{(2.5)}& T=T_{\infty }\text{as}\eta \to \infty \end{array}

(2) The regularity condition:

\begin{array}{}\text{(2.6)}& T_{\text{S}}=O(1)\text{as}\eta \to 0\end{array}

(3) The interface conditions, at η = η_s(ξ, θ), we have

(a) the temperature continuity condition:

\begin{array}{}\text{(2.7)}& T_L=T_{\text{S}}\end{array}

(b) the Gibbs–Thomson condition: The surface tension coefficient can be expressed in the form:

\begin{array}{}\text{(2.8)}& \gamma ={\gamma }_0[1+{\widehat{\alpha }}_4{Q}_4(\mathit{n})]\end{array}

where γ₀ is the coefficient of isotropic surface tension, α₄ is the coefficient of forth fold anisotropic surface tension,n is the unit normal vector of the interface pointing to the liquid phase, and Q_4(\mathit{n})=(n_x^4+n_y^4+n_z^4) is the anisotropic scalar function. The temperature at the interface can be determined by the Herring formula [10, 11]:

\begin{array}{}\text{(2.9)}& T_{\text{S}}=T_{\text{I}}=-{\epsilon }^2\left[1+{\alpha }_4\left(\frac{\partial Q_4}{\partial {\theta }_1^2}+\frac{\partial Q_4}{\partial {\theta }_2^2}\right)\right]{\mathcal{K}}_3\{{\eta }_{\text{s}}(\xi \mathrm{,}\theta )\}\end{array}

where (θ₁, θ₂) respectively represent the orientation angles of the directions of the two principal curvatures; {\mathcal{K}}_3\{{\eta }_{\text{s}}(\xi \mathrm{,}\theta )\} is the local twice mean curvature of the interface η_s(ξ, θ), which is derived as:

\begin{array}{}\text{(2.10)}& \begin{array}{l}{\mathcal{K}}_3\{{\eta }_{\text{s}}(\xi \mathrm{,}\theta )\}\\ =-\frac{1}{{({\xi }^2+{\eta }_{\text{s}}^2)}^{\frac{3}{2}}{\left[{\xi }^2{\eta }_{\text{s}}^2(1+{\eta }_{\text{s}}^{\prime 2})+{\eta }_{\text{s}\mathrm{,}\theta }^2({\xi }^2+{\eta }_{\text{s}}^2)\right]}^{\frac{3}{2}}}\\ \cdot \left\{\right[(1+{\eta }_{\text{s}}^{\prime 2})({\xi }^2+{\eta }_{\text{s}}^2)][-{(\xi {\eta }_{\text{s}})}^2({\eta }_{\text{s}}{\eta }_{\text{s}}^{\prime }-\xi )\\ +2\xi {\eta }_{\text{s}\mathrm{,}\theta }^2({\xi }^2+{\eta }_{\text{s}}^2)+{\xi }^2{\eta }_{\text{s}}{\eta }_{\text{s}\mathrm{,}\theta \theta }({\eta }_{\text{s}}{\eta }_{\text{s}}^{\prime }-\xi )\\ -\xi {\eta }_{\text{s}}{(\xi {\eta }_{\text{s}})}^{\prime }({\eta }_{\text{s}\mathrm{,}\theta }^2+{\eta }_{\text{s}}{\eta }_{\text{s}\mathrm{,}\theta \theta })]-2[{\eta }_{\text{s}\mathrm{,}\theta }{\eta }_{\text{s}}^{\prime }({\xi }^2+{\eta }_{\text{s}}^2)]\\ \cdot [{(\xi {\eta }_{\text{s}})}^{\prime }{\eta }_{\text{s}\mathrm{,}\theta }({\xi }^2+{\eta }_{\text{s}}^2)+\xi {\eta }_{\text{s}}({\eta }_{\text{s}}{\eta }_{\text{s}}^{\prime }-\xi )(\xi {\eta }_{\text{s}\mathrm{,}\theta }{)}^{\prime }\\ -\xi {\eta }_{\text{s}}{({\eta }_{\text{s}}{\eta }_{\text{s}\mathrm{,}\theta })}^{\prime }{(\xi {\eta }_{\text{s}})}^{\prime }]+[{(\xi {\eta }_{\text{s}})}^2+{\eta }_{\text{s}\mathrm{,}\theta }^2({\xi }^2+{\eta }_{\text{s}}^2)]\\ \cdot \left[{(\xi {\eta }_{\text{s}})}^{″}\xi {\eta }_{\text{s}}({\eta }_{\text{s}}{\eta }_{\text{s}}^{\prime }-\xi )+\xi {\eta }_{\text{s}}{(\xi {\eta }_{\text{s}})}^{\prime }(1-{\eta }_{\text{s}}^{\prime 2}-{\eta }_{\text{s}}{\eta }_{\text{s}}^{″})\right]\}\end{array}\end{array}

Especially, for the special axi-symmetric case, η_s = η_s(ξ), we have

\begin{array}{}\text{(2.11)}& \begin{array}{ll}{\mathcal{K}}_3\{{\eta }_{\text{s}}(\xi \mathrm{,}\theta )\}=\hfill & -\frac{1}{\sqrt{{\xi }^2+{\eta }_{\text{s}}^2}}\left[\frac{{\eta }_{\text{s}}^{″}}{{(1+{\eta }_{\text{s}}^{\prime 2})}^{\frac{3}{2}}}-\frac{1}{{\eta }_{\text{s}}{(1+{\eta }_{\text{s}}^{\prime 2})}^{\frac{1}{2}}}\right]\hfill \\ \hfill & +\frac{{\eta }_{\text{s}}^{\prime }({\eta }_{\text{s}}^2+2{\xi }^2)-\xi {\eta }_{\text{s}}}{\xi ({\xi }^2+{\eta }_{\text{s}}^2)(1+{\eta }_{\text{s}}^{\prime 2}{)}^{\frac{1}{2}}}\hfill \end{array}\end{array}

The complete form of Eq. (2.9) in the paraboloidal coordinate system is very complicated. Fortunately, in this investigation, we only need its first order approximation form, which is given in the Appendix. In Eq. (2.9), the interfacial stability parameter ε is defined as:

\begin{array}{}\text{(2.12)}& \epsilon =\frac{\sqrt{{\mathcal{l}}_{\text{c}}/{\mathcal{l}}_{\text{T}}}}{{\eta }_0^2}\mathrm{,}{\mathcal{l}}_{\text{c}}=\frac{{\gamma }_0{c}_{\text{p}}\rho T_{\text{M0}}}{{(\text{Δ}H)}^2}\end{array}

where T_M0 is the melting temperature of the flat interface.

(c) the heat balance condition:

\begin{array}{}\text{(2.13)}& \begin{array}{ll}\hfill & \left(\frac{\partial }{\partial \eta }-{\eta }_{\text{s}}^{\prime }\frac{\partial }{\partial \xi }\right)\left(T-T_{\text{S}}\right)\hfill \\ \hfill & -{\eta }_{\text{s}\mathrm{,}\theta }(\frac{1}{{\xi }^2}+\frac{1}{{\eta }_{\text{s}}^2})\frac{\partial }{\partial \theta }(T-T_{\text{S}})+{\eta }_0^2(\xi {\eta }_{\text{s}}{)}^{\prime }=0\hfill \end{array}\end{array}

(4) The tip-smoothness condition: at ξ = 0,

\begin{array}{}\text{(2.14)}& {\eta }_{\text{s}}(\mathrm{0,}\theta )=\mathrm{1,}{\eta }_{\text{s}}^{\prime }(\mathrm{0,}\theta )=0\end{array}

The problem under study is how to solve the three unknown functions involved: the temperature fields T (ξ, η, θ), T_S(ξ, η, θ), and the interface shape η_s(ξ, θ).

3 Steady state asymptotic expansion solution away from the root region

It is known that the interfacial stability parameter ε is, in practice, very small, whose numerical magnitude is about 0.1–0.2. Thus, in order to examine the steady perturbation induced by the anisotropic surface tension in the region away from the root of dendrite, one may make the regular perturbation expansion (RPE) around the Ivantsov solution in the limit ε→0. Namely, let

\begin{array}{}\text{(3.1)}& \{\begin{array}{l}T=T(\xi \mathrm{,}\eta \mathrm{,}\theta )=T_0(\eta )+{\epsilon }^2{\eta }_0^2{T}_1(\xi \mathrm{,}\eta \mathrm{,}\theta )+\cdots \hfill \\ T_{\text{S}}=T_{\text{S}}(\xi \mathrm{,}\eta \mathrm{,}\theta )=T_{\text{S}0}+{\epsilon }^2{\eta }_0^2{T}_{\text{S}1}(\xi \mathrm{,}\eta \mathrm{,}\theta )+\cdots \hfill \\ {\eta }_{\text{s}}={\eta }_{\text{B}}(\xi \mathrm{,}\theta )=1+{\epsilon }^2{\eta }_1(\xi \mathrm{,}\theta )+\cdots \hfill \end{array}\end{array}

Substitute (3.1) into the system (2.3)–(2.13) and equate coefficients of like powers of ε to zero, then we can successively derive the systems of each order approximation.

3.1 Zeroth-order approximation O({\epsilon }^0)

The zeroth order approximation solution is the Ivantsov’s solution,

\begin{array}{}\text{(3.2)}& \{\begin{array}{l}{T}_0=T_{\infty }+\frac{{\eta }_0^2}{2}{\text{e}}^{\frac{{\eta }_0^2}{2}}{E}_1\left(\frac{{\eta }_0^2{\eta }^2}{2}\right)\hfill \\ T_{\text{S}0}=0\hfill \\ {\eta }_0=1\hfill \\ T_{\infty }=-\frac{{\eta }_0^2}{2}{\text{e}}^{\frac{{\eta }_0^2}{2}}{E}_1\left(\frac{{\eta }_0^2}{2}\right)\hfill \end{array}\end{array}

where E₁(x) is the exponential integral function defined as (see Ref. [13]):

\begin{array}{}\text{(3.3)}& E_1(x)={\displaystyle {\int }_x^{\infty }}\frac{{\text{e}}^{-t}}{t}\text{d}t\end{array}

3.2 First order approximation O({\epsilon }^2)

In the first-order approximation, we derive that in the liquid phase region

\begin{array}{}\text{(3.4)}& \begin{array}{ll}L\left\{T_1\right\}=\hfill & [\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)\frac{\partial }{\partial \xi }\hfill \\ \hfill & +\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)\frac{\partial }{\partial \eta }]T_1+(\frac{1}{{\xi }^2}+\frac{1}{{\eta }^2})\frac{{\partial }^2{T}_1}{\partial {\theta }^2}=0\hfill \end{array}\end{array}

while in the solid phase region

\begin{array}{}\text{(3.5)}& \begin{array}{ll}L\left\{T_{\text{S1}}\right\}=\hfill & [\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)\frac{\partial }{\partial \xi }\hfill \\ \hfill & +\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)\frac{\partial }{\partial \eta }]T_{\text{S1}}+(\frac{1}{{\xi }^2}+\frac{1}{{\eta }^2})\frac{{\partial }^2{T}_{\text{S1}}}{\partial {\theta }^2}=0\hfill \end{array}\end{array}

The boundary conditions are

(1) As η → ∞,

\begin{array}{}\text{(3.6)}& T_1\to 0\text{(exponentially)}\end{array}

(2) It is hard to determine the asymptotic behavior of the solution at the first space, which specifies whether the solution is growing or decaying as ξ → ∞. However, one may very reasonably assume that the solution cannot grow exponentially, if it grows. Thus, we impose that as ξ → ∞,

\begin{array}{}\text{(3.7)}& \begin{array}{l}{T}_1\to 0\text{(or,}\text{at}\text{most,}\text{algebraically}\text{grows)}\hfill \end{array}\end{array}

(3) As η → 0,

\begin{array}{}\text{(3.8)}& T_{\text{S}1}=O(1)\mathrm{,}\text{is regular}\end{array}

(4) The interface conditions are linearized around the interface of Ivantsov paraboloid (η = 1), which can be written as:

(a) Temperature continuity condition:

\begin{array}{}\text{(3.9)}& T_1=T_{\text{S}1}+{\eta }_1\end{array}

(b) Gibbs-Thomson condition: The first order approximation of the interface temperature can be derived in the form (Ref. [11], the details of the derivation is given in the Appendix.):

\begin{array}{}\text{(3.10)}& T_{\text{S}1}=-[K_0+{\widehat{\alpha }}_4(F_0+F_4\cos 4\theta )]\end{array}

where

\begin{array}{}\text{(3.11)}& \begin{array}{ll}{K}_0\hfill & =\frac{2+{\xi }^2}{{(1+{\xi }^2)}^{\frac{3}{2}}}\mathrm{,}{F}_0=\frac{\frac{3}{4}{\xi }^6+\frac{3}{2}{\xi }^4+21{\xi }^2-6}{{(1+{\xi }^2)}^{\frac{7}{2}}}\hfill \\ F_4\hfill & =\frac{-\frac{15}{4}{\xi }^6-\frac{15}{2}{\xi }^4}{{(1+{\xi }^2)}^{\frac{7}{2}}}\hfill \end{array}\end{array}

(c) Heat balance condition:

\begin{array}{}\text{(3.12)}& \begin{array}{l}\frac{\partial }{\partial \eta }(T_1-T_{\text{S}1})+(2+{\eta }_0^2){\eta }_1+\xi \frac{\text{d}{\eta }_1}{\text{d}\xi }=0\hfill \end{array}\end{array}

(5) The tip-regularity condition, at ξ = 0,

\begin{array}{}\text{(3.13)}& {\eta }_1^{\prime }(0)=\mathrm{0,}{\eta }_1(0)=0\end{array}

Note that Eq. (3.10) is the only inhomogeneous condition in the above system. With the superposition principle, one may separate the solution vector q_1\equiv \left\{\begin{array}{c}{T}_1\mathrm{,}{T}_{\text{S}1}\mathrm{,}{\eta }_{\text{s}1}\end{array}\right\} into three parts:

\begin{array}{}\text{(3.13)}& q_1=q_1^{(\text{iso})}+{\widehat{\alpha }}_4{q}_1^{(0)}+{\widehat{\alpha }}_4{q}_1^{(4)}\end{array}

Among them, q_1^{(\text{iso})} is the isotropic component of the solution with {\widehat{\alpha }}_4=0, {\widehat{\alpha }}_4{q}_1^{(0)} is the axi-symmetric component of the solution induced by the anisotropy of surface tension, while {\widehat{\alpha }}_4{q}_1^{(4)} is the non-axi-symmetric component of the solution induced by the anisotropy of surface tension. These components of the solution are subject to the following simplified Gibbs-Thomson conditions at the interface, respectively:

\begin{array}{}\text{(3.14)}& \begin{array}{ll}\hfill & T_{\text{S1}}^{(\text{iso})}=-K_0(\xi )\mathrm{,}{T}_{\text{S1}}^{(0)}=-F_0(\xi )\hfill \\ \hfill & T_{\text{S1}}^{(4)}=-F_4(\xi )\cos 4\theta \hfill \end{array}\end{array}

while other boundary conditions remain unchanged.

The isotropic component of solution q_1^{(\text{iso})} has been obtained by Xu and Yu previously (refer to Ref. [14]). Therefore, the present paper attempts to solve the solutions {\widehat{\alpha }}_4{q}_1^{(0)} and {\widehat{\alpha }}_4{q}_1^{(4)} only.

3.3 The solution for the temperature field in liquid phase

The solution for Eq. (3.4) can be derived by using separation of variables. Let

\begin{array}{}\text{(3.15)}& T_1(\xi \mathrm{,}\eta \mathrm{,}\theta )=X(\xi )Y(\eta )\cos m\theta \end{array}

From Eq. (3.4) it follows that

\begin{array}{}\text{(3.16)}& X^{″}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)X^{\prime }+\left({\eta }_0^2{\lambda }_1^2-\frac{m^2}{{\xi }^2}\right)X=0\end{array}

\begin{array}{}\text{(3.17)}& \begin{array}{l}{Y}^{″}+\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)Y^{\prime }-\left({\eta }_0^2{\lambda }_1^2+\frac{m^2}{{\eta }^2}\right)Y=0\hfill \end{array}\end{array}

By letting

\begin{array}{}\text{(3.18)}& \sigma =\frac{{\eta }_0^2{\xi }^2}{2}\mathrm{,}X=X(\sigma )\end{array}

the equation (3.16) is transformed into the Kummer equation [13]:

\begin{array}{}\text{(3.19)}& \begin{array}{ll}\sigma X^{″}(\sigma )+(1-\sigma )X^{\prime }(\sigma )+\left(\frac{{\lambda }_1^2}{2}-\frac{m^2/4}{\sigma }\right)X=\hfill & 0\hfill \end{array}\end{array}

Furthermore, letting

\begin{array}{}\text{(3.19)}& X(\sigma )={\text{e}}^{\frac{\sigma }{2}}{\sigma }^{-\frac{1}{2}}W(\sigma )\end{array}

Eq. (3.19) is transformed in the Wittaker’s equation:

\begin{array}{}\text{(3.20)}& \frac{{\text{d}}^{2}W}{\text{d}{\sigma }^2}+\left(-\frac{1}{4}+\frac{\kappa }{\sigma }+\frac{1-4{\mu }^2}{4{\sigma }^2}\right)W=0\end{array}

with

\begin{array}{}\text{(3.20)}& \kappa =\frac{1+{\lambda }_1^2}{2}\mathrm{,}\mu =\frac{m}{2}\end{array}

which has the two fundamental solutions:

\begin{array}{}\text{(3.21)}& W(\sigma )={\text{e}}^{-\frac{\sigma }{2}}{\sigma }^{\frac{1}{2}+\frac{m}{2}}\{\begin{array}{c}M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\\ U\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\end{array}\end{array}

Here M and U are the confluent hypergeometric functions [13].

Hence, the fundamental solutions of Eq. (3.19) are

\begin{array}{}\text{(3.22)}& X(\sigma )={\text{e}}^{\frac{\sigma }{2}}{\sigma }^{-\frac{1}{2}}W(\sigma )={\sigma }^{\frac{m}{2}}\{\begin{array}{l}M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\hfill \\ U\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\hfill \end{array}\end{array}

Note that function M and U have the following properties:

• \begin{array}{l}M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\hfill \end{array} is regular at σ = 0;

• \begin{array}{l}U\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\hfill \end{array} has a logarithmic singularity for m = 0 and O(σ^−m/2) singularity for m > 0 at σ = 0. Thus, one must choose

\begin{array}{}\text{(3.23)}& X(\sigma )={\sigma }^{\frac{m}{2}}M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\end{array}

On the other hand, as σ → ∞, \begin{array}{l}M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,1}+m\mathrm{,}\sigma \right)\hfill \end{array} grows algebraically only for

\begin{array}{}\text{(3.24)}& \frac{m-{\lambda }_1^2}{2}=-n=\mathrm{0,}-\mathrm{1,}-\mathrm{2,}\cdots \end{array}

For other values of m-{\lambda }_1^2/2, the function M\left(\frac{m-{\lambda }_1^2}{2}\mathrm{,}1+m\mathrm{,}\sigma \right) grows exponentially, as σ → ∞. From the far field condition of the liquid phase (3.7), one must set {\lambda }_1^2/2=n+m/2 so that we derive

\begin{array}{}\text{(3.25)}& X(\sigma )={\sigma }^{\frac{m}{2}}M(-n\mathrm{,1}+m\mathrm{,}\sigma )=a_0{\xi }^m{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\end{array}

Here, we have denoted a_0={\left(\frac{{\eta }_0^2}{2}\right)}^{\frac{m}{2}}\frac{n!\Gamma (m+1)}{\Gamma (m+n+1)} and applied the formula:

\begin{array}{}\text{(3.26)}& L_n^{(m)}(\sigma )=\frac{\Gamma (m+n+1)}{n!\Gamma (m+1)}M(-n\mathrm{,1}+m\mathrm{,}\sigma )\end{array}

where L_n^{(m)}(x) is the m-th order general Laguerre polynomial.

We now return to solve Eq. (3.17). Let

\begin{array}{}\text{(3.27)}& \tau =\frac{{\eta }_0^2{\eta }^2}{2}\mathrm{,}Y=Y(\tau )\end{array}

So that Eq. (3.17) is transformed to

\begin{array}{}\text{(3.28)}& \tau Y^{″}(\tau )+(1+\tau )Y^{\prime }(\tau )-\left(\frac{{\lambda }_1^2}{2}+\frac{m^2/4}{\tau }\right)Y=0\end{array}

With the new variable Z(τ) defined as

\begin{array}{}\text{(3.29)}& Y(\eta )={\text{e}}^{-\frac{\tau }{2}}{\tau }^{-\frac{1}{2}}Z(\tau )\end{array}

Eq. (3.28) can be transformed to the Whittaker equation:

\begin{array}{}\text{(3.30)}& Z^{″}(\tau )+\left(-\frac{1}{4}+\frac{\kappa }{\tau }+\frac{\frac{1}{4}-{\mu }^2}{{\tau }^2}\right)Z(\tau )=0\end{array}

where

\begin{array}{}\text{(3.30)}& \mu =\frac{m}{2}\mathrm{,}\kappa =-\frac{1+{\lambda }_1^2}{2}\end{array}

Hence, we can obtain the two fundamental solutions of Eq. (3.28):

\begin{array}{}\text{(3.31)}& Y(\tau )={\text{e}}^{-\tau }{\tau }^{m/2}\{\begin{array}{l}M(1+m+n\mathrm{,1}+m\mathrm{,}\tau )\hfill \\ U(1+m+n\mathrm{,1}+m\mathrm{,}\tau )\hfill \end{array}\end{array}

Note that as z → ∞,

\begin{array}{}\text{(3.31)}& \begin{array}{l}M(a\mathrm{,}b\mathrm{,}z)\sim \frac{\Gamma (b)}{\Gamma (a)}{\text{e}}^z{z}^{a-b}\mathrm{,}U(a\mathrm{,}b\mathrm{,}z)\sim z^{-a}\hfill \end{array}\end{array}

In the present case, a = 1 + m + n, b = 1 + m, a−b = n. It is seen that the two basic solutions of Y(η) have the following properties, respectively:

• Y(\eta )={\text{e}}^{-\tau }{\tau }^{m/2}M(n+m+\mathrm{1,1}+m\mathrm{,}\tau ) vanishes or grows algebraically as η → ∞;

• Y(\eta )={\text{e}}^{-\tau }{\tau }^{m/2}U(n+m+\mathrm{1,1}+m\mathrm{,}\tau ) vanishes exponentially as η → ∞.

From the boundary condition (3.6), one must choose

\begin{array}{}\text{(3.32)}& Y(\eta )={\text{e}}^{-\frac{{\eta }_0^2{\eta }^2}{2}}\frac{{\eta }_0^m{\eta }^m}{2^{m/2}}U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2{\eta }^2}{2}\right)\end{array}

Thus, the general solution for the temperature in the liquid region is obtained as:

\begin{array}{}\text{(3.33)}& \begin{array}{ll}{T}_1(\xi \mathrm{,}\eta \mathrm{,}\theta )=\hfill & \frac{{\eta }_0^{2m}{\xi }^m{\eta }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\beta }_n{L}_n\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \\ \hfill & \cdot \frac{{\text{e}}^{-\frac{{\eta }_0^2{\eta }^2}{2}}U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2{\eta }^2}{2}\right)}{{\text{e}}^{-\frac{{\eta }_0^2}{2}}U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)}\cos m\theta \hfill \end{array}\end{array}

where the coefficients β_n are to be determined later.

3.4 The solution for the temperature field in solid phase

The solution for Eq. (3.5) can be derived similarly. Let

\begin{array}{}\text{(3.34)}& T_{\text{S}1}(\xi \mathrm{,}\eta \mathrm{,}\theta )=X(\xi )Y(\eta )\cos m\theta \end{array}

The solution for X(ξ) is entirely the same as that obtained for the liquid phase, namely,

\begin{array}{}\text{(3.35)}& \begin{array}{l}X(\sigma )={\sigma }^{\frac{m}{2}}M(-n\mathrm{,1}+m\mathrm{,}\sigma )=a_0{\xi }^m{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \end{array}\end{array}

The function Y (η) is also subject to Eq. (3.17). However, in the solid phase region, we let

\begin{array}{}\text{(3.36)}& \widehat{\tau }=-\frac{{\eta }_0^2{\eta }^2}{2}\end{array}

so that Eq. (3.17) is transformed into

\begin{array}{}\text{(3.37)}& \widehat{\tau }{Y}^{″}(\widehat{\tau })+(1-\widehat{\tau })Y^{\prime }(\widehat{\tau })+\left(\frac{{\lambda }_1^2}{2}-\frac{m^2/4}{\widehat{\tau }}\right)Y=0\end{array}

The form of Eq. (3.37) is completely the same as Eq.(3.19). Hence, its solution, regular at \widehat{\tau }=\eta =0, can be found as:

\begin{array}{}\text{(3.38)}& Y(\widehat{\tau })={\widehat{\tau }}^{\frac{m}{2}}M(-n\mathrm{,1}+m\mathrm{,}\widehat{\tau })=b_0{\eta }^m{L}_n^{(m)}\left(\frac{-{\eta }_0^2{\eta }^2}{2}\right)\end{array}

where b_0={\left(-\frac{{\eta }_0^2}{2}\right)}^{\frac{m}{2}}\frac{n!\Gamma (m+1)}{\Gamma (m+n+1)}. Thus, the general solution for the temperature in the solid phase is obtained as:

\begin{array}{}\text{(3.39)}& \begin{array}{l}{T}_{\text{S}1}(\xi \mathrm{,}\eta \mathrm{,}\theta )=\frac{{\eta }_0^{2m}{\xi }^m{\eta }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\alpha }_n\\ \begin{array}{ll}\hfill & \cdot \frac{L_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)L_n^{(m)}\left(-\frac{{\eta }_0^2{\eta }^2}{2}\right)}{L_n^{(m)}\left(-\frac{{\eta }_0^2}{2}\right)}\cos m\theta \hfill \end{array}\end{array}\end{array}

where the coefficients α_n are to be determined later.

3.5 The solution of the interface shape function

We expand the interface shape function η₁(ξ, θ) in the following general Laguerre series:

\begin{array}{}\text{(3.40)}& \begin{array}{ll}{\eta }_1(\xi \mathrm{,}\theta )\hfill & ={\eta }_1^{(\text{iso})}(\xi )+{\widehat{\alpha }}_m{\eta }_1^{(0)}(\xi )+{\widehat{\alpha }}_m{\eta }_1^{(m)}(\xi \mathrm{,}\theta )\hfill \\ \hfill & ={\displaystyle \sum _{n=0}^{\infty }}{\gamma }_n^{(\text{iso})}{L}_n\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)+{\widehat{\alpha }}_m{\displaystyle \sum _{n=0}^{\infty }}{\gamma }_n^{(0)}{L}_n\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \\ \hfill & +{\widehat{\alpha }}_m\frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\gamma }_n^{(m)}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\cos m\theta \hfill \end{array}\end{array}

Thus, the form of the solution derived above contains three sets of unknown constants: {α_n, β_n, γ_n} (n =0, 1, 2, ···), which can be determined by the three boundary conditions (3.9), (3.10) and (3.12).

Once the solution η₁(ξ, θ) is fully determined, one may determine the tip radius as well as the Peclet number with the effect of anisotropic surface tension. As a matter of fact, after linearization, Eq. (2.10) is reduced to

\begin{array}{}\text{(3.41)}& \begin{array}{l}{\mathcal{K}}_3\{{\eta }_{\text{s}}(\xi \mathrm{,}\theta \mathrm{,}\epsilon )\}=-\frac{2+{\xi }^2}{{(1+{\xi }^2)}^{\frac{3}{2}}}+\frac{{\epsilon }^2}{{(1+{\xi }^2)}^{\frac{1}{2}}}\\ \begin{array}{lll}\hfill & \hfill & \cdot [\frac{{\xi }^2(2+{\xi }^2)+4}{{(1+{\xi }^2)}^2}{\eta }_1+\frac{1+2{\xi }^2}{\xi (1+{\xi }^2)}{\eta }_1^{\prime }\hfill \\ \hfill & \hfill & +{\eta }_1^{″}+\frac{(1+{\xi }^2)}{{\xi }^2}{\eta }_{\mathrm{1,}\theta \theta }]+O({\epsilon }^4)\hfill \end{array}\end{array}\end{array}

Due to η₁(0) = η′(0) = 0, the twice of the mean curvature at the tip can be calculated as:

\begin{array}{}\text{(3.41)}& {\mathcal{K}}_3\{{\eta }_{\text{s}}(\mathrm{0,}\theta \mathrm{,}\epsilon )\}=-2+2{\epsilon }^2{\eta }_1^{″}(\mathrm{0,}\theta )+O({\epsilon }^4)\end{array}

By the definition of the Peclet number, we derive

\begin{array}{}\text{(3.42)}& \begin{array}{ll}\text{Pe}\hfill & =\frac{{\mathcal{l}}_{\text{t}}}{{\mathcal{l}}_{\text{T}}}=\frac{2{\eta }_0^2}{-{\mathcal{K}}_3(0)}+O({\epsilon }^4)\hfill \\ \hfill & =\frac{{\eta }_0^2}{1-{\epsilon }^2[{\eta }_1^{{(\text{iso})}^{″}}(0)+{\widehat{\alpha }}_4{\eta }_1^{{(0)}^{″}}(0)]}+O({\epsilon }^4)\hfill \end{array}\end{array}

It yields the dependence of the Peclet number of dendrite growth on the anisotropy of surface tension.

3.6 Determination of the unknown constants

By substituting the above obtained solutions (3.39) and (3.40) to the boundary conditions (3.9) and (3.10), one may derive the three sets of algebraic equations for the three sets of unknown constants. The details are given below.

3.6.1 Temperature continuity condition

From Eq. (3.9) it follows that

\begin{array}{}\text{(3.43)}& {\beta }_n={\alpha }_n+{\gamma }_n\end{array}

3.6.2 Gibbs-Thomson condition

From Eq. (3.10) it follows that

\begin{array}{}\text{(3.44)}& \begin{array}{ll}\hfill & \frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\alpha }_n^{(m)}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\cos m\theta \hfill \\ \hfill & =-F_m(\xi )\cos m\theta \hfill \end{array}\end{array}

It is derived that

• for the component with m = 0, we have

\begin{array}{}\text{(3.45)}& {\displaystyle \sum _{n=0}^{\infty }{\alpha }_n^{(0)}{L}_n^{(0)}}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)=-F_0(\xi )\end{array}

• for the component with m = 4, we have

\begin{array}{}\text{(3.46)}& \frac{{\eta }_0^8{\xi }^4}{2^4}{\displaystyle \sum _{n=0}^{\infty }{\alpha }_n^{(4)}{L}_n^{(4)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)=-F_4(\xi )}\end{array}

and

\begin{array}{}\text{(3.47)}& {\displaystyle \sum _{n=0}^{\infty }{\alpha }_n^{(\text{iso})}{L}_n^{(0)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)=-K_0(\xi )}\end{array}

Noting that the associated Laguerre polynomials L_n^{(m)}(x) (n = 0, 1, 2, ···) form the weighted orthogonal set of functions:

\begin{array}{}\text{(3.47)}& {\displaystyle {\int }_0^{\infty }}{x}^m{\text{e}}^{-x}{L}_n^{(m)}(x)L_k^{(m)}(x)=\{\begin{array}{lll}\mathrm{0,}\hfill & n\ne k\hfill & \hfill \\ \frac{(n+m)!}{n!}\mathrm{,}\hfill & n=k\hfill & \hfill \end{array}\end{array}

the coefficients α_n can be determined from the following integral:

\begin{array}{}\text{(3.48)}& {\alpha }_n^{(0)}=-{\displaystyle {\int }_0^{\infty }}{\text{e}}^{-x}{L}_n^{(0)}(x)F_0(\sqrt{2x}/{\eta }_0)\text{d}x\mathrm{,}m=0\end{array}

\begin{array}{}\text{(3.49)}& \begin{array}{ll}{\alpha }_n^{(4)}=\hfill & -\frac{\Gamma (n+1)}{\frac{{\eta }_0^4}{2^2}\Gamma (n+5)}\hfill \\ \hfill & \cdot {\displaystyle {\int }_0^{\infty }}{\text{e}}^{-x}{x}^2{L}_n^{(4)}(x)F_4(\sqrt{2x}/{\eta }_0)\text{d}x\mathrm{,}m=4\hfill \end{array}\end{array}

and

\begin{array}{}\text{(3.50)}& {\alpha }_n^{(\text{iso})}=-{\displaystyle {\int }_0^{\infty }}{\text{e}}^{-x}{L}_n^{(0)}(x)K_0(\sqrt{2x}/{\eta }_0)\text{d}x\end{array}

3.6.3 Heat conservation condition

We now turn to apply the boundary condition (3.12). From Eqs. (3.39) and (3.33), one obtains

\begin{array}{}\text{(3.51)}& \begin{array}{ll}\hfill & \frac{\partial T_{\text{S}1}}{\partial \eta }(\xi \mathrm{,1})\hfill \\ \hfill & =\frac{m{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\alpha }_n^{(m)}\cos m\theta \hfill \\ \hfill & +\frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}(-{\eta }_0^2)L_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \\ \hfill & \cdot \frac{n{L}_n^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)-(m+n)L_{n-1}^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}{-\frac{{\eta }_0^2}{2}{L}_n^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}{\alpha }_n^{(m)}\cos m\theta \hfill \\ \hfill & =\frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}\left[m+2n-2(m+n)\frac{L_{n-1}^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}{L_n^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}\right]\hfill \\ \hfill & \cdot L_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\alpha }_n^{(m)}\cos m\theta \hfill \end{array}\end{array}

and

\begin{array}{}\text{(3.51)}& \begin{array}{ll}\hfill & \frac{\partial T_1}{\partial \eta }(\xi \mathrm{,1})\hfill \\ \hfill & =\frac{m{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\beta }_n^{(m)}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\cos m\theta \hfill \\ \hfill & +\frac{{\eta }_0^{2m}{\xi }^m}{2^m}\{{\displaystyle \sum _{n=0}^{\infty }}{\beta }_n^{(m)}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \\ \hfill & \cdot \frac{-{\eta }_0^2{\text{e}}^{-\frac{{\eta }_0^2}{2}}U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)+2{\text{e}}^{-\frac{{\eta }_0^2}{2}}}{U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right){\text{e}}^{-\frac{{\eta }_0^2}{2}}}\hfill \\ \hfill & \cdot [(n+m+1)(1+n)U\left(n+m+\mathrm{2,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)\hfill \\ \hfill & -(n+m+1)U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)\left]\right\}\cos m\theta \hfill \\ \hfill & =\frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}[m-2(m+n+1)-{\eta }_0^2\hfill \\ \hfill & +2(n+m+1)(1+n)\frac{U\left(n+m+\mathrm{2,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)}{U\left(n+m+\mathrm{1,1}+m\mathrm{,}\frac{{\eta }_0^2}{2}\right)}]\hfill \\ \hfill & \cdot L_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\beta }_n^{(m)}\cos m\theta \hfill \end{array}\end{array}

In the above we have applied the formulas:

\begin{array}{}\text{(3.52)}& \begin{array}{ll}\hfill & a(1+a-b)U(a+\mathrm{1,}b\mathrm{,}x)=aU(a\mathrm{,}b\mathrm{,}x)+x{U}^{\prime }(a\mathrm{,}b\mathrm{,}x)\hfill \\ \hfill & x(L_n^{(m)}{)}^{\prime }(x)=n{L}_n^{(m)}(x)-(m+n)L_{n-1}^{(m)}(x)\hfill \end{array}\end{array}

With the notations

\begin{array}{}\text{(3.53)}& \{\begin{array}{ll}{a}_n\hfill & =(m+2n)-2(n+m)\frac{L_{n-1}^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}{L_n^{(m)}\left(\frac{-{\eta }_0^2}{2}\right)}\hfill \\ b_n\hfill & =-{\eta }_0^2-2n-m-2\hfill \\ \hfill & +2\left(n+m+1\right)\left(n+1\right)\frac{U\left(n+m+2,m+1,\frac{{\eta }_0^2}{2}\right)}{U\left(n+m+1,m+1,\frac{{\eta }_0^2}{2}\right)}\hfill \end{array}\end{array}

one may write

\begin{array}{}\text{(3.54)}& \begin{array}{l}\frac{\partial T_1}{\partial \eta }(\xi \mathrm{,1})=\frac{{\eta }_0^{2m}{\xi }^m}{2^m}\left[{\displaystyle \sum _{n=0}^{\infty }{b}_n{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)}\right]{\beta }_n\cos m\theta \\ \frac{\partial T_{\text{S}1}}{\partial \eta }(\xi \mathrm{,1})=\frac{{\eta }_0^{2m}{\xi }^m}{2^m}\left[{\displaystyle \sum _{n=0}^{\infty }{a}_n{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)}\right]{\alpha }_n\cos m\theta \end{array}\end{array}

Thus, from the boundary condition (3.12), we derive

\begin{array}{}\text{(3.55)}& \begin{array}{ll}\hfill & {\displaystyle \sum _{n=0}^{\infty }}{b}_n{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\beta }_n^{(m)}-{\displaystyle \sum _{n=0}^{\infty }}{a}_n{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\alpha }_n^{(m)}\hfill \\ \hfill & +(2+{\eta }_0^2){\displaystyle \sum _{n=0}^{\infty }}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right){\gamma }_n^{(m)}+{\displaystyle \sum _{n=0}^{\infty }}\{m{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\hfill \\ \hfill & +\left[2n{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)-2(n+m)L_{n-1}^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\right]\}{\gamma }_n^{(m)}\hfill \\ \hfill & =0\hfill \end{array}\end{array}

So that,

\begin{array}{}\text{(3.56)}& \begin{array}{l}{b}_n{\beta }_n^{(m)}-a_n{\alpha }_n^{(m)}+(2+{\eta }_0^2+2n+m){\gamma }_n^{(m)}\\ -2(n+m+1){\gamma }_{n+1}^{(m)}=\mathrm{0,}n=\mathrm{0,1,2,3,}\cdots \end{array}\end{array}

Combining Eqs. (3.43) and (3.56), we derive that

\begin{array}{}\text{(3.57)}& \begin{array}{l}{\gamma }_{n+1}^{(m)}=g_n{\gamma }_n^{(m)}+f_n{\alpha }_n^{(m)}\mathrm{,}n=\mathrm{0,1,2,}\cdots \hfill \end{array}\end{array}

where

\begin{array}{}\text{(3.58)}& g_n=(n+1)\frac{U\left(n+m+\mathrm{2,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}{U\left(n+m+\mathrm{1,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}\end{array}

\begin{array}{}\text{(3.59)}& \begin{array}{ll}{f}_n=\hfill & g_n+\frac{m+n}{m+n+1}\frac{L_{n-1}^{(m)}\left(-\frac{{\eta }_0^2}{2}\right)}{L_n^{(m)}\left(-\frac{{\eta }_0^2}{2}\right)}\hfill \\ \hfill & -2+\frac{m+1-\frac{{\eta }_0^2}{2}}{n+m+1}\hfill \end{array}\end{array}

Eq. (3.57) is a non-homogeneous difference equation for \left\{{\gamma }_n^{(m)}\right\}. For any given initial value {\gamma }_0^{(m)}, the formula (3.57) allows us to generate the sequence:

\begin{array}{}\text{(3.59)}& \{{\gamma }_0^{(m)}\mathrm{,}{\gamma }_1^{(m)}\mathrm{,}{\gamma }_2^{(m)}\mathrm{,}\cdots \mathrm{,}{\gamma }_n^{(m)}\mathrm{,}\cdots \}\end{array}

Thus, the function {\eta }_1^{(m)}({\gamma }_0^{(m)}\mathrm{,}\xi \mathrm{,}\theta ) can be written in the form:

\begin{array}{}\text{(3.60)}& \begin{array}{ll}\hfill & {\eta }_1^{(m)}({\gamma }_0^{(m)}\mathrm{,}\xi \mathrm{,}\theta )\hfill \\ \hfill & =\frac{{\eta }_0^{2m}{\xi }^m}{2^m}{\displaystyle \sum _{n=0}^{\infty }}{\gamma }_n^{(m)}{L}_n^{(m)}\left(\frac{{\eta }_0^2{\xi }^2}{2}\right)\cos m\theta \hfill \end{array}\end{array}

For the case m = 0, one may apply the tip condition (3.13), namely

\begin{array}{}\text{(3.61)}& {\eta }_1^{(0)}({\gamma }_0^{(0)}\mathrm{,0})={\displaystyle \sum _{n=0}^{\infty }{\gamma }_n^{(0)}=0}\end{array}

to determine {\gamma }_0^{(0)}. Once it is done, all the coefficients {\gamma }_n^{(0)}(n=\mathrm{0,1,2,}\cdots ) are determined. However, for the case m = 4, the tip conditions (3.14) are automatically satisfied by any constant {\gamma }_0^{(4)}. Hence, in order to determine the constant {\gamma }_0^{(4)}, one must use the far field condition for η₁(ξ, θ), that is η₁(ξ, θ) → 0, as ξ → ∞.

4 Asymptotic behavior of the needle-like crystal solution as ξ → ∞

The Laguerre polynomial expansion for η₁(ξ, θ) obtained above (3.40) and (3.60) is a convergent eigen-function enpension, valid in the whole interval . However, this expansion is an alternative series, and the numerical computation of this series at the point \xi \gg 1 cannot yield an accurate value of the sum function due to the numerical instability. In fact, the partial summation of the expansion with any large number of terms N, starts to rapidly oscillate when \xi \gg 1, as the partial sum becomes dominated by the past term {\gamma }_N{L}_N^{(m)}({\xi }^2/2) [see Fig. 2(a) and (b)]. In other words, the Laguerre polynomial expansion solution obtained above, numerically, is useful only in a finite region behind the tip (0⩽\xi ⩽{\xi }_{\text{max}}). To obtain the accurate numerical value of the solution in the region \xi \gg 1, one needs to find a different asymptotic form of the solution. For this purpose, we make the following asymptotic expansion for the solution as ξ → ∞:

\begin{array}{}\text{(4.1)}& \begin{array}{ll}\hfill & T_1(\xi \mathrm{,}\eta \mathrm{,}\theta )=\left[\frac{B_1^{(m)}(\eta )}{\xi }+\frac{B_2^{(m)}(\eta )}{{\xi }^2}+\cdots \right]\cos m\theta \hfill \\ \hfill & T_{\text{S}1}(\xi \mathrm{,}\eta \mathrm{,}\theta )=\left[\frac{A_1^{(m)}(\eta )}{\xi }+\frac{A_2^{(m)}(\eta )}{{\xi }^2}+\cdots \right]\cos m\theta \hfill \\ \hfill & {\eta }_1(\xi \mathrm{,}\eta \mathrm{,}\theta )=\left[\frac{{\widehat{C}}_1^{(m)}}{\xi }+\frac{{\widehat{C}}_2^{(m)}}{{\xi }^2}+\cdots \right]\cos m\theta \hfill \end{array}\end{array}

By substituting Eq. (4.1) into Eqs. (3.4) and (3.5), one may derive each order approximation of the solution successively.

4.1 Zeroth order approximation O(ξ⁻¹)

In the zeroth order approximation, we obtain the equations:

\begin{array}{}\text{(4.2)}& \left[\frac{{\text{d}}^2}{\text{d}{\eta }^2}+\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)\frac{\text{d}}{\text{d}\eta }-\left(-{\eta }_0^2+\frac{m^2}{{\eta }^2}\right)\right]\left(\begin{array}{c}{A}_1^{(m)}\\ B_1^{(m)}\end{array}\right)=0\end{array}

\begin{array}{}\text{(4.3)}& X^{″}+\left(\frac{1}{\xi }-{\eta }_0^2\xi \right)X^{\prime }+\left({\eta }_0^2{\lambda }_1^2-\frac{m^2}{{\xi }^2}\right)X=0\end{array}

\begin{array}{}\text{(4.4)}& \begin{array}{l}{Y}^{″}+\left(\frac{1}{\eta }+{\eta }_0^2\eta \right)Y^{\prime }-\left({\eta }_0^2{\lambda }_1^2+\frac{m^2}{{\eta }^2}\right)Y=0\hfill \end{array}\end{array}

with the interface conditions: at η = 1,

(1)

\begin{array}{}\text{(4.5)}& B_1^{(m)}(1)=A_1^{(m)}(1)+{\widehat{C}}_1^{(m)}\end{array}

(2)

\begin{array}{}\text{(4.6)}& \frac{A_1^{(m)}(1)}{\xi }+\frac{A_2^{(m)}(1)}{{\xi }^2}+\cdots =-F_m(\xi )\end{array}

where

\begin{array}{}\text{(4.6)}& F_0=\frac{\frac{3}{4}{\xi }^6+\frac{3}{2}{\xi }^4+21{\xi }^2-6}{{(1+{\xi }^2)}^{\frac{7}{2}}}\mathrm{,}{F}_4=\frac{-\frac{15}{4}{\xi }^6-\frac{15}{2}{\xi }^4}{{(1+{\xi }^2)}^{\frac{7}{2}}}\end{array}

The values of A_1^{(m)}(1) follow from the above, which are

\begin{array}{}\text{(4.6)}& \begin{array}{l}{A}_1^{(4)}(1)=\frac{15}{4}\mathrm{,}{A}_1^{(0)}(1)=-\frac{3}{4}\hfill \end{array}\end{array}

(3)

\begin{array}{}\text{(4.7)}& \frac{\partial }{\partial \eta }(B_1^{(m)}-A_1^{(m)})+(1+{\eta }_0^2){\widehat{C}}_1^{(m)}=0\end{array}

Eq. (4.2) is of similar type as that for Y (η) in Eq.(3.17). Letting

\begin{array}{}\text{(4.7)}& \tau =\frac{{\eta }_0^2{\eta }^2}{2}\text{and}F(\eta )=\frac{{\text{e}}^{-\frac{\tau }{2}}}{{\tau }^{\frac{1}{2}}}Z(\tau )\end{array}

where F(η) represents both A_1^{(m)}(\eta ) and B_1^{(m)}(\eta ), then Eq. (4.2) is transformed to Eq. (3.30) with {\lambda }_1^2=-1. Thus, we obtain the solutions

\begin{array}{}\text{(4.8)}& F(\eta )=\frac{{\eta }_0^m{\eta }^m}{2^{m/2}}{\text{e}}^{-\frac{{\eta }_0^2{\eta }^2}{2}}\{\begin{array}{c}M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2{\eta }^2\right)\\ U\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2{\eta }^2\right)\end{array}\end{array}

Taking into account the far field condition of the liquid phase region and the regular condition at the η = 0 in the solid phase region, we derive the solutions in the liquid region and solid region, respectively as follows:

\begin{array}{}\text{(4.9)}& A_1^{(m)}(\eta )={\widehat{A}}_1^{(m)}\frac{{\eta }_0^m{\eta }^m}{2^{m/2}}{\text{e}}^{-\frac{{\eta }_0^2{\eta }^2}{2}}M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2{\eta }^2\right)\end{array}

\begin{array}{}\text{(4.10)}& B_1^{(m)}(\eta )={\widehat{B}}_1^{(m)}\frac{{\eta }_0^m{\eta }^m}{2^{m/2}}{\text{e}}^{-\frac{{\eta }_0^2{\eta }^2}{2}}U\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2{\eta }^2\right)\end{array}

where {\widehat{A}}_1^{(m)} and {\widehat{B}}_1^{(m)} are constants to be determined. To determine the three sets of constants: \left\{{\widehat{A}}_1^{(m)}\mathrm{,}{\widehat{B}}_1^{(m)}\mathrm{,}{\widehat{C}}_1^{(m)}\right\}, we apply the interface conditions (4.5)–(4.7). From Eq. (4.6), we obtain

\begin{array}{}\text{(4.11)}& {\widehat{A}}_1^{(m)}=-\frac{{\text{e}}^{\frac{{\eta }_0^2}{2}}{2}^{m/2}}{{\eta }_0^{m}M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}{A}_1^{(m)}(1)\end{array}

From Eq. (4.5), we derive

\begin{array}{}\text{(4.12)}& {\widehat{B}}_1^{(m)}{P}_2={\widehat{A}}_1^{(m)}{Q}_2+{\widehat{C}}_1^{(m)}\end{array}

where we have defined

\begin{array}{}\text{(4.13)}& \{\begin{array}{c}{P}_2=\frac{{\eta }_0^m}{2^{m/2}}{\text{e}}^{-\frac{{\eta }_0^2}{2}}U\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)\\ Q_2=\frac{{\eta }_0^m}{2^{m/2}}{\text{e}}^{-\frac{{\eta }_0^2}{2}}M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)\end{array}\end{array}

From Eq. (4.7), we derive

\begin{array}{}\text{(4.14)}& {\widehat{B}}_1^{(m)}{P}_1-{\widehat{A}}_1^{(m)}{Q}_1+(1+{\eta }_0^2){\widehat{C}}_1^{(m)}=0\end{array}

where we have defined

\begin{array}{}\text{(4.15)}& \begin{array}{ll}{P}_1=\hfill & \frac{{\eta }_0^m{\text{e}}^{-\frac{{\eta }_0^2}{2}}}{2^{m/2}}[-U\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)\hfill \\ \hfill & -2U\left(\frac{m}{2}-\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)]\hfill \\ Q_1=\hfill & \frac{{\eta }_0^4{\text{e}}^{-\frac{{\eta }_0^2}{2}}}{2^2}[-M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)\hfill \\ \hfill & +(m+1)M\left(\frac{m}{2}-\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{1}{2}{\eta }_0^2\right)]\hfill \end{array}\end{array}

In terms of formulas (4.11)–(4.14), and making some organizations, one may derive

\begin{array}{}\text{(4.16)}& {\widehat{C}}_1^{(m)}=-\frac{2R+(m+1)S}{{\eta }_0^2-2R}{\widehat{A}}_1^{(m)}(1)\end{array}

where the following notations are used:

\begin{array}{}\text{(4.16)}& \begin{array}{l}R=\frac{U\left(\frac{m}{2}-\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}{U\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}\\ S=\frac{M\left(\frac{m}{2}-\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}{M\left(\frac{m}{2}+\frac{1}{2}\mathrm{,}m+\mathrm{1,}\frac{{\eta }_0^2}{2}\right)}\end{array}\end{array}

The constants {\widehat{B}}_1^{(m)} then follow from Eq. (4.14). Up to this point, the asymptotic solution in the far field is fully determined.

5 Matching the condition of the RPE with the far field asymptotic solution

Since the asymptotic solution (4.1) is applicable in the far field as ξ → ∞, it can thus yield the accurate numerical value of solution η₁(ξ_max, θ) at \xi ={\xi }_{\text{max}}\gg 1. On the other hand, as it has been demonstrated, the Laguerre expansion solution can be applied in the region (0⩽\xi ⩽{\xi }_{\max }) with good accuracy, provided ξ_max is not too large. Thus, in order to determine the parameter {\gamma }_0^{(m)}, we may let the Laguerre expansion solution shoot the value of the asymptotic solution (4.1) at ξ = ξ_max. In other words, we may choose the value of {\gamma }_0^{(m)} properly, so that the Laguerre expansion solution can meet the far field asymptotic solution at the intermediate point ξ = ξ_max. This procedure is shown in Fig. 2 (a) and (b) for a typical case of m = 0 and m = 4, respectively. The shape of a typical needle’s interface for the case {\eta }_0^2/2=0.5, {\widehat{\alpha }}_4=0.5, ε = 0.15 is shown in Fig. 3.

6 Numerical results

We have computed the solutions with various parameters of ε, {\eta }_0^2 and {\widehat{\alpha }}_4. Our numerical results for the isotropic case {\widehat{\alpha }}_4=0 recovered the results previously obtained by Xu and Yu [14]. Our numerical results are summarized as follows.

6.1 Variation of Peclet number with the parameters ε and {\widehat{\alpha }}_4

The Peclet number of needle crystal growth can be calculated with the formula (3.42). The numerical values of Pe for the case with various parameters ε and {\widehat{\alpha }}_4 are shown in Table 1. The numerical values of {\eta }_1^{″}(0) and {\mathcal{K}}_3(0) at the tip are computed through the paraboloid-like tip that best fits the solution for the shape of the interface obtained by using the method of least square approximation.

The variations of relative Peclet number \widehat{\text{Pe}}=\text{Pe}/{\text{Pe}}_0 with anisotropy parameter {\widehat{\alpha }}_4 for the typical case {\eta }_0^2/2 =0.5 and ε = 0, 0.1, 0.12, 0.15, 0.2, 0.3, 0.4 are shown in Fig. 4 (a) and (d). Here Pe₀ is the Peclet number for the case of ε = 0. It is seen that the results that we obtained for the case of {\widehat{\alpha }}_4=0 are consistent with those previously obtained by Xu and Yu [14], that is, the Peclet number decreases as ε increases. However, if {\widehat{\alpha }}_4 is not small, the Peclet number increases as ε increases. Particularly, when {\widehat{\alpha }}_4\approx 0.2, as ε increases from 0.2 to 0.4, the relative Peclet number \widehat{\text{Pe}} increases very rapidly.

The variations of the Peclet number Pe with T_∞ for the cases of anisotropic parameters {\widehat{\alpha }}_4 =−0.5,−0.2,−0.1, 0.1, 0.2, 0.5 are shown in Fig. 5. It is seen again that the anisotropy parameter {\widehat{\alpha }}_4 may significantly affect the Peclet number, especially in the large undercooling regime. When the {\widehat{\alpha }}_4 number is negative or positive but small, the Peclet number decreases with increasing ε. However, when the {\widehat{\alpha }}_4 number is larger than a certain positive value α_c ≈ 0.2, the Peclet number will increase as ε increases.

6.2 Shape of the cross-section of the needle crystal in (x, y)-plane

The contours of the needle crystal on the (x, y)-plane for the cases ε = 0.15, {\widehat{\alpha }}_4=\mathrm{0.2,0.5} are shown in Fig. 6 (a) and (b), respectively. It is seen that the fourth fold symmetric feature of the interface morphology generated by the anisotropic surface tension in steady state of needle crystal growth can be most clearly seen over a certain region of needle interfaces behind the tip. However, such a feature can only be weakly seen in the region far away from the dendrite’s tip. This fact is also consistent with the Gibb-Thomson condition (3.10), which shows that as \xi \gg 1, the interface temperature T_{\text{S}1}=\mathcal{O}(1/\xi ). It can also be seen from the asymptotic expansion form of the solution in the far field (4.1), which shows that as \xi \gg 1, the interface shape {\eta }_1=\mathcal{O}(1/\xi ). The result implies that the correction to Ivantsov’s solution induced by the anisotropic surface tension is rather weak in the root region far away from the needle’s tip.

7 Discussion and conclusions

In the present paper, we considered the steady needle-like crystal growth with full 3-D anisotropic surface tension with m-th fold symmetry. The uniformly valid asymptotic form of the solution is obtained in terms of the general Laguerre function expansion in the finite region around the tip of needle crystal, which matches with the asymptotic expansion solution in the far field behind the tip. On the basis of the solution, we derive the effect of the anisotropy parameter on the Peclet number Pe and the interface’s morphology. The results show that two components of the variation of the interface morphology of the needle crystal are generated due to the presence of the 3-D anisotropy of surface tension: the axi-symmetric component and 4-th fold symmetric component. The numerical computations for the case of m = 4 show that a fourth fold symmetric feature of the interface morphology can be most clearly seen over a certain section of the needle interface behind the tip; such feature, however, is rather weak and faint as \xi \gg 1 in the far field. One of the most important results derived in this work is that the 3-D anisotropy of surface tension may significantly affect the Peclet number Pe. For example, our numerical calculations for the case of {\eta }_0^2/2=0.005, ε = 0.12−0.15 show that, with a small value of anisotropy parameter {\widehat{\alpha }}_4=0.005, we have the relative Peclet number \widehat{\text{Pe}}\approx 1. However, when the anisotropy parameter increases to the value {\widehat{\alpha }}_4\approx 0.2, and ε = 0.3−0.4, one has \widehat{\text{Pe}}=2-6.

The present paper does not attempt to carry out a stability analysis for the needle solution obtained. The detailed analysis and discussions for the selection and determination of the tip-velocity U_tip will be given in future papers.

References

[1] D. A. Kessler and H. Levine, Phys. Rev. A, 1987, 36(2):4123

[2] D. A. Kessler and H. Levine, Phys. Rev. Lett., 1986, 57:3069

[3] F. Liu and N. Goldefeld, Phys. Rev. A, 1990, 42(8): 5052

[4] E. A. Brener, J. Cryst. Growth, 1990, 99: 165

[5] J. J. Xu, Canad. J. Phys., 1991, 69(7): 789

[6] J. J. Xu, Phys. Rev. E, 1996, 53(5): 5031

[7] J. J. Xu, Interfacial Wave Theory of Pattern Formation:Selection of Dendrite Growth and Viscous Fingering in a Hele-Shaw Flow, Berlin: Springer-Verlag, 1997

[8] M. E. Glicksman and S. P. Marsh, “The Dendrite” in: Handbook of Crystal Growth, Volume 1: Fundamentals, Part B: Transport and Stability, edited by D. T. J. Hurle, Amsterdam, North-Holland: Elsevier Science, 1993

[9] K. F. Gurski and Jeff. B. McFadden, Proc. R. Soc. Lond. A, 2003, 459: 1

[10] C. Herring, in: Structure and Properties of Solid Surface, edited by R. Gomer and C. S. Smith, University of Chicago, 1952

[11] Jeff. B. McFadden, Private communication, Nov., 2007

[12] E. A. Brener and V. I. Melnikov, Adv. Phys., 1991, 40: 53

[13] M. Abramovitz and I. A. Stegun (Eds.), Handbook of Mathematical Functions, National Bureau of Standards, 1964

[14] J. J. Xu and D. S. Yu, J. Cryst. Growth, 1998, 187: 314

[15] J. J. Xu, J. Cryst. Growth, 2002, 245: 134

Rights & permissions

Higher Education Press and Springer-Verlag 2008