A level set based model for damage growth: The thick level set approach

In this paper, we introduce a new way to model damage growth in solids. A level set is used to separate the undamaged zone from the damaged zone. In the damaged zone, the damage variable is an explicit function of the level set. This function is a parameter of the model. Beyond a critical length, we assume the material to be totally damaged, thus allowing a straightforward transition to fracture. The damage growth is expressed as a level set propagation. The configurational force driving the damage front is non‐local in the sense that it averages information over the thickness in the wake of the front. The computational and theoretical advantages of the new damage model are stressed. Numerical examples demonstrate the capability of the new model to initiate cracks and propagate them even in complex topological patterns (branching and merging for instance). Copyright © 2010 John Wiley & Sons, Ltd.


INTRODUCTION
It is well known that fracture mechanics alone may not model the full scenario of the degradation of solids under mechanical loading. Crack initiation for instance requires damage mechanics to model the gradual loss of stiffness in a small area [1,2]. When it comes to modelling damage, special care needs to be taken to avoid so-called spurious localizations. Several damage models were proposed in the literature to avoid spurious localization as detailed in [3,4]: • non-local integral damage model: the damage evolution is governed by a driving force which is non-local, i.e. it is the average of the local driving force over some regions [5,6]; • higher order, kinematically based, gradient models through the inclusion of higher order deformation gradient [7][8][9] or additional rotational degrees of freedom [10]; • higher order, damage based, gradient models: the gradient of the damage is a variable as well as the damage itself. This leads to a second-order operator acting on the damage [11][12][13].
Non-local and higher order gradient approaches were compared in [14]. More recently two strategies to avoid spurious localization also appeared. The phase-field approach emanating from the physics community [15][16][17] and the so-called variational approach [18][19][20].
In this paper, we wish to model damage as a propagating level set front. Let us first tour the existing papers using a level set approach to model a degradation front. In report [21] a level 360 N. MOËS ET AL. this zone the crack (although we shall see it is not necessarily of zero thickness). This zone may be of very complex topology, thus dealing easily with branching and merging. In our approach, we thus do not need a transition to a cohesive crack to handle full degradation as in [28][29][30]. The crack appears as a consequence of the damage front motion. Finally, note that in the presented approach, the fully damaged zone is clearly delimited by a level set. It is thus not a set of elements similar to the element deletion approach [31]. Moreover, this zone is still allowed to carry compressive loading if needed.

A NEW LEVEL SET BASED DAMAGE MODEL
It is well known that purely local damage models suffer from spurious localization because no length scale is introduced in the models. In this paper, we introduce a length l c in the wake of the front creating what we call a thick level set (TLS). The damage depends on a signed distance function to the front. This description induces a specific evolution law for the damage, the evolution follows the motion of the level set. The damage itself could be regarded as a level set, but to preserve the generality of the approach, we choose to separate the evolution of the geometry described by the level set and the associated value of damage.
In order to introduce the TLS approach we consider a simple damage model: an elasto-damage model with a scalar damage parameter.

A simple local damage model
The free energy per unit volume depends on the strain and a scalar damage variable denoted by d The state laws are obtained by differentiating the free energy where r is the stress and Y the local energy release rate. For instance, if we consider a symmetric behavior in traction and compression, we may use the potential leading to A potential with closure effect will also be considered later in the numerical simulations. The evolution of the damage is given by a dissipation potential * (Y ) which is a convex function of Y : Alternatively, we may write The above potentials (ḋ) and * (Y ) are dual and related by the Legendre-Fenchel transformation For non-regular functions, the definition of the gradient is replaced by the notion of sub-gradient for all admissible d * . A similar definition is associated to * . Note that from the above we have the following inequality for any couple (Y,ḋ) satisfying or not the damage evolution law * (Y )+ (ḋ)−Yḋ 0 The equality is satisfied if and only if the damage evolution law is satisfied too: Note that the convexity of or * ensure positive dissipation, this is due to (8) considering the properties (ḋ) convex, (ḋ) 0, (0) = 0 (11) The damage model described above suffers from spurious localizations. To regularize this situation, we introduce a length scale in a different way than that proposed in the literature.

