Locking‐free interface failure modeling by a cohesive discontinuous Galerkin method for matching and nonmatching meshes

In this work, modeling of brittle failure of the interface for a linear elastic material is presented. The idea is to integrate a novel extrinsic cohesive zone model into the incomplete interior penalty Galerkin variant of the discontinuous Galerkin (DG) method. As a result, the initial stiffness in the prefailure regime is omitted without having to remesh the crack path during the crack propagation. The interface model is used in combination with different discretization techniques, including matching and nonmatching meshes. This is possible due to the DG method's weak continuity constraint. Moreover, the locking problem in the bulk is cured by the application of a reduced Gaussian integration scheme on the boundary terms. The performance of the new cohesive discontinuous Galerkin elements with different integration schemes is compared with one of the standard intrinsic cohesive models. Due to the elimination of locking, crack initiation at the interface can be realistically displayed.

computation of the material tangent. 15 Over the years, several remedies were proposed to overcome these problems. 16,17 Apart from the numerical problems, a proper definition of the traction-separation relation plays an important role. 18,19 In general, thermodynamical consistency 20 and satisfaction of the main balance laws 21 have to be guaranteed.
Among the different types of CZ models, two main approaches mainstream this field, namely, the intrinsic and the extrinsic methods. Although both manifest the idea of the TSL, there are main differences in terms of the form of the TSL as well as the implementation. In the intrinsic CZ method, elements are inserted between the bulk elements prior to failure. Consequently, the intrinsic CZ model is ideal for parallel computations due to the uncoupled structure of its stiffness matrix. Regarding the form of the TSL, an initial elastic regime results in an artificial compliance (see Figure 1A). Increasing the initial slope may help to reduce this effect. Nonetheless, exaggeerating it can bring about time step problems and ill-conditioning of the stiffness matrix. 22 In contrast, the extrinsic approach ( Figure 1B) excludes the initial elastic regime and is thus more robust especially in dynamic calculations. 23 On the other hand, the parallelization of this method remains a cumbersome task since the mesh structure at the tip of the crack is steadily changing. 24 This is due to the fact that extrinsic CZ elements are inserted into the conforming elements once the failure criterion is met and therefore a remeshing of the fracture process zone is required.
By use of nonconforming FE methods such as the discontinuous Galerkin method (DG), one can benefit from the advantages of the aforementioned CZ approaches simultaneously. Most DG methods permit weak discontinuities. Consequently, DG elements no longer share mutual nodes. Hence, the extrinsic CZ elements can be implemented between the bulk elements prior to the failure and thereafter no regeneration of the mesh is needed. Thanks to the structure of the DG stiffness matrix, parallelization of the method is no longer an obstacle.
DG methods were first introduced by Reed and Hill 25 for the neutron transport problem. A penalty term was added at the element boundaries to stabilize the solution of the discrete problem (see Nitsche 26 ). After its initialization, the DG method was used in the context of different types of partial differential equations including hyperbolic, near hyperbolic and fourth-order problems (see References 27,28). Besides applications in fluid dynamics, 29-31 a number of interesting applications in solid mechanics are found (eg, References 32,33).
One of the first applications of the DG methods in the failure modeling of the interface is seen in the work of Mergheim et al. 34 The hybrid DG/CZ method of Mergheim et al 34 was further developed in the context of composite materials by Prechatel et al 35 and Wu et al 36 to model the debonding of the matrix-fiber interface. Dynamic fracture and fragmentation of solids in three dimensions is presented in the work of Radovitzky et al 23 by applying the DG/CZ model. A detailed implementation of this method is discussed by Nguyen 22 for both elastic and inelastic solids undergoing finite deformation. Versino et al 37 modeled the failure of the multilayered composite shells applying the DG/CZ method. Shear locking is alleviated in their work by application of the assumed natural strain and enhanced assumed strain techniques. Nevertheless, a comprehensive study on the effect of both, shear and volumetric locking, on the interface behavior using DG/CZ method is not conducted in the literature.
DG methods have also been applied in combination with other failure models than CZ approach. For instance, DG methods show good potential in addressing some stability issues when they are combined with other methods to predict failure and damage. 38 Aduloju and Truster 39 developed a DG formulation for modeling dynamic debonding in composite materials. Ghosh et al 40 generalized Nitsche's method to enforce stiff anisotropic cohesive laws. In this way, the problem of oscillations due to spurious tractions could be overcome. Bird et al 41 combined the idea of configurational forces 42 with a DG method to predict arbitrary crack propagation. The transition from diffuse damage to sharp discontinuities is an interesting and challenging topic that helps to model the entire ductile fracture process. Leclerc et al 43  advantages of a nonlocal damage model with a cohesive model in a DG finite element framework to address this problem. Furthermore, the simultaneous evolution of plastic deformation and damage has to be considered for a proper ductile damage model. 1,44 DG formulations are also applied to discretize phase-field fracture propagation 45 or gradient-enhanced schemes. 46 Conforming methods impose a strong continuity constraint on the discretized body. As a consequence, the meshes on opposing sides of the crack must be conforming. This leads to unnecessarily high CPU time effort. The application of the DG method facilitates the incorporation of hanging nodes into the discretization. Nonmatching meshes can be advantageous when it comes to, for example, heterogeneous interfaces. A novel node to segment/node to surface CZ element was developed by Paggi and Wriggers 47 to address nonmatching meshes at the interface. Nguyen et al 48 applied a similar approach in combination with the DG method.
The coupling of the DG method with the CZ model leads to benefits such as ease of implementation as well as increased robustness. Nevertheless, the artificial stiffening of the bulk has to be diminished prior to and after the occurrence of failure. There exist already some variants of the DG methods which do not show locking, for example, an incomplete interior penalty Galerkin (IIPG) variant combined with a reduced integration scheme can alleviate both volumetric and shear locking phenomena. 49,50 Grieshaber et al 51 also showed that, specifically in the case of quadrilateral, bilinear elements, underintegration of selected edge terms in their new interior penalty method results in elimination of volumetric locking. Wulfinghoff et al 52 introduced a hybrid DG (HDG) method for geometrical nonlinearity which is free of locking as well. The latter was improved and compared with a locking-free element formulation 53 in the work of Reese et al. 54 In Reference 55, the concept of control points was introduced into an HDG method and applied to crystal plasticity.
In the present work, a novel extrinsic CZ model with contact, viscosity as well as cyclic loading formulations is integrated into the IIPG version of the DG method. Based on Reference 23, a simple implementation of the cohesive DG method (CDG) is applied. The discretization in case of matching meshes is performed by means of four-node elements. Hanging nodes (nonmatching meshes) have been treated by the application of three-node elements from References 47,48. In addition, the influence of locking on the fracture behavior of the material is studied. To the best of the authors' knowledge, it is the first time that this issue is investigated. The problem is cured by the new locking-free CDG method. This paper is structured as follows: First, the mathematical background of the model including the strong form, weak form, and discretization is clarified. Later, different numerical integration schemes are illustrated. Several validation tests in addition to five benchmark examples are considered to investigate the performance of the CDG method. Finally, the work is summarized as well as conclusions are drawn.

