Numerical experiments in 2D variational fracture

In the present work we present some results of numerical experiments obtained with a variational model for quasi-static Griffith-type brittle fracture. Essentially the analysis is based on a recent formulation by Francfort and Marigo the main difference being the fact that we rely on local rather than on global minimization. Propagation of fracture is obtained by minimizing, in a step by step process, a form of energy that is the sum of bulk and interface terms. To solve the problem numerically we adopt discontinuous finite elements based on variable meshes and search for the minima of the energy through descent methods. We use a sort of mesh dependent relaxation of the interface energy to get out of small energy wells. The relaxation consists in the adoption of a carefully tailored cohesive type interface energy, tending to the Griffith limit as the mesh size tends to zero.


INTRODUCTION
ur research group has developed in the last decade a computer model for fracture nucleation and propagation in 2d brittle solids, based on variational fracture (see [1], [2] and references therein). The present work reports some numerical results obtained with this computer model in quasi-static problems of fracture nucleation and propagation. Essentially the analysis is based on the formulation of Francfort and Marigo [3] the main difference being the fact that we rely on local rather than on global minimization. Nucleation and propagation of fracture is obtained by minimizing in a step by step process a form of energy that is the sum of bulk and interface terms. Recent attempts of producing numerical codes for variational fracture (see for example the works of Bourdin [4] or Del Piero et al [5]) are based on the approximation of the energy, in the sense of Gamma-convergence, by means of elliptic functionals (e.g. Ambrosio-Tortorelli approximations [6]). Here instead we adopt discontinuous finite elements and search for local minima of the energy, close to equilibrium trajectories, through descent methods. In constructing such a numerical model there are two main questions that must be answered. 1. The variational model for fracture requires the ability to accurately approximate the location of cracks, as well as their length. Cracks may not be restricted to propagate along the skeleton of a fixed finite element mesh. To overcome this O difficulty our mesh is made variable in the sense that node positions in the reference configuration of the body are considered as further variables. 2. On adopting a Griffith type interface energy the initiation of fracture in a previously virgin body is always "brutal", that is the system cannot proceed along neutral equilibrium paths or, in other words, descent directions toward local minima do not exist. To reach local minima the system must have the ability to surmount small energy barriers and then proceed through "non-equilibrium" descent paths. Here we use a sort of mesh dependent relaxation of the interface energy to get out of small energy wells. The relaxation consists in the adoption of a carefully tailored cohesive type interface energy, tending to the Griffith limit as the mesh size tends to zero. In the first part of the paper we give the main ingredients of the theory and of the model implementation. In the final part of the paper, some examples concerning fracture nucleation and propagation are presented.

GRIFFITH'S VERSUS VARIATIONAL FRACTURE MECHANICS
n the theory of brittle fracture due to Griffith, cracking is associated to the existence of surface energy and the propagation of cracks to the competition between surface energy and potential energy release during an infinitesimal increase of the crack length "c". A number of authors, since the pioneering work of Ambrosio and Braides [7] in 1995, has proposed variational models of brittle fracture based on global energy minimization with free discontinuities (for early thoughts on minimization with free discontinuities see De Giorgi [8]). The basic idea is to use the competition between bulk and interface energy that is implicit in Griffith's theory, [9], as the driving criterion for optimal crack selection. Global energy minimization remedies to two main shortcomings of Griffith's theory, namely  Initiation: Griffith's is a model for crack growth from a pre-existing crack, not a model for crack nucleation. Indeed the cost of surface energy for the creation of an infinitesimally long crack cannot be compensated by the elastic energy release in an originally "flawless" body.  Crack path selection: Griffith's model is restricted to the study of crack propagation along a given path. In the original formulation the energy balance is expressed in terms of "c" only. The variational formulation based on global minimization which looks at the crack set as one of the sought unknowns, leads naturally to the search of crack paths without any restriction on crack topology, thus allowing, at least in principle, the study of initiation and, possibly, branching. There are two main criticisms to variational fracture: 1. The minimum energy formulation implies reversibility but fracture is an inherently irreversible phenomenon. 2. Global minimality does not seem physically correct. The "infinitesimal" local minimality of the Griffith's theory is presumably not correct either but at least does not lead to paradoxes such as the impossibility of equilibrium under purely tensile loading. In 1998 Francfort & Marigo formulated a consistent variational model of fracture respecting the irreversibility of the phenomenon, but again based on global minimality. The problem is formulated as a quasi-static evolution problem in which the geometry and size of the cracks are limited by their precedessors. A more realistic approach that would investigate local minimizers still lacks of the necessary mathematical apparatus (though the recents attempts of Dal Maso & Toader [10] and Larsen [11] must be aknowledged). Our work, as described in the Introduction, is devoted to assess, through numerical experiments, the ability of a variational model of fracture based on local minimization to model fracture initiation and propagation. From the numerical point of view the main issue with local minimization is how the material moves from an equilibrium state to another. Though it seems natural to look at "gradient flows", so that the material moves downhill in the steepest direction, in studying initiation (that is the onset of cracks in a previously flawless body) the crack cannot open up in any "nearby" optimal direction without having the energy to increase somewhere. A way out to restore physical plausibility is to assume that the system can overcome "small" energy barriers, where the term "small" refers both to the amount of the energy increase before the descent begins and to the length of the path leading to the point of descent. Such a "smallness" could be measured by introducing a convenient barrier-norm to be compared with a given threshold specific of the material considered. Since an increase in energy always requires the work of external active forces, these energy jumps should be actually forbidden in a quasi-static evolution and can be justified only admitting the presence of finite disturbances in the geometry and in the data.

