# Competition between Surface Energy and Elastic Anisotropies in the Growth of Coherent Solid State Dendrites

## Abstract

A new phase field model of microstructural evolution is presented that includes the effects of elastic strain energy. The model’s thin interface behavior is investigated by mapping it onto a recent model developed by Echebarria et al [Kar04]. Exploiting this thin interface analysis the growth of solid state dendrites are simulated with diffuse interfaces and the phase field and mechanical equilibrium equations are solved in real space on an adaptive mesh. A morphological competition between surface energy anisotropy and elastic anisotropy is examined. Two dimensional simulations are reported that show that solid state dendritic structures undergo a transition from a surface dominated [10] growth direction to an elastically driven [11] growth direction by changes in the elastic anisotropy, the surface anisotropy and the supersaturation. Using the curvature and strain corrections to the equilibrium interfacial composition and linear stability theory for isotropic precipitates as calculated by Mullins and Sekerka, the dominant growth morphology is predicted.

## I Introduction

Recent years have seen the increased use of phase field modeling of free-boundary problems related to microstructure growth. This method has the advantage in that it avoids explicit front tracking by making phase boundaries spatially diffuse through the use of one or more order parameters, , which vary continuously across interfaces. This method of simulation has been very successful in modeling dendritic solidification [Boe02, Pro05]. In particular, Phase Field theory has confirmed the predictions of the microscopic solvability theory of dendrite growth [Kessler1988, Ben1984, Meiron1986]. That is, the anisotropy in the surface energy (), although small, controls the growth rate and morphology of the resultant precipitate [Kar96].

More recently the phase field approach has been used to study dendrite growth directions in cases where multiple sources of anisotropy can compete to control the resultant morphological structure. Haxhimali et al. [Haxhimali2006] showed that the dendrite growth direction in Al-Zn alloys can change continuously from (100) to (110) with the addition of solute. The authors explained their observations as an interplay between the four fold and six fold anisotropy parameters [HoytRev], was conjectured to vary as a function of solute. Analogously, Provatas et al. [provatas2003] modeled the change from (100) directed dendirites to a seaweed structure as the surface energy anisotropy competed with the direction of the temperature gradient.

The phase field method has also been used to model the effects of elastic and plastic effects on microstructural evolution [Morin1995, lowengrub1997, Aguenaou1998, Wang2002A, Wang2002B, Yeon2005, Steinbach2006]. In these studies the elastic strain energy was added to the phase field energy functional and the elastic modulus tensor () was assumed to be a function of the phase field . Many of these studies have focused on the interaction of elastic fields on growth and coarsening rates of solid state particles [Chen2004, Zhu2004]. These studies have shown that elastic effects play an important role in controlling growth rates and morphology, consistent with previous analytical predictions of Laraia et al. [Laraia1988] and Voorhees et al [Voorhees1992]. Energy minimization techniques also showed that the elastic free energy drives the interface to become unstable for mismatched elastic coefficients under both applied and self strain [Kassner2001A, Kassner2001]. Other studies have also incorporated the effect of mobile dislocations in phase field modeling [Haataja2001]. They have examined such effects as the role of dislocations on the coarsening rate during spinoidal decomposition [Mah04].

Another phenomenon of both scientific and practical interest, and that has has not received enough theoretical study, concerns the morphology and growth rate of solid state dendrites when two sources of anisotropy can compete. Anisotorpy can result from the surface energy () or the more interesting result of anisotropy emerging in the elastic modulus tensor (). Experimentally it has been found that dendritic structures can appear in the solid state [Malcolm67, Husain1999] provided specific elastic and kinetic relationships between the parent and precipitate phases are met. As discussed by Husain et al. [Husain1999], solid state dendrites are made possible if the precipitating phase and the parent phase have similar atomic lattice structures, there is a high rate of atomic transfer between the two phases and the transformation rate is controlled by diffusion in the presence of an anisotropy. More exotic structures can also be found through competing mechanisms leading to interface instabilities as shown by Yoo [Yoo2005].

