Analytical predictions of delamination threshold load of laminated composite plates subject to flexural loading

Accepted Manuscript Analytical Predictions of Delamination Threshold Load of Laminated Composite Plates subject to Flexural Loading Jiawen Xie, Anthony M. Waas, Mostafa Rassaian PII: DOI: Reference:

S0263-8223(17)31255-2 http://dx.doi.org/10.1016/j.compstruct.2017.07.009 COST 8664

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

20 April 2017 5 July 2017 8 July 2017

Please cite this article as: Xie, J., Waas, A.M., Rassaian, M., Analytical Predictions of Delamination Threshold Load of Laminated Composite Plates subject to Flexural Loading, Composite Structures (2017), doi: http://dx.doi.org/ 10.1016/j.compstruct.2017.07.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Analytical Predictions of Delamination Threshold Load of Laminated Composite Plates subject to Flexural Loading Jiawen Xiea , Anthony M. Waasb,∗, Mostafa Rassaianc a Dept.

of Aerospace Engineering, University of Michigan, Ann Arbor, MI 48109-2140 E. Boeing Dept. of Aeronautics and Astronautics, University of Washington, Seattle, WA 98195-2400 c Boeing Company, Seattle, WA, 98124

b William

Abstract An analytical approach is proposed to determine delamination threshold loads of fiber-reinforced laminated composite plates with arbitrary stacking sequences under transverse loading conditions. Following the concept of cohesive zone modeling, a laminated plate is considered as an assembly of two sub-laminates connected by a virtual elastic-brittle layer with infinitesimal thickness. The problem is formulated and solved by the Rayleigh-Ritz method based on firstorder shear deformation theory. The problem of quasi-static face-on (transverse) indentation test is analyzed as an example. The results, including elastic stiffness of flexural response, traction distributions over the potential crack interface, and threshold loads and initiating locations of delamination, are found to be in very good agreement with finite element simulations using cohesive elements. The modeling strategy, therefore, is useful for aerospace structural engineers at the preliminary design stage of laminated composite aerospace structures. Keywords: laminate, plate, analytical, Rayleigh-Ritz, FSDT, cohesive zone modeling, delamination threshold load, flexural response

∗ Corresponding

author Email addresses: [email protected] (Jiawen Xie), [email protected] (Anthony M. Waas), [email protected] (Mostafa Rassaian)

Preprint submitted to Composite Structures

July 19, 2017

1. Introduction During service, unexpected impacts on laminated composite structures by foreign objects, such as dropped tools, runway debris, and service vehicles, can occur. Damage induced by those impact events is a major concern in designing laminated composite structures because it can result in barely visible damage while altering the functionality of a designed structure in terms of stiffness and strength [1]. Delaminations or interfacial debonding is one common failure mode of impact damage to composite laminates. For preliminary design and sizing, it is beneficial to identify the following: the weakest interface of a laminate subjected to flexural loading; a corresponding critical load level that initiates a delamination (or, namely, a delamination threshold load); extending directions of the delamination; and their dependencies with design parameters, such as geometry and stacking sequences. A simple approach [2], which analyzes growth of an existing delamination at a mid-plane in an axisymmetric plate using linear elastic fracture mechanics (LEFM), was proposed to quickly estimate the critical load to extend a delamination in a quasi-isotropic laminated plate; Pc2 =

8π 2 Eh3 GIIc 8π 2 = Eeff h3 GIIc 2 9(1 − ν ) 9

(1)

where h is the laminate thickness, Eeff is the effective Young’s modulus of the plate, and GIIc is the mode II critical energy release rate. The merit of this approach is that the crack size term vanishes in the final equation, as shown in Eqn. (1), indicating that the critical load for delamination growth is independent of crack size. Therefore, the model is also considered to estimate a load level corresponding to an initial occurrence of damage. Though a 30% difference was found in predicting the delamination threshold load due to its simplicity [3, 4], a proportional relation of the critical load and scaling with h3/2 suggested by this equation have been found to be in good agreement with experimental measurements [5, 6]. The simple fracture mechanics approach has been extended to consider multiple delaminations [7, 8, 9] and to analyze high-velocity impact

2

events [10]. However, this simple approach cannot provide any evaluation of the stacking sequence effects because of the assumption of axisymmetry. Qualitative models have also been proposed to analytically explain a characteristic peanut shape of delaminations, with consideration of a mismatch in bending stiffness [11] or transverse displacement [12], by analyzing a stack of two plies with different material orientation. Compared to LEFM-based solutions and mismatch theories, cohesive zone modeling (CZM) considers additional interactions between adjacent plies at their interface [13, 14, 15] and apply interfacial constitutive laws that combines strength- and energy- based criteria to guide delamination initiation and propagation without prior knowledge of the crack location or growth path. Following Kanninen’s original idea [16], analytical CZM solutions consider the laminate as a stack of two sub-laminates connected by virtual deformable springs in normal and tangential directions. The solutions have been developed for many delamination characterization tests [17, 18, 19, 20, 21, 22, 23]. As a result, a critical failure load, both pre-peak and post-peak flexural responses, and a process zone length are obtained, on which parametric studies of interfacial properties, geometry and loading conditions can be performed. However, analytical CZM solutions presented in past literature are limited to beam configurations with some simplifications on stacking sequences and delamination locations. Analytical CZM solutions are extended in this paper to analyze delamination evolution in quasi-static indentation tests on laminated composite plates with arbitrary stacking sequences. Considering possible extension-shear couplings, the Rayleigh-Ritz method is used to find approximate solutions for the pre-peak flexural response and the delamination threshold load. The RayleighRitz approximation derived from first-order shear deformation theory (FSDT) for an elastic bending problem of a simply supported anisotropic plate, is considered in Section 2. The approach is then applied to analyze an assembly of two anisotropic sub-laminates with an elastic-brittle cohesive layer in-between, shown in Section 2. Elastic stiffness of the pre-peak response obtained by these two models are cross-checked against each other. Results of elastic stiffness, 3

traction distributions over potential crack interfaces, threshold loads and initiating locations of delaminations are further compared against those of finiteelement (FE) simulations using shell elements and cohesive elements. In this study, shear driven delamination, or mode II/III interfacial fracture, is the only failure mode considered motivated by its dominance in laminate flexure. Other occurrences of delamination that may result in mode I opening-crack failure or a state of mixed modes, such as matrix cracking induced delamination [24, 25, 26] or buckling induced delamination growth [27, 28], nor other intra-lamina failure modes [29, 30], are included. Note that failure progression in composite laminates is complicated as a consequence of competition among many different possible failure mechanisms. To determine a dominant failure mode, a strategy to respectively determine a critical load for each mode is necessary. Following that strategy, an analytical study on the delamination threshold load for shear driven delamination is presented.

2. Elastic Bending of a Laminated Plate

q(x, y) SS

z SS y h

SS b x SS a

Figure 1: A simply supported laminated plate under face-on (transverse) loading .

Consider an n-layer laminated plate with simply supported length, a, width, b, and thickness, h, subject to transverse pressure distribution q on its top surface, as shown in Figure 1. The displacement field of FSDT is of the form

4

[31],    U (x, y, z) = u(x, y) + zφx    V (x, y, z) = v(x, y) + zφy      W (x, y, z) = w(x, y)

(2)

where (U , V , W ) are displacement components along the (x, y, z) coordinate axes, respectively, of a point in the laminated plate; (u, v, w) are displacements of a point on the mid-plane (z = 0); φx (x, y) and φy (x, y) describe the rotations of a transverse normal about the y- and x-axes, respectively. The constitutive relations between through-the-thickness force resultants (N , M , Q) and midplane strains (0 , χ, γ) are,      ∂u         Nx  A A12 A16 B11 B12 B16    ∂x     11            ∂ v           N A A A B B B y 12 22 26 12 22 26     ∂y                ∂v ∂u   N  A   A26 A66 B16 B26 B66 + ∂y  xy 16 ∂x  =   ∂ φx         M B B B D D D x 11 12 16 11 12 16     ∂x                ∂ φ y           M B B B D D D y  12 22 26 12 22 26    ∂y              ∂ φy ∂ φx  M    B16 B26 B66 D16 D26 D66 + xy ∂x ∂y      N  A B 0   ⇔ = M  B D χ      ∂w  ts ts  Q  φ K + K x x 11 12  ∂x = ∂w  ts ts  Q  φ + K K y

12

y

22

(3a)

∂y

ts

⇔Q=K γ

(3b)

where A, B and D are the axial, bending-axial coupling, and bending stiffness matrices, respectively. The stiffness matrices are determined by lamina plane¯ k (k = 1, 2, . . . , n), stress stiffness matrix Q h A B

i

D =

n Z X

k=1

zk+1

zk

h ¯k 1 z Q

i z 2 dz

(z1 = −h/2, zk is the z-coordinate of lower surface of the kth lamina, zn+1 = h/2) (4)

5