Free energy and state laws for the level set based damage model
The damage description is depicted in Figure 1. The level set = 0 separates the domain into an undamaged zone and damaged one. In the damaged zone, the damage variable d is an explicit function of the level set. Note that a single level set is needed even if damaged zones are disconnected as shown in Figure 1. The damage increases progressively as the level set value rises. Mathematically this is expressed by the following inequalities: d () 0, 0 l c (13) d() = 1, l c (14) 362 N. MOËS ET AL.
The function d = d() is assumed to be continuous for clarity. Discontinuous functions can be chosen in the proposed approach but the obtained expressions must be modified with additional contributions. At a distance l c from the damage front, the material is assumed completely damaged. A spatial relationship between the damage value along the normal to the front = 0 is then deduced. Consider a body . Along the external boundary * , loading T d (t) is applied on * T and displacements u d (t) are prescribed on the complementary part * u ; * = * T ∪* u , * T ∩* u = ∅. Some parts of may be totally damaged i.e. d = 1. We denote by c , the part of not fully damaged: Potential energy is written as a function of the displacement field and the level set function : where (u) denotes the strain, i.e. symmetric part of the gradient of u. The displacement field must satisfy the kinematic boundary conditions whereas the gradient of the level set field must be of norm 1 since it is a signed distance function. The function d() is assumed to be continuous. The displacement, the free energy and the stress vector are then continuous functions along each iso-curve . Otherwise, a more complex expression involving jump terms must be written. Along the displacement is continuous, then The set of admissible variations is then defined by The condition on indicates that to remain a signed distance the level set variation must be constant if we follow a path along the gradient of . In order to visualize this property, let us consider a plane situation as depicted in Figure 2 with the band of material comprised between the boundaries 0 and c . This band may be parametrized by two curvilinear coordinates and s, where s is the curvilinear abscissa along 0 . The iso-coordinate curves are orthogonal to each other. The condition (23) implies that is a function of s only and must be seen as a normal perturbation of the front location. We may then write = a(s) Note that the level set (x, t) = 0 defines a surface of normal n then the velocity of this surface is v n such that Finally, note that in the curvilinear system of coordinates (, s), the surface element d is expressed as where (, s) is the radius of curvature at (, s) of the iso-curve. The last equality comes from the fact that the radius of curvature is linear along the iso-coordinates and we have used the shortened notation (s) = (0, s). Equation (26) is well defined provided is differentiable. In the 3D case, similar formula may be derived. The directional derivative of the energy E with respect to ( u, a) ∈ A 0 may now be computed. We take into account the state laws (2) Zeroing the directional derivative of the potential energy with respect to ( u, 0) ∈ A 0 d , we obtain the equilibrium equation Let us now analyze the variation of the potential energy with respect to . First, we note that d () is only non-zero in the thick band oc in which the level set varies from 0 to l c Variation of the potential energy with respect to is reduced to:  This relation defines the driving force associated with the moving band. The term is quite similar to the energy release rate introduced in [22,23]. The term g(s) can also be interpreted as a configurational force per unit length on the front In the above, l is either l c or some smaller values as illustrated in Figure 3. The configurational force (33) is an average of the local damage driving force Y weighted by d () and the evolution of the front curvature along the thickness. Note that the equality (33) bears some similarity with the transition from damage to fracture as proposed in [32], whereas here it is relating diffuse damage to a degradation front advance. Finally, let us compute the dissipation. To do so, the virtual front advance a is replaced by the true front velocity v n . And therefore the dissipation reduces to: The last equation defines the global dissipation in terms of an energy release rate along a moving surface. The link between this global formulation and the local variables Y,ḋ inside the front is essentially due to Equation (25).