This paper studies solid state dendritic growth in the presence of competitive interactions between surface energy and elastic anisotropy. It quantifies a transition between dominant growth morphologies controlled by each of these two anisotropic effects. Section II presents the phase field model used in this study. Section III examines the effect of isotropic strains on the equilibrium composition. Section IV provides an approximation to the anisotropic strain field at the interface of the precipitate. Section V illustrates the effect of the surface energy anisotropy, elastic energy anisotropy and supersaturation on the growth morphology of dendritic growth. Finally, the transition between a surface energy controlled morphology and an elastic energy controlled morphology is quantified in Section VI.

## Ii Phase Field Model

This section derives a phase field free energy for solid state transformations, which assumed a dilute alloy free energy for the bulk chemical thermodynamics. To this is added an additional contribution to account for elastic free energy. Corresponding phase field equations are then derived from this free energy. It is shown that the phase parameters can be related according according to the thin interface asymptotics of Ref. [Kar04] to a good approximation. The thin interface limit allows for significant efficiency of simulations of the phase field equations, particularly when combined with novel adaptive mesh refinement algorithms [Pro99, Badri2007].

The bulk chemical free energy is given by

(1) |

, and are the solute concentration, order parameter and total strain fields, respectively. The parameter is the gas constant, is the molar volume, is the temperature and i the melting temperature of component of a dilute binary alloy. The elastic free energy is defined by Hooke’s law, given by

(2) |

where is the phase dependent eigenstrain, given by . Here is the hydrostatic lattice eigenstrain with and is the lattice parameters in each phase. The phase dependent elastic tensor is interpolated across the interface by . A convenient choice for the interpolation function which maintains the bulk phases at is given as

(3) |

The elastic energy can be re-cast as a function of , leaving the elastic free energy in the form of

(4) |

Each pre-factor in the polynomial of (i.e., , , and ) is a function dependent on the elastic modulus tensor. Explicit expressions for these prefactors are given in the appendix (Appendix A) for the cubic elastic modulus tensor in two-dimensions.

It is straightforward to show using the common tangent technique that the equilibrium composition, denoted , is modified by elasticity according to,

(5) |

where is the equilibrium coexistence line corresponding to the parent phase in the absence of elasticity, is the correction to the phase diagram due to a local change in the strain and is the partition coefficient.

### ii.1 A Phase Field Model For Elastically Influenced Phase Transformations

By applying the usual dissipative dynamics for the order parameter,
mass conservation for concentration and strain relaxation for
the strains, the following kinetic equations for the order parameter
, undercooling (where is the partition
coefficient and is defined in below) and the strain fields
are derived:

Phase Mobility:

(6) |

Chemical Diffusion:

(7) |

Strain relaxation:

(8) |

The explicit form of in equation 6 is

(9) |

where is the equilibrium composition corrected for the strain by equation 5, i.e., with being the equilibrium composition in the absence of elasticity ( has been neglected here, as discussed below). The quantities and are the usual phase field parameters setting the dimensions of space and time [Kar96, Kar04]. The function controls the four-fold surface energy anisotropy, where is the angle of orientation between the interface normal and a reference axis.

The dimensionless undercooling in equation 7 is modified from its form in Ref. [Kar04] to incorporate an elastic correction to the equilibrium composition. This results in

(10) |

where is the prefactor to in the elastic free energy. The function modulates the diffusion through the interface correcting for the diffuse nature of and is given by

(11) |

Two-sided diffusion is controlled by the function where modulates the diffusion in the two phases, to simulate equal diffusion coefficients . The dimensional diffusion coefficient is denoted .The value of elastic coefficients used herein is a actually reported as .

The phase field equations above are derived in the limit where , the condition for a small difference in the elastic coefficients in either phase. This condition also holds for all materials in which holds. Generally speaking even in the most extreme disparities of the elastic coefficients, so this assumption is not unreasonable.

Following the procedure for the derivation of the phase field model presented in [Kar04], the constants , and are inter-related by the asymptotic analysis of Ref. [Kar96]. This analysis maps the phase field model above, without elasticity, onto the sharp interface limit governed by the Gibbs-Thomson condition , where is the capillary length, is interface kinetic coefficient and is the concentration jump across the to phase interface. Echebarria et al. showed that for the case of vanishing kinetic coefficient (appropriate for the study of solid state denderites) the following relations must be obeyed: and where is the diffusion coefficient in dimensional units (). These relationships arise by an expansion of the phase () and composition () inside the interface, which is matched to solutions outside the interface such as to emulate the sharp interface boundary conditions.

