Available online at www.sciencedirect.

com

Acta Materialia 60 (2012) 5311–5321

www.elsevier.com/locate/actamat

Phase-field modeling of precipitate evolution dynamics in

elastically inhomogeneous low-symmetry systems: Application
to hydride precipitation in Zr
L. Thuinet ⇑, A. De Backer, A. Legris
Université de Sciences et Technologies de Lille, Unité Matériaux Et Transformations (UMET), CNRS UMR 8207, Université Lille 1,
ENSCL, F-59655 Villeneuve d’Ascq, France

Received 10 April 2012; received in revised form 30 May 2012; accepted 31 May 2012
Available online 24 July 2012

Abstract

A phase-field model was developed within the framework of heterogeneous elasticity theory to study the precipitation of particles with
trigonal symmetry in a hexagonal matrix. The model is first calibrated and successfully compared with previous analytical calculations
performed to explain the effect of symmetry-breaking transformations on precipitate morphology. Secondly, the model was adapted to
study the precipitation of the coherent f hydride phase in zirconium. The results are consistent with the well-established experimental
observation of the existence of acicular precipitates aligned along the dense directions in the basal plane. Moreover, original kinetic path-
ways are implied by the presence of a threefold axis of symmetry, leading to the emergence of original morphological bifurcations not
previously reported and probably related to the inconsistency between the threefold symmetry and the inversion properties of the B func-
tion introduced by Khachaturyan. In spite of its simplicity (only one order parameter is taken into account), the present phase-field
model gives rise to very complex morphological sequences.
Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

Keywords: Precipitation; Phase-field models; Hydrides

1. Introduction easily incorporates strain-induced effects during the nucle-

ation, growth and coarsening. Most particularly, this the-
The equilibrium particle shape during the early stages of ory establishes a very useful criterion to determine the
precipitation is the result of competing effects of the inter- habit plane of well-developed coherent precipitates assum-
facial and elastic energies giving rise to a high diversity of ing that the interfacial energy plays a negligible part.
morphologies. Moreover, the kinetic path followed by the According to this theory, the normal to equilibrium habit
precipitates can go through different transient shapes plane n0 is the direction that minimizes the following func-
before reaching the equilibrium shape. For example, a wide tion B:
range of c0 morphologies (doublets, octets or chaplets of
BðnÞ ¼ r0ij e0ij  ni r0ij Xjm ðnÞr0mn nn ð1Þ
cuboids) are reported in the literature concerning nickel-
based superalloys, some of them clearly induced by kinetic This function includes all the elastic parameters of the sys-
effects, as shown by numerical investigations using finite- tem: the stress-free strain tensor e0kl , which is associated
element [1] or phase-field [2] calculations. with the stress-free volume and shape change of the matrix
Equilibrium shapes can be investigated using the micro- due to the transformation, r0ij ¼ C ijkl e0kl with the elastic con-
elasticity theory developed by Khachaturyan [3], which stants of the system assumed to be elastically homogeneous
(the matrix and precipitate have the same elastic con-
⇑ Corresponding author. stants), and the inverse Green function tensor X1 jm ðnÞ ¼
E-mail address: Ludovic.Thuinet@univ-lille1.fr (L. Thuinet). C jklm nk nl where ni are the Cartesian components of n. This

1359-6454/$36.00 Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.actamat.2012.05.041
5312 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321

criterion remains valid even in the case of inhomogeneous inside the hexagonal matrix. In order to determine whether
elasticity conditions (distinct elastic constants between the such variants can be elastically induced, Eshelby’s equiva-
matrix and the precipitate), by considering only the elastic lent inclusion method was adopted in a previous study
constants of the precipitate. Note that all the crystallo- [6], allowing us to conclude, under the heterogeneous elas-
graphic symmetry information is included in function B. ticity assumption, that the most stable shape was elliptic in
Despite the huge impact of this criterion, it does not the basal plane, with the preferential orientation of the tri-
incorporate the interfacial effect which may induce mor- gonal precipitates along the dense directions, this orienta-
phological bifurcations. Different approaches, classified in tion being dictated by the extra elastic constant C15 of
Ref. [4], have therefore been proposed to study the mor- the precipitate.
phological evolution of a misfitting particle inside a matrix Although useful, this approach is based on the assump-
under the combined elastic and interfacial effects. The first tion that the precipitate is ellipsoidal, though no argument
category consists in restricting the possible precipitate can ascertain whether it is a good approximation of the ex-
shapes to a particular geometrical family (for example, act equilibrium shape. Moreover, this method does not al-
ellipses in two dimensions and ellipsoids in three dimen- low an exhaustive approach of all the possible transient or
sions). It offers the advantage of giving analytical solutions equilibrium morphologies that can be encountered in this
to morphological bifurcation problems, precisely determin- type of system. Indeed, by analogy with homogeneous cu-
ing the influence of the various material parameters [5,6]. bic systems with purely dilatational misfit, it is expected
However, some possible shapes may be ignored by this that two kinds of shape transition will be observed in the
method, which is its main drawback. Other numerical ap- system (Fig. 1): symmetry-conserving (resp. -breaking)
proaches (called nonrestricted methods, or NRM) allow shape transition, for which the symmetry of the precipitate
no restriction on transient or equilibrium particle morphol- shape is equal to or greater (resp. less) than the symmetry
ogies [7–16]. Moreover, they can simulate shape evolution resulting from the intersection of the point group of the
dynamics during growth and coarsening. On the latter as- precipitate and matrix phases. By definition, the circle–
pects, inhomogeneous elastic constants in the system may ellipse bifurcation occurring in the basal plane and studied
have a significant effect on shape instabilities. For example, in Ref. [6] is a symmetry-breaking transition. Indeed,
the application of the discrete atom method [10] to an ini- applying the Curie principle, the resulting shape of a sym-
tially soft circular isotropic inclusion in an isotropic matrix metry-conserving transition
 should contain  at least a three-