Dissipation potentials and evolution laws
We now need to provide the evolution law relating the front speed v n to the configurational force g.
Either we give this relationship in a phenomenological way prescribing directly a relation between v n and g or we somehow homogenize the local damage evolution law (5). For this, we can use the notion of a shock-generating function introduced in [33,34] for an elastic continuum and hence we obtain a global constitutive relation as in plasticity [35].
We propose here to satisfy as much as possible (5) while fulfilling the constrainṫ This relation is defined for any distribution v n along 0 . In general situation, we cannot fulfill at every point the local damage evolution law since for a given Y field, the damage evolution given by (5) will not satisfy the damage constraint (35). We also note from Equation (35) that the front velocity must be positive at all time and all point on the front to prevent a reduction in damage.
Given a field Y , we thus minimize the expression below which may be interpreted as the error on the local constitutive law. The value of the infimum will not be zero and somehow measures the effect of homogenizing the evolution law over the 365 band thickness. Noticing that v n is present in the last two terms of (36) the search of the infimum reads as The above integral may be limited to oc . Using (32) we may also write and we define the homogenized potential through Note that the homogenized potential expression depends explicitly on the location on the front. This potential is defined for any distribution v n given along 0 . It thus ensures a local definition of the potential at point s: The infimum (37) where the dual potential * is obtained through the Legendre-Fenchel transformation of . The local damage evolution law (10) has been homogenized into a non-local evolution law v n = * * (g, s) Note that through the homogenization process, the properties of the local potential (11) are preserved at the homogenized level. The following properties are satisfied at any point on 0 ensuring the positivity of the dissipation gv n along the front. This automatic preservation of duality when creating the non-local model is an advantage over classical non-local models [36]. A typical case. Let us apply this procedure to a power law-type damage evolution law with thresholdḋ with the positive · + (negative · − ) part of a scalar following the usual definition The corresponding potentials are where n = 1/n and I + is the indicator function of R + The homogenized potential is given by (39) ( where The dual potentials and the corresponding evolution law are where we did denote r as part of which may be differentiated in the classical sense. Note that in practice, as it will be seen in the numerical implementation, it is not necessary to obtain the analytical global potentials .

Advantages of the proposed approach
• Compared with non-local approaches in which the duality may be corrupted [36], the proposed approach is thermodynamically sound. • Although the approach is non-local the specific numerical treatment for non-locality is limited to a band. • The scenario of the damage variable getting to 1 is included in the model. • The criterion is maximum at the tip of a crack and not a little away from it as in classical non-local approach [37]. • The issue of the proper boundary conditions in the gradient type is removed. Also, the displacement field is no longer forced to be extra smooth as in the higher order displacement gradient models.

Relationship with fracture mechanics
Let us now consider that the damage front has a shape depicted in Figure 4 and that it progresses so that it keeps its shape with a speedȧ. This means that the normal velocity on the front is where n is the outer normal to the front and e the crack direction. Substituting this representation of v n into (53), (55), (56), we get the following potentials: where To better see the relationship between damage parameters (Y c ,l c , k) and fracture mechanics parameters (G c , K ), let us consider the particular case of a circular damage front and a linear damage to level set relationship. In this case, we get

SPACE AND TIME DISCRETIZATION
In this section, we consider only one front 0 . In the presence of multiple, disconnected fronts i, the same procedure is performed for each individual front i 0 . Regarding the space discretization, we use the X-FEM. This method is well known to handle cracks [38,39] without explicit meshing of the crack. The X-FEM may also handle voids [40] and it is this feature that we will be using to model completely damaged zones. Thus, the boundary c does not need to be meshed explicitly. For a given location of the damage front, once the mechanical fields are evaluated, we need to compute the front velocity. We shall look for the The modes F I are explicit functions of s initially defined on the front 0 , then extended through the band oc . The extension is performed just as a velocity is extended from the iso-zero of the level set when dealing with the level set equation [41,42]. In order to find the velocity coefficients, the discretized expression of the velocity is introduced in (37) inf The integration is restricted over oc since outside this band d () = 0. In the minimization, the velocity v n must be positive. If the velocity is not positive the potential is infinite. In order to impose the positivity of the velocity, introduction of a Lagrange multiplier is useful.
where r is the regular part of the potential which is no longer infinite when the velocity is negative. Note that in the above equation, the term also involves the gradient of the damage. This is important if we do not want to loose the constraint of positive velocity as l c goes to zero. The interpretation of the Lagrange multiplier is as follows: it will be zero in zones where the front is moving and related to the gap Y c −Y in areas where it is not moving. Finally, in order to remove the constraint 0 in (68) the trick, introduced for contact by Ben Dhia et al. [43], is to write where >0 is a parameter that does not affect the solution. Differentiating with respect to v n and , we get in which we did use the notation This system may be solved using a Newton-Raphson technique. Two non-linearities appear. One coming from the constitutive law r and the second coming from the contact condition . It was shown in [43,44] that solving one difficulty at a time was efficient. Basically, there are two embedded loops in the Newton-Raphson procedure.
We now discretize where the modes F I depend only on the curvilinear abscissa s ⎡ ⎢ ⎣ where 0 = ( 0 , v 0 n ). We note that the velocity computation on the front does not require an explicit computation of the configurational forces on the front.
Regarding the loading, different types of control may be used. In the numerical experiments we did use a crack mouth opening displacement (CMOD) for the three-point bending test and a dissipation control-type algorithm for the other experiments. For the dissipation control algorithm we were inspired by the work of Verhoosel et al. [45]. In our case, the dissipation reads as It is tempting to impose the value of D by some prescribed constants. Unfortunately at initiation, the left term goes to zero and the equality may only be satisfied for infinite loading. Hence, we introduce a given parameter and we take into account the volume of the damaged zone by the constraint This constraint is imposed through a Lagrange multiplier . Taking into account the dissipation constraint, the inf-sup principle (69) now reads as where Y 0 is the value of Y for a unit loading. It may be observed in the above that the meaning of the Lagrange parameter is simply the square of the loading. The corresponding Newton-Raphson iteration is given below ⎡ Once the normal velocity of the level set is known, the level set may be updated using classical level set approaches [41,42].

NUMERICAL SIMULATIONS FOR SEVERAL 2D CASES
In order to take into account more general loading cases and the closure effect in the fully damaged zone, we consider an asymmetric tension/compression damage potential with a scalar damage variable. We have considered a strain-based model. Other types of models involving closure effects may be found in [46] ( , d) = (1−d) + : where and are the Lamé elastic coefficients and h 1. The positive (negative) part of the tensor is applied to its spectral decomposition where i , i = 1, 2, 3, are the eigenvalues of the strain tensor and e i , i = 1, 2, 3, the eigenvectors. The potential is expressed as follows using the eigenvalues (and plane strain assumption 3 = 0)) One can easily show that the eigenvectors of the stress and strain tensors coincide. The eigen stresses are obtained from the eigen strains by differentiating (91) Storing the two eigen stresses inr and the two eigen strains in˜   It is important for the Newton-Raphson algorithm to compute also the tangent stiffness. First, we note that the stress is the derivative of the potential finally, the local energy release Y is given by We now discuss some mathematical aspects. The potential is convex with respect to the strain for any damage value (d 1). It is differentiable in the classical sense everywhere, i.e. for any given strain, the corresponding stress is unique. However, the tangent stiffness H, i.e. second derivative of the potential is discontinuous for 1 = 0 or 2 = 0 or 1 + 2 = 0. This can be seen from the expression of the D operator.

