The hodograph equation for slow and fast anisotropic interface propagation.

Using the model of fast phase transitions and previously reported equation of the Gibbs-Thomson-type, we develop an equation for the anisotropic interface motion of the Herring-Gibbs-Thomson-type. The derived equation takes the form of a hodograph equation and in its particular case describes motion by mean interface curvature, the relationship 'velocity-Gibbs free energy', Klein-Gordon and Born-Infeld equations related to the anisotropic propagation of various interfaces. Comparison of the present model predictions with the molecular-dynamics simulation data on nickel crystal growth (obtained by Jeffrey J. Hoyt et al. and published in Acta Mater. 47 (1999) 3181) confirms the validity of the derived hodograph equation as applicable to the slow and fast modes of interface propagation. This article is part of the theme issue 'Transport phenomena in complex systems (part 1)'.

Introduction

Anisotropy of interfaces plays a crucial role in the formation of equilibrium shapes [1], changing of growth direction of crystals to the preferable one at a critical governing parameter [2], selection of stable mode of dendritic growth [3] and facetting with the formation of microscopic defects in the bulk of crystals [4]. Figure 1 shows a region around a dendritic tip, the growth of which has been selected and faceted due to the existence of anisotropy of solid–liquid interface. From the near equilibrium properties, these anisotropic interfaces lead to crystallographically oriented properties such as free energy, magnetization and surface tension [5]. From the kinetic properties, the interface mobility and coefficients of atomic attachment to the interface strongly depend on the crystallographic faces orientation [6,7].

Figure 1

Anisotropic interfaces of crystals in solidified Ni-Cu droplet processed in electromagnetic levitation facility [5]. (a) Dendritic patterns. (b) Region around the tip of one of the dendrites with clearly visible steps. (c) Visible steps at the dendrite surface at higher resolution.

In the last century there has been great success in the description of equilibrium shapes and their motion. Herring [8] generalized the formalism of Wulff on thermodynamics of crystal surfaces. Cahn & Hoffman [9,10] introduced a useful construction, the so-called -vector, to characterize equilibrium shapes and missing orientations of sharp interfaces. Caginalp and then Kobayashi [11-13] were the first to introduce the anisotropy in the diffuse interface formalism with the surface tension anisotropy by allowing the gradient energy factor to be dependent on the orientation of the phase interface. Wheeler & McFadden [14-16] provided description of anisotropic interfaces using the Cahn–Hoffman -vector within the framework of the phase field model. Avoiding unstable interfacial orientations, the growth of highly anisotropic interfaces with facetting and edges was analysed by Eggleston et al. [17] and Debierre et al. [18]. Finally, an overview by Sekerka [19] presented approaches in the description of anisotropic crystal growth.

In the present work, we extend the description of anisotropic interfaces to the highly rapid regimes of their motion with the appearance of local non-equilibrium effects [20]. With this aim, a kinetic phase field model [21,22], which would be reduced to the single equation of motion called ‘the hodograph equation of interface’, is used [23,24]. The hodograph equation predicts non-stationary regimes as well as steady-state regimes of interface motion and it has been applied, for instance, to quantitative estimations of non-stationarity periods of dendrite growth [25]. Using the formalism of the Cahn–Hoffman -vector and the presently developed anisotropic phase field model, we follow the analysis of Wheeler & McFadden [14] who obtained the sharp interface limit corresponding to the case for which the diffuse interface width is small compared to a characteristic macroscopic length scale.

The article is organized as follows. As an advancement of the isotropic phase field model (§2) we introduce the kinetic phase field model with crystalline anisotropy in its hyperbolic formulation (§3). The obtained hodograph equation in the form of Herring–Gibbs–Thomson-type equation is compared with the previous equations of equilibrium, acceleration-velocity-dependent equation as well as Born–Infeld and Klein–Gordon equations describing the propagation of isotropic and anisotropic interfaces (§4). A quantitative comparison of linear and nonlinear equations for interface motion which follow from the obtained kinetic equation for anisotropic crystalline interface is given (§5). A summary of our conclusions is presented (§6). Finally, electronic supplementary material, [26] and appendix A add the material for the derivation of the present hodograph equation applicable to the slow and fast modes of interface propagation.

Isotropic phase field equation

