Application of continuum decohesive finite element to progressive failure analysis of composite materials

Composite Structures 212 (2019) 365–380

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Application of continuum decohesive finite element to progressive failure analysis of composite materials Shiyao Lina, Nhung Nguyenb, Anthony M. Waasa, a b

T

⁎

Composite Structures Laboratory, Aerospace Engineering, University of Michigan, Ann Arbor, MI 48109-2140, United States Department of Industrial and Systems Engineering, International University – Vietnam National University, Ho Chi Minh City, Viet Nam

A R T I C LE I N FO

A B S T R A C T

Keywords: Progressive failure analysis Representative volume element Delamination toughness test Mixed mode law

The continuum decohesive finite element (CDFE) is a novel finite element technique combining continuum and cohesive crack modeling seamlessly. In CDFE, the transition from a continuum to non-continuum is modeled physically by introducing pairs of dummy nodes to account for the crack separations. A static condensation algorithm is applied to solve for and preserve the separations. In this paper, CDFE has firstly been applied to the modeling of transverse crack growth in a representative volume element (RVE) of composite materials. Microscopic cracks initiate separately and coalesce into a macroscopic transverse crack. In this case, mode I cracking is dominant and the maximum principal stress criterion is used for crack angles. The second problem set is the delamination toughness tests, including a double cantilever beam (DCB) test, an end notched flexure (ENF) test and mixed mode bending (MMB) tests of three mixed mode ratios. In these cases, the crack plane is welldefined and mode mixity is encountered. For better numerical stability, the CDFE inner-element discretization scheme is modified and a novel mixed mode cohesive formulation is implemented. Through analyses, the capability of CDFE to deal with general progressive failure problems of composite materials is shown.

1. Introduction Fiber reinforced polymer composites (FRPC) have shown appealing applications in industry, especially in the aerospace-related fields because of desirable properties such as high specific strength and stiffness, thermal stability, conductivity, capability to be specifically designed, etc. However due to the complexity in material composition and damage mechanisms, virtual progressive failure analysis (PFA) with predictive capability is essential in the assessment of the structural integrity and damage tolerance (SIDT) of structures made of FRPC. Various PFA models have been developed with different levels of resolution. In the macroscale failure analyses, a composite laminate is treated as a general orthotropic material with apparent mechanical properties. The crack angle and path are not presumed. In the mesoscale failure analyses, a composite laminate is usually treated as plies of transversely isotropic lamina adhered by resin-rich interfacial layers. Matrix crack and fiber rupture can take place as intra-ply failure while delamination can occur as inter-ply failure. In this scale, for the intraply modeling, crack angle is usually assumed parallel or perpendicular to the fiber orientation while crack path is unknown. For the inter-ply crack (delamination) modeling, both the crack angle and crack path are fixed. In the microscale failure analyses, where matrix and fibers are ⁎

explicitly modeled, a crack can initiate and propagate within the matrix with arbitrary crack angle and location. In this case, the crack angle needs to be determined by an initiation criterion. From a literature survey, the most used criteria in the scope of damage mechanics for arbitrary crack growth in composites include maximum principal stress criterion (also referred to as the Rankin criterion, maximum circumferential stress criterion or maximum hoop tensile stress criterion) [1–4], maximum principal strain criterion [5] and maximum principal shear stress criterion (also referred to as the Mohr–Coulomb criterion) [6,7]. The scenario to which these criteria are applied vary. If the direction of principal normal stress is positive (tensile), maximum principal stress or strain criteria should be used. The maximum principal stress and strain criteria are essentially identical, but differ in the way of implementation. For cases involved with pre-peak nonlinearity [4], it is easier to formulate the initiation criterion based upon the maximum principal strain. If the direction of the principal normal stress is negative (compressive), it is assumed that the cracks within the element is constrained from growing in mode I. Instead, a mode II crack would form and orient along the plane of maximum shear stress. In this case, the Mohr–Coulomb criterion should be used. Numerical modeling techniques for PFA of composite materials rely mostly on the finite element method (FEM). Non-FEM methods such as

Corresponding author. E-mail address: [email protected] (A.M. Waas).

https://doi.org/10.1016/j.compstruct.2019.01.021 Received 22 October 2018; Received in revised form 7 December 2018; Accepted 2 January 2019 Available online 09 January 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 212 (2019) 365–380

S. Lin et al.

enriched laminate element for tensile failure of composites. Based on the FNM idea, an enriched shell element has been developed and applied to the study of low velocity impact damage of composite materials and the analysis time has been reduced compared to traditional, 3D impact modeling which uses layer by layer discretization using FEM. However, finer features of the delamination shapes which can influence the compression after impact strength cannot be captured with such an approach, [36]. ECDM is developed in the framework of XFEM. However, unlike XFEM, enriched degrees of freedom (DOFs) are condensed after crack initiation and the crack opening follows the cohesive behavior without the crack-tip enrichment in XFEM for singular crack-tip stress distribution. After condensation, an equivalent stiffness matrix is obtained such that the method can be implemented with standard FEM codes [32]. Various problems of composite materials have been studied with ECDM including four point bending of a stiffened laminated composite panel, crack propagation in a composite T-joint and mode I crack opening and delamination migration in composite beams [32–35]. Through the studied examples, the effectiveness, robustness and efficiency of the method have been shown. Additionally, it is worthwhile to point out that, in the delamination migration example, compared to the traditional cohesive zone model (CZM) approach, ECDM reduce the CPU time of prediction by more than 90%. Inspired by the VMCM, where a sharp crack is represented with enriched elemental shape functions and static condensation is used to remove the extra DOFs, continuum decohesive finite element (CDFE) has been proposed as a novel FEM based technique to seamlessly bridge the continuum damage model and the discrete damage model [37,38]. In CDFE, a crack is modeled as a physical discontinuity inside an element by introduction of pairs of dummy nodes. Before the crack initiation, the CDFE elemental formulation follows the continuum element scheme, while after the initiation, a cohesive crack following the discrete cohesive zone model (DCZM) is inserted inside the element with the calculated or presumed crack angle through the element centroid. Within each time increment, the equivalent elemental stiffness matrix is obtained by statically condensing the stiffness matrix of the DCZM element representing the cohesive crack into the continuum element, which is inspired by the last step of VMCM. Detailed mathematical and numerical illustration of CDFE will be given in Section 2. After a careful literature study, it should be pointed out that the proposed CDFE method has similarities to the method of ECDM. In both methods, introduced DOFs are condensed into equivalent stiffness matrix such that the methods can be implemented into standard FEM framework and to improve the efficiency of the methods. And, the crack initiation and propagation in both methods are based on cohesive crack growth. The major difference between the two methods are stated as follows. On a methodological level, CDFE is inspired by VMCM while the ECDM is motivated from XFEM. In CDFE, cohesive crack is physically introduced into the element while in ECDM, like XFEM, the crack is represented by enriched DOFs without being physically inserted into the element. In this way, the partition of unity (POU) in CDFE is for cohesive crack insertion while that in ECDM is solely for numerical integration in sub-domains. Also, due to this difference, the damage factor in CDFE is based upon physical crack separation, through the traction-separation law, while that in ECDM, the cohesive law is related to the strain field. On a numerical implementation level, the current CDFE is developed within an explicit framework, using Abaqus subroutine VUEL. The ECDM is developed within an implicit framework, using Abaqus subroutine UEL. However, the transition between implicit and explicit is always achievable. The cohesive law dictates the traction-separation behavior of the cohesive crack. Cohesive traction-separation laws were implemented in node-based user elements and proposed as DCZM by Xie and Waas [39]. Following this work, a thorough parametric study about the effects of the cohesive law on DCZM modeling has been performed by Gustafson and Waas [40]. In these studies, the crack tip is defined through the simultaneous vanishing of normal and tangential traction components.

