Stability Analysis of Variable Stiffness Composite Laminated Plates with Delamination Using Spline-FSM 1

The dynamic behavior of variable stiffness composite laminated (VSCL) plate with curvilinear fiber orientation subjected to inplane end-loads is investigated. A variable stiffness design can increase the laminated composite structural performance and also provides flexibility for trading-offs between various structural properties. In each ply of the VSCL plate, the fiber-orientation angle assumed to be changed linearly with respect to horizontal geometry coordinate. The spline finite strip method based on both classical as well as higher order shear deformation plate theories is formulated to explain the structural behavior. The panel is assumed containing internal square delamination regions with friction and contact conditions at delaminated interfaces are not considered. In order to demonstrate the capabilities of the developed method in predicting the structural dynamic behavior, some representing results are obtained and compared with those available in the literature. The effects of change in curvilinear fiber orientation angles on the structural stability is studied. The obtained results show very good conformity in comparison with those exists in the available literature.


INTRODUCTION
Aeronautical, space and marine structures are among disciplines where the least structural weight besides providing high available strength must be achieved.Thus, thin-walled structures will usually come to play.Every structural component especially those thinner in thickness under in-plane harmonically varying excitation added to a constant mean load, meets situations where the instability conditions may appear.The amplitude of the dynamic instability load may even be lower than the Latin American Journal of Solids and Structures 14 (2017) 528-543 value corresponding to static bifurcation point.These excitation conditions are prevalent in case of mechanical structures as well as fluid-structural interactions.
The traditional composite designs consider the composite lamina properties to be constant throughout the entire ply by the usage of straight and uniformly spaced fibers.This type of construction provides constant unchanged stiffness throughout the whole lamina.A ply with variable mechanical properties could be achieved by changing the fiber orientation angle with respect to the locality.With the automated fiber placement technology, it is possible to fabricate composite plies with variable fiber orientations within their geometrical domain.As a result of changed fiber orientation, the ply achieves variable directional stiffness through the geometry and may be called as a variable stiffness composite laminate (VSCL).On the other hand, a widespread defect of composite laminated structures is the debonding of layers called delamination.Delamination occurrence causes total strength reduction of the structure and activates low energy local instability and failure modes.So it is of high importance to calculate the stiffness reduction of a delaminated structure.Therefore, it is essential to study the various effects of loading and delamination on the dynamic characteristics and response of layered plates under static and periodic in-plane loads.
The first reported studies on curvilinear fiber VSCL plates can be traced back to the works by Hyer et. al (1991) and Gurdal and Olmedo (1993).Parhi et al. (2001) presented dynamic analysis of a squared plate with delamination.The finite element equations were developed based on first order shear deformation theory (FST) for a developed eight-node isoparametric element.Some parametric studies on the dynamic behavior of delaminated plates in case of various boundary conditions, different lay-ups and changing geometries were presented.Hu et al. (2002) studied the vibration of moderately thick laminated plates containing delaminations using FST finite element method.The effects of delamination on changing the vibration behavior of the geometry were inspected.Yang and Fu (2007) examined the parametric instability of a thin-walled laminated cylinder with delamination zone.The Rayleigh-Ritz method besides using of Heaviside-type displacement functions were utilized.The problem governing Mathieu equations were solved through Bolotin's approximation method and the effects of external excitation amplitude, delamination location and size, and the material properties on the natural frequencies and the fundamental instability regions were studied.In more recent years, Akhavan and Ribeiro (2011) studied the free vibration of VSCL plates with curvilinear fibers using the third-order shear deformation theory.Ovesy et al. (2014) developed a layer-wise spline finite strip method (FSM) based on FST and analyzed the parametric instability problem of laminated plates containing through the width delamination.The delamination region simulated using step displacement approximation functions.The results of change in delamination size and position in static buckling, natural frequencies and parametric instability problems were investigated.Tornabene and co-workers (2015) investigated higher-order structural theories for the static analysis of doubly curved laminated composite panels reinforced by curvilinear fibers using the generalized differential quadrature (GDQ) method.They extract the strain as well as stress distribution through the thickness direction.Mohanty et al. (2015) analyzed the dynamic instability of delaminated plates under in-plane harmonic loading.The FST finite element method was utilized.The effects of delamination position, lay-up, material orthogonality, geometry and loading on the instability characteristics were evaluated.The author of the current manuscript in some earlier publications (2010)(2011)(2012)(2013), developed semi-analytical as well as spline finite strip formulations and Latin American Journal of Solids and Structures 14 (2017) 528-543 investigated the parametric instability problem of flat and curved shell panels with and without internal cutout and longitudinal stiffeners under uniform in-plane loadings.
In this paper the static as well as dynamic stability behavior of moderately thick laminated flat panels containing delamination subjected to inplane loadings has been investigated.The laminate assumed to be variable stiffness due to curvilinear fiber placement that changes linear in the panel longitudinal direction.The in-plane loading is assumed to change harmonically with time.A Bspline version of FSM has been developed.The formulations are based on both the classical thin plate theory and the Reddy type higher order shear deformation theory in order to include the transverse shear stresses effect in case of moderately thick structures.The governing equations are derived using full energy concepts on the basis of the principle of virtual work.The dynamic behavior including natural frequencies as well as instability load-frequency margins are extracted utilizing the Bolotin's first order approximation followed by an eigenvalue analysis.Some representative problems are numerically studied and compared to those in the literature wherever available.To the best of the author's knowledge, this is the first application of B-spline FSM to the problem.Moreover, many studies have been devoted to the curvilinear fiber VSCL panels, while there are few published works on the parametric instability of curvilinear fiber variable stiffness composite laminated panels with delamination.

