Global stiffness determination of a cross-ply laminated composite plate with distributed delaminations and matrix cracks and its application

Journal Pre-proofs Global stiffness determination of a cross-ply laminated composite plate with distributed delaminations and matrix cracks and its application In-Bong Kim, U. Song-Hak, Su-Hyok Sin PII: DOI: Reference:

S0263-8223(19)30765-2 https://doi.org/10.1016/j.compstruct.2019.111586 COST 111586

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

1 March 2019 16 September 2019 21 October 2019

Please cite this article as: Kim, I-B., Song-Hak, U., Sin, S-H., Global stiffness determination of a cross-ply laminated composite plate with distributed delaminations and matrix cracks and its application, Composite Structures (2019), doi: https://doi.org/10.1016/j.compstruct.2019.111586

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.

Revised manuscript

Global stiffness determination of a cross-ply laminated composite plate with distributed delaminations and matrix cracks and its application In-Bong Kim, Song-Hak U*, Su-Hyok Sin Department of Mechanics of Composite Materials, Faculty of Mechanics, Kim Il Sung University, Democratic People’s Republic of Korea *

Corresponding author Email: [email protected]

Abstract This paper is focused on a global stiffness method for cross-ply laminated composite plates with distributed delaminations and matrix cracks that makes it possible to find its solution by simply replacing stiffness quantities in the known solution for the undamaged counterpart. Local and zero stiffness model are introduced to determine the local stiffness in damaged regions and using the energy functional expression based on the first-order shear deformation theory(FSDT) and the classic laminate theory(CLT), the global stiffness is analytically determined by combining stiffness terms in undamaged and damaged regions with corresponding area ratio weights. The obtained global stiffness is used to analyse bending and buckling problem of a simply supported plate with distributed delaminations and matrix cracks subjected to in-plane and transverse loads, being validated through ANSYS simulation.

Keywords Global stiffness, Delamination, Matrix crack, Laminated composite plate, First-order shear deformation theory

1. Introduction Because of its excellent physical and mechanical properties, plates and shells fabricated by composite materials are widely used in many fields such as aerospace, turbine, wind turbine blades, nuclear reactor, etc. However, those are prone to delaminations, cracks and defects during manufacturing as well as working. Therefore, it is meaningful for reliable and safe designs to build a proper model taking these factors into account. Recently, the mechanical behaviour of composite plates with delaminations and matrix cracks, or shortly DMCs, is being widely researched. In Allam and Zenkour [1], researchers analysed a composite plate weakened by an oval opening and deformed by forces applied to the middle plane. The stress concentration factor for a fiber-reinforced viscoelastic plate is given using the method of effective moduli. Also, the numerical results on stress concentration factor are given for the problems of tension and pure bending in the various directions of a structurally anisotropic composite plate. 1

In some references including Wang and Tong [2] and Parlapalli et al. [3], the plate or plate-beam with embedded delaminations of certain length is divided into the undamaged region and delaminated region which is considered by dividing it into upper and lower subregions having delamination as boundary. And then the beam theory is applied in each subregion to derive basic equations, regarding the continuous condition at adjacent boundary, which yield the analytical solutions for vibration, bending and buckling problem. In most practical problems, the delamination does not exist individually, but in the combined form of transverse cracks and delamination, especially it is the case when structures are subjected to bending load. Transverse crack fracture is modeled by using continuous damage and plasticity model. Ahmed and Sluys [4] pointed that a drawback of this approach was such that, even though delamination damage was modeled using a discrete fracture approach, the use of a plasticity model for intra-ply damage can still cause problems due to strain localization, resulting in ill-posedness of the governing equations and mesh dependency. They developed the computation framework, considering the influence between DMCs in composite plates under transverse loading. In Kaya et al. [5] damage mechanisms such as debonding, fiber fracture, delamination and matrix cracking within composite plates subjected to tensile loading are analysed using acoustic emission technique. In Zhang et al. [6] an approach to predict the initiation and propagation of damage in laminated composite plates is brought forward by penalty function method. And the potential delamination and matrix cracking areas are considered as cohesive zone and the damage process as contact behaviour between the interfaces. On the other hand, Yang et al. [7] established a multiscale numerical approach to model damages in random glass fiber composites. A representative volume element of a random glass fiber composite is employed to analyze microscale damage mechanisms, such as matrix cracking and fiber-matrix interfacial debonding. In Lee et al. [8] a nonlinear finite element code, DELAM3D, with a three-dimensional layered solid element based on the updated Lagrangian formulation, is developed to simulate the compressive response of a laminated composite plate with multiple delaminations. An analytical model is established to characterize the mechanical behaviour such as postbuckling, contact of the delaminating interface, delamination growth and fiber-matrix failure. In Yan and Yam [9] are studied the detection of crack damage in composite plates using the embedded piezoelectric materials and dynamic responses by wavelet analysis. In this case the changed structural elasticity parameters presented in Talreja [10] for cracks perpendicular to and parallel with the fiber are used. These stiffness parameters can be written in terms of crack damage variable correlating to crack density. Nakamura and Wu [11] studied the effect of delamination and ply-arrangement on multilayered structures. Energy release rate and mixed-mode stress intensity factors are computed to quantify the crack driving force. The results are used to determine dominant failure initiation mode, structural buckling or delamination growth. Cheng et al. [12] devoted to modeling of elastic behaviour of laminated composite shells, with special emphasis on incorporating interfacial imperfection. The conditions of imposing 2

traction continuity and displacement jump across each interface are used to model imperfect interfaces. In Czarnocki [13] the analysis is focused on the extent to which the thermal loading resulting from a difference between the curing and service temperatures and from variations in the thermal expansion coefficients from layer to layer can affect the growth of delaminations. Square laminate plates containing circular delaminations were examined. Wang and Lin [14] considered the buckling of plate with delamination under cylindrical bending, modeling it as plate-beam having delamination that passes through it along width. However, preceding researches focused only the case in which delamination passed through the width of plate-beam. The plate with delamination that does not reach the four sides of plate and embeded inside plate would be investigated by FEM ([15]). Kim and Keith [16] proposed the analytical method for predicting the onset of buckling in delaminated composite plates, with which both global and local sublaminates buckling deformation modes were predicted for rectangular plates of any thickness. Chattopadhyay et al. [17] and Haozhong et al. [18] developed new high-order theories and applied to research of delamination buckling and postbuckling in the cylindrical composite shells. They concluded that the transverse shear effect in delamination buckling is of importance. In Zhu et al. [19] the reference surface element formulation of a four-node quadrilateral membrane-shear-bending element (ZQUA24) was presented and numerical investigations were performed for composite plates and shells with various delamination shapes. The numerical results show that the technique is simple, reliable and able to model delamination buckling and postbuckling behaviour of laminated plates or shells. Elisa Pietropaoli and Aniello Riccio [20] and Endel V. Iarve et al [21] presented methods base on FEM, e.g. x-FEM and a global/local FEA, respectively. D. H. Li et al [22] established an extended layerwise method (XLWM) for the composite laminated beams with multiple DMCs. A new displacement field in the thickness direction was proposed, which contains the linear Lagrange interpolation functions, the one-dimensional weak discontinuous function and strong discontinuous function to model the displacement discontinuity and the strain discontinuity. Summarizing the previous researches above, it is drawn that laminated composite plates with DMCs were analyzed by introducing FE models such as the x-FEM, the global/local FEA and the reference surface element combined with sublayer FEM while distributed crack damages were treated in terms of crack damage variable according to crack density. In this paper we present the global stiffness of cross-ply laminated composite plates with distributed DMCs, which makes it possible to find its solution by simply replacing stiffness quantities in the known solution for the undamaged counterpart. Based on the energy functional expression in the FSDT and CLT, the global stiffness quantities are analytically determined by combining stiffness terms in undamaged and damaged regions with corresponding area ratio weights. Local and zero stiffness models are introduced to determine the local stiffness in damaged regions. We use the global stiffness quantities obtained to analyse the bending and buckling 3