K ts is the transverse shear stiffness matrix, obtained by equating shear strain energy computed by transverse force resultants of laminate with an integration of shear strain energy density of lamina [32] (see Appendix A). Closed-form solutions to the problem described in Figure 1 are available for cross-ply laminates [31] (see Appendix B). However, for laminates containing off-axis angle plies, closed-form solutions are unavailable due to an existence of extension-shear couplings, represented as non-zero 16, 26 components of A, B, D. The Rayleigh-Ritz method is used to determine approximate solutions for laminates with arbitrary stacking sequences. Weak form of governing equations is obtained using the principle of virtual work. Admissible approximate functions of (u, v, w, φx , φy ) that satisfy simply supported boundary conditions are chosen as double sine series with unknown coefficients. The virtual displacements and rotations are assumed to be approximated similarly. Since the the coefficient of each double sine term of virtual displacements and rotations can be arbitrary chosen, the weak form equation becomes a system of equations by which the unknown coefficients of approximate functions can be solved. Thus, the load-displacement response of elastic bending of an anisotropic laminate can be obtained. The detail formulation is provided in Appendix C.

3. Bending of a Laminated Plate with a Potential Crack Interface Suppose it is of interest to study delamination evolution in one interface between two adjacent plies in a laminated plate during bending. The potential crack interface, parallel to the laminate plane, virtually divides the plate into an upper and a lower sub-laminate, with thickness of h1 and h2 , respectively. The subscript α = 1 denotes the upper sub-laminate, while α = 2 denotes the lower one. Following the concept of cohesive zone modeling, the laminated plate is considered as a stack of two sub-laminates that are connected by a zero-thickness virtual deformable layer, referred to as cohesive layer, at their interface, as shown in Figure 2. As a result, the cohesive layer introduces three pairs of equal and opposite traction distributions, σz , τx and τy , to the bottom

6

Figure 2: Modeling of a laminated plate with a potential crack interface.

surface of the upper sub-plate and the top surface of the lower sub-plate. The traction distributions, as well as damage initiation and crack propagation of the potential crack interface, is determined by separation displacements between those two surfaces, i.e., subject to traction-separation laws. Each sub-laminate is considered individually within a framework of FSDT. Displacement is defined in a local coordinate system located at the mid-plane of each sub-laminate, as shown in Figure 2,    U (α) (x, y, zα ) = u(α) (x, y) + zα φ(α)  x   (α) (α) (α) , V (x, y, zα ) = v (x, y) + zα φy      W (α) (x, y, z ) = w(α) (x, y) α

(α = 1, 2)

Laminate-level constitutive relations are      (α)  N (α)  A(α) B (α) 0   = M (α)  B (α) D (α)  χ(α)  Q(α) = K ts(α) γ (α) ,

(α = 1, 2)

(5)

(6a) (6b)

In this study, the upper and lower surface of the potential crack interface are assumed perfectly bonded in the transverse direction, w(x, y)(1) = w(x, y)(2) ≡ w(x, y) 7

(7)

resulting in an unknown reaction traction σz (x, y) at the cohesive surfaces, while relative deformations of the cohesive layer along tangential directions are allowed. Shear tractions are determined by tangent separation displacements,

τx =

∆U τ (∆S) , ∆S

τy =

∆V τ (∆S) ∆S

(8)

where the separation displacements along x- and y-axes are ∆U = U (1) (x, y, −h1 /2) − U (2) (x, y, h2 /2) = u(1) − u(2) −

h1 (1) h2 (2) φ − φx 2 x 2

(9a)

∆V = V (1) (x, y, −h1 /2) − V (2) (x, y, h2 /2) = v (1) − v (2) −

h1 (1) h2 (2) φ − φy 2 y 2

(9b)

and ∆S is their resultant ∆S =

p ∆U 2 + ∆V 2

(10)

The cohesive law and fracture parameters are assumed to be identical for mode II and mode III. The shear traction resultant τ is subject to a mode II/III cohesive law, chosen as a linear elastic-brittle law, as shown in Figure 3 or mathematically τ (∆S) =

  K1 ∆S 

0

∆S ≤ δc

(11)

∆S > δc

The parameters of the cohesive law satisfy GIIc =

1 τc δc , 2

K1 =

τc τ2 = c δc 2GIIc

(12)

A limitation of the suggested approach is that nonlinearity of the pre-failure response cannot be captured by the simple application of a linear cohesive law without softening segments. With one possible cause by accumulations of microscopic damage, the pre-failure nonlinearity of some materials can be significant in experimental measurements. However, previous studies on cohesive zone

8

τ τc K1 GIIc ∆S

0

δc

Figure 3: A linear elastic-brittle traction-separation law of mode II/III fracture.

modeling has shown that the shape of bi-linear cohesive laws has a limited influence on the critical failure load [22], indicating that the peak load is mainly determined by critical energy release rates. Considering that the objective of this study is to analytically predict delamination threshold loads and damage initiating locations of general laminates, the simplification to the cohesive laws appears acceptable, as well as beneficial from cost saving on formulation and computation time. Similar to the technique for solving the elastic bending of a laminate, the Rayleigh-Ritz method is used to determine approximate solutions to bending of the system containing two sub-laminates and one cohesive layer. The detailed formulation is provided in Appendix D. With the series coefficients solved, the displacement field of the two sub-laminates can be obtained under a given loading condition. Interfacial shear traction distribution can be further computed by Eqn. (9), (10), (11). Thus, the solutions can provide answers to the following concerns. 3.1. Initiation of Delaminations on Pristine Interfaces Before any delamination occurs, the interface is pristine. According to the definition of the linear elastic-brittle law, delaminations will initiate at locations where the traction is maximum over the interface when the maximum value reaches the critical strength value τc . The critical load corresponding to delamination initiation is recorded as the load when the maximum traction reaches the critical strength value. Repeating the analysis to consider other interfaces of interest that divide 9

the laminate into different pairs of sub-laminates, the minimum value of the obtained critical loads can be considered as the delamination threshold load of the laminates. 3.2. Propagation of Existing Delaminations It should be pointed out that the proposed method has a limited predictive capability on propagation of existing delaminations. Similar to delamination initiation, propagation is also determined by the interfacial shear traction distribution. However, stress concentration in a narrow zone from the edge of delamination introduces certain singularities and therefore, may require much more terms of series in computations than the case of pristine interfaces. Additionally, with the application of the elastic-brittle cohesive law, the abrupt jump in traction at the edge of delamination brings in a discontinuity. Gibbs phenomenon can be significant and will not die out as more terms are added. With those restrictions, it can be difficult to obtain an acceptable yet computationally efficient interfacial traction distribution for the problem with existing delaminations. Note that the Rayleigh-Ritz method finds approximate solutions by satisfying essential boundary conditions and a principle of minimum total potential energy of a system of two sub-laminates and one cohesive layer with existing delaminations. Though it can be difficult to obtain localized details as mentioned, the proposed method can provide information on system-level parameters. One application can be predictions of stiffness change due to delaminations.

4. Results and Discussions Numerical evaluations of the proposed method have been performed using Matlab to analyze configurations of quasi-static face-on indentation tests of laminated composite plates. As shown in Figure 4, a plate is simply supported on its four edges and loaded at the center of its top surface by a hemispherical indenter. Dimensions are given in Figure 4. Material properties and interfacial

10

fracture properties of IM7/8552 graphite/epoxy lamina [33] were used, as shown in Table 1. If the linear elastic-brittle law of mode II/III fracture is used, as shown in Figure 3, δc = 2.21 × 10−2 mm, K1 =3.17 GPa/mm. SS P, ∆ R=25.4 mm z SS y x

h=3 mm

SS

◦ 90+45 ◦ 0◦

b=101.6 mm

SS a=152.4 mm

Figure 4: The static face-on indentation test configuration.

Table 1: Homogenized lamina properties and interface fracture properties from Ref. [33]. For the interface strength σc and τc , values were assumed from similar materials in literature. E11 E22 = E33

161

GPa

11.38

GPa

ν12 = ν13

0.32

ν23

0.45

G12 = G13

5.2

GPa

G23

3.9

GPa

GIc

0.212

N/mm

GIIc = GIIIc

0.774

N/mm

σc

50

MPa

τc

70

MPa

A parabolic distribution was used to approximate transverse pressure applied by the indenter,    r2 2   2 1− 2 P , rc q(x, y) = πrc   0,

with total external force as F =

RbRa 0

0

r ≤ rc

(13)

r > rc

q(x, y)dxdy = P . rc is the radius of

a small but finite contact area between the indenter and the top surface of the plate, which is taken as rc = R/9 by measuring the contact area in FE 11

simulations, and r is the distance from the center of the plate, s 2   a 2 b r= x− + y− 2 2

(14)

The load coefficient Qpq , shown in Eqn. (C.6), can be computed by Gaussian quadrature in polar coordinates.

Figure 5: FE model of the quasi-static face-on indentation test of a laminated composite plate.

Accuracy of the proposed Rayleigh-Ritz method was first evaluated. Elastic bending of center loaded pristine plates with eight stacking sequences was considered. A different ply thickness was used for the different laminate cases so that the total plate thickness was kept 3 mm. FE simulations using Abaqus/ Standard [34] were performed as a reference. In the FE model, shown in Figure 5, the plate was modeled by conventional composite shell elements (S4) of size 0.508 mm (L) × 0.508 mm (W) while the indenter was modeled as an analytical rigid surface. Simply supported boundary conditions, similar to those of analytical solutions, were applied to the edges of the plate. During the simulations, the plate was loaded by the displacement-controlled indenter via a contact interaction between them. The comparison of the stiffness or the slope of load-displacement response, K = P/∆, is summarized in Table 2. It can be seen that the predictions made by the Rayleigh-Ritz method with 80×80 double sine series agree well with the results obtained by FE simulations within 2.5% error. The predictions are slightly better for plates with the symmetric stacking sequences than the asymmetric ones. For cross-ply laminates, the stiffness are found exactly same as those provided by closed-form solutions shown in Ap12