(lp = lm/2, with lp and lm the shear moduli of the precip- fold axis of symmetry m6 m2 m2 \ 3 m2 ¼ 3 m2 , which is not the
itate and matrix, respectively) reveals that it is unstable and case of an ellipse with distinct semi-axes. It was established
develops topological waves at its surface before reaching its (under assumptions specified in Ref. [6]) that this transition
equilibrium elliptical shape. The corresponding simulation occurs if L > 4 (Fig. 1), where L is given by:
was achieved in Ref. [9] with a phase-field model (PFM)
taking into account elastic inhomogeneities. This feature
was incorporated via the algorithm presented in Ref. [17]
and based on the perturbation method of Ref. [18]. Both
methods, although developed at two distinct scales, give
the same evolution of the precipitate shape, meaning that
the shape destabilization mechanism is not model depen-
dent but has a true physical origin. Moreover, despite the
simple isotropy assumption of the elastic constants in both
matrix and precipitate, these studies show that elastic inho-
mogeneities can induce significant shape instabilities.
A thorough survey of the literature shows that there is a
lack of information on systems with lower symmetry than
cubic. This paper proposes an extension of the study of
shape instabilities, leading to the formation of transient
or equilibrium shapes, to the case of hexagonal systems
in inhomogeneous elasticity conditions. Furthermore, to
illustrate the possible impacts of the study on technological
alloys, an example of great industrial interest is chosen as
the reference case: the precipitation of f hydrides in an
aZr matrix [6,19,20], for which a crystallographic symme-
try break occurs, since f hydrides belong to the trigonal
point group 3m whereas the aZr matrix belongs to the
Fig. 1. Diagram representing the 2-D bifurcation of (a) a misfitting
hexagonal point group 6/mmm. The formation of three dif- trigonal precipitate inside the basal plane of a hexagonal matrix
ferent orientation variants is experimentally observed in the (inhomogeneous elasticity) and (b) a particle in a homogeneous cubic
basal plane during precipitation of the trigonal precipitates system with a purely dilatational misfit.
 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321 5313

e2 Crp F ¼ F chem þ Eelas ð4Þ

L¼ ð2Þ
c Fchem is given by an integral over the material volume Vm:
Z Z Z  
where e is the stress-free strain (isotropic in the basal b 2
plane), rp is the radius of the precipitate, c is the interfacial F chem ¼ fhom ðcÞ þ ðrcÞ dV ð5Þ
Vm 2
energy (assumed to be isotropic) and C is an elastic param-
eter equal to ½3C 
C2
12 þ2C44 2 15
. It must be emphasized that this In practice, the continuous integral is replaced with a
C12 þ2C44 C44
criterion is formally very close to the one proposed in Ref. discrete sum, the order parameters being evaluated at the
[13] for the symmetry-breaking bifurcation from fourfold points of a grid with spacing a. fhom(c) is the homogeneous
symmetric shapes to twofold symmetric shapes in homoge- incoherent free energy density represented by a double-well
neous cubic systems with purely dilatational misfit: polynomial:
 
1 2 1 4
e2 C44 rp fhom ðcÞ ¼ A  ðcðr; tÞ  cÞ þ 2 ðcðr;tÞ  cÞ ¼ Af~ hom ðcÞ
L¼ > 5:6 ð3Þ 2 Dc
c
ð6Þ
The difference is that criterion (2) was analytically estab-
lished, whereas criterion (3) was numerically established A is the free energy scale parameter of the model and
using an NRM. The work in Ref. [6] clearly established f~ hom ðcÞ is the nondimensional homogeneous incoherent
that Eshelby’s equivalent inclusion method is successful free-energy density, c ¼ c1 þc
2
2
and Dc ¼ c2  c1 , where c1
for studying the symmetry-breaking circle–ellipse bifurca- and c2 correspond to the two minima of the double-well
tion in our reference system, but a transient or equilibrium potential, i.e. the noncoherent equilibrium concentrations
shape resulting from a symmetry-conserving transition of the two phases (c1 < c2).
2
cannot be taken into account by this shape-restricting The heterogeneous part fhet ðcÞ ¼ b2 ðrcÞ is strictly posi-
approach. tive and contributes to F only in regions where c(r,t) is sub-
This preliminary comparison justifies the use of an jected to spatial variations, i.e. in the vicinity of interfaces.
NRM to obtain a general overview of all the shape transi- b is called the stiffness parameter and is linked to the inter-
tions in the system. To this end, the PFM was chosen in facial energy. Due to the particular form of fhom (Eq. (6)), it
this paper to investigate in detail the possible transition can be established that the interfacial energy (c) and width
shapes in an inhomogeneous system with a crystallographic of the interface (e) are given by:
0:5
symmetry-break precipitation. Indeed, as mentioned Dc2 ðAb=2Þ
above, the reliability of this method, even with inhomoge- c¼ ð7Þ
3
neous elastic constants, has been successfully tested. 0:5
and e ¼ ne a ¼ 2ð2b=AÞ ð8Þ
The paper is organized as follows. The first part of this
paper briefly recalls the basic assumptions, formulation with ne the number of grid points at the interface.
and calibration of the PFM with inhomogeneous elasticity. From Eqs. (7) and (8), it can be established that:
The second part then tests its reliability in the case of a AeDc2 Ane aDc2
crystallographic symmetry-break precipitation by compar- c¼ ¼ ð9Þ
12 12
ison with previous analytical results. Finally, in the last In this approach, the interfacial energy is isotropic.
part, new results concerning shape instabilities and shape The elastic energy in the absence of any applied stress is:
bifurcations in the trigonal/hexagonal system will be char- Z Z Z
acterized, discussed and compared to previous results Eelas ¼ C ijkl ðrÞðeij ðrÞ  e0ij ðrÞÞðekl ðrÞ  e0kl ðrÞÞdV
reported in the literature. Vm