problem of a simply supported plate with distributed DMCs subjected to in-plane and transverse load and validate the present method through ANSYS simulation. 2. Global stiffness determination of a cross-ply laminated composite plate with distributed DMCs Consider a cross ply laminated plate with sides 𝑎, 𝑏 and thickness ℎ, subjected to inplane loads 𝑃𝑥,𝑃𝑦 and transverse distributed load 𝑞. It is assumed that the plate has several distributed DMCs and their interrelations between neighbouring DMCs can be ignored. Besides, each delamination embedded in the plate is square-shaped and matrix crack that is vertical to it may be placed in one layer or through successive two or more layers, cutting both matrix and fibers. Divide the plate into undamaged and damaged regions according to distributed DMCs and take a cubic element containing the eth DMC tightly as the eth damaged region (Fig.1). Obviously, the length and width of element are equal to the delamination sides 𝑎𝑘 × 𝑏𝑘 that lies at the interface between the 𝑘th and (𝑘 + 1)th layer of plate, while its thickness is equal to the plate thickness itself. The matrix crack of length by depth equal to 𝑙𝑘 × 𝑑𝑘 stands vertically on the delamination and its length is towards y-axis.

Fig. 1. The laminated plate with DMC, its section and the eth damaged (grey) region Let the stiffness of element in damaged region and that of total plate be called the local stiffness and global stiffness, respectively. Quantities that are relevant to damaged regions are distinguished by the superscript “e”. To start with, consider the stiffness determination in damaged regions. For this, we introduce two stiffness models, i.e., the local and zero stiffness model. 2.1. The local and zero stiffness model for damaged layers The first one is a local stiffness model which uses the local stiffness of 𝑘th cracked layer in the damaged region ([23]). The local stiffness constants of 𝑘th layer having DMC in the eth damaged region are (1) (𝐶𝑒𝑖𝑗)𝑘 = (𝐶𝑖𝑗)𝑘 ― (∆𝐶𝑒𝑖𝑗)𝑘 ,(𝑖,𝑗 = 1,2,4,5,6) 𝑒 , where (𝐶𝑖𝑗)𝑘 are the elastic stiffness constants without DMCs and ∆𝐶𝑖𝑗 are the incremental elastic stiffness constants attributed to DMCs. The latter constants are bounded as 0 ≤ (∆𝐶𝑒𝑖𝑗)𝑘 ≤ (𝐶𝑖𝑗)𝑘, being equal to zero (∆𝐶𝑒𝑖𝑗)𝑘 = 0 for no DMC and (∆𝐶𝑒𝑖𝑗)𝑘 = (𝐶𝑖𝑗)𝑘 for the 4

completely damaged element layer by crack growth. The transformed reduced stiffnesses (𝑄𝑒𝑖𝑗)𝑘 corresponding to (𝐶𝑒𝑖𝑗)𝑘 are given as (2) (𝑄𝑒𝑖𝑗)𝑘 = (𝑄𝑖𝑗)𝑘 ― (∆𝑄𝑒𝑖𝑗)𝑘 𝑒 , where (𝑄𝑖𝑗)𝑘 and (∆𝑄𝑖𝑗)𝑘 are the transformed reduced stiffnesses corresponding to (𝐶𝑖𝑗)𝑘 and (∆𝐶𝑒𝑖𝑗)𝑘, respectively. Especially, (∆𝑄𝑒𝑖𝑗)𝑘 are those for the 𝑘th layer in damaged region and determined by the expression ∆𝐶𝑒𝑖𝑗 = 𝑓(𝐴,𝐹,𝐺,𝜑𝑖𝑗,ℎ𝑒), where 𝐴,𝐹,𝐺 - the constants related to material properties, 𝜑𝑖𝑗 - the stiffness relations in the crack coordinate system, ℎ𝑒 - the thickness regarding stress concentration at the crack tip. The second one is a zero stiffness model which ignores the stiffness of the 𝑘th cracked layer in damaged region. The local stiffness model based on (∆𝐶𝑒𝑖𝑗)𝑘 reflects enough the local status of matrix crack, but evaluating (∆𝐶𝑒𝑖𝑗)𝑘 is rather complicated. The zero stiffness model is based on experimental observation ([24]) on damage state of the 𝑘th cracked layer in damaged region. Now, determine the transformed reduced stiffnesses (𝑄𝑒𝑖𝑗)𝑘 (Eq. (2)) of damaged region without using (∆𝐶𝑒𝑖𝑗)𝑘. If a matrix crack grows, producing fiber bridging and fiber/matrix ruptures, then bending and shear resistance of the cracked layer gets decreased remarkably. As seen often in practice, due to repeated or impulsive transverse loading, damages including DMC are initiated and grows through the layer thickness, resulting in fiber bridging, fiber rupture, fiber pull-out and fiber/matrix debonding. Of course, these damages may be initiated by physical-chemical actions. Then, the damaged layer is not able to resist tension/compression as well as bending/shear so that the layer stiffness gets reduced to zero. Therefore, for the 𝑘th cracked layer of damaged region, the following is assumed. (3) (𝑄𝑒𝑖𝑗)𝑘 = 0 This equation means the zero local stiffness of damaged layer, i.e., the complete damage of local part. 2.2. Global stiffness determination based on the FSDT The global stiffness can be obtained by finding stiffnesses in undamaged and damaged regions and combining those with a weight equal to the ratio of each area to total plate area. Here, the stiffnesses refer to extensional 𝐴𝑖𝑗, coupling 𝐵𝑖𝑗 and bending 𝐷𝑖𝑗 stiffness terms while the area means the corresponding middle plane area. The stiffnesses of damaged region are determined by the local and zero stiffness model introduced above. For the sake of brevity, consider the case of symmetric cross ply laminated plates. In such case, 𝐴16 = 𝐴26 = 𝐷16 = 𝐷26 = 0 and [𝐵] = 0. According to the FSDT, the followings are obtained for this case. i) The local stiffness model 𝑁

𝐴𝑒𝑝𝑞 = ∑𝑘 = 1(𝑄𝑒𝑝𝑞)𝑘(𝑧𝑘 ― 𝑧𝑘 ― 1) 1 𝑁

𝐷𝑒𝑖𝑗 = 3∑𝑘 = 1(𝑄𝑒𝑖𝑗)𝑘(𝑧3𝑘 ― 𝑧3𝑘 ― 1)

for the damaged region

(4)

for the undamaged region

(5)

𝑁

𝐴𝑝𝑞 = ∑𝑘 = 1(𝑄𝑝𝑞)𝑘(𝑧𝑘 ― 𝑧𝑘 ― 1) 1 𝑁

𝐷𝑖𝑗 = 3∑𝑘 = 1(𝑄𝑖𝑗)𝑘(𝑧3𝑘 ― 𝑧3𝑘 ― 1)

5

(𝑁- total lay number, 𝑝𝑞 = 11, 12, 22, 44, 55, 66; 𝑖𝑗 = 11, 12, 22, 66) It is noted that, in Eq. (4), the stiffnesses of damaged layer (𝑄𝑒𝑖𝑗)𝑘 are determined from local stiffness quantities (𝐶𝑒𝑖𝑗)𝑘, but those of undamaged layer in the damaged region are the same as one for the undamaged region. Thus, denoting the damaged layer number in the eth damaged region by 𝛾(𝑒) gives (𝑄𝑒𝑖𝑗)𝑘