The isotropic phase field equation and its solution were obtained in the work [23] and analysed in comparison with MD-data in [27]. In this section, we recollect the results of these works necessary to reproduce anisotropic advancement of the hyperbolic phase field model.

The free energy functional for the entire system of the volume in the isotropic case can be written in the form which yields the dimensional isotropic form of the dynamic equation as Here, is the phase field variable defining the phase state as is the relaxation time of the gradient flow for the phase field , is the mobility, is the solute concentration and is the temperature.

The Gibbs free energy change on transformation takes into account direction of transformation and the variety of transformations obtained for enthalpy of fusion, [28], dilute mixtures [29], and functions obtained from thermodynamics databases [30]. Note that () is the Gibbs free energy of the phase, and the indexes and are related to the liquid and solid phases, respectively.

The interpolation function and the double-well function are defined by [31] According to definition (2.3), the interpolation function varies monotonically from to .

The gradient energy coefficient in equation (2.2) is related to the interface energy at equilibrium, , for which equation (2.2) admits one-dimensional steady solution at : with the diffuse interface stationary width and under the boundary conditions as and as . Then, the surface energy of the isotropic interface is given by which is proportional to the product of the interface width and the energy per unit volume associated with the barrier height .

In the dynamics, , equation (2.2) admits one-dimensional travelling-wave solution [23,32] with the boundary conditions as and as , with the constant velocity limited by as a maximum speed of phase field propagation, the velocity-corrected effective interface thickness, and the mobility related to the isotropic interface mobility as First, the particular solution (2.9)–(2.12) with the hyperbolic tangent function follows from the general set of analytical solutions of Allen–Cahn-type equations [33] which is given by equation (2.2). Second, the interface velocity, , cannot exceed the maximum speed of disturbance propagation in the phase field, because the phase field itself dictates the interface shape and its velocity, i.e. in the solutions to equations (2.9)–(2.12). Third, with regard to the effective interface thickness (2.11), one has to note two important issues: (i) with increasing interface velocity, should become smaller than the constant interface width that has been chosen as a reference for the interface thickness in equilibrium state, equation (2.7); (ii) within the limit , one gets , therefore, the phase field variation will be steeper with the tendency to build up a sharp interface as the velocity increases.

In the dimensionless isotropic form, equation (2.2) can be re-written taking into account equations (2.8) and (2.12) as

To obtain the Herring–Gibbs–Thomson-type equation following from the hyperbolic phase field equation (2.13), we shall obtain the sharp interface limit corresponding to the case that the diffuse interface width is small compared to a characteristic macroscopic length scale , which can be taken to represent a typical radius of interfacial curvature in the non-planar case [14]. More specifically, we consider the distinguished limit while maintaining a finite value for the ratio in order that the surface energy remains finite in the sharp interface limit. Using these scales, it is convenient to choose the length and coordinates in units of , the energy density and energy barrier height in unit of latent heat , the time and relaxation time in units of a diffusive time scale consistent with the accompanying energy diffusion processes, the surface energy in units of and the interface kinetics given by in units of . Then, the isotropic form of equation (2.2) follows from equation (2.13) as where plays the role of a small parameter in the subsequent asymptotic treatment and it is defined by This ratio means that the high value of barrier between phases provides a smallness of and the asymptotic limit .

Hyperbolic phase field model with anisotropy

Free energy and scaled phase field equation

Consider a binary system consisting of solvent and solute under isothermal condition, with the temperature being constant in the overall system. The system is undergoing phase transition, solidification/melting, from the undercooled/overheated state. Taking the existence of the anisotropic diffuse solid–liquid interface into account [34,35], the free energy functional in units of for entire system of the volume is described as where length, time, Gibbs potential, and homogeneous functions of surface energy, , are measured in units of , and , respectively, and is the surface energy of the interface depending on the normal vector pointing from solid to liquid, The homogeneous extension of is given by [36] such that the free energy functional (3.1) becomes where the Gibbs potential is described by the expression [37,38] which has the local equilibrium contribution and the local non-equilibrium contribution In addition to the interpolation function and the double-well function given by equation (2.5), the contributions (3.6) and (3.7) include: the phenomenological coefficients proportional to the relaxation time of the gradient flow , the barrier between phases as well as the Gibbs free energies and .