eij(r) is the local strain in the system and e0ij ðrÞ ¼ e00 ij
2. Phase field model ðcðrÞ  c0 Þ, where e00 ij is the tensor of the Vegard coefficients
and c0 is the average composition. To calculate the elastic
Generally, in a PFM, two different types of order strain energy of the inhomogeneous system by the pertur-
parameters are considered: one that has to verify conserva- bation method, the difficulty rests in the determination of
tion laws (conserved parameter), which is the case of the the displacement field uk ðrÞ that obeys the following equa-
molar fraction or composition noted c(r,t), and one that tion [18,17]:
is a nonconserved parameter noted g(r,t), which is related   
to the crystal structure of the solid solution and the precip- @2 @ @
C ijkl þ DC ijkl DcðrÞ uk ðrÞ
itates. In the case considered here, precipitation of only one @rj @rl @rj @rl
type of eigenstrain variant is modeled, and a single field @DcðrÞ @Dc2 ðrÞ
corresponding to the local composition c(r,t) is sufficient ¼ ðC ijkl e00
kl  DC ijkl
ekl Þ þ DC ijkl e00
kl ð10Þ
@rj @rj
to determine the free energy F of the evolving system. F
has two main contributions, which respectively depend with C ijkl ¼ cc02 c 1
C 2 þ cc22 c
c1 ijkl
0
C 1 , DC ijkl ¼ c2 c
c1 ijkl
1
1
ðC 2ijkl  C 1ijkl Þ,
on the local chemistry (Fchem) and nonlocal elastic effects u
Dc(r) = c(r)  c0 and C ijkl are the elastic constants of phase
(Eelas): u.
5314 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321

eij is the homogeneous strain: since the case of an applied 3. Results

stress is not discussed here, eij is taken to be zero due to the
choice of the reference state. The treatment of the nonlinear In this section, we present the results obtained with the
term on the left-hand side of the equation (perturbation PFM. The first part deals with generic cases in order to
term) requires the use of an iterative approach in the Fourier compare the PFM results with previous analytical results.
space to solve the elastic displacement field [17]. The In the second part, the PFM was parameterized to investi-
displacement u is first initialized to the value u(0), which gate the precipitation of the recently characterized coherent
corresponds to the solution of the homogeneous system: hydride phase f.
@2 @DcðrÞ
C ijkl u0k ðrÞ ¼ C ijkl e00
kl ð11Þ 3.1. Comparison of the PFM with analytical results
@rj @rl @rj
The solution is then determined iteratively using the fol- In order to validate the PFM for the system under con-
lowing relation: sideration in this paper, its numerical results are compared
@2 @DcðrÞ to the analytical relation established in Ref. [6] and recalled
C ijkl unþ1 ðrÞ ¼ ðC ijkl e00
kl  DC ijkl
ekl Þ in Eq. (2). First, a numerical study is performed to deter-
@rj @rl k @rj
2
  mine the required number of iterations n in the perturbation
00 @Dc ðrÞ @ @
method presented in Section 2 to ensure the convergence in
þ DC ijkl ekl þ DC ijkl DcðrÞ un ðrÞ ð12Þ
@rj @rj @r1 k displacement by solving Eqs. (11) and (12). The calculation
The evolution of the concentration field is governed by is carried out for an initial circular precipitate of radius
the Cahn–Hilliard equation [21]: rp = 100a placed in the center of a two-dimensional (2-D)
simulation domain of 512a  512a. The tensor of the
@cðr; tÞ dF
¼ Mr2 ð13Þ Vegard coefficient is e00 ij ¼ 0:025dij . The values of the elastic
@t dcðr; tÞ constants adopted for this calculation are reported in
For the sake of numerical efficiency, the evolution equa- Table 1 (they correspond to the experimental or numerical
tion is solved in Fourier space: values of the a-f system, as calculated in Ref. [25]). In the
paper, all the elastic constants are expressed in the basis
@cðq; tÞ dF
¼ Mq2 ðq; tÞ ð14Þ ([1 0 10],[1210],[0 0 0 1]). Fig. 2 shows that an accurate
@t dcðq; tÞ solution is reached for n = 2. This result is probably due
where F(q,t) is the Fourier transform of F(r,t). q is the re- to the small differences in the elastic constants between the
ciprocal space vector and is limited to the first Brillouin matrix and the precipitate, which is at most 23 GPa. So,
zone. The calculation domain is defined by Li and Ni, for all the following simulations reported in this paper,
which are respectively the size and number of cells along n = 3 will be adopted.
the ith Cartesian direction. Because of the periodic bound- After this convergence study, we consider the shape evo-
ary conditions applied to the simulation domain, each lution of circular precipitates under the effect of a trigonal
Cartesian coordinate qi of q along the ith direction is equal perturbation inside a homogeneous isotropic system. The
to 2p Lnii , with –Ni/2 < ni 6 Ni/2 (with Ni even and ni2Z). Eq. initial precipitate radius rp = 43a is adopted, the simulation
(14) is a nondimensionalized by introducing the following domain being 256a  256a. c1 and c2 are arbitrarily set at
reduced quantities (noted by a tilde): values of 0.1 and 0.9. Since a is set at 0.5 nm, rp = 21 nm.
An interface energy of c = 10 mJ m–2 is adopted, which
q ¼ aq
~
implies a value of A calculated from Eq. (9):
~t ¼ t=t0 with t0 ¼ a2 =MA A = 11.43  107 J m–3. The system is chosen to be elasti-
e elas ¼ Eelas =A cally homogeneous and isotropic except in the precipitate,
f~ hom ¼ fhom =Aðsee Eq: ð6ÞÞ and E ð15Þ
for which a trigonal perturbation is introduced through a
~ ¼ b=Aa2
b nonzero value of C15. The elastic constants are specified in
Eq. (14) can then be expressed as follows: Table 2 and correspond to an isotropic approximation of
" # the real elastic constants of pure zirconium. The tensor of
@cð~q;~tÞ ~ ~ elas
dE
2 @ f hom ~ 2
¼ ~
q ð~q;~tÞ þ b~ q;~tÞ þ
q cð~ q;~tÞ
ð~ ð16Þ
@~t @c dc Table 1
Values (in GPa) of the elastic constants of a and f.
Eq. (16) is solved using the semi-implicit algorithm
Elastic constants Matrix Precipitate
described in Ref. [22]. By combining Eq. (8) with the
expression of b~ (Eq. (15)), it is easy to show that b ~ ¼ n2e . c11 155 168
8
c12 67 89
Due to interface mobility problems, three nodes must be c13 65 67
used at least to describe the interface (ne = 3). Conse- c33 173 195
~ is set at 1.5. The values of the nor-
quently, the value of b c44 36 29
malizing parameters will be specified for each example of c66 44 39.5
the following sections. c15 0 23
 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321 5315

