A Discontinuous Model to Study Fracture of Brittle Materials

A variety of computational techniques exist to describe the fracture behaviour of quasi-brittle materials. These models can be classified into two main groups: continuous and discontinuous models. In a continuous model, the displacement and strain fields remain continuous, even after a strong localization of the deformations. Localization of deformation can be triggered by strain softening. A major problem with classical continuum models is that the governing equations lose ellipticity for quasi-static problems and hyperbolicity for dynamic problems if strain softening is introduced. When using finite elements, this results in a strong sensitivity to the mesh. Upon mesh refinement, deformations localize into a band of zero thickness and complete structural failure can occur without dissipation of energy. To regularize the governing differential equations, non-locality or rate dependency can be introduced in the constitutive model. This controls the zone in which the deformations tend to localize. Examples of regularized continuum models are non-local damage models [1], gradient damage models [2], Cosserat continuum models [3] or viscous models [4, 5]. Discontinuous models represent cracks as displacement discontinuities. A discontinuous term can be incorporated into the strain field (weak discontinuity) [6–8] or into the displacement field (strong discontinuity) [9–18]. In this paper, a displacement discontinuity is introduced using a specific property of finite element shape functions. These form a partition of unity, which allows enhancing nodes with additional degrees of freedom. The first section covers the kinematics of a body crossed by a discontinuity. Then, the governing equations are derived. In the third section, implementation aspects are discussed and finally an example is treated.


Introduction
A variety of computational techniques exist to describe the fracture behaviour of quasi-brittle materials.These models can be classified into two main groups: continuous and discontinuous models.In a continuous model, the displacement and strain fields remain continuous, even after a strong localization of the deformations.Localization of deformation can be triggered by strain softening.A major problem with classical continuum models is that the governing equations lose ellipticity for quasi-static problems and hyperbolicity for dynamic problems if strain softening is introduced.When using finite elements, this results in a strong sensitivity to the mesh.Upon mesh refinement, deformations localize into a band of zero thickness and complete structural failure can occur without dissipation of energy.To regularize the governing differential equations, non-locality or rate dependency can be introduced in the constitutive model.This controls the zone in which the deformations tend to localize.Examples of regularized continuum models are non-local damage models [1], gradient damage models [2], Cosserat continuum models [3] or viscous models [4,5].
In this paper, a displacement discontinuity is introduced using a specific property of finite element shape functions.These form a partition of unity, which allows enhancing nodes with additional degrees of freedom.The first section covers the kinematics of a body crossed by a discontinuity.Then, the governing equations are derived.In the third section, implementation aspects are discussed and finally an example is treated.

Kinematics of a displacement jump
Consider a body W crossed by 2 non-intersecting discontinuities, G 1 and G 2 , as shown in Fig. 1.The displacement field is given by: where u is the displacement field of the body W, $ u and ũi are continuous fields and H G i is the Heaviside step function defined as: where W i + and W i -are subvolumes of W such that W W and the discontinuity G i is the border between the two subvolumes.The normal n i to the discontinuity is pointed towards Taking the symmetric gradient of the displacement field (1) results in the infinitesimal strain field: where d G i is the Dirac delta distribution centred on the discontinuity.

Partition of unity concept
A partition of unity is a set of functions N i , such that:

K. De Proft, W. P. De Wilde
In this paper, the partition of the unity property of finite element shape functions is used to introduce displacement discontinuities into finite elements.The discontinuous character of the displacement field is captured with the Heaviside step function.Using the partition of unity concept, the governing equation of the continuum and the discontinuity are separated and are consequently described by different constitutive laws.Inside the discontinuity, a plasticity based constitutive law is used to describe the decrease of tractions in function of the crack opening while the continuum is assumed to remain elastic.The methodology will be described and validated with a comparison between numerical simulations and experimental results.This paper is dedicated to J. Sejnoha, TU Prague, with respect and admiration for his scientific achievement.where n is the number of discrete points.Since finite element shape functions must be able to describe rigid body modes, they must fulfil Eq. ( 4) and consequently they form a partition of unity.Duarte and Oden [19] showed that a field such as the displacement field can be interpolated making use of these partitions of unity, where a i are the regular degrees of freedom, b ij are the enhanced degrees of freedom and g j is an enhanced basis with k basis terms.As can be seen, the classical degrees of freedom are enriched with additional degrees of freedom, when necessary.This enhancement can be done locally.Belytschko and Black [20] and Moës et al. [21] used the partition of unity concept to enrich standard approximations with near-tip asymptotic fields and a discontinuous function for elastic crack growth.Wells and Sluys [17] introduced the partition of unity technique for describing cohesive cracks in a material.
When the Heaviside function is chosen as an enhanced basis, the displacement field can be interpolated according to equation ( 5): where N is a matrix containing the finite element shape functions, a are the regular degrees of freedom and b i are the enhanced degrees of freedom related to discontinuity i.For each crack, a basis term and an additional set of degrees of freedom is added.

Governing finite element equations
The weak form of the virtual work equation without body forces reads: where h is taken from the set of admissible displacement variations and S e is the outer surface where external tractions t are applied.Using a Galerkin approach, the admissible displacement variations can be decomposed in the same manner as the actual displacement field.Inserting the kinematical expressions for multiple (in this case m) non-intersecting discontinuities, the virtual work equation can be rewritten as: Taking first variations of $ h( hi = 0) followed by variations of hi ($ h = 0), a set of m + 1 equations is obtained.
where T i are the traction forces working at the discontinuity G i .To obtain equation ( 9), the enhanced displacement field, ũ, is assumed to be zero where essential boundary conditions are imposed.From equations (9, 10), the set of discretized equations is obtained by introducing the discretized form of the displacement and the strain field.
In these equations, the continuum response and the discontinuity response are completely split.The continuum is assumed to remain elastic during the complete computation.The stress rate in the continuum can be easily obtained from the expression of the strain field in the bulk, given by Eq. ( 3): where C e is the elastic material tensor.The behaviour at the discontinuity is stated in terms of tractions and separations.
The separation of the discontinuity can be computed as: Eq. (13) shows that the separation of the discontinuity is expressed as a function of the enhanced degrees of freedom.The traction rate is defined as: where D is the material tangent for the discontinuity.This material tangent will be further elaborated in the next paragraph.
Eq. ( 11) is further linearized by inserting Eq. ( 12) and Eq. ( 14).After some mathematical manipulations, the linearized set of equations is given by where It can be seen from previous equations that the global system of equations (Eq.( 15)) remains symmetric when C e and D are symmetric.It should be noted that equations (16) are valid when only discontinuity j crosses the element, but the element is also influenced by discontinuity i.Note that all stiffness contributions in equations ( 16) are very similar.The crucial difference between the terms in equations ( 16) is the presence of the Heaviside function.This makes the finite element implementation relatively simple.

Integration of the crossed elements
From Eq. ( 16) it can be seen that the Heaviside step functions should be taken into account for the integration of the stiffness matrix.This means that the integration must only be performed on the side of the element where the Heaviside function equals one, i.e.W i + .
Obviously, enough integration points must be present on that side so that an adequate integration is performed.Since a discontinuity can cross the element arbitrarily, the safest solution is to redefine the integration scheme.When only one discontinuity crosses the considered element, the integration rule is adapted following Wells and Sluys [17].For the triangular quadratic element as the underlying finite element, 23 integration points -21 in the continuum and 2 for the discontinuity -are inserted as shown in Fig. 2a.When 2 non-intersection discontinuities cross the same element, the integration rule is re-adapted again as shown in Fig. 2b.In this case, 15 integration points for the continuum and 4 integration points for the discontinuity are used.Of course, depending on the position of the discontinuities, the integration rule might change.