A stable evolution of the entire system is given by the Lyapunov condition of non-positive change of the total Gibbs free energy. For the functional (3.1), this condition gives the inequality from which one finds the following phase field equation (appendix A) Note that and are measured in units of and , respectively. Here, the -vector of Hoffman & Cahn [9] is described by where the second equality in equation (3.10) defines the -vector as a homogeneous function.

The hodograph equation as a form of the Herring–Gibbs–Thomson-type equation

The anisotropic form of the phase field equation can be written if we assume that the surface energy and the interface mobility are homogeneous functions of the first degree in order to include their anisotropy, i.e. and , as for -vector equation (3.10). Indeed, taking into account these conditions and using the same treatments from §2 for equation (2.2), one can find that the anisotropic equation (3.9) takes the following form: where equation (2.17) holds for the definition of small parameter . If the inertial term is neglected, in equation (3.11), the difference in Gibbs free energy is taken as [28] (with the melting temperature) and the equality (3.11) is accounted for, then the anisotropic phase field equation of Wheeler & McFadden [14] is recovered.

Using asymptotic analysis and some necessary derivations (see electronic supplementary material, [26], especially, eqn (83)), one gets where is the orientation dependent interface mobility, and stand, respectively, for the velocity and acceleration normal to the interface, and the divergency of the -vector is with () being new coordinates depending on curvilinear coordinates from which both principal curvatures () can be obtained (see eqns (90)–(93) in electronic supplementary material, [26]). Note that the velocities and are measured in and the acceleration is measured in .

In the sharp interface description (see eqns (96)–(97) in electronic supplementary material, [26]), equation (3.13) becomes where is the term dependent on the curvatures , surface energy , orientation-dependent interface mobility and their first order derivatives with respect to . In the particular case of a convex shape having negative principal curvatures ( for ), the term plays a role of deceleration.

Equation (3.16) represents the compact form of acceleration-velocity-dependent Herring equation of a moving curved interface for the arbitrary driving force . In a general case, this equation can be considered as the hodograph equation for the equilibrium state, slow and fast propagation of anisotropic interface.

Recovering previously derived model equations

Equilibrium by Herring equation

The classic form of the Herring equation follows directly from equation (3.16) by taking , and [8]: This equation also follows from the zero variational derivative from free energy (3.1) with [4]. Such a condition defines the minimal thermodynamic potential with respect to the transfer of particles from one phase to the other. Therefore, equation (4.1) describes the shape of anisotropic particle in equilibrium, presenting a balance between a volumetric tendency to exchange and a surface tendency to save a shape of coexisting phases.

Acceleration-velocity-dependent Gibbs–Thomson equation

In the case of the isotropic phase field model, the surface energy, the interface mobility and the maximum speed of disturbances are independent of the orientation, i.e. , and have constant averaged values. Equation (3.16) is then simplifying for due to the zero derivatives and

Taking into account equations (2.7), (2.8) and (2.12) established in the case of the isotropic model, the following relations among the model parameters are obtained where we also used equation (4.3). Using solution (2.6), and taking into account equations (2.7), (2.8), (2.12) and (4.4), the time-dependent thickness of the interface is now given by [23] The inequality in equation (4.5) has the same meaning as for solutions (2.9)–(2.12): because the phase field itself dictates the interface shape and its velocity, the interface velocity cannot exceed the maximum speed of disturbance propagation in the phase field, .

By the above considerations in the isotropic scheme, we recover from equation (3.16) our generalized acceleration-velocity isotropic Gibbs–Thomson equation found in [23] and analysed for (100)-Ni-crystal orientation in [27] with the same parameters, with the negative curvature , the phase field assumes the value and in the solid and liquid phases, respectively. At small and moderate velocity, and with the absence of the driving force, , equation (4.6) arrives at In two dimensions, this equation was used by Gurtin & Podio-Guidugli [39] to explain the evolution of interface possessing effective inertia consistently with experimental data on the oscillation of quantum crystals [40].

Velocity-dependent Herring equation

In the case of absence of inertial effects, (, ) and constant normal velocity , equation (3.13) reduces to the Herring–Gibbs–Thomson equation in terms of the Cahn–Hoffman -vector that is found by Wheeler & McFadden [14] Using equation (3.15), this equation leads to the velocity-dependent Herring equation [41] which describes the interface motion of anisotropic particle due to imposed Gibbs free energy change on transformation, , and both principal curvatures, and .

