Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method

Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method

Journal Pre-proof Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method Nikhil Garg Conce...

1MB Sizes 0 Downloads 68 Views

Journal Pre-proof

Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method Nikhil Garg ConceptualizationMethodologySoftwareWriting - Original DraftVisualization , Nilanjan Das Chakladar MethodologySupervisionWriting - Review & Editing , B. Gangadhara Prusty SupervisionWriting - Review & EditingProject administratorFunding acquisition , Chongmin Song ResourcesSupervisionWriting - Review & Editing , Andrew W. Phillips SupervisionWriting - Review & EditingFunding acquisition PII: DOI: Reference:

S0020-7403(19)33614-8 https://doi.org/10.1016/j.ijmecsci.2019.105349 MS 105349

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

24 September 2019 26 November 2019 27 November 2019

Please cite this article as: Nikhil Garg ConceptualizationMethodologySoftwareWriting - Original DraftVisualization , Nilanjan Das Chakladar MethodologySupervisionWriting - Review & Editing , B. Gangadhara Prusty SupervisionW Chongmin Song ResourcesSupervisionWriting - Review & Editing , Andrew W. Phillips SupervisionWriting - Review Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method, International Journal of Mechanical Sciences (2019), doi: https://doi.org/10.1016/j.ijmecsci.2019.105349

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier Ltd.

Highlights    

SBFEM is explored for laminated composite plates with weakly bonded interfaces. Proposed mesh reduces the discretisation difficulties and computational cost. It can handle smaller interface elements along with large adjoining bulk element. Proposed formulation can predict more accurate results than other available models.

1

Modelling of laminated composite plates with weakly bonded interfaces using scaled boundary finite element method Nikhil Garga,*, Nilanjan Das Chakladar b, B. Gangadhara Prustyc, Chongmin Songd, Andrew W. Phillipse a,c

ARC Training Centre for Automated Manufacture of Advanced Composites, School of Mechanical and Manufacturing Engineering, UNSW Sydney, Australia b d e *

Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur 721302, India. School of Civil and Environmental Engineering, UNSW Sydney, Australia Maritime Division, Defence Science and Technology Group (DST Group), Melbourne, Australia Corresponding Author, Ph No. +61-410218735; E-mail address: [email protected]; Postal

Address: Desk No. 12, Room 208, J-17, Ainsworth Building, UNSW, Sydney, Australia, 2033

Abstract Scaled Boundary Finite Element Method (SBFEM), a relatively new approach, is investigated for the analysis of the weakly bonded laminated composite plates, under cylindrical bending. A plane strain formulation, which does not require any assumptions on the displacement or stress field, is developed. The interface between two plies is simulated using a linear spring layer model and incorporated in the numerical modelling with the help of zero thickness interface elements. The model is coded and postprocessed in MATLAB to capture general delamination, i.e. separation of the plies at interfaces in both tangential and normal direction. Due to the semi-analytical solution procedure of SBFEM, the model provides highly accurate results. Use of SBFEM reduces the meshing burden and complexities in the discretisation. Moreover, the method can handle hanging nodes, so one can use small interface elements without decreasing the size of adjoining domains. These qualities make the proposed numerical model computationally cheaper as well as very accurate. The presented formulation is demonstrated for a variety of examples in laminated composites including unidirectional laminates, symmetric cross ply laminates and antisymmetric cross ply laminates. The proposed model’s close agreement as well as its superiority over the other methods is established for future reference. Keywords: Laminated composites, Plates, Cylindrical bending, Weak interface, Scaled boundary finite element method, SBFEM.

2

Nomenclature - Transformed rigidity co-efficient

[

,

- Coefficients describing the extent of the delamination (or imperfection)

,

- Shape functions of a two noded line element in FEM ] - Matrix of the shape functions

- Shear slip coefficient - Normal separation coefficient - Traction in normal and tangential direction l - Length of the plate h - Height of the plate { } - Displacement vector - Displacement jump in normal and tangential direction - Displacement jump in x and z direction - Jump between node i and node j in x and z direction i.e.

,

{ }, { } - Vectors of the nodal co-ordinates - Stress - Strain - Radial co-ordinate in SBFEM - Circumferential co-ordinate in SBFEM - Eigen values { } - Corresponding eigen vector - Derivative with respect to

1.

Introduction

Composite materials are being used in ever increasing array of engineering components due to their many advantageous properties like high stiffness and strength to weight ratios, high fatigue performance and resistance to corrosion. Moreover, exploiting the directional properties, a laminated composite structure can be designed as per the requirement by tailoring the laminate properties in different directions. The high use of composites has also increased the demand for their accurate modelling under stringent conditions. Often exact solutions are not available due to the geometric and material complexities of the composite 3

structures. Approximate methods were proposed to model the mechanical behaviour of laminated composites which includes shear deformation theories [1–4], layer-wise theories [5,6] and zigzag models [7–10]. Readers may refer to detailed reviews on these plate theories [11–14]. Copious literature is available on assumption of a perfectly bonded interfaces between the layers, i.e. interlaminar displacement continuity is retained. However, perfect interfaces do not exist in reality and are only an idealisation of a complex situation. Various defects such as microcracks, inhomogeneities, and voids may get introduced at the interfaces during the manufacturing process of the laminates [15]. Also, in the service life, damage is initiated at interfaces which can escalate the local failure to the loss of structural integrity of the component itself. Accordingly, there is a requirement to study the effect of interfacial damage or delamination on the global response of the laminate, so that the model can more realistically represent the interface behaviour. Delamination can be categorised into two types, local delamination and weakly bonding of layers [16]. The local delamination characterizes a complete separation of the two layers, observed on a localised area. In this case, the interfacial stresses in that area becomes zero and there is no interaction between the delaminated plies. Whereas weakly bonded layers refer to a weak interface, in which the whole interface is only partial de-bonded. In this type of de-bonding, the interfacial stresses are non-zero and generally assumed to be a function of the interfacial displacement jump. One of the ways to model the local delamination is to use the sub-laminate approach [17–23]. However, this paper is limited to the second type of delamination, i.e. imperfectly bonded interfaces, hence the literature covering the local delamination phenomenon is out of scope of this research. A detailed literature review dealing with the modelling of imperfect interfaces is presented in the upcoming segments of this section. Initial theoretical studies investigating the imperfect interfaces started by using shear slip theory [24–26], where a tangential discontinuity was allowed at the interface and separation was related to the interfacial shear stress. In shear slip theory, the opening at the interface was not taken into consideration. Studies were then extended to model general delamination considering both normal as well as shear imperfections [27–29]. However, some authors [30,31] prefer not to use normal separations in order to avoid layer interpenetration. Due to its simplicity, a spring layer model was adopted in most of these studies which assumes that