peridynamics [8–10], smoothed particle hydrodynamics (SPH) [11,12], element-free Galerkin method (EFG) [13] are promising in modeling dynamic fracture and massive crack interaction of composite materials due to incorporating pointwise interaction among material points. However, in the current state of the art, these methods are computationally expensive to be applied to solving large scale engineering problems of interest. FEM based models usually fall under two categories, continuum damage governed or discrete damage governed. In the continuum damage model, the crack separation is smeared over a characteristic length as crack strain and the crack softening behavior is defined through an equivalent stress–strain relationship based upon a tractionseparation law that encompasses the dissipated energy during failure. The representative studies of the continuum damage model include the crack band method [14] and the smeared crack method [15]. The traditional crack band method has been enhanced with Schapery Theory [16], which is a thermodynamically-based work potential theory, to capture the pre-peak nonlinearity caused by micro-damage in matrix [4]. Advanced features such as interaction between in-plane damage and delamination and crack spacing are implemented as the intra-inter crack band model (I2CBM) [17–19]. In the discrete damage model, when the crack angle and crack path are known a priori, such as delamination and adhesively bonded structure modeling, cohesive cracks are modeled using interfacial elements [20,21]. In the cases where crack path is not presumed, more powerful discrete damage methods are developed. Nodal enrichment or elemental enrichment is performed on the finite element and advanced FEM technologies are proposed, such as the eXtended FEM (XFEM) [22,23], phantom node method (PNM) [24], strong discontinuity approach (SDA) [2] and the variational multiscale cohesive method (VMCM) [25,3]. In nodal enrichment methods [22–24], the support for discontinuous displacement modes is one of the given nodal shape functions and in elemental enrichment methods [2,25,3], the support is one of the given elements. A comprehensive comparative study of different discrete damage models has been reported in [26]. It was concluded that, when implemented based on the same mesh, both methods converge to coincident results, while in the multiple-crack cases, the computational cost is more favorable to the elemental enrichment method. In fact, in the examples studied, the elemental enrichment was found to be independent of the number of cracks while the nodal enrichment cost was found to increase linearly with the number of cracks. This observation suggests that elemental enrichment is more suited for problems where damage can include the need to capture several hundreds of cracks. Recently, methods including the augmented FEM (AFEM) [27,1,28,29], floating node method (FNM) [30,31] and extended cohesive damage model (ECDM) [32–35] have emerged as promising techniques for PFA of composite materials. A-FEM firstly inherited the scheme of PNM, in which the crack element has two states, namely the weak discontinuity state and the strong discontinuity state. The elemental integrals are performed on parts of the coincident elements [27]. Then, with the implementation of a novel consistency-check based condensation algorithm, the dummy nodes are introduced to deal with the crack separation and the crack inside the element is physically modeled [1,28,29]. So far, most of the A-FEM studies concentrate on arbitrary crack growth analyses, in both composite materials and general solid materials. The FNM, before crack initiation, is also essentially equal to the PNM. After the initiation, the scheme becomes different, in which the FNM element is split into sub-elements, and the isoparametric mapping and elemental integration are performed separately for the sub-elements [30]. A significant feature of FNM is that the intra-ply damage and inter-ply damage are captured in the same element, such that direct interaction between different types of damage can be modeled. The interaction between different damage modes are of importance in PFA of composite materials. An enriched laminate element based on FNM is developed within which all the plies and interfaces are modeled [31]. High-fidelity results have been acquired with the 366

Composite Structures 212 (2019) 365–380

S. Lin et al.

Following this definition, Nguyen and Waas reexamined different cohesive formulations and proposed a novel mixed mode cohesive formulation [41]. The formulation, in the 2D case, as before, only requires four material parameters, including fracture strengths and toughness of mode I and mode II. In the method, the mode II traction-separation law is scaled to the mode I law, and an effective separation is defined for the scaled laws. Effective tractions are determined on the scaled laws based on the effective separation and finally scaled back to physical values. The traction components of both modes are guaranteed to vanish simultaneously and smoothly by the architecture of the method. For lamina-level analyses, where mode mixity is usually encountered, this characteristic of the mixed mode cohesive formulation would be beneficial to obtain more physically sound and numerically stable results. Therefore, for mesoscale composite damage analysis, CDFE is enhanced with this formulation and illustrated in Section 2. Prior work has been done using CDFE to analyze open hole tension problems, single edge notched bending problems [37] and tensile failure problems of composite materials [38]. Mesh-objectivity has been demonstrated through examples with different mesh densities. However, the capability of CDFE to model arbitrary crack growth has not been demonstrated yet. In this work, the first example applies CDFE to modeling the transverse crack formulation in a RVE of composite material. In this case, mode I cracking is dominant and the maximum principal stress criterion is utilized. There are two major reasons for using this criterion. The first reason is the RVE is cracked under transverse tensile stress, which prevent the usage of the Mohr–Coulomb criterion. The second reason is that the CDFE results are compared with that in [6], where the maximum principal stress criterion was implemented to study the transverse tensile crack of the RVE. Three meshes of the RVE with the bi-linear and exponential traction-separation laws are studied to investigate the sensitivity of the method to mesh sizes and cohesive law types. The second example set deals with problems where mode mixity is encountered. To better capture the mixed mode scenarios, CDFE is enhanced with a new inner-element discretization scheme and the mixed mode cohesive formulation. This problem set is composed of three MMB tests of different mixed mode ratios. Nevertheless, a DCB and an ENF test are also modeled for comparison. Results and a discussion can be found in Section 3.

Fig. 1. CDFE mathematical scheme: (a) before initiation, (b) after initiation.

is as shown in Fig. 1(a), σC is the critical strength for mode I cracking. The crack angle is fixed as being perpendicular to the maximum principal stress direction at initiation.

The PVW of the CDFE element at and after crack initiation is expressed as, 3.

∫Ω

1

∫Ω

2

∫Ω w·bdV + ∫Ω w·bdV + ∫Γ + ∫ w·fdS + ∫ w·t (〈w〉)dS Γ Γ + ∫ w·t (〈w〉)dS Γ

∇w: σ dV =

1

2

w·fdV

t1

c1

c2

(3)

In Eq. (3), Ω1 and Ω2 stand for two sub-domains. Γt1, Γt 2, Γk1 and Γk2 are normal and kinematic boundaries on the corresponding sub-domains. Γc1 and Γt2 represent two crack faces, on which crack traction t is distributed. t is determined according to designated traction-separation law and the crack separation 〈w〉, where 〈·〉 represents the determination of separation from the displacement field. The calculation of t will be explained in detail in Section 2.1.2. 2.1.2. Damage evolution and criteria In this section, the damage evolution scheme together with initiation and failure criteria will be explained in detail.

2.1. CDFE for arbitrary crack analysis

2.1.2.1. Damage evolution. The damage evolution of the crack is dictated by the cohesive laws. Any type of cohesive law can be implemented with CDFE. In this work, the bilinear cohesive laws are adopted. The mode I and II traction-separation laws are shown in Fig. 2. To determine these laws, five input parameters are required, including mode I crack strength σC and fracture toughness GIC , mode II crack strength τC and fracture toughness GIIC and the ratio between the initiation separation δn0, δt 0 and the critical separation δnC , δtC , which is denoted as η. The critical normal and tangential separations are

2.1.1. Mathematical formulation In this section, the scheme of the CDFE method for arbitrary crack problems will be illustrated and the principle of virtual work (PVW) formulation for the CDFE before and after initiation will be given. 2.1.1.1. Before initiation. Before crack initiation, the CDFE element is identical to a continuum element. The PVW then provides, t

∇w: σ dV +

t2

2. Methodology

∫Ω ∇w: σ dV = ∫Ω ∇w·bdV + ∫Γ w·fdS

(2)

σn ⩾ σC

(1)

where w is the virtual displacement vector, σ is the actual stress tensor, b is the body force vector, f is the external surface force vector, Γt is the normal boundary, Γk is the kinematic boundary, and Ω is the continuum domain. An illustrating diagram is shown in Fig. 1(a). 2.1.1.2. At and after initiation. As demonstrated in Fig. 1(b), right at and after the initiation of the crack, a cohesive crack following the DCZM model is embedded into the continuum body through the centroid with calculated crack angle. The CDFE can be viewed as two continuum sub-elements connected by a DCZM element. In this work, for arbitrary crack problems, crack initiation criterion is the maximum principal stress criterion, as in Eq. (2). In the equation, σn is the normal component of the principal stress, the local principal coordinate system

Fig. 2. Cohesive laws (a) mode I and, (b) mode II. 367

Composite Structures 212 (2019) 365–380

S. Lin et al.

calculated as in Eq. (4). The initiation separations, δn0 and δt0 , are determined as in Eq. (5).

δnC =

2GIC , σC

(4a)

δtC =

2GIIC τC

(4b)

δn0 = ηδnC ,

(5a)

δt 0 = ηδtC

(5b)

The traction laws for mode I and mode II are,

tn (δn ) =

σC δn (δn ⩽ δn0), δn0

(6a)

tn (δn ) =

δn − δn0 σC (δn0 ⩽ δn < δnC ), δnC − δn0

(6b)

tn (δn ) = 0 (δnC ⩽ δn )

(6c)

τ tt (δt ) = C δt (δt ⩽ δt 0), δt 0

(7a)

tt (δt ) =

δt − δt 0 τC (δt 0 ⩽ δt < δtC ), δtC − δt 0

Fig. 3. CDFE FEM scheme for arbitrary crack problems: (a) before initiation, (b) after initiation.

relating displacement and strain. D is the material stiffness matrix, N is the matrix of interpolation functions and Fn is the nodal force vector.

(7b) (7c)

tt (δt ) = 0 (δtC ⩽ δt )

2.1.3.2. At and after initiation. After crack initiation, the continuum element is split into two sub-elements connected by a DCZM element. The scheme is shown in Fig. 3(b). The element stiffness matrices for the two sub-elements follow the standard FEM scheme, where i denotes sub-element domain 1 and 2.

where, δn, δt and tn, tt are normal and tangential components of separation and traction vectors. It should be noted that, in CDFE for arbitrary crack problems, only mode I cracking is considered. Mode II damage evolution law is demonstrated here for the completeness and the structure of the article and for later examples.

Ki =

2.1.2.2. Initiation criterion. For CDFE method dealing with arbitrary crack problems, the crack initiation criterion is the maximum principal stress criterion, as in Eq. (2). After the initiation, the crack angle is fixed, being perpendicular to the maximum principal tensile stress.

(8)

δn ⩾ δnC

(9)

R

F ext e =

(10b)

∫Ω N T·bdV + ∫Γ N T·fdS + ∑ N T·Fn t

