Phase-field simulations with inhomogeneous elasticity: Comparison with an atomic-scale method and application to superalloys

Available online at www.sciencedirect.com

Acta Materialia 58 (2010) 4170–4181 www.elsevier.com/locate/actamat

Phase-field simulations with inhomogeneous elasticity: Comparison with an atomic-scale method and application to superalloys Guillaume Boussinot a,b, Yann Le Bouar a,*, Alphonse Finel a a

Laboratoire d’Etude des Microstructures, CNRS/ONERA, 29 avenue de la division Leclerc, B.P. 72, F-92322 Chaˆtillon, France b Institut fu¨r Festko¨rperforschung, Forschungszentrum Ju¨lich, D-52425 Ju¨lich, Germany Received 21 December 2009; received in revised form 2 April 2010; accepted 4 April 2010 Available online 11 May 2010

Abstract We present a 2D and 3D phase-field analysis of microstructure evolution in the presence of a lattice misfit and with inhomogeneous elastic constants. The method is first critically compared with a Monte Carlo modeling at the atomic scale. We then apply the phase-field model to the Ni–Al system under external load along a cubic axis. We find that the microstructure becomes anisotropic and that the situation qualitatively differs depending on the sign of the applied stress. The microstructure evolution operates mainly by shape changes and alignments of precipitates, but also by splitting of precipitates initially elongated along directions perpendicular to the stress-induced, elastically favorable directions. The final microstructure is finally qualitatively analyzed in terms of a mean field theory in which the elastic inhomogeneity is embedded into an effective eigenstrain. This analysis leads to a simple formulation which can be used to easily predict the coherent microstructural anisotropy induced by any external loading condition. Ó 2010 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field models; Phase transformations; Elastic behavior; Coarsening; Mean field analysis

1. Introduction Ni-based superalloys are currently used in the hottest parts of aeroplane engines. Turbine blades, commonly made of a single crystal, are exposed to strong centrifugal forces (of the order of 100 MPa). In industry, much effort is devoted to improving the mechanical properties of such alloys by controlling their microstructure, and to increase the lifetime during service. This mainly consists in optimizing their chemical composition and performing complex thermal treatments. Mechanical properties are controlled by the plastic behavior, which, in turn, is mainly governed by the underlying microstructure. The typical microstructure of Nibased superalloys consists of c0 precipitates surrounded by channels of c matrix. The most important property of this microstructure is to localize the plastic activity in the *

Corresponding author. Fax: +33 146734155. E-mail address: [email protected] (Y. Le Bouar).

thin c channels. Mechanical properties can be further improved by the precipitation of small particles in the c channels using specific thermal treatments [1,2] (leading to bimodal microstructures), or by the addition of heavy elements which produce strong local distortions and modify the lattice parameters of the phases [3]. During service, the alloy experiences large mechanical stresses and thermal constraints that lead to a creep regime during which the microstructure develops directional coarsening known as rafting. These morphological evolutions operate at a scale of a few hundreds of nanometers, and therefore are best described by mesoscopic approaches such as phase-field models [2,4,5]. In this paper, we present our work on the phase-field simulations of Ni–Al alloy, which is often used to model Ni-based superalloys. We first present our phase-field model. In this model, particular attention will be paid to the treatment of the elastic inhomogeneity between the two phases. This is a crucial point in order to take into account external loading in our simulations. However, in

1359-6454/$36.00 Ó 2010 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2010.04.008

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

the presence of inhomogeneous elasticity the solution of the mechanical equilibrium does not reduce to a linear problem and a numerical iterative scheme is required. In addition, it is known that these elastic inhomogeneities often lie at the origin of the development of complex instabilities and morphologies displayed by growing microstructures. It is therefore important to test the reliability of the iterative scheme that will be used to solve the elastic equilibrium in the phase-field approach. This will be done here through a careful comparison between our continuous modeling and an atomistic Monte Carlo scheme [13], in a situation where complex growth instabilities are observed. This overall approach is then used in the specific case of c  c0 in the Ni–Al system. First, 2D simulations are performed in order to investigate the differences between homogeneous and inhomogeneous elasticity on microstructures large enough to be statistically representative, and also to analyze directional coarsening on a sufficiently large-scale. Secondly, we present 3D simulations of c  c0 microstructural evolution under external load in order to investigate the influence of the elastic driving force on the directional coarsening. This is to our knowledge the first 3D simulation taking into account simultaneously the elastic inhomogeneity and the ordered nature of the c0 phase. We will analyze the obtained microstructure under external load in terms of a mean-field approach of inhomogeneous elasticity in order to better understand the underlying physical processes. 2. Phase-field model Phase-field models for the simulation of c  c0 phase transitions have been widely developed [2,6,7]. They incorporate a concentration field c(r) and three long-range order (LRO) parameters gi(r), that may be defined using a concentration-wave description of the atomic configurations. In the case of large c0 volume fraction, it is crucial to introduce these LRO parameters to distinguish the different variants of the L12 ordered structure. The phase-field model relies on the time evolution of these fields. When one wants to reproduce a nucleation regime (a precipitation from a metastable, weakly supersaturated matrix), one has to use a stochastic phase-field model, in order to overcome the nucleation free energy barrier [8]. Here, we want to reproduce a Ni–Al c  c0 microstructure with a large c0 volume fraction. The initial disordered matrix will contain a large supersaturation, and will be essentially unstable. Therefore a simple spatial noise in the initial configuration will be sufficient to initiate the phase separation. As usual, the conserved concentration obeys the Cahn– Hilliard equation:

@gi ðrÞ dF ; ¼ L @t dgi ðrÞ

4171

ð2Þ

where L is a constant. In order to implement these evolution equations, one has to specify the free energy F . This is the goal of the next section. 2.1. Free energy of the system The free energy F of the system is written as follows: Z f GL ðrÞdV þ Est ; F ¼ V

