Estimation of Local Delamination Buckling in Orthotropic Composite Plates Using Kirchhoff Plate Finite Elements

We analyse the buckling process of composite plateswith through-the-width delamination and straight crack front applying uniaxial compression.We are focusing on themixedmode buckling case, where the non-uniform distribution of the in-plane forces controls the occurence of the buckling of the delaminated layers. For the analysis, semi-discrete finite elements will be derived based on the Lèvy-type method. The method of harmonic balance is used for taking into account the force distribution that is generally non uniform in-plane.

Many researchers have studied the field of delaminated composite plate buckling experimentally, analytically using finite element method (FEM).Chai et al. [25] have developed an analytical one-dimensional model for simulating buckling of delaminated beam or plate specimens.The delamination was assumed to be close to the surface (thin-film approach), and mixed buckling modes have been determined.On a beam-plate model it was shown that the transverse shear effect reduces the critical buckling loads [26].The effect of the boundary conditions was investigated on a 1D model, taking into account the large deformations, and the results were compared with the thin-film approximation and experiments [27].Anastasiadis and Simitses [28] presented a modified 1D model for improving the critical buckling loads.They used artificial springs along the crack line for preventing crack opening, which would result in inadmissible mode shapes [29].Ovesy et al. [30] developed a layerwise theory based on the First-Order Shear Deformation Theory (FSDT) for analysing the postbuckling behaviour of multiple delaminated laminates.In order to prevent the inadmissible mode shapes they used contact constraints on the delaminated area.For the double delaminated beam-plate model an analytical solution was introduced by Shu using a constrained model for the global buckling calculation [31].Kim and Kedward [32] used the Classical Laminate Plate Theory (CLPT) for modelling a delaminated rectangular plate and analysed its behaviour with respect to the buckling using the Navier method.For the global stability analysis the delaminated region was treated as a reduced stiffness zone.The local stability analysis was carried out on a clamped plate assuming the same load along the local part as for the global model.The delamination buckling and delamination growth of cross-ply laminates were examined, induced by low velocity impact in [33].The results of the optical measurement were compared with X-ray nondestructive evaluation results.Wang and Lu [34] examined the failure mechanism of near the surface embedded delaminations under compression load using the energy method and experimental testing.Through-the-width delaminated composite laminates subjected to compressive loads were investigated, using the Rayleigh-Ritz method and CLPT theory by Kharazi and Ovesy [35].Experiments on multiple delaminated composite plates and a FEM study using ANSYS were carried out by Aslan and S ¸ahin [36].Kharazi developed a layerwise theory based on the FSDT where the prevention of the penetration of the delaminated layers was solved by using the penalty method and artificial springs [37].The mixed mode buckling was captured with a contact model in the paper of Hwang and Liu [38].A layerwise FEM solution using Abaqus and VICONOPT is given for composite plates with embedded rectangular delamination in Damghani et al. [39].Marjanović and Vuksanović [40] used the Generalized Layerwise Plate Theory of Reddy and built a finite element (FE) model which was capable of handling more embedded delaminations in the composite structure.A 3D FEM model was developed for analysing the buckling behaviour of embedded and through-the-width delaminated plates in [41].
For the global model an analytical solution was presented in Juhász and Szekrényes [42].In that paper it was shown that the in-plane force distribution greatly influences the buckling behaviour of the plate.In this paper a FE model is introduced, which concerns a plate with two simply supported opposite edges so the plate problem cannot be solved using beam theories.The other two boundary conditions (BCs) can be arbitrarily chosen.The plate is loaded with uniaxial compression.The critical in-plane loads for the global buckling are calculated using a constrained model [29,31].Using this approach there is no need for using constraints along the delamination, which results in a very efficient calculation method contrary to the widely used contact models.The continuity between the delaminated and undelaminated parts is maintained using special transitional elements.This model is capable of predicting the global critical buckling loads and determines the corresponding in-plane force distributions.Using the solution of the global model the local stability is analysed separately using a local FE model for each delaminated part.Because of the nonuniformity in the distribution of axial forces the stability analysis was carried out using the method of harmonic balance [43].It was shown that these effects greatly influence the local buckling loads.In this paper the effects of different type of boundary conditions are analysed in a numerical example where the length of the delamination is varied, and the local and global critical amplitudes are calculated, which results in a stability map for the different BCs.It is shown that the different BCs greatly influence the buckling loads and the corresponding mixed mode buckling shapes.Using our method and the resulting stability maps the failure process of a compressed delaminated plate can be estimated, and it can be determined at which point will the plate buckle globally or locally, or in mixed mode case.