(11)

(12)

As in Eq. (12), R is the rotation matrix transforming the stiffness matrix in the local coordinates of the principal direction to the global coordinates. R is written in Eq. (13) in detail, in which θ is the crack angle. K abα means the α component of the stiffness of the spring connecting node a and b. For example, K57t represents the stiffness in the local tangential direction of the spring 5–7. The t and n components of the stiffness are local tangential and normal secant stiffness on the corresponding traction-separation laws, which are calculated in Eq. (14), where tt and tn are determined from Eqs. (6) and (7). It should be noted that, since arbitrary crack problem is mode I dominated, the tangential traction is automatically zero, hence the secant stiffness K abt is also zero. However, the Eqs. (12)–(14) are still written in full form for completeness.

(10a)

∫Ω BTDBdV,

(i = 1, 2)

− K57t 0 0 0 0 0 0 ⎤ ⎡ K57t K57n − K57n 0 0 0 0 0 ⎥ ⎢ 0 ⎢ 0 K 68t − K 68t 0 0 0 0 0 ⎥ ⎥ ⎢ K K − 0 0 0 0 0 0 68n 68n ⎥ = RT ⎢ ⎢− K57t K57t 0 0 0 0 0 0 ⎥ ⎢ 0 K K − 0 0 0 0 0 ⎥ n n 57 57 ⎥ ⎢ K K − 0 0 0 0 0 0 ⎥ ⎢ 68t 68t ⎢ 0 K 68n ⎥ − K 68n 0 0 0 0 0 ⎦ ⎣

2.1.3.1. Before initiation. Before crack initiation, the CDFE element follows the classical finite element formulation. A quadrilateral element is taken as an example in Fig. 3(a). The elemental equations are,

Ke =

B Ti DB i dV

KDCZM

2.1.3. FEM formulation In this section, the FEM scheme of the CDFE method will be explained in detail.

RHS = K e Ue − F ext e ,

i

As shown in Fig. 3(b), the DCZM element is a node-based cohesive element. Two discrete nonlinear elements connect the upper and lower surfaces of the element. The stiffness matrix for the DCZM element 5 − 6 − 7 − 8 is in Eq. (12).

2.1.2.3. Failure criterion. Since only mode I cracking is involved in the arbitrary crack problems, the failure criterion is simply the mode I energy dissipated GI being greater than mode I critical energy GIC , as in Eq. (8). Since no tangential traction is involved, this is essentially equivalent to normal separation δn being greater than the normal critical separation δnC , which is shown in Eq. (9). Eq. (9) is implemented in the CDFE method.

GI ⩾ GIC

∫Ω

(10c)

where K e is the elemental stiffness matrix, Ue is the nodal displacement vector, F ext is the elemental external force vector, B is the matrix e 368

Composite Structures 212 (2019) 365–380

S. Lin et al.

0 0 0 0 0 0 ⎤ ⎡ cosθ sinθ 0 0 0 0 0 0 ⎥ ⎢− sinθ cosθ ⎢ 0 0 cosθ sinθ 0 0 0 0 ⎥ ⎢ 0 0 − sinθ cosθ 0 0 0 0 ⎥ R= ⎢ ⎥ θ θ 0 0 0 0 cos sin 0 0 ⎥ ⎢ ⎢ 0 0 0 0 − sinθ cosθ 0 0 ⎥ ⎢ 0 0 0 0 0 0 cosθ sinθ ⎥ ⎥ ⎢ 0 0 0 0 0 − sinθ cosθ ⎦ ⎣ 0

(n) , refers to the time increment n. Otherwise, the variables are temporary and would not be passed back to the FEM solver or saved. As shown in Fig. 4, at the start of time increment n, the information inX, cluding the nodal coordinates material properties E11, E22, ν12, G12, σC , τC , GIC and GIIC , current nodal displacement U (n) , values from the previous increment including crack status ’cs (n − 1) ’, crack angle θ (n − 1) , displacement vector of the dummy nodes U∧ (n − 1) and normal crack separation are passed into the CDFE program. The variable ′cs′ means the crack status. If ′cs′ is equal to 0, the element is still pristine. If ′cs′ is equal to 1, the crack has initiated in the element. If ′cs′ is equal to 2, the element has failed. At the start of each increment, ’cs (n − 1) ’ is first checked to see which state the element is at. At the end of the increment, the right-hand side, updated crack status ’cs (n) ’, crack angle θ (n) , displacement vector of the dummy nodes U∧ (n) and normal separation δn(n) are output and passed back to the FEM solver.

(13)

K abt =

tt (a = 5, 6;b = 7, 8), δt

(14a)

K abn =

tn (a = 5, 6;b = 7, 8) δn

(14b)

In Eqs. (14), the normal and tangential separations are obtained from the displacement vector U∧ of the dummy nodes 5-6-7-8. This procedure is shown in Eq. (15).

δn = L n RU∧

(15a)

δt = L t RU∧

(15b)

2.2. Enhanced CDFE for composite damage analysis 2.2.1. Enhanced CDFE Unlike arbitrary crack growth problems, damage in composite laminates tend to occur with well-defined crack angles. For intra-ply damage, matrix splitting is usually along the fiber and the fiber rupture is perpendicular to the fiber direction. For inter-ply damage, the delamination is most likely to happen between plies. Therefore, CDFE for composite damage analysis does not require principal stress calculation and can be specially treated for mixed mode scenarios. In this section, CDFE is enhanced with a new inner-element discretization scheme and a novel mixed mode cohesive formulation proposed by Nguyen and Waas [41].

where δn and δt are normal and tangential separation vector, L t and L n are the matrices correlating separations with displacements in local coordinates. L t and L n are shown in Eqs. (16).

− 1 0 0 0 1 0 0 0⎤ Ln = ⎡ ⎣ 0 0 − 1 0 0 0 1 0⎦

(16a)

0 − 1 0 0 0 1 0 0⎤ Lt = ⎡ ⎣0 0 0 − 1 0 0 0 1⎦

(16b)

The next step is to assemble K1, K2 and KDCZM into the global elemental stiffness matrix K , and the nodal displacement and force vectors are, ∗∗ ∗∧ U∗ F∗ K = ⎡ K∧∗ K∧∧ ⎤, U = ⎡ ∧ ⎤, F = ⎡ ∧ ⎤ ⎣F ⎦ ⎣U ⎦ ⎣K K ⎦

2.2.1.1. The new inner-element discretization scheme. At the time increment where initiation occurs, the inner-element discretization scheme would change instantaneously from one continuum element to two continuum sub-elements connected by a DCZM element, namely from Fig. 3(a)–(b). The sudden transition from one element to multiple sub-elements would generate an abrupt stiffness degradation due to the addition of extra DOFs. Such a stiffness degradation would induce an instantaneous drop of the load on the element. After the load drop, the load would continue to pick up following the new elemental stiffness until it degrades according to the traction-separation law. A sudden local stiffness change is also mentioned in [28]. To avoid such issue, since the crack angle is well defined, in enhanced CDFE, the DCZM element is fixed along the fiber direction through elemental centroid from the start of the analysis. In this way, when the crack initiates, the DCZM would change from the initial region to the softening region, while the inner-element discretization is constant throughout the analysis, hence no extra DOF is added and the abrupt change of the stiffness is avoided. The new inner-element scheme is shown in Fig. 5(a) and (b).

(17)

where, ∗ stands for DOFs relative to nodes of the CDFE element before initiation, ∧ stands for the DOFs associated to the dummy nodes on the cohesive crack. To be specific, as illustrated in Fig. 3(b), ∗ represents DOFs relative to node 1, 2, 3, 4 and ∧ represents DOFs relative to node 5, 6, 7, 8. It should be noticed that, since dummy nodes 5, 6, 7 and 8 are not explicit to the FEM solver, no external nodal forces are exerted on them, and hence F∧ = {0} . The final step of the CDFE scheme is to perform static condensation to solve for the inner-element equilibrium and condense the nodal displacements associated with the dummy nodes to generate an equivalent stiffness matrix and the right-hand side RHS . The static condensation is as in Eqs. (18).

U∧ = −(K∧∧)−1K∧∗U∗,

(18a)

∗∗ ∗∧ ∧∧ −1 ∧∗ K eq e = K − K (K ) K ,

(18b)

ext ∗ RHS = K eq e U − Fe

(18c)

2.2.1.2. The novel mixed mode cohesive formulation. When mode mixity is encountered, a popular failure criterion is the first-order power law, as in Eq. (19). In the equation, GI and GII are dissipated energy of mode I and II.

After the static condensation step, the equivalent stiffness matrix K eq e and RHS can be readily implemented into any FEM solver with both implicit and explicit time integration algorithms. For the robustness of the analysis, explicit solver of Abaqus is utilized and CDFE is encoded through the user element subroutine VUEL. The input to the subroutine is the material properties, nodal coordinates, nodal displacement vector of the current time increment and the state variables inherited from the previous increment. The output of the subroutine to the FEM solver is the RHS vector and updated state variables. The details of the implementation can be found in Section 2.1.4.

GI G + II ⩾ 1 GIC GIIC

(19)