where fGL is the usual Ginzburg–Landau free energy density and Est, the residual elastic energy, will be discussed below in the framework of inhomogeneous elasticity. 2.1.1. Ginzburg–Landau free energy density The Ginzburg–Landau contribution consists in an homogeneous free energy density f[c(r),gi(r)] and an heterogeneous term penalizing the variations of the fields: k bX jrgi ðrÞj2 : ð3Þ fGL ðrÞ ¼ f ½cðrÞ; gi ðrÞ þ jrcðrÞj2 þ 2 2 i The coefficients k and b should be chosen in order to reproduce experimental or calculated interfacial energies (in the latter case special attention has to be paid to the nature of the interface with respect to local elastic relaxations [2]). One should notice that the model is isotropic, which is consistent with the observed spherical shape of the c0 precipitates when they are small. Concerning the homogeneous free energy density f, the simplest polynomial adapted to our situation is: P 2 A B g 2 f ðc; g1 ; g2 ; g3 Þ ¼ ðc  cc Þ þ ðcs  cÞ i i 3 2 2 P C D i g4i : ð4Þ  g 1 g2 g3 þ 3 4 3 In this equation, the parameters have to be chosen to reproduce the phase diagram. For a sake of simplicity, we operated in such a way that free energies of c and c0 are equal at c  c0 equilibrium. We also impose that the order parameters gi of equilibrium c0 phase have an amplitude equal to 1. This leads to:   B ¼ 2A cc0  cc   C ¼ 6A cc0  cc ðcs  cc Þ ð5Þ    D ¼ 2A cc0  cc cc0 þ 2cs  3cc

ð1Þ

where ccðc0 Þ is the equilibrium concentration of c(c0 ) phase, cs the concentration below which the c phase is unstable, and A is the energy scale, which remains a free parameter at this point.

where it is assumed that the mobility coefficient M is constant. The non-conserved gi parameters follow a simple Allen–Cahn equation:

2.1.2. Elastic contribution The elastic energy is written in the framework of linear elasticity:

@cðrÞ dF ¼ Mr2 ; @t dcðrÞ

4172

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

Z n o 1 Est ¼  V raijij þ d 3 rkijkl ðrÞ ij ðrÞ  0ij ðrÞ 2 r    kl ðrÞ  0kl ðrÞ

ð6Þ

where the first term is the work corresponding to the applied stress raij on the sample of volume V with ij being the resulting strain. In the second term, 0ij ðrÞ is the local eigenstrain and ij(r) the actual strain. The elastic constants kijkl(r) are space dependent through their concentration dependence:   ðcÞ i 1 h ðc0 Þ ðcðrÞ  cc Þkijkl þ cc0  cðrÞ kijkl ð7Þ kijkl ðrÞ ¼ c c0  c c ðcÞ

ðc0 Þ

where kijkl and kijkl (cc and cc0 ) are elastic constants (equilibrium concentrations) of c and c0 , respectively. One then conveniently defines the average elastic constants by:   ðcÞ i 1 h ðc0 Þ  ðc  cc Þkijkl þ cc0  c kijkl ð8Þ kijkl ¼ cc0  cc with c the average value of c(r), and the elastic inhomogeneities: i 1 h ðc0 Þ ðcÞ kijkl  kijkl ð9Þ k0ijkl ¼ cc0  cc so that kijkl ðrÞ ¼  kijkl þ k0ijkl DcðrÞ

ð10Þ

with DcðrÞ ¼ cðrÞ  c. The eigenstrain 0ij ðrÞ is written in the form of a Vegard law: 0ij ðrÞ

00

¼  dij DcðrÞ:

ð11Þ

where the Vegard coefficient 00 is defined by 00 ¼

d c c0  c c

ð12Þ

In our case of phase-transformation controlled by atomic diffusion, static elasticity (elastic equilibrium) can be assumed since elastic waves will relax instantaneously with respect to the concentration field. This means that the sound velocity is considered infinite. The solution of the elastic equilibrium is then obtained by minimizing the elastic energy with respect to the elastic displacement field together with boundary conditions. If one writes the strain tensor ij(r) as the sum of its average value ij and its deviation dij(r), the elastic equilibrium is written : ð14Þ

The average strain ij is then determined by the boundary conditions. If one fixes the displacement at the boundaries, ij is the corresponding applied deformation aij . If one chooses to impose a stress raij at the boundaries, then ij is solution of : dEst ¼ raij dij

ð17Þ

which reads:

ij ¼ S ijkl k0klmn ð00 dmn hDc2 ðrÞi  hdmn ðrÞDcðrÞiÞ þ rakl ; ð18Þ with S ijkl kklpq ¼ dip djq

and

1 hwðrÞi ¼ V

Z

d 3 rwðrÞ:

ð19Þ

r

In the homogeneous case of equal elastic constants ð kijkl Þ in both phases, the elastic equilibrium is written: 2 kijkl @ uk ðrÞ ¼ r0 @DcðrÞ ij @rj @rl @rj

ð20Þ

with r0ij ¼ kijkl 00 dkl ij ¼ S ijkl rakl

ð21Þ

ð15Þ

ð22Þ

The displacement uk(r) is easily found using the explicit Green function G of linear elasticity in the Fourier space [9] defined through its inverse:  G1 ij ðqÞ ¼ kikjl qk ql :

ð23Þ

In this case, the average strain ij does not appear in the elastic equilibrium equation. Therefore, the heterogeneous strain dij(r) does not depend on ij , and the elastic energy is written as: Z 1 d 3q 2 BðqÞjcðqÞj ; ð24Þ Est ¼ 2 q – 0 ð2pÞ3 where c(q) is the Fourier transform of the concentration field: Z 1 d 3 rcðrÞeiqr ; ð25Þ cðqÞ ¼ 3=2 ð2pÞ r and B(q) is given by: BðqÞ ¼ r0ij 00 dij  qi r0ij Gjk ðqÞr0kl ql :

where ui(r) is defined by 1 dij ðrÞ ¼ ½@ui ðrÞ=@xj þ @uj ðrÞ=@xi Þ: 2

ð16Þ

The boundary condition is

and where the misfit d is linked to the lattice parameters ac and ac0 by: a c0  a c d¼2 : ð13Þ a c0 þ a c

dEst ¼0 dui ðrÞ

Eq. (14) reads:    2 kijkl @ þ k0 @ DcðrÞ @ uk ðrÞ ijkl @rj @rl @rj @rl

 @DcðrÞ @Dc2 ðrÞ þ k0ijkl 00 dkl : ¼ kijkl 00 dkl  k0ijklkl @rj @rj

ð26Þ

For isotropic elasticity, B does not depend on q. For anisotropic elasticity, B only depends on n = q/jqj. The second

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

