Multiscale studies on the nonlinear vibration of delaminated composite laminates–global vibration mode with micro buckles on the interfaces

This paper presents a multiscale approach to study the nonlinear vibration of fiber reinforced composite laminates containing an embedded, through-width delamination dividing the laminate into four sub-laminates. The equations of motion are established from macroscopic nonlinear mechanics for plates and shells and micro-mechanics of composite material to allow for the influences of large amplitude, membrane stretching in the neutral plane, and the interactions of the sublaminates. Analytical solutions obtained in this paper reveal that the interaction penalty at the interfaces plays a coupling effect between sublaminates, which eventually alters the vibration characters of the four-sublaminate lamina in macroscopic and microscopic mechanism. From a macro perspective, sub-laminates above and below the delamination vibrate in exactly the same mode in spite of their different stiffness and the four-sublaminate lamina has a consistent global vibration mode. In accompanying with the macro vibration, micro buckles occur on the interfaces of the delamination with amplitude about 10−3 times of that of the global mode. It is found that the vibration frequency is an eigenvalue of the delaminated lamina determined only by the geometry of the delamination. Authentication of the multiscale study is fulfilled by comparing the analytical solutions with the FEA results.

Carbon-fiber-reinforced composite materials consist of two parts: matrices and reinforcements. The reinforcements are usually carbon fibers which provide the strength and rigidity for the composite materials, measured by the stress and the elastic modulus respectively. The matrix is usually a polymer resin, such as epoxy, to bind the reinforcements together. Unlike isotropic materials, such as steel and aluminum, carbon-fiber-reinforced composite has directional strength and properties. The properties of carbon-fiber-reinforced composite depend on the layouts of the carbon fiber and the proportion of the carbon fibers in the polymer.
During manufacturing process, fiber breakage, matrix cracking, and fiber/matrix interface debonding may exist in the composite materials and, as a result, develop delamination in the composite materials. Delamination is one of the primary failure types in composite laminates and will cause the stiffness degradation and strength loss of the entire laminate, consequentially the premature failure of the composite materials. The failure of fiber reinforced composite structures containing buried delamination is much complex due to presence of various mechanisms. The analysis of the response of the delaminated composite structures under the external loadings needs to account for various factors, e.g., large deformation, nonlinear constitutive behavior, failure, delamination, and interaction penalty effects, etc. Ignoring the interaction penalty effects may result in the underestimation of the load carrying capacity and the penetration of the deflection between the sublaminates above and below the interfaces of delamination 1 . Therefore, researchers have launched various studies to investigate the effects of contact on the mechanical responses of the delaminated composite structures.
Numerous approaches as well as different results have been reported in many literatures. Chai et al. 2 may be one of the first to carry out the study of the delaminated composite laminates. By simplifying the composite plate as a beam-plate, they proposed a one dimensional analytical model to assess the compressive strength of a delaminated composites plate. Cranford 3 investigated the mechanical stability of graphene-based nanoscale multi-layers and derived the critical length scale and required adhesion energy to prevent delamination under buckling conditions. Oyewole et al. 4 presented the results of experimental and theoretical/computational micro-wrinkles and buckling on the surfaces of stretchable poly-dimethylsiloxane (PDMS) coated with nano-scale Gold (Au) layers. The interfacial delamination occurred along with the buckling was also studied using finite element simulations of the interfacial crack growth. Ovesy et al. 5 suggested a novel layerwise theory to evaluate the buckling and post-buckling behavior of delaminated composite plates with multiple through-the-width delaminations based on the first order shear deformation theory (FSDT). Giannakopoulos et al. 6 and Whitcomb 7 investigated the stability by considering the nonlinearity associated with the contact problem. The results reveal that the effects of the contact conditions are significant, especially with regard to the propagation of the debonded region. Hu et al. 8 introduced artificial springs to compute the fictitious contact forces occurred in the buckling of laminates with an embedded delamination based on the Mindlin plate theory. The other research topics include simulation of the multi-failure responses, see Ho et al. 9 ; multiscale method to compute the elastic deformations for general bilayer, see Kumar et al. 10 ; contact of interfacial debonding front region and face matrix cracking, etc. in face/core interface debonded composite sandwich plates, see Chen and Bai 11 ; the non-linear behaviors of delaminated sandwich panels with a compressible core, see Frostig and Thomsen 12 ; fatigue delamination growth for piezoelectric laminated cylindrical shells, see Zhu et al. 13 ; analysis of the mixed mode delamination in composite laminates from point view of fracture mechanics and contact mechanics, see Bruno et al. 14 ; experimental studies on laminated plates with strip-type delamination under pure bending, see Yet et al. 15,16 .
Except the static analysis, the dynamic responses of composite laminate with delamination are also the focus of a number of studies. One of the earliest models for vibration analysis of composite beams with delaminations was proposed by Ramkumar et al. 17 . They modeled a beam with one through-width delamination by simply using four Timoshenko beams connected at delamination edges. Wang et al. 18 improved the analytical solution by including the coupling between flexural and axial vibrations of the delaminated beams which were modeled as split regions. Recurrence equations relating integration constants for adjacent split regions were established by satisfying continuity conditions at junctions of sub-regions. Ramkumar et al. 's and Wang et al. 's models were essentially one-dimensional models. Dynamic responses of laminates with two-dimensional embedded delamination were studied by Dey and Karmakar 19 to examine free vibration of multiple delaminated angle-ply composite conical shells using Mindlin's theory and the multi-point constraint algorithm and by Noh and Lee 20 to analyze the dynamic stability of laminated skew plates by developing a finite element formulation based on HSDT. Such models are called "free model". As pointed out in present published paper 21 , the "free model" result in interpenetration of the delaminated sub-regions, which is physical impossible. To overcome this defect a piecewise linear spring model inserted in the delaminated region was introduced. Several researchers conducted their studies based on this model. Luo and Hanagud 22 carried out the dynamic response analysis of a delaminated beam by taking into account of shear effect, rotary inertia terms and bending-extension coupling. Kargarnovin et al. 23 examined a delaminated Timoshenko beam under the motion of a constant amplitude point force traveling with uniform velocity by accounting for the Poisson's effect, shear deformation and rotary inertia. Chen et al. 24,25 presented a formula of element stiffness and mass matrices for the composite laminates using the first-order shear deformation theory combined with the selecting numerical integration scheme. A virtual interface linear spring element was employed to avoid the overlap and penetration phenomenon between the upper and lower sublaminates at the delamination region. Chattopadhyay et al. 26 presented an FE model using the first-order zigzag theory for delaminated composites and smart composite plates. The linear spring model needs to specify the spring constant before solving the problem, which is basically equivalent to specifying whether the sublaminates are in contact or not beforehand, this model therefore does not reflect the essential feature of the contact problem, that is, the contact region and the degree of interpenetration are all unknown a priori. Oh et al. 27 and Kwon and Aygunes 28 developed finite element models adopting the approach given by Hughes et al. 29 for correcting the velocity, acceleration and contact force values during release-to-contact condition. Their models were able to achieve nonpenetration of interlaminates but the results they present can still be argued upon for its physical correctness, as both the top and bottom laminates experience a kind of total separation. In reality, they must collide more than once during the vibration as both of the sub-laminates are shown to be vibrating by themselves at higher frequencies than the contact occurrence. Recently, Wang and Tong 30 suggested a nonlinear interpenetration constraint model to handle the contact problem. In this nonlinear model, the contact force was considered to be proportional to the relative deflections of the mid-planes of subregions above and below the delamination interfaces. Schwarts-Givli et al. 31 further put forward a step function in the nonlinear constrain model to incorporate "with and without contact" conditions in the governing equations, which reflects the nonlinear nature of the contact behavior.
Although, a vast literature exists on the effects of delamination on the dynamic response of composite structures, they have mostly ignored large deformation effects. When delaminated laminates undergo large deflection, the large deflection induces geometric nonlinearity in the delaminated region and leads to the couple effect between the membrane stretching and transverse bending and the occurrence of the contact phenomenon at the delaminated region. As a result, non-uniform interaction constrains are produced along the interfaces in-between the upper and lower sublaminates. They restrict the two sublaminates to deform penetrating with each other by causing a local compressive stress and a vertical deformation along the interfaces of the delamination [32][33][34] . The bending analysis of the delaminated plates is an interdiscipline and multiscale problem. The research involves two mechanical mechanisms: theories of plates and interaction mechanics. The multiscale composition of the composite lamina causes it to undergo macro deformation as well as micro buckling. The problem becomes even more intricate due to the obscurity of the interaction properties, such as the locations and sizes of the contact area, the intensity of the contact force, etc. For this reason, almost all the analyses and results in above mentioned literatures are accomplished in virtue of computer techniques. Such analyses may lead to an argument that the contact properties, such as the locations and the sizes of the contact area, are unpredictable.
In this article, we develop a general framework to analyze the nonlinear vibration of fiber-reinforced composite laminate with embedded delamination using a multi-scale modeling approach and infer possible consequences of delamination geometry on vibration properties. While the interaction forces at the interfaces of the delamination is investigated from the point view of micro-mechanics, governing equations for the delaminated laminates undergoing nonlinear vibration are derived based on von Karman nonlinear mechanics for plates and shells to include nonlinear geometry deformation and membrane stretching and are solved via Galerkin approach. To demonstrate the validity of our method, we analyze the nonlinear free vibration of a simply supported, symmetrically cross-plied laminate with a though-width delamination. We obtain close-form analytical solutions and provide the macroscopic and microscopic explanation for the outcome of the consistent macro vibration mode and the micro buckles on the interfaces of the delamination. We have compared the vibration frequency as well as the vibration mode predicted by our multi-scale approach with the dynamics simulations from ABAQUS, obtained an excellent agreement between these two approaches.