However, in a mixed mode scenario, at propagation, there are still finite values for both mode I and II traction components. When the tractions are forced to zero, instantaneous changes of the stiffness matrix and RHS would be generated. This instantaneous change would possibly cause spurious numerical oscillations in the analysis, which reduces the robustness of the method in mixed mode cases. In addition, propagation strictly requires the simultaneous vanishing of tractions due to the creation of new surfaces. Nguyen and Waas proposed a novel mixed mode cohesive scheme

2.1.4. Flow chart Flow chart of CDFE for arbitrary crack problems is shown in Fig. 4 to elucidate the method more clearly. The variables with a superscript 369

Composite Structures 212 (2019) 365–380

S. Lin et al.

INPUT: X, U(n), E11, E22, Ȟ12, G12, ıC, IJC, GIC, GIIC, cs(n-1), ș(n-1), U^(n-1), įn(n-1)

2GIC /

nC

IF

C

cs(n-1)==0

,

2GIIC /

tC

C

THEN T

K en

B e DB e dV

ȍe

(n)

(n)

(n) Dİ (n) İe BeUe , ı e e acquire principal direction

RHS (n) e

U ^(n)

U ^(n

IF

n

1)

0, į (n) n

n

į (n n

1)

0

THEN

C (n)

0,

(n)

cs (n)

1,

(n)

cs

and normal pincipal stress

p

(n) K (n) Feext(n) e Ue

( n -1)

0

ELSE p

ENDIF ELSEIF cs(n-1)==1 THEN T

K1

ȍ1

B1 DB1dV , K 2

T

ȍ2

B 2 DB 2 dV

Determine K 57 n , K 57 t , K 68 n , K 68t from equation (14) Determine K DCZM from equation (12) Assemble into global stiffness matrix K n Static condensation : K eq e

RHS (n) e IF į n

K *^

K ^*

K ^^

K ** K *^ (K ^^ ) 1 K ^*

K eq(n) U (n) Feext(n) , U^(n) e e

(K ^^ ) 1 K ^* U (n) , į (n) n

L n RU^(n)

į nC THEN (n)

1,

(n)

( n -1)

cs (n)

2,

(n)

( n -1)

cs

K **

ELSE ENDIF ELSEIF cs(n-1)==2 THEN

K en RHS (n) e

T

ȍe

B e DB e dV /10 20 (n) K (n) Feext(n) e Ue

ENDIF (n) (n) ^(n), į (n) OUTPUT: RHS (n) n e , cs , ș , U Fig. 4. Flow chart of the CDFE for arbitrary crack problems.

will be explained in detail.

[41], where no extra material parameters are required while the tractions are guaranteed to vanish simultaneously and smoothly at propagation. In this work, this novel mixed mode cohesive formulation is implemented in CDFE to enhance its behavior in mixed mode cases. Details of the mixed mode cohesive formulation will be given in Section 2.2.3.

2.2.3.1. Damage evolution. This part will concentrate on introducing the novel mixed mode cohesive formulation proposed by Nguyen and Waas [41]. This formulation requires no extra material parameters other than σC , GIC , τC , GIIC and η, while the tractions are guaranteed to vanish simultaneously and smoothly at propagation The scheme is shown in Fig. 6 and the formulation is shown in Eq. (20)–(24). In the formulation, first, a scale factor β is calculated by stretching the mode II law to mode I, while keeping the stretched mode II critical separation equal to that of mode I. The stretching is represented in Eq. (20) and the critical strength and separation for the stretched mode II laws are βτC and βτnC (=δnC ). With scaled traction-separation law, an effective separation δ , as determined by Eq. (21), is used to obtain the effective tractions t1 and t2 from the scaled traction-separation laws. The calculation of the effective tractions is shown in Eqs. (22) and (23).

2.2.2. Mathematical formulation In enhanced CDFE for composite damage analyses, the DCZM element is fixed inside the element from the start. Therefore, the PVW framework for the enhanced CDFE before and after initiation is identical to Eq. (3). 2.2.3. Damage evolution In this section, the damage evolution together with initiation and failure criteria following the novel mixed mode cohesive formulation 370

Composite Structures 212 (2019) 365–380

S. Lin et al.

4

4

3

3 Continuum sub-element 1

Continuum sub-element 1

6(8) 5(7)

5

Continuum sub-element 2

1

2

(a)

2

(b)

(24a)

tt (δn, δt ) = t2 ×

δt δ

(24b)

2.2.3.2. Initiation criterion. The initiation criterion for the mixed mode cohesive formulation is simply the effective separation being greater than the initiation separation of mode I or the stretched mode II. The equivalence of them is proved in Eq. (25) and the initiation criterion is in Eq. (26).

Continuum sub-element 2

1

δn , δ

8 6

7

tn (δn, δt ) = t1 ×

βδt 0 = βηδtC = ηδnC = δn0

(25)

δ ⩾ δn0 (or δ ⩾ βδt 0)

(26)

Fig. 5. Modified CDFE for composite damage analysis: (a) before initiation, (b) at and after initiation.

2.2.3.3. Failure criterion. The failure criterion for the mixed mode cohesive formulation is simply the effective separation being greater than the critical separation of mode I or the stretched mode II. These two values are equal and hence the simultaneous vanishing of tractions of both modes are guaranteed. The failure criterion is shown in Eq. (27).

tt

tn ı0

IJ0

δ ⩾ δnC (or δ ⩾ βδtC )

įnC t1

įn0

įtC įt0

t2

(a)

ı0

2.2.4. FEM formulation In this section, the FEM scheme of the enhanced CDFE for composite damage analyses will be explained in detail.

(b)

2.2.4.1. Before initiation. Before crack initiation, the DCZM element is already fixed between the two continuum sub-elements along the fiber angle. The stiffness matrices are calculated as in Eqs. (11) and (12). The secant stiffness used by KDCZM is obtained from Eq. (14). The tractions are determined from Eqs. (22a), (23a) and (24). As before, the element equations are in Eq. (17) and the static condensation step is in Eq. (18).

IJ0

t1

0

t2

įn0

įnC*

į

(c)

t0

į

(d)

tC

įtC

2.2.4.2. At and after initiation. Since in the enhanced CDFE, there is no change of the inner-element scheme throughout the analysis, the element formulation at and after the initiation is identical to that before initiation. The only difference is the traction calculation. Before initiation, the calculation of the tractions follow Eq. (22a), (23a) and (24). At and after initiation, the calculation follows (22b), (22c), (23b), (23c) and (24).

Fig. 6. (a) Mode I cohesive law (b) mode II cohesive law (c) scaled mode I cohesive law (d) scaled mode II cohesive law.

The last step is to scale back the effective tractions to physical values, as shown in Eqs. (24).

β=

δ=

δnC δtC

2.2.5. Flow chart Flow chart of CDFE for arbitrary crack problems is shown in Fig. 7 to elucidate the method more clearly. The variables with a superscript (n) means at the time increment n. Otherwise, the variables are temporary and would not be passed back to the FEM solver or saved. As shown in Fig. 7, at the start of time increment n, the information inX, cluding the nodal coordinates material properties E11, E22, ν12, G12, σC , τC , GIC and GIIC , current nodal displacement U (n) , values from the previous increment including crack status cs (n − 1) , crack angle θ , displacement vector of the dummy nodes U∧ (n − 1) and normal, tangential and effective separations δnn − 1, δtn − 1, δ n − 1 are passed into the CDFE program. It should be specifically pointed out that, the input θ is a constant in enhanced CDFE, since the DCZM is fixed from the start of the program. The variable cs means the crack status. If cs is equal to 0, the element is still pristine. If cs is equal to 1, the crack has initiated in the element. If cs is equal to 2, the element has failed. At the start of each increment, cs (n − 1) is first checked to see which state the element is at. At the end of the increment, the right-hand side RHSe(n) , updated crack status cs (n) , displacement vector of the dummy nodes U∧ (n) and normal, tangential and effective separations δnn, δtn, δ n , are output and passed back to the FEM solver.

(20)

δn2 + (βδt )2

(21)

t1 (δ ) =

σC δ (δ ⩽ δn0), δn0

(22a)

t1 (δ ) =

δ − δn0 σC (δn0 ⩽ δ < δnC ), δnC − δn0

(22b)

t1 (δ ) = 0 (δnC < δ )

(22c)

t2 (δ ) =

βτC δ (δ ⩽ βδt 0), βδt 0

(23a)

t2 (δ ) =

δ − βδt 0 βτC (βδt 0 ⩽ δ < βδtC ), βδtC − βδt 0

(23b)

t2 (δ ) = 0 (βδtC < δ )

(27)

(23c) 371

Composite Structures 212 (2019) 365–380

S. Lin et al.

INPUT: X, U(n), E11, E22, Ȟ12, G12, ıC, IJC, GIC, GIIC, ș, cs(n-1), U^(n-1), įn(n-1), įt(n-1), į(n-1)

2GIC /

nC

K1

į

,

tC

2GIIC /

ȍ1

į

(n 1) 2 n

cs(n-1)==0

C

,

nC

/

tC

T

B1 DB1dV , K 2

(n 1)

IF

C

T

ȍ2

B 2 DB 2 dV

(n 1) 2 t

THEN