FORMULATION
The assumed model includes a flat square laminated panel with a squared through-the-width or embedded delamination zone.The panel laminates are made from curvilinear (linearly changed) fiber orientations.The geometry is divided to a number of longitudinally adjacent finite strips.Figure 1 shows a sample numerical mesh of a geometry of width b, length L and total thickness t.The geometry is made from number of strip elements of length L and width bs.A typical embedded single delamination zone is also indicated in the figure .A uniform loading as prestress is assumed on finite strips.The loading is consisted from a non-changing (static) component and a harmonically changing (dynamic) component which are indicated using S and D superscripts, respectively.The loading as a function of the static buckling load could be expressed as, With w , S a and D a as the excitation frequency, static part coefficient and dynamic part coefficient, respectively.
The assumed model displacement field based on Reddy-type third order shear deformations in thickness direction (zero shear at top and bottom surfaces), may be expressed as, Latin American Journal of Solids and Structures 14 (2017) 528-543 Where , , u v w are the displacement components of any arbitrary point, 0 0 0 , , u v w are the corresponding displacement components at the strip mid-surface, and , x y b b are the rotations around y and x axis, respectively.The pointer C is set to 1 and may be set to 0 in case of classical plate theory assumptions.The mid.-surface displacement field is approximated as multiplication of longitudinal-directional approximation functions.In longitudinal direction of a finite strip, the summation of a series of B3splines are employed while in the strip width, inplane linear Lagrange functions in conjunction with out of plane third order Hermitian ones are chosen (Fazilati and Ovesy, 2013).Any type of boundary constraints (free, hinged, clamped) may be implemented according to the approximation displacement functions chosen.
The linear strains for flat geometry are calculated through, , .( ) where ',' defines a differentiation operator.Substituting the displacement functions (equation 2) into the strain equations (equation 3) leads to the strain field as:  Solids and Structures 14 (2017) 528-543 where the strain coefficients are defined as: -, + , , ,-, -, , , 0,0,0 , , ., -, - Solution of the instability problems is sought through the principle of virtual work.The total energy of a strip is defined as summation of kinetic (T), pre-stress (Ug) and elastic strain (Ue) energy components: Where the energy terms could be defined as, .
e e e e g g g g With denoting of the material mass density as , differentiation with respect to time as a upper dot and a matrix transpose operator as superscript T .The force resultants (N, ,O,P,Q,R,T,U) can be related to the strain terms via the curvilinear fiber laminated material equivalent stiffness matrices.i.e.: (2) Substituting the strain and force resultant equations in energy integrals (equation 7), minimizing the energy equilibrium equation 6, factorizing with respect to the degrees of freedom vectors, and some further handlings including assembling the strip equations and implementing of necessary boundary conditions, a Mathieu-type differential governing equation is obtained as, Latin American Journal of Solids and Structures 14 (2017) 528-543 Where M, K, Kg S and Kg D are the global structural matrices corresponding respectively to mass, elastic strain, static stress and dynamic initial stress energies.d is the global vector of uncon- strained degrees of freedom.By implementing the Bolotin's first order approximation corresponding to the period twice the loading period, which is more critical (Bolotin, 1964), the time varying vector, d , is approximated as: A and B are time-independent coefficient vectors called degrees of freedom vectors.Substitution of equation 10 into equation 9, factorization of harmonic terms and setting their coefficients to zero leads to a set of homogenous equations.For a non-trivial solution of unknown vectors A and B, the determinants of the coefficient matrices should be set to zero.The governing equations are reduced to two subsequent eigenvalue problems as, The equations corresponding to vectors A and B are not coupled and can be solved separately which reveals the two boundaries of the instability region for the structure in terms of loading parameter sets of ( , , )   S D a a w .