FORMULATIONS
Within this section, the unsymmetric IIPG variant of the DG methods as well as our new CZ model are presented. In the following, the strong form along with its unified weak form of the linear momentum balance is introduced, followed by the discretization techniques for matching and nonmatching meshes. Next, the linearization of the residuals will be explained. Finally, numerical integration schemes are clarified and discussed.

Strong form
The body Ω with the outer boundary Ω is considered as shown in Figure 2. The outer boundary Ω is divided into Ω t and Ω u with the prescribed traction t p and the prescribed displacement u p acting on them, respectively. In general, the relations Ω u ∪ Ω t = Ω and Ω u ∩ Ω t = 0 hold for Ω. Within the body Ω, the discontinuities Γ divide the body into a finite number N e of subdomains Ω e . The relations Ω = ∪ N e e=1 Ω e and Γ = ∪ N e e=1 Ω e hold for the union of the subdomains and their boundaries, where e is the subdomain index. On each side of the discontinuity Γ, the subdomains are expressed by − and + with their corresponding outer normal vectors n − and n + , respectively. In case of weak discontinuities, in order to simplify, we define the unique normal vector pointing from the negative to the positive side. This does not hold in the presence of the strong discontinuities since the direction of the normal vectors may not be necessarily the same. Weak discontinuities denote the discretization error in the DG method and not any physical cracks. This is due to the fact that continuity constraint in most DG methods is fulfilled only in a weak sense. On the contrary, strong discontinuities represent a physical interface (see for instance Reference 34). In Figure 2, the tractions are drawn only at a point (later known as integration point) for the sake of simple illustration. Obviously, tractions exist on the entire interface. The equilibrium of the inner and outer forces obtained from the quasi-static linear momentum balance along with its boundary conditions reads as follows: In the given strong form, represents the Cauchy stress tensor and f the body force vector. Note that in case of weak discontinuities (prefailure), the continuity condition of the displacements and the tractions on the inner borders of the subdomains is satisfied in the continuous solution A detailed definition of the boundary conditions is described in Reference 54.
It is necessary to differentiate between the negative and positive quantities on the sides of the discontinuities Γ. To this end, the average and jump operators are defined below In contrast, in case of strong discontinuities, the displacement jumps are allowed. Nevertheless, the continuity of the tractions remains valid. Note that these assumptions hold for "cohesive interface models." For more information, refer to Reference 56. At the interface, the tractions are a function of the displacement jump: More details regarding the interface model will be given later. This work is restricted to linear elasticity. Hooke's law is given by where C represents the fourth-order elasticity tensor and is the infinitesimal strain tensor defined by where x is the material position vector.