RESULTS
We present here four test cases to validate and illustrate the applicability of the TLS approach.
The first test case is the three-point bending problem with an imposed CMOD. Both the mesh sensitivity and the speed effect, i.e. modifications of the load-displacement curves depending on the  ratio between the crack mouth opening velocity and the damage front velocity [47], are observed and the results are compared with a reference crack solution. The three remaining test cases illustrate the advantages of the method by considering the coalescence, branching and merging of cracks with an imposed dissipation rate to represent severe snap-backs on the load-displacement curves.
All computations are performed using a linear relationship d() and the non-linear potential (88). Note that this non-linearity does not have much impact on the solution for the three-point bending problem since the damage zone essentially experiences traction. The dissipation potential used is the one corresponding to the damage evolution law (44), using n = 1. The velocity v n is discretized using polynomial (Hermite) modes providing a partition of unity. Each mode is defined (i.e. non-zero) on a length of about eight elements in the curvilinear coordinate along the front. The total number of modes is therefore increasing as the front length increases. For instance, 50 modes are used at the end of the coalescence test case (Figure 11).
Before the damage front may be grown it must be initiated. It is initiated at the location and for the load for which the critical local energy release Y c is reached. A level set is placed with a radius of about two element sizes.

Three-point bending
We consider the beam in the three-point bending configuration depicted in Figure 5. Crack growth for this specimen has been studied by Carpinteri and Colombo [48] using finite elements. These results will be used as the reference results for the comparison. The material characteristics are the Young modulus E = 36.5×10 9 Nm −2 and the Poisson ratio of = 0.1. Other parameters are chosen to match approximately the order of magnitude of the reference solution: Y c = 100 Nm −2 , l c = 7.5×10 −3 m and k = 2.5×10 −2 s −1 while we ensure a constant crack mouth opening velocityU of 2×10 −7 ms −1 . We therefore impose in this test case a CMOD instead of a given dissipation rate.
First, we consider three different grid refinement levels to observe the mesh sensitivity. Figure 6 is a close-up view of the three displacement fields and damaged zones obtained for the same physical time t = 103.25 s on three embedded regular meshes with an element size h =l c /3,l c /6 and l c /12. The corresponding load-displacement curves are depicted in Figure 7 for both the deflection and the CMOD. The same order of magnitude is reached for the three computations and the reference solution (continuous line). In Figure 8 is depicted the time evolution of the crack and the norm of the displacement field. The major loss of stiffness of the specimen happens between states B and C.
Then, a family of load-displacement curves can be obtained with the variation of the speed effect non-dimensional numberU kl c representing the ratio between the loading velocity and the front velocity v n ∝ kl c . We consider here variations of k to observe the impact of the speed effect. Three load-displacement curves are depicted in Figure 9 for k = 2.  with the observations in [47]: the increase in the speed effect tends to scale the curves. Indeed, for a same damaged zone, the specimen stiffness F/U remains the same, but larger load and displacement are reached.