term on the right-hand side of Eq. (26) is the configurationdependent contribution to the elastic interactions. One can then define a configuration-dependent elastic energy: Z 1 d 3q 2 q r0 Gjk ðqÞr0kl ql jcðqÞj ð27Þ Econf ¼  2 q – 0 ð2pÞ3 i ij The minimum of elastic energy is reached when the microstructure shows parallel plate-like precipitates (with a characteristic length in one direction much smaller than in the two others) which are perpendicular to the q-vector nmin that minimizes B(q) [9]. In the case of a negative anisotropy coefficient n = C11  C12  2C44, the soft directions correspond to the three cubic directions. Here, an important conclusion has to be drawn. The microstructure response to elastic stresses depends on $2dEst/d c(r). According to Eq. (6) the elastic driving force in the homogeneous case is dEst =dcðrÞ ¼ kijkl 00 dij ½kl þ dkl ðrÞ  00 dkl DcðrÞ. In the present situation (Vegard’s law) and if homogeneous elasticity is considered, dij(r) does not depend on ij . As a consequence, $2 dEst/ dc(r) does not depend on ij either. Then the microstructure evolution will not depend on fixed boundary conditions

ij ¼ aij or fixed applied load (see Eq. (22)). In the case of inhomogeneous elasticity, the heterogeneous strain field dij(r) and the average strain ij are coupled by Eqs. (16) and (18). In the homogeneous case (k0 = 0), the differential operator in front of uk(r) in Eq. (16) does not depend on space and is local in the Fourier space. This is why ui(r) can be found directly using the closed-form Green function Gij(q). For the inhomogeneous case, the differential operator in Eq. (16) is space dependent and the equation has to be solved numerically. One technique is to use an iterative algorithm based on the Green function of the homogeneous elasticity Gij [11,12]. One initializes u to u(0), the analytical solution of the homogeneous elasticity equilibrium:  @2 @Dc ð0Þ  uk ðrÞ ¼  kijkl 00 dkl ðrÞ; ð28Þ kijkl @rj @rl @rj ð0Þ

The average strain  is initialized to ij ¼ aij if a strain is ð0Þ applied, or to ij ¼ S ijkl rakl if a stress is imposed. Then, one finds iteratively the solution of the elastic equilibrium with inhomogeneous elasticity by solving: ðnþ1Þ  @ 2 uk ðrÞ  00 ðnÞ @DcðrÞ  ¼ kijkl  dkl  k0ijklkl kijkl @rj @rl @rj @Dc2 ðrÞ  k0ijkl @rj ! ðnÞ @uk ðrÞ DcðrÞ ; @rl

þ

k0ijkl 00 dkl



@ @rj

ð29Þ

ðnÞ

with ij ¼ aij if a strain is applied, but with an iteratively calculated average strain if a stress rakl is imposed:

a ijðnþ1Þ ¼ S ijkl k0klmn ð00 dmn hDc2 ðrÞi  hdðnÞ mn ðrÞDcðrÞiÞ þ rkl : ð30Þ

4173

In Ref. [11], the boundary condition consisted of ij ¼ 0, and consequently the last equation was not solved. The number of iterations n needed to converge obviously depends on the differences between the elastic constants of the two phases. In the case of Ni–Al (or Ni-based superalloys), these differences are quite small, of the order of 10% ðk0 =k ’ 10%Þ. Then, only one iteration is usually needed to obtain a reasonable convergence. As a summary, we have seen in this section that the investigation of external loading in c  c0 microstructure evolution necessitates taking into account the inhomogeneous elasticity within the system. Whereas the elastic driving force has an analytic expression in the case of homogeneous elasticity, it has to be found numerically for the inhomogeneous case. This is done through an iterative scheme using the Green function of the homogeneous elasticity problem associated with the average elastic constants. In the next section, we investigate the ability of this procedure to reproduce morphological evolutions of precipitates presenting large differences of elastic constants with respect to the matrix. This is done through comparison with Monte Carlo simulations. 3. Inhomogeneous elasticity: Comparison with Monte Carlo In order to verify the validity of the inhomogeneous elasticity algorithm introduced in the phase-field model, we make a comparison with Monte Carlo simulations. In Ref. [13], Lee developed a Monte Carlo algorithm in which elastic interactions are introduced via a quadratic potential energy between atoms. We refer to Ref. [13] for more details of the numerical method. The comparison concerns the 2D evolution of a single precipitate toward its equilibrium shape [14]. In the phase-field model used here, we simulate a simple phase separation and a single concentration field is defined. Hence only the Cahn–Hilliard equation (Eq. (1)) is solved. In the Cahn–Hilliard equation, a small noise is added in order to overcome possible energy barriers during the evolution. This noise has to be locally conserved and white. The total free energy of the system is written as:   Z k 2 3 d r f ½cðrÞ þ jrcðrÞj þ Est ð31Þ F ¼ 2 r where

"  2  4 # 1 5 1 f ðcÞ ¼ A  c  þ : c 2 2 2

ð32Þ

The elastic energy of the system without applied stress is written as: Z n o  1 d 3 rkijkl ðrÞ ij ðrÞ  0ij ðrÞ kl ðrÞ  0kl ðrÞ ð33Þ Est ¼ 2 r The elastic interactions are introduced in the Cahn–Hilliard equation through the functional derivative of Eq. (6) with respect to c(r). For numerical purposes, we intro-

4174

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

duce a dimensionless free energy F =A, and we discretize with a grid step d. Therefore the experimental interface ene through ergy C is linked to its dimensionless counterpart C e C ¼ CAd. Hence a choice of A fixes d (or vice versa). In the first simulation, matrix and precipitate are elastically isotropic. Numerical inputs are chosen to reproduce the simulations presented in Ref. [14]: 00 = 0.05, lM = 5.97  1019Ja2, lP = lM/2 (the precipitate is soft) and C = 2.5  1021J a1 with a the interatomic distance in the atomistic simulation. li ¼ C i12 ¼ C i44 are the shear modulii of the matrix (i = M) and the precipitate (i = P) and C i11 ¼ 3li (note that since we are in 2D, elastic constants are in units of J a2 and the interface energy is in units of J a1). The elastic inhomogeneity between precipitate and matrix is significantly large and four iterations are needed to achieve convergence within the algorithm used for mechanical equilibrium. Finally, as we want to compare our model with an atomistic approach, we use a small grid step d = 1.4a. Initially, the precipitate has a circular shape, and its radius equals 35a. The results of the first simulation are shown in Fig. 1, where we present the evolution according to the two methods: the phase-field is shown in the top row, and the Monte Carlo counterpart in the bottom row. Even though the two methods are completely different (a continuous description for the phase-field method and a discrete description for the atomistic method), we observe a very good agreement for the evolution of the precipitate shape. Indeed, in both cases, the circular shape is unstable and a wave appears at the surface of the precipitate. The wavelength is the same in both simulations, which suggests that this destabilization mechanism is an intrinsic property of the system. Furthermore, a few growing lobes appear and finally the precipitate stabilizes in an elliptical shape. One should note that the ellipse can be oriented along any direction. Here the selection of the orientation is a result of the finite size of the box and the periodic boundary conditions. In the second comparison also, a soft precipitate is embedded in the matrix but the system is elastically aniso-