It expected that for the phase field model presented in this work the thin interface relationships between , and will be, to lowest order, the same as those for the model of Echebarria et al. based on two observations. The first is drawn from the work of Yeon et al. [Dong2005], who showed that, to first order, the inclusion of elasticity is decoupled from the curvature and kinetic effects in the Gibbs-Thomson condition, which reads . In their asymptotic matching procedure the derivation of and remained unaffected by the presence of the contribution.

The second observations validating the use of the thin interface analysis of Ref. [Kar04] for the phase field model presented here is that for our simulations it was foind that in Equation 7 and in equation 6. The first quantity was found numerically to be at least an order of magnitude smaller than the other terms in its equation. It was likewise also found that (which includes ). Both and are zero in the limit of zero strain, but are also small for small lattice eigenstrains. For larger strains these terms are still valid as long as the interfacial velocity () is small. In effect a small velocity eliminates any excess kinetic effects that arise at the interface. In this work only small eigenstrains () and small growth rates () are considered.

It should be noted that the small variable
accounts for the effect of a variable strain field on the
development of the precipitate. It plays an analogous
role to the temperature correction of that is used
when simulating solidification of alloys with *non-linear*
coexistence phase boundaries [Tong2007].

## Iii Growth of an Isotropic Second Phase Precipitate with Coherent Interfaces in an Isotropic Parent Phase

This section examines the interface compositions of an isolated second phase precipitated into a parent phase, which is grown under isotropic conditions in both the parent and precipitate phases. The surface energy is made isotropic by setting the surface energy anisotropy coefficient . The elasticity equations are formulated in terms of cubic tensor coefficients, ie , and as described in appendix A. For isotropic linear elastic coefficients the cubic elastic terms are related by,

(12) |

where for this simulation and ( is dimensionalized by ) and the coherent hydrostatic eigenstrain is set to . The convergence constant of the phase field equations is . The equilibrium composition is with an initial alloy composition of and the solute partition coefficient is set to a value of . The diffusion coefficients are set equal in both phases thus removing the phase dependency in the diffusion coefficient (ie. ) in equation 7.

The total domain size simulated is on a side with periodic boundary conditions. The phase field and diffusion equations are solved using an explicit time stepping algorithm on an adaptive mesh with a grid spacing of at the lowest level of refinement and a time step of . The displacement field is solved by direct Gauss-Seidel iteration, which is found to require () operations per time step on an adaptive mesh [MikeThesis2008].

The precipitate particle is grown and a cross section of the composition field solved by Equation 7 is shown in Figure 1. This figure indicates the corresponding value of the equilibrium interfacial composition as calculated by equation 5. is plotted as a function of the phase by interpolating it through the interface to its corresponding precipitate side value by . The composition field is interpolated to the interface described by the point where both from inside the precipitate bulk and from outside the precipitate in the parent phase. The points are denoted and and are found to have excellent agreement with the interpolations to the center of the interface.

## Iv Approximation of the Strain Field Around a Precipitate with Cubic Elastic Coefficients

In this section the elastic field is analyzed around a circular precipitate where the cubic elastic coefficients are equal in both the precipitate and matrix phases. The anisotropy is entered into the cubic elastic coefficients by introducing a deviation from the isotropic relation in equation 12 as defined here by,

(13) |

where is the deviation from the isotropic elastic coefficients. and remain unchanged. While the analytical solution to the isotropic strain field has been derived for elliptical inclusions under a hydrostatic eigenstrain [Eshelby1957], a solution to an anisotropic precipitate under the same conditions is mathematically cumbersome do deal with [Ru2003] and is instead solved here numerically.

A circular precipitate in a parent phase with identical elastic coefficients in both phases is considered. The dimensionless elastic coefficients are set to values of and ( is dimensionalized by , ie. ), the hydrostatic eigenstrain is set to and the concentration field is made constant (The elasticity here is not influenced by compositional effects and therefore any concentration field will produce similar results). The deviation from elastic isotropy is studied by two controls, the particle radius and the strength of the deviation from isotropic elasticity ().