Fig. 2. (a) The xx component and (b) yy component of the elastic stress along the x direction across the center of a trigonal precipitate obtained for
different iteration numbers of the perturbation method.

the Vegard coefficients e00 ij was arbitrarily set at 0.01dij. C crit

15 = ± 23 GPa. As a conclusion, the critical value of C15
Fig. 3 represents the morphological evolution of the initial deduced from the phase-field simulations agrees with the
precipitate for three different values adopted for C15 = 24, analytical simulation quite well. This comparison turns out
26 and 28 GPa. For C15 = 26 and 28 GPa, a circle–ellipse to be a good quantitative validation of the PFM applied to a
bifurcation is observed, with the ellipse corresponding to crystallographic symmetry-breaking phase transformation.
the equilibrium shape of the precipitate. The ellipse is The greater the value of C15, the fewer the number of
aligned along one dense direction of the system. As a result, time steps required to detect the occurrence of the bifurca-
the analytical approach reported in Ref. [6] and the phase- tion. If C15 is close to C crit
15 , the number of time steps
field simulations give the same qualitative results. required to detect a bifurcation is very high, which makes
The circle–ellipse bifurcation is observed for C15 = 26 the determination of C crit
15 difficult to achieve by the PFM.
and 28 GPa (Fig. 3b and c), but not for C15 = 24 GPa Surprisingly, the PFM results reported in Fig. 3 show
(Fig. 3a). This means that the critical value C crit15 of C15 that, whatever the value of C15, no symmetry-conserving
below which no bifurcation is observed is between 24 and shape transition is detected in the different simulations,
26 GPa. For a given value of rp, it is possible to analytically which makes this system qualitatively different from the
determine C crit
15 from Eq. (2): homogeneous cubic system with a dilatational misfit shown
2 31=2 in Fig. 1.
4cC44 The simulation in Fig. 4 was obtained with the same
6 7
C crit
15 ¼ 4 2 5 ð17Þ parameters apart from the elastic constants, which are
3C12 þ2C44
C12 þ2C44
e2 r p reported in Table 1, in order to study the effect of a trigonal
perturbation inside an inhomogeneous anisotropic system.
The stress-free strain between both phases is A transient morphology is observed in Fig. 4 for ~t = 27,500
e0ij = 0.01dij(c2 – c1) = 0.8%, but this value must be slightly and 275,000, with the development of two branches corre-
corrected. Indeed, the coherency strain energy induces a sponding to two different orientation variants among the
decrease in the two-phase stability domain and a Dccoh/incoh three possible variants, giving birth to a “bone” morphol-
shift of the equilibrium concentrations. With a two-phase ogy. Then there is the selection of only one variant during
elastically isotropic and homogeneous system, it can be coarsening. The equilibrium morphology is still elliptical,
shown [23,24] that Dccoh/incoh = B(c2 – c1)/(2A + 2B) with but the variant selected during the simulation is different
B = 2l[e00]2(1 + m)/(1  m), l and m being, respectively, the from that in Fig. 3. It was suggested in Ref. [6] that the hex-
system’s shear modulus and Poisson’s ratio. Consequently, agonal heterogeneities of the Voigt matrix promote the
ccoh coh
1 = c1 + Dccoh/incoh = 0.15 and c2 = c2 – Dccoh/incoh = deformation from the circle to an ellipse, which is con-
0.85, inducing a new value of e0ij = 0.01dij firmed by the simulation reported in Fig. 4: for the same
(ccoh coh
2  c1 ) = 0.7%. By integrating this result in Eq. (17), value of C15 (24 GPa), no bifurcation was observed in the
isotropic approximation.
Despite the emergence of new morphologies, none of
Table 2
them exhibits the threefold symmetry characteristics of a
Values (in GPa) of the elastic constants in the isotropic
approximation. symmetry-conserving transition. Indeed, since the bone
morphology has a twofold symmetry, the corresponding
c11 157.5
c12 67.5 shape bifurcation is symmetry breaking.
c13 67.5 The next section aims at establishing whether the previ-
c33 157.5 ous conclusions concerning shape transitions of an isolated
c44 45 precipitate are still valid for a population of interacting
c66 45
precipitates.
5316 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321