pendix B. Convergence studies of the stiffness on number of terms in the series considered in computation have been performed for laminates with stacking sequences (0) and (+45/-45/0/90)2 . As shown in Figure 6, convergence errors of the stiffness K between series, KM ×N − K(M −1)×(N −1) , econv = KM ×N

M = N = 20, 21, . . . , 80

(15)

for both laminates are less than 0.1% starting from 40×40 terms. Table 2: Stiffness predicted by the Rayleigh-Ritz method (80×80 series) and FE simulations. Stacking Sequences

Stiffness K [N/mm]

(from bottom to top)

FE

R.-R.

Error

(0)

382.34

386.96

1.21%

(0/90)s

490.01

496.69

1.36%

(+45/-45/0/90)s

938.13

956.44

1.95%

(-45/+45/90/0)s

987.14

1006.70

1.98%

(0/90)2

657.89

668.37

1.59%

(+45/-45/0/90)2

869.15

888.18

2.19%

(-45/+45/90/0)2

842.23

860.38

2.15%

(0/90/+45/-45/0/-45/+45/90/0)

732.24

743.56

1.55%

Delamination initiation has been studied by the proposed Rayleigh-Ritz approximation with cohesive zone modeling for two laminate cases: a unidirectional laminate (0)2 with a ply thickness of 1.5 mm, and a quasi-isotropic laminate (0/90/+45/-45/0/-45/+45/90/0) with a ply thickness of 0.33 mm. A corresponding FE model is shown in Figure 7. A single layer of cohesive elements (COH3D8) of a small thickness, 0.01h, that are properly tied to the two sub-laminates, was used to model the interface. A quadratic stress criterion for crack initiation 

hσz i σc

2

+



τx τsc

2

+



τy τtc

2

=1

(16)

with τsc = τtc = τc =70 MPa, and a linear power law of fracture energy for crack propagation GII GIII GI + + =1 GIc GIIc GIIIc

13

(17)

0.08 FEM FSDT Error

388

0.06

386

0.04

384

0.02

382 0 20x20 30x30 40x40 50x50 60x60 70x70 80x80

Convergence Error between Series [%]

Stiffness P/∆ [N/mm]

390

(0) (IM7/8552) ctr. dist. loading Convergence Study of Stiffness

Number of Series MxN (+45/-45/0/90) (IM7/8552) ctr. dist. loading

Stiffness P/∆ [N/mm]

910

0.2 FEM FSDT Error

900

0.15

890 0.1 880 0.05

870

860 0 20x20 30x30 40x40 50x50 60x60 70x70 80x80

Convergence Error between Series [%]

2

Convergence Study of Stiffness

Number of Series MxN

Figure 6: Convergence studies of the stiffness on number of series considered in the RayleighRitz method.

were assigned in FE analyses. Initial elastic stiffness of cohesive elements were σc2 τc2 set as KI = = 5.90 GPa/mm, KII = KIII = = 3.17 GPa/mm, 2GIc 2GIIc following linear elastic-brittle laws. The FE model developed here was to verify the analytical solutions in a similar setting. The stiffness values used here are smaller to those commonly used in bi-linear quasi-brittle laws (Ref. [35] as an example). The delamination was assumed initiated at locations of the interface that first reaches the quadratic stress criterion while the load was also recorded as the critical load. The flexural stiffness and the critical loads for the two laminates are summarized in Table 3 and Table 4, respectively. Good agreement is found between 14

Figure 7: FE model with a cohesive layer.

analytical solutions and FE simulations, especially for the stiffness. The stiffness also agree well with the values shown in Table 2, indicating that the overall stiffness of the plate will not be significantly affected by inserting a thin and relatively soft interfacial layer for both analytical solutions and FE analyses. The analytical solutions slightly over-estimate the critical loads yet within an average error of just 8%. The delamination threshold load of laminates is predicted as 2.553 kN and 2.985 kN for the uni-directional laminate (0/0) and the quasi-isotropic laminate (0/90/+45/-45/0/-45/+45/90/0), respectively, with an over-estimation of 5% compared with the value found in FE simulations. The predicted delamination threshold loads were further compared with the critical threshold load given by the simple fracture mechanics model [2] (dEqn. (1)). For specially orthotropic plates, the effective plate stiffness D∗ can be approximated as [9], ∗

D =

s

1 D11 D22 2



D12 + 2D66 √ +1 D11 D22



(18)

Thus, the effective Young’s modulus of the uni-directional laminates is Eeff = 12D∗ = 35.1 GPa. The simple fracture mechanics model provides a good h3 15

estimation on the threshold load Pc = 2.537 kN, with an error of 4.5% compared with the value obtained by the FE model. For the quasi-isotropic laminate, the effective Young’s modulus was obtained as a range, Eeff = 63.2 ∼ 98.1 12D11 12D22 GPa, since it can be estimated by: Eeff = or . The threshold h3 h3 load Pc = 3.403 ∼ 4.241 kN, given by Eqn. (1), has at least a 20% overprediction above the FE analyses results. Overall, the proposed method, with an over-prediction of 5%, improves the analytical predicability on the delamination threshold load for laminates of arbitrary stacking sequences. Moreover, the ranking of weakness of the interfaces is captured correctly. As shown in Table 4, the fourth (-45/0) and fifth (0/-45) interface of the laminate with stacking sequence (0/90/+45/-45/0/-45/+45/90/0) are the weakest ones. Because of the assumption of perfectly bonded in the transverse direction between two sub-laminates and the nature of FSDT, the proposed method does not distinguish the relative vertical position of the two sub-laminates, nor the through-the-thickness location where transverse loading and interfacial traction are applied. In other words, the analytical solutions provide exactly the same answers when a configuration, with an order of lower-sublaminate, interface and upper-sublaminate from bottom to top, is flipped. Therefore, for symmetric laminates, the flexural stiffness and the delamination threshold load are the same for the interfaces that are symmetric about the mid-plane, reflected in the analytical results shown in Table 4. Table 3: Stiffness and critical loads predicted by the Rayleigh-Ritz method (80×80 series) and FE simulations of laminate with stacking sequence (0)2 . Interface (from bottom to top) #1: 0/0

Stiffness K [N/mm]

Critical Load Pc [kN]

FE

R.-R.

Error

FE

R.-R.

Error

377.94

376.92

0.27%

2.428

2.553

5.15%

Convergence studies of the stiffness and the delamination threshold load on the number of terms in the double sine series have been conducted for the uni-directional laminate. As shown in Figure 8,the results are converged after 60×60 terms with both convergence errors of the stiffness and the delamination

16

Table 4: Stiffness and critical loads predicted by the Rayleigh-Ritz method (80×80 series) and FE simulations of laminate with stacking sequence (0/90/+45/-45/0/-45/+45/90/0). Interface

Stiffness K [N/mm]

Critical Load Pc [kN]

(from bottom to top)

FE

FE

R.-R.

Error

R.-R.

Error

#1: 0/90

735.44

737.66

0.30%

5.742

6.247

8.79%

#2: 90/+45

725.17

728.97

0.52%

4.289

4.590

7.03%

#3: +45/-45

716.03

722.19

0.86%

3.115

3.378

8.45%

#4: -45/0

713.42

719.41

0.84%

2.825

2.985

5.66%

#5: 0/-45

713.42

719.41

0.83%

2.848

2.985

4.79%

#6: -45/+45

715.85

722.19

0.89%

3.004

3.378

12.46%

#7: +45/90

725.35

728.97

0.50%

4.298

4.590

6.80%

#8: 90/0

736.77

737.66

0.12%

6.652

6.247

6.09%

threshold load less than 1%. Comparisons of the normalized interfacial shear traction distribution computed by the Rayleigh-Ritz approximations and FE simulations are shown in Figure 9 for the uni-directional laminates, and Figure 10, Figure 11 for the quasi-isotropic one. The initiating locations of delamination of each interface can be directly seen in these plots as red concentration areas, as an early form of peanut shaped crack patterns. Good agreements between the analytical solutions and the numerical simulations can be found. As mentioned earlier, the traction distribution of the analytical solutions are the same for a pair of interfaces symmetric about the mid-plane of laminates with stacking sequence (0/90/+45/-45/0/-45/+45/90/0). The symmetry can be observed in the results of the FE simulations also.

5. Conclusions Analytical CZM solutions have been developed to determine the pre-peak flexural response and the delamination threshold load in a simply supported laminated composite plate subject to transverse loading, without any restriction on laminate stacking sequences and locations of a potential or existing delamination. The laminated plate is modeled as two anisotropic sub-laminates

17

0.1

FEM FSDT Error

380

0.08

379

0.06

378

0.04

377

0.02

376 0 20x20 30x30 40x40 50x50 60x60 70x70 80x80

Convergence Error between Series [%]

Stiffness P/∆ [N/mm]

381

(0/0) (IM7/8552) ctr. dist. loading Convergence Study of Stiffness