tropic. The anisotropy coefficient n = 1.13 C44 is the same in both phases (i = P,M). The elastic constants are given 19 M J a2 , and C Pij ¼ C M by: C M 44 ¼ C 12 ¼ 8:36  10 ij =2. Again, four iterations are needed for convergence of the inhomogeneous elasticity algorithm. The interface energy is the same as in the previous simulation. The phase-field grid step is chosen equal to the interatomic distance of the atomistic model (d = a). The precipitate is initially circular and its radius is 40d. In Fig. 2, the evolution of the precipitate towards equilibrium is presented for the phase-field simulation (top row) and for the Monte Carlo simulation (bottom row). The circular shape of the precipitate evolves towards a 4-fold symmetry shape, with concave edges. The 4-fold symmetry remains for a certain time before being destabilized, and the precipitate adopts an elliptical shape, aligned with the [1 0] direction. Of course, the [0 1] direction is equivalent, and the noise added to the Cahn–Hilliard equation is responsible for the symmetry breaking. Note that at the end of the simulation, the largest dimension of the precipitate becomes comparable with the size of the simulation box. The further evolution would be partly determined by the boundary conditions, and we stopped the simulation before the precipitate reached the limit of the box. However, again, the evolution of the precipitate shows a very good agreement between the phase-field model and the atomistic model. In conclusion, we note that the continuous phase-field model reproduces remarkably well the kinetic path and the elastic behavior observed in the Monte Carlo simulation even though the continuous modeling has been used down to the atomic scale. This is a very satisfactory result for the validation of the mechanical equilibrium, especially considering that the elastic inhomogeneity was very large (matrix twice as hard as the precipitate). 4. Inhomogeneous elasticity in Ni–Al In this section, we present an analysis of the microstructural evolution in a model Ni–Al superalloy under external

Fig. 1. Evolution of an isotropic precipitate in an isotropic matrix from phase-field simulation (top) and Monte Carlo simulation (bottom) [14]. The elastic constants of the precipitates (P) are related to the ones of the matrix (M) by: C Pij ¼ C M ij =2. The phase-field simulations correspond to the dimensionless times 0, 2000, 4000, 8000, 12000 and 48000. Monte Carlo steps (MCS) are expressed in millions (M).

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

4175

Fig. 2. Evolution of a system with anisotropic elasticity with phase-field method (top) and Monte Carlo simulation (bottom) [14]. The anisotropy coefficient n = 1.13 (see text for definition) for both phases, and the inhomogeneity is given by C Pij ¼ C M ij =2. The time is given in dimensionless units in the phase-field simulation and the Monte Carlo steps (MCS) are expressed in millions (M).

load. The numerical simulation are based on the phasefield model presented in Section 2. We first detail how we link the parameters of the model to the material constants. We then present the simulation results in 2D and 3D. The first step is to select proper values for the parameters of the Landau functional. For T = 1073 K, the equilibrium concentrations are cc = 13.5%Al and cc0 ¼ 23:4%Al [15]. We need also to select the value of cs which controls the stability of the c phase. Indeed, from Eq. (4), it is obvious that the c phase is unstable for c > cs if B > 0. In the case of high c0 phase fraction, the precipitation mechanism is of spinodal type. Hence, no temporal noise is required to overcome free energy barriers, and a simple spatial noise is sufficient to initiate the precipitation. The early evolution (exponential amplification of waves) depends on the curvatures of the metastable parts of the free energy and the characteristic time for this first regime then depends on the calibration of the functional. However, non-linear effects rapidly come into play, and this first precipitation regime is very short. The subsequent evolution is essentially controlled by the diffusion between precipitates. Therefore it is not crucial to reproduce quantitatively the early evolution during the precipitation regime. We will then simply impose cs < c for all the simulations. More precisely, we choose, cs = 15% of Al. The Ginzburg term (penalizing the interfaces) in Eq. (3) is controlled by k and b which in turn enter into the interfacial energies. An important feature of Ni-based superalloys is that the anti-phase boundary (APB) energies between different variants are usually much larger than twice the c/c0 interfacial energy, and therefore no APB are observed. To this regard, we choose k/ad2 = 0.28 and b/ad2 = 0.0013, which leads to sufficiently diffuse interfaces and to APB energies larger than twice the c  c0 interface energy. Concerning the elastic parameters, we choose the following set of elastic constants: C c11 ¼ 190 GPa; C c12 ¼ 138

0

GPa; C c44 ¼ 110 GPa for c phase, and C c11 ¼ 210 GPa; c0 c0 C 12 ¼ 150 GPa; C 44 ¼ 124 GPa for c0 phase. These values are typical of Ni–Al alloys [16,17]. The inhomogeneity is approximately 10%, and within the iterative algorithm described above, one iteration is usually sufficient to achieve mechanical equilibrium. At T = 1073 K, the lattice misfit is 0.37%, which leads to a Vegard coefficient 00 = 0.037. Finally for the timescale, according to the formalism in Section 2, we have to link M and the diffusion coefficient D in the disordered phase using:  2 @ f : ð34Þ D¼M @c2 c We chose D = 1017 m2 s1 which is very close to the impurity diffusion coefficient of Al and Ni at T = 1073 K in pure Ni [18], and of Al in Ni3Al at the same temperature [19]. One important effect that is automatically incorporated in our model is the influence of elasticity on the equilibrium concentrations in the two-phase region of the coherent phase diagram. Specifically, in the case of homogeneous elasticity and of Vegard’s law, the equilibrium in the presence of coherency effects corresponds to plates with a residual elastic energy: Z 1 2 d 3 rBmin ðDcðrÞÞ : ð35Þ 2 r This elastic energy will compete with the incoherent Landau energy (Eq. (4)), and leads to new equilibrium concentrations. Obviously, the shift of the new equilibrium concentrations will depend on the second derivatives of the Landau polynomial. However, in the way we parameterize the present model, we do not have any precise information about these second derivatives, and the shifts that will result cannot be linked to any precise experimental data. Therefore, in order to suppress these uncontrolled effects, we subtract Bmin(Dc)2/2 from Eq. (4), so that the

4176

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

