Thermal fracture as a framework for quasi-static crack propagation
aa r X i v : . [ c ond - m a t . m t r l - s c i ] S e p Thermal fracture as a framework for quasi-static crackpropagation
F. Corson † , M. Adda-Bedia † , H. Henry ‡ and E. Katzav †† Laboratoire de Physique Statistique, ENS, Paris VI, Paris VII, CNRS,24 rue Lhomond, 75005 Paris, France. ‡ Laboratoire de Physique de la Mati`ere condens´ee, CNRS UMR 7643,Ecole Polytechnique, 91128 Palaiseau Cedex, France.November 8, 2018
Abstract
We address analytically and numerically the problem of crack path prediction in themodel system of a crack propagating under thermal loading. We show that one can explainthe instability from a straight to a wavy crack propagation by using only the principle oflocal symmetry and the Griffith criterion. We then argue that the calculations of the stressintensity factors can be combined with the standard crack propagation criteria to obtainthe evolution equation for the crack tip within any loading configuration. The theoreticalresults of the thermal crack problem agree with the numerical simulations we performedusing a phase field model. Moreover, it turns out that the phase-field model allows to clarifythe nature of the transition between straight and oscillatory cracks which is shown to besupercritical.
Crack path prediction is one of the main challenges in the field of fracture mechanics. The reasonis simple - a satisfactory equation of motion of a crack tip is associated with a fundamental under-standing of material separation mechanisms (Freund, 1990; Broberg, 1999; Fineberg and Marder,1999; Adda-Bedia et al., 1999; Leblond, 2003). From a more general point of view and alongsideits importance in many applications, the determination of crack propagation laws is necessaryto describe the fracture phenomenon as a pattern formation process induced by mechanicalstresses. Within the framework of Linear Elastic Fracture Mechanics (LEFM), the propagationof a crack is mainly governed by the singular behavior of the stress field in the vicinity of itstip (Freund, 1990; Broberg, 1999; Leblond, 2003). For a two-dimensional quasi-static crack,which is the main purpose of the present work, this behavior is given by σ ij ( r, φ ) = K I √ πr Σ I ij ( φ ) + K II √ πr Σ II ij ( φ ) + O ( r ) , (1)where Σ I ij ( φ ) and Σ II ij ( φ ) are universal functions describing the angular variation of the stressfield, and K I and K II are the stress intensity factors (SIFs). The evolution of the crack tip isgoverned by the Griffith energy criterion (Griffith, 1920; Freund, 1990; Broberg, 1999; Leblond,2003), which states that the intensity of the loading necessary to induce propagation is given by G = Γ, where G is the energy release rate and Γ is the fracture energy of the material, that is1he energy needed to create new free surfaces. This criterion can be rewritten using the stressintensity factor, in which case it is referred to as the Irwin criterion (Irwin, 1957) K I = K Ic ≡ p µ Γ , (2)where K Ic is the toughness of the material and µ is the shear modulus. The Griffith criterion wasoriginally formulated for quasi-static crack propagation (Griffith, 1920), and later generalizedto rapidly moving cracks (Freund, 1990). While this criterion is very useful in predicting crackinitiation, it cannot predict the direction of the crack tip, and therefore in most cases it is notsufficient to determine the actual path of the crack. In order to achieve this, several suggestionshave been made. Among them, the Principle of Local Symmetry (PLS) states that the crackadvances in such a way that in-plane shear stress vanishes in the vicinity of the crack tip, orexplicitly K II = 0 , (3)This rule was proposed for in-plane quasi-static cracks (Gol’dstein and Salganik, 1974; Leblond,1989), and then generalized to rapidly moving cracks (Adda-Bedia et al., 1999). Note thatfrom a historical point of view, the first criterion for crack path selection, based on symmetryarguments similar to the PLS, was first formulated for antiplane (mode III) crack propagation(Barenblatt and Cherepanov, 1961). Recently, the dynamics of a rough crack in two dimensionshas been successfully described by an equation of motion derived using the PLS (Katzav et al.,2007a). Another suggestion, based on symmetry arguments, and which recovers the principle oflocal symmetry in a certain limit, was proposed by Hodgdon and Sethna (1993), who formulatedan equation of motion of the crack tip. However, this formulation introduces additional lengthscales, which are not known a priori. A different approach, known as the maximum energyrelease rate criterion (Erdogan and Sih, 1963) states that the crack advances in a direction thatmaximizes its energy release rate. Interestingly, even though this approach is quite differentfrom the PLS, it produces very similar results to those obtained using the PLS to such an extentthat they were even conjectured to coincide (Bilby and Cardew, 1975). However, as shownin (Amestoy and Leblond, 1992) for kinked cracks and in (Katzav et al., 2007b) for branchedcracks, the results are not exactly the same and a clear distinction can be made between thetwo. Actually, it has been shown that the PLS is the only self-consistent criterion (Leblond,1989, 2003). A theoretical challenge is to explain the current rich crack path phenomenology,including the various known instabilities, using a pure LEFM approach combined with theGriffith criterion and the PLS. This has been recently argued for roughening instabilities ofcracks in (Katzav et al., 2007a), as well as for the branching instability in (Katzav et al., 2007b)and here we wish to contribute further to this effort.Trying to address the problem of crack path prediction from the experimental side, the ther-mal crack problem has attracted a lot of attention. The propagation of cracks induced by ther-mal gradients has been widely studied (Ronsin et al., 1995; Ronsin, 1996; Yuse and Sano, 1997;Ronsin and Perrin, 1998; Yang and Ravi-Chandar, 2001; Deegan et al., 2003; Yoneyama et al.,2006; Sakaue et al., 2008; Yoneyama et al., 2008) since the work of Yuse and Sano (1993). In atypical experiment, a glass strip with a notch at its end is pulled at a constant velocity froman oven into a cold bath (see Fig. 1). The control parameters in this experiment are mainlythe width of the strip, the temperature gradient between the oven and the cold bath, and thepulling velocity. If the velocity is small enough, a crack does not propagate. Above a firstcritical velocity, the crack starts propagating following a straight centered path, and above asecond critical velocity, the crack begins to oscillate with a well defined wavelength. However,the nature of the transition from straight to wavy path is not fully understood. For example,the experimental results (Ronsin, 1996; Yuse and Sano, 1997) are not sufficient to determine2hether the bifurcation is super-critical or sub-critical (also referred to as continuous or dis-continuous transition) . At higher pulling velocities, irregular oscillations and branching areobserved. Interestingly, these are very reproducible regimes, which makes the thermal crack anideal model experiment as it allows to study the slow propagation of cracks under well controlledconditions.In an attempt to provide a theoretical explanation, Marder (1994) successfully determinedthe propagation threshold by calculating K I for a straight crack propagating at the center of theplate. He also proposed a description of the oscillatory instability using the stability criterionderived by Cotterell and Rice (1980) for a crack in an infinite plate (the celebrated T -criterion),but that implied the improbable requirement that the fracture energy should depend stronglyon the velocity. Eventually, this result was found to be incompatible with the experimentalevidence (Ronsin et al., 1995), and serves as a classical example for the violation of the T -criterion. Sasa et al. (1994) attempted to predict the threshold of the oscillatory instability andits wavelength by using the PLS under the infinite-plate approximation. Although they gave areasonable qualitative description, their results deviated systematically from the experimentalmeasurements. The next step was taken by Adda-Bedia and Pomeau (1995), who calculated K II directly for a sinusoidal perturbation in a finite strip, in the configuration where the crack tip isat the center of the strip. However, in order to predict the oscillatory instability they resortedto an additional criterion, namely that the stability with respect to this perturbation dependson the sign of K II . More recently, Bouchbinder et al. (2003) wrote an amplitude equation fora sinusoidal perturbation using the propagation criterion of Hodgdon and Sethna (1993). Oneshould note that all these approaches rely on additional criteria when determining the instabilitythreshold and cannot do so using only the classical principles of crack propagation. In addition,theoretical results, except an attempt in (Bouchbinder et al., 2003), are mostly limited to thestudy of the linear stability of a straight crack or to the study of transient regimes. Thissituation calls for a reexamination of the thermal crack problem in order to explain it usingsolid arguments.In this work, we begin with the study of quasistatic cracks in a quenched glass plate usinga phase-field approach. Phase field models were originally introduced by Caginalp and Fife(1986); Collins and Levine (1986) to describe the propagation of solidification fronts. Ex-tensions of this work allowed quantitative modeling of dendritic growth (Karma and Rappel,1998) and since then, phase field methods have been applied to many solidification prob-lems such as growth of binary alloys (Etchebarria et al., 2004) and formation of polycristallinesolids (Kobayashia and Warren, 2005). Furthermore, they have been extended to other freeboundary problems including viscous fingering (Folch et al., 2000) and crack propagation (Karma et al.,2001; Karma and Lobkovsky, 2004; Henry and Levine, 2004; Hakim and Karma, 2005, 2009).In the latter case, the phase field approach turns out to be very similar to the variationalapproach of fracture which uses the description of free discontinuity problems (Braides, 1998;Acerbi and Braides, 1999) using the Γ-convergence (Ambrosio and Tortorelli, 1990, 1992). Thenumerical results of the phase-field model we derive allow to clarify the nature of the transitionbetween straight and oscillatory cracks.We then provide a theoretical analysis of the experiment of Yuse and Sano (1993). As inprevious works (Adda-Bedia and Pomeau, 1995; Bouchbinder et al., 2003), the analysis is basedon the calculation of the SIFs. However, we calculate them to first order with respect to a straightcentered crack for an arbitrary trajectory in a strip of finite width. We are then able to determinethe instability threshold and the wavelength of the oscillations using only the PLS without In (Yuse and Sano, 1997), the bifurcation is claimed to be supercritical. Nonetheless, data points close tothreshold are not accurate enough to rule out any other hypothesis (see Fig. 5 on page 372 in (Yuse and Sano,1997)). Results presented in (Ronsin, 1996) suffer from the same lack of data points close to threshold. xblv bT + ∆ TT OvenCold bath
Figure 1: The experimental setup. A glass plate of width 2 b is moved slowly at a constantvelocity v from an oven maintained at temperature T + ∆ T to a cold bath of fixed temperature T . The crack tip is steady in the laboratory frame (moving in the plate frame) and its positionwith respect to the cold bath ( bl ) is controlled by the temperature field.introducing any additional criterion. Results agree very well with the numerical simulationsusing the phase field model. And comparing them with the results of Adda-Bedia and Pomeau(1995), we find a small correction for the instability threshold, and a significant one for thewavelength. Also, unlike Bouchbinder et al. (2003), we argue that the calculation of the SIFs canbe combined with the standard crack propagation criteria (that is with the Irwin criterion andthe PLS) to obtain an evolution equation for the crack tip. Our main conclusion is that the PLSprovides a good description of crack paths without any additional criteria. Recently, a similarapproach for the thermal crack problem has been proposed and solved numerically (Pham et al.,2008) which is in agreement with both our theoretical findings and with experimental resultsof Ronsin and Perrin (1998).The paper is organized as follows. We first provide a description of the thermal crackproblem. We then describe the phase field approach for this problem, and compare its resultswith those of a theoretical analysis we perform. An important result of the theoretical analysisis an evolution equation for the crack path which allows quantitative predictions. We concludeby discussing the relevance of this work to path prediction of cracks in general. We first introduce the thermal crack problem and set the basic definitions and notations. Here,we consider the idealized problem of an infinitely long strip containing a semi-infinite crack. Thecoordinate system is chosen so that the axis of symmetry of the strip corresponds to y = 0 andthe tip of the crack lies at x = 0 (see Fig. 1). We assume that the temperature field is uniformacross the thickness of the strip, and that the strip is under plane stress conditions, so that theproblem is actually two-dimensional. We will also use dimensionless variables, measuring thelengths in units of the half-width of the strip b , the temperature in units of ∆ T , the strains inunits of α T ∆ T /b , and the stresses in units of Eα T ∆ T , where E is the Young’s modulus and α T the coefficient of thermal expansion. As will be shown, the behavior of the system is mainlygoverned by two dimensionless parameters : the Peclet number P = bv/D , where D is thethermal diffusion coefficient and the dimensionless material toughness ˆ K Ic = K Ic / ( Eα T ∆ T √ b ).Provided the velocity of the plate is not too low, the temperature field, which is assumed notto be affected by the presence of the crack, can be described by the advection-diffusion equation ∂ xx T + P ∂ x T = 0, with the boundary conditions: T ( x ) = 0 in the cold bath and T ( x ) = 1 for4 → ∞ (Marder, 1994). The solution of this equation is T ( x ) = (cid:16) − e − P ( x + l ) (cid:17) Θ( x + l ) , (4)where x = − l corresponds to the location of the surface of the cold bath (see Fig. 1) and Θ( x ) isthe Heaviside function. Recall that the temperature field as given by Eq. (4) is an approximationof the experimental one. For a quantitative comparison with experiments, the real temperaturefield should be used (Pham et al., 2008).We consider the problem of a crack propagating through the heated strip within the frame-work of linear elastic fracture mechanics. Under plane stress conditions, the two-dimensionalstress tensor σ is related to the two-dimensional strain tensor ǫ by σ ij = 11 − ν [(1 − ν ) ǫ ij + νǫ kk δ ij − (1 + ν ) T ( x ) δ ij ] , (5)where ν is the Poisson ratio and the strain tensor ǫ is related to the displacement field ~u by ǫ ij = 12 (cid:20) ∂u i ∂x j + ∂u j ∂x i (cid:21) . (6)The strip is free from external traction, thus within any approach one should always satisfy theboundary conditions σ yy ( x, ±
1) = σ xy ( x, ±
1) = 0 . (7) The classical theory of crack propagation, where the crack behavior is determined by the singu-larities of the stress field at the crack tip leads to difficult numerical issues when considering themovement and interaction of many cracks or to track crack motion in three dimensional systemswhere even the dynamics of crack fronts is not well understood. In this context, the phasefield approach to crack propagation (Aranson et al., 2000; Karma et al., 2001; Eastgate et al.,2002; Marconi and Jagla, 2005) is very useful as it allows to study problems in arbitrary geome-tries and go beyond simple crack paths. It has succeeded in reproducing qualitative behaviorof cracks such as branching (Karma and Lobkovsky, 2004) and oscillations under biaxial strain(Henry and Levine, 2004). More recently, it has been shown that the Griffith criterion and thePLS are embedded in the phase field model of crack propagation (Hakim and Karma, 2009).The general idea of phase field modeling is to introduce an additional field (the phase field)that describes the state of the system. The main advantage of this approach is that one does notneed to treat the interface explicitly since it is defined implicitly as an isosurface of the phasefield. This advantage becomes clear when studying the evolution of complex shapes since thealgorithmic cost of using an additional field is much smaller than the cost of dealing explicitlywith complex moving surfaces. In the case of fracture, this phase field φ indicates whether thematerial is intact ( φ = 1) or broken ( φ = 0). In the usual sharp interface representation, acrack surface is an infinitely thin boundary between a region where the material is intact and anempty region that can not sustain any stress. In contrast, the equations of the phase field modelare such that φ varies continuously in space, and a crack surface is represented by a region offinite thickness. The thickness can be chosen freely and does not affect numerical results as longas it is much smaller than the characteristic length scale of the studied phenomenon. It can beshown that through an appropriate coupling between the phase field φ and the elastic fields,the model can have the desired properties : no stresses are transmitted across a crack (in the5imit where the system size is much larger than the width of the diffuse interface) and one canassociate a finite surface energy with the interface, so that the fracture energy is well-defined.Here, we use this model to describe a quasistatic crack propagating due to a thermal gra-dient. We begin by presenting the model, including the adjustments needed to incorporatethermoelastic effects. The dimensionless elastic energy is written as (Adda-Bedia and Pomeau,1995) E el = Z Z dxdy g ( φ ) " − ν ) ( ǫ kk − T ( x )) + 12(1 + ν ) (cid:18) ǫ ij − δ ij ǫ kk (cid:19) , (8)where the function g ( φ ) = φ (4 − φ ) describes the coupling between the phase field and the elas-tic fields. This function is such that in regions where the material is broken ( φ = 0), the contri-bution to the elastic energy is zero, while in regions where the material is intact, the contributionto the elastic energy recovers the one prescribed by linear elasticity. In (Karma et al., 2001), it isshown that the particular choice of g does not affect the results as long as g (0) = g ′ (0) = g ′ (1) = 0and lim φ =0 g ( φ ) ∼ φ α , with α >
2. The additional term T ( x ) corresponds to the thermal ex-pansion. This effect is assumed to be the same in the intact and in the partially broken partsof the material.Since we are considering a quasi-static crack growth, in which the elastic body is at me-chanical equilibrium at all times, the elastic energy should be an extremum with respect to thedisplacement field ~u ( x, y ). That is δE el δu i = 0 . (9)We now need to specify the evolution of the phase field. This is done by associating with thephase field a dimensionless free energy given by (Henry and Levine, 2004) E φ = Z Z dxdy (cid:20) D φ ∇ φ ) + V ( φ ) + g ( φ ) ( E φ − E c ) (cid:21) , (10)where E φ is defined by E φ = − ν ) ( ǫ kk − T ( x )) + ν ) (cid:16) ǫ ij − δ ij ǫ kk (cid:17) if ( ǫ ii − T ( x )) > ν ) (cid:16) ǫ ij − δ ij ǫ kk (cid:17) if ( ǫ ii − T ( x )) < , (11)where V ( φ ) = 1 . h (1 − φ ) φ is a double well potential, with an energy barrier of height ∝ h , thatensures that the preferred states of the homogeneous system are either φ = 1 (intact) or φ = 0(completely broken). E φ coincides with the standard elastic energy density when the materialis under tension and includes only the contribution of shear when the material is under com-pression. This choice avoids the propagation of a crack under compression (Henry and Levine,2004). Taking into account the coupling through the function g ( φ ), when E φ is larger than athreshold value E c , the broken phase is favored while when E φ < E c , the intact phase is favored.The evolution equation for φ is relaxational with a kinetic bias that ensures that the materialwill not be driven from a broken state to an intact state : τ ˙ φ = min (cid:18) − δE φ δφ , (cid:19) , (12)where τ is a constant. The irreversibility of Eq. (12) means that a crack cannot heal (whichwould happen in our set up at the crack tail) and ensures that the total elastic energy stored inthe material cannot increase when φ varies. 6he elastic equation (9) is solved on a rectangular domain [ − L , L ] × [ − ,
1] with L ≫ x = ± L/
2. Herewe have chosen ǫ ij ( − L/ , y ) = 0 , (13) σ xy ( L/ , y ) = σ xx ( L/ , y ) = 0 . (14)It is clear that the final solution is insensitive to these boundary conditions as long as L ≫ /P .The phase field equation (12) is solved using an Euler forward scheme (Press et al., 2007). Thevalues of the grid spacing dx and of the diffusion coefficient D φ are chosen so that the diffuseinterface occupies few grid points and is much smaller than the characteristic wavelength ofthe oscillating crack, while the characteristic time τ is chosen so that the phase-field relaxesfast enough to its equilibrium (or more precisely, multiplying τ by a factor 2 does not affectquantitatively the results.). Numerical checks were performed to ensure that neither changes inthe convergence criterion of the Gauss-Seidel iterative method nor changes in the grid resolutionaffect the results. The size of the domain [ L ] × [2] was typically 1500 ×
300 grid points, andnumerical checks showed that doubling L did not affect the results. Finally, the fracture energyΓ = K / (2 µ ) was computed using the expression (Karma et al., 2001)Γ = Z dφ q D φ ( E c − ( E c g ( φ ) − V ( φ )) . (15)leading to the dimensionless material toughness ˆ K Ic = √ µ Γ / ( Eα T ∆ T √ b ).Before turning to the description of numerical results, we find it necessary to summarize thedifferent parameters used here and to give some rationale for their choice. In addition to thephysical parameters associated to the material and to the temperature field, the pertinent phasefield parameters are τ , D φ , h and E c . In dimensionless units the values of τ = 1, D φ = 4 . × − , h = 1 and E c = 1 have been chosen. Notice that the three latter parameters determine boththe fracture energy and the interface width of the phase field. Since we are considering thedimensionless material toughness ˆ K Ic , the only effect we are concerned with, is the dependanceof the width of the diffuse interface (the broken surface) on D φ . More precisely, in the spiritof the phase field approach, the value of D φ , for a given value of ˆ K Ic should not affect thebehaviour of the crack. Within the present simulations, we have checked that using a value of D φ corresponding to an interface which is twice thinner does not affect quantitatively the finalresults.We now summarize the results of the simulations. By tuning the control parameters P andˆ K Ic , we could observe the different regimes observed in experiments : no crack propagation,propagation of a straight crack at the center of the strip or a crack following a wavy path (seeFig. 2). However, the main purpose of the present simulations was to determine the nature of thebifurcation from straight to oscillatory cracks, since no direct experimental results are available.Namely in experiments, the amplitude of oscillations varies almost linearly with the controlparameter, but since experimental results do not allow to explore the behavior of the crackclose to the threshold it was impossible to characterize the crack path in the close vicinity of thetransition point. Furthermore, it should be emphasized that apart from attempts describing thistransition by a Landau-Ginzburg expansion (Bahr et al., 1995; Bouchbinder et al., 2003), theavailable theoretical work has been mostly limited to linear stability analysis. Using the phasefield model and without resorting to any assumption, we were able to compute the steady-statepaths and to determine the nature of the bifurcation.First, in Fig. 3 we present some typical crack paths close to the threshold for the sake ofcomparison with the analytical results presented below and also in order to clarify the nature7igure 2: Typical crack paths obtained using the phase field model. The white region corre-sponds to the broken surface and the color code corresponds to the density of elastic energy inthe material (red: high elastic energy density, blue: zero elastic energy). (a) A straight centeredpropagating crack. (b)
A wavy crack above the instability threshold. This picture shows thatthe crack width is much smaller than the wavelength of the observed oscillations. -0.15 0 0.15 0 2 4 -0.25 0 0.25 0 1.35 2.7
Figure 3: (a) : Crack tip trajectory below the threshold of linear instability for a crack whichwas initially off center. P = 7 . / ˆ K Ic = 7 . (b) : Crack tip trajectory for an initiallyslightly off center crack above the threshold of linear instability. P = 7 . / ˆ K Ic = 8 . slightly off center (shifted by 2 grid points away from the center) and slightly above thethreshold. One can see oscillations of the path that are amplified as the crack advances. Theseoscillations eventually saturate and the crack path becomes oscillating with a finite constantamplitude as observed in experiments. Using the model we were able to compute the amplitudeof the oscillations when a steady state was reached. The simulations showed that the steadystate is independant of initial conditions, such as the initial position of the crack, and is a singlevalued function of the control parameters. Typical results are shown in Fig. 4 where the squareof the amplitude is plotted as a function of 1 / ˆ K c for a given value of P . One can see thatabove a threshold value of 1 / ˆ K c the amplitude of oscillations is no longer zero and scales as thesquare root of the deviation from the threshold. The same behavior holds when using the valueof P as a control parameter. The amplitude curves together with the damped oscillations belowthe threshold and the amplified oscillations above the threshold indicate that the bifurcation issupercritical or continuous.For values of the control parameters well above the threshold, we were not able to re-cover complex paths similar to those observed in experiments (Yuse and Sano, 1993). This maybe due to the fact that well above the threshold the quasistatic approximation is no longervalid and dynamic effects through the temperature field as given by Eq. 4 become relevant.It is nonetheless interesting to note that the observed crescent-like path in Fig. 5 exhibits a8 A (a) 0 0.1 0.2 11 13 15 A (b) Figure 4: (a)
Square of the amplitude A of the oscillations for P = 7 .
89 as function of 1 / ˆ K Ic ,which is related to the amplitude of the temperature gradient. (b) : same as (a) with P = 3 . A in the neighborhood of the transition shows that the transition fromstraight to oscillatory propagation is supercritical (continuous). -0.33 0 0.33 0 1.35 2.7 Figure 5: Crack tip trajectory for a crack deep in the nonlinear regime, that is far above thethreshold of the instability. The values of the control parameters are P = 7 . / ˆ K Ic =11 .
6. Note the striking resemblance to crack paths appearing in (Ghatak and Mahadevan, 2003;Audoly et al., 2005; Sendova and Willis, 2003).9 scillating modestraight modeno propagation(a) 1 /P / ˆ K Ic (b) λ /P Figure 6: Main theoretical and numerical results. (a) The phase diagram of the different statesof the system as a function of the two control parameters P and ˆ K Ic . (b) The wavelength λ as a function of P at the transition from straight to wavy crack propagation. In both figures,the symbols + and × correspond to the numerical results obtained from the phase field model,the solid lines correspond to the analytical calculation presented in Sec. 4 and the dashed linescorrespond to the results reported in (Adda-Bedia and Pomeau, 1995).striking similitude with paths obtained when cutting a thin elastic sheet with a moving thickobject (Ghatak and Mahadevan, 2003; Audoly et al., 2005) and with crack paths obtained dur-ing drying of silicate sol-gel films fabricated using spin-coating techniques (Sendova and Willis,2003).The phase diagram in Fig. 6(a) shows the threshold for the propagation of a centered straightcrack and for the transition to a wavy crack propagation in the 1 /P –1 / ˆ K c phase space. Itcan be seen that the theoretical predictions presented in Sec. 4 are in quantitative agreementwith numerical results of the phase field model without any adjustable parameters. The sameagreement is reached when considering the dimensionless wavelength at the threshold of thetransition from straight to wavy crack path (see Fig.6(b)). For this latter quantity, one shouldalso note the good agreement with experimental values for small values of 1 /P (linear behaviorwith λ ≃ .
28 + 2 . /P ) (Ronsin et al., 1995). Moreover, the phase field computation and ourtheoretical approach lead to the same deviation from the linear behavior for 1 /P > .
2, wherethe discrepancy between the present results and those reported in (Adda-Bedia and Pomeau,1995) become significant. This feature has also been observed in experiments (Ronsin, 1996).
We now turn to a theoretical study of the transition from straight to oscillatory crack propa-gation. As in previous theoretical treatments (Marder, 1994; Adda-Bedia and Pomeau, 1995;Bouchbinder et al., 2003), the analysis is based on the calculation of the stress intensity factors(SIFs). However, here, they are computed to first order for an arbitrary trajectory withoutusing the infinite plate approximation and without the restriction that the crack tip lies at thecenter of the sample. Within the framework of linear elastic fracture mechanics, the propagation10f a crack is governed by the singular behavior of the stress field in the vicinity of its tip, asdefined in Eq. (1). The propagation laws we use are the Irwin criterion defined in Eq. (2) andthe principle of local symmetry (PLS) defined in Eq. (3). To calculate the SIFs, one must solvethe equilibrium equations ∂σ ij ∂x j = 0 , ∇ σ ii = −∇ T ( x ) , (16)with the boundary conditions (7) at the sides of the plate and the additional boundary conditionat the crack surface σ ij n j = 0 , (17)where −→ n is the normal to the crack faces.To investigate the transition from a straight to an oscillating regime, we consider a smalldeviation from a straight centered crack, y ( x ) = Ah ( x ) , (18)where h ( x ) is defined for x ≤ | A | ≪
1. We expand the displacement andstress fields to first order as u i = u (0) i + Au (1) i + O (cid:0) A (cid:1) , (19) σ ij = σ (0) ij + Aσ (1) ij + O (cid:0) A (cid:1) . (20)The SIF K I (resp. K II ) is an even (resp. odd) function of A . Therefore their expansions havethe form K I = K (0)I + O ( A ) , (21) K II = AK (1)II + O ( A ) . (22)Note that in particular, K II = 0 for a straight centered crack. The different terms of the aboveexpansion can be determined by solving the equilibrium equations using repeatedly the Wiener-Hopf method. The details of these calculations are presented in Appendix A. Here we simplysummarize the calculation of K (1)II , which is novel. K (1)II can be written as K (1)II [ h ( x )] = κh (0) + K (1)II [ h ( x ) − h (0)] , (23)which amounts to decomposing the general problem into that of a straight off-center crack andthat of an oscillating crack with a centered tip. The latter problem was solved by Adda-Bedia and Pomeau(1995), yielding K (1)II [ h ( x ) − h (0)] = K (0)I h ′ (0) + Z −∞ dx ∂∂x h σ (0) xx ( x, h ( x ) − h (0)) i p + ( − x ) , (24)where p + ( x ) is a weight function that does not depend on the physical parameters. We have gen-eralized the results of Adda-Bedia and Pomeau (1995) with the calculation of the contribution κh (0) to obtain K (1)II [ h ( x )] = (cid:20) cK (0)I − ∂∂l K (0)I (cid:21) h (0) + K (0)I h ′ (0) + Z −∞ dx ∂∂x h σ (0) xx ( x, h ( x ) i p + ( − x ) , (25)where c ≃ .
272 is a numerical constant (see Appendix A for details).Having calculated the SIFs to first order, we can now examine the particular form taken bythe crack propagation criteria for a small perturbation. The Irwin criterion reads K (0)I = ˆ K Ic .11 K (0)I Figure 7: K (0)I as a function of l , for P = 5. Recall that l is the distance of the crack tip fromthe boundary of the cold bath. The threshold of incipient growth of a centered straight crack isgiven by ˆ K Ic = K (0)I ( l c ), where l c is the location of the maximum of K (0)I ( l ).The transition to a propagating straight crack is thus governed by the existence of l such that K (0)I ( l ) = ˆ K Ic . This occurs when ˆ K Ic = K (0)I ( l c ), where l c is the location of the maximumof K (0)I ( l ) (see Fig. 7). This condition determines the transition curve in Fig. 2 from no crackpropagation to a propagation of a centered straight crack. Also, above the propagation thresholdEq. (21) shows that to first order in A , the distance between the cold front and the crack tip,which is determined by K I [ h ( x )], does not depend on the crack path y ( x ) = Ah ( x ). Thissimplifies the problem, as we can parameterize the “time” evolution of the crack path by theposition of its tip in a coordinate system attached to the plate. If h ( x ) now describes the crackpath in such a coordinate system, the value of K (1)II when the tip is at the point x is simplygiven by K (1)II ( x ) = (cid:20) cK (0)I − ∂∂l K (0)I (cid:21) h ( x ) + K (0)I h ′ ( x ) + Z −∞ du ∂∂u h σ (0) xx ( u, h ( x + u ) i p + ( − u ) . (26)Therefore, the PLS which for a small perturbation reads K (1)II ( x ) = 0, yields0 = (cid:20) cK (0)I − ∂∂l K (0)I (cid:21) h ( x ) + K (0)I h ′ ( x ) + Z −∞ du ∂∂u h σ (0) xx ( u, h ( x + u ) i p + ( − u ) . (27)Eq. (27) is a linear equation for a crack that deviates slightly from the straight centered path.To determine the stability of the straight centered crack propagation, we must examine theevolution of a small perturbation, which amounts to solving the following problem : given a crackpath h ( x ) defined for x ≤
0, which does not necessarily satisfy K (1)II ( x ) = 0 for x <
0, find acontinuation of h such that K (1)II ( x ) = 0 for x >
0. Now, since σ (0) xx ( x,
0) decreases exponentiallyas x → −∞ , one can expect the long-term behavior of the crack to be the same as that of acrack path satisfying K (1)II ( x ) = 0 for all x , and this is indeed observed in the simulations (seeFig. 3(b) in the previous section, and the discussion below). Therefore, examining crack pathsthat satisfy K (1) II = 0 from the beginning is sufficient to determine the stability threshold, andwe will use such paths for that purpose. 12e look for solutions of the form h ( x ) = Re [ δ exp( iqx )] , (28)where Re stands for the real part, δ and q are in general complex numbers. Then using Eq. (26), K (1)II ( x ) takes the form K (1)II ( x ) = Re [ δ K II ( q ) exp( iqx )] , (29)where K II ( q ) = (cid:20) cK (0)I − ∂∂l K (0)I (cid:21) + iq K (0)I Z −∞ du ∂∂u h σ (0) xx ( u, e iqu i p + ( − u ) . (30)A solution K (1)II ( x ) = 0 is obtained if K II ( q ) = 0, which is the dispersion relation of the problem.Since q and K II ( q ) are in general complex numbers, we can expect the equation K II ( q ) = 0to admit a discrete set of solutions for a given set of the dimensionless parameters P and ˆ K Ic .By examining the behavior of K II ( q ), we find that the stability is determined by solutionscorresponding to a pair of conjugate oscillating modes, which are shown on Fig. 8(a) for a givenvalue of P . When Im( q ) >
0, the perturbation (28) is damped which implies that the straightand centered crack propagation is stable. When Im( q ) <
0, the perturbation is amplified, leadingto a wavy crack propagation regime. Thus the instability threshold corresponds to Im( q ) = 0,and is denoted by q ≡ q c . Using this result, another way to determine the instability threshold isillustrated in Fig. 8(b). It consists in supposing from the beginning that q ≡ ω is a real number and in looking for the specific value of ω = q c for which Im[ K II ( ω )] = 0 and Re[ K II ( ω )] = 0simultaneously.By repeating the above computation for different values of P and ˆ K Ic , we obtained in Fig. 6(a)the phase diagram and in Fig. 6(b) the wavelength of the oscillations at the threshold, whichis given by λ = 2 π/q c . These figures show that a good quantitative agreement is reached withthe numerical results of the phase field model for both the wavy instability threshold and thewavelength at the transition.Our results for the instability threshold differ from those of Adda-Bedia and Pomeau (1995),in which the threshold was determined by the existence of a real solution q c such that K (1)II [sin ωx ] = dK (1)II dω [sin ωx ] = 0 for x = 0 and the stability criterion was based on the sign of K (1)II [sin ωx ]. Ingeneral, such a solution has no reason to satisfy K (1)II ( x ) = 0 for x = 0. Indeed, the wavelengthspredicted by these two criteria (i.e., the present one and the one of Adda-Bedia and Pomeau(1995)) differ significantly, especially for small values of P (see Fig. 6(b)). The thresholds pre-dicted also differ, but less notably (see Fig. 6(a)). As can be seen in Fig. 8(b), the deviationsbetween the estimates of ˆ K Ic and λ obtained by the two methods are related to the vertical andhorizontal distances between the minimum of Im( K II ) and its intersection with Re( K II ), andthe vertical distance scales as the square of the horizontal one. In this paper we address the problem of crack path prediction for quasi-static crack propagation.We used the framework of the thermal crack experiment since it is a well controlled experiment,and the best known to us. The reason is that stresses are induced internally by an imposed we use ω instead of q in order to emphasize the fact that it is a real number m( q )Re( q )(a) 1 / ˆ K Ic Im( K II )Re( K II )(b) ω/ π Figure 8: (a) Solutions of K II ( q ) = 0 versus 1 / ˆ K Ic for P = 5. The solution branch disappearswhen Im( q ) exceeds a certain value, as the integral appearing in the calculation of K II ( q ) ceasesto converge. (b) The real and imaginary part of K II ( ω ), with ω real, at the instability thresholdfor P = 5 ( ˆ K Ic ≃ . The straight off-center crack
In the following, we present briefly the method of resolution for a straight off-center crack in athermal gradient. The method of resolution follows the same steps as in (Adda-Bedia and Pomeau,1995). Consider a straight crack which is slightly above/below the center of the strip, namelyits location is given by y ( x ) = Ah (0) ≡ δ . Due to the absence of symmetry with respect to thecenter of the strip y = 0, let us split the displacement, strain and stress fields in the strip asfollows f ( x, y ) = 12 (cid:0) f + ( x, y )Θ( y − δ ) + f − ( x, y )Θ( δ − y ) (cid:1) , (31)where Θ is the heaviside function and f stands for any component of the elastic fields. Let usalso define [ f ]( x ) = 12 (cid:0) f + ( x, δ ) − f − ( x, δ ) (cid:1) , (32)Note that all the elastic fields are continuous in the unbroken regions. However, along the cracksurfaces the displacements might be discontinuous while σ yy and σ xy are continuous there. Thus,for this specific problem, the boundary conditions (7,17) can be rewritten as σ ± yy ( x, ±
1) = σ ± xy ( x, ±
1) = 0 , (33)[ σ yy ] ( x ) = [ σ xy ] ( x ) = 0 , (34)[ u x ] ( x ) = [ u y ] ( x ) = 0 for x > , (35) σ ± yy ( x, δ ) = σ ± xy ( x, δ ) = 0 for x < , (36)Now, we Fourier transform along the x direction all the elastic fields and the temperaturefield and plug them into the equilibrium equations (16) and into the boundary conditions (33,34).After some algebraic manipulations, this yields relations of the type (cid:18) [ˆ u x ] ( k )[ˆ u y ] ( k ) (cid:19) = (cid:18) C xx ( k ) C xy ( k ) D x ( k ) C yx ( k ) C yy ( k ) D y ( k ) (cid:19) ˆ σ xy ( k, δ )ˆ σ yy ( k, δ )ˆ T ( k ) , (37)Performing the first-order expansions in δ , as given by Eqs. (19)-(20):ˆ u i ( k, δ ) = ˆ u (0) i ( k,
0) + δ ˆ u (1) i ( k, , (38)ˆ σ ij ( k, δ ) = ˆ σ (0) ij ( k,
0) + δ ˆ σ (1) ij ( k, , (39)and inverting (37) gives the final resultˆ σ (0) yy ( k,
0) = − F ( k ) h ˆ u (0) y i ( k,
0) + D ( k ) ˆ T ( k ) , (40)ˆ σ (0) xx ( k,
0) = H ( k )ˆ σ (0) yy ( k,
0) + S ( k ) ˆ T ( k ) , (41)ˆ σ (1) xy ( k,
0) = − P ( k ) (cid:16)h ˆ u (1) x i ( k, − ik h ˆ u (0) y i ( k, (cid:17) + ik ˆ σ (0) xx ( k, , (42)where F ( k ) = k sinh k − k sinh 2 k + 2 k , D ( k ) = 2 (1 − cosh k )(sinh k − k )sinh 2 k + 2 k ,H ( k ) = k + sinh k sinh k − k , S ( k ) = sinh k − k sinh k + k ,P ( k ) = k sinh k − k sinh 2 k − k . (43)15he first step of resolution which involves only the zeroth-order expansion, i.e. the case of astraight centered crack, is the calculation of K (0)I . Using the boundary conditions (35,36) to lead-ing order and applying the Wiener-Hopf method to equation (40) yields (Adda-Bedia and Pomeau,1995) K (0)I = Z + ∞−∞ D ( k ) ˆ T ( k ) F + ( k ) dk π , (44)where F ( k ) = F − ( k ) F + ( k ) , (45)and F + ( k ) ( F − ( k )) is analytic in the upper (lower) half-plane, respectively.To complete the calculation of the SIFs, we return to the case of a straight off-center crackand determine the coefficient κ appearing in Eq. (23). Using the boundary conditions (35,36)to first order in δ and applying again the Wiener-Hopf method to Eq. (42) yields h ˆ u (1) x i ( k ) = ik h ˆ u (0) y i ( k ) − P − ( k ) Z + ∞−∞ ik ′ P + ( k ′ )ˆ σ (0) − xx ( k ′ , k ′ − k + iǫ dk ′ iπ + αP − ( k )= − kF − ( k ) Z + ∞−∞ D ( k ′ ) ˆ T ( k ′ ) F + ( k ′ ) k ′ − k + iǫ dk ′ π − P − ( k ) Z + ∞−∞ ik ′ P + ( k ′ )ˆ σ (0) − xx ( k ′ , k ′ − k + iǫ dk ′ iπ + αP − ( k ) , (46)where P ( k ) = P − ( k ) /P + ( k ) with P + ( k ) ( P − ( k )) analytic in the upper (lower) half-plane,and ǫ is a vanishingly small positive number. notice that the solution of h ˆ u (0) y i ( k ) obtainedfrom the zeroth order calculations has been used, and that σ (0) xx ( x,
0) can be obtained fromEq. (41). For both F ( k ) and P ( k ) the factorization into F ± ( k ) and P ± ( k ) can be achievedusing a Pad´e decomposition, as previously explained in (Bouchbinder et al., 2003; Katzav et al.,2007b). The unknown constant α is introduced because the decomposition involved in theWiener-Hopf method is not unique. The value of α can be determined by noting that the h ˆ u (1) x i ( k ) behaves as O ( k − / ) as k → ∞ . Canceling out the k − / term on the r.h.s. yields α = − Z + ∞−∞ D ( k ) ˆ T ( k ) F + ( k ) dk π , (47)and thus one has h ˆ u (1) x i ( k ) = (cid:18) F − ( k ) − P − ( k ) (cid:19) Z + ∞−∞ D ( k ′ ) ˆ T ( k ′ ) F + ( k ′ ) dk ′ π + 1 F − ( k ) Z + ∞−∞ − ik ′ D ( k ′ ) ˆ T ( k ′ ) F + ( k ′ ) k ′ − k + iǫ dk ′ iπ . (48)Finally by looking at asymptotic behavior of h ˆ u (1) x i ( k ) as k → ∞ one has κ = cK (0)I − ∂∂l K (0)I − Z −∞ ∂∂x σ (0) xx ( x, p + ( − x ) dx , (49)where p + ( x ) is the inverse Fourier transform of P + ( k ) and c is a numerical constant given by c = lim k →∞ ( ik ) √ (cid:20) F − ( k ) − P − ( k ) (cid:21) ≃ . . (50)Using this value of κ in Eq. (23) yields Eq. (25).16 eferenceseferences