Theoretical modelling
Material modeling. This article mainly focuses on investigating the nonlinear vibration of a delaminated composite laminate. The study includes investigation of the overall dynamic responses of the composite laminate and determination of the local transverse deformation of sub-laminates above and below the delamination due to the interaction forces between them. The overall dynamic responses of the composite laminate are mainly achieved by virtue of shell theories in which the in-plane material properties is of much more significance. Figure 1a-f illustrate the material modeling of a symmetrically cross-plied composite laminate. The laminate consists of unidirectional fiber reinforced composite plies (Fig. 1a). By performing analysis shown in Fig. 1b-e, the FRCP is homogenized as anisotropic plies with the following equivalent material properties of Young's modulus E 11 and E 22 , Poisson's ratio μ 21 and Shear modulus G 12  12 where E m , μ m , G m and v m are the Young's modulus, Poisson's ratio, Shear modulus and volume fraction of the matrix, and E f , μ f , G f and v f are the corresponding properties of the fiber. The laminate is manufactured by stacking the fiber reinforced composite plies in an orientation sequence, as shown in Fig. 1f. Let θ be the angle of the orientation of m-th ply with respect to the horizontal axis in the global coordinate systems, the stiffness coefficient matrix of the m-th ply Q ( ) jk m in the global coordinate system is given by 35  where c = cosθ, s = sinθ, and To investigate the influence of the delamination on the dynamic behaviors of the delaminated composite laminates, we consider a rectangular composite lamina of length a, width b, and height h. The composite lamina consists of K layers of fiber-reinforced composite ply with a through-width delamination located at depth h 2 from its top surface. The delamination divides the spatial region of the laminate into four sub-laminates of length l i and height h i , separately referred to as sub-laminate Ω (i) (i = 1, 2, 3, and 4), as shown in Fig. 1g. The mechanical properties of the membranous stiffness A (i) , the tension-bending coupling stiffness B (i) and the bending stiffness D (i) of sub-laminate Ω (i) (i = 1, 2, 3, and 4) are given by the following formula 35 In this article, we are interested in analytical investigation on the interaction mechanism between the sub-laminates and its influence on the nonlinear vibration properties. Such analysis requires solving eight nonlinear partial differential equations and thirty-two auxiliary conditions. To clearly present our analytical approach without hindered by tedious calculation, we consider that all the four sub-laminates Ω (i) are symmetrically cross-plied. In this way, the in-plane resultant stresses N (i) (x, y, t) and the bending moment M (i) (x, y, t) of sub-laminate Ω (i) are related to the in-plane strain components ε (i) (x, y, t) and the change of curvatures κ (i) (x, y, t), respectively, through the following constitutive relations: Because of much small thickness dimension comparing to the in-plane dimensions, the transverse displacement at any point in the laminate is represented by the deflection at the corresponding point on the middle plane in the shell theories, i.e., the unit elongation in the transverse direction of the composite laminate ε z is ignored. For this reason, the material property in the transverse direction of the composite laminate E 33 is not required in shell theories. Such approximation won't produce obvious error in analyzing the overall responses of the composite laminate. Nevertheless, to scrutinize the nonlinear vibration of the delaminated composite laminate, our discussion is on the conduct of local interaction effects between the sub-laminates.
Prior to the commencement of local deformation, we need to ascertain the material properties of the four sub-laminates in transverse direction. As shown in Fig. 1f, a composite laminate is manufactured by stacking the composite plies in an orientation sequence in which the angle of orientation of each composite ply is inclined in-plane with respect to x -or y -axis. Therefore, the vertical stacking of the composite plies does not affect the material properties in transverse direction. In other words, the material properties in transverse direction are Scientific RepoRts | 7: 4468 | DOI:10.1038/s41598-017-04570-3 sensitive to the material properties and the volume frictions of the matrix and the fiber but are insusceptible to the orientation sequence. Following this observation, we draw a very important conclusion that regardless of the stacking sequence and the number of plies, the equivalent Young's modulus in transverse direction of the four sub-laminates are identical and equal to that of the single ply provided that the composition for each single ply remains unchanged, i.e.