Determine K 57 n , K 57 t , K 68 n , K 68t from equation (14) Determine K DCZM from equation (12) Assemble into global stiffness matrix K Static condensation : K en RHS (n) e į (n) n

K *^

^*

K ^^

K

K ** K *^ (K ^ ^ ) 1 K ^*

(n) K (n) Feext(n) , U^(n) e Ue

L n RU^(n) , į t(n)

K **

(K ^ ^ ) 1 K ^*U (n)

L t RU^(n) , į (n)

į (n) n

2

(n) 2 t

IF į(n)<įn0 THEN ELSE

cs (n)

0

cs (n)

1

ENDIF ELSEIF cs(n-1)==1 THEN

Determine K 57 n , K 57 t , K 68 n , K 68t from equation (14) Determine K DCZM from equation (12) Assemble into global stiffness matrix K Static condensation : K en RHS (n) e į (n) n

K *^

K ^*

K ^^

K ** K *^ (K ^ ^ ) 1 K ^*

(n) K (n) Feext(n) , U^(n) e Ue

L n RU^(n) , į t(n)

K **

(K ^ ^ ) 1 K ^*U (n)

L t RU^(n) , į (n)

į (n) n

2

(n) 2 t

IF į(n)<įnC THEN

cs (n) ELSE

cs (n)

1 2

ENDIF ELSEIF cs(n-1)==2 THEN

K en RHS (n) e

T

ȍe

B e DB e dV / 10 20 (n) K (n) Feext(n) e Ue

ENDIF OUTPUT: RHS (n) , cs(n), U^(n), įn(n), įt(n), į(n) e

Fig. 7. Flow chart of the enhanced CDFE for composite damage analysis.

3. Applications

microscope (SEM) [42]. Initiation of micro-cracks and formation of macro-cracks were observed. Several numerical tools based on FEM have been developed for the problem. Vaughan and McCarthy modeled the transverse damage by seeding discrete layers of interfacial cohesive elements on the fiber–matrix interfaces [43]. However, in this model, the crack cannot propagate into the matrix from the interfaces. This feature makes the model not suitable for modern composite materials due to the technology of fiber sizing, plasma surface modification and carbon nanotube coating, et al. [44] which elevates the interface toughness significantly. The fiber/matrix adhesion is greatly increased by these technologies such that the adhesion energy of the interface can

3.1. Micromechanical crack growth The micromechanical crack growth analysis of a representative volume element (RVE) with 13 randomly packed fibers is performed to show the capability of the CDFE to model crack propagation in cases without well-defined crack angle and path. The micromechanical transverse crack growth in composite materials has been studied by many. Gamstedt and Sjogren performed in situ inspection of tensile failure of a unidirectional composite using scanning electron

372

Composite Structures 212 (2019) 365–380

S. Lin et al.

Table 1 Elastic and fracture properties of the glass fiber and epoxy matrix used in the RVE [6]. Property

Value

Property

Value

Ef (GPa) νf

74.0

Gm (GPa)

1.79

0.2

εCm

0.0143

Gf (GPa)

30.8

γCm

0.0237

Em (GPa)

4.65

m (N/mm) GIC

0.000563

νm

0.35

m (N/mm) GIIC

0.00385

Extends to be larger than the cohesion energy of matrix [45], hence the fiber/ matrix debonding described in [43] is not likely to happen. Therefore, the matrix itself should be discretized with crack elements capable of modeling arbitrary crack initiation and propagation. Following this idea, Pineda et al. implemented the crack band model within the framework of high-fidelity generalized method of cells (HFGMC) and FEM to study a RVE containing 13 randomly arranged fibers under transverse tension, compression and shear loadings [6]. Kaleel et al. utilized the component-wise (CW) approach, which is an extension of the Carrera unified formulation (CUF) [46], combined with crack band method to study the progressive failure of the same RVE [47]. Detailed modeling results of crack propagation was reported. In this paper, transverse tensile crack formation in the RVE studied in [6] is modeled with CDFE. The size of the RVE is 21.25 μm × 21.25 μm . The material properties for the fiber and matrix are listed in Table 1. Same as in [6], plane strain material stiffness matrix is implemented in CDFE. As mentioned before, the maximum principal stress criterion is used to decide crack orientation and initiation. For the purpose of comparison, and to validate the effectiveness of CDFE, periodic boundary condition is imposed on the RVE, as in [6,47], using kinematic constraint of the DOFs of the corresponding node pairs and a specified reference node, as in Eq. (28). However, whether periodicity can be applied to this problem is an open question.

ui1 − ui2 = uiRF

Fig. 8. Comparison of the bi-linear and exponential traction-separation laws.

In Fig. 10, the x-direction displacement discontinuity is clearly observed. Among three meshes, the crack discontinuity paths are different. However, for the same mesh with bi-linear and exponential traction-separation laws, the paths are essentially identical. A more detailed inspection can be conducted by explicitly plotting the dummy crack nodes, which are illustrated in Fig. 3 as node 5, 6, 7 and 8. The storage of the crack dummy node information throughout the analysis is one of the major advantages of the CDFE method. The CDFE results at the end of loading step with dummy nodes plotted are in Fig. 11. It should be noted that only coordinates of regular nodes and dummy nodes are utilized for plotting, without any field variable, element deletion or coloring scheme. From Fig. 11, macroscopic cracks cutting through the RVE are clearly seen. The macroscopic cracks are formed with dummy cracks inside each cracked element. At some locations, the dummy cracks among elements are quite continuous. Nevertheless, in each case, the continuity is interrupted by several dummy cracks not aligned well with the macroscopic crack path to various extents. A more quantified comparison can be drawn in terms of stress–strain responses of three RVE meshes with two traction-separation laws and Pineda’s results [6], as shown in Fig. 12. From Fig. 12, a good coincidence is found among all curves. Critical parameters including initial stiffness, peak stress and failure strain are listed and compared in Table 2. Root-mean-square deviation (RMSE) and maximum error of the values are measured to investigate the variance of the data. The RMSEs for three critical parameters are 8.9%, 11.4% and 4.8%. The maximum errors are 8.9%, 17.4% and 10.2%. It is worthwhile to observe the dummy crack nodes with more details. As a representative the result of the mesh size 0.6 μm with bilinear traction-separation law is plotted with several locations magnified in Fig. 13. In Fig. 13(b), the dummy cracks of the cracked elements are well aligned, forming an clear crack path. In Fig. 13(c), although crack path is still visible, a saw-tooth-like crack interface pattern is seen. This is due to the fact that no crack continuity is enforced across elemental boundaries. In Fig. 13(d), the macroscopic crack path is interrupted with a cracked element having dummy cracks not well aligned. This interrupting cracked element will be illustrated more in combination with a careful observation of the crack development history of the RVE. For the case of the mesh size 0.6 μm with bi-linear traction-separation law, historical data of the dummy crack nodes is preserved and plotted with the stress–strain history to thoroughly understand the mechanical behavior of the cracking RVE in Fig. 14. The progressive failure of the RVE is as follows. First, as seen in Fig. 14(a), no dummy crack is present in the RVE and from 14(f), before the point (a), the stress–strain curve is linear, meaning that there is no damage

(28)

In Eq. (28), the superscripts 1 and 2 represent the first and second node of a certain node pair. RF represents the reference node. The subscript i represents the ith DOF of the node. To study the sensitivity of CDFE to mesh sizes and the types of the traction-separation law, three meshes are seeded on the same RVE and each is studied with bi-linear and exponential traction-separation law. The bi-linear traction-separation law is as shown in Eqs. (6), (7) and Fig. 2. Only the mode I exponential traction-separation law is written in Eq. (29), since the arbitrary crack growth in RVE is a mode I dominated problem. The required input parameters are the same as that for the bilinear traction-separation law. The comparison of bi-linear and exponential traction-separation laws is shown in Fig. 8.

tn (δn ) =

σC δn (δn ⩽ δn0), δn0

tn (δn ) = σCexp (

−σC (δn − δn0)) (δn0 ⩽ δn ⩽ δnC ), GC

tn (δn ) = 0 (δnC ⩽ δn )

(29a) (29b) (29c)

The three meshes of the RVE are shown in Fig. 9. The mesh sizes are 1 μm, 0.6 μm and 0.2 μm. There are three types of the element, including the reduced integration plane strain element CPE4R for the fiber, CDFE user element for the matrix, and dummy elements for the post-processing purpose. Extremely low stiffness is assigned to the dummy elements such that no contribution is made to the crack analysis from these elements. A global strain of 0.004 in the global 2nd direction is applied to the RVE. The displacement fields in the x direction of the three meshes with bi-linear and exponential traction-separation laws are shown in Fig. 10. 373

Composite Structures 212 (2019) 365–380

S. Lin et al.

Fig. 9. Meshes of the RVE: (a) mesh 1, (b) mesh 2, and (c) mesh 3.

