Phase Field Modeling of Fracture and Stress Induced Phase Transitions

0.28 MB·English

by  R. Spatschek

Phase Field Modeling of Fracture and Stress Induced Phase Transitions
R. Spatschek, C. Mu¨ller-Gugenberger, and E. Brener
Institut fu¨r Festk¨orperforschung, Forschungszentrum Ju¨lich, D-52425 Ju¨lich, Germany
B. Nestler
Institute of Computational Engineering, University of Applied Sciences D-76133 Karlsruhe, Germany
7 (Dated: February 6, 2008)
0
0 Wepresentacontinuumtheorytodescribeelasticallyinducedphasetransitionsbetweencoherent
2 solidphases. Inthelimitofvanishingelasticconstantsinoneofthephases,themodelcanbeusedto
describefractureonthebasisofthelatestageoftheAsaro-Tiller-Grinfeldinstability. Startingfrom
n
a sharp interface formulation we derive the elastic equations and the dissipative interface kinetics.
a
J Wedevelopaphasefieldmodeltosimulatetheseprocessesnumerically;inthesharpinterfacelimit,
it reproduces the desired equations of motion and boundary conditions. We perform large scale
5
simulationsoffractureprocessestoeliminatefinite-sizeeffectsandcomparetheresultstoarecently

developed sharp interface method. Details of the numerical simulations are explained, and the
]
i generalization to multiphase simulations is presented.
c
s
- PACSnumbers: 62.20.Mk,46.50.+a,81.40.Np
l
r
t
m
I. INTRODUCTION of motion. But they suffer from inherent finite time sin-
. gularities which do not allow steady state crack growth
t
a unless the tip radius is again limited by the phase field
The propagationof cracks is very important for many
m interfacewidth[10]. Verydifferentapproacheswhichare
applications and a centraltopic in physics and materials
- not based on a phase field as order parameter introduce
science. The most fundamental basis of understanding
d
a tip scale selection by the introduction of complicated
n fracture traces back to Griffith [1]; according to his find-
nonlinear terms in the elastic energy for high strains in
o ings,thegrowthofcracksisdeterminedbyacompetition
the tip region [11], requiring additional parameters.
c of a releaseof elastic energyand a simultaneous increase
[ of the surface energy if a crack extends. Although much Recently, we developed a minimum theory of fracture

1 progresshasbeenmadeinunderstandingthestrikingfea- [12] which is only based on well-established thermody-
v tures of cracks [2], the mechanisms which determine the namical concepts. This is also motivated by experimen-
1 dynamics of crack propagation are still under heavy de- tal results showing that many features of crack growth
0 bate. Atypicaldescriptionofcracksstartsontheatomic are rather generic [13]; among them is the saturation of
1 levelandinterpretsthe propagationby successivebreak- the steady state velocity appreciably below the Rayleigh
1
ingofbonds;itisobviousthatthetheoreticalpredictions speed and a tip splitting for high applied tension. This
0
significantly depend on the underlying empirical models theory describes crack growth as a consequence of the
7
0 of the atomic properties (see for example [3]). A rather Asaro-Tiller-Grinfeld(ATG)instability[14]intheframe-
/ complementary approach takes into account effects like work of a continuum theory. Mass transport at the ex-
t
a plasticity,whichcanleadtoextendedcracktips(finitetip tended crack tips can be either due to surface diffusion
m radiusr )[4]. Recentexperimentalinvestigationsoffrac- or a phase transformation process. The latter has been
0
- ture in brittle gels [5] possibly revealmacroscopicscales. investigated numerically by phase field simulations [15]
d It is obvious that under these circumstances a full mod- and sharp interface methods [16]. It turned out that the
n
elingofcracksshouldnotonlydeterminethecrackspeed phase field simulations were still significantly influenced
o
butalsotheentirecrackshapeandscaleself-consistently. byfinite sizeeffects andinsufficientseparationofthe ap-
c
: During the past years, phase field modeling has pearing length scales, and therefore the results did not
v
coincide. One central aim of the current paper is there-
emerged as a promising approach to analyze fracture
i
X by continuum methods. Recent phase field models go fore to carefully extrapolate new phase field results ob-
tained by large-scalecomputations. As we will show, we
r beyond the microscopic limit of discrete models, and
a then get a very convincing agreement of the approaches.
encompass much of the expected behavior of cracks
Also, we explain details of the phase field method and
[6,7,8,9];However,asignificantfeatureofthesedescrip-
the underlying sharp interface equations in more detail.
tionsisthatthescaleofthegrowingpatternsisalwaysset
by the phase field interface width, which is a purely nu- The paper is organized as follows: First, in Section
mericalparameterandnotdirectlyconnectedtophysical II, we introduce the sharp interface equations to de-
properties; therefore these models do not possess a valid scribe crack propagation as a phase transformation pro-
sharpinterfacelimit. Alternativedescriptions,whichare cess. The basic selection mechanisms for crack growth
intendedtoinvestigatetheinfluenceofelasticstresseson are reviewed in Section III. We introduce a phase field
the morphological deformation of surfaces due to phase descriptiontosolvethearisingmovingboundaryproblem
transitionprocesses,arebasedonmacroscopicequations inSection IV. We demonstrate the numericalseparation
2
of lengthscales and obtain results which are in excellent
agreement with sharp-interface predictions (Section V).
This part contains important new results concerning the V ( 2 )(t)
v
underlying continuum theory of fracture. In Section VI n n
webrieflyexplainhowthemodelcanbeextendedtosys- n=n(1)
(2)
tems consisting of multiple phases. Detailed derivations n
τ
of the sharp interface equations are given in Appendix
A(t)
A. Since coherencyof two solidphases leads to an unex-
pectedexpressionforthechemicalpotential,itsrelevance (1)
V (t)
is analyzed specifically for the motion of a planar inter-
face (Appendix B). In Appendix C it is demonstrated
that the presented phase field model recovers this effect
inthesharpinterfacelimit. Finally,detailsofthenumer-
icalimplementationofthe phase fieldmodel aregivenin
Appendix D. FIG. 1: Geometry of the phase transition scenario. Phase
transitions between phases 1 and 2 are possible and lead to
interface motion with local normal velocity vn. The volumes
of the two phases and the interface A(t) are therefore time-
II. SHARP INTERFACE DESCRIPTION
dependent.
The fracture process in [15] is interpreted as a first
order phase transition from the solid to a “dense gas
phase”, driven by elastic effects. More generally, we Here, the index n denotes the normal direction of the
investigate the transition between two different solid interface, with the perpendicular tangential direction τ
phases. In the limiting case that one phase is infinitely (see Fig. 1). We mention that this equation holds only
soft, crackpropagationcanbe studied. Acentralsimpli- for the specific conditions of equal mass densities and
ficationis due to the assumptionofequalmassdensity ρ coherency at the interface. In other cases one obtains
in both phases and the condition of coherency at the in- moregeneralrelationsformomentumconservationatthe
terface,i.e.the displacementfieldu is continuousacross interface,whichalsoinvolvetheinterfacenormalvelocity
i
the phase boundary, v [2].
n
The elastic contribution to the chemical potential at
u(1) =u(2), (1)
i i the interface for each phase is
where the upper indices indicate the different phases.
1 1
The strain ǫ is related to the displacement field u by µ(α) =Ω σ(α)ǫ(α)− σ(α)ǫ(α)−σ(α)ǫ(α) . (6)
ik i el 2 ττ ττ 2 nn nn nτ nτ
(cid:18) (cid:19)
1 ∂u(α) ∂u(α)
ǫ(α) = i + k . (2) Since the chemical potential has the dimension energy
ik 2 ∂xk ∂xi ! per particle, we introduced the atomic volume Ω. It is
quiteremarkablethatthenormalandshearcontributions
Strainandstressσ areconnectedthroughHooke’slaw;
ik enter into the expression with negative sign, in contrast
for the specific case of isotropic materials, it reads
to the natural expectation µ˜ =Ωσ ǫ /2, which is the
el ik ik
E(α) ν(α) potentialenergydensity. Thereasonforthismodification
σi(kα) = 1+ν(α) ǫ(ikα)+ 1−2ν(α)δikǫ(llα) (3) is the coherency constraint which has to be fulfilled at
(cid:18) (cid:19) the interface. An illustrativeexampleto understandthis
for each phase. Here, E and ν are elastic modulus and unexpectedexpressionforthechemicalpotentialisgiven
Poissonratio,respectively. Wenotethateigenstraincon- in Appendix B. However, this effect is only important
tributionsduetodifferentunitcellsofthephasesarenot for solid-solid transformations. For crack propagation,
considered here for brevity, but they can easily be intro- where we assume that the new phase inside the crack is
duced [17]. For simplicity, we assume a two-dimensional infinitely soft, the normal and shear stresses vanish at
plane-strain situation. the interface, and therefore the discrepancy between the
Allfollowingrelationsareobtainedinaconsistentway chemicalpotential(6)andthenaiveguessµ˜ disappears.
el
fromvariationalprinciples,andthisisdescribedindetail We would also like to mention that the equality of the
inAppendixA. Theequationsofdynamicalelasticityare massdensityandthecoherencyleadstotheabsenceofki-
neticenergycontributionstothechemicalpotential. The
(α)
∂σ
ik =ρu¨ , (4) reason is that such a contribution is continuous across
i
∂xk the interface; finally, only the chemical potential differ-
for eachphase. Onthe interface,we obtainthe expected enceµ(1)−µ(2) entersintotheequationofmotionforthe
continuity of normal and shear stresses interface, and therefore such a term would cancel. Nev-
ertheless, we note that the kinetic energy contribution
σ(1) =σ(2). (5) would enter into the chemical potential with negative
in in
3
sign, i.e. −ρu˙2, because in the Lagrangian the kinetic and branching for high driving forces are rather generic
i
and potential energy contributions appear with opposite and are similarly valid for models with conserved order
sign. Kinetic contributions may play a role if instead of parameters. In [12], we derived the similar equations of
a phase transformationprocess,surface diffusionalong a motionifinsteadofaphasetransition,masstransportis
free boundary drives the evolution. due to surface diffusion along the free crack boundary.
For the motion of the interface, surface energy is also In the latter case, the elastic boundary conditions are
takenintoaccount. Sinceitdoesnotcoupletotheelastic replaced by
terms, it simply gives an additional contribution to the
chemical potential, σ(sd) =−ρu˙ v , (8)
in i n
µ =Ωγκ, and the chemical potential becomes
s
with the local interface curvature κ and the surface en- 1 1
µ(sd) =Ω σ ǫ − ρu˙2+γκ . (9)
ergy density γ. The curvature is positive if phase 1 is 2 ik ik 2 i
(cid:18) (cid:19)
convex. Then the motion of the interface due to a phase
transitionprocessisdescribedbythe localnormalveloc- It differs from the expression (6) first by the elastic en-
ity(Disakineticcoefficientwithdimension[D]=m2/s) ergydensity,becausenocoherencyconstraintshavetobe
fulfilled here. Second,the kinetic energydensityappears
D here, because it does not cancel in the derivation from
v =− (µ(1)−µ(2)+γκ), (7)
n γΩ el el theLagrangian. Theequationofmotionfortheinterface
is replaced by
which is positive if phase 1 grows. The set of equations
(1)-(7) describes the dynamics of the system. We point D(sd)∂2µ(sd)
v(sd) = (10)
outthatitleadstoacomplicatedfreeboundaryproblem, n γΩ ∂τ2
and the arising interfacial patterns are self-consistently
selected during this nonequilibrium process if external with the surface diffusion coefficient D(sd).
forces are applied to the system. For an initially almost In both cases, stresses on the boundary of the crack
flat interface between a soft and a hard phase, which tip with finite radius r scale as
0
is nonhydrostatically strained, first the ATG instability
develops: Longwavemorphologicalperturbationsleadto σ ∼Kr−1/2, (11)
0
a decrease of the total energy and the formation of deep
notches, similar to cracks (see e.g. [10, 18]). As we have and the curvature behaves as κ ∼ 1/r0. Therefore,
shownin[12,15,16],itisessentialtoincludetheinertial all contributions to the chemical potentials scale like
contributions,becauseotherwiseasteadystategrowthof µ∼r−1/2, and this is ultimately the reasonfor the cusp
0
these cracks is impossible, and the system collapses into singularityofthe Grinfeldinstability andthe impossibil-
the finite time cusp singularity of the ATG instability. ityofasteady-statecrackgrowth,ifonlystaticelasticity
is taken into account: Then, the equations of motion (7)
or(10)canberescaledtoanarbitrarytipradiuswhichis
III. CRACK PROPAGATION: SELECTION not selected by the dynamical process. The explanation
PRINCIPLES is that the linear theory of elasticity and surface energy
define only one lengthscale, the Griffith length, which
Inthecasethatonephaseisinfinitelysoft,crackprop- is macroscopic, but do not provide a microscopic scale
agationcanbestudied. Here,growthofthecrackisbased which allows the selection of a tip scale. Formally, the
on a phase transition of the solid matrix to a “dense gas equations of motion depend only on the dimensionless
phase” which has the same density as the solid. In this combinations vr /D for the phase transition dynamics
0
sense,it is similar to othermodels of fracturebasedona andvr3/D(sd) forsurfacediffusion;theradiusr andthe
0 0
non-conservedorder parameter [6, 7, 8]. The crucial dif- steady state velocity v therefore cannot be selected sep-
ferenceisthatthecurrentmodelisbasedonwell-defined arately; any rescaling which maintains the value of the
sharp interface equations, and therefore the predictions product would therefore describe another solution. The
do not depend on inherently numerical parameters like situationchangesifinertialeffectsaretakenintoaccount,
a phase field interface width. However, numerically, it whichis reasonablefor fast crackpropagation. Then ad-
requires a tedious separation of scales to obtain these ditionally the ratio v/v (v is the Rayleigh speed) ap-
R R
results; this is described in the next section. pears in the equations of motion, and therefore a rescal-
Understanding fracture as a phase transition process ing is no longer possible. Instead, D/v for the phase
R
offers many numerical advantages, as phase field models transition dynamics and (D(sd)/v )1/3 for surface diffu-
R
can be derived to solve this moving-boundary problem. sionsetthetipscale. Thus,weconcludethatfaststeady
We point out that the underlying selection principles state growth of cracks is possible if inertial effects are
which allow a steady state crack growth with propaga- taken into account. More formal analyses also including
tionvelocitieswellbelowtheRayleighspeed,tipblunting rigorous selection mechanisms due to the suppression of
4
growing crack openings far behind the tip are given in The elastodynamic equations are derived from the free
[15] for phase transition processes and in [12] for surface energybyvariationwithrespecttothedisplacementsu ,
i
diffusion.
δF
Thisanalysisshowsthattheselectionprincipleswhich
ρu¨ =− , (19)
i
allow a fast steady state growth of cracks are similar (cid:18)δui(cid:19)φ=const
for the simple phase transition process studied here and
surface diffusion. The latter does not require the intro- and the dissipative phase fields dynamics follows from
duction of a dense gas phase inside the crack and obeys
∂φ D δF
conservation of the solid mass itself. Even more, surface =− . (20)
∂t 3γξ δφ
diffusioncanbeunderstoodinageneralizedsenseasplas- (cid:18) (cid:19)ui=const
tic flow in a thin region around the extended tip which
Ithasbeenshownin[10]thatinthequasistaticcase,the
can be effectively described in the spirit of a lubrication
aboveequationsleadtothesharpinterfaceequations(4)-
approximation. Therefore, many general statements ob-
(7) if the interface width ξ is significantly smaller than
tainedforthephasetransitiondynamicscanalsobeused
all physical lengthscales present in the system. In Ap-
for crack growth propelled by surface diffusion. The lat-
pendix C, it is illustrated that this model also correctly
ter is more tedious to implement numerically, since the
incorporates the modification of the chemical potential
equation of motion (10) is of higher order [19].
(6) due to the coherency constraint. Details of the nu-
For both mechanisms, a tip splitting is possible for
merical implementation are given in Appendix D.
high applied tensions due to a secondary ATG instabil-
ity: Since σ ∼ Kr−1/2 in the tip region and the local
0
ATG length is L ∼ Eγ/σ2, an instability can occur,
G V. PHASE FIELD MODELING OF CRACK
provided that the tip radius becomes of the order of the PROPAGATION
ATG length. In dimensionless units, this leads to the
prediction ∆ ∼O(1).
split Thecentralpredictionofthistheoryoffractureisthat
a well-defined steady-state growth with finite tip radius
and velocities appreciably below the Rayleigh speed is
IV. PHASE FIELD MODEL
possible. This also cures the problem of the finite-time
cusp singularity of the Grinfeld instability. These pre-
To describe systems with moving boundaries accord-
dictions have been confirmed by phase field simulations
ing to the equations of motion developed above, we im-
[15]andsharpinterfacemethods [16]whicharebasedon
plemented a phase field model. Let φ denote the phase
a multipole expansion of the elastodynamic fields. Sur-
field with values φ = 1 for phase 1 and φ = 0 for phase
prisingly, it turned out that the obtained results seem
2. The energy density contributions are
to differ significantly: In particular, the sharp interface
method predicts a range of driving forces inside which
f =µ(φ)ǫ2 +λ(φ)(ǫ )2/2 (12)
el ij ii the velocity of the crack is a monotonically decreasing
function. Here, we demonstrate that the discrepancy of
for the elastic energy,with the interpolatedshear modu-
results is due to finite size effects of the previous phase
lus and Lam´e coefficient
field results [15], and that by careful extrapolation of
µ(φ) = h(φ)µ(1)+(1−h(φ))µ(2), (13) large-scalesimulations,acoincidingbehaviorisobtained.
We investigate crack growth in a strip geometry with
λ(φ) = h(φ)λ(1)+(1−h(φ))λ(2) (14)
fixed displacements at the upper and lower grip. The
multipole expansion technique [16] is designed to model
where
a perfect separation of the crack tip scale D/v to the
R
h(φ)=φ2(3−2φ) (15) strip width L: In most real cases, crack tips are very
tiny,andthereforeitistheoreticallydesirabletodescribe
interpolatesbetweenthephases,andthesuperscriptsde- this limit. For the phase field method, however, a finite
note the bulk values. The surface energy is strip width L is necessary, and a good separation of the
scales therefore requires time-consuming large-scale cal-
fs(φ)=3γξ(∇φ)2/2 (16) culations. We typically use strip lengths 2L and shift
the system such that the tip remains in the horizontal
with the interface width ξ. Finally,
center. This allows to study the propagation for long
times until the crack reaches a steady state situation.
f =6γφ2(1−φ)2/ξ (17)
dw Apart from this finite size restriction, we had to intro-
ducetheinterfacewidthξ asanumericalparameter,and
is the well-known double well potential. Thus the total
the phase field method delivers quantitative results only
free energy is given by
in the limit that all physical scales are much larger than
this lengthscale. The latter has to be noticeably larger
F = dV (f +f +f ). (18)
el s dw thanthe numericallatticeparameter∆x, butthe results
Z
5
2 0.9
D/v ξ = 2.33 D/ξv =3.09
R R
D/v ξ = 4.64 D/ξv =4.64
R R
D/vRξ = 9.27 0.8 D/ξvR=6.18
1 D/v ξ = 18.55 D/ξv =9.27
R R
D/ξv =18.55
R
D vR0.7
/R 0 /ξ
v L,
y v 0.6
-1
0.5
0.4
-2 0 0.05 0.1 0.15 0.2
-5 0
xv /D (ξ/L)1/2
R
FIG. 2: Crack shapes for different scale separations D/vRξ FIG. 3: First step of the extrapolation procedure for the di-
andfixedratioLvR/D=11.03;theaspectratioofthesystem mensionless velocity v/vR. The system size L/ξ is increased
is 2 : 1. The driving force is ∆ = 1.4. By improvement of and the ratio D/ξvR is kept fixed for each curve. For each
the separation, the crack opening is reduced, and finally the ratio, an extrapolated velocity vξ/vR corresponding toan in-
boundaries become straight parallel lines. finitesystemsizeisobtained,asindicatedbythedashedlines.
The drivingforce is ∆=1.4.
show that the choice ξ =5∆x is sufficient. We therefore
length scale dependencies). In the first step, we ex-
have to satisfy the hierarchy relation
tend the simulations to an infinite system size. There-
D fore, we decrease the ratio ξ/L → 0 for fixed tip scale
ξ ≪ ≪L, (21)
vR ratio D/ξvR. This step is demonstrated in Fig. 3 for
∆ = 1.4. Here, the dimensionless propagation velocity
which is numerically hard to achieve. We developed a v˜ =v /v isplottedasfunctionoftheinversesquare
L,ξ L,ξ R
parallel version of the phase field code which is run on root of the system size (ξ/L)1/2. In this representation,
up to 2048 processors, with system sizes up to 8192× the data for the larger systems can be extrapolated lin-
4096·(∆x)2. All computations are performed onthe su- early to infinite system sizes, since we numerically get a
percomputersJUMPandJUBLoperatedattheResearch scaling
Center Ju¨lich.
For the strip geometry, a dimensionless driving force D L D ξ 1/2
v˜ (∆, , )=v˜ (∆, )+α (23)
∆ is defined as L,ξ ξ
v ξ ξ v ξ L
R R (cid:18) (cid:19)
δ2(λ+2µ)
∆= , (22) for large systems, ξ/L ≪ 1, with a constant α > 0 for
4γL each curve. Since the separation of D/v to ξ is still
R
imperfect, the extrapolated values v˜ (∆,D/v ξ) do not
with δ being a fixed displacement by which the strip is ξ R
yet cumulate to a single point, and therefore a second
elongated vertically. The elastic constants of the new
extrapolation step is necessary.
phase inside the crack are zero. The value ∆ = 1 cor-
Hence, in Fig.4, the dependence of the velocity v /v
responds to the Griffith point. All calculations are done ξ R
ontheseparationparameterv ξ/Dfor∆=1.4isshown.
with Poisson ratio ν =1/3. R
The extrapolated values from Fig. 3 are used, and we
In [16], it was shown that close to the Griffith point,
obtain a scaling
dissipation free solution exists in the framework of the
model: Inthis regime1<∆<1.14anadditionalmicro- D v ξ
R
v˜ (∆, )=v˜(∆)−β (24)
scopiclengthscaleisneededtoselectthesmalltipradius ξ
v ξ D
R
which is no longer determined by the ratio D/v . This
R
can here be mimicked by the phase field interface width withaconstantβ >0andthe dimensionlesssharpinter-
andwasalreadydonein[15]. Here,wefocusonthemore face limit velocity v˜=v/v .
R
interestingregimeofhigherdrivingforces,butstillbelow Thistediousprocedurewasperformedforseveraldriv-
the threshold of instability. Typical crack shapes in the ing forces, and in Fig. 5 the comparison to the multi-
vicinity of the tip are shown in Fig. 2. pole expansion method [16] is shown. The agreement
To fulfill the scale separation (21), we perform a dou- of the results which are obtained from completely differ-
ble extrapolation of the obtained steady state velocities ent methods is very convincing. The small deviation for
v (the subscripts indicate the additional non-resolved ∆=1.8is due to the fact that this value is alreadyclose
L,ξ
6
0.7
R0.6
v
/
ξ
v
0.5
0 0.1 0.2 0.3 0.4
v ξ/D
R
FIG. 4: Second extrapolation step to obtain the sharp in-
terface velocity v. The extrapolated velocities obtained from
FIG. 6: Irregular tip splitting scenario for ∆=1.9. Weused
Fig. 3 are plotted as function of the scale separation param-
eter vR/ξD. In thisexample ∆=1.4 is used. sLyvsRte/mDis=2:414..2Taimnde Dis/gviRveξn=in9u.3n;itsthDe/avs2p.ecTt hreattiohicokfntehses
R
oftheinterfaceisthephasefieldinterfacewidth,indicatinga
good separation of thescales.
0.8
multipole expansion
0.7 phase field
0.6
withthe conjecture thatbranchingoccursas soonas the
0.5
steady state tip curvature becomes negative, leading to
R
v 0.4 the prediction ∆split ≈1.8 [16].
/
v The numerical determination of a characteristic crack
0.3
width scale in the sharp interface limit is more difficult,
0.2
and therefore we refrain from performing a double ex-
0.1 trapolation procedure. The explanation is that if the
soft phase inside the crack still possesses small nonva-
0
nishing elastic constants, the equilibrium situation far
1 1.5 2 2.5 3
behind the crack tip correspondsto a full opening of the
∆
crack, instead of the opening being of the order D/v :
R
As it is shown in Appendix B, the elastic energy is min-
FIG. 5: Comparison of the steady-state crack velocity ob-
tained from the multipole expansion technique [16] and the imized if the hard phase completely disappears. Small
extrapolated value from phase field simulations. remaining elastic constants can be due to an insufficient
separation of the scales D/v and ξ, since according to
R
Eqs. (13) and (14), the elastic constants decay only ex-
ponentially inside the crack, even if this soft phase has
tothethresholdofthetip-splittinginstabilitywhichcan- nominally vanishing elastic coefficients. Therefore, the
not be captured by the multipole expansion method. In crack opening is a weakly growing function of the dis-
particular, we find evidence for the prediction that the tance from the crack tip, and this slope becomes smaller
steady state velocitydecaysweaklywithincreasingdriv- withbetterscaleseparation,seeFig.2. Wepointoutthat
ing force, which might be an artefact of the model. thisopeningissolelyduetothephasetransitionprocess,
For higher driving forces, we observe tip-splitting in and the shapes are drawn without elastic displacements
the phase-field simulations (see Fig. 6). The onset of which should be added to obtain the real shape under
this irregularbehaviordependssensitivelyonthe system load. For example, the vertical displacement obeys the
size, because in relatively small systems, the branches usual scaling uy ∼ |x| for large distances |x| from the
of the crack cannot separate since they are repelled by tip.
p
the boundaries. Therefore,the steadystate growthis al- Thesameeffectcanbeseenifweinvestigatesolid-solid
ways stabilized by finite size effects. On the other hand, transformations towards a soft phase with small elastic
initial conditions can trigger an instability, and then a constants. The Poisson ratios in both the surround-
long transient is required to get back to steady state so- ing solid and the new inner phase are chosen equally,
lutions. Despite these restrictions, we are still able to ν = 1/3, but the bulk moduli differ by many orders of
makethe predictionthatthe thresholdofsplitting obeys magnitude. The softer the inner phase becomes, the less
∆ . 1.9 in the phase field model. It is in agreement theopeningofthe“crack”growswithincreasingdistance
split
7
2 N
λs/λh = 10-3 Fs = 3ξ γij (φi∇φj −φj∇φi)2dV, (28)
λ/λ = 10-3.5 4
s h i,j=1 Z
1 λs/λh = 10-4 X
λ/λ = 10-4.5
s h
N
D F = 3 γ φ2φ2dV. (29)
/R 0 dw ξ ij i j
v i,j=1 Z
y Xi6=j
Here, the interpolated elastic constants are
-1
N N
µ(Φ)= h(φ )µ(i), λ(Φ)= h(φ )λ(i), (30)
i i
-2
-5 0 Xi=1 Xi=1
xv /D
R with µ(i) and λ(i) being the elastic constants of the indi-
vidual bulk phases. Also, we have mutual interfacial en-
FIG. 7: Solid-solid transformation in a strip. A very soft
ergy coefficients γ = γ for phase boundaries between
phase grows at the expense of a harder phase. Parameters ij ji
iandj. Noticethatthird-phasecontributionsappearin-
are LvR/D = 11.03 (vR is the Rayleigh speed of the hard
side two-phase boundaries which lead to a renormaliza-
phase), D/ξvR =9.27, the aspect ratio is 2:1. The Poisson
tion of the bare interfacialenergies. This canbe avoided
ratioν =1/3isequalinbothphases,andthedrivingforceis
∆=1.4. by the addition of higher order terms to the multiwell-
potential (29). See Refs. [20, 21] for a discussion of this
issue.
Variation gives explicitly:
from the tip, see Fig. 7. Only very far away, the new
phase fills the whole channel. δF
el = h′(φ ) µ ǫ2 +λ (ǫ )2/2 ,
δφ k k ij k ii
k
N (cid:0) (cid:1)
VI. MULTIPHASE MODELING δF
s
= 3ξ γ 2(φ ∇φ −φ ∇φ )·∇φ
ik k i i k i
δφ
k
A simple approach to derive multiphase equations to Xi=1 h
describe systems consisting of more than two phases −φi(φi∇2φk−φk∇2φi) ,
starts again from variational principles. The volume N i
δF 12
fraction of each phase is described by a field variable dw = γ φ2φ .
φk, k = 1,...,N with N being the number of phases, δφk ξ i=1 ik i k
Φ=(φ ,...,φ ). Inthe sharpinterfacelimit, onephase Xi6=k
1 N
field variable has the value one inside the bulk phases,
Similarly, the elastic equation of motion is just the gen-
andtheothersarezero. Theirtemporalevolutionisgiven
eralization of Eq. (19)
by
∂φk =−2D˜ δF −Λ , (25) ρu¨i = ∂σ∂ikx(Φ) (31)
∂t 3ξ δφ k
(cid:18) k (cid:19)
with
whereweintroduceaLagrangemultipliertomaintainthe
phase conservation, N φ =1. We redefine the diffu-
k=1 k σ (Φ)=2µ(Φ)ǫ +λ(Φ)ǫ δ . (32)
sion coefficient D/γ → D˜, which is more appropriate to ik ik ll ik
P
generalize to arbitrary interfacial energy coefficients γik ThisdescriptionisadirectgeneralizationofRef.[20],and
between phases i and k. The additional factor 2 in the therefore it is known that it leads to appropriate sharp
equation of motion above is chosen to recover the previ- interface equations similar to Eqs. (1)-(7) and contact
ous phase field model in the case N =2. The expression angles as predicted by Young’s law [21].
for the Lagrange multiplier is given by
N
1 δF APPENDIX A: DERIVATION OF THE SHARP
Λ= . (26)
N δφ INTERFACE EQUATIONS
k
i=1
X
The free energy F = F +F +F has the following
el s dw Hereweexplainindetailhowthe equationsofmotion,
contributions:
the appropriate boundary conditions and the chemical
potential, which is responsible for interface motion, are
F = µ(Φ)ǫ2 +λ(Φ)(ǫ )2/2 dV, (27)
el ij ii derived in a unique way from variational principles. We
Z
(cid:0) (cid:1)
8
assume that the two phases are coherent, i.e. the dis- ∂σ(2)
− σ(1) δu dτ + ik δu dV
placement field is continuous across the interface, and in(1) i ∂x i
the mass densities are equal in both phases. Since we AZ(t) V2Z(t) k
do not consider lattice strains or surface tension here,
all elastic stresses arise from external forces. For sim- − σ(2) δu dτ .
in(2) i
plicity, we assume a two-dimensional plane-strain situa- #
Z
A(t)
tion. Contributions due to surface energy do not couple
to the elastic fields. Ultimately, they give only an addi-
The first integral is integrated by parts, assuming as
tivecurvature-dependenttermtothe chemicalpotential,
usual that the variations δu vanish for t and t . Since
i 0 1
as incorporatedin Eq. (7). We note that this expression
also the normal vectors of both phases are antiparallel,
was already obtained in [22], but here we take also in- n:=n(1) =−n(2), thus σ =−σ , we get
in(2) in
ertial effects into account, since the velocity of cracks is
typTihcaellkyinoeftticheenoerrdgeyrdoefntshietysoisunindbsoptehedpsh.ases δS = t1dt ∂σi(k1) −ρu¨ δu dV
i i
∂x
1 Z " Z k !
T = ρu˙2, (A1) t0 V1(t)
2 i
∂σ(2)
and the potential energy density reads + ik −ρu¨i δuidV
∂x
Z k !
U(α) = 1σ(α)ǫ(α). (A2) V2(t)
2 ik ik
− (σ(1)−σ(2))δu dτ .
Notice that certain components of the stress and strain in in i #
Z
tensors are in general discontinuous at the interface, as A(t)
will be elaborated below.
Demanding vanishing variation δS gives in the bulk the
We assume the total volume V of the entire system
usual equations of motion
to be constant in time and to be decomposed into two
subvolumes V(1)(t) and V(2)(t) of different solids. Up-
∂σ(α)
per indices discriminate between the phases (see Fig. 1). ik =ρu¨ , (A5)
i
The common interface A(t) := ∂V(1)(t)∩∂V(2)(t) with ∂xk
normal n and tangential τ is moving in time due to
and on the interface we obtain the continuity of normal
phase transitions, and consequently, the phase volumes
and shear stresses
aretime-dependentaswell. However,wedonotyetspec-
ify a concrete dynamical process here.
(1) (2)
σ =σ . (A6)
The Lagrangianis defined as in in
The next step is to calculate the change of the total
L(t)= TdV − U(1)dV − U(2)dV, (A3)
energy when the interface moves in the course of time.
Z Z Z
V V1(t) V2(t) Thisisdoneinthreesteps: First,wecalculatethechange
ofenergyduetothetimeevolutionoftheelasticfieldsfor
and the action is
fixed interface position. Second, we calculate the change
t1 of elastic energy due to the motion of the interface for
S = L(t)dt, (A4) fixed elastic fields in the bulk phases. After this second
Z step,the coherencyconditionatthe interfaceisviolated.
t0 In the last step, we therefore have to do additional work
with arbitrary beginning and end times t0 and t1. to adjust the displacements appropriately.
We obtain the usual elastic equations by varying the The first contribution is
action (A4) with respect to the displacement field for
fixed interface positions. Thus we get dW1 = ∂ (T +U(1))dV
dt ∂t
t1 V(1Z)(t)
δS = dt ρu˙ δu˙ dV − σ(1)δǫ(1)dV
" i i ik ik + ∂ (T +U(2))dV
Z Z Z
t0 V V1(t) Z ∂t
V(2)(t)
− σ(2)δǫ(2)dV
ik ik = ρu˙ u¨ dV + σ(1)ǫ˙(1)dV
Z # i i ik ik
V2(t) Z Z
V V(1)(t)
t1 ∂σ(1)
= dt ρu˙ δu˙ dV + ik δu dV + σ(2)ǫ˙(2)dV.
i i ∂x i ik ik
Z "Z Z k Z
t0 V V1(t) V(2)(t)
9
We note that the kinetic energy density is continuous ∂ u ,∂ u ,ǫ ,ǫ are discontinuous across the inter-
n n n τ nn nτ
across the interface. Furthermore, by the equations of face.
motion (A5) In the second step of energy calculation, we extended
smoothly the fields into the receding domain. The inter-
dW1 = u˙ ∂σi(k1)dV + u˙ ∂σi(k2)dV face at this new time t+∆t is now locatedat a different
dt i ∂x i ∂x position. This leads to discontinuities of the displace-
Z k Z k
V(1)(t) V(2)(t) ments,e.g.forthenormalcomponentatthenewposition
∂ ∂σ(1) of the interface
+ σ(1)u˙ dV − ik u˙ dV
∂x ik i ∂x i
V(1Z)(t) k (cid:16) (cid:17) V(1Z)(t) k ∆un = [∂nun](1)−[∂nun](2) vn∆t=(ǫ(n1n)−ǫ(n2n))vn∆t,
(2) (cid:16) (cid:17)
∂ ∂σ
+ σ(2)u˙ dV − ik u˙ dV where [...](α) denotes the evaluation of a probably dis-
∂x ik i ∂x i
V(2Z)(t) k (cid:16) (cid:17) V(2Z)(t) k continuous expression at the previous interface position,
taken for the phase α. Similarly, for the tangential com-
= σ(1) u˙ dτ + σ(2) u˙ dτ =0, ponent
in(1) i in(2) i
Z Z
A(t) A(t) ∆u = 2ǫ(1)−[∂ u ](1)+[κu ](1)−2ǫ(2)+[∂ u ](2)
τ nτ τ n τ nτ τ n
whereweassumedforsimplicitythatu˙i =0onallbound- −(cid:0) [κuτ](2) vn∆t
ariesapartfromA(t), i.e.no externalworkis exertedto
the solids. In the last step, we used the boundary condi- = 2(ǫ(n1τ)−ǫ(n(cid:1)2τ))vn∆t.
tions (A6), σ(1) =σ(2) =−σ(2) ; also, by definition, the
in in in(2) To zeroth order in ∆t, the stresses at the new inter-
displacement rate u˙ is continuous across the interface.
i face positionareequalonboth sides andidenticalto the
The above result is quite clear since the elastodynamic
stresses at the previous interface position. To reconnect
time evolution is purely conservative.
the displacements, we have to apply the coherency work
Thesecondcontributionarisesduetothemotionofthe
rate
interface for fixed elastic fields. We extend the elastic
state of the growing phase analytically into the newly dW
acquired region. This assures that the bulk equations dt3 = vn −(ǫ(n1n)−ǫ(n2n))σnn−2(ǫ(n1τ)−ǫ(n2τ))σnτ dτ.
remain fulfilled in both phases even after the forward AZ(t) h i
motion of the interface. Thus this contribution to the
energy change rate reads Altogether, the change of the energy is given by
dW
dt2 = vn U(1)−U(2) dτ. (A7) dW = d(W1+W2+W3)
AZ(t) (cid:16) (cid:17) dt dt
1 1
Theinterfacenormalvelocityispositiveifthephase1lo- = vn 2στ(1τ)ǫ(τ1τ)− 2σn(1n)ǫ(n1n)−σn(1τ)ǫ(n1τ)
cally extends. Here, we immediately used the continuity Z "(cid:18) (cid:19)
A(t)
of the kinetic energy density, which therefore cancels.
1 1
Afterthephasetransformationinthissecondstep,the − σ(2)ǫ(2)− σ(2)ǫ(2)−σ(2)ǫ(2) dτ.
displacements are no longer continuous at the interface. (cid:18)2 ττ ττ 2 nn nn nτ nτ(cid:19)#
Thusextraworkhastobeinvestedtoremovethismisfit.
In the local coordinate system n and τ (see Fig. 1) the We can therefore define an appropriate chemical poten-
strain tensor becomes tial for each phase at the coherent interface
ǫnn = ∂nun, (A8) µ(α) =Ω 1σ(α)ǫ(α)− 1σ(α)ǫ(α)−σ(α)ǫ(α) . (A11)
ǫττ = ∂τuτ +κun, (A9) el 2 ττ ττ 2 nn nn nτ nτ
(cid:18) (cid:19)
1
ǫ =ǫ = (∂ u +∂ u −κu ). (A10)
nτ τn τ n n τ τ Noticethat, incontrastto afreesurface,the normaland
2
shearcontributionsappearwithnegativesign. Then,the
Here,κisthe interfacecurvature,whichispositiveifthe
energy dissipation rate can be written as
phase 1 is convex.
At this point, a few comments concerning the conti- dW 1
nuity of various fields across the coherent interface are dt = Ω (µ(e1l)−µ(e2l))vndτ. (A12)
in order. Since the displacement field has to be con- Z
A(t)
tinuous across the interface, also its tangential deriva-
tives are continuous, but the normal derivatives are Due to the coherency condition and the requirement of
not. Consequently, the following quantities are con- equal mass density, the kinetic energy density does not
tinuous: ∂ u ,∂ u ,κu ,κu ,ǫ . On the other hand, appear.
τ τ τ n n τ ττ
10
δ δ strip (per unit length) as function of L(1):
1 1
y U(L(1)) = L(1)σ(1)ǫ(1)+ L(2)σ(2)ǫ(2)
(1) 2 ik ik 2 ik ik
L
L 1 δ2(2µ(1)+λ(1))
x
= . (B3)
2L(1)+(L−L(1))2µ(1)+λ(1)
2µ(2)+λ(2)
(2)
L
ThisfunctionismonotonicinL(1). Assumingthatphase
1 is harder than phase 2, 2µ(1)+λ(1) >2µ(2)+λ(2), the
energy is minimized if L(1) = 0, i.e. if the hard phase
FIG. 8: Motion of a planar interface in a strip geometry. disappears. Notice that on the other hand, the elastic
Vertically, a constant displacement δ is applied to the grips energy density is higher in the softer phase, σ(1)ǫ(1)/2<
of the strip which consists of two solid phases with different ik ik
σ(2)ǫ(2)/2. We get from Eq. (B3)
elastic constants. ik ik
dU(1) 1 1 1
=− σ(1)ǫ(1)+ σ(2)ǫ(2) = µ(1)−µ(2) .
dL(1) 2 yy yy 2 yy yy Ω el el
APPENDIX B: MOTION OF A PLANAR PHASE
(cid:16) (cid:17)
BOUNDARY
Here we clearly see that the negative elastic energy den-
sityforthenormaldirectionentersintotheenergychange
The elastic contribution to the chemical potential (6)
rateandthusintothechemicalpotential;inmoregeneral
attheinterfacebetweentwosolidphaseshastheremark-
cases one can easily verify that this is also true for shear
ablepropertythatthenormalandshearcontributionsare
contributions.
negative definite. This means that growth of the phase
with higher elastic energy density at the expense of the
phasewithlowerenergydensitycanstillreducethetotal
energy. In order to illuminate this point, we consider a APPENDIX C: PHASE FIELD MODELING OF
simple example. COHERENT SOLIDS
Here, two strips of different solid materials are coher-
ently connected (see Fig. 8), and the interface can move
The aim of this section is to show that the phase field
due to a phase transition. We assume for simplicity
model presented in Section IV leads to the correctsharp
that the process is slow and inertial effects can be ne-
interfacelimitevenforsolid-solidtransformations,where
glected. We apply a fixed displacement δ at the upper
the chemical potential differs from the usual expression
end of this layer structure and set the displacement at
[23]. Since the treatment of the surface energy contribu-
the lower grip to zero. The total strip width L is dis-
tioniswell-knownandentersadditivelyintothechemical
tributed among the two layers, L=L(1)+L(2) =const.
potential, we focus on the elastic fields here.
We have a homogeneous strain situation in each of the
phases (α = 1,2), i.e. u(α) = 0, u(α) = u(α)(y) with For smooth interfaces, all stresses remain finite even
x y y inthe sharpinterfacelimit. We note that this statement
(α)
a pure linear dependence of u on y. Strains are
y holds even for fracture processes, where usually stresses
(α) (α) (α)
ǫxx = ǫxy = 0, ǫyy = ∂yuy = const. The elongation candivergeatsharpcornersandtipsintheframeworkof
of each phase is δ(α) = L(α)ǫ(α). In sum, they are equal thelineartheoryofelasticity;nevertheless,inthecurrent
yy
to the prescribed total opening δ = δ(1) +δ(2) = const. description, the tip radius r0 is always finite and there-
Stresses and strain are connected through Hooke’s law fore stresses are limited to values σ ∼ Kr−1/2, where
0
for isotropic materials, thus σ(α) = λ(α)ǫ(α), σ(α) = 0, K is a stress intensity factor. Consequently, the dis-
xx yy xy
σ(α) = (2µ(α)+λ(α))ǫ(α). At the interface, the equality placement field must be continuous in the sharp inter-
yy yy
face limit, because a finite mismatch δu would lead to
of normal stresses, σ(1) =σ(2), leads to
yy yy diverging stresses and strains ǫ∼σ/E ∼δu/ξ, and thus
a divergentelastic energy. Then, obviously,the equation
2µ(1)+λ(1)
ǫ(2) = ǫ(1). (B1) of motion (19) leads to the usual elastic bulk equations
yy 2µ(2)+λ(2) yy
and also to the continuity of normal and shear stresses
at the interface in the limit ξ →0.
The strain can be computed in terms of the given total
opening, Next, we calculate the total elastic energy change due
totheinterfacemotion. Wedeterminetheenergycontri-
δ
ǫ(1) = . (B2) butions and changes in a line perpendicular through the
yy L(1)+(L−L(1))2µ(1)+λ(1) interface, assuming the interface curvature to be small,
2µ(2)+λ(2)
κξ ≪1. Todoso,weintroducealocalcoordinatennor-
Hence, we can calculate the total elastic energy in the maltotheinterface(pointingfromphase1intophase2).