coherent equilibrium concentrations are cc and cc0 . Of course, this procedure is only approximate in the case of inhomogeneous elasticity. However here, since the inhomogeneity is small, this approach is still sufficient to control equilibrium concentrations at coherent equilibrium. Concerning the c  c0 interface energy, we use, for the 3D calculations, C = 5 mJ m2. This is consistent with experimental measurements [20] and with preliminary ab initio calculations reported in Ref. [2], for which a procedure of local relaxation close to the interface was performed. This relaxation procedure accounts for the displacements that are not included in the linear elasticity theory introduced in the phase-field simulations. An important feature of the c  c0 microstructural evolution is the competition between surface energy, which favors spherical precipitates when this energy is isotropic, and elastic energy, which favors precipitates with large interfaces perpendicular to the elastically soft directions. This competition leads to a characteristic length scale lc, which may be estimated as the typical size of a precipitate when it undergoes a transition from sphere-to-cube. This critical size depends, of course, on the dimensionality. More precisely, we have the following estimation for the sphere-to-cube transition [9]1: l3D c 

K C 0 0:0601 þ 0:0479D C 44 d2

ð36Þ

where the dimensionless constants K and D0 depend on the elastic constants through n D0 ¼ C 11 þ 2C 12 þ 4C 44 C 11 C 44 ðC 11 þ C 12 þ 2C 44 Þ K¼ ; ð37Þ jnjðC 11 þ 2C 12 Þ2 where we have used the anisotropy constant n ¼ C 11  C 12  2C 44 : We use a similar approach to evaluate the circle-to-square transition for a plain-strain condition in two dimensions; this gives 00 l2D c  2:33K

C C 44 d2

given in Ref. [20] and the one inherited from the preliminary ab initio calculations reported in Ref. [2]. We have seen in Section 2.1.2 that the elastic response to an external load requires taking into account the difference in elastic constants between c and c0 . Before analyzing the influence of external load on microstructural evolution, we first focus on the influence of inhomogeneous elasticity on a microstructure without external load. In Fig. 3, we present a c/c0 microstructure, with homogeneous elasticity (left) and inhomogeneous elasticity (right). The different variants of the L12 structure of the c0 phase are represented by different colors, and the c matrix is represented in grey. The grid step is equal to d = 64 nm, and we use a 512  512 simulation box, corresponding to a linear size of 32.77 lm. In the homogeneous elasticity simulation, the elastic constants of c and c0 are equal to the weighted average of the elastic constants presented above (see definition of  k in Eq. (8)). The initial configuration for the two simulations are the same, and the two configurations displayed in Fig. 3 correspond to the same evolution time. The c  c0 interfaces (and thus the matrix channels) are perpendicular to the cubic directions according to the negative elastic anisotropy n. The c0 phase fraction is 75% (according to the lever rule) and only small channels of c matrix are present. One observes qualitative differences between the two simulations. Generally [11,21], when the precipitates are harder than the matrix, they have a tendency to display compact and regular shapes. Together with elastic anisotropy, this should lead to regular alignments between precipitates and to a slow down of the coagulation mechanism. Those effects are indeed observed in Fig. 3 where we have identified different locations where elastic inhomogeneity leads to qualitative difference between the two simulations. These simulation results show that an elastic inhomogeneity as low as 10% is sufficient to significantly modify the shape of the c0 precipitates. We now present the microstructural evolution under external load in an elastically inhomogeneous c  c0 microstructure. We first perform 2D simulations to characterize the evolution of a microstructure containing a large num-

ð38Þ

where the dimensionless constant K00 is K 00 ¼

ðC 11 þ C 12 þ 2C 44 ÞC 44 C 11 nðC 11 þ C 12 Þ

2

Using the elastic constants given above, we get 3D l2D c  0:24lc :

ð39Þ

Therefore, in order to reproduce the 3D critical length in our 2D simulations, we have used an interface energy equal to 20 mJ m2, which is about four times larger than the one

1 The numerical constants that we have obtained in the expression of the critical length given by Eq. (36) differ from the ones given in Ref. [9].

Fig. 3. Comparison of c/c0 microstructure in the homogeneous elasticity case (left) and in the inhomogeneous elasticity case (right). The c0 precipitates are represented in blue, green, red and white (four translational variants) and the c is shown in grey. The c0 elastic constants are 10% higher than c elastic constants (see text for precise values). The highlighted areas show qualitative differences between the two simulations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

Fig. 4. Evolution of the c(black)–c0 (white) microstructure under uniaxial strain of 1.86% in the [1 0] direction. The time is given in dimensionless units d2/D. The system, which is initially (t = 0) the results of a stress-free evolution, evolves towards a rafted microstructure, where the symmetry between the two cubic directions is broken.

Average phase fraction

1.0

0.8

0.6

0.4

0.2

0.0 0

128

256

384

512

Position along (10) 1.0

Average phase fraction

ber of precipitates. The simulation box consists in 512  512 grid points with a grid spacing d = 20 nm. The initial configuration (see Fig. 4 at t = 0) is a microstructure resulting from a stress-free evolution with a phase fraction of c0 precipitates equal to 35%. As in Fig. 3, the precipitates exhibit interfaces perpendicular to the cubic directions and are strongly aligned along those directions. This is more pronounced than in Fig. 3 because of the smaller precipitate phase fraction. A large positive uniaxial strain app is then applied in the horizontal direction [1 0] and maintained until the end of the simulation. We used a fairly large applied strain in order to accelerate the symmetry breaking and thus reduce the simulation time. The resulting microstructural evolution is presented in Fig. 4 using shades of gray according to the value of the concentration field (c0 phase appears in white and c phase in black). During the microstructural evolution, the smallest precipitates disappear according to classical Ostwald ripening. Precipitates initially elongated in the horizontal direction, or having more or less square shape, elongate and align along the [1 0] direction. Precipitates initially strongly elongated in the [0 1] direction split into different precipitates. The splitting occurs because these precipitates are initially elongated along a direction perpendicular to the stress-induced, elastically favorable directions. Note that, as it is a local mechanism, this splitting will be observed also for a misaligned single precipitate [14]. The resulting microstructure consists of rafts containing the loading direction as observed experimentally when the lattice misfit between c and c0 and the load are positive. At long evolution times, some rafts exhibit “snaky” shapes due to the coalescence of the microstructure. This 2D simulation with a large applied strain allows us to see on a large-scale the c  c0 microstructural evolution under external load. In Fig. 5, we can observe the symme-

4177

0.8

0.6

0.4

0.2

0.0 0

128

256

384

512

Position along (01) Fig. 5. Plot of the average c0 phase fraction along the [1 0] direction (top) and along the [0 1] direction (bottom). The dashed lines correspond to the initial configuration of Fig. 4 and solid lines correspond to the final one.