attention is paid to the interrupting cracked element highlighted in Fig. 15b) and the dummy crack development inside this element through the analysis is presented in Fig. 15(c). In Fig. 15(c), the red frame represent the four edges of the quadrilateral CDFE element. In stage I, there is no crack initiated, which correspond to Fig. 14(a). In stage II, the crack initiate inside the element and the crack separation is so small that the edges of the sub-elements are almost overlapping. This stage correspond to Fig. 14(b) and (c), where the drastic stress drop happens. From stage III to IV, one of the sub-element is sheared outside the CDFE elemental boundary (the red frame) and causes the interruption of the continuity of the macroscopic crack. The crack initiation and deformation of the sub-elements imply that the stress state of the cracked element varies constantly through the analysis, due to the crack initiation and propagation at other locations, and this would lead to shearing of the sub-elements after crack initiation since once the CDFE element is initiated, the crack orientation is fixed until the end of the analysis. This phenomenon is due to the fact that no constraint is applied currently to refrain the sub-element from reaching into the

happening to lower the stiffness of the RVE. From Fig. 14(a)–(b), cracks initiate and open, although with hardly visible separations, at several locations which are mostly between closely spaced fibers. These places are where stress concentration is most likely to occur. Fig. 14(b) correlate to the point (b) in Fig. 14(f), which is the peak of the stress–strain curve. Between point (a) and (b), the nonlinearity is seen, due to the initiation and opening of multiple cracks. Right after the point (b), there is a drastic drop of the stress, and this corresponds to the connection of the multiple cracks from Fig. 14(b)–(c). As seen in Fig. 14(c), the connected macroscopic crack almost cut through the RVE and this marks the failure of the RVE. Different from the dummy crack plots in Fig. 13, at stage (c), the dummy cracks in the cracked elements align with the macroscopic path quite well. However, from Fig. 14(d)–(e), some interrupting cracked elements are seen, and in Fig. 14(f), the stress becomes oscillating around a trivial value. The reason behind the interruption of the macroscopic crack path deserves more attention. Again, the result of the mesh size 0.6 μm with bi-linear tractionseparation law is magnified locally in Fig. 15(a) and (b). Special

Mesh: 1 ȝm

Mesh: 0.6 ȝm

(a)

(b)

(c)

(d)

(e)

(f)

Mesh: 0.2 ȝm

Bi-linear

Exponential

Fig. 10. Displacement fields of the RVE: (a) Mesh: 1.0 μm/Bi-linear, (b) Mesh: 0.6 μm/Bi-linear, (c) Mesh: 0.2 μm/Bi-linear, (d) Mesh: 1.0 μm/Exponential, (e) Mesh: 0.6 μm/Exponential and (e) Mesh: 0.2 μm/Exponential. 374

Composite Structures 212 (2019) 365–380

S. Lin et al.

Fig. 11. Crack plots with dummy nodes of the RVE: (a) Mesh: 1.0 μm/Bi-linear, (b) Mesh: 0.6 μm/Bi-linear, (c) Mesh: 0.2 μm/Bi-linear, (d) Mesh: 1.0 μm/ Exponential, (e) Mesh: 0.6 μm/Exponential and (e) Mesh: 0.2 μm/Exponential.

(b)

(d)

(c) (a) Fig. 13. Crack plots with dummy nodes: (a) of the whole mesh, (b) where cracks are forming a continuous pattern, (c) where cracks are forming a sawtooth-like pattern, and (d) where the macroscopic crack is interrupted.

Fig. 12. Comparison of the stress strain curves obtained with CDFE and HFGMC.

Table 2 Comparison of the critical parameters. CDFE 1 μm

Initial stiffness (GPa) Peak load (MPa) Failure strain

0.6 μm

RMSE

Max error

8.9% 11.4% 4.8%

8.9% 17.4% 10.2%

0.2 μm

Bi.

Exp.

Bi.

Exp.

Bi.

Exp.

18.3 36.1 2.78e−3

18.3 34.4 2.82e−3

18.3 42.1 2.75e−3

18.3 37.8 2.86e−3

18.3 36.0 3.03e−3

18.3 33.2 2.67e−3

375

16.8 40.2 2.75e−3

Composite Structures 212 (2019) 365–380

S. Lin et al.

(b)

(a)

(c) (b)

(a) (c)

(d)

(e)

(d)

(e)

(f)

Fig. 14. Transverse crack growth in the RVE.

through the RVE. To evaluate the CDFE method, three meshes of the mesh size 1 μm, 0.6 μm and 0.2 μm are studied, with two types of traction-separation law implemented, including the bi-linear law and exponential law. The results are compared against published results in [6]. According to the CDFE results, the crack paths predicted by CDFE is sensitive to the mesh, since the crack orientation prediction is highly dependent on the stress state of each element. However, the CDFE method has been proved to be mesh objective in terms of stress–strain responses, and the stress–strain curves obtained with CDFE agree well with the HFGMC results in [6]. RMSEs and maximum errors are measured for three critical parameters, including the initial stiffness, peak stress and the failure strain. The RMSEs are within 12% and the maximum errors are within 18%. It is noticed that the macroscopic crack is usually interrupted by some cracked elements, and the cause is found to be the sub-elements of the cracked CDFE element shearing outside the element boundary. This is due to the fact that no crack continuity across elemental edges is enforced in the CDFE at the present time. Future studies will evaluate the efficiency of implementing such a feature.

(b)

(a)

I

II

III

IV

(c) 3.2. Delamination toughness tests Fig. 15. Detailed history of element cracking in the RVE.

In Section 3.1, the crack growth was dominated by mode I cracking. No mode mixity is present in the problem. As illustrated in Section 2.2, enhanced CDFE is developed to model the mixed mode crack growth. Therefore, in this section, delamination toughness tests with well-defined crack plane and mode mixity are modeled with CDFE and enhanced CDFE to validate the improvement of numerical stability in the results. The delamination toughness tests are already equipped with benchmark results obtained using virtual crack closure technique (VCCT) [49], analytical solutions [50,51], simulations [52,53] and experimental results [54], The tests include a double cantilever beam bending (DCB) test for pure mode I cracking, an end notched flexure (ENF) test for mode II cracking and the mixed mode bending (MMB) test for mixed mode cracking. In this paper, delamination toughness tests in [49] are modeled with CDFE and enhanced CDFE. The results

bordering elements. For example, in Fig. 3, applying constraints on edge 1–5-7–4 and 2-6-8-3 to constrain the node 1, 5, 7, 4 and 2, 6, 8, 3 to be on the same straight edges can prevent the shearing-out issue. The constraints can be implemented based on the penalty method or the Lagrangian multiplier method [48]. After obtaining the stiffness matrices from the applied constraints, the constraint stiffness matrices can be assembled into the global stiffness matrix as in Eq. (17) and condensed into the standard stiffness matrix as in Eqs. (18). From the CDFE analyses of the tensile failure of the RVE, the mechanism of the transverse crack initiation and propagation has been investigated in detail. The drastic nature of failure is captured and the mechanism of the failure is found out as the connection of multiple cracks at different locations into one microscopic crack that cuts

376

Composite Structures 212 (2019) 365–380

S. Lin et al.

Fig. 16. Set-up for the delamination toughness tests: (a) DCB, (b) ENF and (c) MMB.

are also verified with existing benchmark solutions. All the simulations use the same time increment size of 5e−7 s . It should be specifically pointed out that, in the pure mode I (DCB) and mode II (ENF) cases, the results of the enhanced CDFE should be identical to that of CDFE, since the mixed mode formulation illustrated in Section 2.2.3 would converge to single mode formulation in single mode cases. However, in mixed mode (MMB) cases, differences in the results obtained by CDFE and enhanced CDFE are expected and should be highlighted.

Table 4 Material properties of T300/1076 [55].

60

Force (N)

50

40

30

20 Krueger (2015) CDFE with Power Load CDFE with the Mixed Mode Formulation

10

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

End Displacement (mm) Fig. 17. Load–displacement curve of the DCB case.

Table 5 Geometry parameters for the ENF case [49].

3.2.2. ENF The set-up of the ENF case is shown in Fig. 16(b). The geometry

a L h b

Table 3 Geometry parameters for the DCB case [49]. Initial crack length Half support distance Half thickness Out-of-plane width

4.6 GPa 3.54 GPa 0.170 N/mm 0.494 N/mm

G12 = G13 G23 GIC GIIC

70

3.2.1. DCB The set-up of the DCB case is shown in Fig. 16(a). The geometry parameters and material properties follow [49] and listed in Table 3 and 4. To enable CDFE explicit dynamic analyses, material density ρ , mode I and mode II crack strengths are introduced from [51], as ρ = 1.495 g/cm3, σC = 50 MPa , τC = 70 MPa . The load–displacement curves obtained with CDFE and enhanced CDFE are compared and validated with the benchmark solution provided in [49] in Fig. 17. As seen in Fig. 17, in the pre-peak region, the initial stiffnesses agree very well with the benchmark solution. The peak load of the benchmark solution is 62.0 N, and of CDFE and enhanced CDFE results are both 61.9N. The error is 0.16%. In the post-peak region, both CDEF and enhanced CDFE results are oscillatory, which is due to the nature of explicit algorithm. In general, for the DCB case, the CDFE and enhanced CDFE results agree very well with the benchmark solution. The curves obtained by CDFE are almost right on top of the curve obtained with enhanced CDFE. This is expected, since there is only mode I cracking in the DCB test, the mixed mode cohesive formulation would converge to the single mode formulation.

a L h b

139.4 GPa 10.16 GPa 0.3 0.436

E11 E22 = E33 ν12 = ν13 ν13