Fig. 3. Influence of a trigonal perturbation inside an isotropic homogeneous system on the morphological evolution of the precipitate for (a)
C15 = 24 GPa, (b) C15 = 26 GPa and (c) C15 = 28 GPa. The values of the elastic constants are reported in Table 2.

Fig. 4. Evolution of a trigonal precipitate morphology inside a close-packed hexagonal matrix for C15 = 24 GPa. The values of the other elastic constants
are reported in Table 1.

3.2. Studies of the morphological evolution dynamics of f composition c0 = 4 at.%. The corresponding equilibrium
zirconium hydrides in the basal plane volume fraction of hydride at the end of the simulation is
5%. As shown in Fig. 6, the hydride volume fraction
All the simulations in this section are performed with the increases rapidly at the beginning of the simulation, which
following parameters. Since the reference case under study is typical of the growth stage: as the supersaturation progres-
corresponds to the a–f system, c2 = 0.33 (due to the partic- sively decreases inside the matrix due to the growth of the
ular stoichiometry of f hydrides) and c1 is taken to be equal precipitates, the growth speed decreases as well until the
to the solubility limit css of hydrogen in an a solid solution. supersaturation completely vanishes. The hydride volume
css(T) is calculated from the following experimental rela- fraction then reaches a second stage that looks very similar
tion [26]: to a plateau (very slow linear evolution with time), which
is characteristic of the coarsening stage. A dotted line has
C ss ðT Þ ¼ 16:075 expð4562=T Þ ð18Þ been added in Fig. 6 as a guide to show the coarsening stage.
T is given in kelvin. At T = 600 K, css = 0.89%. e00ij =
The point at which the curve representing the hydride vol-
0.025dij and the elastic constants are given in Table 1. ume fraction evolution differs from this line graphically dif-
The grid spacing is fixed at a = 1 nm, and the interface ferentiates these two stages, as shown by the vertical lines in
energy is 20 mJ m–2. The domain is meshed in the basal Fig. 6: for c0 = 4 at.%, the growth stage lasts until approxi-
plane of zirconium with 1024  1024 cells, so the size of mately ~t = 30,000. During this stage, there is no significant
the domain is approximately 1 lm. shape evolution of the precipitates: they remain almost cir-
The simulation represented in Fig. 5a was obtained for 20 cular except for one, which becomes elliptical (top-left-hand
nuclei randomly distributed inside the matrix and an initial corner of Fig. 5a) for ~t = 5000 and ~t = 20,000. The coarsening
 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321 5317

Fig. 5. Morphological evolution of interacting precipitates for an initial composition of (a) c0 = 4 at.% (corresponding final volume fraction of hydride of
5%) and (b) c0 = 5 at.% (final volume fraction of hydride of 8%).

as the one applied for c0 = 4 at.%, Fig. 6 shows that the

growth stage lasts until approximately ~t = 50,000. During
the growth stage and contrary to the previous simulation,
the precipitates undergo a circle–triangle transition, reveal-
ing the existence of two distinct variants of triangular pre-
cipitates (Fig. 5b) for ~t = 20,000. As shown in Fig. 7a, the
normal to the faces of the triangular shape correspond to
the h1 0 1 0i directions: they are in fact the minima of the
B function defined by Eq. (1) and calculated with the elastic
constants of the trigonal precipitate, so they will be referred
to as the elastically soft directions of the system. Fig. 7b
shows that the circle–triangle bifurcation does not always
Fig. 6. Evolution of volume fraction of hydride for c0 = 4 at.% (red empty
circles) and c0 = 5 at.% (green crosses). The coarsening stages for c0 = 4 occur for the same volume of precipitate but depends on
at.% and c0 = 5 at.% are represented by the red and green dotted lines the local configuration of the neighboring precipitates.
respectively. The vertical red and green solid lines separate the growth and Moreover, single triangular precipitates are rarely isolated,
coarsening stages for c0 = 4 at.% and c0 = 5 at.% respectively. (For but interact with one (or several) other triangular precipi-
interpretation of the references to colour in this figure legend, the reader is
tate(s). The elastic interactions between precipitates pro-
referred to the web version of this article.)
mote the “spike–spike” configuration represented in
Fig. 8a, but destabilize the opposite “face–face” configura-
stage proceeds mainly by coagulation and the disappearance tion, as shown in Fig. 8b: the right precipitate develops a
of smaller precipitates. It induces an increase in the average spike towards the face of the other one, then the left precip-
size of the remaining precipitates, which gradually follow a itate develops its own spike towards the right precipitate.
circle–ellipse transition. This bifurcation mainly occurs dur- The last stage of this sequence shows the coagulation of
ing the coarsening stage, not the growth stage. At both precipitates.
~t = 200,000, three orientation variants of the elliptic precip- During coarsening, the spike–spike configuration is fol-
itates are present in the simulation domain and aligned along lowed by a coagulation stage that gives birth to a bone
the dense directions of the basal plane. It is worth noting that morphology characterized by the development of two
they all correspond to the same eigenstrain variant, since branches, quite similar to the one observed in Fig. 4. Only
only one order parameter was introduced in the PFM (Sec- one of these two branches survives to give the final elliptical
tion 2). The shape evolutions observed in this simulation shape. Whereas the circle–triangle bifurcation occurs dur-
are mainly symmetry breaking and similar to the one pre- ing the growth stage, the triangle–ellipse bifurcation occurs
sented in Fig. 3. during the coarsening stage.
The simulation of Fig. 5b is identical to the previous one This simulation shows that the system can undergo
except that c0 = 5 at.%. The corresponding equilibrium another shape transition, the circle–triangle bifurcation,
volume hydride fraction is 8%. Using the same procedure which is symmetry conserving. This seems a priori to be
5318 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321