Model Creation
Let us consider a layered plate with several orthotropic plies and a closed through-the-width delamination (see Figure 1).The plate length is , the length of the delamination is , and the thickness of the plate is .The plate is composed of two equivalent single layers (ESL) with thicknesses   and   , respectively.
The edges parallel to the -axis are simply supported, and the constraint of the other two can be arbitrarily chosen.The plate is subjected to uniaxial compression in the midplane.We are considering thin structures; therefore, the solution is based on CLPT [44][45][46][47]: where û0 and V0 are constant parts in the in-plane displacement functions of the undelaminated portions, which arise because of the kinematic coupling on the interface as Figure 2 shows.  0 and V  0 are constant through the thickness in the top and bottom parts, where  takes either "" or "," respectively (see Figure 2).Moreover,  0 is the deflection of the plate parts, respectively.Along the undelaminated parts we assume perfect adherence and no crack propagation at the crack tips.Therefore, the continuity of the undelaminated portions is ensured using the system of exact kinematic conditions (SEKC) [48]: And the in-plane displacements in the global reference plane are From the literature it is well known that a free model allows the intersection of the delaminated portions into each other, which results in kinematically inadmissible mode shapes and wrong critical loads [29,37,[49][50][51].Avoiding this a constrained model is used, where the deflection of the delaminated top and bottom portions is common.Using these assumptions the strain fields of the undelaminated portions can be given as And for the delaminated portions the strain field will be where 0 means the constant and 1 means the linear part of the strain field in terms of the  coordinate [52,53].
For expressing the stresses of the composite laminate with the strain fields we have to use the constitutive equation of composite laminates (CEL) [52][53][54][55]: where the matrices A, B, and D can be calculated based on the literature.Using these the strain energy density of the delaminated and undelaminated portions can be given as [17,52,56] where the in-plane force and moment resultants are depending on the strain fields based on (6).

Finite Element Discretization
The finite element discretization concept can be seen in Figure 3.
As we have written in the previous section the in-plane displacement fields are continuous on the undelaminated portions and independent of each other in the delaminated top and bottom portions (see Figure 3).Therefore, between the delaminated and undelaminated portions special transition elements were used on sections 2 and 2 which capture the crack tips.These elements ensure the continuity of the displacements [57].Because of the opposite simply supported edges the Lèvy-type method is applicable and the displacements fields in the  direction can be given by the terms of Fourier series, Hosseini-Hashemi et al. [58], Bodaghi and Saidi [55], Thai and Kim [59], Szekrényes [48,60], Thai and Kim [61], and Nguyen et al. [62]: where  = / and  0 (),  0 (),  0 () are the amplitudes in the  direction.As the displacement in the  direction is given it is enough to discretize the model along , so we can write the nodal displacements as

Crack tip Crack tip Delamination
where u  is the nodal displacement vector and R  is a diagonal matrix containing the corresponding trigonometric terms from (8).This results in a semidiscrete finite element, which is capable of modelling plate-like buckling.For the analysis we need to derive the material and the geometric stiffness matrices of the different elements of Figure 3.
Integrating the strain energy density the material stiffness matrix of an element can be derived based on the Hamilton principle [63][64][65][66]: where K  is the material or general stiffness matrix and Ω  is the element domain.Based on the literature the geometric stiffness matrix with respect to the axial compression can be given as [63,67] As the in-plane   forces along the delaminated region will not be zero, for the local stability analysis we also need the geometric stiffness matrix with respect to the   load: In (( 11)-( 12)) the  0 is the applied load on the element, N  is the vector of interpolation functions of the element, and B  is the derivative of the vector of interpolation functions.

Element of the Nondelaminated Parts.
For the undelaminated parts an 8DoF element is used: where   is the rotation of the cross-section in the node.The matrix of the trigonometric coefficients can be given as For the in-plane displacements linear interpolation functions were used, whereas for the transverse deflection a third-order function was applied: where   ,   , and   are constant coefficients and  varies between 0 and 1.Their value can be calculated using the CLPT from the following equations: 1 where the derivation is carried out with respect to the dimensionless  coordinate.By solving (( 15)-( 21)) the vector interpolation functions can be obtained [63]: where  can be , V, or .Substituting the discretized displacement fields into (10) the strain energy can be given as Carrying out the integration over  the element material stiffness matrix can be obtained.
Using the interpolation functions from ( 11) the geometric stiffness matrix can be derived.

Element of the Delaminated Parts.
Because of the constrained model the transverse deflection is common, but the in-plane displacements are independent in the delaminated portion, which results in a 12DoF element: The R  matrix can be composed based on the nodal displacement vector.For the in-plane displacement the same linear interpolation functions were used, given by (( 15)-( 16)), and for the transverse deflection the third-order function was applied in accordance with (17).
As the in-plane displacements are independent, the potential energy has to be evaluated for both the top and bottom parts and the material stiffness matrix can be derived from the sum of the potential energies based on (23).

The Transition Elements.
In accordance with Figure 3, the elements of sections 2 and 2 ensure the kinematic continuity between the delaminated and undelaminated portions.The vector of displacements of the 2 element is And for the element denoted by 2 we have Based on the kinematic continuity the following 4 equations can be written based on the applied plate theory for the 2 section: The equations take similar form for the 2 section: Using the equations above and ( 20) and ( 21) the vector of interpolation functions can be obtained.The stiffness matrices can be calculated on the same way as it was shown before.

Stability Analysis
Based on Section 3 the structural matrices of the global model can be obtained.After applying the selected BCs on the  = 0 and  =  edges the critical loads and the corresponding eigenvectors can be calculated as The corresponding global mode shapes and the resultant inplane force distributions can be obtained using the vector of interpolation functions by (22) and the CEL given by (6).
Because of the in-plane resultant forces at the crack tips the plate is able to buckle locally along the delamination.The local stability is analysed individually for the top and bottom delaminated plate portions, assuming plates with built-in end BCs along the crack tips.The local stability is affected by the distribution of the in-plane forces (see Figure 4).For the local FE model we derived the elements of the individual top and bottom layers, using the same method as for the elements of the global model.The nonuniform resultant in-plane forces of the global model are evaluated for every element at the middle.These values were normed with the value at the crack tip, and the element geometric stiffness matrices were multiplied with these values, taking into consideration the distribution along .Because of the simply supported edges the plate will have a half wave shape along the width, which results in the fact that the load along the crack tip will not be uniform (see Figure 4).Taking this aspect of the problem into consideration we applied the method of harmonic balance and wrote the Fourier series of the nodal displacements [43]: where   are constant coefficients and  is the vector of displacement values.Taking this back into (29) and applying some trigonometric identities we can obtain a system of equations in matrix form: The critical values and the corresponding mode shapes can be calculated from ( 31) and (29).
For validation purposes the model was solved using Abaqus.The plate is made by carbon/epoxy material using the following layup order: [±45 ∘ ; 0 ∘ ; ±45 ∘ ].Engineering constants of the layers are detailed in Table 1.The series expansion in (30) was carried out for two terms.Along the  direction the plate was discretized using 14 elements, to capture the higher order mode shapes.The obtained critical values from (31) are ∼40% higher than the loads of the problem with constant distribution along .The top ESL of the example in Section 6 was checked assuming constant force distribution along .The width of the plate was 100 mm, and the length of the plate was 105 mm.The S4R shell element was used for the analysis with 1 mm element size.The results show good agreement with the present calculations (see Table 2).

Boundary and Continuity Conditions
In this paper the process of loss of stability is determined by using a displacement controlled model based on Section 2. For solution, the Lèvy-type method is used with the statespace approach [52].From (7) using Hamilton's principle, the governing PDEs of each section can be derived [52].Applying (8) the obtained ODEs can be rearranged into the state-space model [52,68]: where Z is the state vector.The general solution of ( 32) is [52,68,69] where K  is the vector of constants {  }.At the crack tips we have to define 10-10 continuity conditions (CCs) between the plate portions A, B and B, C. Because of the closed delamination (see Figure 1) the so-called Mujumdar conditions have to be used for fitting the   moment and the Kirchhoff equivalent shear force [29]: where  can take either,  1 or ( 1 + ) respectively, and {  ,   ,   ,   ,   ,   } depends on { 0 ,  0 ,  0 } and their derivatives. can take 1, 2, or 3 depending on the sections, which will be fit.In the BCs an  0 axial displacement at  = 0 has to be prescribed, and there is no other load.Substituting the solution of the state-space model into the BCs and CCs a system of inhomogeneous equations can be obtained: which can be solved for the  all  constants.Using (33) we can get the displacement functions, and the in-plane forces can be calculated using (6).
Using this model we calculated the arising forces at the edge of the plate and at the crack tips with respect to the axial displacement  0 .The critical values of the global and local stability analysis were compared with these results.

Criterion of Constant Arc
Length.All of the mode shapes were calculated with a maximum amplitude of 1 mm and scaled to fit the physical requirements.The amplitudes of the global and local modes were controlled using an arc length criterion [57].This means that the arc length of the superimposed eigenshapes minus the axial displacement has to be equal to the length of the plate or the delamination: where    is the scale factor for the mode shapes and    () is the buckled shape of the th buckling mode.For global mode shapes Δ =  0 .For local mode shapes it is the signed sum of the axial displacements at the left and right crack tips.In case of mixed mode buckling Δ is where    () is the axial displacement of the global model for the th mode scaled with the  Amp. ( 0 ) amplitude and  Static  ( 0 , ) is the static axial displacement.

Superposition of the Mode Shapes.
Based on the linear relationship between the in-plane normal force and the axial displacement of the displacement controlled model, critical axial displacements can be calculated for the critical loads.From these values the amplitudes are rising linearily: We assume that the dominant in-plane force distribution is determined by the first global mode.Therefore, local buckling in mixed mode case occurs only if the critical values of the local modes, which belong to the first global mode, are reached.

Numerical Example
In this section we adopt the method on a carbon/epoxy layered plate.The ply order of the plate is [±45 ∘ ; 0 ∘ ; ±45 ∘ ].The plate consists of 9 layers.The corresponding material data can be found in Table 1.The plate is symmetrically delaminated, and its geometric data is presented in Table 3.The stiffness matrices of each single layer were determined based on the elastic properties given by Table 1.The analysis was carried out with  = 1 condition in (8).The order of the matrix in (31) was set to 2. The plate was discretized using 12 elements in all sections and 1-1 additional transitional elements were used at the crack tips.The position of the delamination was set above the 5th layer.At the edges  = 0 and  =  the same simply supported (S-S) or built-in (B-B) BCs were used.The length of the delamination was varied from 10 mm to 100 mm.The global critical forces with respect to the delamination length can be seen in Figure 5.
It can be seen that the obtained critical loads of the builtin plate are higher, but as the length of the delamination increases the effect of BCs gets less significant.The critical amplitudes of the local top and bottom delaminated portions can be seen in Figures 6-7.
As it can be seen the local critical values are higher in the simply supported cases.This is because different eigenshape belongs to the different BCs which results in different inplane force distribution.Again as the delamination length increases the effect of the BCs gets less significant.Using the displacement controlled model the critical axial displacements can be calculated for each critical amplitude.Based on this calculation stability diagrams can be obtained with respect to the axial displacement and the delamination length (see Figures 8-9).
On both pictures below the blue line the plate is stable.In the orange region the plate buckles globally, in the green region it buckles globally and the crack opens as the local top plate loses its stability, and above the green line the delaminated bottom portion buckles too.It has to be remarked that in the B-B case the bottom buckles only at    higher axial compression; therefore, the green line is outside the range shown in Figure 9.The maximal critical amplitude was set to 2 mm.It can be seen that the built-in end plate is more stable, and its bottom part does not lose its stability up to the maximal axial displacement, whereas the simply supported plate loses its stability on smaller amplitudes.It can be noticed that as the delamination length increases the point of the global and local stability loss of the top plate gets close to each other.The presented critical loads are the first critical amplitudes.But if the plate is weak against uniaxial compression, higher order mode shapes are also feasible.These mode shapes can be superimposed using the arc length criterion.In the following we will show the process of stability loss of the simply supported and built-in end plates with 100 mm delamination length.The global critical amplitudes for the two types of BCs are listed in Table 4.For these values the critical axial compressions can be determined based on the displacement controlled model.The resulting forces with respect to the axial displacement for the simply supported case are shown in Figure 10.On the same way the critical axial displacements of the built-in end plate can be determined.Whereas the critical loads are higher than in case of simply supported BCs the critical axial displacements of the first 2 modes are smaller and only the third mode appears at higher displacement: 0.14 mm, 0.15 mm, and 0.33 mm.The maximal axial compression was chosen in both cases for the 120% of the third mode.The critical values of the delaminated N Amp.
N Amp.portions were calculated for the local buckling case where the nonuniform distribution of the in-plane forces does not count, but the calculated critical axial displacements were higher than the critical axial displacement of the first global mode; therefore, the plate loses its stability first globally.From Figures 11 and 12 it can be seen that because of the different BCs different mode shapes appear.For the mixed mode buckling the local critical values were calculated for both cases using the nonuniform force distribution of the global modes.Here we present only the critical loads of the realizing local modes (see Table 5).As it can be seen only the first two local modes appear in both cases.The third mode would only appear at higher axial compression.At the built-in end case the local modes, calculated with the force distribution of the second global mode, are not present during the stability loss, because the critical values of these modes are much higher.The plate was also examined for the   forces, but according to the results no stability loss occurs with respect to the   forces at the crack tip.In accordance with Figures 8 and 9 the delaminated bottom part does not lose its stability at the selected maximal axial displacement.
The shapes of these modes were calculated with the inplane force distribution resulting from the corresponding global modes and were superimposed using the arc length W Amp.

Figure 12:
The global mode shapes of the built-in end case.Note that the distributions involve a half sine wave in the  direction.On the superimposed shapes it can be seen that the dominant part of the solution is always the global and local first modes, but the higher order modes influence the shape slightly.

Conclusion
In this paper the buckling process of a delaminated layered plate was investigated.The formulation of the problem is based on the system of exact kinematic conditions (SEKC) by cutting the plate in the plane of the delamination and forming the continuity conditions.The problem was solved using FEM with self-developed semidiscrete finite elements.The model contains special transitional elements which ensure the kinematic continuity between the delaminated and undelaminated portions.The delaminated region was modelled as a constrained section in the global model; therefore, there is no need for using contact along the delaminated area, which results in a calculation efficient and simple method for the estimation of the global critical buckling loads and the corresponding shapes.The local behaviour of the delaminated portion was analysed by a separate FE model.For the consideration of the nonuniform in-plane force distribution the method of harmonic balance was used.On a numerical example the effects of the simply supported and built-in end BCs were determined with respect to the delamination length.It was shown that the BCs are influencing not only the critical loads but also the corresponding global mode shapes.Because of the different global mode shapes the local behaviour of the delaminated portions is different, as the inplane force distributions differ significantly.This results in the fact that whereas the simply supported plate buckles globally W Amp. W Amp.
W Amp. W Amp. W Amp. at lower values this configuration is more stable locally than the built-in end configuration.It was also shown that this effect is more significant if the delamination length is small.Stability diagrams with respect to the axial displacement and the delamination length were given where the global and mixed mode stability loss cases were shown.At one delamination length the process of stability loss was presented for both BCs.Here the effect of the BCs and the nonuniform in-plane force distribution can be seen.This nonuniform distribution was not observed with respect to the different type of BCs in the literature, and we can state that it greatly alerts the buckled shape of the delaminated layers.

Figure 1 :
Figure 1: Simply supported layered plate with delamination subjected to uniaxial compression.

Figure 2 :
Figure 2: Cross-sections and deformations of the top and bottom elements of an unsymmetrically delaminated orthotropic plate.

Figure 3 :
Figure 3: Finite element discretization and nodal DOFs of the layered plate with delamination.

Figure 4 :
Figure 4: Model of the local stability analysis (top or bottom part) (a).Example for the distribution of the in-plane normal force along the  direction (b).

Figure 5 :Figure 6 :
Figure 5: The global critical amplitudes with respect to the delamination length.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: The local critical amplitudes of the bottom plate portion with respect to the delamination length.

Figure 10 :
Figure 10: The static   and   curves and global critical forces and the corresponding axial displacements of the simply supported case.

Figure 11 :
Figure 11: The global mode shapes and the corresponding in-plane force distributions of the simply supported case.Note that the distributions involve a half sine wave in the  direction.

Table 1 :
Elastic properties of single carbon/epoxy composite laminates.

Table 2 :
Difference between the critical amplitudes of the constant and sine loaded plate, and the difference between the two types of loads of each method (Δ), and the difference between the results of the ABAQUS model with sinusoidal loading and the results of the present method (Δ sin ).The dimensions of the results are in N mm −1 .

Table 3 :
Geometric parameters of the plate modelled for the numerical examples.

Table 4 :
The global critical buckling loads in N mm −1 .