(4) 33
Interaction kinematics of sub-laminates. When the laminate undergoes vibration, it experiences transverse deflection. Once the laminate deviates from its equilibrium position, it starts to bend first at the center of sub-laminate Ω (2) . Sub-laminate Ω (2) continues to deflect by pushing sub-laminate Ω (3) to deform downwards simultaneously until both Sub-laminates Ω (2) and Ω (3) reach their equilibrium states. During this process, a reaction force q * (x, y, t), called interaction force, is produced between Sub-laminates Ω (2) and Ω (3) . The vibration modes represent the overall response of the delaminated laminate to the external load, thus the amplitudes and the patterns of the external forces, the supporting conditions at the boundary and the geometries of the laminate determine the modes of the vibration. However, the general modes are not the prior to determine the interaction constraint. On the contrary, they depend on the interaction constraint. At any particular moment of the vibration, each of the sub-laminate Ω (i) experiences a particular deflection w (i) (x, y, t) (for i = 1, …, 4). Depending on the amplitudes of w (2) and w (3) , contact regions (where w (2) ≥ w (3) ) and voids (where w (2) < w (3) ) may form alternately along the interfaces of the delamination between Sub-laminates Ω (2) and Ω (3) . Although the contact effects, such as the positions where the contact may occur, the size of the contact areas and the magnitude of the contact force, are still unknown, it does not prevent us from further analyzing. It sounds quite rational to consider that the contact may occur at several regions along the interface of sub-laminates Ω (2) and Ω (3) . Denoting the coordinates of one of such regions as R ∈ [x, x + Δx], a contact force q * (x, y, t) is produced at this place since the deflection of sub-laminates Ω (2) is larger than that of Ω (3) , as shown in Fig. 2a. The following analysis is based on Contact Region R, but we will show later in Subsection Nonlinear dynamic stability analysis that the analysis is valid for the entire laminate. As shown in Fig. 2a, the interaction force q * (x, y, t) developed at the contact regions produces thickness reductions by amounts of δ (2) and δ (3) for Sub-laminates Ω (2) and Ω (3) , respectively (2) 2 33 where h 2 and h 3 are original thickness of Ω (2) and Ω (3) , and E 33 (2) and E 33 (3) are the equivalent Young's modulus of Ω (2) and Ω (3) in the z-direction. The relative movement of the mid-plane of sub-laminates Ω (2) and Ω (3) is approximately Put Eq. (6) into Eqs (7) and (8) and combine with Eq. (9). The kinematic relation between the relative movement w (2) -w (3) of Sub-laminates Ω (2) and Ω (3) and their interaction forces q * is found to be where k is the interaction factor given by: 33 Figure 2. Illustration of the deformation mechanism in one of the contact regions R∈[x, x + Δx] between Sub-laminates Ω (2) and Ω (3) during vibration. (a) Local interaction kinematics at the interfaces of Sub-laminates Ω (2) and Ω (3) . (b) A compatiable deformation between macro deflection and micro local deformation of Sublaminates Ω (2) and Ω (3) .
Scientific RepoRts | 7: 4468 | DOI:10.1038/s41598-017-04570-3 The interaction brace plays a pivotal role in our problem. It determines both the general response (macro behavior) and the local buckling (micro behavior) of the laminate subjected to the action of the external load.
Nonlinear theory of dynamic analysis. We consider the composite laminate under the action of an external excitation. The external excitation is resolved into two components: a horizontal one acting in the middle plane of the composite laminate and a transverse one normal to the middle plane. The horizontal component causes the laminate to undergo in-plane or membrane deformation and the transverse component induces out-of-plane deflection or bending in the composite laminate. Therefore, the deformation mechanisms of the composite laminate are analyzed via both membrane and flexural mechanisms of plates. When the in-plane deformations and the flexural deflection of the laminate are small, the membrane and the flexural shell models are analyzed independently. Under such circumstance, we say that the problem is linear. On the contrary, the membrane and the flexural analyses are coupled with each other. We refer to such case as nonlinear problem. In this section, we derive a theory of nonlinear vibration for delaminated composite laminate.
We begin with the nonlinear membrane kinematics. When the delaminated composite laminate is subjected to high levels of acoustic pressure, it undergoes large amplitude vibrations. The large amplitude produces in-plane membrane stretching in the composite laminate which is considered as in-plane force resultants applied to the middle planes of each of the sub-laminate and parallel to their middle planes. In these conditions, the problems are reduced to plane stress problems. Ignoring the body forces and the in-plane inertial forces, the differential equations of equilibrium for membrane deformation of sub-laminate Ω (i) are obtained from theory of elasticity as follows: where a comma denotes partial differentiation with respect to the corresponding coordinates. Equation (12) are satisfied by introducing a stress function ψ (i) (x, y, t) such that The in-plane strain components ε (i) (x, y, t) are related to the in-plane displacements and the transverse deflection through von Karman nonlinear geometric kinematics as follows 36,37 : x Eliminating u (i) and v (i) from Eq. (14) and using Eqs (5) and (13) we have where A (i)* is the inverse of the membranous stiffness A (i) . Equation (15) is the compatibility equation for in-plane membrane deformation of Sub-laminate Ω (i) in terms of the stress function ψ (i) and the transverse deflection w (i) . Next, we examine the nonlinear flexural mechanism of the composite laminate. From the point view of shell theory, the nonlinear problem differs from the linear problem mainly on the treatment of the membrane forces. In linear problem, the membrane forces in the middle plane of the laminate are fixed, while in nonlinear problem they change with the increasing of the amplitude of the deflection and the influence of the membrane forces on the force equilibrium in transverse direction cannot be ignored. For this reason, the differential equation of equilibrium in transverse direction and the moment equilibrium for nonlinear bending of sub-laminate Ω (i) take the following forms 37 : where Q j (i) (x, y, t) is the transverse shear forces of sub-laminate Ω (i) , m (i) is the mass moment of inertia, q (i) (x, y, t) is the distributed load in transverse direction applied to Sub-laminates Ω (i) , and the changes of curvatures κ (i) (x, y, t) are given by Substituting the last two of Eq. (16) into the first equation and introducing Eqs (5), (13) and (17) we have In Eq. (18), the term m (i) w ,tt (i) represents the component of inertial force along with nonuniform motion. Thus, Eq. (18) is also known as the equation of motion for nonlinear vibration of Sub-laminates Ω (i) . Now, we need to consider the influence of the interaction penalty between sub-laminates. Due to the interaction forces q * (x, y, t) between sub-laminates, the distributed load acting on each sub-laminate q (i) (x, y, t) is not the same. The interaction forces transmit the external force from sub-laminate Ω (2) to Ω (3) , causing a reallocation of the external force q(x, y, t) between sub-laminates Ω (2) and Ω (3) : that is a net force of q(x, y, t) -q * (x, y, t) on sub-laminates Ω (2) and q * (x, y, t) on Ω (3) . With the available interaction forces q * (x, y, t) given in Eqs (10) and (11) the distributed load on each sub-laminate q (i) is reapportioned as follows: Delamination integration. Although we have established the governing equations of nonlinear vibration for each particular sub-laminate Ω (i) , the four sub-laminates are still independent with each other. Except the interaction effect between Ω (2) and Ω (3) , we need to impose constrains to make the four sub-laminates be the four parts of a single delaminated composite laminate. Therefore, we enforce the equilibrium conditions and compatibility conditions between the sub-laminates to unite them as a monolithic entity with a delamination integrated.
The equilibrium conditions require that the moments, the shear forces, the in-plane forces must be balanced at the fronts of delamination. The continuity conditions mean that the transverse deflection, the slope, and the in-plane displacements must be compatible. They are At x = l 1 (1) At x = l 1 + l 2 In Eqs (20) and (21) the moment M x (i) (x, y, t), the shear force Q x (i) (x, y, t), the in-plane displacements u (i) (x, y, t) and v (i) (x, y, t) are related to the transverse deflection w (i) and the stress function ψ (i) through the following equations, respectively, Equations (15), (18)(19)(20)(21)(22) constitute the general theory for nonlinear vibration of fiber-reinforced composite laminates containing an embedded delamination. The theory accounts for large geometric deformation in transverse direction, the in-plane membrane stretching, and the coupling interaction penalty at the interfaces of the delamination. It is highly nonlinear-not only because of the coupling of membrane mechanism and flexural mechanism within each sub-laminate, but also due to the coupling of the transverse deflection between sub-laminates.
Nonlinear dynamic stability analysis. As an application of the proposed theory, we analyze the free vibration of a delaminated composite laminate simply supported at its four edges. The laminate is made of graphite/epoxy composite plies with equivalent material properties of the longitudinal modulus E 11 = 140 GPa, the transverse modulus E 22 = 10 GPa, the shear modulus G 12 = 5 GPa, and the Poisson's ratios μ 12 = 0.3. It contains a though-width delamination symmetrically located at its center and is referred to a global coordinate system (x, y, z) centered at the left corner of the middle surface of the sublaminate, in which x and y are in the longitudinal and breadth directions and z is in the direction of the inward normal to the middle surface. In a real problem, mechanical behaviors depend on the dimensions of the physical quantities. To establish a general model with performances insusceptible to the dimensions of the physical quantities, we introduce the following non-dimensional parameters: where ω is the frequency of free vibration of the composite lamina. The free vibration of the composite laminate is an eigenvalue problem, thus the frequency of the free vibration ω is a characteristic value of the composite laminate and unaffected by external forces, such as the transverse load q on the surfaces or in-plane forces N ij at the edges. Under this circumstance, the membrane forces induced by the transverse deflection do not need to be considered, i.e. Equation (15) is not required 38 . Substituting Eq. (23) into (18) and setting ψ ,ij (i) and q to be zero, the governing equation for free vibration of the delaminated composite laminate is rewritten from Eqs (18) xxxx xxy y y y y y tt As the membrane forces and deformations are not included for the problem discussed herein, we need to remove these terms from equilibrium and continuity conditions in Eqs (20) and (21). The remaining formulas are the equilibrium of bending moment M x and shear force Q x and the continuity of the deflection w and it's slope w ,x , which are expressed in nondimensional forms as follows: (1, , ) (1, , ) (1, , ) where G j (i) and L j (i) are the coefficient constants of the global vibration mode and local buckling mode, λ G (i) and λ L are the modal parameter of global vibration and local buckling. The modal parameters λ G (i) and λ L are found from the governing equation as follows (See section Method for detail): where a i , b i and c i are constants determined by the geometric parameters and material properties of Sublaminate Ω (i) as follows: The coefficient constants for Sub-laminates Ω (2) and Ω (3) are also found to be related to each other as follows: where Equations (32)(33)(34)(35)(36)(37)(38) are the solutions of the vibration mode. We have mentioned in sub section of Interaction kinematics of sub-laminates that the solutions are valid for Region R. However Region R is just arbitrarily chosen. It denotes any one of the contact regions. If we use different contact area to analysis, we obtain the same solutions given in (32)(33)(34)(35)(36)(37)(38). It means that Eqs (32)(33)(34)(35)(36)(37)(38) are valid for all the contact regions. We know that the contact regions and the voids are sequentially distributed along the interfaces of the delamination. Then when we plot out Eqs (32) and (33) along the entire length of the laminate, we can identify the contact regions where w (2) > w (3) and the voids where w (2) < w (3) . Equations (32) and (33) manifests that the nonlinear free vibration of the delaminated composite laminate consists of global vibration and local buckling. In particular, each of the four sublaminate undergoes a global vibration, but the local buckling only occurs between Ω (2) and Ω (3) . In the rest of this section we will explain how the performances of the global mode are different from those of the local buckling.
Going back to Eqs (1-4), we note that Sub-laminates Ω (1) and Ω (4) possess the same material properties, that is A (1) = A (4) , D (1) = D (4) . We have mentioned at the beginning of this section that the composite laminate has a symmetrically located delamination at its center, in geometrical expressions l 1 = l 4 , l 2 = l 3 . Under such circumstance we have λ G (1) = λ G (4) from Eqs (34) and (36). It illuminates that the modal parameters of the global vibration for Sub-laminates Ω (1) and Ω (4) are the same. So do those for Sub-laminates Ω (2) and Ω (3) , as shown in Eq. (35). We now observe that the coefficient constants of the global vibrations for Sub-laminates Ω (2) and Ω (3) are almost equal to each other, as given in the first of Eq. (37). All the relationships rationalize that Sub-laminate Ω (2) and Sub-laminate Ω (3) exhibit an identical global vibration. This holds for Sub-laminate Ω (1) and Sub-laminate Ω (4) as well because of symmetrical delamination.
Based on above observation, we estimate the approximate values of λ G Eq. (34) are estimated to be proportional to K i 2 (h 0 /l i ) 4 in which K i is the number of composite ply of Sub-laminate Ω (i) , therefore they are roughly of the same magnitude order. Depending on the geometries of a particular delamination, the magnitudes of a i , b i and c i range from 10 −12 to 10 −7 . On the other hand, the nondimensional interaction factor k for composite laminate made of graphite/epoxy composite plies with defined material properties in the first paragraph of this section is approximated to be 10 −4 . With such properties, the numerical values of global vibration λ G (i) and local buckling λ L in Eq. (35) are approximated to be λ G (i) ∝ b i /a i and λ L ∝ (k/a 2 ) 1/4 , respectively. Considering that a i , b i , c i ∝ K i 2 (h 0 /l i ) 4 and k ≫ a 2 , it is easy to see that λ L ≫ λ G (i) . Our further investigations indicate (except very limit conditions when l i < 0.05a) that As the global vibration and local buckling are the effective components of the total free vibration mode, their leading order terms in Eq. (33) must be of the same order, i.e. ∝ λ λ . Using Eq. (39), the following relation is established Equation (40) illustrates that the amplitudes of global vibration are about 10 3 times of those of local buckling, that is to say, the global vibration is macroscopic and the local buckling is microscopic. In addition, it can be seen from Eqs (34)(35)(36) that the global vibration modes λ G (i) and the local buckling mode λ L are related to the nondimensional frequency ω and the nondimensional interaction factor k, respectively. Referring to Eqs (11), (23) and (35), it is found that the value of λ L is determinate for a particular case in which material properties and geometric parameters of the laminate and the delamination are all known. Determination of λ G (i) calls for the application of the equilibrium conditions and the continuity conditions in Eqs (27)(28)(29) and the boundary conditions in Eq. (30) from which a set of homogeneous equations of the indeterminate coefficients are derived and are expressed in a concise manner as follows: where ω R [ ( )] is 16-by-16 coefficients matrix containing the nondimensional frequency ω as the only unknown and Χ [ ] is a collection of the sixteen indeterminate coefficients given by The premise of solving the homogeneous Eq. (41) is the determinant of ω R [ ( )] to be zero, i.e.