4

interface behaves in a linearly elastic manner. In these spring layer models, the traction was assumed to be directly proportional to the displacement jumps at the interfaces. It can be understood as a specific case of cohesive zone modelling (CZM), which does not consider the softening of the interface. Attempts were made to extend plate theories to incorporate weak interfaces. Shear deformation theories predict the global response with sufficient accuracy however they are not able to provide good insights into the state of the interfacial stresses. Layer-wise theories can overcome this drawback and hence are exploited to study the weak interfaces. Lu and Liu [26] extended the layer-wise theory so that it can capture the interfacial jumps in the tangential directions. Later, the theory was further developed to incorporate the normal separation [29] along with the shear slip. A generalized model for delamination, which includes the shear as well as normal separation, is also presented by Williams and Addessio [27]. Drawback of layer-wise models is that they are computationally very heavy. Zigzag models reduce the computational cost without compromising the accuracy much and hence, are extensively used to incorporate displacement jumps at interfaces. Cheng and colleagues [30,31] presented the results for cylindrical [31] and rectangular plate bending [30] using the third order zigzag theory. A mixed formulation (stress and displacement) was also provided along with the results from a displacement-based model. Massabo and Campi improved upon the available zigzag theories [32] and applied the proposed efficient theory for imperfect interfaces [33]. The theory was later used to develop the results for the complex plate problems other than simply supported boundary conditions [34]. Kim et al. [35] attempted to incorporate the effect of higher order zigzag model with an efficient first order shear deformation theory. They calculated the effective shear modulus of the first order theory using the strain energy transformation of the zigzag theory. Zigzag theories have also been used to study the imperfect interface in shells [16,36] and geometrically nonlinear plates [30,37–40]. Chen et al. [41] presented the exact solutions for the rectangular cross ply laminates with imperfect interfaces using the state space approach. The solution for cylindrical bending were presented as a special case, where the width of the plate is assumed to be very large. Results of the third order zigzag model [30] were compared with the exact results in the study. It was found that as imperfection increases, the zigzag theory, which is one of the best models to

5

describe laminate behaviour in perfectly bonded conditions, is not able to predict the static response of the plate correctly. Using the similar state space approach, exact solutions were developed for angle ply laminates [42] and laminated composite beams [43]. For the composite beams, plane stress conditions were used to develop the results, i.e. width of the beam is taken to be very small as compared to the thickness. This assumption is not very correct, and plane strain conditions would have represented the model more accurately. Various other studies dealing with imperfect interfaces with the use of finite element method (FEM) are also available in the literature. Bui et al. [44,45] provided FEM solutions for the problem and studied the effect of imperfect interfaces on bending, buckling [44], interlaminar stresses and strain energy release rate [45]. Chakraborty and Sheikh [46] attempted to formulate a new element type for plate analysis which can incorporate displacement jumps at the interfaces. Nevertheless, use of FEM could be computationally intensive for real structures. The present study attempts to model the effect of imperfect interface on wide laminated plates using a relatively new technique, Scaled Boundary Finite Element Method (SBFEM). The method was first developed by Song and Wolf [47] and has various advantages over traditional FEM. SBFEM is a semi-analytical technique which solves the boundary of the domain (circumferential direction) numerically whereas an analytical solution is obtained in the radial directions. It provides an upper edge over the accuracy of the results with fewer nodes. Moreover, the method only requires the discretisation of the boundary, hence the discretisation process is reduced by one spatial dimension. It greatly saves the required computational effort as well as eases meshing of a complex domain. SBFEM has proven to be very advantageous in fracture mechanics studies [48–51]. However, its utility in delamination studies of composite materials has not been explored yet. As per the author’s knowledge, there are no published results available which deals with the implementation of SBFEM for the analysis of laminated composite plates with imperfect interfaces. In this paper, SBFEM is used to study cylindrical bending of laminated composite plates with interlaminar bonding imperfections. A plane strain formulation is used which does not need any plate theory based approximations on the displacement field. The results obtained are compared against best practice layer wise, zigzag and FEM model approaches. The proposed model is able to capture the behaviour of general delamination, i.e. the effect of normal separation can be incorporated along with shear slip. Furthermore, presented SBFEM

6

formulation showed very accurate results even for the thick laminates where shear deformation theories generally fail. This high accuracy does not come at an expense of the high computational cost as meshing requirements are minimal due to the semi-analytical nature of the method. The paper is set out as follows. First the mathematical formulation is introduced including the meshing strategy, constitutive equations, SBFEM and interface element formulations. In the following section the material properties and boundary conditions of the cylindrical bending problem is described. Finally, the developed SBFEM approach is compared against best practice modelling approaches for a range of different laminate sequences. 2.

Mathematical formulations

The following section presents the mathematical formulation developed for studying the effect of weak interfaces in composite laminates subjected to cylindrical bending. The coordinate system adopted for modelling the laminated plates is shown in Figure 1. The plane strain formulation is implemented in the xz plane. The origin is chosen to be at the left edge of the plate at mid height. The length of the plate is denoted with l and total thickness with h. z

h/2

x

h/2 l

Figure 1: Schematic showing the cylindrical bending problem and adopted coordinate system (The dark shaded regions are individual lamina while the light shaded region represents the zero-thickness interface) 2.1

Meshing

In traditional methods like FEM, the most common elements are triangular or quadrilateral, whereas SBFEM can have polygonal elements with any number of sides. Furthermore, each side of the polygon can be discretised into any number of 1D line elements. Hence, even if small interface elements are used, structures can have large elements adjoining them (Figure 2). Traditional FEM cannot solve the problem with this type of meshing due to the presence of hanging nodes (vertex of one element lies at the interior of the edge of another element). Hanging nodes usually occur in FEM when adaptive meshing is used. Due to this 7

requirement of mesh conformity in FEM, the size of interface element will determine the aspect ratio of the other regular elements too. The meshing strategy used in the present study is represented in Figure 2. Meshing is shown for a two-layered laminate with l/h = 2. Each lamina is divided into rectangular domains (aspect ratio is kept close to one). The boundary of these domains is discretised into line elements. In the figure, each side of the rectangle is shown to have four line-elements. This type of mesh is represented as ‘4×4’ mesh in further sections of the document. Interface elements (shown in grey) connects the individual layers. A four noded zero thickness interface element is adopted in the modelling. Due to zero thickness, the nodes at the top edge of the bottom layer and the bottom edge of the top layer overlap each other. For the sake of clarity, interface elements are shown with some thickness in the figure. Lamina-1