(0/0) (IM7/8552) ctr. dist. loading Convergence Study of Threshold Load of Interface #1 6

FEM FSDT Error

4

c

Threshold Load P [kN]

4.5

5 4

3.5 3 3 2 2.5 2 20x20

1

30x30

40x40

50x50

60x60

70x70

0 80x80

Convergence Error between Series [%]

Number of Series MxN

Number of Series MxN

Figure 8: Convergence studies of the stiffness and the delamination threshold load on number of series considered in the Rayleigh-Ritz method.

Figure 9: Normalized interfacial shear traction distribution over the interface of laminate with stacking sequence (0)2 by the Rayleigh-Ritz method (left) and FE simulations (right).

that are perfectly bonded in the transverse direction while separable in the tangential directions in which the cohesive interactions are subject to a linear elastic-brittle traction-separation law. In this sense, only mode II/III driven delaminations are considered. The problem is formulated within the framework of first order shear deformation theory, and further solved by the Rayleigh-Ritz

18

approximations that satisfy kinematic boundary conditions and the principle of minimum total potential energy of the system. A critical load corresponding to delamination initiation on each interface of a laminate can be obtained as a load level when a maximum traction on that interface reaches a critical strength value. The delamination threshold load of the laminate is considered as the minimum value of all critical loads on all possible interfaces. The delamination initiating locations are where the traction is maximum over the weakest interface. Convergence of the results shown is ensured by the results of the convergence studies. Results of flexural stiffness of the pre-peak response obtained by the models without and with a cohesive layer agree well with the values obtained in FE simulations. Though the delamination threshold loads are on average overpredicted by 8% above that of FE simulations, the proposed method makes a considerable improvement compared with the simple fracture mechanics model, in terms of accuracy, and capability of identifying the weakest interface among all interfaces and delamination initiating locations. The proposed method is easy to implement using Matlab or other computer programs with all equations given and solution techniques as described in this paper. With calculations done in minutes, the proposed method can be a design tool to preform parametric studies and explore optimization possibilities on laminate stacking sequences and geometry. The proposed analytical solutions can be used with confidence to predict delamination evolution in laminated composite plates and other general multilayered structures. The proposed method can also be formulated following classical lamination theory (CLT) (see Appendix E). The method is extendable to considering more complicated cohesive laws, multiple interfaces in a single model, in-plane loading conditions, and other boundary conditions. For some of the boundary conditions, admissible displacement functions are provided in Ref. [31].

19

Acknowledgments The authors are grateful for financial support from the Boeing company.

Appendix A. Transverse Shear Stiffness of Laminated Plates Following Ref. [32], consider only bending and shear of a laminated plate with an assumption of cylindrical bending along the width direction. Therefore, force resultants are zero except Mx and Qx , and all derivatives with respect to y are zero. Assume the reference plane associated with pure bending in the x-direction, denoted as z = zx0 , is not necessarily the mid-plane (z = z0 = 0). Thus, in-plane strain,  = 0 + (z − zx0 ) χ

(A.1)

where mid-plane strain can be represented by force resultants as the inverse of the laminate constitutive relation in Eqn. (3a)   0  χ



A

= B

  −1   N  N   ≡H M  M  D B

(A.2)

where H is the laminate compliance matrix. Substituting Eqn. (A.1) and (A.2) into lamina stress-strain relation ¯ k σ=Q

(A.3)

the normal stress in x-direction for the kth lamina can be written as

where


k k Mx σxk = Bx1 + (z − zx0 ) Bx2 k Bx1 = Qk11 H14 + Qk12 H24 + Qk16 H34 k Bx2 = Qk11 H44 + Qk12 H54 + Qk16 H64

(A.4)

(A.5)

Further combined Eqn. (A.4) with a stress equilibrium in the x-direction ∂ σxk ∂τk + xz = 0 ∂x ∂z 20

(A.6)

and a moment equilibrium about the y-axis Qx =

∂ Mx ∂x

(A.7)

one will get a relation between laminate shear force resultant and lamina transverse shear stress k
∂ τxz k k = − Bx1 + (z − zx0 ) Bx2 Qx ∂z

(A.8)

Integrating Eqn. (A.8) through the thickness of the laminate, and using boundary conditions and continuities at z = z1 = −h/2 :

1 τxz =0

n =0 at z = zn+1 = h/2 : τxz

(A.9)

k k+1 τxz = τxz

at z = zk :

gives the transverse shear stress for the kth lamina   1 2 k k k k ˜ τxz = − (z − zk ) Bx2 + (z − zk ) Bx1 + Bx0 Qx 2

(A.10)

where k k k ˜x1 = Bx1 + (zk − zx0 ) Bx2 B  k−1 X 1 k i i ˜ ti Bx2 + Bx1 Bx0 = ti 2 i=1
PN k k tk Bx1 + 12 (zk+1 + zk ) Bx2 zx0 = k=1
PN k k=1 tk Bx2 A11 H14 + A12 H24 + A16 H34 + B11 H44 + B12 H54 + Q16 H64 = =0 A11 H44 + A12 H54 + A16 H64 (A.11)

and tk is the thickness of the kth lamina. Eqn. (A.11) indicates that the reference plane associated with pure bending in the x-direction is the reference plane chosen for formulation, which is the mid-plane in this paper 1 . 1 It

can also be proven that zx0 = z0 is true even if the reference plane (z = z0 6= 0) is

offset from the midplane.

21

Similarly, assuming pure bending in the y-direction gives   1 2 k k k k ˜y1 =− τyz (z − zk ) By2 + (z − zk ) B + By0 Qx 2

(A.12)

Consider that shear strain energy computed by transverse force resultants of laminate is equivalent to an integration of shear strain energy density of lamina,     n Z zk+1 n Q  1 X τ  o o 1n x xz ts ¯k = (A.13) Qx Qy F τxz τyz S Q  2 τ  2 zk y

yz

k=1

where lamina transverse shear compliance matrices are  ¯k Q ¯ k =  44 S ¯k Q 45

¯k Q 45 ¯k Q 55

−1 

,

k = 1, 2, . . . , n

(A.14)

¯ 55 = G23 , Q ¯ 45 = 0. Substituting Eqn. (A.8) and (A.12) ¯ 44 = G12 , Q For 0-ply, Q to Eqn. (A.13) and performing a through-the-thickness integration gives the transverse shear compliance component of laminate " n X
ts k k 2 k ˜k + tk Bx0 S¯44 tk Bx0 Bx1 F11 = k=1

ts F22

=

n X

k=1

ts F12

k tk S¯55

"

#  
1 3 ˜k k 1 2  ˜ k 2 1 4 k k k 2 Bx1 + Bx0 Bx2 + tk Bx1 Bx2 + tk Bx2 + tk 3 4 2

k By0


2

k ˜k + tk By0 By1

#  
1 2  ˜ k 2 1 4 1 3 ˜k k k k k 2 + tk By1 + By0 By2 + tk By1 By2 + tk By2 3 4 2 "  n  1 X 1  k ˜k k k k k k k ˜k ˜x1 ˜x1 = tk Bx0 By0 + tk Bx0 By0 S¯45 By1 + B By1 + t2k B 2 3 k=1 #   1
1 1 3  ˜k k k k k k k ˜k 4 k k + + tk Bx1 By2 + Bx2 By1 + tk Bx2 By2 Bx0 By2 + Bx2 By0 2 8 2 (A.15)

Transverse shear stiffness can be further obtained as  −1 K ts = F ts 22

(A.16)

For a specially orthotropic (or homogeneous orthotropic) plate, Bx0 = By0 = 0. The plate is naturally symmetric about its mid-plane, B = 0, and the laminate compliance matrix is H=



1 −1 Q h

resulting in Bx1 = By1 = 0, Bx2 = By2 =

12 −1 h3 Q 12 h3 ,

ts ts ts K11 = 56 G12 h, K22 = 56 G23 h, K12 = 0, where



(A.17)



˜x1 = B ˜y1 = − 62 . Therefore, B h

5 6

is the shear correction factor

known for orthotropic plates.

Appendix B. Closed-form Solutions For Bending of Cross-ply Laminates For cross-ply laminates, ¯ k16 = Q ¯ k26 = Q ¯ k45 = 0 (k = 1, 2, . . . , n) Q ⇒ A16 = A26 = 0 ,

B16 = B26 = 0 ,

D16 = D26 = 0 ,

ts K12 =0

(B.1)

Governing equations are  ∂ Nx ∂ Nxy   +   ∂x ∂y     ∂ N ∂ Ny  xy  +    ∂x ∂y    ∂Q ∂ Qy x +  ∂x ∂y     ∂ M ∂ M x xy   +   ∂x ∂y       ∂ Mxy + ∂ My  ∂x ∂y

=0 =0 = q(x, y)

(B.2)

= Qx = Qy