Curvilinear Fiber Simulation Considerations
Using machine fiber placement technologies, it is possible to change the fiber orientation through the geometry.The resulting production is called as a curvilinear fiber placement.In this research it is assumed that the fiber angle changes linearly just in the longitudinal direction of the geometry.As a result, the lamination layup in the panel width is constant while changing in the longitudinal direction.The changing fiber angle is denoted by a two-angle set <T0,T1> where the former one and the latter are the fiber angle at the plate middle section and plate longitudinal ends, respectively.In other words, the fiber angle in a single ply is symmetric with respect to the plate middle section.So, the fiber angle at every arbitrary point in the geometry is given by the following equation, Diagram of Figure 2 depicts the changing fiber orientation in a sample strip coordinate frame.Due to the changing nature of lay-up properties in the longitudinal direction, the calculations for laminate equivalent stiffness matrices is performed in each integration point, independently.

Delamination Modeling Technique
In the delamination region, the plate is actually a set of two thinner plates.To bring a single delamination effects into consideration, the main idea is to use double strips in the thickness direction.This means that the whole plate is composed of two similar layers of strip meshes with special layups.Inside the delamination zone, the two layers are independent and has no connections (note that the probable contact between two adjacent layers is ignored at present.).Out of delamination zone, all of degrees of freedom of the two layer must rigidly linked to each other via knots' merging process.The corresponding strips in upper and lower layers have the same geometrical and numerical characteristics.According to Figure 3, the strip knots of the same planar positions are merged together in all the perfect panel areas and also in the delamination zone edges.It is important to have knots on the edges of the delamination region.This approach could also be generalized for the case geometry with N delamination in thickness direction (N+1 layers is needed to be defined).Every strip layer has the same geometry properties but are different in bending stiffness.To fulfill the true bending properties of every layer in a strip with respect to the plate mid.-surface, every layer lay-up is considered similar to the whole plate layup with the redundant layers' material changed to a null stiff-less and weight-less one.(see Figure 3) These considerations provide the physical conditions at the edges of the delamination zone.