The behaviour of the strain trace due to changes in the strength of the elastic anisotropy () is studied by holding the precipitate radius constant while increasing the value of . The amplitude of the trace of the strain tensor ( is found to have a strong sensitivity to variations in . This is illustrated in Figure 2 for , where the individual strain components and their sum at the precipitate/matrix interface are plotted as a function of the angle , with zero representing the (10) direction. Notice the four fold symmetry of the trace of the strain. Although the strain energy depends strongly on the value of , the precipitate radius is found to have little effect on the amplitude of any perturbation to the strain trace.

The dependence of the strain trace can be approximated to lowest order by a single Fourier mode defined by,

(14) |

To measure the functional form for the amplitude of the strain trace (), the anisotropic strength () is varied and the amplitude of the strain trace () is measured for the corresponding waveform at the interface of the precipitate. These measured values are fitted to a functional form, given by the equation,

(15) |

The functional form of equation 15 shows good agreement for all values of in the regime where as shown in figure 3. This range of anisotropies is well within the limits of the relative anisotropic strengths that are studied here.

## V Conditions Influencing the Morphology of Precipitates

This section characterizes three distinct controlling influences on the selection of a dominant morphology of precipitated dendrites. These are the anisotropies of the elastic tensor and surface energy and the supersaturation. The effect of each of these parameters are systematically tested by increasing the strength of each parameter while holding the other parameters constant.

In the following phase field simulations, and, as required by the sharp interface analysis [Kar04], the dimensionless diffusion coefficient is set to . The temperature is set such that the equilibrium composition is and the partition coefficient used is . The elastic coefficients were converted to units of the model by the elastic modulating factor and in these units the elastic coefficients are set to a value of and , where .

Precipitate structures are grown in a system with periodic boundaries, where the system size is set to ( being the interface width, determined along with from the asymptotic analysis used). The precipitates grew to sizes of at most and the solution to the displacement field drops off as . This justifies the claim that purely isolated precipitates are studied while using the periodic boundaries. The diffusion coefficients and elastic coefficients have no phase dependence. The grid spacing is set to and the explicit time step is set to .

### v.1 Elastic Anisotropy ()

The elastic anisotropy emerges from the elastic tensor (, and ). The anisotropy of the tensor is varied by holding , constant and varying through changes in . Figure 4 (a-c) shows the effect of increasing the strength of the elastic anisotropy by increasing while holding surface energy anisotropy and supersaturation (controlled via the average composition ) constant. From top to bottom the values of used are respectively. As can be seen in this figure, a small elastic anisotropy causes the surface energy to dominate and the dendrite grows in the direction ( Figure 4 (a) ). When the elastic anisotropy is increased to sufficient strength ( Figure 4 (c) ) the dendrite grows in the direction. When the anisotropies effectively destructively interact the resultant structure leads to an almost isotropic growth morphology ( Figure 4 (b) ).

(a) , and |

(b) , and |

(c) , and |

*(Colour Online) Morphologies of growth precipitates vs. variation in the elastic anisotropy in the modulus tensor through an alteration in . The composition field is plotted on the adaptive mesh, blue being low concentration and yellow is a higher concentration. The total system size is , plotted is an area of on one side. .*

### v.2 Surface Energy Anisotropy ()

(a) , and |

(b) , and |

(c) , and |

*(Colour Online) Morphologies of growth precipitates for a variation in the surface energy anisotropy by modification of . The composition field is plotted on the adaptive mesh, blue being low concentration and yellow is a higher concentration. The total system size is , plotted is an area of on one side.*

The surface energy anisotropy is entered into the model using the simple form for 4-fold surface energy . The effect of surface energy anisotropy is examined by varying and holding both (or equivalently the super saturation ) and the elastic anisotropy () constant. Figure 5 (a-c) shows the effect of increasing the strength of the surface anisotropy under these conditions. From top to bottom in Figure 5 the values of the used are respectively. Analogously with the effect shown in section V.1, increasing the strength of causes the morphology to transform from a preferential growth along the [11] direction ( Figure 5 (a) ) to that of the [10] direction ( Figure 5 (c) )with a transition region where the precipitate structure is isotropic ( Figure 5 (b) ).

### v.3 Supersaturation()

(a) , and |

(b) , and |

c) , and |

*(Colour Online) Morphologies of growth precipitates for a variation in the supersaturation with both anisotropic strengths, and , held constant. The composition field is plotted on the adaptive mesh, blue being low concentration and yellow is a higher concentration. The total system size is , plotted is an area of on one side.*