Weak form
By the use of the Gauss theorem and considering the symmetry property of the Cauchy stress, we obtain the weak form for the standard continuous Galerkin method where no internal boundaries exist. In the present context, the surface integral must not only be applied to the outer boundaries Ω with prescribed tractions t p , but also to the inner boundaries Γ with t + and t − . Since the inner boundaries differ based on the existence of either the weak or the strong discontinuity, the formulation of the weak form varies accordingly.

Weak discontinuities
Taking the weak discontinuities first into consideration, the final weak form for the DG method is given by Recalling Equations (1) and ( 5), the weak form is further simplified to Applying the relation 34 and by considering the traction continuity condition (see Equation (4)), the final weak form of the DG is given by In order to obtain a stable solution, a penalty term is added in the spirit of Nitsche's method. 26 Please note that the addition of this term does not violate the consistency of the weak form with the strong form since the jump of the displacement field in the exact solution vanishes. Consequently, one obtains where = E∕h[N∕m 3 ] is a penalty parameter and depends on Young's modulus E, the mesh size h and a sufficiently large positive value . 50,57

Strong discontinuities
Once failure occurs, strong discontinuities take place and the DG method cannot be applied any more as displacement jumps occur in the continuous solution. Instead, a CZ model is applied to prescribe the tractions t as a function of displacement jumps [[u]] on the discontinuities Γ.
Pursuing Equation (9) in view of its corresponding boundary conditions (Equation (6)), the weak form of the CZ becomes In analogy to the derivation of the weak form in DG (Equation (11)), we obtain the final weak form of the CZ as follows:

Unified weak form
Here, in order to incorporate the strong discontinuities into the DG framework, the parameter is utilized to specify which weak form is active. Inserting the CZ term into the DG weak form, one obtains where = 0 denotes the prefailure as default. Once the failure criterion is met, = 1 is set and remains as a history variable at the integration point of the interface element. The failure criterion as well as the choice of are clarified in the next section.

CZ model
The interface behavior, introduced in the CZ model, is described by means of a traction-separation relation. In this work, isothermal and isotropic interface behavior is assumed. Furthermore, it is well-known to introduce a local coordinate system at the interface to differentiate between the normal and shear directions. To this end, displacement jumps [[u]], defined on global coordinates, are substituted by the gap vector which are computed of the crack tip. Here, g s and g n are the gap components in the shear and normal directions, respectively. Clearly, they are orthogonal to each other. Unlike the weak discontinuity in linear elasticity, the opposing sides of the CZ may not stay parallel to each other during separation. Therefore, a mid-plane is defined in the case of a four-node element, with respect to which the shear and normal directions are defined ( Figure 3B). Denoting the angle between the local and global coordinates as and defining the rotation matrix we can compute the gap vector based on the displacement jumps as follows: In order to differentiate between an opening and a closure of the crack especially in the case of contact, we define the effective quantities by the use of the Macaulay brackets Accordingly, the effective gap vector, as well as the effective traction vector, read as follows: .
Here, assigns different weights to the normal and tangent separations. The effective separation and the effective traction values are given by respectively. Finally, we define the TSL where max represents the maximum separation that the interface has reached so far and f is the elongation at full failure. In addition, the maximum strength of the interface is denoted by t 0 . The power m is a material parameter that controls the convexity of the drop of the cohesive tractions whereas the power n is mainly introduced to obtain numerical stability in case of contact. The penalty term includes the contact force (see, eg, Reference 58) to avoid penetration while the last term represents a viscous force with the viscosity parameter . The latter has been introduced to consider viscous effects at the interface as well as to avoid possible numerical instabilities such as snap-back (refer to References 15,59). For the contact force, the same penalty values as for the DG part is used. Note that max is an internal variable at each integration point and can be related to a scalar damage parameter D. In fact, the relation holds with D = max ∕ f . Nonetheless, we derive our equations applying the notation max in the following. On the local coordinates, once the effective tractions t eff resulting from the DG part t DG = { } n on Γ reach the maximum strength t 0 , the cohesive traction comes into play and depending on the separation and its history, different forms of cohesive tractions are applied. In other words, in case t eff DG > t 0 , the traction vector will take one of the following forms: Bear in mind that the first case shows the loading case whereas the second one could be either unloading or reloading up to max . In contrast, the last case serves as the ultimate failure with no cohesive tractions. In all three cases, as soon as the normal separation begins to obtain negative values, it gets penalized via the contact contribution ⟨−g n ⟩. This is done in a similar manner to the DG penalty term with the same penalty parameter. Nevertheless, the shear contribution remains unpenalized.
Two different two-dimensional (2D) fracture modes, namely I and II, are illustrated in Figure 4 for the extrinsic CDG TSL. Note that in the following figures, the initial slope is plotted in an over-magnified manner to demonstrate the presence of the DG in the pre-failure regime.

