A stabilized finite element method for delamination analysis of composites using cohesive elements
AA stabilized finite element method for delamination analysis ofcomposites using cohesive elements
Gourab Ghosh a , Ravindra Duddu a, ∗ , Chandrasekhar Annavarapu b a Department of Civil and Environmental Engineering, Vanderbilt University, Nashville, Tennessee. b Department of Civil Engineering, Indian Institute of Technology Madras, Chennai, India.
Abstract
We demonstrate the ability of a stabilized finite element method, inspired by the weighted Nitscheapproach, to alleviate spurious traction oscillations at interlaminar interfaces in multi-ply multi-directional composite laminates. In contrast with the standard (penalty-like) method, the stabilizedmethod allows the use of arbitrarily large values of cohesive stiffness and obviates the need forengineering approaches to estimate minimum cohesive stiffness necessary for accurate delamina-tion analysis. This is achieved by defining a weighted interface traction in the stabilized method,which allows a gradual transition from penalty-like method for soft elastic contact to Nitsche-likemethod for rigid contact. We conducted several simulation studies involving constant strain patchtests and benchmark delamination tests under mode-I, mode-II and mixed-mode loadings. Ourresults show clear evidence of traction oscillations with the standard method with structured andperturbed finite element meshes, and that the stabilized method alleviates these oscillations, thusillustrating its robustness.
Keywords:
Nitsche method, Mixed-mode delamination, Traction oscillations, Numericalstability, Cohesive zone model
1. Introduction
In laminated fiber-reinforced composites, delamination is one of the most dominant failuremechanisms, which involves progressive damage accumulation and fracture along interlaminar in-terfaces [1]. Delamination under static and cyclic-fatigue loading has been widely studied in theliterature over the past 40 years [2, 3], because it causes localized damage that is hard to detect andmay lead to sudden structural collapse. The cohesive zone modeling approach has been extensivelyused to analyze and predict mixed-mode delamination propagation, despite its drawbacks and lim-itations. A particular drawback of the standard finite element implementation of cohesive zonemodels (CZMs) is its occasional numerical instability, which causes spurious traction oscillationsat delamination/crack interface [4]. Simple engineering solutions [5] may mitigate numerical is-sues with the standard finite element method (FEM) on a case-by-case basis, but they are not robustand introduce parametric uncertainty [6]. We recently illustrated that a Nistche-based stabilized ∗ Corresponding author
Email address: [email protected] (Ravindra Duddu) a r X i v : . [ c s . C E ] A ug EM is robust and accurate for enforcing stiff cohesive laws and simulating fracture propagationin isotropic, homogeneous elastic materials [7]. The purpose of this article is to investigate thenumerical stability and accuracy of standard FEM and weighted Nitsche-based stabilized FEMfor delamination analysis of composites using cohesive elements, especially at anisotropic anddissimilar interlaminar interfaces.The cohesive zone modeling approach, which is based on continuum damage mechanics [8],has been widely used to simulate mixed-mode delamination of composites [9–11], including thegrowth of multiple delamination cracks extending into various ply interfaces [12, 13]. Despite itssuccess, the standard FEM implementation of CZMs comes with certain outstanding challenges,including mesh dependence, parametric uncertainty, computational efficiency, and numerical in-stability. Mesh dependence of the predicted crack path or directional mesh bias can be an issuewith CZMs, as cracks can only propagate along finite element edges. To allow the propagationof arbitrary cohesive cracks and/or the inclusion of intra-element cohesive interfaces, various ap-proaches based on the partition of unity concept (including G/XFEM) or virtual/phantom nodeswere proposed [14–21]. The performance of the standard FEM to simulate cohesive cracks canbe poor with distorted or low quality meshes, so mesh free methods were proposed to alleviatesuch difficulties [22]. Recently, phase-field damage models have also been proposed to simulatedelamination of orthotropic laminates [23, 24].Parametric uncertainty or sensitivity of CZMs to the choice of model parameters can affectthe accuracy of delamination analysis. Depending on the choice of cohesive stiffness, cohe-sive/interface strength, and fracture toughness parameters, the accuracy of load and interface trac-tion prediction can vary significantly [6, 25]. Out of these three parameters, only fracture toughnesscan be determined or well-constrained from experiments for a given composite material. Notably,in the case of mixed-mode delamination at dissimilar interfaces in bi-materials, fracture tough-ness is also dependent on the mode mixity (loading) [26]. In contrast, the cohesive strength isdifficult to determine from experiments, and is usually assumed based on mesh size considera-tions [5, 8, 27], as a compromise between computational efficiency and numerical accuracy. Thecohesive stiffness is generally regarded as a numerical penalty parameter in intrinsic CZMs thatassume an initially elastic response, and in extrinsic CZMs that assume initially rigid response andelastic unloading/reloading response. The FEM implementation of extrinsic CZMs requires ad-vanced algorithms [28–30], which increases computational complexity; whereas, intrinsic CZMsare relatively straightforward to implement within a legacy/commercial finite element code, butthe choice of cohesive stiffness can affect numerical stability and/or convergence. It is noteworthythat extrinsic CZMs under fatigue and compressive loading scenarios also suffer from numericalinstabilities observed in intrinsic CZMs with large values of cohesive stiffness.In quasi-static fracture/delamination analysis, the standard FEM implementations using intrin-sic CZMs can exhibit spurious traction oscillations along the cohesive interface, especially nearcrack tips, if a large initial cohesive stiffness is specified [5, 31–36]. Even with potential-based in-trinsic CZMs, it is important to control the elastic behavior through initial slope indicators to avoidinstability [37]. Past studies indicate that using the Gaussian full integration scheme [1, 4, 31, 33]can cause spurious oscillations in the traction profile along cohesive interfaces. Although theNewton-Cotes integration scheme was suggested as an alternative, some studies reported that it2an result in an odd wrinkling mode in thin laminates undergoing delamination process [38–40].Improved cohesive stress integration schemes for continuum CZMs [41] were developed to furthertackle the issues of stability and robustness. Alternatively, discrete CZMs were developed [42–45] that use spring-like elements to connect the finite element nodes at delamination interfaces,instead of element edges, to alleviate numerical instability and/or convergence issues. However,in discrete CZMs, when using non-uniform meshes, or modeling interfacial kinks, relating force-displacement relations for the spring like elements with interface traction is not straightforward.The issue of numerical instability in the standard FEM implementation of intrinsic continuumCZMs, when using large cohesive stiffness and full/reduced Gauss integration schemes, arises dueto ill-conditioning of discrete systems, as typical with penalty-like formulations [46]. AlthoughLagrange-multiplier-based mixed or two-field formulations for intrinsic/extrinsic CZMs [47, 48]can overcome numerical instability, they can be computationally expensive or cumbersome owingto the difficulty in determining a stable Lagrange multiplier space.To broadly address the numerical instability issues with penalty-like formulations and standardFEM for interface/contact problems, discontinuous Galerkin (dG) methods or Nitsche-based meth-ods were proposed in the last two decades [49, 50]. Several novel dG approaches were developedfor fracture problems, including dG interface [51–53], space-time dG approaches [54, 55], andhybrid dG-CZM approaches [56–59], where the dG method is generally used before fracture initi-ation and CZMs are used for simulating fracture propagation. A limitation, however, for the wideruse of dG approaches for fracture is the inherent complexity associated with their implementation,especially in legacy/commercial finite element codes/software. Conversely, Nitsche methods wereadvocated by [60, 61] for modeling strong/weak discontinuities and elastic interface problems.Nitsche methods overcome the issue of numerical instability affecting penalty-like formulations,by adding consistency terms [62, 63]. Nitsche’s method has been extended for modeling frictional-sliding on embedded interfaces [64, 65] and small-sliding contact on frictional surfaces, includingstick?slip behavior [66]. Inspired by the work of [67], we recently extended the Nitsche’s methodto cohesive fracture problems, and developed a stabilized FEM that alleviates traction oscillationswith stiff, anisotropic cohesive laws [7].In this paper, we illustrate the ability of the stabilized FEM of Ghosh et al. [7] in alleviatingtraction oscillations at interlaminar interfaces in multi-directional orthotropic composite laminatesunder different loading conditions. A specific aim is to illustrate its robustness for composite de-lamination analysis, with regard to the choice of the cohesive stiffness and the structure of the finiteelement mesh (e.g. uniform structured versus perturbed or semi-structured meshes), which has notbeen addressed before. The rest of this paper is organized as follows: in Section 2, we brieflydescribe the governing equations of the cohesive fracture problem and the weak forms correspond-ing to the standard and stabilized methods. In Section 3, we discuss the salient aspects of thenumerical implementation and the selected model/material parameters. Specifically, we focus onthe implementation via user element subroutines in commercial finite element software ABAQUS[68], so that it can be utilized by the broader composite modeling community. In Section 4, wepresent several benchmark numerical examples to compare the standard and stabilized methods,with a particular emphasis on the accuracy of the interface traction field and load–displacementcurves. Finally, in Section 5, we conclude with a summary and closing remarks.3 . Governing equations and weak formulations
In this section, we briefly review the Nitsche-inspired stabilized finite element method orig-inally proposed in [7] for enforcing stiff cohesive laws. We will begin with a description ofthe strong form of the governing equations followed by the anisotropic bilinear cohesive law formixed-mode loading. Subsequently, we will discuss the weak form for the standard and stabilizedmethods, and the choice of stabilization parameters and weights.
We define an initial domain Ω ⊂ R , which is partitioned into two non-overlapping bulk do-mains Ω (1) and Ω (2) separated by a pre-defined internal cohesive interface Γ ∗ , such that Ω = Ω (1) ∪ Ω (2) (see Fig. 1). Throughout this paper, we use the notation that numbers within parentheses inthe superscript identify the domain partitions. Dirichlet and Neumann boundary conditions are en-forced on two disjointed parts of the domain boundary Γ ≡ ∂ Ω in such a way so that ∂ Ω = Γ D ∪ Γ N with Γ D ∩ Γ N = ∅ . The outward unit normal to the boundary ∂ Ω is denoted by n e and unit normalvector associated with the interface boundary Γ ∗ is denoted by n and points from Ω (2) to Ω (1) (thus n = − n (1) = n (2) ). ̅𝐭𝜕Ω = Γ ! ∪ Γ " Γ " 𝐧 (𝐗) Γ ∗ Ω = Ω (&) ∪ Ω (() x ! x X Ω (&) Ω (() 𝐧Γ ! , u ̅𝐭 𝐧Ω (&) Ω (() Γ ∗ , u Γ " Γ ! X ! X X 𝐧 Ω = Ω (&) ∪ Ω (() 𝜕Ω = Γ ! ∪ Γ " Figure 1: A schematic of the undeformed domain for the quasi-static delamination/debonding problem. We chooseX and X as the in-plane coordinates for the laminate material, and X is the out-of-plane coordinate. For delamination and debonding at sharp interfaces, we assume the bulk domains consist of ahomogeneous anisotropic linearly elastic material, thus intralaminar damage is neglected and onlyinterlaminar damage is considered. Assuming small displacements, the Cauchy stress tensor canbe defined in Ω (1) and Ω (2) as σ ( m ) = D ( m ) : (cid:15) ( m ) , m = { , } , (1)where D denotes the fourth-order anisotropic elasticity tensor and the small strain tensor (cid:15) = ( ∇ u + ( ∇ u ) T ) is defined by the symmetric part of the displacement gradient tensor. The gov-4rning elasto-static equilibrium equations in the absence of body forces are given by: ∇ · σ ( m ) = in Ω ( m ) , m = { , } , (2) u = ¯ u on Γ D , (3) σ · n e = ¯ t on Γ N , (4) t c = − σ (1) · n (1) = σ (2) · n (2) on Γ ∗ , (5)where ¯ u is the prescribed displacement vector on the Dirichlet boundary Γ D and ¯ t is the prescribedtraction on the Neumann boundary Γ N . The interface traction t c is related to the Cauchy stresstensor evaluated in the sub-domains Ω (1) and Ω (2) as given by (5), and is continuous across thecohesive interface Γ ∗ to satisfy the force equilibrium. The cohesive traction ˆ t c is the Newton’sthird law pair to the interface traction t c on a given delamination/crack surface and can be definedas a function of the interface separation or displacement jump as ˆ t c = − t c = α ( δ ) δ , (6) δ = [[ u ]] = u (2) − u (1) , (7)where the cohesive stiffness matrix α is usually a nonlinear function of the interface separation. Remark 1.
In the case where the interlaminar interface delineates two dissimilar materials (i.e., D (1) (cid:54) = D (2) ), the Cauchy stress tensor evaluated at interface from either side can be different (i.e., σ (1) (cid:54) = σ (2) ), but the traction field must be continuous across the interface.2.2. Anisotropic bilinear cohesive law We consider an intrinsic bilinear traction–separation or cohesive law that consists of an initialelastic region followed by a softening region. For mixed-mode delamination under quasi-staticloading in two-dimensions, we cast the bilinear cohesive law in the damage mechanics frameworkas detailed in [6]. The tangential t τ and normal t n components of the interface traction vector t c are related to the tangential δ τ and normal δ n components of the interface separation δ as t c = (cid:26) t τ t n (cid:27) = − (1 − d s ) α τ (cid:18) − d s (cid:104) δ n (cid:105) δ n (cid:19) α n (cid:26) δ τ δ n (cid:27) , (8)where α n and α τ represent the initial cohesive stiffness in the normal and the tangential directions,respectively, and the scalar damage variable d s is given by d s = if δ e < δ ce ,δ ue ( δ e − δ ce ) δ e ( δ ue − δ ce ) if δ ce ≤ δ e < δ ue , if δ ue ≤ δ e , (9)and the equivalent separation δ e = (cid:112) (cid:104) δ n (cid:105) + δ τ . In the above equations, (cid:104)·(cid:105) denotes Macaulaybrackets, so that (cid:104) δ n (cid:105) = max (0 , δ n ) , which ensures that there is no damage growth or damage5ffect on the normal cohesive stiffness response under compression or contact. The critical andultimate interface separation parameters δ ce and δ ue , respectively, defined as [10]: δ ce = (cid:115)(cid:18) α n cos Iσ max (cid:19) + (cid:18) α τ cos IIτ max (cid:19) , (10) δ ue = (cid:18) α n δ ce (cos I ) G IC (cid:19) + (cid:18) α τ δ ce (cos II ) G IIC (cid:19) , (11)where the direction cosines cos I = δ n /δ e and cos II = δ τ /δ e , σ max and τ max are the pure mode Iand mode II cohesive strengths, and G IC and G IIC are the pure mode I and mode II critical fractureenergies. For illustration, the normal and tangential traction profiles as a function of the normaland tangential interface separations are shown in Fig. 2.
Remark 2.
Although the anisotropic bilinear cohesive law is phenomenological, it can be relatedto a potential function defined as
Ψ = − (cid:18) (1 − d s ) α τ δ τ + (cid:18) − d s (cid:104) δ n (cid:105) δ n (cid:19) α n δ n (cid:19) . (12) The crack surface traction components can be defined based on the above potential function as t τ = ∂ Ψ ∂δ τ and t n = ∂ Ψ ∂δ n . The expressions for (8) are approximations obtained by neglecting thenonlinearity due to interface damage d s , which is a function of the interface separation. Remark 3.
The anisotropic bilinear cohesive law of Jiang et al. [10] has six independent param-eters, namely initial cohesive stiffness, maximum cohesive strength and critical fracture energy ofpure mode I and II loadings, to describe the traction-separation relationship. If these parame-ter values are chosen to be the same for both normal and shear modes, then we get an isotropicbilinear cohesive law with only three independent parameters.2.3. Standard weak form
We apply the Galerkin procedure of weighted residuals to derive the standard weak form, which isdetailed in [7]. By weighting the equilibrium equation in (2) by a test function w , integrating byparts, applying the divergence theorem, using the traction continuity condition at the interface in(5), and the constitutive relation in (1), we obtain the weak form as: (cid:88) m =1 (cid:90) Ω ( m ) ∇ s w ( m ) : D ( m ) : ∇ s u ( m ) d Ω − (cid:90) Γ ∗ [[ w ]] · t c d Γ = (cid:90) Γ N w · ¯ t d Γ , (13)where the jump in the test function is defined as [[ w ]] = w (2) − w (1) . Substituting the traction-separation relation in (6) into the weak form in (13) we obtain the standard weak form as6 t n (a) n t (b) Figure 2: Traction-separation relations defined by the anisotropic, intrinsic bilinear cohesive law of [10]: (a) normaltraction field and (b) tangential traction field. (cid:88) m =1 (cid:90) Ω ( m ) ∇ s w ( m ) : D ( m ) : ∇ s u ( m ) d Ω + (cid:90) Γ ∗ [[ w ]] · α ( δ ) δ d Γ = (cid:90) Γ N w · ¯ t d Γ . (14)Because the cohesive traction and separation components are defined in the normal and tangentialdirections, the standard weak form is implemented as, (cid:88) m =1 (cid:90) Ω ( m ) ∇ s w ( m ) : D ( m ) : ∇ s u ( m ) d Ω + (cid:90) Γ ∗ (1 − d s ) (cid:0) [[ w n ]] α n δ n + [[ w τ ]] α τ δ τ (cid:1) d Γ = (cid:90) Γ N w · ¯ t d Γ . (15)Thus, in the standard weak formulation, the cohesive traction is simply enforced as a mixedboundary condition on the interface. Remark 4.
If the initial cohesive stiffness parameters α n and α τ are taken to be large enough, thestandard weak form resembles the penalty method for enforcing displacement continuity acrossthe interface. However, for stiff cohesive laws, where cohesive stiffness is several orders of mag-nitude greater than the elastic modulus, the standard weak form becomes ill-conditioned leadingto numerical instability and/or convergence issues. In the limiting case of a non-interpenetration(contact) constraint or an extrinsic cohesive law, where α n → ∞ and/or α τ → ∞ , the standardweak form is not well defined.2.4. Stabilized weak form The stabilized FEM developed in [7] extended Nitsche’s method [67, 69] to cohesive fractureproblems, which we review here for the sake of clarity. The key idea is to evaluate the interface7raction in terms of a weighted average stress in the bulk material across the interface and thetraction in the cohesive interface as t c = ( I − S ) (cid:104) σ (cid:105) γ · n − S α δ , (16)where I is the second-order identity matrix, S is the stabilization matrix defined as, S = β τ α τ (1 − d s ) + β τ β n α n (1 − d s ) + β n , (17) β τ , β n are the stabilization parameters, and (cid:104) σ (cid:105) γ is the weighted average of the stress tensors onboth sides of the interface defined as (cid:104) σ (cid:105) γ = ( γ (1) σ (1) + γ (2) σ (2) ) ∀ γ (1) + γ (2) = 1 , γ (1) > , γ (2) > . (18)Substituting (16) into the weak form in (13) we obtained the stabilized weak form as (cid:88) m =1 (cid:90) Ω ( m ) ∇ s w ( m ) : D ( m ) : ∇ s u ( m ) d Ω − (cid:90) Γ ∗ [[ w ]] · ( I − S ) (cid:104) σ (cid:105) γ · n d Γ+ (cid:90) Γ ∗ [[ w ]] · S α ( δ ) δ d Γ = (cid:90) Γ N w · ¯ t d Γ . (19)In the above equation, the second and third terms on the left hand side ensure consistency andstability, respectively. If the stabilization matrix is taken as the identity matrix, then the stabilizedweak form in (19) becomes identical to the standard weak form in (14). Remark 5.
The stabilized method presented here is unsymmetric and resembles the incompleteinterior penalty method [70, 71]. It can be proved that the displacement solution u of the strongform equations (2) – (5) is satisfied by the solution to the weak form equation (19), which estab-lishes consistency for any value of cohesive stiffness; the mathematical procedure for proving thisis similar to that described in [67, Lemma 2.1]. The stabilized weak form with the interface traction and separation components expressed inthe normal and tangential coordinates is given by (cid:88) m =1 (cid:90) Ω ( m ) ∇ s w ( m ) : D ( m ) : ∇ s u ( m ) d Ω − (cid:90) Γ ∗ [[ w ]] · ( I − S ) (cid:104) σ (cid:105) γ · n d Γ+ (cid:90) Γ ∗ (cid:18) [[ w n ]] (1 − d s ) α n β n α n (1 − d s ) + β n δ n + [[ w τ ]] (1 − d s ) α τ β τ α τ (1 − d s ) + β τ δ τ (cid:19) d Γ = (cid:90) Γ N w · ¯ t d Γ . (20)8 emark 6. As (1 − d s ) α n , (1 − d s ) α τ → ∞ , the stabilized weak form in (20) resembles the Nitschestabilized finite element method for frictional contact presented in [65]. Thus, the stabilized weakform remains well-defined in the limiting case of a non-interpenetration (contact) constraint or anextrinsic cohesive law, unlike the standard weak form in (15).2.5. Choice of stabilization parameters and weight factor The stabilization parameters β τ , β n and the weights γ (1) , γ (2) play a key role in the numericalperformance of a Nitsche-based stabilized FEM [72]. For instance, taking a small value for thestabilization parameters may undermine the positive definiteness of the global linear system ofequations; whereas, taking a large value for the stabilization parameters will essentially lead to apenalty-like method for enforcing the interface constraints [73]. Moreover, taking equal weightsfor dissimilar material interfaces may hamper numerical performance, or may not alleviate tractionfield oscillations. For constant strain triangular and tetrahedral elements, Annavarapu et al. [72]provided estimates for the stabilization parameters using a local coercivity analysis as given by β n = β τ = 2 (cid:32) | D (1) | ( γ (1) ) meas (Ω (1) ) + | D (2) | ( γ (2) ) meas (Ω (2) ) (cid:33) meas (Γ ∗ ) . (21)Here, we simply use the above estimate of stabilization parameters for bilinear quadrilateral finiteelements and conduct parametric sensitivity studies to illustrate their adequacy. Equation (21)establishes a functional dependency between the stabilization parameter and interface weights, andthe weights are chosen so that it minimizes the stabilization parameter while ensuring coercivity.Here, we simply use expression for interface weights derived in [72]: γ (1) = meas (Ω (1) ) | D (1) | meas (Ω (1) ) | D (1) | + meas (Ω (2) ) | D (2) | ; γ (2) = 1 − γ (1) , (22)where | D | denotes the two-norm of the elasticity tensor, meas (Ω) denotes the area of neighboringbulk element in 2D, and meas (Γ ∗ ) is the length of the interface element. Remark 7.
For the weak form in (19) , precise estimates for the stabilization parameter can bederived for constant strain elements following the procedure described in [72]. For higher-orderelements, closed-form analytical estimates for the stabilization parameters are yet to be derived.However, the stabilization parameters can be specified by solving a local eigenvalue problem [74].
3. Numerical Implementation and Model Parameters
We implemented the stabilized FEM in the commercial software ABAQUS through user-defined subroutines. In this section, we briefly present key details of ABAQUS implementationand list the model parameters that are specific to delamination analysis. The full details of the nu-merical implementation (omitted here), including the finite element approximation, discretization9nd linearization of the standard and stabilized weak forms, and the expressions for continuum andinterface element force vectors and matrices can be found in [7].
We use bilinear quadrilateral four-noded plane stress continuum elements with four-point Gaussintegration scheme and four-noded linear zero-thickness interface elements with two-point Gaussintegration scheme. Although our method can be implemented in 3-D, we use the 2D plane stressapproximation as it has been extensively used in prior studies and has been validated with experi-mental data [11, 27, 37, 52]. We chose the finite element mesh and interface element size so thatthere are at least three interface elements within the estimated cohesive process zone; this is nec-essary for an accurate representation of the numerical stress distribution within the process zone atthe point of initial crack propagation, as elaborated in [27]. User-element subroutines in ABAQUStypically require the user to provide the stiffness matrix (AMATRX) and the right hand side (RHS)force vector. In our implementation, we utilize the UELMAT subroutine to provide the continuumelement force vector and stiffness matrix, and the UEL subroutine to provide the cohesive elementforce vector and stiffness matrix. The UELMAT subroutine allows the user to access some of theinbuilt material models through utility subroutines
MATERIAL_LIB_MECH , unlike the UMATsubroutine. Using global modules, we store and share the stress and shape function derivativematrices calculated in the UELMAT subroutine to the UEL subroutine for computing interfaceforce vector and stiffness matrix. A detailed description of these computations can be found in [7]for isotropic elasticity; whereas, here we use anisotropic elasticity for composites. For the sakeof verification or comparison in some simulation studies, we used ABAQUS in-built 4-noded 2Dcohesive elements (COH2D4) along with 4-noded plane stress 2D continuum elements (CPS4).
The anisotropic bilinear cohesive law described in Section 2.2 is highly nonlinear, because thescalar damage variable is a complex nonlinear function of the normal and tangential separations.The consistent tangent stiffness corresponding to this cohesive law can be derived from (8) as K tan = − ∂t τ ∂δ τ ∂t τ ∂δ n ∂t n ∂δ τ ∂t n ∂δ n = (1 − d s ) α τ − α τ δ τ ∂d s ∂δ τ − α τ δ τ ∂d s ∂δ n − α n δ n (cid:104) δ n (cid:105) δ n ∂d s ∂δ τ (cid:18) − d s (cid:104) δ n (cid:105) δ n (cid:19) α n − α n δ n (cid:104) δ n (cid:105) δ n ∂d s ∂δ n , (23)Deriving and implementing a closed-form expression of the above consistent tangent is arduous,and using it can cause numerical issues, as the diagonal terms become negative in the softeningportion of the cohesive law. In our previous work [6, 44, 75], we argued that the secant stiffness ismore advantageous with this cohesive law and demonstrated its accuracy and convergence with animplicit scheme. The simpler secant stiffness matrix can be derived from (8) as K sec = (1 − d s ) α τ (cid:18) − d s (cid:104) δ n (cid:105) δ n (cid:19) α n , (24)10lthough we include the simple, linearized secant stiffness terms in the AMATRX defined inthe UEL subroutine, the nonlinearity of the cohesive fracture problem is handled in ABAQUS/Stan-dard outside of the user subroutines. As detailed in the user manual [68, Chapter 7: Analysis Solu-tion and Control], ABAQUS/Standard combines incremental and iterative (Newton-Raphson) pro-cedures for solving nonlinear problems. Using the secant stiffness necessitates smaller load/disp-lacement increments (pseudo-time steps) to gradually reach the final applied load/displacement.The user typically suggests the maximum and minimum increment size and the size of the firstincrement, and ABAQUS/Standard automatically chooses the size of the subsequent increments.Within each increment, ABAQUS/Standard automatically performs iteration to find an equilibriumsolution based on a user-defined criteria for residual force and displacement correction. In all ofour simulations, we use sufficiently small increment size that most displacement steps converge ina single iteration, except at certain time steps where more than one cohesive element fails. Also,in cohesive fracture simulations, the ABAQUS/Standard default criteria for residual force toler-ance may be too small that numerical convergence may not be attainable as cohesive elementsfail. Especially, using the secant stiffness it becomes necessary to increase these tolerances ap-propriately (sometimes by two orders of magnitude than the default value) to attain convergence.Therefore, we compare the predicted load–displacement curves against analytical solutions, exper-imental data, other numerical model results to ensure the accuracy of our simulations. We consider HTA/6376C unidirectional carbon-fiber-reinforced epoxy laminate as the genericcomposite material. Here, we only consider delamination or debonding along 2D straight inter-faces between two HTA/6376C lamina with the fibers aligned either in the X (in-plane horizontal)or X (out-of-plane) direction, as indicated in Figure 1. According to the standard notation usedto define stacking sequences in multi-ply composites laminates [76], “[0/0]” laminate denotes thetwo-ply specimen with fibers in the top and bottom lamina oriented in X direction. Similarly,“[0/90]” laminate denotes the cross-ply specimen with the fibers in the top lamina oriented in the X direction and fibers in the bottom lamina oriented in the X direction. Usually, this notationimplies that each laminate layer has the same thickness and made of the same composite mate-rial (HTA/6376C in our study). Unfortunately, experimental data from delamination tests is onlyavailable for [0/0] laminates. The anisotropic (transversely isotropic) linear elastic material prop-erties for the 0 ◦ ply HTA/6376C laminate are listed in Table 1, which is directly obtained fromexperiments [77]. Using coordinate transformation relations we can easily obtain these materialproperties for the 90 ◦ ply HTA/6376C laminate. Table 1: Material properties of carbon fiber/epoxy laminated composite HTA/6376C obtained from [77] . E E = E G = G G ν = ν ν (N/mm ) (N/mm ) (N/mm ) (N/mm ) . × . × . × . × .4. Cohesive zone model parameters The CZM parameters chosen in our simulation studies are listed in Table 2. The mode I andmode II fracture energies for the 0 ◦ ply HTA/6376C laminate are taken from [77]. For ensuringthe accuracy and convergence of delamination analysis using the FEM with cohesive elements,two conditions must be satisfied [5]: (1) the element size must be less than the cohesive (process)zone length, which is determined by fracture toughness and cohesive strength; and (2) the cohesivestiffness must be large enough to avoid the introduction of artificial compliance. The cohesivestrength is often chosen based on cohesive zone length and mesh size considerations [see Eq. (7) inRef. 27], owing to computational cost or limitations. Choosing a small cohesive strength may yielda poor peak load prediction, but beyond a certain value choosing a larger cohesive strength will notimprove model fit with the load–displacement data from quasi-static delamination tests, but willseverely restrict the mesh size and increase computational cost. Here, we use the cohesive strengthvalues suggested in [27] to ensure numerical accuracy and efficiency. Due to the unavailabilityof experimental data for the 90 ◦ ply HTA/6376C laminate, we assume the same CZM parameterslisted in Table 2.The cohesive stiffness is generally considered to be a penalty parameter and various guide-lines have been proposed in the literature for selecting the stiffness. Although the purpose of thecohesive stiffness in intrinsic/extrinsic cohesive zone models is to account for the elastic loading,unloading and reloading response of the fracture/delamination interface, it can contribute to theglobal deformation response and introduce artificial compliance or numerical instability issues.Based on 1D laminate model, the cohesive stiffness necessary to avoid artificial compliance issuecan be estimated as [5] α = E M/t, (25)where E is the Young’s modulus of the material along the laminate thickness direction, t is thesub-laminate thickness , and M is a non-dimensional number that is to be chosen much largerthan one. For cohesive interfaces in non-laminates t is not defined, so in (25) it can be replacedby a certain length measure h of the bulk material [78] or the finite element mesh size. Taking M = 100 , E = E = 1 . × N/mm and t = 1 . mm for the delamination tests (see Section4), we estimate α ≈ N/mm . To demonstrate the performance of the standard and stabilizedmethods, we assume three values of cohesive stiffness in our studies, including values that are twoorders of magnitude smaller and larger than the above estimate. However, we note that for thin-plylaminates with t < . mm or for small values of length measure h < . mm in non-laminates,the estimated cohesive stiffness α > N/mm . Table 2: Cohesive zone model parameters for the carbon fiber/epoxy laminated composite HTA/6376C are taken from[27], except the cohesive stiffness values. α n α τ G IC G IIC σ max τ max (N/mm ) (N/mm ) (N/mm) (N/mm) (N/mm ) (N/mm ) { , , } { , , } . Numerical Examples In this section, we present several examples to demonstrate the ability of the Nitsche-inspiredstabilized formulation in alleviating oscillations in interface traction using constant strain patchtests, and pure mode I, mode II and mixed mode delamination tests. Through these tests, wespecifically examine numerical stability at similar and dissimilar laminate interfaces defined byanisotropic and isotropic cohesive laws using perturbed, structured and unstructured meshes.
We assess the ability of the standard and stabilized formulations in alleviating traction oscil-lations at horizontal and inclined straight interfaces using the constant strain patch test. Undercompressive loading, we assume a stiff elastic response in the normal direction to enforce contactand a weak elastic sliding response in the tangential direction, which is captured by the anisotropicCZM. We assign the square plate with a side length L = 1 mm and the horizontal delaminationinterface at mid-height. To apply the compressive load, we constrain both vertical and horizon-tal displacements at the bottom edge of the plate, and prescribe a uniform vertical displacement ∆ = − . mm at the top edge of the plate. We specify traction-free conditions at the left and rightedges of the square plate. The cohesive parameters and material properties assumed for this testare listed in Tables 1 and 2, respectively. To examine mesh sensitivity, we generate a × structured square mesh with element lengthof 0.1 mm (Fig. 3a) and perturb the interface nodes by ≈ of the element length (Fig. 3b). In Fig.4, we show the normal traction profile versus the horizontal coordinate along the [0/0] laminateinterface obtained from the standard and stabilized methods. We consider high stiffness-contrastwith the anisotropic CZM, where α n = 10 N/mm and α τ = 10 N/mm . According to (21), wetake the stabilization parameters β n = β τ = 3 × N/mm . As shown in Fig. 4a, if the interfaceis perfectly flat, then both standard and stabilized methods yield smooth traction profiles withoutany spurious oscillations. However, we note that the standard FEM exhibits spurious traction os-cillations even with the unperturbed mesh, if the normal stiffness α n = α τ ≥ N/mm (resultsnot shown here). From Fig. 4b, we can see that the standard FEM suffers from severe numericalinstability even with a slightly perturbed interface mesh, as evident from the large amplitude trac-tion oscillations. In contrast, our stabilized FEM is stable and alleviates traction oscillations foranisotropic CZMs with high stiffness contrast between normal and tangential directions. We consider a straight interface inclined at an initial angle of 140.4 ◦ with the horizontal X axis embedded within the square plate of side length L = 1 mm (Fig. 5a). We discretize thedomain using a × semi-structured mesh with quadrilateral element, so that the interfaceis divided into 13 elements (Fig. 5b). According to (21), we take the stabilization parameters β n = β τ = 5 × , × N/mm for [0/0] and [0/90] laminate interfaces, respectively. Weconsider the anisotropic CZM with high stiffness contrast, where α n = 10 N/mm and α τ = 10 N/mm . In Fig. 6, we show the normal traction profile versus the horizontal X coordinate of13 a) (b) ∆ x x x x x x x x x x x ∆ x x x x x x x x x x x X X Figure 3: Boundary conditions and mesh used for the square plate with horizontal interface: (a) straight interface; (b)perturbed interface. The nodes are perturbed by ≈ % of the element size and a zoom of the interface undulations isshown in the inset. (a) (b)StandardStabilized Figure 4: Traction profiles obtained from the standard and stabilized methods with the anisotropic CZM ( α n = 10 and α τ = 10 N/mm ) for the square plate made of [0/0] laminate: (a) horizontal interface (b) perturbed interface. the integration points along the inclined interface for [0/0] and [0/90] laminates obtained from thestandard and stabilized methods. While the standard FEM exhibits instability evident from thelarge amplitude traction oscillations, the stabilized FEM is able to alleviate oscillations and yieldsa smooth traction profile. This study illustrates the drawback of the standard formulation for semi-structured meshes and potentially unstructured meshes when using anisotropic CZMs with highstiffness contrast. Overall, the two patch tests highlight the robustness of our stabilized FEM, butin the following sections we will investigate its performance for benchmark delamination tests. We investigate the ability of the stabilized FEM in recovering oscillation-free interface tractionduring mode-I delamination crack growth, using the double cantilever beam (DCB) test. Thespecimen geometry and test set-up along with the finite element mesh are shown in Fig. 7. To14 ∆ (a) (b) x x x x x x x x x x x x x x Figure 5: Boundary conditions and mesh used for the square plate with inclined interface: (a) schematic diagram; (b)structured finite element mesh with quadrilateral (non-rectangular) elements. (a) (b)StandardStabilized
Figure 6: Traction profiles obtained from the standard and stabilized methods with the anisotropic CZM ( α n = 10 and α τ = 10 N/mm ) for the square plate with an inclined interface: (a) [0/0] laminate (b) [0/90] laminate. initiate the delamination process, a pre-crack is placed at the left end of the beam, and equal andopposite vertical displacements ( ∆ ) are applied on the upper and lower nodes. The correspondingload (P) is determined from the simulation using the reaction force at the corresponding node. Thefixed boundary condition is applied at the right end of the beam.In Fig. 8a, we compare the load–displacement curves obtained from the stabilized formula-tion for [0/0] laminate with experimental data from [79], to check that our choice of strength anddisplacement increment is appropriate and to determine the sensitivity to cohesive stiffness. Evi-dently, there is good agreement between the numerically predicted load-displacement response andthe experimental data. The slight mismatch in the predicted peak load and the initial slope of theload–displacement curves is plausibly due to the idealization of boundary and loading conditions inthe model, compared to the experimental test setup. In Fig. 8b, we compare the load–displacementcurves obtained from the stabilized FEM for [0/90] laminate with those obtained from ABAQUSinbuilt cohesive elements [68], as experimental data is unavailable. The perfect match of load–displacement curves obtained with our user-defined and ABAQUS inbuilt (COH2D4) elements15erifies the correctness of our implementation. The load–displacement curves obtained from thestandard FEM exactly match with those from the stabilized FEM, so we do not show them here.As the two methods differ primarily in their ability to recover interface traction fields, we examinetraction fields along the delamination interfaces obtained from standard and stabilized methods for[0/0] and [0/90] laminates. ∆, P∆, PH a L Figure 7: Geometry and boundary conditions for the double cantilever beam test. The dimensions are: L = 150 mm,H = 3.1 mm and a = 35 mm. Displacement (mm) L oa d ( N ) Experimental α n = 10 N/mm α n = 10 N/mm α n = 10 N/mm (a) Displacement (mm) L oa d ( N ) ABAQUS Inbuilt α n = 10 N/mm α n = 10 N/mm α n = 10 N/mm (b) Figure 8: Load versus displacement curves for the double cantilever beam test obtained from the stabilized FEM: (a)[0/0] laminate and (b) [0/90] laminate.
From Fig. 9a, we find that for smaller values of stiffness α n = 10 , N/mm the normaltraction profile from the standard FEM is smooth, but for the larger value of α n = 10 N/mm significant oscillations are observed in the traction profile. From Fig. 9b, it is evident that thestabilized FEM is able to alleviate traction oscillations for the large value of cohesive stiffness.We also notice that the location of the region of traction oscillations (i.e. X ≈ mm) coincideswith the transition from tensile to compressive normal traction. To further examine the effect ofthis instability on interface damage, we plotted the respective damage profiles obtained from both16ethods in Figs. 9c and 9d. For α n = 10 N/mm at the location of traction oscillations, wealso find an oscillation in the damage profile with the standard FEM; whereas, there is no suchoscillation in the damage profile with the stabilized FEM. This study shows evidence of instabilitywith the standard FEM even with a perfectly flat interface and structured rectangular mesh. (a) (b)(c) (d) ↵ n = 10 ↵ n = 10 ↵ n = 10 N/mm N/mm N/mm Figure 9: Normal traction and damage versus interface length curves from the double cantilever beam test ([0/0] plyorientation) for different cohesive stiffness values: (a), (c) standard FEM and (b), (d) stabilized FEM.
We next examine the interface traction along the delamination interface recovered with a per-turbed rectangular mesh, where we change the coordinates of interface nodes by ≈ of theelement length. In Fig. 10a and 10b, we show the normal traction profiles obtained from the stan-dard FEM and a zoomed-in image showing evidence of spurious oscillations in the compressionregion of the traction. From Fig. 10c it is evident that the stabilized FEM is able to alleviate thetraction oscillations. This study demonstrates that numerical instability with the standard FEM ismore pronounced under contact conditions (i.e. in the regions where the normal traction is neg-ative) when using perturbed or unstructured finite element meshes, although the amplitude of thespurious oscillations in the mode I DCB test is generally small.17 a) (b) StraightPerturbed (c) Figure 10: Normal traction versus interface length for α n = 10 N/mm obtained from DCB test ([0/0] ply orienta-tion): (a) standard FEM, (b) zoom in of standard FEM near the crack tip and (c) stabilized FEM. We now compare the interface traction fields recovered from standard and stabilized methodsat dissimilar material interfaces. Recall that the 0 ◦ laminate exhibits anisotropic material behaviorin the X and X plane; whereas, the 90 ◦ laminate exhibits isotropic material behavior in this planewith an order of magnitude contrast in elastic modulus ( E /E = 11 . ). In Fig. 11, we showthe normal traction along the delamination interface recovered from the standard and stabilizedmethods for different values of initial cohesive stiffness. We do find a minor oscillation in normaltraction and damage profiles in Figs. 11a and 11c with the standard FEM for α n = 10 N/mm ,which is not the case with the stabilized FEM. Notably, there are minor differences in the tractionand damage curves for different stiffness values, unless stiffness is taken greater than N/mm . We investigate the ability of the stabilized FEM in recovering oscillation-free interface tractionduring mode-II delamination crack growth, using the end notch flexure (ENF) test. The specimengeometry and test set-up are shown along with the finite element mesh in Fig. 12. To initiatethe delamination process, a pre-crack is placed at the left end of the beam and downward verticaldisplacement is applied at the middle of the simply supported beam. In Fig. 13a, we comparethe load–displacement curves obtained from the stabilized FEM for [0/0] laminate specimen alongwith the experimental data and the corrected beam theory solution from [27] for different valuesof cohesive stiffness. In Fig. 13b, we compare the load–displacement curves obtained from thestabilized FEM for [0/90] laminate with those obtained from ABAQUS inbuilt cohesive elements[68], due to the unavailability of experimental data or analytical solutions. The good match ofload–displacement curves, including the predicted peak load, obtained with our user-defined andABAQUS inbuilt (COH2D4) elements verifies the correctness of our implementation. The load–displacement curves obtained from the standard FEM exactly match with those from the stabilizedFEM, so we do not show them here. 18 a) (b)(c) (d) ↵ n = 10 ↵ n = 10 ↵ n = 10 N/mm N/mm N/mm Figure 11: Normal traction and damage versus interface length curves from the double cantilever beam test ([0/90] plyorientation) for different cohesive stiffness values: (a), (c) standard FEM and (b), (d) stabilized FEM.
From Fig. 14a, it is evident that the tangential traction profile obtained from the standard FEMis oscillation-free for smaller cohesive stiffness values α τ = 10 , N/mm , but for the largervalue α τ = 10 N/mm a significant oscillation can be seen. Also, the traction field for α τ = 10 from the standard FEM does not match with that for larger stiffness values, which emphasizesthe importance of choosing a large enough cohesive stiffness. The traction profile obtained fromthe stabilized FEM shown in Fig. 14b does not exhibit any oscillations for all cohesive stiffnessvalues, thus demonstrating stability and robustness. We also plot the respective damage profilesobtained from both the methods in Figs. 14c and 14d, which show a corresponding oscillation inthe damage profile for the standard FEM for α τ = 10 N/mm , but no such oscillation exists in thedamage profile with the stabilized FEM. These results indicate that even though the standard andstabilized methods capture the load displacement curves well, numerical instability can corrupt thetangential traction and damage field at the interlaminar interface.We next examined the tangential traction along the delamination interface recovered with a19 , PH a L Figure 12: Geometry and boundary conditions for the end notch flexure test. The dimensions are: L = 100 mm, H =3.1 mm and a = 35 mm. Displacement (mm) L oa d ( N ) ExperimentalCorrected Beam Theory α τ = 10 N/mm α τ = 10 N/mm α τ = 10 N/mm (a) Displacement (mm) L o a d ( N ) ABAQUS Inbuilt = 10 N/mm = 10 N/mm = 10 N/mm (b) Figure 13: Load versus displacement curves for the end notch flexure test obtained using the stabilized FEM: (a) [0/0]laminate and (b) [0/90] laminate. perturbed rectangular mesh, where we change the coordinates of interface nodes by ≈ of theelement length. In Fig. 15a and 15b, we show the tangential traction profiles obtained from thestandard FEM and a zoomed-in image showing clear evidence of spurious oscillations. From Fig.15c it is evident that the stabilized FEM is able to alleviate the large amplitude oscillations. Thisstudy demonstrates that numerical instability with the standard FEM can corrupt tangential traction(i.e. shear stress at the interface) when using perturbed or unstructured finite element meshes, andthe amplitude of some spurious oscillations in the mode II ENF test can be large (as much as 50%local error). Therefore, we recommend using the stabilized FEM to improve accuracy and avoidconvergence issues with stiff cohesive laws. We now examine the interface traction fields obtained from both standard and stabilized meth-ods for the [0/90] laminate interface. Although we expect a similar response to the [0/0] laminateinterface case, we present the results here for completeness. In Fig. 16, we show the tangentialtraction along the delamination interface recovered from the standard and stabilized methods fordifferent values of initial cohesive stiffness. We do find large amplitude oscillations in the tangen-20 a) (b)(c) (d) ↵ ⌧ = 10 N/mm ↵ ⌧ = 10 N/mm ↵ ⌧ = 10 N/mm Figure 14: Tangential traction and damage versus interface length curves from the end notch flexure test ([0/0] plyorientation) for different cohesive stiffness values: (a), (c) standard FEM and (b), (d) stabilized FEM. tial traction for α τ = 10 N/mm in Fig. 16a, and in the corresponding damage profile in Fig. 16cwith standard FEM, whereas no traction oscillations are seen with the stabilized FEM. We also investigated the ability of the stabilized FEM in recovering normal and tangentialtraction fields during mixed-mode delamination crack growth using the fixed ratio mixed mode(FRMM) test. The specimen geometry and test set-up along with the finite element mesh areshown in Fig. 17. To initiate the delamination process, a pre-crack is placed at the left end andvertical displacement is applied only at the upper corner node. The fixed boundary condition isapplied at the right end of the beam. In Fig. 18a, we compare the load–displacement curves fromthe stabilized FEM with the corrected beam theory solution from [27] for [0/0] laminate, due tothe lack of experimental data. In Fig. 18b, we compare the load–displacement curves from thestabilized FEM with the ABAQUS inbuilt cohesive elements for [0/90] laminate, due to the lack21 a) (b) StraightPerturbed (c)
Figure 15: Tangential traction versus interface length for α τ = 10 N/mm obtained from ENF test ([0/0] ply orienta-tion): (a) standard FEM , (b) zoom in near the crack tip for standard FEM and (c) stabilized FEM. of experimental data or analytical solutions. Because the ABAQUS inbuilt cohesive elements usea different mixed-mode criteria proposed in [9], the predicted peaks loads and softening portionsof the load–displacements are different in Fig. 18b. However, by choosing different cohesivestrengths with different mixed-mode criteria it is possible to match the load–displacement curvesobtained from stabilized FEM and inbuilt ABAQUS elements. In Figs. 19a and 19b, we show the normal and tangential traction profiles for [0/0] laminatefrom the standard FEM. We find that for the large cohesive stiffness α n = α τ = 10 N/mm ,large amplitude oscillations appear in the tangential traction, but a few small amplitude oscillationappears in the normal traction profile. The stabilized FEM alleviates oscillations in both normaland tangential traction fields, as evident from Figs. 19c and 19d. Notably, for the smaller stiffnessvalue of N/mm , the peak tangential traction value is considerably smaller than the other twocohesive stiffness cases (i.e., N/mm and N/mm ) from both methods, which also happenswith ABAQUS in-built cohesive elements. This suggests that it is important to take the cohesivestiffness large enough to accurately recover the interface traction for mixed mode loading. In Figs. 20a and 20b, we show the normal and tangential traction profiles for [0/90] laminatefrom the standard FEM. We find that for α n = α τ = 10 N/mm , large amplitude oscillationsappear in both normal and tangential traction profiles. The stabilized FEM alleviates these oscil-lations, as evident from Figs. 20c and 20d. We also notice for the smaller stiffness value of N/mm , the peak tangential traction value is considerably smaller than the other two cohesive stiff-ness cases (i.e., N/mm and N/mm ). Once again, this indicates the importance of takinga larger value for cohesive stiffness to accurately recover the interface traction.22 a) (b)(c) (d) ↵ ⌧ = 10 N/mm ↵ ⌧ = 10 N/mm ↵ ⌧ = 10 N/mm Figure 16: Tangential traction and damage versus interface length curves from the end notch flexure test ([0/90] plyorientation) for different cohesive stiffness values: (a), (c) standard FEM and (b), (d) stabilized FEM.
5. Conclusion
In this paper, we compared and contrasted the performance of standard and stabilized finite el-ement methods for composite delamination analysis using cohesive elements. First, we illustratedthrough constant strain patch tests that numerical instability with standard FEM causes spuriousoscillations in the interface traction fields at both similar and dissimilar composite material inter-faces. We demonstrated that perturbed, semi-structured and potentially unstructured meshes canaggravate the numerical instability with standard FEM, especially in anisotropic CZMs where thecontrast between normal and tangential stiffness is high. For constant strain patch tests, we then il-lustrated the ability of a stabilized method inspired by the weighted Nitsche approach in alleviatingspurious traction oscillations at the delamination interface.We next examined the traction and damage profiles recovered from the standard and stabilizedmethods using delamination tests ( i.e.,
DCB, ENF and FRMM tests). Once again, spurious os-23 , PH a L Figure 17: Geometry and boundary conditions for the fixed ratio mixed mode test. The dimensions are: L = 50 mm,H = 3.1 mm and a = 35 mm. Displacement (mm) L oa d ( N ) Corrected Beam Theory α τ = 10 N/mm α τ = 10 N/mm α τ = 10 N/mm (a) Displacement (mm) L o a d ( N ) ABAQUS Inbuilt n = 10 N/mm n = 10 N/mm n = 10 N/mm (b) Figure 18: Load versus displacement curves for the fixed ratio mixed mode test obtained from the stabilized FEM: (a)[0/0] orientation and (b) [0/90] orientation. cillations were observed in normal and tangential traction fields for a moderately large choice ofthe cohesive stiffness parameter ( i. e., three orders of magnitude larger than the bulk stiffness).However, we found that the load–displacement curves obtained from both standard and stabilizedFEM match well with theory and experiments, implying that the numerical instability observedin the standard method does not corrupt it’s load predictions. The key point is that the numericalinstability seems to only corrupt the traction fields recovered from the standard FEM when thecohesive stiffness parameter is taken to be moderately large. However, simply taking a small cohe-sive stiffness parameter is not appropriate because it diminishes the accuracy of the traction field,especially under mixed-mode loading, where the interface damage occurs before the maximum co-hesive strength is reached. In contrast, our stabilized FEM ensures stability even with large valuesof cohesive stiffness, and enables accurate recovery of the interface traction field. We found thatstandard FEM with structured uniform meshes is less prone to instability in the three delaminationtests. Typically, we only find isolated oscillations in normal and tangential traction profiles nearthe crack tip region, where the traction and damage fields varies sharply. However, the standardFEM with perturbed meshes exhibits spurious oscillations in compressive/contact regions of thedelamination interface, although they are smaller in amplitude.24 a) (b)(c) (d) ↵ n = 10 N/mm ↵ n = 10 N/mm ↵ n = 10 N/mm Figure 19: (a) Normal and tangential traction fields versus interface length from the fixed ratio mixed mode test ([0/0]ply orientation) for different cohesive stiffness values: (a), (b) standard FEM and (c), (d) stabilized FEM.
In conclusion, with the standard FEM the choice of cohesive stiffness can be non-trivial andinvolves a precarious trade-off – prescribing a value that is larger than necessary results in numer-ical instabilities, whereas prescribing a value that is smaller than necessary results in inaccuratetraction recovery. A “correct” choice of cohesive stiffness is not always apparent and can only befound by conducting sensitivity studies based on 1D model estimates. In contrast, the proposedstabilized FEM is robust and stable for a wide range of cohesive stiffness values. Furthermore, asdiscussed in Ghosh et al. [7] the stabilized FEM does not increase the computational cost com-pared to the standard FEM, so it is advantageous for fracture and composite delamination analysisusing cohesive/interface elements. a) (b)(c) (d) ↵ n = 10 N/mm ↵ n = 10 N/mm ↵ n = 10 N/mm Figure 20: Normal and tangential traction fields versus interface length from the fixed ratio mixed mode test ([0/90]ply orientation) for different cohesive stiffness values: (a), (b) standard FEM and (c), (d) stabilized FEM.
Acknowledgements
Gourab Ghosh and Ravindra Duddu gratefully acknowledge the financial support of the Of-fice of Naval Research – award
References [1] Y. Mi, M.A. Crisfield, G.A.O. Davies, and H.B. Hellweg. Progressive delamination usinginterface elements.
Journal of Composite Materials , 32(14):1246–1272, 1998.262] J.A. Pascoe, R.C. Alderliesten, and R. Benedictus. Methods for the prediction of fatigue de-lamination growth in composites and adhesive bonds–a critical review.
Engineering FractureMechanics , 112:72–96, 2013.[3] L. Banks-Sills. 50th anniversary article: review on interface fracture and delamination ofcomposites.
Strain , 50(2):98–110, 2014.[4] J.C.J. Schellekens and R. De Borst. A non-linear finite element approach for the analysis ofmode-I free edge delamination in composites.
International Journal of Solids and Structures ,30(9):1239–1253, 1993.[5] A. Turon, C.G. Davila, P.P. Camanho, and J. Costa. An engineering solution for mesh sizeeffects in the simulation of delamination using cohesive zone models.
Engineering FractureMechanics , 74(10):1665 – 1682, 2007.[6] S. Jim´enez and R. Duddu. On the parametric sensitivity of cohesive zone models for high-cycle fatigue delamination of composites.
International Journal of Solids and Structures ,82:111–124, 2016.[7] G. Ghosh, R. Duddu, and C. Annavarapu. A stabilized finite element method for enforc-ing stiff anisotropic cohesive laws using interface elements.
Computer Methods in AppliedMechanics and Engineering , 348:1013–1038, 2019.[8] G. Alfano and M.A. Crisfield. Finite element interface models for the delamination analysisof laminated composites: mechanical and computational issues.
International Journal forNumerical Methods in Engineering , 50(7):1701–1736, 2001.[9] P. P. Camanho, C. G. Davila, and M.F. De Moura. Numerical simulation of mixed-mode pro-gressive delamination in composite materials.
Journal of Composite Materials , 37(16):1415–1438, 2003.[10] W.G. Jiang, S.R. Hallett, B.G. Green, and M.R. Wisnom. A concise interface constitutive lawfor analysis of delamination and splitting in composite materials and its application to scalednotched tensile specimens.
International Journal for Numerical Methods in Engineering ,69(9):1982–1995, 2007.[11] A. Turon, P.P. Camanho, J. Costa, and J. Renart. Accurate simulation of delamination growthunder mixed-mode loading using cohesive elements: definition of interlaminar strengths andelastic stiffness.
Composite Structures , 92(8):1857–1864, 2010.[12] M. May and S.R. Hallett. A combined model for initiation and propagation of damage underfatigue loading for cohesive interface elements.
Composites Part A: Applied Science andManufacturing , 41(12):1787–1796, 2010.[13] R. Higuchi, T. Okabe, and T. Nagashima. Numerical simulation of progressive damage andfailure in composite laminates using XFEM/CZM coupled approach.
Composites Part A:Applied Science and Manufacturing , 95:197–207, 2017.2714] G.N. Wells and L.J. Sluys. A new method for modelling cohesive cracks using finite elements.
International Journal for Numerical Methods in Engineering , 50(12):2667–2682, 2001.[15] N. Mo¨es and T. Belytschko. Extended finite element method for cohesive crack growth.
Engineering Fracture Mechanics , 69(7):813–833, 2002.[16] J.J.C Remmers, R. de Borst, and A. Needleman. A cohesive segments method for the simu-lation of crack growth.
Computational Mechanics , 31(1-2):69–77, 2003.[17] J. Mergheim, E. Kuhl, and P. Steinmann. A finite element method for the computationalmodelling of cohesive cracks.
International Journal for Numerical Methods in Engineering ,63(2):276–289, 2005.[18] J.H. Song, P.M.A. Areias, and T. Belytschko. A method for dynamic crack and shear bandpropagation with phantom nodes.
International Journal for Numerical Methods in Engineer-ing , 67(6):868–893, 2006.[19] S. Maiti, D. Ghosh, and G. Subhash. A generalized cohesive element technique for arbitrarycrack motion.
Finite Elements in Analysis and Design , 45(8-9):501–510, 2009.[20] X. Zhang, D.R. Brandyberry, and P.H. Geubelle. IGFEM-based shape sensitivity analysis ofthe transverse failure of a composite laminate.
Computational Mechanics , 64(5):1455–1472,2019.[21] J. Zhi, B.Y. Chen, and T.E. Tay. Geometrically nonlinear analysis of matrix cracking and de-lamination in composites with floating node method.
Computational Mechanics , 63(2):201–217, 2019.[22] T. Rabczuk and G. Zi. A meshfree method based on the local partition of unity for cohesivecracks.
Computational Mechanics , 39(6):743–760, 2007.[23] B. Dhas, M. Rahaman, K. Akella, D. Roy, and J.N. Reddy. A phase-field damage model fororthotropic materials and delamination in composites.
Journal of Applied Mechanics , 85(1),2018.[24] F.A. Denli, O. G¨ultekin, G.A. Holzapfel, and H. Dal. A phase-field model for fracture ofunidirectional fiber-reinforced polymer matrix composites.
Computational Mechanics , pages1–18, 2020.[25] X. Lu, M. Ridha, B.Y. Chen, V.B.C. Tan, and T.E. Tay. On cohesive element parameters anddelamination modelling.
Engineering Fracture Mechanics , 206:278 – 296, 2019.[26] Y. Freed and L. Banks-Sills. A new cohesive zone model for mixed mode interface fracturein bimaterials.
Engineering Fracture Mechanics , 75(15):4583–4593, 2008.[27] P.W. Harper and S.R. Hallet. Cohesive zone length in numerical simulations of compositedelamination.
Engineering Fracture Mechanics , 75(16):4774–4792, 2008.2828] A. Mota, J. Knap, and M. Ortiz. Fracture and fragmentation of simplicial finite ele-ment meshes using graphs.
International Journal for Numerical Methods in Engineering ,73(11):1547–1570, 2008.[29] R. Espinha, K. Park, G.H. Paulino, and W. Celes. Scalable parallel dynamic fracture simula-tion using an extrinsic cohesive zone model.
Computer Methods in Applied Mechanics andEngineering , 266:144–161, 2013.[30] R.R. Settgast, P. Fu, S.D.C. Walsh, J.A. White, C. Annavarapu, and F.J. Ryerson. A fullycoupled method for massively parallel simulation of hydraulically driven fractures in 3-dimensions.
International Journal for Numerical and Analytical Methods in Geomechanics ,41(5):627–653, 2017.[31] J.C.J. Schellekens and R. De Borst. On the numerical integration of interface elements.
International Journal for Numerical Methods in Engineering , 36(1):43–66, 1993.[32] M.G.G.V. Elices, G.V. Guinea, J. Gomez, and J. Planas. The cohesive zone model: advan-tages, limitations and challenges.
Engineering Fracture Mechanics , 69(2):137–163, 2002.[33] R. De Borst. Numerical aspects of cohesive-zone models.
Engineering Fracture Mechanics ,70(14):1743–1757, 2003.[34] A. Simone. Partition of unity-based discontinuous elements for interface phenomena: com-putational issues.
International Journal for Numerical Methods in Biomedical Engineering ,20(6):465–478, 2004.[35] Q. Yang and B. Cox. Cohesive models for damage evolution in laminated composites.
Inter-national Journal of Fracture , 133(2):107–137, 2005.[36] T. Q. Thai, T. Rabczuk, and X. Zhuang. Numerical study for cohesive zone model in delam-ination analysis based on higher-order B-spline functions.
Journal of Micromechanics andMolecular Physics , 2(01):1750004, 2017.[37] K. Park, G. H. Paulino, and J. R. Roesler. A unified potential-based cohesive model of mixed-mode fracture.
Journal of the Mechanics and Physics of Solids , 57(6):891–908, 2009.[38] C.G. D´avila, P.P. Camanho, and A. Turon. Effective simulation of delamination in aero-nautical structures using shells and cohesive elements.
Journal of Aircraft , 45(2):663–672,2008.[39] C.G. D´avila, C.A. Rose, and P.P. Camanho. A procedure for superposing linear cohesivelaws to represent multiple damage mechanisms in the fracture of composites.
InternationalJournal of Fracture , 158(2):211–223, 2009.[40] Q.D. Yang, X.J. Fang, J.X. Shi, and J. Lua. An improved cohesive element for shell delamina-tion analyses.
International Journal for Numerical Methods in Engineering , 83(5):611–641,2010. 2941] B.C. Do, W. Liu, Q.D. Yang, and X.Y. Su. Improved cohesive stress integration schemes forcohesive zone elements.
Engineering Fracture Mechanics , 107:14–28, 2013.[42] W. Cui and M.R. Wisnom. A combined stress-based and fracture-mechanics-based model forpredicting delamination in composites.
Composites , 24(6):467–474, 1993.[43] D. Xie and A.M. Waas. Discrete cohesive zone model for mixed-mode fracture using finiteelement analysis.
Engineering Fracture Mechanics , 73(13):1783–1796, 2006.[44] X. Liu, R. Duddu, and H. Waisman. Discrete damage zone model for fracture initiation andpropagation.
Engineering Fracture Mechanics , 92:1–18, 2012.[45] Y. Wang and H. Waisman. Progressive delamination analysis of composite materials usingXFEM and a discrete damage zone model.
Computational Mechanics , 55(1):1–26, 2015.[46] E. Svenning. A weak penalty formulation remedying traction oscillations in interface ele-ments.
Computer Methods in Applied Mechanics and Engineering , 310:460–474, 2016.[47] E. Lorentz. A mixed interface finite element for cohesive zone models.
Computer Methodsin Applied Mechanics and Engineering , 198(2):302–317, 2008.[48] F. Cazes, M. Coret, and A. Combescure. A two-field modified lagrangian formulation forrobust simulations of extrinsic cohesive zone models.
Computational Mechanics , 51(6):865–884, 2013.[49] P. Hansbo. Nitsche’s method for interface problems in computational mechanics.
GAMM-Mitteilungen , 28(2):183–206, 2005.[50] V.P. Nguyen. Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementationcaveats and applications in computational fracture mechanics.
Engineering Fracture Me-chanics , 128:37–68, 2014.[51] T.J. Truster and A. Masud. A discontinuous/continuous Galerkin method for modeling ofinterphase damage in fibrous composite systems.
Computational Mechanics , 52(3):499–514,2013.[52] D. Versino, H. M. Mourad, C. G. Davila, and F. L. Addessio. A thermodynamically consistentdiscontinuous Galerkin formulation for interface separation.
Composite Structures , 133:595– 606, 2015.[53] S.C. Aduloju and T.J. Truster. A variational multiscale discontinuous galerkin formulationfor both implicit and explicit dynamic modeling of interfacial fracture.
Computer Methods inApplied Mechanics and Engineering , 343:602 – 630, 2019.[54] R. Abedi and R. B. Haber. Spacetime simulation of dynamic fracture with crack closureand frictional sliding.
Advanced Modeling and Simulation in Engineering Sciences , 5(1):22,2018. 3055] R. Abedi and P.L. Clarke. A computational approach to model dynamic contact and fracturemode transitions in rock.
Computers and Geotechnics , 109:248–271, 2019.[56] J. Mergheim, E. Kuhl, and P. Steinmann. A hybrid discontinuous Galerkin/interface methodfor the computational modelling of failure.
International Journal for Numerical Methods inBiomedical Engineering , 20(7):511–519, 2004.[57] A. Seagraves and R. Radovitzky. Large-scale 3D modeling of projectile impact damage inbrittle plates.
Journal of the Mechanics and Physics of Solids , 83:48–71, 2015.[58] V.P. Nguyen, C.T. Nguyen, S. Bordas, and A. Heidarpour. Modelling interfacial crackingwith non-matching cohesive interface elements.
Computational Mechanics , 58(5):731–746,2016.[59] H. R. Bayat, S. Rezaei, T. Brepols, and S. Reese. Locking-free interface failure modeling bya cohesive discontinuous galerkin method for matching and non-matching meshes.
Interna-tional Journal for Numerical Methods in Engineering , 121(8):1762–1790, 2019.[60] A. Hansbo and P. Hansbo. A finite element method for the simulation of strong and weak dis-continuities in solid mechanics.
Computer Methods in Applied Mechanics and Engineering ,193(33-35):3523–3540, 2004.[61] J.D. Sanders, J.E. Dolbow, and T.A. Laursen. On methods for stabilizing constraints over en-riched interfaces in elasticity.
International Journal for Numerical Methods in Engineering ,78(9):1009–1036, 2009.[62] F. Liu and R.I. Borja. A contact algorithm for frictional crack propagation with the ex-tended finite element method.
International Journal for Numerical Methods in Engineering ,76(10):1489–1512, 2008.[63] P. Wriggers and G. Zavarise. A formulation for frictionless contact problems using a weakform introduced by Nitsche.
Computational Mechanics , 41(3):407–420, 2008.[64] C. Annavarapu, M. Hautefeuille, and J.E. Dolbow. A Nitsche stabilized finite element methodfor frictional sliding on embedded interfaces. Part II: Intersecting interfaces.
Computer Meth-ods in Applied Mechanics and Engineering , 267:318 – 341, 2013.[65] C. Annavarapu, M. Hautefeuille, and J.E. Dolbow. A Nitsche stabilized finite element methodfor frictional sliding on embedded interfaces. Part I: single interface.
Computer Methods inApplied Mechanics and Engineering , 268:417–436, 2014.[66] C. Annavarapu, R.R. Settgast, S.M. Johnson, P. Fu, and E.B. Herbold. A weighted Nitschestabilized method for small-sliding contact on frictional surfaces.
Computer Methods in Ap-plied Mechanics and Engineering , 283:763–781, 2015.[67] M. Juntunen and R. Stenberg. Nitsche’s method for general boundary conditions.
Mathemat-ics of Computation , 78(267):1353–1374, 2009.3168]
Abaqus Analysis User’s Guide, Version 6.14 . Dassault Systemes Simulia Corp., 2014.[69] C. Annavarapu, M. Hautefeuille, and J.E. Dolbow. Stable imposition of stiff constraints inexplicit dynamics for embedded finite element methods.
International Journal for NumericalMethods in Engineering , 92(2):206–228, 2012.[70] D.N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discontinuousGalerkin methods for elliptic problems.
SIAM Journal on Numerical Analysis , 39(5):1749–1779, 2002.[71] R. Liu, M.F. Wheeler, and C.N. Dawson. A three-dimensional nodal-based implementation ofa family of discontinuous Galerkin methods for elasticity problems.
Computers & Structures ,87(3):141–150, 2009.[72] C. Annavarapu, M. Hautefeuille, and J.E. Dolbow. A robust Nitsche’s formulation for in-terface problems.
Computer Methods in Applied Mechanics and Engineering , 225:44–54,2012.[73] J.D. Sanders, T.A. Laursen, and M.A. Puso. A Nitsche embedded mesh method.
Computa-tional Mechanics , 49(2):243–257, 2012.[74] A. Embar, J. Dolbow, and I. Harari. Imposing Dirichlet boundary conditions with Nitsche’smethod and spline-based finite elements.
International Journal for Numerical Methods inEngineering , 83(7):877–898, 2010.[75] S. Jim´enez, X. Liu, R. Duddu, and H. Waisman. A discrete damage zone model for mixed-mode delamination of composites under high-cycle fatigue.
International Journal of Frac-ture , 190(1-2):53–74, 2014.[76] J.N. Reddy and A. Miravete.
Practical analysis of composite laminates . CRC press, 2018.[77] L.E. Asp, A. Sj¨ogren, and E.S. Greenhalgh. Delamination growth and thresholds in a car-bon/epoxy composites under fatigue loading.
Journal of Composites Technology & Research ,23(2):55–68, 2001.[78] S.H. Song, G.H. Paulino, and W.G. Buttlar. A bilinear cohesive zone model tailored forfracture of asphalt concrete considering viscoelastic bulk material.
Engineering FractureMechanics , 73(18):2829–2848, 2006.[79] S. R. Hallett, B. G. Green, W.G. Jiang, and M.R. Wisnom. An experimental and numeri-cal investigation into the damage mechanisms in notched composites.