The supersaturation is varied by changing the initial (i.e. average) alloy composition . Figure 6 (a-c) shows the effect of decreasing the supersaturation while holding the strength of the surface anisotropy and elastic anisotropy constant. From top to bottom the values of the initial alloy composition used are respectively. For large super saturations (Figure 6 (a)) the growth direction is dominated along directions preferred by the surface energy, i.e. the [10] directions. As the supersaturation is decreased a transition from the [10] growth direction to the [11] direction is observed ( Figure 6 (b) and (c) ).

## Vi Characterization of the Morphological Transition

In this section we discuss a technique by which the dominant precipitate growth direction can be predicted. First, the point at which this transition occurs is defined and measured by examining the envelope of the precipitate tips in space. For a given value of the precipitate radius, analysis of the amplitude of the envelope makes it possible to determine the critical surface energy anisotropy (), for a specified elastic anisotropy (), where a morphological transition from to growth directions occurs. The elastic strain energy is, however, proportional to the precipitate area (in 2D) and the surface energy contribution varies as the particle perimeter. As such the critical is also a function of the precipitate size. The particle size dependency of vs. is then found by balancing the Gibbs-Thomson corrections corresponding to surface energy vs. elastic anisotropy. This critical radius scale is found to be proportional to the Mullins-Sekerka instability radius (). Finally, a condition relating as a function of and at the transition is proposed.

### vi.1 Defining the Transition Point

*The transition point that characterizes the controlling
mechanism of growth morphology is defined as the point at which all competing
anisotropies exactly cancel.* Under this condition an isolated
precipitate will grow (ideally) as a circle (a sphere in three
dimensions) until the interface becomes unstable by the
Mullins-Sekerka instability. While in this case the interface will
become unstable, the envelope around the particle will continue to
grow as a spheroid. It is this envelope that allows the point of
transition between anisotropically controlled directions to be characterized.

The concept of the precipitate envelope is illustrated in Figure 7 where the interface for 4 precipitates with values of , and in space are plotted, where . The dashed line shows the envelope surrounding the interface. As the magnitude of approaches the transition point the amplitude of the envelope decreases, approaching zero.

### vi.2 Measurement of the Critical Surface Energy Anisotropy -

The critical surface energy anisotropic coefficient (denoted ) is defined as the value of at a given supersaturation () and elastic anisotropy () which results in an envelope amplitude of zero. is interpolated from the amplitudes of the precipitate envelopes obtained by varying for given values of the supersaturation and elastic anisotropy .

The envelope amplitude is approximated by measuring the difference of the total growth distance from the center of the precipitate along the [10] direction to the growth distance along the [11] direction. The transition point is the interpolated value for such that these amplitudes approach zero. Seven different supersaturations are considered here, where is varied between and and is varied from to . The envelope amplitudes are measured at arbitrary times, chosen in each case, however, such that the precipitate has outgrown any initial transients.

For each value of at each alloy composition , the envelope amplitude is plotted vs . The inset to Figure 8 illustrates this for an alloy with an average composition of , and the deviation from elastic isotropy is characterized for values of . For each value of , the data is fitted linearly and is interpolated to the transition line to extract the critical surface anisotropic value .

### vi.3 Critical Tip Radii at the Transition Point

The competing anisotropic effects controlling morphology cancel when all correction terms that are dependent on the interface normal angle () in the interfacial equilibrium composition ( ) exactly cancel. Here and are the isotropic corrections to interfacial equilibrium composition and and represent the strength of the surface energy and elastic energy anisotropies. Assuming a linear fourier expansion with a 4-fold symmetry in both and (ie. ) the terms in the interface solute correction can be grouped by order of the fourier expansion.

(17) |

where and are the relative anisotropic strengths of the surface energy correction and the elastic energy correction respectively. The factor of in the capillarity term comes from the stiffness of the capillarity, . The elastic anisotropy strength, , is linked to the strength of the elastic anisotropy through and is derived below (see Equation 21).

For an isotropic morphology to emerge the coefficient of in Equation 17 is required to vanish, ie.

(18) |

In Equation 17 the capillary term () contains a curvature correction, while the elastic term () does not. When these terms balance each other a curvature is selected, which will be associated with a critical radius (ie. ). This is determined next.

The capillarity correction is , which is used to calculate as

(19) |