Discretization
In this section, we only consider the terms on the discontinuity (interface) Γ and the rest of the terms in the bulk and on the Neumann boundaries are treated in a standard approach. Furthermore, the quantities are transferred into Voigt notation due to the symmetric structure of the Cauchy stresses. The 2D discretized form of the unified CDG weak form is derived for two different linear elements, namely the four-node quadrilateral and the three-node triangular elements. The former is applied on matching meshes while the latter is exploited in the presence of hanging nodes. On each discontinuity (interface) element Γ e , the element position x h , displacement u h , test function u h , as well as the stresses h , are interpolated from their nodal values (denoted in capital letters) by the use of linear shape functions N. These quantities read as where n en is the number of nodes per element and is the reference element coordinate. Note that the nodal stresses on Γ are directly imported from the neighboring bulk elements. Henceforth, we neglect the superscript h in the approximated quantities. Depending on the type of element, n en can vary between 3 and 4. As shown in Figure 2A, the normal vector to the weak discontinuity Γ (prefailure regime) is given in a matrix form as follows: .
In the following, the quantities which differ based on the type of discretization, are presented separately.

Four-node quadrilateral element (matching meshes)
Recalling n en = 4, the discretized forms of the displacement jumps and stress averages on Γ read as where the indices 1 and 2 refer to the negative side of the discontinuity/interface while 3 and 4 refer to its positive side. In addition, N J and N A are the matrices of all shape functions for the jump and average operators, correspondingly. For a detailed definition of the above-mentioned matrices, please refer to Appendix A1.

Three-node triangular element (nonmatching meshes)
For this element with n en = 3, the discretized forms of the jump and average terms read as As it is seen, the interpolation takes place only between the nodes 1 and 2. More details will be given in Section 2.6.

Linearization
Once the spatial quantities are discretized, the final discretized weak form regarding only the terms on Γ on the element level can be written as Here, remember that the gap and its variational form are already expressed in terms of the jump as follows: Due to the arbitrariness of U, the residual forces for the terms on Γ are given as follows: Prior to the failure t eff DG ≤ t 0 , that is, = 0, the residual forces consist of only the DG terms. The linearization of them leads to the DG stiffness matrix where the linearization of the first DG internal residual force is neglected, based on the work of Reference 22. As a result, the cumbersome implementation of the DG contribution is reduced. Moreover, the neighboring bulk elements do not need to be evaluated in the discontinuous element. Nonetheless, the contribution of the bulk elements have to be computed prior to that of the DG elements. For a detailed assembly process, refer to Reference 22. As soon as the failure initiates t eff DG > t 0 , that is, = 1, the residual forces and the stiffness matrices will solely result from the cohesive tractions of Equation (27) where K is the 2 × 2 cohesive stiffness matrix. As a remark,̇g in the discrete form is given by Δg∕Δt with Δt being the size of the time step for the corresponding increment Δg = g t n − g t n−1 . Accordingly, after a lengthy but straightforward derivation, we obtain the stiffness matrices K for different cases based on the state of the effective separation consistent with Equation (27). We have to differentiate between the two regimes. The first regime ( max < f ) represents the TSL prior to the full failure. The second regime ( max ≥ f ) represents the full failure, where no tractions exist anymore.
• Regime one (before full failure) max < f : Within this regime, again two different cases can show up depending on either the system is been loaded (case a) or unloaded/reloaded (case b).
-Case a: during loading = max : where I is the identity matrix of second-order and sgn(x) is the sign function, which extracts the sign of the real number x. -Case b: during un-and reloading < max : • Regime two (after the full failure) max ≥ f : .

Numerical integration
The integrals of the weak form (Equation (17)) are evaluated numerically by application of the quadrature points. Taking only the left-hand side of this equation into account, one denotes the bulk and the discontinuity Γ terms as surface and line integrals in a 2D configuration, respectively. The first term is fully integrated on the surface using four Gaussian F I G U R E 5 Different integration schemes on the Q1 and cohesive discontinuous Galerkin element quadrature points with linear shape functions in a quadrilateral element (Q1). The rest of the terms are integrated on Γ in the form of four-and three-node elements explained in the following.