A couple of important particular cases can be outlined, which follow from equation (4.9). First, if the driving force is absent (i.e. with the absence of supersaturation or supercooling), equation (4.9) predicts the interface motion by both principal curvatures The latter effect can be seen as mean curvature flow described by the Allen–Cahn equation [36,42,43] if both principal curvatures are equal to the mean curvature of the surface, i.e. () and (). Therefore, equation (4.10) reduces to Second, in the case of planar interface, , a simplest equation of motion, can be found from equation (4.9). This equation follows from the classic theory of irreversible processes and from traditional phase field model [11,34].

Born–Infeld equation

One of the specific cases of nonlinear wave propagation can be found from equation (3.16) under the absence of driving force, , and with the infinite orientation-dependent interface mobility, Equation (4.13) is the anisotropic Born–Infeld equation, which in its isotropic form, in equation (3.17), is described by with standing for the negative curvature. Equation (4.14) can be directly obtained from equation (4.6) and it is used in nonlinear electrodynamics [44,45]. Such an analogy of nonlinear waves and the fast phase interface propagation can be clarified in forthcoming works.

Klein–Gordon equation

Accepting and in equation (3.9), one gets the undamped Klein–Gordon equation extended to the propagation of anisotropic interface. As is known, to describe peculiarities of the first moments of matter formation, the Klein–Gordon equation is used for the inflation stage of the Universe [46]. In this sense, equation (4.15) can be considered as the Klein–Gordon equation for the anisotropic field propagation of the inflation stage of the matter. Using the definition of the -vector of Hoffman & Cahn (3.10), equation (4.15) can be transformed to the isotropic version of the undamped Klein–Gordon equation [23].

A simplest form of the damped Klein–Gordon equation can also be found from equation (3.9) with but within the limit : This equation describes the oscillatory motion of matter within the action of the potential function . Such a type of front propagation is used for qualitative analysis of the inflation stages of various fields [46].

Comparison with data of atomistic simulations

In this section, we shall transform from the dimensionless variables and functions (as described in §2, 3 and 4) to the dimension equations and expressions as is given by table 1.

Table 1

Notations and dimensions.

parameter	dimension
Latin:
D_ϕ, phase field diffusion coefficient	m^2 s^−1
G, Gibbs free energy density	J m^−3
ΔH_f, enthalpy of fusion	J m^−3
ℓ, velocity-dependent interface thickness	m
M_ϕ, mobility	m^3 J^−1 s^−1
ℜ, spatial scale	m
ΔS_f, entropy of fusion	J mole^−1 K^−1
T_m, melting temperature	K
t, time	s
V_n, normal interface velocity	m s^−1
A_n, normal interface acceleration	m s^−2
V_ϕ, maximum speed of the phase field	m s^−1
v_0, volume of the entire system	m^3
W_ϕ, energetic barrier between phases	J m^−3
Greek:
γ, interfacial free energy	J m^−2
δ, equilibrium interface thickness	m
ε_ϕ, gradient energy factor	(J m^−1)^1/2
κ, interfacial mean curvature	m^−1
κ_i (i=1,2), principal curvature	m^−1
μ(n), orientation-dependent interface mobility	m^4 J^−1 s^−1
μ_klm, interface kinetic coefficient	m K^−1 s^−1
ϕ, phase field variable	—
τ, time scale	s
τ_ϕ, relaxation time of the gradient flow	s

At small driving force of solidification or melting, interface velocity has a linear dependence on undercooling or overheating, respectively [6] that is also predicted by equation (4.6) in comparison with MD-data [23,24,27]. At large driving force, MD-simulation predicts interface velocity with non-linearity of two types: velocity-undercooling relationship with saturation [47] and with maximum by the crystal growth velocity [48,49]. These nonlinear dependencies in growth kinetics of crystals are well described by the solution of equation (4.6) and by the hyperbolic phase field equation (2.2) [22]. The kinetic undercooling obtained from the steady-state mode predicted by equation (4.6) included in a general undecooling balance at the tip of dendrite allows us to describe the solidification kinetics in glass-forming alloys [50].

Velocity-driving force relationship