try breaking observed in Fig. 4 by plotting the average c0 phase fraction along the [1 0] direction and along the [0 1] direction (bottom and top of the Fig. 5, respectively). The dashed lines correspond to the initial configuration of the evolution shown in Fig. 4 and the solid lines to the final one. For the initial state, the amplitude of the oscillations along the [1 0] and [0 1] directions are similar, showing the equivalence of the two cubic directions. In the final state, oscillations with large amplitude are observed along the [0 1] direction, whereas an almost constant value of the c0 phase fraction is observed along the [1 0] direction. Indeed, as compared to the initial configuration, the standard deviation of the phase fraction is twice as large in the [0 1] direction, whereas it is one order of magnitude lower in the [1 0] direction. The maxima of the c0 phase fraction along the [0 1] direction correspond to the location of the rafts that extend over the whole simulation box in the [1 0] direction. The further evolution of c0 phase fraction along the [0 1] direction should lead to a broadening of the peaks corresponding to the rafts and an increasing distance between them. Finally, note that the c0 phase fraction decreases significantly (about 10%) between the initial stress-free configuration and the final rafted microstructure presented in Fig. 4. Even if the phase fraction is expected to evolve under an applied stress, the evolution may be overesti-

4178

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

mated here, mainly due to the large value of the applied stress used in the present simulations. We present now our 3D simulations of c  c0 microstructure when an external stress is applied to the system. Under this condition, Eqs. (16) and (18) need to be solved simultaneously within the iterative algorithm summarized in Eqs. (29) and (30). The initial configuration consists of an unstable homogeneous c phase with concentration c ¼ 20%, and a stress is applied during the whole simulation. We use a 64  64  64 simulation box and a grid spacing d = 4 nm, which leads to a linear size of 256 nm. The situation in Fig. 6 corresponds to 13 h of evolution under a compression of 600 MPa along the [1 0 0] direction. The precipitates are displayed in gray, and two different views are presented. According to the left panel in Fig. 6, we can clearly observe a rafted microstructure with platelike groups of precipitates perpendicular to the loading axis. In the right panel in Fig. 6, a view perpendicular to the loading axis is presented. Because the transverse sizes of the precipitates in the plane perpendicular to the loading axis are only slightly above the elastic crossover length between sphere and cube l3D c  20 nm , we observe a moderate alignment perpendicular to the elastic soft directions. We now apply a tension of 600 MPa along the [1 0 0] direction, starting again from an unstable homogeneous c phase with c ¼ 20%. In Fig. 7, the microstructure is displayed after 13 h of evolution. We clearly observe rod-like precipitates (with a characteristic length in one direction much larger than in the two others) parallel to the loading axis. We note that rod sections display circular shapes when they are small enough and square-like shapes for larger sections. This kind of rod-like microstructure has been observed experimentally [22]. The 3D simulation results of c  c0 microstructural evolution in the Ni–Al system under an external stress (Figs. 6 and 7) have clearly shown a qualitative difference resulting from the sign of the applied stress (compression or tension). In the case of compression, the c0 precipitates are elongated perpendicular to the loading axis, while they

Fig. 7. Microstructure obtained after 13 h of evolution under a 600 MPa uniaxial tension along the direction [1 0 0]. The microstructure exhibits rod-like precipitates parallel to the loading axis [1 0 0].

are elongated parallel to the loading axis in the case of tension. This anisotropic coarsening of the microstructure is often referred as type N (P) rafting when the precipitates are perpendicular (parallel) to the loading axis. In the next section, we discuss the origin of the anisotropy of this anisotropic coarsening. 5. Discussion The microstructure evolution observed during external loading presented in Figs. 6 and 7 is the complex result of the competition between several driving forces: the lattice misfit, the elastic inhomogeneity and the applied stress. In this section, we show that this problem can be drastically simplified in the limit of small elastic inhomogeneity. We apply the mean-field approach developed in Ref. [10], which provides an analytical treatment of the problem posed in Section 2.1.2 for small elastic inhomogeneity. This

Fig. 6. Microstructure obtained after 13 h of evolution under a 600 MPa uniaxial compression along the direction [1 0 0]. The microstructure exhibits plate-like groups of precipitates (rafts) perpendicular to the loading axis [1 0 0].

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

4179

approach allows a qualitative discussion for the results obtained in the previous section. In the mean-field approach, the space-dependent operator that operates on uk in Eq. (16):    @2 @ @ 0  ð40Þ þ kijkl DcðrÞ kijkl @rj @rl @rj @rl

By analogy with homogeneous elasticity, one can use Eq. (46) to define an effective eigenstrain ij by:

is approximated by its homogeneous expression:

In this approximation where the Green’s function of homogeneous elasticity is used to solve the elastic equilibrium, the coupling between the microstructure and the external load appears through an effective stress–free strain * which includes the elastic inhomogeneity k0 . As mentioned in Section 2.1.2, in the case of homogeneous elasticity, the minimization of elastic energy corresponds to plate-like precipitates perpendicular to the soft directions. These soft directions are functions of the elastic constants and of the components of the eigenstrain. As seen in Eq. (49), the effective stress–free eigenstrain ij depends on the average strain ij . In the case of an applied strain aij , ij should be replaced by aij . On the other hand, in the case of an applied external stress ra, the average strain ij is coupled to the inhomogeneous displacement fields uk(r), and should be identified by solving simultaneously Eqs. (18) and (42). However, if the condition   ð50Þ rakl  k0klmn 00 dmn hDc2 ðrÞi  hdmn ðrÞDcðrÞi

2

 kijkl

@ ; @rj @rl

ð41Þ

which is the Green’s function of the homogeneous elasticity problem (Eq. (23)). On can then obtain explicitly the displacement uk(r) in Fourier space by solving: 2

 @DcðrÞ kijkl @ uk ðrÞ ¼  kijkl 00 dkl  k0ijklkl @rj @rl @rj þ k0ijkl 00 dkl

@Dc2 ðrÞ : @rj

ð42Þ

We show now that, within this mean-field approximation and as long as the configuration-dependent elastic energy is involved, our problem of inhomogeneous elasticity under external load can be mapped into a problem of homogeneous elasticity without external stress and characterized by an effective eigenstrain that will be identified in the following. In order to write the elastic energy in the form of Eq. (24), we need to find B*(q), the equivalent of B(q) given by Eq. (26). Eq. (42) should then be written in a similar form as Eq. (20). For this purpose, one makes the change of variable: /ðrÞ ¼

