Table Of Content

Phase Field Modeling of Fracture and Stress Induced Phase Transitions R. Spatschek, C. Mu¨ller-Gugenberger, and E. Brener Institut fu¨r Festkörperforschung, 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 obtained 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 interface (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 instability: 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é 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 interface 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 parameter 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 obtained 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ï = ∂σ∂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ï δ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 derivatives 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).