For the planar interface () propagating with constant velocity (, , ) the hodograph equation (3.16) together with equation (3.14) yields where the phase field diffusion coefficient, , and the maximum phase field speed, , are considered as two parameters which define the relaxation time of the gradient flow of the phase field by With the local equilibrium limit, i.e. with (), equation (5.1) reduces to the linear relation , which is obtained from the parabolic phase field equation and which can exhibit nonlinear behaviour only due to nonlinearity of the function (see discussion in [22]).

Anisotropic functions

By employing cubic harmonics [51,52] and adopting established conventions [53-55], the interfacial free energy can be written as with where is the average value of the interfacial free energy, the coefficients are the strengths of anisotropy, and the basic variables and in the context of the normalization take the following form of Fehlner & Vosko [51]: where are the Cartesian components of the unit normal vector .

Hoyt et al. [56] defined the values of stiffness for different crystallographic orientations in the solid–liquid interface energy needed for parametrization of . Following the results of their work, the anisotropic function and its corresponding stiffness can be evaluated around the vanishing angle between the normal to the interface and the normal to the crystal face, i.e. for different orientations. Tables 2 and 3 summarize the obtained results for different orientations, where the first column gives the Miller indices for orientation and parallel direction. Note that the linear system of four interfacial stiffness equations (the middle column in table 3) can be solved using the computed values of interfacial stiffness in MD-simulations from ref. [57].

Table 2

Interfacial free energy for different orientations obtained analytically from equation (5.3).

orientation	γ(n)/γ_0
100[010];100[001]	1+25ϵ_1+47ϵ_2+413ϵ_3
110[11¯0];110[001]	1−110ϵ_1−1314ϵ_2+952ϵ_3
111[1¯10];111[112¯]	1−415ϵ_1+6463ϵ_2+32351ϵ_3

Table 3

Interfacial stiffness for different crystallographic orientations obtained analytically from equation (5.3) (middle column) and the value of linear fitted interfacial stiffness in MD simulations [57] (right column) for Ni.

orientation	(γ+(d^2γ/dθ^2))/γ_0	γ+(d^2γ/dθ^2)
100[010];100[001]	1−185ϵ_1−807ϵ_2−14013ϵ_3	0.177
110[11¯0]	1+3910ϵ_1+15514ϵ_2−354ϵ_3	0.405
110[001]	1−2110ϵ_1+36514ϵ_2+17552ϵ_3	0.228
111[1¯10];111[112¯]	1+125ϵ_1−128063ϵ_2+1120351ϵ_3	0.386

Application to growth of nickel crystal

Free energy for different crystal orientation

The linear system of four equations for interfacial stiffness (the middle column in table 3) can be solved using the computed values of stiffness in MD-simulation due to Rozas & Horbach [57] (right column in table 3). One gets then the following values of averaged interface energy, and anisotropic strength coefficients, , These values can be used to compute the interfacial free energies with respect to , and crystal orientation for Ni, i.e. , and , respectively, from the equations shown in table 2. The comparison with the results of works [57,58], and the values of the interfacial free energies obtained by using the average interfacial energy, from [59,60] is shown in table 4.

Table 4

Comparison between the obtained values of interfacial free energy, , and , and averaged interfacial free energy, , for Ni and those found in the literature.

γ_0	γ_100	γ_110	γ_111	source
0.302	0.325	0.310	0.304	[57]
0.2186	0.2322	0.2134	0.2126	[58]
0.3402	0.3846	0.3229	0.3216	[58]
0.25	0.269	0.257	0.252	[59] and present work
0.278	0.299	0.285	0.280	[60] and present work
0.312	0.336	0.321	0.314	present work

Velocity-undercooling relationship

In the case of a pure (chemically one component, elemental) system, the driving force, , can be given by equation (3.12). Therefore, equation (5.1) yields the following anisotropic interface ‘velocity-undercooling’ relationship from which one can identify the expression of coefficients of crystal growth kinetics, , as where , and are the Miller indices.

Results of computations

Using equation (5.7) with equations (5.3)–(5.5) and (5.8), the growth of Ni-crystal in , and 111 crystallographic directions can be characterized as follows:

Knowing the values of interfacial stiffness computed in MD-simulations, the development of the interfacial free energy and its stiffness in cubic harmonics with three strengths of anisotropy, , shown in table 2 and table 3, yield a good agreement in comparison with [57,58] and those obtained using the averaged interface energy of [59,60], see table 4.