Substituting Eqn. (3), (C.2) and (B.1) into the governing equations, one will get five equations for each combination of p = 1, 2, . . . , M and q = 1, 2, . . . , N :  h  pπ 2  qπ 2  pπ qπ i + A66 bpq A11 apq + (A12 + A66 ) a b a b    pπ 2  qπ 2 h pπ qπ i dpq + (B12 + B66 ) + B11 + B66 epq = 0 (B.3a) a b a b 23

  pπ 2  qπ 2  pπ qπ i + A22 bpq apq + A66 a b a b   pπ 2  qπ 2  h pπ qπ i dpq + B66 epq = 0 + (B12 + B66 ) + B22 (B.3b) a b a b   2  2  4 ts pπ ts qπ ts pπ ts qπ + K22 dpq + K22 epq = − Qpq (B.3c) cpq + K11 K11 a b a b ab  h  pπ 2  qπ 2  pπ qπ i ts pπ + B66 B11 bpq + K11 cpq apq + (B12 + B66 ) a b a b a   h  pπ 2  qπ 2 pπ qπ i ts + D11 + D66 + K11 epq = 0 dpq + (D12 + D66 ) a b a b h

h

(A12 + A66 )

(B.3d)   pπ 2  qπ 2  ts qπ apq + B66 cpq bpq + K22 (B12 + B66 ) + B22 a b a b b    pπ 2  qπ 2 h pπ qπ i ts + (D12 + D66 ) dpq + D66 + D22 + K22 epq = 0 a b a b pπ qπ i

(B.3e)

The equations above are de-coupled between different terms in double sine series. For each combination of p and q, five coefficients, apq , bpq , cpq , dpq , epq can be solved by five equations directly using Eqn. (B.3).

Appendix C. The Rayleigh-Ritz Approximate Solutions For Bending of Laminates Weak form of governing equations can be obtained by the principle of virtual work, b

a

"

    ∂ δv ∂ δv ∂ δw ∂ δw + Ny + Qy + qδw + Nxy + Qx ∂x ∂y ∂x ∂y 0 0    # ∂ δφx ∂ δφx ∂ δφy ∂ δφy + Mxy + Qx δφx + Mxy + My + Qy δφy + Mx dxdy = 0 ∂x ∂y ∂x ∂y

Z

Z

∂ δu ∂ δu + Nxy Nx ∂x ∂y



(C.1)

24

Admissible approximation functions that satisfy simply supported boundary conditions are double sine series      M X N h i i X h  iπx jπy   cos = sin aij dij  u(x, y) φx (x, y)  a b   i=1 j=1         N h M X i i X h iπx jπy sin = cos v(x, y) φy (x, y) bij eij a b   i=1 j=1         N M X  X  jπy iπx   sin w(x, y) = cij sin   a b i=1 j=1

(C.2)

where aij , bij , cij , dij , eij , (i = 1, 2, . . . , M , j = 1, 2, . . . , N ) are coefficients to be solved; M and N are numbers of terms related to x and y, respectively. Virtual displacements and rotations (δu, δv, δw, δφx , δφy ) are assumed to be approximated similarly to Eqn. (C.2) with virtual coefficients δapq , δbpq , δcpq , δdpq , δepq , (p = 1, 2, . . . , M and q = 1, 2, . . . , N ). Substituting Eqn. (3) and approximate functions into the weak form, Eqn. (C.1) becomes a form of " # M X N X (. . . ) δapq + (. . . ) δbpq + (. . . ) δcpq + (. . . ) δdpq + (. . . ) δepq = 0 p=1 q=1

(C.3)

Since the virtual coefficients can be arbitrarily chosen, expressions inside each parenthesis need to be zero to satisfy Eqn. (C.3). Therefore, one will get five equations for each combination of p and q, ( M,N X  iπ pπ ip jπ pπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − A16 (b) A11 a a b a i,j

 iπ qπ ip jπ qπ ip jq jq −A16 I (a)Isc (b) + A66 I (a)Icc (b) aij a b sc b b cc  jπ pπ ip iπ pπ ip jq jq + A12 I (a)Iss (b) − A16 I (a)Ics (b) b a ss a a cs  jπ qπ ip iπ qπ ip jq jq −A26 I (a)Isc (b) + A66 I (a)Icc (b) bij b b sc a b cc  iπ pπ ip jπ pπ ip jq jq I (a)Iss I (a)Ics + B11 (b) − B16 (b) a a ss b a cs  iπ qπ ip jπ qπ ip jq jq I (a)Isc (b) + B66 I (a)Icc (b) dij −B16 a b sc b b cc 25

 jπ pπ ip iπ pπ ip jq jq + B12 I (a)Iss (b) + B16 I (a)Ics (b) b a ss a a cs

 ) jπ qπ ip iπ qπ ip jq jq I (a)Isc (b) + B66 I (a)Icc (b) eij = 0 −B26 b b sc a b cc (C.4a)

( M,N X  i,j

− A16

iπ pπ ip jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + A66 (b) a a sc b a cc

 iπ qπ ip jπ qπ ip jq jq (b) − A26 (b) aij Iss (a)Iss Ics (a)Ics a b b b  jπ pπ ip iπ pπ ip jq jq I (a)Isc I (a)Icc + − A26 (b) + A66 (b) b a sc a a cc  jπ qπ ip iπ qπ ip jq jq (b) − A26 (b) bij +A22 Iss (a)Iss Ics (a)Ics b b a b  iπ pπ ip jπ pπ ip jq jq + − B16 I (a)Isc (b) + B66 I (a)Icc (b) a a sc b a cc  iπ qπ ip jπ qπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − B26 (b) dij +B12 a b b b  jπ pπ ip iπ pπ ip jq jq I (a)Isc I (a)Icc + − B26 (b) + B66 (b) b a sc a a cc  ) jπ qπ ip iπ qπ ip jq jq +B22 I (a)Iss (b) − B26 I (a)Ics (b) eij = 0 b b ss a b cs +A12

(C.4b) ( M,N X  ts jπ pπ ip jq ts iπ pπ ip jq K12 Isc (a)Ics (b) + K11 Icc (a)Iss (b) b a a a i,j

 jπ qπ ip jq ts iπ qπ ip jq (b) + K12 (b) cij Iss (a)Icc Ics (a)Isc b b a b   ts pπ ip jq ts qπ ip jq + K11 Icc (a)Iss (b) + K12 Ics (a)Isc (b) dij a b   ) ts pπ ip jq ts qπ ip jq + K12 Isc (a)Ics (b) + K22 Iss (a)Icc (b) eij = −Qpq a b

ts +K22

(C.4c) ( M,N X  iπ pπ ip jπ pπ ip jq jq Iss (a)Iss Ics (a)Ics B11 (b) − B16 (b) a a b a i,j −B16

 iπ qπ ip jπ qπ ip jq jq Isc (a)Isc Icc (a)Icc (b) + B66 (b) aij a b b b 26

 jπ pπ ip iπ pπ ip jq jq + B12 I (a)Iss (b) − B16 I (a)Ics (b) b a ss a a cs  jπ qπ ip iπ qπ ip jq jq Isc (a)Isc Icc (a)Icc (b) + B66 (b) bij −B26 b b a b   ts jπ ip jq ts iπ ip jq Isc (a)Ics Icc (a)Iss + K12 (b) + K11 (b) cij b a  iπ pπ ip jπ pπ ip jq jq + D11 (b) − D16 (b) I (a)Iss I (a)Ics a a ss b a cs iπ qπ ip jπ qπ ip jq jq −D16 Isc (a)Isc I (a)Icc (b) + D66 (b) a b b b cc 

ts ip jq +K11 Icc (a)Iss (b) dij

 jπ pπ ip iπ pπ ip jq jq I (a)Iss I (a)Ics + D12 (b) − D16 (b) b a ss a a cs jπ qπ ip iπ qπ ip jq jq −D26 Isc (a)Isc I (a)Icc (b) + D66 (b) b b a b cc ts ip jq Isc (a)Ics (b) +K12

 ) eij = 0 (C.4d)

( X  ij

− B16

iπ pπ ip jπ pπ ip jq jq Isc (a)Isc I (a)Icc (b) + B66 (b) a a b a cc

 iπ qπ ip jπ qπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − B26 (b) aij a b b b  jπ pπ ip iπ pπ jq jq + − B26 (b) + B66 (b) I (a)Isc I ip (a)Icc b a sc a a cc  jπ qπ ip iπ qπ ip jq jq +B22 I (a)Iss (b) − B26 I (a)Ics (b) bij b b ss a b cs   ts jπ ip jq ts iπ ip jq + K22 Iss (a)Icc (b) + K12 Ics (a)Isc (b) cij b a  pπ iπ pπ ip jπ jq jq + − D16 I (a)Isc (b) + D66 I ip (a)Icc (b) a a sc b a cc iπ qπ ip jπ qπ ip jq jq +D12 I (a)Iss (b) − D26 I (a)Ics (b) a b ss b b cs  +B12

ts ip jq +K12 Ics (a)Isc (b) dij

 jπ pπ ip iπ pπ ip jq jq Isc (a)Isc I (a)Icc + − D26 (b) + D66 (b) b a a a cc jπ qπ ip iπ qπ ip jq jq (b) − D26 (b) +D22 I (a)Iss I (a)Ics b b ss a b cs 27

ts ip jq +K22 Iss (a)Icc (b)

 ) eij = 0 (C.4e)

where the integration terms are mn (l) = Iss

mn Icc (l) =

mn Ics (l)

=

Z

l

Z

l

Z

l

sin

0

cos

0

cos

0



mπξ l



sin



mπξ l



cos



mπξ l



sin



nπξ l



dξ =



nπξ l



dξ =



nπξ l