I ANALYTICAL AND NUMERICAL FORMULATION OF THE PROBLEM
ur approach is reported in details in [2] and it will be briefly summarized for the reader's convenience.

Preliminaries
Consider a two-dimensional body B occupying in the original configuration a bounded plane domain  and undergoing small deformations. Let  be the infinitesimal strain tensor and u the displacement field defined over . The boundary of is partitioned into two parts where displacements u and tractions p are given. We admit that u may be discontinuous on a set  assumed to belong to an admissible set of cracks S(). The crack pattern is the most relevant unknown in variational fracture problems and the theorems of Ambrosio [12], valid in particular for Griffith's type variational formulations based on global minimizers, ensure that  is sufficiently regular to be approximated by sets composed by regular arcs, such as line elements. Therefore  is a 1d closed set and the unit normal n and tangent t exist a.e. along it. The jump of u across , denoted [[u]] = u + −u -, is resolved into two components relative to n, t If a is positive the corresponding points on the two opposite sides of the interface  separate and determine a vacancy. Finally we assume that u belongs to H 1 (\), that is u has the regularity required in linear elasticity away from .

Volume and surface energy densities
The material is assumed to be elastic and isotropic, that is characterized by the elastic energy density a scalar field defined over \. In (2)  and  are the elastic moduli for generalized plane stress.
To model brittle fracture, we introduce on  the interface energy density where G c is the thoughness of the material and a,b are the components of the displacement jump along .

Typical Boundary Value Problems
We consider quasi-static BVP's. A typical example is represented by the plane strip of height H and width B shown in Fig. 1 to which we refer for notations. On the lateral sides of the rectangle the traction condition p = 0 is enforced and on the bases a combination of relative displacements u is given. We assume u to be time dependent. At any time t of the loading process we seek equilibrium states of the strip as points of local minima of the functional under the boundary conditions described above.

Space discretization
To discretize the problem we split the domain into triangles, as shown in Fig. 1, and identify  with the part of the skeleton of the triangulation that is not at the boundary of . In order to duplicate nodes and edges of the skeleton of the mesh, we introduce special interface elements with zero thickness placed along the edges of the continuous elements. A topological view of these elements is depicted in Fig. 13a for a plane strip. In order to give to the skeleton the ability to choose optimal crack paths, the positions of the nodes of the triangulation in the reference configuration , are taken as further variables and can be moved to minimize the total energy.

