Modeling of failure at the interface of ductile materials by applying the cohesive discontinuous Galerkin method

In this study, the failure behavior at the interface of ductile materials is investigated. In order to capture the degradation of the tractions at the interface, a cohesive zone (CZ) model is applied. The choice of the type of the CZ approach, i.e. either intrinsic or extrinsic, brings about different drawbacks. The former includes an elastic regime at the interface prior to the failure, which can result in numerical difficulties whereas the latter necessitates the re-meshing of the structure during crack propagation. In order to overcome these problems, the incomplete interior penalty Galerkin variant of the discontinuous Galerkin (DG) method is applied both at the interface and in the bulk instead of the standard conforming finite element method. In addition, the application of the DG method enables to use nonmatching meshes in the discretized model. To treat the bulk, an elastoplastic material model with isotropic hardening as well as different hardening rules for small strains is incorporated into the DG framework. Two numerical examples are computed to study the convergence behavior of the new cohesive discontinuous Galerkin (CDG) method in comparison to that of the conventional models. The new CDG method outperforms the conventional CZ continuous Galerkin elements in the presence of locking effects as well as hanging nodes.


Intr Introduction oduction
In the world of engineering, the prediction of failure plays an important role when it comes to material modeling.
Several methods [1] for the modeling of failure have been developed in the last decades depending on the material, geometry, scale, etc. For instance, to model the failure within the bulk, sophisticated advanced damage models based on continuum damage mechanics [2,3,4] are of significance while cohesive zone models [5] mainly address failure at the interfaces.
On the one hand, cohesive zone models are a simple tool to model failure at the interfaces. On the other hand, the simulation results based on CZ models usually agree well with experimental findings [6]. The embedment of CZ models into standard conforming finite elements usually results in some restrictions at the interface. These can be the occurrence of the initial elastic regime prior to failure as in intrinsic CZ models or the necessity of re-meshing of the structure during the crack propagation as in extrinsic approaches [7]. In addition, extra effort has to be made to include hanging nodes in nonmatching meshes [8,9]. One possible solution is to utilize the CZ model together with the DG method at the interface as well as in the bulk [10,11].
Discontinuous Galerkin methods are an alternative to standard continuous Galerkin (CG) methods, in which the continuity of the displacements and tractions are usually only weakly satisfied. Due to various formulations of the DG methods [12], they have found different applications from fluid [13] to solid mechanics [14]. Recently has been a growing trend in the application of DG methods in combination with cohesive zone models in fracture mechanics [15,16,17]. One significant motivation of this popularity is the fact that adjacent elements in DG methods do not share mutual nodes in the discretized structure. Consequently, cracks (also known as strong discontinuities) can propagate through the existing inherent weak discontinuities in the DG elements either within the bulk [18] or at the interface [19]. Another advantage of the use of the discontinuous Galerkin method lies in the removal of locking effects [20,21,22,23]. In case of volumetric as well as shear locking problems, DG methods with reduced integration perform remarkably faster than the standard continuous Galerkin element formulations [10,24,25] in terms of mesh convergence rate. A locking-free modeling of the fracture at the interface of brittle materials with non-matching meshes is established in [26].
In the present work, (quasi-) ductile materials [27,28] are taken into account whereas our previous work aimed at the modeling of failure at the interface of brittle materials [26]. To this end, a plasticity model with isotropic hardening is considered for small strain problems using the incomplete interior penalty Galerkin (IIPG) variant of the DG methods [10]. Hanging nodes are incorporated into the discretized model as well. The performance of the cohesive discontinuous Galerkin method in comparison to the standard conforming CZ method is investigated in the problems prompt to volumetric locking as well as including different materials.
The structure of this paper is as follows: First, the strong form of the linear momentum balance is introduced followed by its weak form. Next, the traction separation law (TSL) is given. Afterwards, the elastoplastic model based on the Helmholtz free energy is defined. The corresponding state relations and evolution equations are derived obeying the second principle of thermodynamics. The discretization as well as the numerical integration are given briefly in the framework of the discontinuous Galerkin method. The next chapter is devoted to numerical examples to study the convergence behavior of the CDG method in comparison to that of the conventional models. Finally, conclusions are drawn and possible improvements are suggested.