30.5 mm 75.0 mm 1.5 mm 25.0 mm

Initial crack length Half support distance Half thickness Out-of-plane width

25.4 mm 50.8 mm 2.25 mm 25.4 mm

parameters and material properties are listed in Table 5 and 6. To enable the CDFE explicit analyses, material density ρ , mode I and mode II crack strengths are introduced from [51], as 377

Composite Structures 212 (2019) 365–380

S. Lin et al.

Table 6 Material properties of IM7/8552 [56]. 161 GPa 11.38 GPa 0.32 0.45

E11 E22 = E33 ν12 = ν13 ν13

Table 7 Geometry parameters for the MMB cases [49]. 5.2 GPa 3.9 GPa 0.212 N/mm 0.774 N/mm

G12 = G13 G23 GIC GIIC

a L h b

Initial crack length Half support distance Half thickness Out-of-plane width

c

Length of lever arm

25.4 mm 50.8 mm 2.25 mm 25.4 mm R = 20%: 92.9 R = 50%: 41.3 mm R = 80%: 27.3

1800 Krueger (2015) CDFE with Power Law CDFE with the Mixed Mode Cohesive Formulation

1600 1400

Force (N)

1200

R=80%

1000 800 600

R=50%

400 200

R=20%

0 0

0.5

1

1.5

2

2.5

3

Mid Displacement (mm) Fig. 18. Load–displacement curve of the ENF case.

ρ = 1.495 g/cm3, σC = 50 MPa , τC = 70 MPa . The load–displacement curves obtained with CDFE and enhanced CDFE are compared with the benchmark solution in Fig. 18. As seen in Fig. 18, the initial stiffnesses agree well. The peak load of the benchmark solution is 1526.0 N, and of CDFE and enhanced CDFE are all 1381.0 N. The error is 9.50%. In both CDFE and enhanced CDFE results, nearly vertical drops are seen right after the peak point. This usually happens in explicit dynamic modeling of the ENF test, since the ENF test has the snap-back behavior at initiation, which is an unstable response and a snap-back response cannot be captured by an explicit algorithm. In the post-peak region, both the CDFE and enhanced CDFE results show oscillations. The magnitude of oscillations is larger than that of the DCB test, since the snap-back behavior makes the ENF test unstable compared to the stable growth in the DCB test. Although the ENF test is pure mode II cracking, from Fig. 18, it is still found that the CDFE result is more oscillatory than that of the enhanced CDFE. This is due to several factors. First, in the ENF modeling, a slight mode mixity is usually expected. Second, as illustrated in Section 2.2.1, at initiation, extra DOFs are introduced to the CDFE element, and hence a sudden change of elemental stiffness is induced. However, in an enhanced CDFE element, the number of DOFs remain constant throughout the analysis. The sudden elemental stiffness degradation would make the CDFE result more oscillatory than that of the enhanced CDFE.

Fig. 19. Load–displacement curve of three MMB cases. Table 8 Load–displacement curve of three MMB cases. R 20% 50% 80%

Benchmark (N)

CDFE (N)

Enhanced CDFE (N)

Max difference

128.5 385.0 751.0

130.3 370.0 750.5

129.2 363.4 738.1

1.40% 5.45% 1.72%

compared in Table 8. As can be found in Table 8, the maximum difference 5.45% happens at the case of the mixed mode ratio equal to 50%. The difference between CDFE and enhanced CDFE results in the post-peak regions should be paid special attention to. Generally, it is shown that, the CDFE results are more oscillatory than that of the enhanced CDFE results. As the case gets more inclined to the pure-mode crack growth problem (R = 20% and 80%), the differences between the results become smaller. This can be explained by noting that CDFE and enhanced CDFE differ in their damage evolution formulation when mode mixity is encountered. In CDFE, when the first-order power law is satisfied, the tractions are forced to be zero even though finite remaining values exist. However, in enhanced CDFE, the implementation of the mixed mode formulation guarantees that tractions of both mode I and mode II vanish simultaneously and smoothly. Therefore, the enhanced CDFE results are more numerically stable. And, when the MMB case gets more inclined to the single mode case, the difference between CDFE and enhanced CDFE should be smaller. Also, as discussed in the ENF case, the sudden elemental stiffness variation at crack initiation in CDFE should also contribute to the differences between the results, however, compared with the mixed mode cohesive formulation, this effect is negligible, [57]. From the delamination toughness tests analyzed with CDFE and enhanced CDFE, for all cases, the modeling results agree well with the benchmark solutions obtained with VCCT. The maximum error of the peak load values from the benchmark solution is below 10%. For the DCB case, the curve obtained with CDEF is right on top of that with

3.2.3. MMB The set-up of the MMB cases is shown in Fig. 16(c). There are three mixed mode ratios studied. The definition of the mixed mode ratio is R = GII /(GI + GII ) , and the relation between geometry parameters and R can be found in the appendix of [52]. The geometry information can be found in Table 7, and the material properties are listed in Table 6. To enable the CDFE explicit analyses, material density ρ , mode I and mode II crack strengths are introduced from [51], as ρ = 1.495 g/cm3, σC = 50 MPa , τC = 70 MPa . The load–displacement curves obtained with CDFE and enhanced CDFE are compared with the benchmark solutions in Fig. 19. As seen in Fig. 19, for all three mixed mode ratios, initial stiffnesses agree very well. The peak loads of the benchmark solutions, CDFE and enhanced CDFE results are listed and 378

Composite Structures 212 (2019) 365–380

S. Lin et al.

enhanced CDFE. For the ENF case, in the post-peak region, around the load drop point, the CDFE result is more oscillatory than that of enhanced CDFE. The difference is mainly credited to the instantaneous elemental stiffness variation in CDFE. For the MMB cases, in the postpeak regions, the CDFE results are more oscillatory than that of enhanced CDFE. As the cases gets more inclined to the pure mode case, the difference gets smaller. These differences are mainly caused by the implementation of the mixed mode cohesive formulation in enhanced CDFE. The MMB cases should be highlighted and regarded as validation in the improvement in numerical stability of CDFE analyses. The improvement is contributed by the mixed mode cohesive formulation in the enhanced CDFE.

analysis of composite materials with efficiency and accuracy. In the future, a full 3D version of CDFE can be developed to model general laminate progressive failure problems under different loadings such as low-velocity impact. The capability of CDFE to model micromechanical failures can be further investigated and a multiscale framework can be developed based on CDFE. Acknowledgments The study reported here was completed at the University of Washington, Seattle, when all of the authors were with the William E Boeing Department of Aeronautics and Astronautics. We gratefully acknowledge the financial support and the support of the HYAK computational cluster in carrying out the study reported here.