Approximate interface (surface) energy density
In the numerical applications we approximate the equilibrium trajectory of the system by considering crack propagation as based on critical points of the energy. To get out of possible small energy wells, either numerical (due to the finite element mesh) or physical (due to fracture initiation) on introducing the limit tensile stress °, the limit shear stress ° the thoughness =G c , the shear stiffness k and the functions ( where  = 1.90086, we adopt the following approximate relaxed form of the surface energy density: In (6)° represents the slope at the inflection point of the curve depicted in Fig. 2a, ° the slope of the dotted lines in Fig.  2b.
Notice that while  and k are constants, °, °, are actually mesh dependent, that is their values must be tuned depending on the current mesh size h and should tend to infinity as h tends to zero, linearly with1/h. Although the density (13) has a shape close to that of the approximate surface energy density defined in [2], the main difference between them is that all parameters in the function (6) have a clear physical meaning and no shape parameter is needed.

FRACTURE NUCLEATION IN 1D
he study of a simple 1d example can clarify the main ingredients of our approach. Consider the long strip of brittle material depicted in Fig. 3, modeled as a 1d continuous bar, fixed at the left end and subject to a given displacement u at the other end. The bar has been discretized into two finite elements of length L/2. Inside the elements the material is linearly elastic and characterized by the elastic energy density =1/2 EA  2 . The spirit of our approximation is to consider that fractures, that is discontinuities in the displacements, may occur at the interface between elements (the skeleton of the mesh), that is at point 2. To allow for the interface energy at node 2 we must duplicate node 2 into two spatially coincident nodes in the reference configuration : 2', 2". Therefore we are reduced to a two degrees of freedom structure. The two degrees of freedom we choose are {u,a} and are defined in terms of the displacements u(2'), u(2") as follows: u= u(2'), a = u(2")-u(2'). Then the total energy is The graph of the energy (u,a) for a particular value of the given displacement u is depicted in Fig. 4. In the figure are visible two minima: one is ½ EA u 2 L, the "elastic" minimum for a=0 and u=u/2 (point A), that is the one for which there is no discontinuity; the other one is G c , the minimum for a=u and u=0 (point B), that is the bar is completely fractured. For u small the "elastic" minimum is smaller than the "completly fractured" minimum. As the displacement u is gradually increased from zero, with the system sitting in the "elastic" minimum, the energy minimum rises and becomes larger than G c . As it is evident from Fig. 4, wathever large be the elastic minimum with respect to G c , there is always an energy barrier to sormount to go from the local to the global minimum, that is there is no descent direction from A to B. The numerical method we propose to get out of the "elastic pocket" is to adopt for the interface energy the relaxed form depicted in Fig. 2a. With such cohesive type energy the system will access a descent direction at u= L °/E. The stress thresold is set again to infinity as soon as the system accesses the gradient flow from A to B; the fractured minimum is T reached through steepest descent. In this example the mesh skeleton is fixed, therefore the crack can open up only at point 2. A way to make the crack move is to consider the mesh 1, 2, 3 variable in the sense that the coordinate x(2) of point 2 is taken as a further variable in the energy: and the minimum is searched with respect to u, a, x (2). Actually all the fractured minima are equivalent if the bar is homogeneous and there is no reason for the bar to move the crack from the midpoint. This can be easily verified by preminimizing (u,a,x(2)) with respect to u (that gives u°=u x(2)/L) and substituting u° for u in the energy: that is actually idependent of x(2). To make the node movements worthwile we consider that the bar is weackened at x=0.366L (Figure: 5), that is the thoughness Gc is not homogeneous and has a minimum at x=0.366L. The bar is discretized into 10 elements in such a way that there is no node at the weackest point. The energy is then a function of 27 variables: 9 node displacements, 9 displacement jumps, 9 node positions: {u(i), a(i), x(i)}, i=2,...,10. A value of the displacement u, large enough to have Gc<1/2EA/L u 2 , is assigned and the minimum state is searched with a descent method using as starting state the displacement field depicted in Fig. 6a.  Through descent the system evolves toward the minimum as showni in Fig. 6. The evolution of the domain and of the node positions at five steps of the evolution is depicted in Fig. 7. In particular in Fig. 7 it can be seen how node 5 moves from the original position to the "weackest point" by creating an interface at x=0.366L.

VALIDATION WITH CLASSICAL LINEAR FRACTURE MECHANICS (CLFM): PROPAGATION OF A STRAIGHT CRACK IN MODE I
oving from 1d to 2d we enter a much more complex world. The ability of the method to mimic the onset of cracks in brittle materials in 2d becomes a difficult problem of approximation of the real quasi-static trajectory of the system with a step-wise sequence of states, and depends on a sensible choice of the parameters °, °. Such parameters are mesh dependent, in the sense that must be tuned according to the mesh size h, in order to make the approximate solution mesh independent. We must say that the story of our "battle" to obtain sound numerical results from our experiments on fracture simulation is full of "casualties". In scientific papers most often appears that the proposed numerical models are good enough to catch the essence of the problems at hand, but there are scant of no reports on the endless numbers of failures. Actually no one is interested in a bunch of meaningless, junk results and we do not want to change this tradition. Here we want just to emphasize that, except from very large and trivial bounds, we were not able to identify any simple rule for the choice of the parameters °, °, that up to now has proceeded in a case by case fashion. In our effort to establish the ability of simulating crack nucleation and propagation in 2d brittle solids by adopting a numerical model based on variational fracture, the first attempt has been to verify the capacity to reproduce the results of CLFM in cases in which the transition to fracture occurs through stationary states of the system. Such is the case of the propagation of a straight crack in mode I, into a strip of homogeneous, isotropic, linearly elastic material. The geometry of the problem and the main notation is reported in Fig. 8. M In CLFM it is assumed that the crack propagates along the line x2=0 when the given displacement U() reaches a critical value. This assumption, largely confirmed by a vast number of experimental results, makes in this case the question not "the how" but "the when". In Fig. 9 the outcome of an experimental test on an isotropic polycarbonate is reported into two photoelastic images, the first taken before and the second after the onset of the crack. In Fig. 10b the quasi-static path of meta-stable states (local minima of the energy) followed by the system in the {c,U} plane, is represented with a solid line. The shaded line represents the part of the path that is made of stationary points of the energy that are not local minima. The values taken by the total energy and by the elastic energy along the stationary path are depicted in Fig. 10a. In Fig. 10c the total energy at increasing values of the given displacement U is reported as a function of c. The dots indicate the stationary points of . The qualitative results of our numerical simulation are reported in Fig. 11, where 12 steps of the discretized evolution are shown. Notice that these are stable equilibrium states of the system, that is local minima, reached through our descent minimization algorithm at twelve values of the imposed displacement (U=0.00,0.005,0.010,0.015,…,0.065). Denoting M the total number of triangles in the mesh and N the total number of nodes, the energy to be minimized is a function of (2) 3 M nodal displacements u1(n), u2(n) (n=1,2,..,3M) and of 2 N nodal positions in the reference configuration x 1 (m), x 2 (m) (m=1,2,..,N). The energy is a complex, non convex function of its arguments with a lot of local minima. From Fig. 11 we see that the fracture propagates at a critical value of the displacement whose value depends on the ratio °/h. In this example we focus on the tuning of the parameter °. In Fig. 12 the values taken by the elastic energy at some steps of the discrete evolution are compared with the graph of Fig. 10a. The two polylines a and b refer to two values of the parameter °: 246Ncm -2 and 100Ncm -2 . For the first value the structures follows closely the exact energy path, even if with a sequence of snaps about the real trajectory.

NUCLEATION IN 2D. RUPTURE OF A STRETCHED AND SHEARED STRIP
he strip of base B and length L=2B, represented in Fig. 13 is considered. In (a) a topological representation of the original mesh is given: it is a structured mesh based on a uniform triangulation of the domain of mesh size h=L/20 2 . In grey the 100 physical triangles of the mesh and in white the 270 interface elements are reported.
Therefore the free node displacements, before considering the geometric constraints at the bases of the strip, are 2(3)100= 600 instead of the 2 (66)=132 of a classical FEM mesh. The node positions in the reference configuration are also taken as variables. Therefore there are 2(66)=132 more variables to be considered. The energy is a complex functions of these variables with a lot of local minima. T Such energy is minimized by fixing a given combination of stretch and shear at the constrained bases, and evolving an initial state compatible with this datum (the elastic solution, Fig. 13b) through a descent algorithm (conjugate gradient). Fig. 13a,b shows the final geometry of the strip and of the corresponding mesh in the original configuration, at the end of the minimization process. The result we show resemble the outcome of a number of fracture tests on concrete but we refrain the reader from being too excited by this. We must report indeed the rather frustrating results that one can obtain by just changing a little the starting state or the parameters (especially important is the interplay of ° and °). Besides we must point out that the result we show in Fig. 13 was obtained by fixing at the boundary a large displacement in a single step. Any attempt to try to follow a sensible quasi-static trajectory has had always miserable results. We do not have a reasonable answer to this but feel that these difficulties may be remedied by introducing dynamics, that is by alternating a variational step and a dynamical step during the discretized evolution.

PROPAGATION: KINKING OF A STRAIGHT CRACK IN MIXED MODE I&II
he propagation of a straight crack in a panel under mixed mode boundary conditions is considered. The geometries and boundary conditions considered are shown in Fig. 1 and Fig. 14, to which we refer for notations. Propagation of the crack in CLFM is still governed by the Griffith's criterion. Consensus is breached when addressing how a path is selected by the crack. In particular, competing criteria have been proposed for the kinking of a straight crack. In an isotropic setting, a well regarded criterion is the Principle of Local Symmetry, PLS-criterion, [14] which states that the crack always propagates in mode I, that is with inplane tractions that remain perpendicular to the crack in a small neighborhood of the crack tip [15]. The maximum tangential stress, MTS-criterion, (Erdogan and Sih [16], 1963) is the simplest criterion and it states that the direction of crack initiation coincides with the direction of the maximum tangential stress along a constant radius around the crack tip. Another possible alternative is the Gmax-Criterion which states that the crack will kink along a direction that maximizes the release rate of its potential energy among all kinking angles. It was shown [17] that the three criteria generically yield different kinking angles. There is in truth scant evidence that would support one criterion over another, and even less so in anti-plane shear because the crack is in mode III, so that the notion of mode I, or mode II propagation is rendered meaningless. In very recent paper, [18 ], Chambolle, Franfort & Marigo claim that the kinking of a straight crack is possible only if one admits discontinuity in time, that is the crack propagation is brutal. Based on their theorem there are only two possibilities to propagate a straight crack: 1. The crack kinks and the system jumps dynamically from a local minimum to another local minimum. 2. Crack branching is observed. In both cases CLFM does not work since the Griffith's criterion does not apply in his classical form, whilst variational fracture does. With the MTS-criterion the angle of kinking, in terms of the Stress Intensity Factor (SIF) ratio =K2/K1, is given by With the GMax-criterion the angle of kinking, in terms of the SIF ratio =K 2 /K 1 , is given by A plot of the kink angle as a function of the SIF ratio is reported in Fig. 16, where also a comparison with some experimental results, taken from the literature and collected by Al Bin Mousa in [13], is shown. Comparisons of the theoretical curves with experimental results (b). The graph (b) is taken from Al Bib Mousa [13] and "present work" in the legend (c) refers to [13].
The case represented in Fig. 1, with H=B=L, S H =(1) and V H =(1), is then solved with our computer model. The elastic domain is discretized with a non-structured triangulation, as shown in Fig. 17a where also the stress level is represented in a gray scale. The given displacement is discretized into 100 steps and the deformation at each step is determined by minimizing the energy with a descent algorithm. Denoting M the total number of triangles in the mesh and N the total number of nodes, the energy to be minimized is a function of (2) 3 M nodal displacements u 1 (n), u 2 (n) (n=1,2,..,3M) and of 2 N nodal positions in the reference configuration x 1 (m), x 2 (m) (m=1,2,..,N). The energy is a complex, non convex function of its arguments with a lot of local minima. Fig. 17a shows the deformation and stress at the onset of propagation. The propagation of the crack just before dynamical instabilization is shown in Fig. 17b. The kink angle predicted by our model is 30.46°. From the deformed state of Fig. 17 the SIF's K1, K2 can be derived approximately by computing the J-integral. The values obtained for K 1 , K 2 are K 1 = 1.3176, K 2 = 0.4288 For these values the kink angle predicted by MTS and GMax criteria can be computed from (7) and (8) We see that our analysis reproduces closely the GMax-criterion. The last results we present are concerned with the convergence analysis relative to the kinking angle assessment. The problem is essentially the same of Fig. 17, except for the signs of the given displacements: S H =(-1) and V H =(1). The "loading" history is again discretized into loading steps but the analysis is stopped as soon as the crack propagates. The domain  is discretized with structured meshes of decreasing mesh size h: h= L/12 2 , L/16 2 , L/20 2 , L/48 2 .
The results of these three analysis are reported in Figs. 18, 19, 20, 21. From these pictures we see that the kinking angle goes from 38.30° degrees for the coarser mesh to the value of 31.18° of the finest one. The value obtained converges toward the value predicted by GMax. Notice that a star shaped region of different size, centered at the crack tip, is always created during the four analyses. Such regions develop before crack propagation, that is in the elastic phase, and are formed by a fixed number of elongated triangles (6). A similar behaviour is detected also in the case of Fig. 17, also if it is less evident. Elongated triangles are usually hill-behaved in FEM analysis, but the point here is that such shapes are chosen through minimization and are the best possible in order to approximate the large stress that develops around the crack tip.