and the radius of the circular precipitate is 200 nm: in this

simulation, the initial size of both precipitates is much
larger than the average size of the precipitates observed
in Fig. 5 at the end of the growth stage.
Whereas the elliptical precipitate grows but remains
elliptical, the circular precipitate follows a much more com-
plex morphological evolution, which can be divided into
the following steps:

– The initial circular shape evolves towards a hexagonal

shape (sixfold symmetry). Whereas all the morphologies
obtained previously (ellipses and triangles) exhibit faces,
the normals of which correspond to the elastically soft
direction of the trigonal system, the exact opposite
applies for this hexagonal shape, i.e. the normal to the
hexagonal faces are along hard directions.
– The faces of the hexagonal shape become concave, and
the resulting new morphology looks like a star, i.e. a six-
fold symmetric shape with 12 faces that respect the sys-
tem’s elastically soft directions.
Fig. 7. (a) The normal to the faces of the triangles correspond to the – One summit of the star shape develops and the others
system’s elastically soft directions. (b) The circle–triangle transition radius shrink, giving rise to the elliptical final morphology
strongly depends on the configuration of the precipitates.
(twofold symmetric shape).

the counterpart of the circle–square transition in cubic

systems; however, the resulting morphology is not an equi- As a result, if the precipitate is allowed to reach a suffi-
librium morphology, since it always evolves towards the cient size during the growth or coarsening stage, it can go
elliptical shape. through more complex shapes than the triangular shape
Since the elastic effects are amplified with the size of the before reaching the equilibrium elliptical shape.
precipitate, other types of shape transitions may be
observed during the evolution of larger inclusions. In order 4. Discussion
to study this effect in our system, two distinct precipitates,
with a different initial morphology (circular and elliptical in The analytical study presented in Ref. [6] and based on
the basal plane), are put in the matrix with a low supersat- the Eshelby’s equivalent inclusion method establishes that
uration of approximately 1 at.%. The two main axes of the the extra elastic constant C15, typical of the trigonal
elliptical precipitate are, respectively, 400 nm  100 nm, point-group system to which f hydrides belong, causes the

Fig. 8. (a) Illustration of the spike–spike configuration. (b) Instability of the face–face configuration.
 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321 5319

circle–ellipse bifurcation observed in the basal plane. The ants is related to the eigenstrain tensors that enter in the
demonstration assumed that the precipitate is ellipsoidal, mechanical Eqs. (11) and (12) at the zeroth order, whereas
which is quite restrictive. Indeed, at this point, no argument the elastic inhomogeneities have contributions only from
allows us to exclude other types of geometrical forms as the first order. Therefore, for crystallographic reasons,
possible equilibrium shapes of f hydrides. The phase-field modeling the fully coherent f hydride and c requires the
model, since it does not rely on any a priori assumption con- use of a different set of order parameters: since there is only
cerning the equilibrium morphology, makes it possible to one eigenstrain variant for f, one order parameter can
obtain a more general result and still confirms that the describe the microstructural evolution provided that heter-
shape of the precipitate minimizing the sum of the elastic ogeneous elasticity assumptions are taken into account,
and interfacial energies is an ellipse aligned with the dense while homogeneous elasticity is sufficient to model c, but
direction of the matrix, emphasizing the significance of three different order parameters are needed. One drawback
the symmetry break occurring during the transformation. of the latter approach resides in the full coherency hypoth-
It is important to emphasize that the alignment of the pre- esis for the c phase, which can be questioned from experi-
cipitates with dense directions of the basal plane is fully con- mental facts (observed misfit dislocations around the c
sistent with the experimental observations and originates precipitates [30,31]).
from the elastic heterogeneities induced by the symmetry Considering systems containing threefold symmetry
break rather than by the existence of three distinct eigen- axes leads to the emergence of new precipitate shapes
strain variants, as suggested in Refs. [27,28]. Indeed, in which are difficult to predict using shape-restricted meth-
the present model, the entire morphological evolution is ods. Indeed, the phase-field results allow one to exhibit
described by means of only one order parameter, which the circle–triangle symmetry-conserving shape transition
implies that only one eigenstrain variant is considered. in the hexagonal/trigonal system, which cannot be simu-
The three orientation variants, a product of the same eigen- lated by the analytical approach developed in Ref. [6]. At
strain, result from the hexagonal/trigonal symmetry break first sight, the appearance of triangles is somewhat at odds
of the transformation. Furthermore, both approaches (the with the B function criterion, which favors morphologies
phase field and the Eshelby’s equivalent inclusion method) with twofold symmetry properties (from Eq. (1), B(n) =
provide approximately the same value for rc, which proves B(-n)). This last symmetry argument may explain why the
that the algorithm proposed in Ref. [17] combined with a particular circle–triangle transition was not detected in
phase-field model is perfectly well adapted to the study of the simulation illustrated in Fig. 3, for which only a single
shape bifurcations induced by elasticity, even in the case isolated precipitate was considered, while it is observed in
of symmetry-breaking transformation. the simulations described in Fig. 5b. Indeed, in the latter
In previous studies [27,28], the precipitation of the case, the triangular precipitates may interact and thus
tetragonal c hydride phase in zirconium was modeled in adopt a collective configuration, like the spike–spike con-
the context of the homogeneous elasticity theory assuming figuration, recovering the twofold symmetry. It must be
full coherency between the matrix and the precipitate. In noted that the contradictions with the twofold symmetry
this case, the phase transformation implies a symmetry properties associated with function B have not been
change from a hexagonal to a tetragonal point group. If reported in the literature, probably because most of the sys-
the sixfold symmetry axis of the hexagonal group coincided tems studied to date exhibited transitions between the cubic
with the fourfold symmetry axis of the tetragonal group, the or tetragonal symmetry, which is consistent with the sym-
eigenstrain tensor would have been isotropic in the basal metry properties of function B.
plane of zirconium, as in the present case, and in principle The circle–triangle bifurcation was observed for c0 = 5
this transformation could have been studied using the meth- at.% and not for c0 = 4 at.%. The reason for such distinct
odology presented here. However, the crystallographic ori- behaviors is not clearly understood. Nevertheless, the fol-
entation (1 1 1)c//(0 0 0 1)Zr and [ 1
1 2 0]c//[1 1 
2 0]Zr implies lowing arguments can be proposed: when the particle
the existence of three anisotropic eigenstrain tensors [29] grows, the anisotropic strain energy contributes more sig-
of the form e011 = 0.551%, e022 = 5.70% and e033 =5.70%, the nificantly and the circular shape, which minimizes the iso-
other e0ij components being equal to 0 on the tropic interfacial energy, gradually transforms to
([1 1 
0],[1 
1 0 0],[0 0 0 1]) hexagonal basis. The other two vari- polyhedron to relax this strain energy, the facets of the
ants are obtained by writing the same tensor in the previous polyhedron usually being the elastically soft planes whose
basis rotated by 120° and 240° around the c axis. In this orientations minimize function B. In both simulations,
case, it is necessary to introduce three nonconserved order the same number of initial nuclei were introduced in the
parameters, each of them associated with a particular eigen- simulation domain, so, at the end of growth, the average
strain. From this point of view, the situation is similar to a radius of the remaining precipitates is larger for c0 = 5
cubic–tetragonal transformation described by three eigen- at.% than for c0 = 4 at.%. In this case, it can be assumed
strain tensors corresponding to the three orientation vari- that the precipitates do not reach the critical size allowing
ants [1 0 0], [0 1 0] and [0 0 1]. For these transformations, the strain energy to control the morphology. As a conse-
most of the microstructural features are grasped using the quence, they remain almost circular, whereas the opposite
homogeneous elasticity theory since the origin of the vari- is true for c0 = 5 at.%. However, this cannot be the only
5320 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321