cðrÞ  cc : cc0  cc

ð43Þ

In a sharp-interface limit, one has /2(r) = /(r), since / = 1 in c0 and / = 0 in c. Then Eq. (42) gives: 2 n   kijkl @ uk ðrÞ ¼ ð kijkl 00 dkl  k0ijklkl Þ cc0  cc @rj @rl o @/   ðrÞ: þk0ijkl 00 dkl cc0  cc cc0 þ cc  2c @rj ð44Þ

ij ¼ S ijkl rkl

ð48Þ

This leads to:



ij ¼ 00 dij  S ijkl k0klmn mn  00 dmn ðcc0 þ cc  2cÞ

ð49Þ

is fulfilled, Eqs. (18) and (42) may be decoupled and the average strain ij may be approximated by Eq. (22). As both terms in the right-hand side of Eq. (50) are of the same order, this inequality is valid if:  2 ð51Þ ra  k0 00 cc0  cc For the Ni–Al system under study here, this means ra  10 MPa, which is consistent with the applied stress of 600 MPa in the 3D simulations of previous section. Following this line, the effective eigenstrain ij may be calculated through Eq. (49) with the approximation: ij  S ijkl rakl :

ð52Þ

In the case of uniaxial loading along a cubic direction, say [1 0 0], and taking into account the cubic anisotropy of the elastic constants, it is easy to show that the effective eigenstrain ij takes a simple tetragonal form, i.e.: 3 2  11 0 0 7 6 ð53Þ c ¼ 4 0 22 0 5

In terms of concentration field, the last equation is equivalent to: @ 2 uk ðrÞ n  00  ¼ ðkijkl  dkl  k0ijklkl Þ kijkl @rj @rl  o @Dc ðrÞ ð45Þ þk0ijkl 00 dkl cc0 þ cc  2c @rj which is of a similar form as Eq. (20) with r0ij replaced by:   kijkl 00 dkl  k0 kl Þ þ k0 00 dkl cc0 þ cc  2c : ð46Þ r ¼ ð

11  22 ¼ ðS 11  S 12 Þ ðk011  k012 Þra

The configuration-dependent elastic energy is then written as: Z 1 d 3q 2 q r Gjk ðqÞrkl ql jcðqÞj ð47Þ Econf ¼  2 q – 0 ð2pÞ3 i ij

where ra is the external stress along the [1 0 0] direction and where the Voigt notation has been used for the compliance S ijkl associated with the average elastic constants and for the elastic inhomogeneity k0ijkl . This result shows that, under external stress along a cubic direction, the increase in

ij

ijkl

ijkl

0

0

22

and that the difference between the diagonal components 11 and 22 is given by: 2

ð54Þ

4180

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

the degeneracy for the three elastically soft cubic directions depends only on k011  k012 , i.e. the difference between the shear components C0 = (C11  C12)/2 of precipitate and matrix, and not on the precipitate volume fraction, in agreement with previous studies [23–26]. More specifically, if the [1 0 0] external stress ra is not too large, the diagonal components of ij will stay positive (we recall that the misfit is positive for Ni–Al), and in that case the new elastically soft directions will be given by:   if ra k011  k012 < 0 ½1 0 0 direction   if ra k011  k012 > 0 ½0 1 0 and ½0 0 1 directions In the present situation, we have k011  k012 = 2(C0 (c0 )  C0 (c)) = 8 GPa > 0. Therefore, the elastically soft direction is [1 0 0] under a [1 0 0] compressive loading, and is degenerated between [0 1 0] and [0 0 1] under a [1 0 0] tensile loading. In order to be more specific and to quantitatively analyze the 3D microstructures presented above, we present now a quantitative analysis of the effective eigenstrain ij in a NixAl1x alloy under uniaxial loading. In the case of ra = 600 MPa in compression (Fig. 6), the effective stress–free strain tensor ij can be calculated by Eq. (49) with mn given by Eq. (52): 3 2 0:046 0 0 7 6 ð55Þ c ¼ 4 0 0:032 0 5 0 0 0:032 This effective stress–free strain tensor corresponds to a cubic-to-tetragonal transformation, and in this case where c c c xx > yy ; zz , [1 0 0] is the unique soft direction [9]. Thus, minimization of elastic energy requires plate-like precipitates perpendicular to the [1 0 0] direction. This is consistent with our result where a rafted microstructure made of plate-like c0 precipitates perpendicular to [1 0 0] axis is obc c served in Fig. 6. One can consider c xx ; yy and zz as an effective lattice misfit along the cubic directions [1 0 0], [0 1 0] and [0 0 1], respectively. The effective lattice misfit is then larger along the [1 0 0] direction than along the [0 1 0] and [0 0 1] directions, and the system tends to avoid interfaces parallel to the [1 0 0] axis. In the case of a tensile stress ra = 600 MPa (Fig. 7), the effective stress–free strain is: 2 3 0:022 0 0 6 7 ð56Þ t ¼ 4 0 0:037 0 5 0

0

0:037 t xx

t In this case where < t yy ; zz , [0 1 0] and [0 0 1] are the soft directions . Therefore, the minimum of elastic energy is reached for plate-like precipitates perpendicular to [0 1 0] or [0 0 1]. In both cases, the precipitates are parallel to [1 0 0], due to an effective lattice misfit t xx smaller than in the other directions. Before the minimization of elastic energy is reached, one can expect the precipitates to adopt both soft directions and to be rod-like with a square or rectangular section in the (y z) plane. For even shorter evolution (or for small characteristic lengths), one can ex-