The total elastic correction for cubic coefficients is defined by Equation 5 as . In the absence of elastic anisotropy the strain trace () is zero and therefore

(20) |

is calculated by consideration of the anisotropy in the strain field by substituting Equation 14 for , giving

(21) |

where and is defined by Equation 15. This results in becoming

(22) |

Substituting Equations 19 and 22 into Equation 18, the critical radius of curvature required to maintain isotropic conditions is given by,

(23) |

A fit to a selected critical radii is attained by substituting the fitted equation for the critical surface anisotropy coefficient (Equation 16) and by choosing a reference point of . This results in a relationship for the magnitude of the transition radius that depends on concentration through given as,

(24) |

The next subsection discusses a separate method by which this radius is estimated without the need for measured values of .

### vi.4 Calculation of the Transition Tip Radius using Linear Stability Approximation

In the previous section a selected precipitate radius is derived based on the interpolated values for the critical surface energy anisotropy. However, this method relates to only once the value of is measured. A method to approximate the selected radius by theoretical consideration of the Mullins-Sekerka linear stability analysis on an isotropic particle is now shown.

Mullins and Sekerka in 1963 [Mul63, Lan80] predicted a critical particle size (, is the instability mode and is the supersaturation) after which the particle interface becomes unstable. A particle at the transition point in our study can be considered to behave similarly to an isotropic particle in a supercooled matrix and the tip’s radius of curvature is assumed to be proportional to this critical radius. Here the supersaturation of the precipitate under elastic strain is modified by the elasticity according to , where is the equilibrium interface composition modified due to elasticity by Equation 5. This supersaturation is used in the linear stability analysis result to predict a minimum critical radius of,

(25) |

### vi.5 Calculation of the Critical Transition Point

With the scale of the critical tip radius determined by linear stability theory as a function of supersaturation, the required measurement of the relationship between and can be eliminated in Equation 23. This is done by substituting the linear stability prediction of the critical radius (Equation 26) into the equation for the selected transition radius (Equation 23). The resulting relationship is solved for the critical surface energy anisotropy pre-factor (). This relationship is

(27) |

Equation 27 defines a morphological transition line as a function of supersaturation (), elastic anisotropy () and the anisotropy of the capillarity (). The transition lines for are plotted in figure 10. Precipitates grown above the transition line will grow in the [10] directions while growth for conditions below the line will grow along [11] directions. Some morphologies are overplotted above and below the transition line for in figure 11 to further illustrate the utility of Equation 27.

## Vii Summary

We have introduced a phase-field model for the study of the morphological development of elastically stressed solid state precipitates. We considered particles that have coherent interfaces and are under elastic self stress by a lattice mismatch eigenstrain. We used a new finite difference based adaptive mesh refinement algorithm to solve the phase-field and strain relaxation equations thereby allowing for very rapid solution times. By investigating the effects of supersaturation, elastic anisotropy in the elastic tensor and anisotropy in the capillarity we developed a scaling relationship to predict which anisotropy will be dominant in the morphological evolution of the precipitate. It is interesting note is the effect of supersaturation on the selected precipitate morphology and growth directions.

We would like to thank the National Science and Engineering Research Council of Canada (NSERC) for financial support and SHARC-NET for supercomputing support.

## Appendix A Cubic Elastic Free Energy Coefficients

In the generalized elastic portion of the phase field free energy, as described by equation 4 , several unknown terms , , and are introduced. These coeffiecients are dependent on the particular values of the elastic modulus tensor in either of the precipitate or matrix phases. Presented here are the explicit forms for these functions for two sided cubic modula and a hydrostatic elastic eigenstrain of the form,

The zeroth order component() has no dependence on the phase at all and is calculated to be

(28) |

This pre-factor has no dependence of phase (nor concentration) and therefore it does not appear in either the phase mobility equation (equation 6) or the chemical diffusion equation (Equation 7) since the growth kinetics are dependent on differences of energy. It does however appear in the static elasticity equation (Equation 8).

The first order component is the most prominent term in the model equations and is calculated to be

(29) |

The second order coefficient to the elastic energy in terms of the phase is calculated to be

(30) |

and in the presence of equal elastic coefficients this term becomes a constant.

The third order component is calculated to be

(31) |

and has no dependence on the dynamic strain field. In the presence of equal elastic coefficients this term vanishes completely.