Interface

4

3

1

2

Interface Element Lamina-2

Figure 2: Adopted meshing strategy under framework of SBFEM 2.2

Constitutive equation

A generalised Hooke’s law for any transversely isotropic material can be expressed as follows

⌈ ⌈ ⌈ ⌈ ⌈

⌉ ⌉ ⌉ ⌉ ⌉

⌈ ⌈ ⌈ ⌈ ⌈ ⌈ [

⌉ ⌉⌈ ⌉⌈ ⌉⌈ ⌉⌈ ⌉⌈ ]

⌉ ⌉ ⌉ ⌉ ⌉

Plane strain formulation is adopted in the xz plane, hence, then reduces to 8

(1)

Equation (1)

[

]

[

2.3

][

]

[ ]

[ ][ ]

(2)

SBFEM formulations

In the present study, two-noded line elements are used to describe the boundary mesh. Only the key equations necessary for the implementation are summarised below. Details of the derivation of the SBFEM formulations can be found in [52]. As per the SBFEM formulations, the co-ordinates can be expressed as (Figure 3) [

]{ }

(3)

[

]{ }

(4)

𝜉 Line Element

𝜂 Scaling centre

Figure 3: A typical polygon SBFEM element with one line element per side The degrees of freedom can be expressed in terms of shape functions as follows {

}

[

]{

}

(5)

Based on linear strain-displacement relations, an operator [L] can be defined as

[ ]

⌈ ⌈ ⌈ ⌈ [

⌉ ⌉ ⌉ ⌉ ]

(6)

Using the Jacobian relationships,

9

[ ]

[

]

[

]

(7)

where the matrices

[

]

[

]

Here, | |

|

|

|

|

[

]

[

(8)

]

(9)

| represents the modulus of the jacobian matrix, and can be expressed as

|

(10)

Strains can be represented as {

}

[ ]{

{

}

[

} ]{

(11) }

[

]{

}

(12)

where, [

]

[

][

]

(13)

[

]

[

][

]

(14)

In a similar manner, stresses are expressed as {

}

[ ]{

}

[ ][

]{

}

[ ][

]{

}

(15)

Using the principle of virtual work, the equation of equilibrium for SBFEM discretisation can be written as [

] {

}

[

]

[ ]

[

]

{

}

[

]{

}

(16)

where, [

]

∫ [

] [ ][

]|

|

(17)

[ ]

∫ [

] [ ][

]|

|

(18)

10

[

]

∫ [

] [ ][

]|

|

(19)

The equation obtained by virtual work principal can be transformed into the following differential equation {

{ {

} } }

{ { {

} } }

(20)

where, [ ]

[

[

[

and { {

]

] [ [ ][

] ] [

[

] ] [ ][ ]

]

(21)

} represents the nodal force function, which is expressed as }

[

] {

}

[ ] {

}

(22)

The differential equation can be solved by using the Eigen value solution by solving the equation [ ]

[ ] { }

(23)

defining, { }

{

}

(24)

and selecting the eigen vectors corresponding to the eigenvalues with positive real part

[ ]

, [ ,

- ,

-

-

-

,

]

* [ *

+ ] +

(25)

The stiffness matrix for each S-element is calculated as [ ]

*

+*

+

(26)

Assembly of the stiffness matrices for all the S-elements and interface elements (described in section 2.4) is performed and nodal displacements are obtained using a similar process as FEM. Displacements inside an S-elements are represented by

11

{

}

*

{ }

+

(27)

represents the matrix of the eigen values with positive real part and { }

where,

represents vector of integration constants. The value of

and { } are found out by the

eq (28) and eq (29) respectively. (28) { }

+ {

*

} when

(boundary nodes)

(29)

Stresses are obtained by { } ,

,

∑ [ ]

-

2.4

[

]

(30)

[

] ,

-

(31)

Interface elements formulation

Interface elements are formulated using the conventional FEM approach [53,54]. A linear uncoupled model is adopted for the interface behaviour. It is assumed that the traction developed is directly proportional to the displacement jumps, hence, [ ]

[

][

]

(32)

It can be observed that when

, the traction becomes zero and represents the

completely de-bonded layers. Similarly,

represents perfectly bonded laminate. The

above equation relates the traction with the displacement jumps and sometimes referred as constitutive equation of the interface elements. To represent the displacement discontinuity in global axes (x,z), rotation of the local axes (normal (n) and tangential (t) directions) is required to be performed. [

]

*

+[

]

[ ][

]

(33)

Displacement jumps can be expressed in terms of the nodal values with the help of shape functions as follows

12

[

]

⌈ ]⌈ ⌈ [

[

⌉ ⌉ ⌉ ]

⌈ [ ]⌈ ⌈ [

⌉ ⌉ ⌉ ]

(34)

Here, superscript (1,4) and (2,3) are local node numbers in an interface elements (Figure 2). Now, by defining an operator matrix, these nodal displacement jumps can be related to the displacement vector.

⌈ ⌈ ⌈ [

⌉ ⌉ ⌉ ]

⌈ ⌈ ⌈ [

⌉ ⌉ ⌉ ]

][ ]