Results
In order to study the performance of the nonlinear vibration of the delaminated laminates, we conduct numerical analysis. Based on the results presented above, we develop MATLAB programs to perform numerical calculations. In the implementations, the laminate has a geometry size of 2 m × 2 m × 0.12 m and is composed of 60 plies made of T300/QY8911 carbon fiber reinforced composite material with the density ρ = 2150 kg/m 3 stacked in a sequence of [0°/90°/0°] 20 with each single-layer thickness h 0 = 0.002 m. The equivalent material properties of each composite ply are still the same as those in Nonlinear dynamic stability analysis, i.e. E 11 = 140 GPa, E 22 = 10 GPa, G 12 = 5 GPa, and μ 12 = 0.3. The data obtained enable us to clearly reveal the nature of the free vibration of composite laminates containing a buried delamination. The profile of the vibration modes of the delaminated laminate at the cross section of x = a/2 is shown in Fig. 3a obtained from MATLAB programs for analytical solution given the non-dimensional delamination length η = l 2 /l = 0.5 and the non-dimensional delamination depth τ = h 2 /h = 0.3. As a comparison, we also perform finite element simulation using ABAQUS software. Unlike the shell theory in theoretical analysis, we use solid elements and General Contact Interaction Module in ABAQUS analysis in order to incorporate the anti-penetration penalty at the interfaces of sublaminates Ω (2) and Ω (3) -above and below the delamination-to prevent them from penetrating with each other during vibration. Figure 3b is the free vibration of the delaminated laminate from ABAQUS analysis for the same η and τ used in Fig. 3a. The profile of the vibration mode in Fig. 3b makes it clear to us that the delaminated composite laminate vibrates as a single whole-sublaminates above and below the delamination vibrate identically. This is exactly the same as what we predict in theoretical analysis given in Fig. 3a. The observation of consistent global vibration mode in ABAQUS justifies our analytical solutions and makes our analysis effective.
We have discussed that except the consistent global vibration the nonlinear vibration of the delaminated composite laminate includes a local buckling. Figure 3c shows the buckling modes in x-direction at cross section y = 1 m for both Sublaminates Ω (2) and Ω (3) . Distinct from their consistent global vibrations, their local buckling modes of Sublaminates Ω (2) and Ω (3) are not identical-not only the buckling amplitude of Sublaminates Ω (2) is much larger than that of Sublaminate Ω (3) due to its weaker stiffness, but also their amplitudes are completely reversal with each other. It is also found from Fig. 3c that the two sublaminates buckle in modes of equal but Scientific RepoRts | 7: 4468 | DOI:10.1038/s41598-017-04570-3 quite short wavelength. With these features, once the delaminated laminate undergoes vibration, local voids are formed on the interfaces at places where upper Sublaminate Ω (2) and lower Sublaminate Ω (3) buckles upwards and downwards respectively. Between the voids the two sublaminates buckle in the opposite directions. Interaction forces are produced in these regions. With these characteristics, voids and interaction areas are developed alternately at the interfaces of the delamination, as illustrated in Fig. 3c. It is remarkable that the magnitudes of the local buckles are only about 10 −3 of the magnitudes of the vibration modes, so the voids are microscopic, invisible to the naked eyes.
For a more comprehensive understanding of the situation, we further investigate cases of different length of delamination with/without the interaction penalty. Figure 4 shows a comparison of the free vibration frequency of the delaminated laminate with/without interaction penalty enforced at the interfaces of the delamination from ABAQUS analysis. The results reveal that the delaminated laminate vibrates in modes of much high frequency when the penalty conditions are enforced. The reason for such difference is that the penalty conditions link the separate sublaminates together as a unitary one, thus greatly enhance the stiffness of the delaminated laminate.
We further carry out a parametric study to investigate the influence of the delamination properties on the vibration frequency. Variation of the nondimensional frequency ω with respect to the nondimensional delamina-  tion length η are given in Fig. 5a given the nondimensional delamination depth τ = 0.1 and τ = 0.2. As a comparison, the results from ABAQUS analysis are drawn in the diagram too. Both the analytical solutions and ABAQUS results show that the nondimensional frequency ω decreases with the increasing of the nondimensional delamination length η. It is a common sense that the smaller the stiffness of a laminate, the larger its amplitude and the lower its frequency. The delaminated laminate is actually a single whole with reduced stiffness due to the existence of the delamination. As the length of the delamination enlarges, the stiffness of the laminate continues deteriorating, thus leading to the decreasing of the frequency.
We also plot variation of the nondimensional frequency ω with respect to the nondimensional delamination depth τ in Fig. 5b given the nondimensional delamination length η = 0.6 and η = 0.8 for both analytical solutions and ABAQUS results. Likewise, the nondimensional frequency ω reduces as the value of nondimensional delamination depth τ rises. Such trend is still owing to the degradation of the stiffness. We are clear that the stiffness of the delaminated laminate depends on the stiffness and the geometries of the four sublaminates. For a specific length of the delamination, the geometries of the four sublaminates as well as the stiffness of Sublaminates Ω (1) and Ω (4) are invariable. Under such circumstance, the stiffness of the delaminated laminate is determined either by Sublaminate Ω (2) or Sublaminate Ω (3) , whichever has a larger stiffness. As the delamination moves from the surface to the middle of the laminate, the stiffness of Sublaminate Ω (3) continues to decrease but still exceeds that of Sublaminate Ω (2) , thus the stiffness of the delaminated laminate is continuously reduced. At the middle plane, Sublaminates Ω (2) and Ω (3) possess equal stiffness and the delaminated laminate has the lowest stiffness. Once the delamination is in the lower part of the laminate, the stiffness of the delaminated laminate is mainly controlled by Sublaminates Ω (2) instead of Ω (3) and increases with the delamination depth.
Comparison of the nondimensional frequency ω with ABAQUS results is presented in Fig. 5b, too. One can see from Fig. 5 that the results from finite element simulation and theoretical analysis have some deviations. The nondimensional frequencies ω from ABAQUS simulation are numerically slightly larger than the analytical predictions. Besides, the declining trends of ω with respect to η from the two approaches are somewhat different, as shown in Fig. 5a. These deviations are due to the facts that we use the two-dimensional plate and shell theory to carry on the theoretical analysis but a three-dimensional solid model to carry on the finite element analysis. To prevent the penetration between sublaminates Ω (2) and Ω (3) in ABAQUS simulation, hard contact property is defined in interaction module to impose constrains on the bottom of Ω (2) and on the top of Ω (3) . However such application is practicable only to solid model in which shear modulus in transverse direction are necessary. As the shear deformation in transverse direction is omitted in theoretical analysis, inevitably, there are deviations between the theoretical and finite element results. Nevertheless, the comparisons are still adequate in validating our theoretical predictions.