dξ =

0

0

a

q(x, y) sin

m=n m 6= n

 l/2 m = n   0 m + n even

 pπx  a

m 6= n

 l/2   0

 

mn nm Isc (l) = Ics (l)

and the load coefficient is Z bZ Qpq =

  0

2n l n2 −m2 π

sin

 qπy  b

(C.5)

m + n odd

dxdy

(C.6)

Therefore, the coefficients, aij , bij , cij , dij , eij , with a total number of 5×M ×N , can be solved by assembling 5 × M × N equations. Appendix D. The Rayleigh-Ritz Approximate Solutions For Bending of Laminates with a Cohesive Layer The weak form can be derived from the principle of virtual work of the system of two sub-laminates and one cohesive layer,    Z b Z a ( (1) (2) ∂ δu(1) ∂ δu(2) (1) ∂ δu (2) ∂ δu Nx(1) + Nxy + τx δu(1) + Nx(2) + Nxy − τx δu(2) ∂x ∂y ∂x ∂y 0 0     (1) (1) (2) ∂ δv ∂ δv (2) (1) ∂ δv (2) ∂ δv + Ny(1) + τy δv (1) + Nxy + Ny(2) − τy δv (2) + Nxy ∂x ∂y ∂x ∂y    ∂ δw   ∂ δw (2) (2) + Q(1) + qδw + Q(1) x + Qx y + Qy ∂x ∂y " #   (1) (1) h1 (1) ∂ δφx (1) ∂ δφx (1) (1) + Mx + Mxy + Qx − τx δφx ∂x ∂y 2 28

+

"

(2) ∂ δφx Mx(2)

∂x

"

∂ δφy ∂x

"

∂ δφy ∂x

+

(2) (2) ∂ δφx Mxy

∂y

(1)

(1) + Mxy

(1)

+ My(1)

∂ δφy ∂y

+ My(2)

∂ δφy ∂y

(2)

(2) + Mxy

(2)

#  h2 (2) + − τx δφx 2 #   h 1 + Q(1) τy δφ(1) y − y 2 #)   h2 (2) (2) dxdy = 0 + Qy − τy δφy 2 

Q(2) x

(D.1) (α)

(α)

Substituting admissible approximation functions of u(α) , v (α) , w, φx , φx

and

virtual terms according to Eqn. (C.2) for the two sub-laminates α = 1, 2 and Eqn. (6) into Eqn. (D.1), one will get " N M X X (2) (1) (2) (. . . ) δa(1) pq + (. . . ) δapq + (. . . ) δbpq + (. . . ) δbpq + (. . . ) δcpq p=1 q=1

(. . . ) δd(1) pq

+

+

(. . . ) δd(2) pq

+

(. . . ) δe(1) pq

+

(. . . ) δe(2) pq

#

=0

and since the virtual coefficients are arbitrarily chosen, Eqn. (D.2) gives nine equations for each combination of p and q: (" M,N X (1) iπ pπ ip (1) jπ pπ ip jq jq A11 Iss (a)Iss (b) − A16 Ics (a)Ics (b) a a b a i,j

qπ ip (1) jπ qπ ip jq jq Isc (a)Isc (b) + A66 Icc (a)Icc (b) a b b b # " # (1) (2) ip jq ip jq ˜ ˜ ˜ ˜ +K1 I (a)I (b) a + − K1 I (a)I (b) a

(1) iπ

−A16

cc

"

ss

ij

cc

ss

ij

pπ ip (1) iπ pπ ip jq jq I (a)Iss (b) − A16 I (a)Ics (b) b a ss a a cs

(1) jπ

+ A12

# qπ ip (1) iπ qπ ip (1) jq jq I (a)Isc (b) + A66 I (a)Icc (b) bij b b sc a b cc

(1) jπ −A26

"

pπ ip (1) jπ pπ ip jq jq I (a)Iss I (a)Ics (b) − B16 (b) a a ss b a cs

(1) iπ

+ B11

qπ ip (1) jπ qπ ip jq jq I (a)Isc Icc (a)Icc (b) + B66 (b) a b sc b b # " # h1 ˜ip ˜jq h2 ˜ip ˜jq (1) (2) −K1 Icc (a)Iss (b) dij + − K1 Icc (a)Iss (b) dij 2 2 (1) iπ

−B16

29

"

(1) jπ

+ B12

pπ ip (1) iπ pπ ip jq jq (b) − B16 (b) Iss (a)Iss I (a)Ics b a a a cs ) # qπ ip qπ iπ (1) (1) jq jq =0 I (a)Isc (b) + B66 I ip (a)Icc (b) eij b b sc a b cc

(1) jπ −B26

(D.2a) (" M,N X

#

"

(1) ip jq ip jq − K1 I˜cc (a)I˜ss (b) aij + K1 I˜cc (a)I˜ss (b)

i,j

(2) iπ

pπ ip (2) jπ pπ ip jq jq (b) − A16 (b) I (a)Iss I (a)Ics a a ss b a cs # qπ jπ (2) iπ qπ ip (2) (2) jq jq −A16 I (a)Isc I ip (a)Icc (b) + A66 (b) aij a b sc b b cc

+A11

"

(2) jπ

+ A12

pπ ip (2) iπ pπ ip jq jq I (a)Iss I (a)Ics (b) − A16 (b) b a ss a a cs

# iπ qπ ip qπ (2) (2) jq jq (b) + A66 (b) bij I (a)Isc I ip (a)Icc b b sc a b cc " # " h1 ˜ip ˜jq h2 ip ˜jq (1) + K1 Icc (a)Iss (b) dij + K1 I˜cc (a)Iss (b) 2 2 (2) jπ −A26

pπ ip (2) jπ pπ ip jq jq I (a)Iss (b) − B16 I (a)Ics (b) a a ss b a cs # qπ jπ (2) iπ qπ ip (2) (2) jq jq I (a)Isc I ip (a)Icc (b) + B66 (b) dij −B16 a b sc b b cc (2) iπ

+B11

"

pπ ip (2) iπ pπ ip jq jq Iss (a)Iss I (a)Ics (b) − B16 (b) b a a a cs

(2) jπ

+ B12

) # qπ ip (2) iπ qπ ip (2) jq jq =0 I (a)Isc (b) + B66 I (a)Icc (b) eij b b sc a b cc

(2) jπ −B26

(D.2b) (" M,N X i,j

(1) iπ

− A16

pπ ip (1) jπ pπ ip jq jq (b) + A66 (b) I (a)Isc I (a)Icc a a sc b a cc # qπ ip qπ jπ (1) (1) jq jq I (a)Iss (b) − A26 I ip (a)Ics (b) aij a b ss b b cs

(1) iπ +A12

"

pπ ip (1) iπ pπ ip jq jq Isc (a)Isc I (a)Icc (b) + A66 (b) b a a a cc

(1) jπ

+ − A26

30

(1) jπ

+A22

qπ ip (1) iπ qπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − A26 (b) b b a b # " # (1) (2) ip jq ip jq +K1 I˜ss (a)I˜cc (b) b + − K1 I˜ss (a)I˜cc (b) b ij

"

ij

pπ ip (1) jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + B66 (b) a a sc b a cc

(1) iπ

+ − B16

# qπ ip (1) jπ qπ ip (1) jq jq I (a)Iss (b) − B26 I (a)Ics (b) dij a b ss b b cs

(1) iπ +B12

"

pπ ip (1) iπ pπ ip jq jq I (a)Isc I (a)Icc (b) + B66 (b) b a sc a a cc

(1) jπ

+ − B26

(1) jπ

qπ ip (1) iπ qπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − B26 (b) b b a b # " # ) h1 ˜ip ˜jq h2 ˜ip ˜jq (1) (2) −K1 Iss (a)Icc (b) eij + − K1 Iss (a)Icc (b) eij =0 2 2

+B22

(D.2c) (" M,N X i,j

pπ ip (2) jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + A66 (b) a a sc b a cc

(2) iπ

− A16

# qπ ip (2) jπ qπ ip (2) jq jq I (a)Iss (b) − A26 I (a)Ics (b) aij a b ss b b cs " # " (1) ip jq ip jq + − K1 I˜ss (a)I˜cc (b) b + K1 I˜ss (a)I˜cc (b) (2) iπ +A12

ij

pπ ip (2) iπ pπ ip jq jq Isc (a)Isc (b) + A66 I (a)Icc (b) b a a a cc # (2) jπ qπ ip (2) iπ qπ ip (2) jq jq +A22 I (a)Iss (b) − A26 I (a)Ics (b) bij b b ss a b cs (2) jπ

−A26

"

pπ ip (2) jπ pπ ip jq jq I (a)Isc (b) + B66 I (a)Icc (b) a a sc b a cc

(2) iπ

+ − B16

# qπ ip (2) jπ qπ ip (2) jq jq I (a)Iss (b) − B26 I (a)Ics (b) dij a b ss b b cs " # " h1 ˜ip ˜jq h2 ip ˜jq (1) + K1 Iss (a)Icc (b) eij + K1 I˜ss (a)Icc (b) 2 2 (2) iπ +B12

pπ ip (2) iπ pπ ip jq jq Isc (a)Isc I (a)Icc (b) + B66 (b) b a a a cc

(2) jπ

−B26

31

# ) qπ ip (2) iπ qπ ip (2) jq jq =0 I (a)Iss (b) − B26 I (a)Ics (b) eij b b ss a b cs

(2) jπ +B22

(D.2d) M,N,2 X i,j,α

("

pπ ip ts(α) iπ pπ ip jq jq I (a)Ics (b) + K11 I (a)Iss (b) b a sc a a cc

ts(α) jπ

K12

# qπ ip ts(α) iπ qπ ip jq jq I (a)Icc (b) + K12 I (a)Isc (b) cij b b ss a b cs " # ts(α) pπ ip ts(α) qπ ip (α) jq jq + K11 I (a)Iss (b) + K12 I (a)Isc (b) dij a cc b cs ) " # ts(α) pπ ip ts(α) qπ ip (α) jq jq = −Qpq I (a)Ics (b) + K22 I (a)Icc (b) eij + K12 a sc b ss

ts(α) jπ +K22

(D.2e) M,N X i,j

("

pπ ip (1) jπ pπ ip jq jq (b) − B16 (b) I (a)Iss I (a)Ics a a ss b a cs

(1) iπ

B11

qπ ip (1) jπ qπ ip jq jq I (a)Isc (b) + B66 Icc (a)Icc (b) a b sc b b # " # h1 ˜ip ˜jq h1 ˜ip ˜jq (1) (2) −K1 Icc (a)Iss (b) aij + K1 Icc (a)Iss (b) aij 2 2

(1) iπ

−B16

"

(1) jπ

+ B12

pπ ip (1) iπ pπ ip jq jq (b) − B16 (b) Iss (a)Iss I (a)Ics b a a a cs

# qπ ip (1) iπ qπ ip (1) jq jq I (a)Isc (b) + B66 I (a)Icc (b) bij b b sc a b cc " # ts(1) jπ ip ts(1) iπ ip jq jq + K12 I (a)Ics (b) + K11 I (a)Iss (b) cij b sc a cc

(1) jπ −B26

"

pπ ip (1) jπ pπ ip jq jq I (a)Iss I (a)Ics (b) − D16 (b) a a ss b a cs

(1) iπ

+ D11

(1) iπ

qπ ip (1) jπ qπ ip jq jq (b) + D66 (b) Isc (a)Isc I (a)Icc a b b b cc #  2 h1 ts(1) ip (1) jq ip jq ˜ ˜ +K11 Icc (a)Iss (b) + K1 Icc (a)Iss (b) dij 2 " # " h1 h2 ˜ip ˜jq (2) ts(1) ip jq I (a)Iss (b) dij + K12 Isc + K1 (a)Ics (b) 2 2 cc −D16

32

(1) jπ

pπ ip (1) iπ pπ ip jq jq Iss (a)Iss I (a)Ics (b) − D16 (b) b a a a cs ) # (1) jπ qπ ip (1) iπ qπ ip (1) jq jq =0 I (a)Isc (b) + D66 I (a)Icc (b) eij −D26 b b sc a b cc

+D12

(D.2f) (" M,N X i,j

− K1

#

"

h2 ˜ip ˜jq h2 ip ˜jq (1) I (a)Iss (b) aij + K1 I˜cc (a)Iss (b) 2 cc 2

pπ ip (2) jπ pπ ip jq jq I (a)Iss I (a)Ics (b) − B16 (b) a a ss b a cs # (2) iπ qπ ip (2) jπ qπ ip (2) jq jq −B16 I (a)Isc (b) + B66 I (a)Icc (b) aij a b sc b b cc (2) iπ

+B11

"

pπ ip (2) iπ pπ ip jq jq I (a)Iss (b) − B16 I (a)Ics (b) b a ss a a cs

(2) jπ

+ B12

# qπ ip (2) iπ qπ ip (2) jq jq I (a)Isc (b) + B66 I (a)Icc (b) bij b b sc a b cc " # iπ ts(2) jπ ip ts(2) jq jq I (a)Ics I ip (a)Iss + K12 (b) + K11 (b) cij b sc a cc " # "   2 h1 h2 ˜ip ˜jq h2 (1) ip jq + K1 Icc (a)Iss (b) dij + K1 (a)I˜ss (b) I˜cc 2 2 2 (2) jπ −B26

(2) iπ

pπ ip (2) jπ pπ ip jq jq (b) − D16 (b) I (a)Iss I (a)Ics a a ss b a cs (2) iπ qπ ip (2) jπ qπ ip jq −D16 I (a)I jq (b) + D66 I (a)Icc (b) a b sc # sc b b cc "

+D11

ts(2) ip jq Icc (a)Iss (b)

+K11

(2)

ts(2) ip jq Isc (a)Ics (b)

dij + K12

pπ ip (2) iπ pπ ip jq jq Iss (a)Iss I (a)Ics (b) − D16 (b) b a a a cs ) # (2) jπ qπ ip (2) iπ qπ ip (2) jq jq =0 I (a)Isc (b) + D66 I (a)Icc (b) eij −D26 b b sc a b cc (2) jπ

+D12

(D.2g) (" M,N X i,j

(1) iπ

− B16

pπ ip (1) jπ pπ ip jq jq (b) + B66 (b) I (a)Isc I (a)Icc a a sc b a cc # qπ ip (1) jπ qπ ip (1) jq jq I (a)Iss (b) − B26 I (a)Ics (b) aij a b ss b b cs

(1) iπ +B12

33

"

(1) jπ

+ − B26

pπ ip (1) iπ pπ ip jq jq (b) + B66 (b) Isc (a)Isc I (a)Icc b a a a cc

qπ ip (1) iπ qπ ip jq jq I (a)Iss (b) − B26 I (a)Ics (b) b b ss a" b cs # # h1 ˜ip ˜jq h1 ˜ip ˜jq (1) (2) −K1 Iss (a)Icc (b) bij + K1 Iss (a)Icc (b) bij 2 2 " # ts(1) jπ ip ts(1) iπ ip jq jq I (a)Icc (b) + K12 I (a)Isc (b) cij + K22 b ss a cs (1) jπ

+B22

"

pπ ip (1) jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + D66 (b) a a sc b a cc

(1) iπ

+ − D16

(1) iπ

+D12

qπ ip (1) jπ qπ ip jq jq Iss (a)Iss Ics (a)Ics (b) − D26 (b) a b b b # "

ts(1) ip jq Ics (a)Isc (b)

+K12

(1)

ts(1) ip jq Iss (a)Icc (b)

dij + K22

pπ ip (1) iπ pπ ip jq jq I (a)Isc (b) + D66 I (a)Icc (b) b a sc a a cc (1) jπ qπ ip (1) iπ qπ ip jq jq I (a)Iss Ics (a)Ics (b) − D26 (b) +D22 b b ss a b ) # " #  2 h1 h h 1 2 (1) (2) ip jq jq I˜ss +K1 =0 I˜ip (a)I˜cc (a)I˜cc (b) eij + K1 (b) eij 2 2 2 ss (1) jπ

−D26

(D.2h) (" M,N X i,j

pπ ip (2) jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + B66 (b) a a sc b a cc

(2) iπ

− B16

# qπ ip qπ jπ (2) (2) jq jq I (a)Iss (b) − B26 I ip (a)Ics (b) aij a b ss b b cs " # " h2 ˜ip ˜jq h2 ip ˜jq (1) + − K1 Iss (a)Icc (b) bij + K1 I˜ss (a)Icc (b) 2 2 (2) iπ +B12

(2) jπ

pπ ip (2) iπ pπ ip jq jq (b) + B66 (b) I (a)Isc I (a)Icc b a sc a a cc # qπ iπ (2) jπ qπ ip (2) (2) jq jq +B22 I (a)Iss I ip (a)Ics (b) − B26 (b) bij b b ss a b cs " # ts(2) jπ ip ts(2) iπ ip jq jq I (a)Icc (b) + K12 I (a)Isc (b) cij + K22 b ss a cs

−B26

"

pπ ip (2) jπ pπ ip jq jq I (a)Isc I (a)Icc (b) + D66 (b) a a sc b a cc

(2) iπ

+ − D16

34

(2) iπ

qπ ip (2) jπ qπ ip jq jq (b) − D26 (b) Iss (a)Iss Ics (a)Ics a b b b # " # h1 h2 ˜ip ˜jq ts(2) ip (2) (1) jq +K12 Ics (a)Isc (b) dij + K1 I (a)Icc (b) eij 2 2 ss

+D12

"

pπ ip (2) iπ pπ ip jq jq I (a)Isc I (a)Icc (b) + D66 (b) b a sc a a cc

(2) jπ

+ − D26

qπ ip (2) iπ qπ ip jq jq Iss (a)Iss I (a)Ics (b) − D26 (b) b b a b cs # )  2 h2 ts(2) ip (2) jq ip jq ˜ ˜ Iss (a)Icc (b) eij +K22 Iss (a)Icc (b) + K1 =0 2

(2) jπ

+D22

(D.2i) According to the definition of the linear elastic-brittle law, the delaminated region of the interface has zero shear traction and zero secant modulus while the remaining uncracked region still follows the elastic traction-separation relation. Therefore, the integrations related to the cohesive elastic stiffness K1 are only performed over the uncracked region of the interface     Z Z  pπx   qπy  iπx jπy ip jq (a)I˜ss (b) = cos I˜cc cos sin sin dxdy a a b b uncracked     Z Z  pπx   qπy  iπx jπy ip jq (a)Iss (b) − cos cos sin sin dxdy = Icc a a b b delam.     Z Z  pπx   qπy  iπx jπy ip jq (a)I˜cc (b) = sin I˜ss sin cos cos dxdy a a b b uncracked     Z Z  pπx   qπy  iπx jπy ip jq (a)Icc (b) − cos cos sin sin dxdy = Iss a a b b delam. (D.3) For a pristine interface, i.e., a interface without any delamination, ip jq ip jq I˜cc (a)I˜ss (b) = Icc (a)Iss (b),

ip jq ip jq I˜ss (a)I˜cc (b) = Iss (a)Icc (b)

(D.4)

For a given shape of delamination, the integrations, as shown in Eqn. (D.3), can be analytically carried out by approximating the delamination area as a collection of small rectangles, or computed by numerical integration methods. (α)

(α)

(α)

(α)

Above all, the coefficients, aij , bij , cij , dij , eij , with a total number of 9 × M × N , can be solved by 9 × M × N coupled equations. 35

Appendix E. The Rayleigh-Ritz Method with the Classical Lamination Theory (CLT) CLT assumes the transverse normals remain perpendicular to the mid-plane during deformation, namely, φx = −

∂w ∂w , φy = − ∂x ∂y

(E.1)

For elastic bending of an anisotropic plate, as shown in Section 2, the weak form of governing equations are    Z b Z a" ∂ δu ∂ δu ∂ δv ∂ δv + Nxy Nx + Nxy + Ny ∂x ∂y ∂x ∂y 0 0  # ∂ 2 δw ∂ 2 δw ∂ 2 δw + My − Mx + 2Mxy − qδw dxdy = 0 ∂x2 ∂x∂y ∂y 2

(E.2)

Additionally considering a cohesive layer at a potential crack interface of a plate, as described in Section 2, the weak form becomes    Z b Z a ( (1) (1) (2) (2) (1) ∂ δu (1) ∂ δu (1) (2) ∂ δu (2) ∂ δu (2) + Nx Nx + Nxy + τx δu + Nxy − τx δu ∂x ∂y ∂x ∂y 0 0     (1) (1) (2) (2) (1) ∂ δv (1) ∂ δv (1) (2) ∂ δv (2) ∂ δv (2) + Nxy + Ny + τy δv + Ny − τy δv + Nxy ∂x ∂y ∂x ∂y "    ∂ 2 δw  2 (1) (2) ∂ δw + 2 M + M − Mx(1) + Mx(2) xy xy ∂x2 ∂x∂y #)    2  h ∂ δw ∂ δw (1) (2) ∂ δw + τy τx − qδw dxdy = 0 + My + My − ∂y 2 2 ∂x ∂y (E.3) References References [1] S. Abrate, Impact on composite structures, Cambridge university press, 2005. [2] G. A. O. Davies,

P. Robinson,

Predicting failure by debond-

ing/delamination, in: AGARD: 74th Structures & Materials Meeting, Debonding/Delamination of Composites, no. 5, Patras, Greece, 1992. 36

[3] D. J. Elder, R. S. Thomson, M. Q. Nguyen, M. L. Scott, Review of delamination predictive methods for low speed impact of composite laminates, Composite Structures 66 (1) (2004) 677–683. [4] J. Xie, A. M. Waas, Predictions of delamination growth for quasi-static loading of composite laminates, Journal of Applied Mechanics 82 (8) (2015) 081004. [5] R. Olsson, A review of impact experiments at FFA during 1986 to 1998, The Aeronautical Research Institute of Sweden. FFA TN 8. [6] G. A. Schoeppner, S. Abrate, Delamination threshold loads for low velocity impact on composite laminates, Composites Part A: Applied Science and Manufacturing 31 (9) (2000) 903–915. [7] H. Suemasu, O. Majima, Multiple delaminations and their severity in circular axisymmetric plates subjected to transverse loading, Journal of Composite Materials 30 (4) (1996) 441–453. [8] H. Suemasu, O. Majima, Multiple delaminations and their severity in nonlinear circular plates subjected to concentrated loading, Journal of Composite Materials 32 (2) (1998) 123–140. [9] R. Olsson, Analytical prediction of large mass impact damage in composite laminates, Composites Part A: Applied Science and Manufacturing 32 (9) (2001) 1207–1215. [10] R. Olsson, M. V. Donadon, B. G. Falzon, Delamination threshold load for dynamic impact on plates, International Journal of Solids and Structures 43 (10) (2006) 3124–3141. [11] D. Liu, Impact-induced delamination – a view of bending stiffness mismatching, Journal of Composite Materials 22 (7) (1988) 674–692. [12] A. J. Lesser, A. G. Filippov, Mechanisms governing the damage resistance of laminated composites subjected to low-velocity impacts, International Journal of Damage Mechanics 3 (4) (1994) 408–432. 37

[13] G. I. Barenblatt, The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks, Journal of Applied Mathematics and Mechanics 23 (3) (1959) 622–636. [14] G. I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Advances in Applied Mechanics 7 (1962) 55–129. [15] D. S. Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids 8 (2) (1960) 100–104. [16] M. F. Kanninen, An augmented double cantilever beam model for studying crack propagation and arrest, International Journal of Fracture 9 (1) (1973) 83–92. [17] R. Olsson, A simplified improved beam analysis of the DCB specimen, Composites Science and Technology 43 (4) (1992) 329–338. [18] J. G. Williams, H. Hadavinia, Analytical solutions for cohesive zone models, Journal of the Mechanics and Physics of Solids 50 (4) (2002) 809–825. [19] A. B. de Morais, Novel cohesive beam model for the end-notched flexure (ENF) specimen, Engineering Fracture Mechanics 78 (17) (2011) 3017– 3029. [20] S. Bennati, P. Fisicaro, P. S. Valvo, An enhanced beam-theory model of the mixed-mode bending (MMB) test-part I: Literature review and mechanical model, Meccanica 48 (2) (2013) 443–462. [21] S. Bennati, P. Fisicaro, P. S. Valvo, An enhanced beam-theory model of the mixed-mode bending (MMB) test-part II: Applications and results, Meccanica 48 (2) (2013) 465–484. [22] J. Xie, A. M. Waas, M. Rassaian, Closed-form solutions for cohesive zone modeling of delamination toughness tests, International Journal of Solids and Structures 88 (2016) 379–400.

38

[23] J. Xie, A. M. Waas, M. Rassaian, Estimating the process zone length of fracture tests used in characterizing composites, International Journal of Solids and Structures 100 (2016) 111–126. [24] Y. Hirai, H. Hamada, J.-K. Kim, Impact response of woven glass-fabric compositesi.: Effect of fibre surface treatment, Composites Science and Technology 58 (1) (1998) 91–104. [25] L. Ma, D. Liu, Delamination and fiber-bridging damage analysis of angleply laminates subjected to transverse loading, Journal of Composite Materials 50 (22) (2015) 3063–3075. [26] S. I. Thorsson, J. Xie, J. Marek, A. M. Waas, Matrix crack interacting with a delamination in an impacted sandwich composite beam, Engineering Fracture Mechanics 163 (2016) 476–486. [27] H. Chai, C. D. Babcock, W. G. Knauss, One dimensional modelling of failure in laminated plates by delamination buckling, International Journal of Solids and Structures 17 (11) (1981) 1069–1083. [28] D. Bruno, F. Greco, An asymptotic analysis of delamination buckling and growth in layered plates, International journal of solids and structures 37 (43) (2000) 6239–6276. [29] G. A. O. Davies, R. Olsson, Impact on composite structures, Aeronautical Journal 108 (1089) (2004) 541–563. [30] R. Olsson, Analytical prediction of damage due to large mass impact on thin ply composites, Composites Part A: Applied Science and Manufacturing 72 (2015) 184–191. [31] J. N. Reddy, Mechanics of laminated composite plates and shells: theory and analysis, CRC Press, Boca Raton, FL, 2004. [32] Dassault Syst`emes Simulia Corp., Johnston, RI, Abaqus 6.14 documentation: Abaqus theory guide (2016). 39

[33] R. Krueger, Development and application of benchmark examples for mixed-mode I/II quasi-static delamination propagation predictions, Tech. Rep. NASA/CR-2012-217562, NIA Report No. 2012-01, NF1676L-14133, NASA Langley Research Center, Hampton, Virginia, USA (2012). [34] Dassault Syst`emes Simulia Corp., Johnston, RI, Abaqus 6.14 documentation (2016). [35] X. Zhang, F. Bianchi, H. Liu, Predicting low-velocity impact damage in composites by a quasi-static load model with cohesive interface elements, The Aeronautical Journal 116 (1186) (2012) 1367–1381.

40

Figure 10: Normalized interfacial shear traction distribution over lower four interfaces of laminate with stacking sequence (0/90/+45/-45/0/-45/+45/90/0) by the Rayleigh-Ritz method (left) and FE simulations (right).

Figure 11: Normalized interfacial shear traction distribution over upper four interfaces of laminate with stacking sequence (0/90/+45/-45/0/-45/+45/90/0) by the Rayleigh-Ritz method (left) and FE simulations (right).

41