Four-node quadrilateral element (matching meshes)
For matching meshes, we consider a four-node element with two line-integrals for each side of the discontinuity. They are integrated using three different integration schemes, namely nodal, full and reduced Gaussian (see Figure 5). The idea of the application of reduced integration was investigated in the work of Reference 60. In case of plasticity, spurious penetration modes were observed for the interface element between two quadratic tetrahedra when they were under-integrated. In the present work, this problem did not occur in the context of linear elasticity. However, the solution to this problem when it comes to plastic deformations could be achieved by the use of a mixed (selective) integration of the interface/boundary terms (see Bayat et al 50 ). In the aforementioned paper, the penalty term was fully integrated while the other boundary term ("DG term") was under-integrated. This possible remedy for the aforementioned unstable modes in plasticity remains to be shown and is not investigated in the context of this paper.

2.6.2
Three-node triangular element (nonmatching meshes) In contrast, there is only one integration point located at p in the local coordinates in the presence of hanging nodes (see Figure 6). The corresponding tractions and displacements are interpolated on p using the information of the neighboring nodes 1 and 2 where N 1 = 1∕2 (1 − ) and N 2 = 1∕2 (1 + ) are the linear shape functions on Γ − . Thereafter, the displacement jump [[u]] (or the gap g in CZ) and average tractions {t} (or the prescribed tractions t cz in CZ) between node 3 and X 3 ′ are evaluated.
In the reference configuration, one obtains the position of p as where x c = (x 1 + x 2 )∕2 is the x coordinate in the middle of the nodes 1 and 2 and Γ represents the length of the side of the triangle, lying on the discontinuity. For the three-node CDG element, the definite integrals on Γ of the residual and stiffness matrices are computed numerically as follows: where l el is the area of the contribution of node 3 on the finer mesh (see Reference 47). This area can be calculated depending on whether the discretized bodies are meshed regularly as pictured in Figure 7. For regular meshes, l el = l∕n is used for the inner nodes whereas l el = l∕2n is used for the outer nodes. Here, n refers to the number of elements on the finer mesh. In the case of irregular mesh, l el = (b + c)∕2 holds for the inner nodes whereas l el = a∕2 is applied for the outer nodes. Nonetheless, once the DG elements are inserted within the lower and upper elements as well, the area l el will be simply half of the length of each element, that is, for the regular mesh l el = l∕2n whereas in case of irregular mesh, for example, here l el = a∕2, l el = b∕2, and l el = c∕2 (see Figure 8).

IMPLEMENTATION
In the assembly, the bulk elements are prior to the CDG elements since the stresses from the bulk should be transferred to the elements on Γ. To this end, the stresses ( xx , yy , and xy ) at the Gauss points of the bulk elements are first extrapolated to the nodes and stored in a global matrix. During the assembly of the interface elements, the stresses from the adjacent nodes are called and interpolated at the integration points of the CDG. This is shown in Figure 9. The implementation of the CDG element is clarified in Algorithm 1. Due to the fact that the TSL is evaluated at each integration point, also the failure is evaluated on the level of integration point. This means that within an element, one integration point may not have experienced failure while the other one has already failed, see Equation (17). Note that we neglect the separation 0 at the maximum strength t 0 due to its insignificant contribution. These gaps result from the permissible displacement jumps in DG prior to the failure. Furthermore, max as well as (both initially set to zero) are variables to be stored at the integration point in order to differentiate between loading and un-or reloading and the state of failure, respectively. F I G U R E 9 Extrapolation of the stresses from bulk Gauss points to the mutual nodes as well as interpolation of them on Γ The mesh for both the three-and the four-node elements (CDG) as well as the bulk elements (Q1) are generated using MATLAB. In the case of matching meshes, first the bulk elements are generated and then the four-node CDG elements are placed between all adjacent volume elements in a straight forward manner.  (39)) end else error calculating t eff DG end end R Γ = R DG + R CZ and K Γ = K DG + K CZ end As clarified before, the placement of the three-node elements depends on whether the bulk is discontinuous or conforming. In the case of a continuous bulk, each node of the triangular elements on the finer side of the mesh is checked to find its two directly neighboring nodes on the coarser side. 47 Thereafter, the triangular element is constructed using these three nodes. If a node on the finer mesh coincides with that of the coarser mesh, it will be connected to both the left and right nodes of the coarser mesh to make two triangles. This is illustrated in Figure 7.
In contrast, once the bulk is nonconforming, there will be double nodes at the same point (in the reference configuration). The mutual overlying nodes of the finer mesh will be connected to the closest nodes of the coarser mesh making two overlapping triangular elements. Those two nodes of the finer mesh which coincide with the other two nodes of the coarser mesh make two nonoverlapping triangles on each side for the coarser mesh as illustrated in Figure 8.