W Weak f eak form orm
The unified weak form of equation (1) considering both weak and strong discontinuities is defined by where = /ℎ is the penalty parameter. The positive free parameter is a sufficiently large value, is Young's modulus and ℎ is the mesh size. The parameter is used to switch between intact ( =0) and damaged interfaces ( =1). In addition, is the gap vector with components in shear ( ) and normal direction ( ) with respect to the discontinuity. In fact, the gap is the displacement jump in a local coordinate system ( , ) and is given by = ⟦ ⟧. Here, is a rotation matrix. For more details, refer to [26].

T Tr raction separ action separation la ation law (TSL) w (TSL)
Finally, the tractions on the strong discontinuity are calculated by In the above equation, 0 is the maximum strength of the interface right before failure, =√(〈 〉 2 + 2 2 ) is the effective separation, is the maximum separation reached so far and is the separation at full failure (see Fig.   2). To control the contribution of the shear separation, the parameter is defined. Here, is a material parameter to control the convexity of the fading tractions in the traction separation law while is mainly introduced to avoid numerical instabilities during contact. The parameter is the same as that of the DG weak form in equation 4. The kinematics of the problem is specified for elastoplasticity with small strains. As a result, an additive decomposition of the strains follows: The elastic part of the strain is denoted by =sym (grad ( )) whereas is the plastic strain. In order to calculate the stresses and the plastic strain, we start with the Helmholtz free energy It is additively divided into an elastic and plastic part. The elastic part is given by Finally, associative plasticity is assumed, leading to the following evolution equations: Here,̇is the plastic multiplier. Practically, two cases can occur: either the material is elastic (Φ<0 and=0) or it ), we take the standard four-node quadrilateral elements with polynomials of first order, namely Q1 elements. The discontinuity Γ is discretized by four-node DG elements between the bulk elements or by CDG elements on the interface. In addition, three-node (C)DG elements are placed where hanging nodes appear in the structure. A detailed description of the aforementioned procedures is out of the scope of this work and can be found in [26].

The discretization of the evolution equations (Equation 12) in time is established by means of a fully implicit backward
Euler scheme. Within the time interval [ , +1] in a Newton-Raphson iteration at the global level, equation (12) is solved iteratively. In the local plasticity iterations, for a given , the quantities , and Δ have to be determined. Here, Δ =Δ ̇ is the plastic multiplier. After solution of the local iterations, the consistent tangent modulus is derived. Fig. 3. Int Fig. 3. Integr egration schemes f ation schemes for bulk (Q1) and discontinuous ((C)DG) elements. or bulk (Q1) and discontinuous ((C)DG) elements.
The integration of the integrand is carried out numerically by applying the Gaussian quadrature rule. To this end, the bulk element is integrated in a full manner, i.e. using four Gauss points in each element. The terms on the discontinuity are integrated in different schemes. The standard procedure is the full integration besides the nodal integration.
However, to eliminate locking effects, reduced integration is applied as in [26]. This is pictured in Fig. 3.

Numerical e Numerical ex xamples amples
This section includes two different examples of elastoplastic structures. The first one is a single-edge notched specimen from [26] with matching meshes. This example aims at the study of the volumetric locking effects and their influence on the crack propagation behavior. The second example is a fiber composite structure from [29]. The debonding between the fiber and the matrix is modeled enabling hanging nodes in the discretized structure.
1.3.1 Sing 1.3.1 Single-edge not le-edge notched specimen ched specimen In this example, we examine a plane strain sample with an inherent notch. It is fixed at the bottom and a displacement Modeling of failure at the interface of ductile materials by applying the cohesive disc...