RESULTS AND DISCUSSION
The first two natural frequencies and fundamental buckling strength of a laminated composite rectangular plate are investigated.The lamina material properties and the model geometry are characterized as, 134.4 , / 13, 0.25 0.127 , / 10, / 125 , / 0;0.2;0.4;0.6 / 0.5, 1480 / The laminate has 8 layers with fixed orientation lay-up [0/90/0/90]s.A central through-thewidth delamination in the mid layer is considered as shown in Figure 4. Table 1 shows the results for the first two natural frequencies as well as fundamental critical buckling load.The results from higher order deformation theory, layerwise theory and experimental tests are also provided from the literature.A good agreement could be seen between the present FSM calculations and the reference ones.The first eight natural frequencies of three-layer square VSCL plates of diverse thickness, for simply-supported and clamped boundary conditions are extracted and presented in comparison to higher order FEM results of Akhavan and Ribeiro (2011) The results show the very good consistency of the FSM higher order as well as classical formulations.
Square three layer VSCL plate of lay-up [<0,45>//<-45,-60>/<0,45>] with an embedded central square delamination between first and second layers is assumed.The geometrical properties are the same as presented in equation ( 14) and Figure 4.The length to thickness of the plate (L/t) is 10.A uniformly distributed longitudinal stress is assumed everywhere on the plate.Effect of different out of plane edge constraint are studied.Seven different constraint sets including CCCC, SSSS, CSCS, CFCF, CCCF, SFSF, SSSF are concerned where C, S and F stand for clamped, simply supported and free conditions starting from a longitudinal end, respectively.No inplane movements are permitted at plate edges.Table 5 and 6 represent the results of the first four critical buckling loads as well as natural frequencies of perfect and delaminated geometries (corresponding to delamination areas of 0%, 4%, 16% and 36%), respectively.The results show that by growing the delamination area, the VSCL plate critical fundamental natural frequencies meet extreme reductions such that in the case of plate with 36% delamination area experiences a 14% fundamental frequency drop with respect to the perfect plate in CCCC conditions.The table also indicates that the higher the models are constrained, the more significant reduction in the natural frequencies will take place.The reductions fall in the interval of 4 to 14 percent in cases under study (see Figure 5).The delamination in VSCL plate leaves more notable effects on reduction of static buckling strength of the panel as shown in Figure 5.In this contest also the more constraint leads to higher softening of the structure due to probable delamination.Figure 7 depicts the change in dynamic instability regions map of the perfect and delaminated VSCL panel with respect to different loading amplitude and frequencies.A pure longitudinal harmonic loading (a S =0.0) is assumed.The results show that the burst of delamination shifts the base instability frequencies (a D =0.0) toward slightly more critical lower ones while shrinks the instability region.A 4% delamination shows very ignorable results while higher delaminations make more significant changes.Clamped VSCL plate containing a central 36% delamination region is considered with the reference lay-up [<0,0>//<90,90>/<0,0>].The central fiber angles in all layers (T0) is changed and the instability regions under inplane longitudinal uniform loading is extracted.Figure 8 depicts the instability region limits of various lay-ups.A pure dynamic loading is assumed.The results show that with change in central lay-up, the instability frequencies shift initially to the higher ones and then to the lower ones.A lay-up with fiber angles of [<15,0>//<-75,90>/<15,0>] benefits with higher instability frequencies that may interpreted to higher stable lay-up of the panel.It is also shown that in case of unchanged central fiber angles, the lay-up [<0,15>//<90,-75>/<0,15>] presented higher instability frequencies.The study demonstrates a slightly higher effect of changing the central fiber angle in comparison with change in end fiber angles.Clamped VSCL plate with and without central 36% delamination region is considered with layup [<0,45>/<-45,-60>/<0,45>]. A longitudinal inplane loading with static preload coefficient (a S =N S /Ncr) of 0.0 to 1.0 is considered and the dynamic instability regions of the VSCL plate are derived.The instability regions are shown in Figure 9.According to the figure, with growing the static loading coefficient, the instability region shifts toward lower loading frequencies while its size increases.In case where the static preload equals to the plate's buckling critical load (a S =1.0), there is already stable plate conditions for sufficiently high loading frequencies.The results also show that the difference between perfect and delaminated plate is reduced for higher static preloads.Moreover, the instability regions of delaminated plate shrinks with respect to the ones corresponding to the perfect geometry.

CONCLUSIONS
The static as well as dynamic stability behavior of moderately thick variable stiffness composite laminated (VSCL) plates containing inter-ply square delamination region is investigated by using of a higher order B-spline finite strip formulation.A harmonically time-varying in-plane longitudinal loading is assumed.The friction and contact effects at delaminated interfaces are neglected.The governing equations are derived using full energy concepts on the basis of the principle of virtual work.The dynamic parametric instability load-frequency margins are extracted utilizing the Bolotin's first order approximation followed by an eigenvalue analysis.Comparisons imply that the spline formulation is a reliable tool in calculation of delamination effects as well as variable stiffness laminated plates stability problems.To the best of the author's knowledge, this is the first application of B-spline FSM to the problem.Some representative results with change in lay-up, boundary conditions and delamination size are provided.The results show that the growth of the delamination region made more significant destabilizing effects for the most constrained cases.It is also shown that proper fiber angle design may lead to minimize the delamination effects on the stability of the plate.The instability of the panel shifts to occur at lower loading frequencies as a delamination growing but the severity of changes depends on the type of out-of-plane boundary conditions.

Figure 1 :
Figure 1: Typical delaminated plate finite strip mesh, finite strip geometry and loading.

Figure 2 :
Figure 2: Curvilinear variable fiber orientation in a single lamina.

Figure 4 :
Figure 4: The geometry and delamination position of the cantilever plate (left) and square plate (right).

Figure 6
Figure 6 depicts the fundamental free vibration mode shape in different boundary condition sets for the case of three-ply VSCL plate with central 36% delamination area.

Figure 5 :Figure 6 :
Figure 5: Change in static and dynamic instabilities due to delamination area growth for different constraints.

Figure 7 :
Figure 7: Dynamic instability region for VSCL plate with variable delamination areas under different boundary conditions.

Table 1 :
Buckling critical load and natural frequencies of thin clamped beam-plate with delamination.

Table 7 :
Natural frequencies of square three-ply VSCL plate with central embedded delamination (Hz).