NUMERICAL EXAMPLES
Unless otherwise stated, a geometrically linear formulation as well as the extrinsic CZ approach without viscosity ( = 0) are used. The Finite Element Analysis Program 61 is utilized for the computations whereas the mesh is generated separately using a user-defined code in MATLAB. First, a 2D example is given to verify different aspects of the new formulation. Afterward, five benchmark problems are considered, three of which are chosen to investigate the locking behavior whereas the rest exemplifies the application of nonmatching meshes. The above-mentioned examples are as follows: • peeling test: mode I delamination (matching meshes), • peeling test: mixed-mode delamination (nonmatching meshes), • mode II delamination: end notched flexure test (matching meshes), • single-edge notched specimen (matching meshes), • fiber pull-out in composite structures (nonmatching meshes).

Two-dimensional validation test
In order to validate the performance of the new formulation, we test it for different loading cases in two simple examples, namely, modes I and II separations. The tractions, as well as the reaction forces for unloading and reloading including contact as well as the damage power (m), are illustrated. In the end, a comparison of the three-node vs four-node element formulations is carried out.

Normal separation
A rectangular block (1 × 2 mm) is fixed on its lower edge in y-direction while being pulled on its upper edge in y-direction as pictured in Figure 10A. To avoid rigid body motions, both left corners are fixed in x-direction additionally. The corresponding parameters are shown in Table 1. The problem is displacement driven (u = 5 mm) and the response of the reaction forces as well as the traction separation diagram (for an integration point) are plotted in Figure 11, respectively. The elements at the interface plotted in Figure 10C and d do not represent any material at the interface.

Shear separation
The same block from the above example is considered here with different boundary conditions. At the bottom, it is completely fixed in both directions while at the top only the y-direction is constrained. A horizontal displacement is applied to the upper nodes (see Figure 10B) and the response of the reactions and tractions for different damage powers (m) is illustrated in Figure 12. Note that m is a material parameter depending on how the material behaves during failure. The parameters are found in Table 1.

Three-node triangular element (nonmatching meshes)
The same block with the very same parameters is investigated using the three-node elements (see Figure 10D). The traction separation diagrams for both modes I and II are plotted in Figure 13. It is noticeable that the performance of the three-node elements is the same as that of the four-node elements.

Peeling test
A benchmark problem of the peeling test is considered in this part for matching and nonmatching meshes. The former is used for the mode I fracture while the latter is applied to capture mixed-mode behavior. Firstly, for a single beam, we demonstrate a mesh convergence study with different element formulations applying only the DG part of the formulations. Afterward, the cohesive behavior of the peel test is studied in terms of mesh convergence and numerical integration schemes.