Crack coalescence
Let us consider the configuration depicted in Figure 10. The crack coalescence is a standard benchmark studied extensively in, for instance [49,50]. We consider the ratios a/b = 0.2, e/b = 0.025 as in [50], with b = 40 mm, E = 200 GPa and = 0.3. Other parameters are l c = 1 mm, Y c = 1 kPa and k = 2.5×10 −6 s −1 , with an imposed dissipation rate of 1 Wm −2 . The load-displacement curve is depicted in Figure 10. The stress and displacement fields of points A, B, C and D can be observed in Figure 11. Oscillations appearing on the load-displacement plot seems to be due to the control in dissipation and will be investigated in a future work. Note however that small dissipation increments have been used in this computation, but there is no particular restriction on these imposed increments. We eventually obtain the expected crack growth pattern of complete failure (state D) in good agreement with experimental observations [49,50].

Crack branching
This test case consists in a specimen made of two materials presenting different Young's moduli. A crack is initiated in the material presenting the lower modulus, as depicted in Figure 12. The damage evolution is presented in Figure 13. The crack propagates normally until approaching the Young modulus transition zone (point A). Then, since branching requires less energy than continuing into the second material, the fully damaged zone becomes larger (point B) and the crack eventually branches into two parts (point C) as expected. This branching process naturally appears as the level set grows, and would clearly be much more difficult to represent with classical models.

Crack percolation
The last test case configuration is depicted in Figure 14. We consider a random distribution of small defects with a maximum initial size of 1 mm with l c = 3.5 mm and the specimen length b = 40 mm. We choose Y c = 100 Pa, k = 2.5×10 −4 s −1 and an imposed dissipation rate of 100 Wm −2 , for the Young modulus and the Poisson ratio of 37.7 GPa and 0.2, respectively. Note that crack percolation was also studied using the X-FEM in [51].
Comparing Figures 14 and 15, we observe that the specimen stiffness is not much affected when the defects grow, until they reach a maximum damage of d = 1. Then, cracks appear as shown in point B, and the stiffness decreases much more quickly. We also observe the merging of damage zones and cracks. Unlike classical fracture mechanics methods, this process appears naturally as the level set function grows, and the cracks merging does not require any particular implementation. Moreover, it is not necessary to explicitly determine which crack should grow. Indeed, the front velocity is already computed as the velocity minimizing the energy while ensuring a given dissipation rate.

CONCLUSIONS
In this paper, we have introduced a new way to model damage growth in solids. The model is based on a so-called Thick Level Set (TLS). The damage variable is tied to the distance to the damage front (iso-zero of the level set). Beyond a critical distance the material is considered completely damaged (but may still sustain compression if it is prescribed in the model). The model thus introduces a length scale and it regularizes local damage models. A set of numerical experiments did demonstrate the capability of the model to handle initiation, growth, branching and coalescence of crack-like patterns.