1856/6
load is prescribed at its upper edge. Geometry, boundary conditions, loading, discretization as well as the evolution of the plastic zones are illustrated in Fig. 4. Material parameters are given in Table. 1. There is a cohesive layer along the notch, which is discretized by CDG elements. Both the notch and the cohesive layer are of zero thickness. The bulk is discretized by square elements with the same number of elements in each direction for each mesh refinement. The aim is to investigate the convergence behavior of different element formulations for near incompressible elastoplastic materials in the presence of a crack.  The reaction force-displacement curves for the given example are plotted in Fig. 5. Due to the incompressibility of the material both in elastic and in plastic regimes, volumetric locking effects are observed (see Fig. 5a). Nonetheless, by the use of the reduced integration scheme on the boundary terms, this problem is remarkably diminished. Figure 5b shows that converged results are achieved with almost eight elements in each direction. A better comparison of the convergence rate for the same mesh refinement level (16 elements) is shown in Fig. 5c, where full integration shows overestimated reaction forces and even with 128 elements in each direction, it does not reach the converged solution completely. Finally, in order to illustrate the presence of the plastic strains, the specimen is unloaded at 40% of the total given displacement to its origin and then reloaded thoroughly (refer to Fig. 5d).   [29] including a fiber with a circular cross-section surrounded in matrix is studied here. Geometry, boundary conditions and loading of the composite structure are pictured in Fig. 6a for a quarter of the geometry due to its double symmetry conditions. The distribution of the accumulated plastic strains in its discretized form together with the propagation of the crack at the interface of the fiber and matrix are depicted in Fig. 6b. In addition, material parameters for the matrix, fiber and the cohesive layer between these two materials are given in Table. 2. The structure is being pulled apart at its left and right edges through a prescribed displacement.
The convergence behavior of the CDG elements in terms of reaction force-displacement curves with and without nonmatching meshes is compared to that of a standard continuous intrinsic cohesive zone element formulation [30] with matching meshes (see Fig. 7). The DG elements in the fiber structure are integrated in a reduced manner while the CDG elements can be integrated using different schemes (here: full, reduced and nodal).
Modeling of failure at the interface of ductile materials by applying the cohesive disc...   Figure 7a shows that by application of the non-matching meshes, a smaller number of elements is required to capture the same behavior. This is due to the fact that the fiber is 200 times stiffer than the matrix. Thus, it does not require a fine mesh as that of the matrix. The numbers indicate the sum of the elements in each direction for the matrix and the fiber, respectively. Figure 7b illustrates the influence of the hardening isotropic plasticity on the behavior of the bulk as well as the interface for the same given parameters. Having a pure elastic material results in a completely different interface behavior. r refinement le efinement lev vel of 32+32 elements. el of 32+32 elements.

Conclusion Conclusion
A recently developed novel extrinsic cohesive discontinuous Galerkin method was extended to model the failure at the interface of ductile materials. To this end, a cohesive zone model was implemented into the incomplete interior penalty Galerkin method. The discretization of the structure enables both matching as well as nonmatching meshes.
The material model is elastoplastic for small deformations. Isotropic hardening with linear and Voce type variants were applied to investigate two different examples. Using nonmatching meshes were advantageous in the discretization of a fiber composite structure. This led to a reduction of the computation costs. Furthermore, the method showed to be free of volumetric locking due to the reduced integration of the boundary terms in the DG method. To avoid shear locking in elastoplastic materials, a selective integration scheme on the discontinuities could be of significance. This is currently under investigation.

A Ackno cknow wledgements ledgements
Financial support of this work, related to the subproject A6, "Multiscale modelling of the damage and fracture behavior in nanostructured laminates", of the Transregional Collaborative Research Center SFB/TRR 87 "Pulsed high power plasmas for the synthesis of nanostructured functional layers" as well as the project "Hybrid discretizations in solid