Mode I delamination (matching meshes)
A structure consisting of two beams with adhesive in between is fixed at its left end while being pulled apart in y-direction at its right end. A precrack on the right-hand side of the beam is considered for the initiation of the crack. The geometry, as well as the boundary conditions, are given in Figure 14. The discretization of each beam is carried out in the same manner simultaneously with the following meshes 2 × 10, 4 × 20, 8 × 40, etc. Material parameters exemplify a harder bulk, here steel, with a weaker adhesive in between as in composites to model delamination between the laminae. The material parameters are given in Table 2. In this example, the mode I fracture is illustrated to study the behavior of the cohesive model. A constant element aspect ratio of 10 ∶ 1 is set. This is of interest to show that shear locking can be overcome. 50 Artificial stiffening or the so-called locking phenomenon can lead to incorrect reaction forces and thus an unrealistic behavior of not only the prefailure regime 49,50 but also the failure regime. The former is investigated by considering only the upper beam whereas for the latter, the entire geometry is taken into account. Different discretization levels as well as different element formulations are compared. In order to apply the DG method in the prefailure regime, the correct parameter, introduced in Section 14, is specified for each numerical integration scheme. A single force 10 −4 N acts on the upper right corner of the upper single beam. The corresponding vertical displacement of the same corner is investigated for different levels of mesh refinements. The standard four-node linear rectangular Q1 is used for the bulk in combination with our 4-node DG elements utilized with full and reduced Gaussian integration schemes as defined in Figure 5. Note that the DG elements are spread out between all adjacent Q1 elements. Moreover, a Q1 element with the element aspect ratio of 1 ∶ 1 is considered as a reference solution. The relative vertical displacement (u y,relative ) of the upper right corner of the beam with respect to the number of elements in y-direction and additionally once with respect to the total number of degrees of freedom is given in Figure 15.
It is clearly seen that the DG method with reduced integration overcomes locking whereas the full integration does not eliminate the influence of the artificial stiffening. In addition, the influence of the stabilization (penalty) parameter on the convergence of the DG elements is studied in Figure 16. The relative vertical displacement (u y,relative ) of the upper right corner of the beam is plotted with respect to the penalty parameter . Two different mesh refinement levels for both cases of reduced as well as full integration schemes are considered. By varying (see Equation (14)) from very small to very large values, the response of the elements varies from too soft to converged and later to unrealistic slightly stiffer values. Choosing too high values for results in numerical instabilities.
Within the failure regime, the application of the DG elements in the bulk and CDG elements at the interface with different integration schemes similar to the previous example is studied. To better illustrate the differences in the performance of the element formulations in the presence of shear locking, an additional intrinsic CZ model from Reference 15 is employed in combination with Q1 elements. In this variant, no DG elements are introduced in the bulk and the CZ elements are evaluated by a Gaussian full integration scheme. Figure 17 shows the convergence rate of the aforementioned elements in terms of the reaction force-displacement curve. Note that the number of elements in the graph refers to the upper beam of the peel test for a better comparison with the intact beam example. Obviously, the complete peel needs twice the number of elements.
As clearly seen in Figure 17, the application of the reduced integration scheme in the CDG formulations results in realistic (here, lower) reaction forces. On the other hand, the false (here higher) reaction forces will bring about a premature tearing of the interface.
In the case of crack propagation of the brittle materials, depending on the material as well as the interface properties such as fracture energy, there need to be enough integration points on the crack path, leading to a smaller distance between the adjacent integration points. This is better seen in Figure 18 for the mesh of 2 × 10 elements. However, by refining the mesh only one level higher (4 × 20 elements), a smooth force-displacement curve is obtained. The same converged solution in the case of full integration is obtained at the earliest with 128 × 640 elements (see Figure 17A). Therefore,

Mixed-mode delamination (nonmatching meshes)
Peeling can very often occur when one layer is considerably thicker/stiffer than the other layer. As a result, discretization of the two layers with the same mesh refinement level is computationally unnecessarily costly. Here, we consider a peeling example in which the stiffer layer possesses a higher thickness than the softer material. This implies that the stiffer bulk does not necessitate a fine mesh. Thus, the three-node CDG elements are applied. The geometry along with the boundary conditions is depicted in Figure 19. The stress states yy for the matching as well as the nonmatching meshes for different time steps are plotted in Figure 20 against each other for a better visual comparison. The corresponding parameters are shown in Table 3 for the lower, upper, and middle layers. The introduction of the DG elements within the bulk is possible as explained in Figures 7 and 8.
It must be noted that unlike other matching discretization techniques, the level of the refinement in this type of the nonmatching discretization 47 must be adequate on both sides of the interface. Otherwise, an unsymmetric stress pattern is observed at the interface. In Figure 20, this discretization error is because the coarse mesh is not fine enough.
The reaction force displacement of the upper-right edge of the geometry is plotted in Figure 21 for both meshes once with and once without the DG elements in the bulk. The converged response is obtained with much less number of elements in the nonmatching meshes.

Mode II delamination: End notched flexure test (matching meshes)
In order to investigate the response of the model for the mode II delamination, a beam (see Figure 22) with a notch on its right end similar to the DIN ENF 6034 standard is considered. Geometry, boundary conditions, and the given displacement  Figure 22A. In addition, the discretization scheme and the distribution of the stress ( 22 ) are depicted in Figure 22B. The element ratio is fixed to 1 ∶ 10, that is, the geometry is discretized by 4 × 10, 8 × 20, 16 × 40,…, 256 × 640 elements to assure the occurrence of shear locking. The parameters for the extrinsic CDG as well as the intrinsic CZ 15 elements are listed in Table 4. It is noticeable that due to the existence of the contact algorithm, the intrinsic CZ model needs a penalty value to be defined. In this way, the extrinsic CDG benefits from a fewer number of parameters.
The reaction force-displacement response of the specimen on its upper mid-point is shown in Figure 23. The aforementioned mesh refinement levels are compared for the upper beam of the specimen in the vertical direction.
The CDG elements with reduced integration scheme result in an astonishingly fast convergence rate in terms of the number of elements in the vertical direction (only eight elements). At this level of mesh refinement, other element formulations show over 30% error.