explanation: indeed, even for c0 = 4 at.%, the remaining in a homogeneous cubic material with isotropic misfit: for
precipitates grow significantly during coarsening, clearly this particular system, shape instabilities start to develop at
allowing them to reach the critical size for the circle–trian- the corners of a square precipitate and its facets become
gle transition. Instead, they evolve directly to the elliptical concave. Common to the case depicted in Fig. 9 and the
equilibrium shape. The circular-triangle shape occurs only example studied in Ref. [2], these kinetic transitions result
during the growth stage, which means sufficient supersatu- in concave facets, which are schematically represented in
ration is needed for its development. To summarize, for Fig. 10b and e. At this stage, it must be emphasized that
sufficiently high supersaturation, a double bifurcation cir- there is a significant difference between the two cases.
cle–triangle, triangle–ellipse is observed, while a direct cir- Before the development of the concavity, the facets of the
cle–ellipse transition is predicted for lower supersaturation. square precipitate in Ref. [2] are oriented along the soft
For larger initial precipitate sizes (radius greater than direction. When the supersaturation vanishes, this kinetic
200 nm), a new circle–hexagon transition is predicted, as effect disappears and the shape recovers the facets which
shown in Fig. 9 for ~t = 30,000. This intermediate morphol- minimize its strain energy, i.e. those it had before shape
ogy infringes the B function rule since, even if the morphol- instability (Fig. 10c). In this case, elastic and kinetic effects
ogy is twofold, the normal to the precipitate facets have opposite consequences on the shape development of
corresponds to the maximum instead of the minimum of the precipitate. In contrast, the facets of the hexagon in
B. Different simulations were attempted to see whether a Fig. 9 are oriented along the hard directions. Since diffu-
circle of the same volume as the one initially entered in sion is enhanced at its corners, each face becomes concave,
the simulation in Fig. 9 could evolve towards a hexagon though the concavity is also promoted by the elastic strain,
oriented with the normal to its facets along the elastically since the resulting faces have their normal coincident with
soft directions (minimum of B), but, surprisingly, no simu- the elastically soft directions, as shown in Fig. 9 for ~t =
lation resulted in this shape. Despite this unusual faceting, 70,000. In this particular case, the concavity is enhanced
it is worth emphasizing that this transition decreases the by both elastic and kinetic effects, and hence there is no
strain energy and is thus permissible. driving force to recover the previous facets of the hexagon.
The hexagon in Fig. 9 evolves towards an ellipse going The transition from a facet oriented along the hard to the
through a star shape. This type of transition may stem soft direction through the concavity development, as
from kinetic effects. Indeed, due to the larger exposure area observed in Fig. 9, is schematically represented in
to the oncoming diffusional flux of solute atoms at the cor- Fig. 10d–f.
ners between the facets of the hexagon, the growth rate is To compare experimental observations with the predic-
locally enhanced and causes a deviation of the interfacial tions of the model, it must be noted that rather unusual
front from these facets, which become concave. This kinetic shapes were always observed for sufficiently large initial
effect has been analyzed in Ref. [2] for a growing precipitate circular precipitates and suitably strong supersaturations.