[

[ ][ ]

(35)

So, the displacement jumps can be expressed in terms of nodal displacements using the following equation [

]

[ ][ ][ ][ ]

[ ][ ]

(36)

Hence, using constitutive equation, tractions can be expressed as [

[ ][

]

]

[ ][ ][ ]

The potential energy

(37)

for the interface element can be expressed in terms of the path

integral, ∫[

] [

]

(38)

Here, represents the width of the plate. Using eq (36), (37) and (38) [ ] (∫[ ] [ ][ ]

)[ ]

(39)

Hence, the stiffness matrix for the interface elements (considering unit width) can be expressed in similar terms to the typical finite element stiffness ∫[ ] [ ][ ]

(40)

13

3.

Modelling details

The following section describes the material properties, boundary and loading conditions and normalising procedures used in the numerical examples. Material Properties: In the examples described below, all the layers of the laminates are assumed to be made up of same orthotropic material. Elastic properties of the material [46] are described as follows which are suitably normalised with respect to

. (41)

Loading and boundary conditions: Sinusoidal loading conditions are applied at the top surface of the plate as described in the following equation. (

)

(42)

Simply supported boundary conditions are considered. The plate is restrained to move along direction at the edges, i.e.

and

required to be restrained in the asserted that the point ( ⁄

. For the numerical considerations, the model is

direction for at least one point. Due to symmetry, it can be ⁄ ) will not undergo a displacement in

hence, the same point is selected for restraining the displacement in the

direction (Figure 4), direction.

z

x

z

x l/2

l/2

Figure 4: Boundary conditions for cylindrical bending used in numerical examples 14

Non-dimensional parameters: The results presented in further sections are shown in the non-dimensional form. The parameters used to get the normalised form are shown as follows

̅

(

)

(

̅

(

̅

⁄ As the coefficient

)

) ̅

̅ (43)

̿

̿

( )



; and

are inversely proportional to

the perfect interface whereas

4.

,

signifies

represents the fully de-bonded layers. In the

numerical computations, when results are calculated for assigned. Similarly, for

and

or

or

, a value of 10-12 is

, a value of 1012 is used.

Numerical results and discussion

The SBFEM methodology is used to study and characterise five laminates with different stacking sequences. The results are compared with the current best practice modelling approaches. First, results for a (0/90/0) plate are presented along with the convergence study. Following that, unidirectional plate with one interface (0/0) is considered. Next, unidirectional plate with two interfaces (0/0/0) is considered. Results for equal as well as unequal individual ply thickness are included. After that, response of a seven layered ((0/90)3/0) symmetric cross ply laminates is studied. Finally, a two layered antisymmetric laminate (0/90) is analysed.

15

Table 1: Non-dimensional deflection for various extents of shear imperfections with the convergence study for (0/90/0) laminate Non-dimensional deflection

along

Stacking sequence

l/h

Reference

0/90/0

4

SBFEM (2×2) SBFEM (4×4) SBFEM (8×8) SBFEM (16×16) SBFEM (32×32) Chen et al. [41] Pagano [55]

2.8465 2.8703 2.8786 2.8867 2.8886 2.8872 2.8872

3.1493 3.1793 3.1839 3.1915 3.1944 3.1956 -

3.4209 3.4537 3.4621 3.4704 3.4698 3.4727 -

3.6696 3.6987 3.7130 3.7195 3.7197 3.7231 -

0.31 0.08 0.06 -

SBFEM (8×8) SBFEM (16×16) Chen et al. [41] Pagano [55]

0.9351 0.9392 0.9316 0.9316

0.9906 0.9887 1.0022 -

1.0697 1.0672 1.0712 -

1.1478 1.1342 1.1385 -

0.62 0.73 -

10

4.1

Avg. % Error

Three layered symmetric cross ply laminate (0/90/0)

A three layered symmetric cross ply laminate (0/90/0) is considered for the initial analysis. All the layers in the laminate are considered to possess equal thickness. At first, a convergence study on the presented formulation is performed. Results for non-dimensional deflection ( ̅ are presented in Table 1, for two different values of span to thickness ratio. For this study, the interface is taken to be perfect in the normal direction i.e. . Non-dimensional deflection is calculated for various mesh sizes varying from 2×2 to 32×32. It can be observed that converged results are obtained for a 16×16 mesh size, hence, in all the other numerical examples solved in this paper, a mesh size of 16×16 is adopted throughout. The results obtained are compared with the exact solutions provided by the Chen et al. [41]. It is observed that SBFEM formulations are very accurate and able to predict the results within 1% error. Also, the results of the perfectly bonded laminate, i.e.

are

very close to the exact solutions provided by the Pagano[55] for the case of cylindrical bending.

16

4.0

͞w

5.0

͞w

SBFEM

SBFEM

4.0

Cheng et al.

3.0

Bui et al.

Cheng et al. Bui et al.

3.0

2.0 2.0 1.0

Rn = 0.0 Rt = 0.2

Rn = 0.0 Rt = 0.6

1.0

0.0

0.0 0

4

8

12

16

20

24

0

Span to thickness ratio (l/h)

(a)

4 8 12 16 20 24 Span to thickness ratio (l/h)

(b)

Figure 5: Variation of non-dimensional deflection of (0/90/0) laminate with span to thickness ratio for (a) small imperfection (b) large imperfection in comparison to results of Cheng et al. [31] and Bui et al. [44] Non-dimensional deflection ̅ for the laminate at

for various span to thickness

ratios is recorded and plotted in Figure 5. Results are obtained for both small large

and

imperfections and are shown in Figure 5 (a) and Figure 5 (b) respectively.

Results are obtained by considering shear slip only and the interface is assumed to be perfect in the normal direction. Results of third order zigzag theory provided by the Cheng et al. [31] and FEM solutions for plane strain formulations provided by Bui et al. [44] are also included to present a comparison. It can be observed that the result of the present SBFEM formulations coincide with the results of Bui et al. [44] as similar plane strain formulation is used in both works. However, observations of Cheng et al. [31] deviate from these results slightly. The deviation is more pronounced in the case of larger imperfections as zigzag theory starts to behave less accurately in those conditions.

17

z/h

z/h

0.5

0.5

0.3

0.3

0.1 -0.1 0.0

Rn = 0 Rt = 4

͞τxz

2.0

3.0

1.0

Rn = 0 Rt = 4

0.1 -60

-30

-0.1 0

-0.3

-0.3

-0.5

-0.5

(a)

30

͞σxx 60

(b) 0.5

Rn = 0 Rt = 4

z/h

SBFEM Massabo and Campi Exact Early Models

0.3 0.1

͞u -3.0

-1.5

-0.1 0.0

1.5

3.0

-0.3 -0.5

(c) Figure 6: Variation of (a) ̅ (b) ̅ and (c) ̅ along the thickness for (0/90/0) laminate with l/h=4; Rn =0, Rt=4 in comparison to results of Massabo and Campi [33], Exact solutions [33], and Early Models [33] Moreover, through the thickness variation of the non-dimensional shear stress longitudinal stress ̅

̅

,

and in-plane deflection ̅ , for the (0/90/0) laminate with span to

thickness ratio of 4 is shown in Figure 6 (a), (b), and (c) respectively. The interfaces are assumed to be perfect in the normal direction and imperfections are considered in shear with . SBFEM simulation results are provided along with the results of the generalised zigzag model of Massabo and Campi [33], the exact 2D solutions [33], and results of early models [33]. Massabo and Campi [33] obtained exact 2D solutions by using the formulation of Williams and Addessio [27], and results of early models by using the first and third order zigzag theories. The results of earlier models is seen to deviate from the exact 2D solutions 18

significantly for all three quantities. For longitudinal stress ̅

and in-plane displacement

̅ , the results of both present SBFEM formulation and Massabo and Campi are very close to the exact 2D results. The largest difference is observed in case of the shear stresses. This is not very surprising as transverse shear stresses are very critical in laminated plates and most of the theories are developed and improved upon so that the effect of transverse shear stresses can be more accurately captured. From Figure 6 (a), it can be observed that only SBFEM formulation are able to predict the results as accurately as the exact 2D solution. z/h 0.50

(a) 0.25 Rt = 5x10-5

̅σxx 0.00 -30

-15

0

15

30

45

SBFEM (Rn = 0) SBFEM (Rn = 10) SBFEM (Rn = 5x10⁶) Liu et al. (Rn = 0) Liu et al. (Rn = 10) Liu et al. (Rn = 5x10⁶)

-0.25

-0.50

z/h 0.50

(b) 0.25 Rt = 5x104

̅σxx 0.00 -60

-40

-20

0

20

-0.25

-0.50

40

60

SBFEM (Rn = 0) SBFEM (Rn = 10) SBFEM (Rn = 5x10⁶) Liu et al. (Rn = 0) Liu et al. (Rn = 10) Liu et al. (Rn = 5x10⁶)

Figure 7: Variation of ̅ along the non-dimensional thickness for different degrees of imperfections in normal direction in (0/0) laminate with l/h=4 (a) (b) in comparison to results of Liu et al. [29] 4.2

Unidirectional laminate (0/0) with one interface

In this example, the behaviour of a unidirectional laminated plate is studied which consist of two layers with equal thickness adjoined by one interface. Fibre orientation in both lamina is 19

0°, i.e. along x direction. A thick plate is considered with span to thickness ratio of 4. The analysis is performed for general delamination considering displacement jumps both in normal and tangential directions. Variation of non-dimensional longitudinal stress ̅

along the thickness (at

depicted in Figure 7. Curves are plotted for three different values of . As mentioned earlier, results for

) is

, namely

, i.e. perfect bond in normal

direction, are calculated by numerically assigning the value of 10 -12. Results are obtained for small

as well as large

shear imperfections and presented in

Figure 7(a) and Figure 7(b), respectively. Results of the present formulation are compared with the results of Liu et al. [29], which are obtained by using the layer-wise laminate theory. It can be observed that results of SBFEM are consistent with those of Liu et al. [29] across the entire thickness for both small and large shear imperfections. 0.5 0.3

z/h

0.5

Rt = 5x10-5

0.3

0.1 -0.3 -0.1 0.0

0.6

0.9

1.2

Rt = 5x104

0.1

̅σzz 0.3

z/h

̅σzz

-0.3 -0.1 0.0

-0.3

-0.3

-0.5

-0.5

0.3

0.6

0.9

1.2

SBFEM (Rn = 0) SBFEM (Rn = 10) SBFEM (Rn = 5x10⁶) Liu et al. (Rn = 0) Liu et al. (Rn = 10) Liu et al. (Rn = 5x10⁶)

(b)

(a)

Figure 8: Variation of ̅ along the non-dimensional thickness for different degrees of imperfections in normal direction in (0/0) laminate with l/h=4 (a) (b) in comparison to results of Liu et al. [29] In a similar manner, the variation of ̅

and

̅

along the non-dimensional thickness is

presented in Figure 8 and Figure 9 respectively, and a comparison is made with the results of layer-wise formulation obtained by Liu et al. [29]. One can easily note that the SBFEM results show excellent agreement with the referenced results. It should be noted that in case of very high imperfections both in normal and shear directions

,

the lower lamina does not experience any type of stress. This is due to the high amount

20

imperfections, which cause both laminae to behave almost independently. As only the upper lamina is loaded externally, the stresses are only apparent in that lamina. z/h

z/h 0.5

0.5

Rt = 5x10-5

0.3

0.3

̅τxz

0.1 -1.0 -0.1 0.0

Rt = 5x104

1.0

2.0

3.0

4.0

0.1

̅τxz

-1.0 -0.1 0.0

-0.3

1.0

2.0

3.0

4.0

SBFEM (Rn = 0) SBFEM (Rn = 10)

-0.3

SBFEM (Rn = 5x10⁶)

-0.5

Liu et al. (Rn = 0)

-0.5

(b)

(a)

Figure 9: Variation of ̅ along the non-dimensional thickness for different degrees of imperfections in normal direction in (0/0) laminate with l/h=4 (a) (b) in comparison to results of Liu et al. [29] 4.3

Unidirectional laminate (0/0/0) with two interfaces

Next, a three-layered unidirectional laminate is chosen to study the effect of weak interfaces. All three layers have equal thickness and fibre orientation of 0°. A thick laminate with simply supported edges is considered for the analysis. Results for through the thickness profile of the non-dimensional shear stress

̿

at the edge (x = 0) of the above described

laminate is presented in Figure 10. Interfaces are assumed to be perfect in the normal direction, i.e.

value is taken to be

for all cases. Curves are plotted for two

extreme conditions i.e. perfectly bonded

and fully de-bonded

In addition, one intermediate value of

interfaces.

is also taken into account. To compare the

SBFEM outcomes, results obtained by Massabo and Campi [33] for the same problem are also plotted. Massabo and Campi [33] used the generalised zigzag theory for their study and compared the results with the exact 2D solutions. Their results for this example show a very high similarity to the exact ones. To avoid any confusion by numerous overlaps, exact solutions are omitted from the plots. It can be observed that the results of the present formulation matches with those of Massabo and Campi [33] (and hence, to the exact solutions). The study suggests that the present simple SBFEM formulation is able to predict the results as accurately as those predicted by very refined zigzag models. 21

z/h

0.3

0.3

̿τxz 0.1

0.2

0.3

0.4

SBFEM Massabo and Campi

0.5

Rn = 0 Rt = 0

0.1 -0.1 0.0

z/h

SBFEM Massabo and Campi

0.5

Rn = 0 Rt = 4

0.1

0.5

̿τxz

-0.1 0.0

-0.3

-0.3

-0.5

-0.5

0.1

0.2

(a)

0.3

0.4

0.5

(b)

SBFEM Massabo and Campi

z/h 0.5 0.3

Rn = 0 Rt = ∞

0.1 -0.1 0.0

̿τxz 0.1

0.2

0.3

0.4

0.5

-0.3 -0.5

(c) Figure 10: Variation of ̿ along the thickness at the edge (x=0) of (0/0/0) laminate for various degrees of tangential imperfection (a) perfectly bonded interfaces (b) imperfect interfaces (c) fully de-bonded interfaces for l/h=5 in comparison to results of Massabo and Campi [33] The effect of interfacial imperfections on the deflection is studied in Figure 11. Deflection values are calculated at

and normalised with the maximum value. The maximum

value of deflection is said to be achieved when full de-bond (in the tangential direction) is present. Hence, wmax is calculated with

. Results of the present SBFEM

formulation are compared with the results of the generalised zigzag theory of Massabo and Campi [33]. In addition, Massabo and Campi generated the exact 2D elasticity solutions by extending the Pagano’s model [55] to incorporate the imperfections, using the strategy described by Williams and Addessio [27]. These exact results are also included in the figure. 22

It can be seen that the results of the zigzag formulation of Massabo and Campi [33] have significant deviation from the exact solutions. The deviation is greater when the imperfections are high. This supports the observation made by Chen et al. [41], that the zigzag model loses accuracy when large imperfections are present at the interfaces. Whereas, results obtained by the present SBFEM formulations are almost identical to the exact results. This demonstrate the benefits of the present SBFEM formulation in terms of accuracy in thick plates with large disbonds. w0 / wmax 1.2 Rn = 0 0.9

0.6 SBFEM Massabo and Campi Exact

0.3

Rt 0.0 0

10

20

30

40

50

Figure 11: Variation of normalised central displacement (w0 / wmax) with shear slip coefficient for (0/0/0) laminate with l/h=5 in comparison to results of Massabo and Campi [33] and Exact 2D solutions [33] The same three-layered unidirectional laminate is also analysed considering different thicknesses for individual lamina. The thickness of the top and bottom plies is ⁄ whereas the middle ply is

⁄ . Non-dimensional shear stress

are recorded for thick laminate

̿

and in-plane displacement ̿

and presented in Figure 12 (a) and (b),

respectively. Curves are obtained for a perfectly bonded interfaces in the normal direction and an imperfect interface in shear

. Results are compared with the exact

solutions and the solutions of zigzag theory, both given by Massabo and Campi [33]. It can be observed that the results of Massabo and Campi deviates from the exact solutions significantly, whereas the results of SBFEM formulation are almost as accurate as the exact ones.

23

z/h

-0.1

0.1

̿τxz 0.0

0.2

0.4

0.6

Rn = 0 Rt = 4

0.3

0.3 0.1

z/h

0.5

SBFEM Massabo and Campi Exact

0.5

-0.4

0.8

-0.2

-0.1

-0.3

-0.3

-0.5

-0.5

(a)

̿u 0.0

0.2

0.4

0.6

(b)

Figure 12: Along the thickness variation of (a) non-dimensional shear stresses ̿ and (b) non-dimensional in-plane displacements ̿ at the edge (x=0) of (0/0/0) laminate with unequal ply thickness (h/6, 2h/3, h/6) for l/h=5; in comparison to results of Massabo and Campi [33] and Exact 2D solutions [33] 4.4

Seven layered symmetric cross ply laminate [(0/90) 3/0]

The problem of a seven layered laminate is considered having stacking of 0º and 90º plies alternate to each other. Both of the outer lamina, i.e. top and bottom lamina are of 0º. Results of non-dimensional deflection for a thick ply with span to thickness ratio of 4 are presented in Table 2. Results are presented with various degrees of imperfection in the shear direction, whereas the normal imperfection is neglected, i.e.

. The presented results are

compared with the exact solutions provided by the Chen et al. [41] as well as the results of the Pagano [55] for perfectly bonded laminates are also included. Relative error calculated with the results of Chen et al. is shown within the table. The present formulation can predict the results within 1% error for such a thick ply with six interfaces. Table 2: Non-dimensional deflection for various extents of shear imperfections ((0/90)3/0) laminate Stacking sequence

l/h

Reference

(0/90)3/0

4

SBFEM Chen et al. [41] Pagano [55]

Non-dimensional deflection 3.1247 3.1110 3.1110

24

4.0409 4.0361 -

4.8965 4.9184 -

5.7566 5.7612 -

for

Avg. % Difference 0.27 -

4.5

Two layered anti-symmetric cross ply (90/0) laminate

The presented formulation is also applied for the two layered antisymmetric laminate. Fibre orientation in the top lamina is at 90º and in bottom one is 0º, i.e. the load is applied on the 90º lamina. A thick laminate plate with

is analysed to obtain longitudinal stress and

shear stress. z/h

z/h

SBFEM (Rn = 0) 0.50 SBFEM (Rn = 100) SBFEM (Rn = 1x10⁶) Liu et al. (Rn = 0) 0.25 Liu et al. (Rn = 100) Liu et al. (Rn = 1x10⁶)

-50

-30

0.50

0.25

Rt = 5x10-5

Rt = 5x104

̅σxx

0.00 -10

10

30

̅σxx -50

-30

0.00 -10

-0.25

-0.25

-0.50

-0.50

(a)

10

30

50

(b)

Figure 13: Variation of ̅ along the non-dimensional thickness for different degrees of imperfections in normal direction in (90/0) laminate with l/h=4 (a) (b) in comparison to results of Liu et al. [29] Through the thickness variation of longitudinal ̅

and shear stress

̅

are depicted in

Figure 13 and Figure 14 respectively. Results are presented for different values of normal separation and tangential slip coefficients. All the results are compared with the results of layer-wise model presented by the Liu et al. [29]. Results of the present SBFEM formulation are consistent with the available results in the literature. It can be noted from the figures that longitudinal stress always experiences a jump at the interface whereas the shear stress profile is always continuous, irrespective of the degree of imperfections. 5.

Conclusions

In this paper, a relatively new semi-analytical approach of scaled boundary finite element method (SBFEM) is employed to study the cylindrical bending of laminated composite plates with interlaminar bonding imperfections. These imperfections develop displacement jumps at the interfaces which are simulated using the spring layer model i.e. a linear traction separation rule. A plane strain formulation is adopted for the modelling hence, there is not a 25

requirement of making any assumptions on the displacement or stress fields like plate theories. The ability of SBFEM to easily handle hanging nodes helps to incorporate the small interface elements without reducing the size of adjoining regular elements. This type of meshing strategy enhances the computational efficiency as compared to traditional FEM. z/h

SBFEM (Rn = 0)

Rt = 5x10-5

0.5

SBFEM (Rn = 100) SBFEM (Rn = 1x10⁶)

0.3

Liu et al. (Rn = 0)

Liu et al. (Rn = 100)

0.1

̅τxz

Liu et al. (Rn = 1x10⁶)

(a) -1.5

-0.1

0.5

2.5

4.5

-0.3

-0.5

z/h 0.5

0.3

SBFEM (Rn = 0)

Rt = 5x104

SBFEM (Rn = 100)

0.1

̅τxz

(b)

SBFEM (Rn = 1x10⁶) Liu et al. (Rn = 0) Liu et al. (Rn = 100)

-1.0

-0.1

0.0

1.0

2.0

3.0

4.0

Liu et al. (Rn = 1x10⁶)

-0.3

-0.5

Figure 14: Variation of ̅ along the non-dimensional thickness for different degrees of imperfections in normal direction in (90/0) laminate with l/h=4 (a) (b) in comparison to results of Liu et al. [29] Numerical examples have been solved to study the effect of span to thickness ratio, lamination sequence, number of layers, layer thickness and degree of imperfection in normal and tangential directions on global response of the cross-ply laminates and on interlaminar

26

stresses. From the case studies, following conclusions can be drawn about the potential of the novel formulation, that is: a. SBFEM is able to accommodate shear slip as well as general delamination considering both shear and normal separation. b. SBFEM predicts results with very good accuracy, i.e. within 1% error. c. The number of lamina or interface count does not affect the accuracy of SBFEM, as reflected by the results obtained for ((0/90)3/0) laminate. d. SBFEM can be a suitable predictor for the case of anti-symmetric plies. Hence, it is observed from the results that the proposed formulation is able to predict the results accurately, even more accurate than previously developed highly efficient zigzag theories. The plane strain formulation along with the semi-analytical nature of the SBFEM leads to the accuracy of the results. Furthermore, the adopted meshing strategy makes it computationally cheaper. Hence, the proposed model is recommended as a very accurate and efficient model to predict the mechanical response of cross ply laminated composite plates with imperfect interfaces undergoing cylindrical bending. As well as being applicable to cylindrical bending the plane strain formulation used here is also applicable for the bending of beams [56]. Future work will include extending SBFEM for 3D elements which will be able to incorporate the more general layup sequence with angle plies. Author statement Nikhil Garg – Conceptualization, Methodology, Software, Writing - Original Draft, Visualization Nilanjan Das Chakladar – Methodology, Supervision, Writing - Review & Editing B. Gangadhara Prusty – Supervision, Writing - Review & Editing, Project administrator, Funding acquisition Chongmin Song – Resources, Supervision, Writing - Review & Editing Andrew Phillips – Supervision, Writing - Review & Editing, Funding acquisition

Declaration of interests

27

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement This project was conducted within the ARC Training Centre for Automated Manufacture of Advanced Composites (IC160100040), supported by the Commonwealth of Australia under the Australian Research Council’s Industrial Transformation Research Program. The authors express their regards and gratitude to Dr. Garth Pearce (School of Mechanical and Manufacturing Engineering, UNSW, Sydney) for his valuable suggestions while conducting the research work. References [1] [2] [3] [4]

[5] [6] [7] [8]

[9]

[10] [11] [12] [13] [14]

E. Reissner, On a variational theorem in elasticity, J. Math. Phys. 29 (1950) 90–95. R.D. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates, J. Appl. Mech. 18 (1951) 31–38. J.N. Reddy, A Simple Higher-Order Theory for Laminated Composite Plates, J. Appl. Mech. (1984). doi:10.1115/1.3167719. N. Grover, B.N. Singh, D.K. Maiti, New nonpolynomial shear-deformation theories for structural behavior of laminated-composite and sandwich plates, AIAA J. 51 (2013) 1861–1871. J.G. Ren, Bending theory of laminated plate, Compos. Sci. Technol. 27 (1986) 225– 248. H. Murakami, Laminated composite plate theory with improved in-plane responses, J. Appl. Mech.(Trans. ASME). 53 (1986) 661–666. S. Kapuria, G.G.S. Achary, An efficient higher order zigzag theory for laminated plates subjected to thermal loading, Int. J. Solids Struct. 41 (2004) 4661–4684. H.D. Chalak, A. Chakrabarti, M.A. Iqbal, A.H. Sheikh, An improved C0 FE model for the analysis of laminated sandwich plate with soft core, Finite Elem. Anal. Des. 56 (2012) 20–31. N. Garg, R.S. Karkhanis, R. Sahoo, P.R. Maiti, B.N. Singh, Trigonometric zigzag theory for static analysis of laminated composite and sandwich plates under hygrothermo-mechanical loading, Compos. Struct. 209 (2019) 460–471. R. Sahoo, B.N. Singh, A new inverse hyperbolic zigzag theory for the static analysis of laminated composite and sandwich plates, Compos. Struct. 105 (2013) 385–397. R. Khandan, S. Noroozi, P. Sewell, J. Vinney, The development of laminated composite plate theories: a review, J. Mater. Sci. 47 (2012) 5901–5910. E. Carrera, Historical review of zig-zag theories for multilayered plates and shells, Appl. Mech. Rev. 56 (2003) 287–308. Y.X. Zhang, C.H. Yang, Recent developments in finite element analysis for laminated composite plates, Compos. Struct. 88 (2009) 147–157. Y.M. Ghugal, R.P. Shimpi, A review of refined shear deformation theories for 28

[15] [16] [17]

[18]

[19]

[20]

[21]

[22]

[23]

[24] [25] [26] [27] [28] [29]

[30] [31] [32] [33] [34]

isotropic and anisotropic laminated beams, J. Reinf. Plast. Compos. 20 (2001) 255– 272. R. Talreja, Manufacturing defects in composites and their effects on performance, in: Polym. Compos. Aerosp. Ind., Elsevier, 2015: pp. 99–113. K.P. Soldatos, X. Shu, Modelling of perfectly and weakly bonded laminated plates and shallow shells, Compos. Sci. Technol. 61 (2001) 247–260. C.K. Hirwani, S.K. Panda, S.S. Mahapatra, S.K. Mandal, L. Srivastava, M.K. Buragohain, Flexural strength of delaminated composite plate–An experimental validation, Int. J. Damage Mech. 27 (2018) 296–329. C.K. Hirwani, S.K. Panda, T.R. Mahapatra, S.S. Mahapatra, Numerical study and experimental validation of dynamic characteristics of delaminated composite flat and curved shallow shell structure, J. Aerosp. Eng. 30 (2017) 4017045. C.K. Hirwani, S.K. Panda, T.R. Mahapatra, Thermomechanical deflection and stress responses of delaminated shallow shell structure using higher-order theories, Compos. Struct. 184 (2018) 135–145. C.K. Hirwani, S.K. Panda, Numerical and experimental validation of nonlinear deflection and stress responses of pre-damaged glass-fibre reinforced composite structure, Ocean Eng. 159 (2018) 237–252. C.K. Hirwani, S.K. Panda, T.R. Mahapatra, S.K. Mandal, S.S. Mahapatra, A.K. De, Delamination effect on flexural responses of layered curved shallow shell panelexperimental and numerical analysis, Int. J. Comput. Methods. 15 (2018) 1850027. S.S. Sahoo, S.K. Panda, T.R. Mahapatra, C.K. Hirwani, Numerical analysis of transient responses of delaminated layered structure using different mid-plane theories and experimental validation, Iran. J. Sci. Technol. Trans. Mech. Eng. 43 (2019) 41–56. C.K. Hirwani, S.K. Panda, Nonlinear finite element solutions of thermoelastic deflection and stress responses of internally damaged curved panel structure, Appl. Math. Model. 65 (2019) 303–317. A. Toledano, H. Murakami, Shear-deformable two-layer plate theory with interlayer slip, J. Eng. Mech. 114 (1988) 604–623. H. Murakami, A laminated beam theory with interlayer slip, J. Appl. Mech. 51 (1984) 551–559. X. Lu, D. Liu, Interlayer shear slip theory for cross-ply laminates with nonrigid interfaces, AIAA J. 30 (1992) 1063–1073. T.O. Williams, F.L. Addessio, A general theory for laminated plates with delaminations, Int. J. Solids Struct. (1997). doi:10.1016/S0020-7683(96)00131-X. J. Aboudi, Damage in composites—modeling of imperfect bonding, Compos. Sci. Technol. 28 (1987) 103–128. D. Liu, L. Xu, X. Lu, Stress analysis of imperfect composite laminates with an interlaminar bonding theory, Int. J. Numer. Methods Eng. (1994). doi:10.1002/nme.1620371608. Z. Cheng, D. Kennedy, W.F. Williams, Effect of interfacial imperfection on buckling and bending behavior of composite laminates, AIAA J. 34 (1996) 2590–2595. Z. Cheng, A.K. Jemah, F.W. Williams, Theory for multilayered anisotropic plates with weakened interfaces, J. Appl. Mech. 63 (1996) 1019–1026. R. Massabo, F. Campi, Assessment and correction of theories for multilayered plates with imperfect interfaces, Meccanica. 50 (2015) 1045–1071. R. Massabò, F. Campi, An efficient approach for multilayered beams and wide plates with imperfect interfaces and delaminations, Compos. Struct. 116 (2014) 311–324. R. Massabò, Influence of boundary conditions on the response of multilayered plates with cohesive interfaces and delaminations using a homogenized approach, Frat. Ed 29

Integrità Strutt. 8 (2014) 230–240. [35] J.-S. Kim, J. Oh, M. Cho, Efficient analysis of laminated composite and sandwich plates with interfacial imperfections, Compos. Part B Eng. 42 (2011) 1066–1075. [36] Z.-Q. Cheng, L.-H. He, S. Kitipornchai, Influence of imperfect interfaces on bending and vibration of laminated composite shells, Int. J. Solids Struct. 37 (2000) 2127– 2150. [37] M. Di Sciuva, Geometrically Nonlinear Theory of Multilayered Plates with Interlayer Slips, AIAA J. (1997). doi:10.2514/2.23. [38] M. Di Sciuva, U. Icardi, L. Librescu, Effects of interfacial damage on the global and local static response of cross-ply laminates, Int. J. Fract. (1999). doi:10.1023/A:1018649824301. [39] Z.-Q. Cheng, W.P. Howson, F.W. Williams, Modelling of weakly bonded laminated composite plates at large deflections, Int. J. Solids Struct. 34 (1997) 3583–3599. [40] M. Di Sciuva, U. Icardi, L. Librescu, On modeling of laminated composite structures featuring interlaminae imperfections, in: Stud. Appl. Mech., Elsevier, 1997: pp. 395– 404. [41] W.Q. Chen, J.B. Cai, G.R. Ye, Exact Solutions of Cross-Ply Laminates with Bonding Imperfections, AIAA J. (2003). doi:10.2514/2.6817. [42] W.Q. Chen, K.Y. Lee, Three-dimensional exact analysis of angle-ply laminates in cylindrical bending with interfacial damage via state-space method, Compos. Struct. 64 (2004) 275–283. [43] W.Q. Chen, J. Ying, J.B. Cai, G.R. Ye, Benchmark Solution of Laminated Beams with Bonding Imperfections, AIAA J. (2004). doi:10.2514/1.4776. [44] V.Q. Bui, E. Marechal, H. Nguyen-Dang, Imperfect interlaminar interfaces in laminated composites: bending, buckling and transient reponses, Compos. Sci. Technol. 59 (1999) 2269–2277. [45] V.Q. Bui, E. Marechal, H. Nguyen-Dang, Imperfect interlaminar interfaces in laminated composites: interlaminar stresses and strain-energy release rates, Compos. Sci. Technol. 60 (2000) 131–143. [46] A. Chakrabarti, A.H. Sheikh, Behavior of laminated sandwich plates having interfacial imperfections by a new refined element, Comput. Mech. 34 (2004) 87–98. [47] C. Song, J.P. Wolf, The scaled boundary finite-element method - Alias consistent infinitesimal finite-element cell method - For elastodynamics, Comput. Methods Appl. Mech. Eng. (1997). doi:10.1016/S0045-7825(97)00021-2. [48] C. Song, J.P. Wolf, Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method, Comput. Struct. (2002). doi:10.1016/S0045-7949(01)00167-5. [49] C. Song, F. Tin-Loi, W. Gao, A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges, Eng. Fract. Mech. (2010). doi:10.1016/j.engfracmech.2010.04.032. [50] E.T. Ooi, C. Song, F. Tin-Loi, Z. Yang, Polygon scaled boundary finite elements for crack propagation modelling, Int. J. Numer. Methods Eng. (2012). doi:10.1002/nme.4284. [51] C. Song, E.T. Ooi, S. Natarajan, A review of the scaled boundary finite element method for two-dimensional linear elastic fracture mechanics, Eng. Fract. Mech. (2018). doi:10.1016/j.engfracmech.2017.10.016. [52] C. Song, The Scaled Boundary Finite Element Method: Introduction to Theory and Implementation, John Wiley & Sons, 2018. [53] J.A. Nairn, Numerical implementation of imperfect interfaces, Comput. Mater. Sci. (2007). doi:10.1016/j.commatsci.2007.02.010. 30

[54] M. Alfano, F. Furgiuele, A. Leonardi, C. Maletta, G.H. Paulino, Fracture analysis of adhesive joints using intrinsic cohesive zone models, in: Convegno IGF XIX Milano 2007, 2007. [55] N.J. Pagano, Exact solutions for composite laminates in cylindrical bending, J. Compos. Mater. 3 (1969) 398–411. [56] V.A. Jairazbhoy, P. Petukhov, J. Qu, Large deflection of thin plates in cylindrical bending–Non-unique solutions, Int. J. Solids Struct. 45 (2008) 3203–3218.

Graphical Abstract

31