Single-edge notched specimen (matching meshes)
A 2D plane strain specimen with a zero-thickness notch is considered here. It is pulled apart in y-direction on its upper edge while being fixed at the bottom. The geometry, boundary conditions, loading as well as the discretization are given in Figure 24. Additionally, a cohesive layer is laid along the notch direction for the evolution of the crack. Different discretization levels are exploited to investigate the convergence behavior of the specimen in the presence of near incomprehensibility. The corresponding parameters are given in Table 5.
In order to investigate the influence of volumetric locking on the convergence rate of the block without any notches, an intact specimen is considered first (see Figure 25). Two different discretization schemes, namely continuous with  only Q1 elements and discontinuous with DG elements among bulk elements (Q1) are utilized. The latter is exploited in two variants, that is, reduced and full integration on Γ.
As clearly seen in Figure 25, the converged reaction forces are obtained with a considerably less number of elements/degrees of freedom once employing DG elements with reduced integration.
Next, the reaction forces are plotted against the prescribed displacements once the notch and the cohesive layer exist (see Figure 26). For the same mesh, the maximum reaction forces are reached at much lower displacements (around 20%) in the CDG element with reduced integration scheme compared to other element formulations. This is due to the elimination of the volumetric locking.
In this example, the reaction force-displacement response of the extrinsic CDG with nodal integration scheme for the mesh of 32 elements is shown in Figure 26D the application of the nodal integration does not improve the convergence rate of the elements with respect to the mesh refinement.

Fiber pull-out in composite structures (nonmatching meshes)
The final example demonstrates the benefit of the application of the nonmatching meshes in composite structures. A notorious instance of failure is observed once the fibers undergo deformation and are finally pulled out. Depending on TA B L E 6 Parameters of the fiber pull-out test the interface properties, they may fail after the matrix does. In this example, due to a precrack in the matrix, once the fibers are pulled, the notch in the matrix propagates and reaches the interface. Thereafter, due to the fixation of the matrix on the left-hand side, the failure starts at the interface. It is noticeable that due to the material properties of the bulk (see Table 6), a fine mesh for the fiber is redundant. Geometry, boundary conditions, and loading are given in Figure 27A. Furthermore, the discretization for both matching and nonmatching meshes are provided in Figure 27B.
For the same time, the stress contour of the matching and nonmatching meshes look alike (see Figure 28). Thus, the application of the nonmatching CZ interface, denoted by red line in Figure 27A on the mid-layer saves a considerable amount of computation time.
In order to shed more light on the differences between these two discretization schemes, the reaction force-displacement response of the structure is plotted in Figure 29. To this end, the reaction forces on the right-hand side of the fiber are plotted against the displacement on the lower-left corner of the specimen. As illustrated, the nonmatching mesh delivers the same behavior as that of the meshing meshes. In addition, by the use of reduced integration, a smoother transition of the pre-to postfailure regimes is observed.

CONCLUSION
A novel extrinsic cohesive formulation was embedded in the framework of the IIPG variant of the DG method. The method was applied in the framework of 2D linear elasticity to model material brittle failure. In addition, two different discretization techniques, namely a four-node quadrilateral element formulation and a three-node rectangular element formulation were exploited to account for matching as well as nonmatching meshes, respectively. In order to circumvent the artificial stiffening within the bulk and to obtain realistic interface behavior, a reduced integration scheme was employed on the boundary terms. A series of validation tests were run to evaluate the performance of the CDG method. Furthermore, mode I, mode II and mixed-mode delamination tests along with a single-edge notched specimen and a fiber pull-out test were investigated to convey a study on the influence of locking on the interface behavior as well as the application of nonmatching meshes. It was found that not only shear locking but also volumetric locking effects were minimized. Thanks to the inherent nonconforming characteristic of DG methods, hanging nodes could be tackled on the discontinuities/interfaces. as a result, the calculation cost was reduced. In addition, the choice of such an extrinsic formulation of the CZ model allowed to avoid unrealistic compliance at the interface. Unlike the conventional extrinsic methods, our CDG method does not require a remeshing of the structure at the interface by virtue of the existence of the DG framework. Determination of the crack path only by means of CZ methods remains debatable. CZ models are able to qualitatively specify the crack path (see Reference 62). This is due to the fact that these methods impose the requirement that the crack path is biased by the direction of the mesh lines. 63 Therefore, a combination of the CZ models with other methods such as nonlocal damage models can capture arbitrary crack paths (eg, refer to References 43,64). However, performing such simulations lies beyond the scope of the current work and will be considered in a forthcoming paper. In addition, it is planned to extend the model for ductile fracture especially in a three-dimensional setting. Furthermore, the use of a locking-free formulation for the bulk is expected to improve the performance of the method in problems prone to locking. The jump and the average shape function operators for the four-node elements are given as: whereas the aforementioned operators for the three-node element are given by (A6)