Fig. 9. Morphological sequence for larger precipitate sizes: initial configuration (~t = 100), intermediate hexagonal shape (~t = 30,000), evolution to a star
(~t = 70,000) and then to an elliptical shape (~t = 100,000 and 300,000).
 L. Thuinet et al. / Acta Materialia 60 (2012) 5311–5321 5321

Fig. 10. Shape instability referenced in Ref. [2]. (a) Facet of the precipitate oriented along the elastically soft direction ns. (b) Evolution of the facet of the
precipitate to a concave shape due to diffusional effect during growth. (c) Recovering of the initial facet during coarsening. Shape instability observed in
Fig. 9. (d) Facet of the precipitate oriented along the elastically hard direction nH. (e) Evolution of the facet of the precipitate to a concave shape due to
diffusional effect during growth. (f) Formation of the facets oriented along the elastically soft directions ns.

Such complicated shapes could be reproduced experimen- Acknowledgements

tally, for instance by quenching a sample that was previ-
ously maintained long enough just below the solubility The authors thank the Centre de Ressources Informa-
limit. tiques of the Université de Sciences et Technologies de Lille
for computational facilities.
5. Conclusions
References
Shape transitions occurring during growth and coarsen-
ing of an elastically inhomogeneous two-phase system con- [1] Hazotte A, Grosdidier T, Denis S. Scripta Mater 1996;34:601.
taining a threefold symmetry axis is studied by a phase-field [2] Wang Y, Khachaturyan AG. Acta Metall Mater 1995;43:1837.
[3] Khachaturyan AG. Theory of structural transformations in sol-
approach without any a priori assumption on the particu- ids. New York: Wiley; 1983.
lar shape of the precipitate. This method confirms the [4] Voorhees PW, Johnson WC. Solid State Phys 2004;59:2.
experimental results obtained for f hydrides in a zirconium [5] Johnson WC, Cahn JW. Acta Metall Mater 1984;32:1925.
alloys: f precipitates are acicular in the basal plane, aligned [6] Thuinet L, Legris A. Acta Mater 2010;58:2250.
[7] Le Bouar Y, Loiseau A, Khachaturyan AG. Acta Mater
along the dense directions. In particular, by comparison
1998;46:2777.
with the analytical results based on Eshelby’s previous [8] Boussinot G, Finel A, Le Bouar Y. Acta Mater 2009;57:921.
modified equivalent inclusion method, it can be concluded [9] Boussinot G, Le Bouar Y, Finel A. Acta Mater 2010;58:4170.
that the influence of the f extra elastic constant C15 on the [10] Lee JK. Scripta Metall Mater 1995;32:559.
bifurcation radius of the precipitate is quantitatively well [11] Lee JK. Mater Sci Eng A 1997;238:1.
reproduced by the phase-field method. Due to its reliability [12] Choy JH, Lee JK. Mater Sci Eng A 2000;285:195.
[13] Thompson ME, Su CS, Voorhees PW. Acta Metall Mater
in the study of crystallographic symmetry-breaking trans- 1994;42:2107.
formations, the phase-field model has been applied to the [14] Thornton K, Akaiwa N, Voorhees PW. Acta Mater 2004;52:1353.
dynamics of shape evolution in systems less symmetric than [15] Thornton K, Akaiwa N, Voorhees PW. Acta Mater 2004;52:1365.
cubic systems. The presence of a threefold symmetry axis [16] Li X, Thornton K, Nie Q, Voorhees PW, Lowengrub JS. Acta Mater
leads to a high diversity of original shape transitions that 2004;52:5829.
[17] Hu SY, Chen LQ. Acta Mater 2001;49:1879.
have not been predicted numerically until now. Besides [18] Khachaturyan AG, Semenovskaya S, Tsakalakos T. Phys Rev B
the simple circle–ellipse bifurcation, transition proceeds 1995;52:15909.
via more complex morphological sequences, including the [19] Zhao Z, Blat-Yrieix M, Morniroli JP, Legris A, Thuinet L, Kihn Y,
circle–triangle–ellipse and circle–hexagon–star–ellipse et al. J ASTM Int 2008;5:101161.
paths. These transitions exhibit some unexpected features, [20] Zhao Z, Morniroli JP, Legris A, Ambard A, Kihn Y, Legras L, et al. J
Microsc 2008;232:410.
namely the development of transient shapes that are not [21] Cahn JW, Hilliard JE. J Chem Phys 1958;28:258.
twofold during growth, or the faceting of transient shapes [22] Chen LQ, Shen J. Comput Phys Commun 1998;108:147.
along elastically hard directions. The present model, [23] Lee JK, Tao W. Acta Metall Mater 1994;42:569.
although containing only one order parameter, gives rise [24] Boussinot G. Thesis, Université de Paris 6; 2007 [in French].
to three orientation variants based on the same eigenstrain. [25] Thuinet L, Besson R. Intermetallics 2012;20:24.
[26] Zuzek E, Abriata JP, San-Martin A, Manchester FD. Bull Alloy
In spite of its simplicity, it reveals a very rich kinetic evolu- Phase Dia 1990;11:385.
tion and allows one to model systems of great technological [27] Ma XQ, Shi SQ, Woo CH, Chen LQ. Comput Mater Sci 2002;23:283.
interest, such as the one covered in this work. By assuming [28] Ma XQ, Shi SQ, Woo CH, Chen LQ. Scripta Mater 2002;47:237.
coherent precipitation, this study suggests a plausible [29] Carpenter GJC. J Nucl Mater 1973;48:264.
[30] Carpenter GJC, Watters JF, Gilbert RW. J Nucl Mater 1973;48:267.
explanation to the long-standing problem of the existence
[31] Carpenter GJC. Acta Metall 1978;26:1225.
of hydride variants in the basal planes of zirconium alloys.