Enhanced nodes
Another implementation issue is the enhancement of the nodes.Only nodes whose support is crossed by a discontinuity should be enhanced.Furthermore, the enhanced degrees of freedom of the nodes on the support of the crack tip remain constrained.Consequently, the separation of the crack tip is zero.An overview of which nodes should be enhanced is given in Fig. 3.The enhanced nodes are represented by squares, while circles represent the nodes at the support of the crack tip.
Since the crack can freely run through the finite elements, it is possible that a discontinuity runs close to a node.As a result, a small proportion of the support of the node lies in either W i + or W i -.In this case, the global stiffness matrix is not necessarily well conditioned.For this reason, an extra condition is introduced: where W s is the volume of the support of a node.The tolerance is dependent on the precision of the solver.When condition ( 17) is not met, the considered node is constrained.The influence will be spread out over the other nodes.

Initiation and propagation rules
An important issue is the initiation of the discontinuity.A criterion is needed to decide when a discontinuity should be initiated and should propagate.A typical example of an initiation criterion is a plasticity yield surface: whenever the stress state in one integration point of the considered element is situated outside the elastic domain bounded by the yield surface, a discontinuity is initiated, as will be explained in section 7.
Another important item is the direction of the discontinuity.
Oliver [16] stated that the direction of the discontinuity could be found from the stress state at the moment of initiation.When using a plasticity yield surface as the initiation criterion, the direction can be obtained by means of a bifurcation analysis.Furthermore, a discontinuity can only grow from a previous crack tip to introduce path continuity, and discontinuities can only grow at the end of a time step.In this case, no discontinuities are introduced in non-equilibrium states and the quadratic convergence of the Newton-Raphson process is ensured [17].A discontinuity crosses each time through the whole element.When allowing multiple discontinuities to grow without intersection, specific interaction rules should be defined.Different cases can be considered, namely a) only one crack tip touches the considered element, b) two crack tips touch the considered element, c) three crack tips touch the considered element, d) one crack tip touches an element that is already crossed.The first case is the most common case.The stress state in the considered element is checked and the initiation criterion decides whether the discontinuity should propagate along the computed direction or not.
The second case implies several possibilities: l the two crack tips link and form one crack, l one crack propagates while the other stops, l both cracks grow, without intersecting.
To decide which case occurs, the normal n comp is computed according to the present stress state in the considered element.
Then the normal n con to the connecting line, i.e. the line that connects the crack tips, is computed.From this normal, a The behaviour within the discontinuity is described by a plasticity based cohesive zone model.The adopted plasticity model was proposed by Carol et al. [22] for use in interface elements.Consequently, the plastic yield function is given in the traction space instead of the stress space.A hyperbolic yield surface is introduced, where , are the normal and tangential component of the traction vector, c is the cohesion, f t the tensile strength and f the internal friction angle of the material.For tension, an associative flow rule is adopted.The evolution of the yield surface is governed by the decrease of tensile strength and cohesion throughout the computation: where f t0 and c 0 are the initial values for the tensile strength and the cohesion, G fI is the mode-I fracture energy, G fII is the mode-II fracture energy and W cr is the energy dissipated during fracture processes.The energy is defined as: ponent of the plastic separation vector.The decrease of tensile strength and cohesion is coupled to the energy dissipated during the fracture processes.Moreover, the choice of Eq. (19a)/(19b) ensures that the total mode-I fracture energy/mode-II fracture energy is dissipated when the tensile strength/cohesion vanishes.Furthermore, the decrease of tensile strength and cohesion is coupled: when a material is damaged due to tensile loading, the tensile strength but also the cohesion decreases.
The tangential stiffness and the stress update are obtained with classical elasto-plastic equations.The tractions are defined through: The plastic deformation rate is defined as: where & l is the plastic multiplier rate.For tension, an associative flow rule is adopted.The plastic deformation rate can be introduced in Eq. ( 22): The plastic multiplier rate can be obtained through the consistency equation: where h is the hardening/softening modulus.Inserting the result for the plastic multiplier into equation ( 23 The tangential stiffness can be inserted in the finite element equations (see Eq. 14).
The elastic stiffness is chosen very high, in theory infinite, in order to suppress the artificial elastic part of the solution.Since a discontinuity is only inserted when the yield surface is violated, the jump is completely inelastic.

Numerical example
Nooru-Mohamed [23] examined double edge notched specimens of different sizes (200×200×50, 100×100×50, 50×50×50 mm), subjected to different loading paths.All specimens were placed in a special loading frame, allowing a combination of shear and tensile loading.For the numerical simulations presented in this section, only one loading path is considered (path 4 in the experiments).
Double-edge notched specimens, shown in Fig. 4, are first loaded in shear until a certain value of the lateral force P s is reached.Afterwards, the lateral shear force is kept constant Acta Polytechnica Vol.44 No. 5 -6/2004 tolerance is computed.The obtained limit values for the normal are n tol1 and n tol2 .When the normal computed from the stress field lies within the zone defined by n tol1 and n tol2 , the cracks are connected.In the other case, one crack propagates along the normal computed from the stress field while the other crack is stopped.The propagating crack can be chosen arbitrarily.Another possibility is to allow the crack with the highest energy dissipation (main crack) to propagate.Another option is to let both cracks propagate.In this case, the cracks may or may not intersect.The case of intersecting cracks is not considered in this paper.The third case can be solved in the same way, only more linking possibilities should be considered here.The last case is simply solved by arresting the crack tip.It would also be possible to let the crack grow, but only if the discontinuities inside the element do not intersect.This is however not considered.It is clear that the defined interaction rules are more or less arbitrary.Nevertheless, the rules are straightforward and relatively easy to implement.If necessary, the adopted rules will be refined.l c 0 = 20 MPa, l G fII = 0.1 N/mm, l tan f = 0.5.
For the continuum, the material paramters are : Young's modulus E = 25000 MPa and Poisson's ratio n = 0.2.
The obtained load-deformation curve is shown in Fig. 6a.When compared with the experimental curve, a good agreement is found.Especially the peak load is simulated remarkably well.The post peak response in the finite element simulation is more brittle.The computed crack path is compared with the experimentally obtained path.As can be seen in Fig. 6b, the computed crack path is in good agreement with the experimentally observed crack path.
During the experiments, Nooru-Mohamed also connected LVDTs to the specimen in order to study local deformations.The position of these additional LVDTs is visualized in Fig. 7.The recorded deformations are plotted versus the shear deformation d s , and compared with the computed values in Fig. 8a-e.Examining Fig. 8a and Fig. 8e, the calculated deformations for LVDT 1 and LVDT 5 show a good agreement with the measured values.The deformation for LVDT 2, Fig. 8b, is not captured with the calculations.The computed crack path runs outside the measuring range of LVDT 2. As a result, the experimentally observed increase in deformations is not observed in the computations.For LVDT 3 and LVDT 4, the measured and calculated deformations show the same tendency.First a small increase is noticed.When the crack has passed the location of the LVDT, the deformations start to decrease because the crack passes outside the measuring range of the LVDT.
The overview in Fig. 8a-e shows that apart from the load deformation curve and the final crack path, also local information is captured in a reasonable way.
Finally, the same model and model parameters are used to simulate load case (a).In this case, the lateral force is increased until P s = 5 kN.Then, the lateral force is kept constant and a tensile load is applied.The comparison of the computed load deformation curve with the experimental obtained curve is shown in Fig. 9a while the experimental and computed crack path are visualized in Fig. 9b.
Again, the peak load is captured remarkably well.The computed post peak response is slightly more brittle.When the crack path is studied, it can be observed that the crack which grows from the left notch is in agreement with the experimentally observed crack.The crack growing from the right notch is more curved than the experimentally observed crack.

Fig. 4 :
Fig.4: Geometry of the specimen (all dimensions in mm)