pect the section of the rod-like precipitates to be circular, as a result of a predominant isotropic interface energy. This is fully consistent with the results presented in Fig. 7. The above results and discussion shows that the determination of the soft directions corresponding to the effective stress–free strain tensor ij (see Eq. (49)) is an efficient tool to predict the morphology of the c  c0 microstructure resulting from our 3D phase-field calculations. This strategy relies on the hypothesis of small elastic inhomogeneity, an hypothesis clearly fulfilled in Ni-based superalloys where the elastic inhomogeneity usually does not exceed 10%. Note also that this approach is valid only if the applied stress ra is not too small, more precisely ra  10 MPa in the present situation, which corresponds to the vast majority of situations experienced by superalloys during service. As mentioned in the Introduction, our simulations and the previous prediction of rafted patterns are only based on elastic interactions. It is well known that plasticity has an important role in the driving force for rafting [27–29]. Recent efforts to model the microstructural evolution in superalloys are based on phase-field models incorporating a plastic activity, either using discrete dislocations [30,31], or an effective plastic strain [32], or viscoplastic laws similar to the ones used in crystalline plasticity [33,34]. Finally, we would like to mention that the effective mean field theory presented above, which is based on the identification of an effective eigenstrain ij incorporating the external loading (see Eqs. (49) and (52)), can be used in the presence of plastic activity. More precisely, this approach consists in using Eq. (49), where the original stress–free eigenstrain 00 is replaced by a “plastically relaxed” eigenstrain ~00 ¼ 00  p , where p is an effective plastic strain whose symmetries are dictated by the geometry of the activated slip systems. 6. Conclusion In this paper, we present a phase-field approach to model microstructural evolutions in materials where the elastic inhomogeneity between coexisting phases has to be taken into account. We first compare the phase-field model for inhomogeneous elasticity with atomistic Monte Carlo simulations from the literature. We show that phase-field model reproduces remarkably well the kinetic path and the elastic behavior observed in the Monte Carlo simulation even though the continuous modeling has been used down to the atomic scale. Moreover, considering that the elastic inhomogeneity was very large (matrix twice as hard as the precipitate), these simulation results can be considered as a stringent test for the accuracy of the fixed-point algorithm used to solve the mechanical equilibrium. In the second part of the paper, we apply the phase-field model to study the microstructure evolution in the Ni–Al system, using both 2D and 3D simulations. The effect of elastic inhomogeneity is first demonstrated with a stressfree simulation where we find that precipitates harder than

G. Boussinot et al. / Acta Materialia 58 (2010) 4170–4181

the matrix tend to adopt more compact shapes. We then study the influence of an imposed strain on a large-scale 2D microstructure where an anisotropic coarsening is observed. This anisotropic rafting is also predicted in 3D simulations when a stress is applied to the microstructure. The situation qualitatively differs depending on the sign of the applied stress, revealing the complexity of the phenomenon of rafting in 3D. We finally show that the role of an applied external stress on the microstructural evolution can be understood in terms of a mean field theory in which the elastic inhomogeneity is embedded into an effective eigenstrain. We show in particular that, under uniaxial loading along a cubic direction, the anisotropy of the microstructure depends only on the difference between the shear elastic constants C0 = (C11  C12)/2 of precipitates and matrix and not on the volume fraction. Note finally that the mean field theory presented here could be applied to qualitatively predict the microstructure evolution during traction or compression tests along any direction, and that an homogeneous plastic activity could also easily be included. Acknowledgement The authors are grateful to Dr. A. Gaubert for her help when performing the 3D phase-field simulations. References [1] Babu SS, Miller MK, Vitek JM, David SA. Acta Mater 2001;49:4149. [2] Boussinot G, Finel A, Le Bouar Y. Acta Mater 2009;57:921. [3] Epishin A, Bruckner U, Portella PD, Link T. Scripta Mater 2003;48:455. [4] Le Bouar Y, Loiseau A, Khachaturyan AG. Acta Mater 1998;46:2777. [5] Singer-Loginova I, Singer HM. Rep Prog Phys 2008;71:106501. [6] Wang Y, Banerjee D, Su CC, Khachaturyan AG. Acta Mater 1998;46:2983. [7] Zhu JZ, Wang T, Ardell AJ, Zhou SH, Liu ZK, Chen LQ. Acta Mater 2004;52:2837. [8] Bronchart Q, Le Bouar Y, Finel A. Phys Rev Lett 2008;100:015702.

4181

[9] Khachaturyan AG. Theory of structural transformations in solids. New York: Wiley; 1983. [10] Khachaturyan AG, Semenovskaya S, Tsakalakos T. Phys Rev B 1995;52:15909. [11] Hu SY, Chen LQ. Acta Mater 2001;49:2001. [12] Gururajan MG, Abinandanan TA. Acta Mater 2007;55:5015. [13] Lee JK. Scripta Mater 1995;32:559. [14] Choy JH, Hackney SA, Lee JK. In: Chen LQ et al., editors. Mathematics of microstructure evolution. EMPMD monograph series 4. TMS/SIAM Publications; 1996. p. 101. [15] Wang J, Osawa M, Yokokawa T, Harada H, Enomoto M. Superalloys 2004. TMS 2004. [16] Yoo MH, Daw MS, Baskes MI. In: Vitek V, Srolovitz DJ, editors. Atomistic simulation of materials: beyond pair potential. New York: Plenum Press; 1989. p. 401. [17] Yoo MH. Acta Metall 1987;35:1559. [18] Landolt-Bo¨rnstein. In: Mehrer H, editor. Numerical data and functional relationships in science and technology, vol. 26. Berlin: Springer; 1990. [19] Fujiwara K., Watanabe M., Horita Z., Nemoto N., Noumi K., Simozaki T., International Conference on Solid–Solid Transformations, JIMIC-3, Koiwa M., Otsuka K., Miyazaki T., 1999. [20] Ardell AJ, Ozolins V. Nature Mater 2005;4:309. [21] Onuki A. In: Chen LQ et al., editors. Mathematics of microstructure evolution. EMPMD monograph series 4. TMS/SIAM Publications; 1996. p. 87. [22] Nabarro FRN. Metal Mater Trans A 1996;27:513. [23] Pineau A. Acta Metall 1976;24:559. [24] Johnson WC, Berkenpas MB, Laughlin DE. Acta Metall 1988;36(12):3149. [25] Nabarro FRN, Cress CM, Kotschy. Acta Mater 1996;44(8):3189. [26] Schmidt I, Gross D. Proc Roy Soc Lond A 1999;455:3085. [27] Fahrmann M, Fahrmann E, Paris O, Fratzl P, Pollock TM. In: Superalloys 1996. Champion, PA. Warrendale, PA: Minerals, Metals & Materials Society; 1996. [28] Buffiere JY, Ignat I. Acta Metall Mater 1995;43:1791. [29] Ve´ron M, Bre´chet Y, Louchet F. Scripta Mater 1996;34:1883. [30] Rodney D, Le Bouar Y, Finel A. Acta Mater 2003;51:17. [31] Zhou N, Shen C, Mills MJ, Wang Y. Acta Mater 2007;55:5369. [32] Zhou N, Shen C, Mills MJ, Wang Y. Acta Mater 2008;56:6156. [33] Gaubert A, Finel A, Le Bouar Y, Boussinot G. In: Jeulin D, Forest S, editors. Continuum models and discrete systems (CMDS 11). Paris: Les Presses de l’Ecole des Mines; 2007. p. 161. [34] Gaubert A, Le Bouar Y, Finel A. Phil Mag B 2010;90:375.