Using material parameters of Ni from table 5, figure 2 shows very good predictions of MD-data due to Hoyt et al. [47] by the hodograph equation, (3.16), for planar interface, , and constant normal velocity given by equation (5.1), where the driving force, , is given by equation (3.12).

Table 5

Comparison of the obtained values of interface kinetic coefficient , equation (5.8), from ‘anisotropic interface velocity-undercooling’ relationship, equation (5.7), for different orientations, i.e. , -and , for Ni and those found in the literature [6,47].

μ_100	μ_110	μ_111	source
0.672	0.586	0.409	[6]
0.719	0.507	0.356	[6]
0.45	0.32	0.18	[47]
0.418	0.253	0.170	present work

Figure 2

Predictions of nonlinear and linear anisotropic ‘interface velocity-undercooling’-relationship, equation (5.7), i.e. with phase field relaxation time, , and without phase field relaxation time, , respectively, (continuous line)-and (dashed line) compared with MD-data due to Hoyt et al. [47] () for Ni. Comparison is made for: (a) crystal direction, (b) crystal direction and (c) crystal direction. The material parameters of Ni used in calculations are shown in table 6.

Table 5 shows the agreement for the obtained values of interface kinetic coefficient, equation (5.8), in comparison with those found in the literature [6,47].

Using expansion of the interface energy (5.3) and (5.4), we obtained values of kinetic coefficients in three crystallographic directions by the fitting of equations (5.7) and (5.8) to data of MD-simulation. In this respect, one should note that the values of these kinetic coefficients correlate well with the values obtained by the authors of the work [47] for the same MD-data, see table 5. This confirms the adequacy of the application of the kinetic equation (5.7) together with (5.8) to the description of anisotropic interfaces in a wide range of undercooling.

Table 6

Material parameters of Ni used in calculations.

parameter		value	source
T_m (K), melting temperature		1706	[6]
ΔS_f (J mole^−1 K^−1), entropy of fusion		10.73	[60]
ΔS_f=ΔH_fT_m (J m^−3 K^−1), entropy of fusion ^a		1.66×10^6	present work
ΔH_f (J m^−3), enthalpy of fusion		2.84×10^9	present work
parameters for (100), (110)-and (111) orientations
γ (n) (J m−2), interfacial free energyb
γ100 =0.336	γ110 =0.321	γ111 =0.314	present work
D_ϕ(n) (m^2 s^−1), phase field diffusion coefficient ^c
D_ϕ,100=8.44×10^−8	D_ϕ,110=4.87×10^−8	D_ϕ,111=3.2×10^−8	present work
V_ϕ(n) (m s^−1), maximum phase field speed ^d
V_ϕ,100=247.3	V_ϕ,110=152.1	V_ϕ,111=136.8	present work

the values of interfacial free energy used in fit for different orientations fromtable 4.

the values of phase field diffusion coefficient, , obtained in this work from fit for different orientations.

the enthalpy of fusion, , is calculated through the entropy of fusion, , due to Jian et al. [60] taking into account the melting temperature, , due to Mendelev et al. [6].

Note that for each , and orientations with the values of interfacial free energies , and , the melting temperature, and the enthalpy of fusion, , respectively (table 6), the phase field diffusion , and , and the maximum speeds of phase field propagation, , and , which are considered as parameters, are obtained by fitting MD-data due to Hoyt et al. [47]. These latter are shown in the 10 and 12 rows of table 6, respectively.

Conclusion

The present work has been devoted to the derivation and analysis of the slow and fast interface propagation described by the hodograph equation following from the phase field model. The obtained solution of the hodograph equation describes the amplitude, width and velocity of the anisotropic interface in the slow and fast regimes of dynamics. We show that the obtained hodograph equation is consistent with classical equations of the Herring–Gibbs–Thomson-type as well as the velocity-dependent Herring equation for crystal growth, Klein–Gordon equation for the front’s propagation of the inflation stage of matter, and Born–Infeld equation of the nonlinear electrodynamics. The first benchmarks of the derived hodograph equation of anisotropic interface motion show the consistency of its solution with the steady-state growth of nickel crystals obtained for different crystallographic directions using molecular dynamics simulation of Hoyt et al. [47].