Discussion
In summary, we present an analysis of the nonlinear vibration of composite laminates containing a buried delamination through the newly developed multiscale framework. The delaminated laminate is treated as four independent sublaminates with granted restrains of forces and deformations at the delamination fronts to reunite them as single whole. By incorporating the anti-penetrating interaction kinematics at the interfaces of the delamination into the nonlinear theory for vibration of composite laminate, the nonlinear vibration of delaminated composite laminates is explored. We consider the multiscale compositions of the laminate and identify two failure modes of global vibration and local buckling. The modes of the global vibration are macroscopic and compatible among the four sublaminates with exactly the same frequency. In particular, sublaminates above and below the delamination undergoes a fully consistency vibration. On the other hand, the local buckling is microscopic, occurring only between sublaminates above and below the delamination with completely reverse amplitudes that are about 10 −3 times smaller than those of the vibration mode. By developing numerical MATLAB program, we investigate the effects of the delamination on the intrinsic behaviors of composite laminate over a wide range of the delamination length and the delamination depth. The free vibration frequency is shown to monotonic decrease upon increasing of the delamination length and to decrease then increase as the delamination depth increases. We perform finite element simulations for the free vibration of delaminated laminates with/without interaction penalty enforced at the interfaces of the delamination, which suggest that the interaction penalty influence not only the vibration frequency but also the vibration mode. For laminates of the same geometry and delamination parameters, the one without the interaction penalty has much lower vibration frequency than that with the interaction penalty. The difference is reflected on the vibration mode in which sublaminate above and below the delamination penetrate with each other when the interaction penalty are not enforced. The characteristic behaviors of the delaminated laminate are consistent in both theoretical solutions and finite element results, thereby demonstrating the effective of the analysis. Now Eq. (53) can be easily solved. The multiscale solutions of λ (2) are given in Eq. (35). With the available solutions of λ (2) , we can formulate the vibration mode in x direction f x ( ) (2) of Sublaminate Ω (2)  Substituting the obtained f x ( ) (2) into Eq. (48), we can determine the solution of f x ( ) (3) : it is quite similar to Eq. (55) and has a close relation with solution of f x ( ) (2) : λ (3) = λ (2) (see Eq. (35)) and G (3) ≈ G (2) and L (3) = −ιL (3) (see Eqs (37) and (38)).