{

(𝑄𝑒𝑖𝑗)𝛾 :𝑘 = 𝛾(𝑒) = (𝑄𝑖𝑗)𝑘 :𝑘 ≠ 𝛾(𝑒)

(6)

ii) The zero stiffness model According to the zero stiffness model, equating the stiffness of damaged layer to zero and regarding that the stiffness of undamaged layer is the same as one in the undamaged region produces the following. (𝑄𝑒𝑖𝑗)𝑘 = (1 ― 𝛿(𝑒) 𝑘𝛾 ) ∙ (𝑄𝑖𝑗)𝑘,

𝛿(𝑒) 𝑘𝛾 =

{

1:𝑘 = 𝛾(𝑒) 0:𝑘 ≠ 𝛾(𝑒)

(7)

Hence, the stiffness quantities in the damaged region are given as 𝑁

𝐴𝑒𝑝𝑞 = ∑𝑘 = 1(𝑄𝑝𝑞)𝑘(1 ― 𝛿(𝑒) 𝑘𝛾 )(𝑧𝑘 ― 𝑧𝑘 ― 1) = ∑𝑘 ≠ 𝛾(𝑄𝑝𝑞)𝑘(𝑧𝑘 ― 𝑧𝑘 ― 1) 1 𝑁

1

3 3 3 3 𝐷𝑒𝑖𝑗 = 3∑𝑘 = 1(𝑄𝑖𝑗)𝑘(1 ― 𝛿(𝑒) 𝑘𝛾 )(𝑧𝑘 ― 𝑧𝑘 ― 1) = 3∑𝑘 ≠ 𝛾(𝑄𝑖𝑗)𝑘(𝑧𝑘 ― 𝑧𝑘 ― 1)

(8)

Those in the undamaged region are the same as one (Eq. (5)) obtained by the local stiffness model and not listed again. Also, the subscripts 𝑝𝑞, 𝑖𝑗 take the same values as those in Eq. (5), otherwise pointed. As mensioned, the global stiffness of plate can be obtained by combining stiffnesses in undamaged and damaged regions with corresponding area ratio weights. Denoting the area of plate middle surface by 𝑆, the areas of undamaged and damaged regions by 𝑆𝑢 and 𝑆𝑑, the area of 𝑒th damaged region by 𝑆𝑒 and the number of damaged regions by 𝑀, respectively gives 𝑀

(9)

𝑆𝑑 = ∑𝑒 = 1𝑆𝑒, 𝑆𝑢 + 𝑆𝑑 = 𝑆 According to the FSDT, the displacements are represented by 𝑢(𝑥,𝑦) = 𝑢0(𝑥,𝑦) +𝑧𝜓𝑥(𝑥,𝑦), ν(𝑥,𝑦) = 𝜈0(𝑥,𝑦) +𝑧𝜓𝑦(𝑥,𝑦), 𝑤 = 𝑤0(𝑥,𝑦)

, where 𝜓𝑥,𝜓𝑦 are the rotations of planes that were perpendicular to the undeformed middle plane about 𝑦,𝑥 axis, respectively. The energy functional of cross ply laminated plates is used to derive the global stiffness. The total energy of plate is expressed in terms of the strain energy 𝑈 and the energy of external loads 𝑉 as ([25]) (10.a) 𝛱=𝑈+𝑉 1

1

𝑀

𝑈 = 2∬𝑠𝛬(𝐴𝑖𝑗,𝐷𝑖𝑗)𝑑𝑥𝑑𝑦 = 2(∬𝑆 𝛬𝑢𝑑𝑥𝑑𝑦 + ∑𝑒 = 1∬𝑆𝑒𝛬𝑒𝑑𝑥𝑑𝑦) 𝑢

1

𝑉 = ― 2∬𝑠[𝑃𝑥(𝑤,𝑥)2 + 2𝑃𝑥𝑦𝑤,𝑥𝑤,𝑦 + 𝑃𝑦(𝑤,𝑦)2 + 2𝑞𝑤]𝑑𝑥𝑑𝑦 6

(10.b) (10.c)

, where 3

3

𝛬𝑢 = ∑𝑠 = 1𝛬𝑢𝑠(𝐴𝑖𝑗,𝐷𝑖𝑗), 𝛬𝑒 = ∑𝑠 = 1𝛬𝑒𝑠(𝐴𝑒𝑖𝑗,𝐷𝑒𝑖𝑗)

(11)

𝛬𝑢1 = 𝐴11(𝑢0,𝑥)2 + 2𝐴12𝑢0,𝑥𝑣0,𝑦 + 𝐴22(𝑣0,𝑦)2 + 𝐴66(𝑢0,𝑦 + 𝑣0,𝑥)2 𝛬𝑢2 = 𝐷11(𝜓𝑥,𝑥)2 + 2𝐷12𝜓𝑥,𝑥𝜓𝑦,𝑦 + 𝐷22(𝜓𝑦,𝑦)2 + 𝐷66(𝜓𝑥,𝑦 + 𝜓𝑦,𝑥)2 𝛬𝑢3 = 𝑙𝐴55(𝑤,𝑥 + 𝜓𝑥)2 + 𝑙𝐴44(𝑤,𝑦 + 𝜓𝑦)2 𝛬𝑒1 = 𝐴𝑒11(𝑢0,𝑥)2 + 2𝐴𝑒12𝑢0,𝑥𝑣0,𝑦 + 𝐴𝑒22(𝑣0,𝑦)2 + 𝐴𝑒66(𝑢0,𝑦 + 𝑣0,𝑥)2 𝛬𝑒2 = 𝐷𝑒11(𝜓𝑥,𝑥)2 + 2𝐷𝑒12𝜓𝑥,𝑥𝜓𝑦,𝑦 + 𝐷𝑒22(𝜓𝑦,𝑦)2 + 𝐷𝑒66(𝜓𝑥,𝑦 + 𝜓𝑦,𝑥)2 𝛬𝑒3 = 𝑙𝐴𝑒55(𝑤,𝑥 + 𝜓𝑥)2 + 𝑙𝐴𝑒44(𝑤,𝑦 + 𝜓𝑦)2

(12)

In Eq. (12), 𝑙 = 5/6 is the shear correction factor. It is well known that the variation method is available using total energy Eq. (10). Alternatively, introducing the global stiffness based on area ratio weight, it is possible to find solution easily and simply for plates with DMC by using the known solution for plates without DMC. As seen in Eq. (11) and Eq. (12), the energy density function 𝛬𝑢,𝛬𝑒 are continuous and linear with respect to stiffness components 𝐴𝑖𝑗, 𝐷𝑖𝑗, 𝐴𝑒𝑖𝑗, 𝐷𝑒𝑖𝑗. Regarding Eq.(12), represent 𝛬𝑢𝑠 and 𝛬𝑒𝑠 in Eq.(11) as the following form. 𝐴𝑖𝑗 :𝑠 = 1, 3 (13) 𝛬𝑢𝑠 = ∑𝑟𝐹𝑖𝑗𝛬𝑢𝑠𝑟, 𝛬𝑒𝑠 = ∑𝑟𝐹𝑒𝑖𝑗𝛬𝑒𝑠𝑟 , 𝐹𝑖𝑗 = 𝐷 :𝑠 = 2 𝑖𝑗

{

,where 𝛬𝑢𝑠𝑟,𝛬𝑒𝑠𝑟 are coefficients of 𝐴𝑖𝑗,𝐷𝑖𝑗, 𝐴𝑒𝑖𝑗, 𝐷𝑒𝑖𝑗 in 𝛬𝑢𝑠,𝛬𝑒𝑠 and 𝑟 is a summation index. Taking consideration of Eq. (13) and the continuity of terms in the corresponding regions, the strain energy expression is as follows.

(

1

3

𝑀

3

𝑈 = 2 ∬𝑆 ∑𝑠 = 1∑𝑟𝐹𝑖𝑗𝛬𝑢𝑠𝑟𝑑𝑥𝑑𝑦 + ∑𝑒 = 1∬𝑆𝑒∑𝑠 = 1∑𝑟𝐹𝑒𝑖𝑗𝛬𝑒𝑠𝑟𝑑𝑥𝑑𝑦 1

𝑢

3

)

𝑀

(14)

= 2∑𝑠 = 1∑𝑟(∬𝑆 𝐹𝑖𝑗𝛬𝑢𝑠𝑟𝑑𝑥𝑑𝑦 + ∑𝑒 = 1∬𝑆𝑒𝐹𝑒𝑖𝑗𝛬𝑒𝑠𝑟𝑑𝑥𝑑𝑦) 𝑢

Giving a weight 𝑆𝑢/𝑆 to the first term and 𝑆𝑒/𝑆 to each damaged region of the 𝑢

second term and then integrating over total area and denoting 𝛬𝑢𝑠𝑟, 𝛬𝑒𝑠𝑟(𝛬𝑒𝑠𝑟 = 𝛬𝑠𝑟) by 𝛬𝑠𝑟 gives 𝑀

𝑀

∬𝑆 𝐹𝑖𝑗𝛬𝑢𝑠𝑟𝑑𝑥𝑑𝑦 + ∑𝑒 = 1∬𝑆𝑒𝐹𝑒𝑖𝑗𝛬𝑒𝑠𝑟𝑑𝑥𝑑𝑦 = ∬𝑆[𝐹𝑖𝑗𝑆𝑢/𝑆 + (∑𝑒 = 1𝐹𝑒𝑖𝑗𝑆𝑒)/𝑆]𝛬𝑠𝑟𝑑𝑥𝑑𝑦 𝑢

(15)

Now, introduce the quantities 𝐴𝑝𝑞,𝐷𝑖𝑗 that are called the global extensional and bending stiffnesses, respectively. 𝑀

𝐴𝑝𝑞 = 𝐴𝑝𝑞𝑆𝑢/𝑆 + (∑𝑒 = 1𝐴𝑒𝑝𝑞𝑆𝑒)/𝑆

(16)

𝑀

𝐷𝑖𝑗 = 𝐷𝑖𝑗𝑆𝑢/𝑆 + (∑𝑒 = 1𝐷𝑒𝑖𝑗𝑆𝑒)/𝑆 Then, Eq. (15) is reduced into 𝑀

∬𝑆 𝐹𝑖𝑗𝛬𝑢𝑠𝑟𝑑𝑥𝑑𝑦 + ∑𝑒 = 1∬𝑆𝑒𝐹𝑒𝑖𝑗𝛬𝑒𝑠𝑟𝑑𝑥𝑑𝑦 =∬𝑆𝐹𝑖𝑗𝛬𝑠𝑟𝑑𝑥𝑑𝑦, (𝐹 = 𝐴,𝐷) 𝑢

and also the strain energy (Eq. (10.b)) and total energy (Eq. (10.a)) into

7

(17)

1 3

1

3

𝑈 = 2∑𝑠 = 1∑𝑟(∬𝑆𝐹𝑖𝑗𝛬𝑠𝑟𝑑𝑥𝑑𝑦) = 2∬𝑆(∑𝑠 = 1∑𝑟𝐹𝑖𝑗𝛬𝑠𝑟)𝑑𝑥𝑑𝑦 1

= 2∬𝑆𝛬(𝐴𝑖𝑗,𝐷𝑖𝑗)𝑑𝑥𝑑𝑦 = 𝑈 𝛱 = 𝑈 + 𝑉 = 𝑈 +𝑉 = 𝛱 As seen in Eq. (16) and Eq. (18.a), Eq. (18.b), the global stiffnesses 𝐴𝑝𝑝,𝐷𝑖𝑗 are

(18.a) (18.b)

the equivalent stiffnesses of undamaged plate which correspond to the stiffnesses of plate with DMCs, and 𝛱 is the total energy of equivalent plate. As a result, the expressions have been obtained that make it possible to replace the plate with DMCs by that without DMC. Therefore, using these stiffnesses, the analysis of plate with distributed DMCs can be undertaken through laminated plate theories including FSDT and CLT. Once the stiffnesses in damaged regions are determined through either local (Eq. (4)) or zero stiffness model (Eq. (8)), Eq. (16) gives the global stiffnesses. For example, when using the zero stiffness model, the global stiffness quantities are found as 𝑁

( (𝑄𝑖𝑗)𝑘(1 ―

) )(𝑧3𝑘 ― 𝑧3𝑘 ― 1)

1 𝑀

𝑒 𝐴𝑝𝑞 = ∑𝑘 = 1(𝑄𝑝𝑞)𝑘 1 ― 𝑆∑𝑒 = 1𝛿(𝑒) 𝑘𝛾 𝑆 (𝑧𝑘 ― 𝑧𝑘 ― 1)

𝐷𝑖𝑗 =

1 𝑁 ∑ 3 𝑘=1

1 𝑀 (𝑒) 𝑒 ∑ 𝑆 𝑒 = 1𝛿𝑘𝛾 𝑆

(19)

If the damaged layer numbers 𝛾(𝑒) are the same for all damaged regions, then 𝛿(𝑒) 𝑘𝛾 = 𝛿𝑘𝛾 and, regarding Eq.(7), the above relation Eq. (19) is reduced to the followings. 𝑁

(

𝑆𝑑

)

(20.a)

𝐴𝑝𝑞 = ∑𝑘 = 1(𝑄𝑝𝑞)𝑘 1 ― 𝛿𝑘𝛾 𝑆 (𝑧𝑘 ― 𝑧𝑘 ― 1)

(

1 𝑁

𝑆𝑑

)

𝐷𝑖𝑗 = 3∑𝑘 = 1(𝑄𝑖𝑗)𝑘 1 ― 𝛿𝑘𝛾 𝑆 (𝑧3𝑘 ― 𝑧3𝑘 ― 1)

(20.b)

Similarly, the global coupling stiffnesses 𝐵𝑖𝑗 can be defined. To be exact, because of DMCs, 𝐵𝑖𝑗 are not equal to zero even for symmetric cross ply laminates, but those are not derived here, assuming to be zero for their weak effect on deflection. 2.3. Global stiffness determination based on the CLT In general, the plane stress state is assumed that ignores stress components 𝜎𝑧,𝜏𝑥𝑧,𝜏𝑦𝑧 in laminated plate theory. However, in thin plates with DMCs, these stresses get increased more than those in plates without DMC because of stress concentration around crack tips and thus the availability of CLT should be estimated. As noted in Ref. [23], for a semi-elliptical matrix crack, at the deepest point at which stress concentration is dominant, 𝜎𝑧 is scarcely affected by the crack while 𝜏𝑥𝑧,𝜏𝑦𝑧 are about one tenth times 𝜎𝑥,𝜎𝑦 and two fifth times 𝜏𝑥𝑦 in magnitude. The local stiffness is determined by SERRs (Strain Energy Release Rates) through SIFs (Stress Intensity Factors). For a given crack size, it is known that SIF is linearly related to stress and SERR is proportional to a square root of SIF. Therefore, it can be drawn that the stress components 𝜎𝑧,𝜏𝑥𝑧,𝜏𝑦𝑧 have no remarkable effect on the local stiffness and further global stiffness, compared to 𝜎𝑥,𝜎𝑦,𝜏𝑥𝑦, so that the thin laminated plate theory based on the assumption of plane stress is available even for 8

plates with DMCs. This is, however, not more than an assumption to make consideration easy. Obviously, results based on this assumption may differ somewhat from the exact solutions obtained by taking account of 𝜎𝑧,𝜏𝑥𝑧,𝜏𝑦𝑧, but it is beneficial in analysis procedure, time and cost, etc. It is validated later through the simulations for buckling and initial problem. Now, determine the global stiffnesses using the CLT, in which the stiffnesses 𝐴44,𝐴55 relevant to shear strains are ignored. The energy functionals 𝛱,𝑈 and 𝑉 are given by Eq. (10.a) - Eq. (10.c). Since 𝐴44,𝐴55 are ignored, Eq. (11) and Eq .(12) are simplified: 2

2

𝛬𝑢 = ∑𝑠 = 1𝛬𝑢𝑠(𝐴𝑖𝑗,𝐷𝑖𝑗), 𝛬𝑒 = ∑𝑠 = 1𝛬𝑒𝑠(𝐴𝑒𝑖𝑗,𝐷𝑒𝑖𝑗)

(21.a)

𝛬𝑢1 = 𝐴11(𝑢0,𝑥)2 + 2𝐴12𝑢0,𝑥𝑣0,𝑦 + 𝐴22(𝑣0,𝑦)2 + 𝐴66(𝑢0,𝑦 + 𝑣0,𝑥)2 𝛬𝑢2 = 𝐷11(𝜓𝑥,𝑥)2 + 2𝐷12𝜓𝑥,𝑥𝜓𝑦,𝑦 + 𝐷22(𝜓𝑦,𝑦)2 + 𝐷66(𝜓𝑥,𝑦 + 𝜓𝑦,𝑥)2 𝛬𝑒1 = 𝐴𝑒11(𝑢0,𝑥)2 + 2𝐴𝑒12𝑢0,𝑥𝑣0,𝑦 + 𝐴𝑒22(𝑣0,𝑦)2 + 𝐴𝑒66(𝑢0,𝑦 + 𝑣0,𝑥)2 𝛬𝑒2 = 𝐷𝑒11(𝜓𝑥,𝑥)2 + 2𝐷𝑒12𝜓𝑥,𝑥𝜓𝑦,𝑦 + 𝐷𝑒22(𝜓𝑦,𝑦)2 + 𝐷𝑒66(𝜓𝑥,𝑦 + 𝜓𝑦,𝑥)2

(21.b)

Based on these relations, the global stiffnesses can be defined in the same way as by the FSDT. It is noted that the global stresses obtained are the same as Eq. (16), but both subscripts 𝑝𝑞 and 𝑖𝑗 take the such values as 11, 12, 22, 66. Once stiffnesses in damaged region are determined, the global stiffnesses 𝐴𝑖𝑗, 𝐷𝑖𝑗 are obtained using Eq. (16). As a result, it can be seen that plates with DMCs can be analysed by simply replacing 𝐴𝑖𝑗, 𝐷𝑖𝑗 in solutions for plates without DMC with 𝐴𝑖𝑗, 𝐷𝑖𝑗. For example, a cracked laminated plate under uniformly distributed load can be solved using the Navier’s solution ([26]). Similarly, other problems of plate with DMCs can be solved as well. 3. Results and Discussions 3.1. Bending of the plate with DMCs based on the FSDT Using the global stiffness quantities determined above, analyse a simply supported square laminated plate subjected to a combination of biaxial in-plane load and transverse pressure and compare the results with those by ANSYS simulation. The differential equation for the symmetric cross ply laminated plate is obtained from a stationary condition of the energy functional 𝛱 (Eq. (18.b)). As in the plate without DMC, the equations of equilibrium become uncoupled into two equations for in-plane displacements and three equations for the deflection and rotations. The latter equation concerning the deflection is given by ―𝑙𝐴55𝑤,𝑥 + 𝐷11𝜓𝑥,𝑥𝑥 + 𝐷66𝜓𝑥,𝑦𝑦 ― 𝑙𝐴55𝜓𝑥 + (𝐷12 + 𝐷66)𝜓𝑦,𝑥𝑦 = 0 ― 𝐴44𝑤,𝑦 + (𝐷12 + 𝐷66)𝜓𝑥,𝑥𝑦 + 𝐷66𝜓𝑦,𝑥𝑥 + 𝐷22𝜓𝑦,𝑦𝑦 ― 𝑙𝐴44𝜓𝑦 = 0 𝑙(𝐴55𝑤,𝑥𝑥 + 𝐴44𝑤,𝑦𝑦) + 𝑃𝑥𝑤,𝑥𝑥 + 𝑃𝑦𝑤,𝑦𝑦 + 𝑙𝐴55𝜓𝑥,𝑥 + 𝑙𝐴44𝜓𝑦,𝑦 = ―𝑞

(22)

Representing the unknowns and transverse pressure in double Fourier series as below 9

𝑤(𝑥,𝑦) = ∑𝑚∑𝑛𝑊𝑚𝑛𝑠𝑖𝑛𝛼𝑚𝑥𝑠𝑖𝑛𝛽𝑛𝑦 , (𝛼𝑚 =

𝑚𝜋 𝑎 ,

𝛽𝑛 =

𝑛𝜋 𝑏)

𝜓𝑥(𝑥,𝑦) = ∑𝑚∑𝑛Ф𝑚𝑛𝑐𝑜𝑠𝛼𝑚𝑥𝑠𝑖𝑛𝛽𝑛𝑦, 𝜓𝑦(𝑥,𝑦) = ∑𝑚∑𝑛𝛹𝑚𝑛𝑠𝑖𝑛𝛼𝑚𝑥𝑐𝑜𝑠𝛽𝑛𝑦

(23.a) (23.b) (23.c)

𝑞 = ∑𝑚∑𝑛𝑞𝑚𝑛𝑠𝑖𝑛𝛼𝑚𝑥𝑠𝑖𝑛𝛽𝑛𝑦 produces a system of three algebraic equations with respect to the amplitudes of the corresponding harmonics for each couple (m, n), according to the FSDT. [𝐿𝑖𝑗]{𝑊𝑚𝑛Ф𝑚𝑛𝛹𝑚𝑛}𝑇 = {0 0 𝑞𝑚𝑛}𝑇

(24)

, where

𝐿11 = 𝑙𝐴55𝛼𝑚, 𝐿12 = 𝐷11𝛼2𝑚 + 𝐷66𝛽2𝑛 + 𝑙𝐴55, 𝐿13 = (𝐷12 + 𝐷66)𝛼𝑚𝛽𝑛 𝐿21 = 𝑙𝐴44𝛽𝑛 , 𝐿22 = 𝐿13, 𝐿23 = 𝐷66𝛼2𝑚 + 𝐷22𝛽2𝑛 + 𝑙𝐴44 𝐿31 = (𝑙𝐴55 + 𝑃𝑥)𝛼2𝑚 + (𝑙𝐴44 + 𝑃𝑦)𝛽2𝑛, 𝐿32 = 𝐿11, 𝐿33 = 𝐿21

(25)

It is evident that the solutions 𝑊𝑚𝑛, Ф𝑚𝑛, 𝛹𝑚𝑛 of Eq. (24) are those obtained by the global stiffnesses 𝐴𝑝𝑝,𝐷𝑖𝑗. As an illustrative example, evaluate the maximum deflection and the maximum equivalent (Von Mises) stress of a plate with the following properties. Length by width: 𝑎 = 0.3𝑚, 𝑏 = 0.2𝑚 Plate thickness: ℎ = 4𝑚𝑚 (8 layers) with each layer of 𝑡𝑘 = 0.5𝑚𝑚 Layup: [0°/90°/0°/90°]s Lay material (E-glass/Epoxy): 𝐸11 = 38.6𝐺𝑃𝑎, 𝐸22 = 8.27𝐺𝑃𝑎, 𝜈12 = 𝜈13 = 0.26, 𝜈23 = 0.35, 𝐺12 = 𝐺13 = 4.14𝐺𝑃𝑎, 𝐺23 = 3.06𝐺𝑃𝑎 Delamination sides: 𝑎𝑘 = 30𝑚𝑚, 𝑏𝑘 = 20𝑚𝑚 Length by depth of the matrix crack: 𝑙𝑘 = 20𝑚𝑚, 𝑑𝑘 = 0.5𝑚𝑚 Biaxial in-plane loads: 𝑃𝑥 = 1200𝑁/𝑚, 𝑃𝑦 = 600𝑁/𝑚, Transverse pressure: 𝑞 = 20𝑘𝑃𝑎

Fig. 2. The distribution pattern of DMC couples For simplicity, suppose that DMCs are distributed in a matrix form and number each DMC as shown in Fig. 2. Each DMC is centred in the corresponding element region and their total number ( M ) has its own pattern as 𝑀=1-⑤, 𝑀 =3-④⑤⑥, 𝑀=5-①③⑤⑦⑨, 𝑀=7-①③④⑤⑥⑦⑨. Thus, according to this pattern, it is centred with the number ⑤ as the plate has a single 10

DMC (M=1) and those are distributed in the row of ④⑤⑥ as it has three DMCs(M=3) and so on. This pattern refers to a uniform distribution of DMCs over the plate. In general, the layer number 𝛾(𝑒) having DMC varies with the damaged region (𝑒). Here, without loss of generality and for convenience, assume that the positions of DMC are the same for all damaged regions (𝛾(𝑒) = 𝛾). Hence, taking the top as the first layer and the bottom as the last layer, the delamination is positioned between the 7th and 8th layer while the matrix crack inside the 7th layer (𝛾 = 7), likewise in all damaged regions. The global stiffnesses are determined by Eq. (16) based on Eq. (4) and Eq. (5), while the plate deflection is obtained from Eq. (23.a) using the solution 𝑊𝑚𝑛 of Eq. (24). The first five terms in series are taken in the calculation. For ANSYS simulation, the plate is meshed with 30,350 elements, using an 8-node 3-D Layered Structural Solid Shell Element (SOLSH190). Nodes at interfaces are merged for the undamaged region but not for the DMC layer in damaged regions to create DMCs (Fig. 3).

Fig. 3. Representative DMC model - merged and unmerged nodes Table 1, 2 show the analytical solution by the global stiffness and the ANSYS simulation result by the DMC model with corresponding relative errors, where “GS” refers to the former while “FEM” to the latter. Table 1 is about the non-dimensional deflection (𝑊 = 𝑊 ℎ) at the plate centre. As shown, the deflection increases with 𝑀. The increasing rate is so monotonic, being 0.44 % (𝑊 = 0.7200 and 0.7232) for 𝑀 = 1 and 𝑀 = 3 and 0.45% (𝑊 = 0.7264 and 0.7297) for 𝑀 = 5 and 𝑀 = 7, respectively. For 𝑀 = 1, the relative error between GS and DMC is 3.1% and it is bounded, decreasing as 𝑀 increases in number. Table 1. The deflection (𝑊) at the plate centre with 𝑀 𝑀

0

GS

0.7184 0.7200 0.7232 0.7264 0.7297 0.7329

FEM

0.7431 0.7432 0.7433 0.7434 0.7435 0.7436

Error (%) 3.3

1

3.1

3

2.7

5

2.3

7

1.8

9

1.4

Table 2 shows the equivalent stress values at the bottom centre (𝑧 = 2𝑚𝑚) of plate with 𝑀. Similar to deflection, the analytical value is close to simulation result with the relative error smaller than 2% for each value of 𝑀. Table 2. The equivalent stress 𝜎𝑒 (𝑀𝑃𝑎) at the bottom centre with 𝑀 11

𝑀

0

1

3

5

7

9

GS

21.936

21.987

22.091

22.196

22.302

22.408

FEM

22.357

22.422

22.450

22.438

22.396

22.351

Error (%)

1.88

1.94

1.59

1.08

0.42

0.26

Next, evaluate the effect of DMC distribution pattern Ω on deflection and stress state using two variants, i.e. Ω1-①④⑦ and Ω2-④⑤⑥. For Ω1, three DMCs are left-sided from the plate centre. Thus, it can be seen that Ω1 is unsymmetrical, but Ω2 is symmetric with respect to the plate centre. Table 3 lists the deflection and equivalent stress comparatively for both patterns. The position of DMC in damaged regions remains the same as above (𝛾 = 7). As obvious in the expressions for global stiffness, the analytical deflections are the same for both patterns. But, the simulation results are different each other, with the relative errors to GS equal to 2.66% for Ω1 and 2.70% for Ω2, respectively. Table 3. The deflection and equivalent stress for different patterns 𝜎𝑒(GS) 𝜎𝑒(FEM) Ω 𝑊(GS) 𝑊(FEM) Ω1 0.7232 0.7430 22.091 22.349 Ω2 0.7232 0.7433 22.091 22.450 Concerning equivalent stress values at the bottom centre (𝑧 = 2𝑚𝑚) of plate for both patterns, the relative errors between GS and FEM solutions are 1.15% (Ω1) and 1.59 %(Ω2). It is worth noting that results are closer each other for unsymmetrical rather than symmetric patterns and thus the present method gives more realistic results for unsymmetrical patterns, although the difference is very small. The local stiffness model reflects different status of DMCs from one the zero stiffness model does, so avoiding to compare results by two models. Instead, buckling and initial deflection problem of plate are analysed below by using the zero stiffness model.

3.2. Buckling and initial imperfection of the plate with DMCs based on the CLT Firstly, consider buckling behaviour of a simply supported cross ply laminated plate with DMCs. According to the global stiffness method, the buckling load can be found by replacing stiffness quantities in the known buckling load formula for its counterpart without DMC by the global stiffnesses. Assuming the buckling mode as Eq. (23.a), the inplane internal forces are 𝑁𝑥 = ― λ𝑃𝑥, 𝑁𝑦 = ― λ𝑃𝑦, 𝑁𝑥𝑦 = 0 for the plate subjected to biaxial edge loads 𝑃𝑥,𝑃𝑦 and the buckling load factor 𝜆 is given as ([25]) 𝜆𝑐𝑟 =

𝐷11𝛼4𝑚 + 2(𝐷12 + 2𝐷66)𝛼2𝑚𝛽2𝑛 + 𝐷22𝛽4𝑛 𝑃𝑥𝛼2𝑚 + 𝑃𝑦𝛽2𝑛

12

(26)

, where 𝐷𝑖𝑗 are the global stiffnesses defined by Eq. (20.b). When the plate is subjected to uniaxial compression (𝑃𝑥 ≠ 0,𝑃𝑦 = 0), buckling occurs with a half wave perpendicular to compression and thus the lowest critical load (𝑛 = 1) is 𝜋 2

2

𝑎 4

(𝑎𝑏)

𝑃𝑥,𝑐𝑟 = (𝑚𝑎) [𝐷11𝑚4 +2(𝐷12 + 2𝐷66)𝑚2

+ 𝐷22(𝑏) ]

(27)

The value of 𝑚 corresponding to the buckling load, is specified dependent on the known 𝐷

1/4

11 inequalities between aspect ratio 𝑎 𝑏 and stiffness ratio 𝜇 = (𝐷22)

𝑎

𝑎

𝑚 = 1;𝑏 ≤ 2𝜇, 𝑚 = 2; 2𝜇 ≤ 𝑏 ≤ 6𝜇, 𝑚 = 3;

. 𝑎

6𝜇 ≤ 𝑏 ≤ 12𝜇,⋯

(28)

Therefore, it can be seen that the value of stiffness ratio (𝜇) for plates with DMCs is dependent not only on their layup and material properties as in its perfect counterpart, but the position (𝛾) and the number (𝑀) of DMCs. As an example, find the buckling load of the square plate (𝑎 = 𝑏 = 0.2𝑚). The other geometric and material properties including material, layup, DMC size, etc. are the same as in 3.1. However, the finite element model is different as the zero stiffness model is used. In order to model a completely damaged layer that can not resist any longer, that particular layer in the damaged region is meshed with elements having so small isotropic materal properties as if it were a weak and soft core of sandwich plate. In other words, it is replaced by a virtual “local core”. In simulation 𝐸 = 1000𝑁 𝑚2,𝜈 = 0 are used as material properties of the local core. As well, the plate is meshed with 1,490 elements, using a 4-Node Structural Shell (SHELL181). The global stiffnesses are determined by Eq. (20.b) and also the buckling load as well as the corresponding half-wave number by Eq. (27) and Eq. (28), using the given aspect ratio 𝑎 𝑏 = 1. Table 4 lists the analytical buckling loads by the global stiffness model (zero stiffness) and the ANSYS simulation results by the DMC model (virtual local core) for the case 𝛾 = 7, respectively. It is noted that the simulation results are the lowest global buckling loads which correspond to one halfwave mode shape and appear after some local ones. Table 4. The buckling load 𝑃𝑥,𝑐𝑟(10𝑘𝑁/𝑚) with 𝑀 𝑀

0

1

3

5

7

9

GS

9.0119 8.9985 8.9718 8.9450 8.9183 8.8915

FEM

8.9895 8.8072 8.8879 8.8552 8.9974 8.9412

Error (%) 0.25

2.13

0.94

1.00

0.89

0.56

As seen, the relative error between GS and FEM is fluctuating but fairly small and bounded as 𝑀 increases. Fig. 4 shows the dimensionless buckling load 𝑃𝑥,𝑐𝑟(𝑃𝑥,𝑐𝑟 = 𝑃𝑥,𝑐𝑟/𝑃0𝑥,𝑐𝑟) with 𝑀 for three different DMCs (𝛾 = 5, 6, 7), where 𝑃0𝑥,𝑐𝑟 is the analytical buckling load for 𝑀 = 0, i.e. the 13

perfect plate with no DMC, being equal to 90.119𝑘𝑁 𝑚. As shown in the figure, the buckling load decreases linearly for each DMC as 𝑀 increases. The bigger 𝛾 is, the more rapidly the buckling load decreases with 𝑀, which means that it is much more affected by outer DMC than inner one. In detail, comparing the buckling loads for 𝑀 = 3, 5, it can be seen that the decreasing rate is equal to 0.02 % (𝑃𝑥,𝑐𝑟 = 0.9996 and 0.9994) as 𝛾 = 5 and 0.16 % (𝑃𝑥,𝑐𝑟 = 0.9975 and 0.9959) as 𝛾 = 6, while 0.44 % (𝑃𝑥,𝑐𝑟 = 0.9933 and 0.9889) as 𝛾 = 7. Fig.5 shows the variation of dimensionless buckling load 𝑃𝑥,𝑐𝑟 as a function of the aspect ratio 𝑎/𝑏 for 𝑀 = 5, 𝛾 = 7. It is noted that the DMC size, i.e. 𝑎𝑘, 𝑏𝑘 are also changed by the same ratio as 𝑎/𝑏.

Fig.4. 𝑃𝑥,𝑐𝑟-𝑀 relation

Fig.5. 𝑃𝑥,𝑐𝑟-𝑎/𝑏 relation

At first, the buckling load decreases with 𝑎/𝑏, by 47.4 % (𝑃𝑥,𝑐𝑟 = 3.04 and 1.60) at 𝑎/𝑏 = 0.4 and 𝑎/𝑏 = 0.6 , 28.75% (𝑃𝑥,𝑐𝑟 = 1.14) at 𝑎/𝑏 = 0.8, 13.15% (𝑃𝑥,𝑐𝑟 = 0.99) at 𝑎/𝑏 = 1 and then increases after the turning point (𝑎 𝑏 ≈ 1.2). Thus, it can be seen that, as 𝑎/𝑏 ≤ 1.2, the buckling load decreases rapidly for the early stage and then gets monotonous. After that, it gradually increases as 1.2 ≤ 𝑎/𝑏 ≤ 1.6 and decreases again with different mode shape (𝑚 = 2) as 𝑎/𝑏 ≥ 1.6. This buckling behaviour is very similar to that of plates without DMC and also verified by the simulation. Finally, consider a simply supported cross ply laminated plate with an initial deflection 𝑤0, subjected to biaxial compression 𝑃𝑥,𝑃𝑦. As known, initial deflections neither introduce a coupling between inplane displacements and transverse deflections, nor affect boundary conditions. Accordingly, the solution can be sought representing both initial deflections (𝑤0) and the deflections (𝑤) from the initially deflected surface in double Fourier series satisfying boundary conditions. 𝑤0(𝑥,𝑦) = ∑𝑚∑𝑛𝑊0𝑚𝑛𝑠𝑖𝑛𝛼𝑚𝑥𝑠𝑖𝑛𝛽𝑛𝑦 , 𝑤(𝑥,𝑦) = ∑𝑚∑𝑛𝑊𝑚𝑛𝑠𝑖𝑛𝛼𝑚𝑥𝑠𝑖𝑛𝛽𝑛𝑦

(29)

Substituting these series into the linear equation of equilibrium and applying the global stiffness method yields the amplitude of the mnth harmonics: ([25])

𝑤𝑚𝑛 = 𝐷

(𝑃𝑥𝛼2𝑚 + 𝑃𝑦𝛽2𝑛)𝑊0𝑚𝑛 4 11𝛼𝑚

+ 2(𝐷12 + 2𝐷66)𝛼2𝑚𝛽2𝑛 + 𝐷22𝛽4𝑛 ― 𝑃𝑥𝛼2𝑚 ― 𝑃𝑦𝛽2𝑛 14

(30)

, where 𝐷𝑖𝑗 are the global stiffnesses and 𝛼𝑚, 𝛽𝑛 are given as Eq.(23.a). Regarding Eq. (27) and Eq. (30) in case of the uniaxial compression (𝑃𝑥 ≠ 0,𝑃𝑦 = 0) gives the total deflection (𝑊 = 𝑤0 +𝑤) as 𝑤0𝑚𝑛

𝑊𝑚𝑛 = 1 ― 𝑃

(31) 𝑥

, where 𝑃𝑥 = 𝑃𝑥 𝑃0𝑥,𝑐𝑟, 𝑤0𝑚𝑛 = 𝑤0𝑚𝑛 ℎ,𝑊𝑚𝑛 = 𝑊𝑚𝑛 ℎ are the non-dimensional inplane load, initial and total deflection, respectively. Table 5 shows the load-deflection relation 𝑃𝑥 ― 𝑊(𝑊/ℎ) of the uniaxially compressed plate. The plate properties including plate size, material, layup, DMC size, etc. are the same as taken in 3.1. The position of DMC is 𝛾 = 6 and the initial deflection at the plate centre is 𝑤0 = 0.125 (𝑚 = 𝑛 = 1). The global stiffnesses are given in Eq. (20.b) by the zero stiffness model. As shown in the table, the deflection increases together with 𝑀 for the given value of compressive load. In detail, the deflection (𝑊) increases by 0.47% equally for 𝑃𝑥 = 0.5 and by 4.3%, 4.6% for 𝑃𝑥 = 0.9. Table 5. The deflection (𝑊) at the plate centre as a function of 𝑃𝑥 0.1

0.3

0.5

0.7

0.9

3

0.3190

0.1791

0.2517

0.4234

1.3312

5

0.1390

0.1794

0.2529

0.4280

1.3922

7

0.1391

0.1798

0.2541

0.4328

1.4597

𝑀

On the other hand, the increasing rate of deflection with 𝑃𝑥 is getting high. For 𝑀 = 3, the rate is 22.4% at first between 𝑃𝑥 = 0.1 and 𝑃𝑥 = 0.3 while rises upto 68.2% at the end between 𝑃𝑥 = 0.7 and 𝑃𝑥 = 0.9. Accordingly, the deflection to be primarily affected by 𝑃𝑥, although it increases with 𝑀. As is obvious from Eq. (31), the deflection increases unlimitedly as the compressive load approaches the critical value (𝑃𝑥→1). A similar behaviour can be found through Table 6 in which the effect of initial deflection on total deflection is evaluated. Table 6. The deflection (𝑊) at the plate centre as a function of 𝑤0 0.1

0.2

0.3

0.4

0.5

3

0.2014

0.4027

0.6041

0.8055

1.0068

5

0.2023

0.4046

0.6069

0.8092

1.0115

7

0.2032

0.4065

0.6097

0.8130

1.0162

𝑀

Hence, it is evident that the deflection increases dependent not only of the initial deflection, but the number as well as the position of DMCs. The behaviour of initially deflected plate with DMCs is similar to that of its perfect counterpart. A plate begins to bend 15

as soon as the load is applied. The deflection increases slowly at early stage but rapidly as the load approaches the critical load. The more the initial deflection and the DMC number, the more deflection increases. However, as the load approaches the critical load, it gets unbounded regardless of the initial deflection and the DMC number, meaning that DMCs give the same affection to the plate as initial imperfections. 4. Conclusions The global stiffness method makes it possible to analyse laminated plates with DMCs by simply replacing the stiffnesses in known formulas for their undamaged counterparts. The global stiffnesses of laminated plates with distributed DMCs are analytically determined and then applied to bending and buckling problem. The comparison shows fair agreements between solutions by the global stiffness and the DMC model (ANSYS). Especially, the zero stiffness model is useful to estimate the ultimate effect by DMCs since it assumes layers to be fully cracked and collapsed. Furthermore, the effect of DMC distribution patterns on deflection, stress and buckling behaviour is investigated, which shows that the present method is simple and easy to apply, representing local damage status of plates reasonably. The present method is available not only for symmetric cross ply laminates, but for any other layups. Acknowledgements The authors are grateful to the anonymous reviewers for their helpful comments which have improved the manuscript significantly.

References [1] Allam MNM, Zenkour AM. Stress concentration factor of a structurally anisotropic composite plate weakened by an oval opening. Composite Structures 2003; 61: 199–211. [2] Wang J, Tong L. A study of the vibration of delaminated beams using a nonlinear antiinterpenetration constraint model. Composite Structures 2002; 57: 483–88. [3] Parlapalli MSR, Wenge T, Shu D. Buckling analysis of tri-layer beams with enveloped delaminations. Composites: Part B 2005; 36: 33–9. [4] Ahmed A, Sluys LJ. A three-dimensional progressive failure model for laminated composite plates subjected to transverse loading. Engineering Fracture Mechanics 2013; 114: 69–91. [5] Kaya C, Kaya F, Butler EG, Boccaccini AR, Chawla KK. Development and characterisation of high-density oxide fiber-reinforced oxide ceramic matrix composites with improved mechanical properties. Journal of the European Ceramic Society 2009; 29:1631-9. [6] Zhang Y, Zhu P, Lai X. Finite element analysis of low-velocity impact damage in composite laminated plates. Materials and Design 2006; 27:513-9. [7] Yang W, Pan Y, Pelegri AA. Multiscale modeling of matrix cracking coupled with interfacial debonding in random glass fiber composites based on volume elements. Journal of Composite Materials 2013; 47(27): 3389-99. 16

[8] Lee YJ, Lee CH, Fu WS. Study on the compressive strength of laminated composite with through the width delamination. Composite Structures1998; 41:229-41. [9] Yan YJ, Yam LH. Online detection of crack damage in composite plates using embedded piezoelectric actuators/sensors and wavelet analysis. Composite Structures 2002; 58: 29–38. [10] Talreja R. Damage mechanics of composite materials. Composite materials series 9. Amsterdam: Elsevier; 1994. [11] Nakamura T, Wu LC. Effects of ply-arrangement on compressive failure of layred structures. Egineering Fracture Mechanics 2000; 67: 421-443 [12] Cheng ZQ, He LH, Kitipornchai S. Influence of imperfect interfaces on bending and vibration of laminated composite shells. Internatinal Journal Solids and Structuers 2000; 37:2127-2150 [13] Czarnocki P. Circular delaminations-the effect of thermal strains on the strain energy release rate. Mechanics of Composite Materials 2001; 37(1): 67-70 [14] Wang JS, Lin CC. Engineering analysis of delaminated beam plates. Composite Structures 1996; 34: 397-407. [15] Whitcomb JD. Analysis of a laminate with a postbuckling embedded delamination, including contact effects. J Compos Mater 1992; 26(10):1523–35. [16] Kim HN, Keith TK. A method for modeling the local and global buckling of delaminated composite plates. Composite Structures 1999; 44: 43-53. [17] Chattopadhyay A, Haozhong G. Modelling of delamination in composite cylindrical shells with a new higher-order theory. Compos Sci Technol 1995; 54(2):223–35. [18] Haozhong G, Chattopadhyay A. Delamination buckling and postbuckling of composite cylindrical shells. AIAA J 1996; 34(6):1279–86. [19] Zhu JF, Zhang G, Howson WP, Williams FW. Reference surface element modelling of composite plate/shell delamination buckling and postbuckling. Composite Structures 2003;61: 255–64. [20] Elisa Pietropaoli, Aniello Riccio. A Global/Local Finite Element Approach for Predicting Interlaminar and Intralaminar Damage Evolution in Composite Stiffened Panels Under Compressive Load. Appl Compos Mater (2011) 18:113–125 [21] Endel V. Iarve, Mark R. Gurvich, David H. Mollenhauer, Cheryl A. Rose and Carlos G. Dávila. Mesh-independent matrix cracking and delamination modeling in laminated composites. Int. J. Numer. Meth. Engng 2011; 88:749–773 [22] D. H. Li, Y. Liu and X. Zhang. An extended Layerwise method for composite laminated beams with multiple delaminations and matrix cracks. Int. J. Numer. Meth. Engng 2015; 101:407–434 [23] In-Bong Kim, Song-Hak U, Tong-Hwan Ji. Determination of local stiffness of a crossply composite plate with delaminations and matrix cracks. Composite Structures 2018; 196: 127–134. [24] L. Zubillaga, A. Turon a, J. Renart, J. Costa, P. Linde. An experimental study on matrix crack induced delamination in composite laminates. Composite Structures 2015; 127: 10–17 17

[25] Birman V. Plate Structures. Springer, ISSN 0925-0042, 2011. [26] Eduard Ventsel, Theodor Krauthammer. Thin plates and shells- theory, analysis, and applications, 2001.

18

Conflict of Interest No.

19