4. Conclusions Continuum decohesvie finite element (CDFE) has been applied to progressive failure analysis of composite materials. The first problem studied is the micromechanical transverse crack growth in an RVE randomly packed with 13 fibers. In this problem, the crack angle is detected by the maximum principal stress criterion and mode I crack is dominating. The second problem set is the delamination toughness tests, including a DCB, ENF and three MMB tests. In these problems, crack plane is well defined but mode mixity can be encountered. To acquire better numerical stability in the analysis results in mixed mode situations, CDFE is enhanced with a new inner-element discretization scheme and a novel mixed mode cohesive formulation to guarantee that both mode I and mode II tractions vanish simultaneously and smoothly at crack growth. Mathematical formulation, damage evolution scheme, criteria and FEM scheme for both CDFE and enhanced CDFE are provided in Section 2. From the micromechanical transverse crack growth analysis, the mechanism of multiple cracks initiating at different sensitive locations and connecting into a macroscopic transverse crack is modeled and the drastic nature of the process is also captured. The stress–strain responses are objective to the mesh size and type of the traction-separation law. The CDFE results agree well with published results obtained with HFGMC. The maximum error in critical parameters including the initial stiffness, peak stress and failure strain is within 18%. From the crack plots with dummy nodes, right after the drastic drop of the stress, the cracks are aligned well and forming a transverse macroscopic crack path cutting through the RVE. However, with the RVE getting further loaded, sub-elements of some already cracked CDFE elements tend to shear out of the CDFE elemental boundaries, and this would cause the interruption of the macroscopic crack path. This phenomenon is caused by the fact that the current CDFE method does not enforce crack continuity across elemental boundaries. A remedy for this issue is to apply constraints on elemental edges to prevent the cohesive crack from shearing out of the elemental boundaries. The constraints implementation can be based on penalty stiffness method or Lagrangian multiplier method. After obtaining the constraint stiffness matrices, assembly into the elemental global stiffness matrix and condensation into standard elemental stiffness matrix is to be performed. From analyses of the delamination toughness tests, all the results obtained with CDFE and enhanced CDFE agree very well with the benchmark solutions acquired with VCCT. The maximum error of the peak load is within 10%, happening at the ENF test. In the pure mode I test (DCB), the CDFE result is right on top of the enhanced CDFE result. In the pure mode II test (ENF), right after the peak point on the load–displacement curve, the CDFE result is more oscillatory than the enhanced CDFE result. This is believed to be caused by the difference in the inner-discretization scheme in the two versions of CDFE. In the mixed mode tests (MMB), from all three groups of modelings, the CDFE results are more oscillatory than the enhanced CDFE results. Also, as the test gets more inclined to a pure-mode case, the difference between the results become smaller. Through these two sets of examples, CDFE has been proved to be capable of dealing with general progressive failure

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.compstruct.2019.01.021. References [1] Liu W, Yang Q, Mohammadizadeh S, Su X, Ling D. An accurate and efficient augmented finite element method for arbitrary crack interactions. J Appl Mech 2013;80(4):041033. [2] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 2: numerical simulation. Int J Numer Methods Eng 1996;39(21):3601–23. [3] Rudraraju S, Salvi A, Garikipati K, Waas AM. Predictions of crack propagation using a variational multiscale approach and its application to fracture in laminated fiber reinforced composites. Compos Struct 2012;94(11):3336–46. [4] Pineda EJ, Waas AM. Numerical implementation of a multiple-isv thermodynamically-based work potential theory for modeling progressive damage and failure in fiber-reinforced laminates. Int J Fract 2013;182(1):93–122. [5] Meyer P, Waas AM. Mesh-objective two-scale finite element analysis of damage and failure in ceramic matrix composites. Integr Mater Manuf Innov 2015:1–18. [6] Pineda EJ, Bednarcyk BA, Waas AM, Arnold SM. Progressive failure of a unidirectional fiber-reinforced composite using the method of cells: discretization objective computational results. Int J Solids Struct 2013;50(9):1203–16. [7] Totry E, González C, LLorca J. Prediction of the failure locus of c/peek composites under transverse compression and longitudinal shear through computational micromechanics. Compos Sci Technol 2008;68(15–16):3128–36. [8] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17–18):1526–35. [9] Hu W, Ha YD, Bobaru F. Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites. Comput Methods Appl Mech Eng 2012;217:247–61. [10] Sun C, Huang Z. Peridynamic simulation to impacting damage in composite laminate. Compos Struct 2016;138:335–41. [11] Medina D, Chen J. Three-dimensional simulations of impact induced damage in composite structures using the parallelized sph method. Compos Part A: Appl Sci Manuf 2000;31(8):853–60. [12] Shintate K, Sekine H. Numerical simulation of hypervelocity impacts of a projectile on laminated composite plate targets by means of improved sph method. Compos Part A: Appl Sci Manuf 2004;35(6):683–92. [13] Liew KM, Zhao X, Ferreira AJ. A review of meshless methods for laminated and functionally graded plates and shells. Compos Struct 2011;93(8):2031–41. [14] Bažant ZP, Oh BH. Crack band theory for fracture of concrete. Matér Constr 1983;16(3):155–77. [15] Rots JG, et al. Smeared crack approach and fracture localization in concrete. HERON 1985;30(1). [16] Schapery RA. A theory of crack initiation and growth in viscoelastic media. Int J Fract 1975;11(1):141–59. [17] Joseph AP, Davidson P, Waas AM. Intra-inter crack band model (i2cbm) for progressive damage and failure analysis of joints. 2017. p. 0200. [18] Joseph AP, Davidson P, Waas AM. Progressive damage and failure analysis of single lap shear and double lap shear bolted joints. Compos Part A: Appl Sci Manuf 2018. [19] Joseph AP, Davidson P, Waas AM. Open hole and filled hole progressive damage and failure analysis of composite laminates with a countersunk hole. Compos Struct 2018;203:523–38. [20] Xie D, Waas AM, Shahwan KW, Schroeder JA, Boeman RG. Computation of energy release rates for kinking cracks based on virtual crack closure technique. Comput Model Eng Sci 2004;6(6):515–24. [21] Xie D, Waas AM, Shahwan KW, Schroeder JA, Boeman RG. Fracture criterion for kinking cracks in a tri-material adhesively bonded joint under mixed mode loading. Eng Fract Mech 2005;72(16):2487–504. [22] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46(1):131–50. [23] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69(7):813–33.

379

Composite Structures 212 (2019) 365–380

S. Lin et al. [24] Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer methods in applied mechanics and engineering 2004;193(33–35):3523–40. [25] Rudraraju SS, Salvi A, Garikipati K, Waas A. Mixed mode in-plane fracture analysis of laminated fiber reinforced composites using the variational multiscale cohesive method. 2010. p. 2860. [26] Oliver J, Huespe AE, Sánchez PJ. A comparative study on finite elements for capturing strong discontinuities: E-fem vs x-fem. Comput Methods Appl Mech Eng 2006;195(37–40):4732–52. [27] Ling D, Yang Q, Cox B. An augmented finite element method for modeling arbitrary discontinuities in composite materials. Int J Fract 2009;156(1):53–73. [28] Liu W, Schesser D, Yang Q, Ling D. A consistency-check based algorithm for element condensation in augmented finite element methods for fracture analysis. Eng Fract Mech 2015;139:78–97. [29] Naderi M, Iyyer N. 3d modeling of arbitrary cracking in solids using augmented finite element method. Compos Struct 2017;160:220–31. [30] Chen B, Pinho S, De Carvalho N, Baiz P, Tay T. A floating node method for the modelling of discontinuities in composites. Eng Fract Mech 2014;127:104–34. [31] Chen B, Tay T, Pinho S, Tan V. Modelling the tensile failure of composites with the floating node method. Comput Methods Appl Mech Eng 2016;308:414–42. [32] Li X, Chen J. The implementation of the extended cohesive damage model for multicrack evolution in laminated composites. Compos Struct 2016;139:68–76. [33] Li X, Chen J. An extended cohesive damage model for simulating multicrack propagation in fibre composites. Compos Struct 2016;143:1–8. [34] Li X, Chen J. An extended cohesive damage model for simulating arbitrary damage propagation in engineering materials. Comput Methods Appl Mech Eng 2017;315:744–59. [35] Li X, Chen J. A highly efficient prediction of delamination migration in laminated composites using the extended cohesive damage model. Compos Struct 2017;160:712–21. [36] McElroy M, Jackson W, Olsson R, Hellström P, Tsampas S, Pankow M. Interaction of delaminations and matrix cracks in a cfrp plate, part i: A test method for model validation. Compos Part A: Appl Sci Manuf 2017;103:314–26. [37] Prabhakar P, Waas AM. A novel continuum-decohesive finite element for modeling in-plane fracture in fiber reinforced composites. Compos Sci Technol 2013;83:1–10. [38] Nguyen N, Waas AM. Continuum decohesive finite element modeling of in-plane fracture: Mesh-objectivity and sensitivity studies. 2017. [39] Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Fract Mech 2006;73(13):1783–96. [40] Gustafson PA, Waas AM. The influence of adhesive constitutive parameters in cohesive zone finite element models of adhesively bonded joints. Int J Solids Struct

2009;46(10):2201–15. [41] Nguyen N, Waas AM. A novel mixed-mode cohesive formulation for crack growth analysis. Compos Struct 2016;156:253–62. [42] Gamstedt E, Sjögren B. Micromechanisms in tension-compression fatigue of composite laminates containing transverse plies. Compos Sci Technol 1999;59(2):167–78. [43] Vaughan T, McCarthy C. Micromechanical modelling of the transverse damage behaviour in fibre reinforced composites. Compos Sci Technol 2011;71(3):388–96. [44] Sharma M, Gao S, Mäder E, Sharma H, Wei LY, Bijwe J. Carbon fiber surfaces and composite interphases. Compos Sci Technol 2014;102:35–50. [45] Paredes J, Martínez-Alonso A, Tascon J. Oxygen plasma modification of submicron vapor grown carbon fibers as studied by scanning tunneling microscopy. Carbon 2002;40(7):1101–8. [46] Carrera E, Cinefra M, Petrolo M, Zappino E. Finite element analysis of structures through unified formulation. John Wiley & Sons; 2014. [47] Kaleel I, Petrolo M, Waas A, Carrera E. Micromechanical progressive failure analysis of fiber-reinforced composite using refined beam models. J Appl Mech 2018;85(2):021004. [48] Belytschko T, Liu WK, Moran B, Elkhodary K. Nonlinear finite elements for continua and structures. John wiley & sons; 2013. [49] Krueger R. A summary of benchmark examples to assess the performance of quasistatic delamination propagation prediction capabilities in finite element codes. J Compos Mater 2015;49(26):3297–316. [50] Xie J, Waas AM, Rassaian M. Closed-form solutions for cohesive zone modeling of delamination toughness tests. Int J Solids Struct 2016;88:379–400. [51] Xie J, Waas AM, Rassaian M. Estimating the process zone length of fracture tests used in characterizing composites. Int J Solids Struct 2016;100:111–26. [52] Camanho PP, Dávila CG. Mixed-mode decohesion finite elements for the simulation of delamination in composite materials. 2002. [53] Alfano G, Crisfield M. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int J Numer Methods Eng 2001;50(7):1701–36. [54] Reeder JR, REWS JH. Mixed-mode bending method for delamination testing. AiAA J 1990;28(7):1270–6. [55] Krueger R. An approach to assess delamination propagation simulation capabilities in commercial finite element codes. 2008. [56] Krueger R. Development and application of benchmark examples for mixed-mode i/ ii quasi-static delamination propagation predictions. 2012. [57] Lin S, Waas AM. Using the continuum decohesive finite element for crack growth analysis in fiber reinforced composites. 2018. p. 0731.

380