Buckling analysis of laminated sandwich beam with soft core

Stability analysis of laminated soft core sandwich beam has been studied by a C 0 FE model developed by the authors based on higher order zigzag theory (HOZT). The in-plane displacement variation is considered to be cubic for the face sheets and the core, while transverse displacement is quadratic within the core and constant in the faces beyond the core. The proposed model satis(cid:12)es the condition of stress continuity at the layer interfaces and the zero stress condition at the top and bottom of the beam for transverse shear. Numerical examples are presented to illustrate the accuracy of the present model.


INTRODUCTION
Laminated composite structures are widely used in aerospace, naval, civil, and mechanical industries due to its superior properties such as high strength/stiffness to weight ratio and greater resistance to environmental degradation over the conventional metallic materials.The use of high strength materials in such construction leads to slender sections, thus making buckling as a primary mode of failure of the member subjected to axial compressive forces.These buckling modes depend upon the material properties and relative stiffness of the face sheets and core.
The most important feature of the laminated composites (e.g., GFRP, CFRP etc.) is that these are weak in shear due to their low shear modulus compared to extensional rigidity.In order to fulfill the weight minimization, a laminated sandwich construction having low strength core and high strength face sheets may be used.As the material property variation is very large between the core and face layers in case of sandwich construction, the effect of shear deformation becomes more complex.
Based on assumed displacement field, the theories proposed for the accurate analysis of composite laminated structures are grouped as Single layer theory and Layer-wise theory.In single layer theory the deformation of the plate is expressed in terms of unknowns at the reference plane, which is usually taken at the middle plane of the plate.This theory is also known as first-order shear deformation theory (FSDT).Based on the first-order shear deformation theory, Allen [5] presented a three-layered model for the analysis of sandwich beams and plates wherein the zigzag deformation pattern was considered.B-spline functions based on FSDT have been used by the Dawe and Wang [18,19] as trial functions for Rayleigh-Ritz analyses and finite strip analyses.Aiello and Ombres [2] developed a model using FSDT to assess the optimal arrangement of hybrid laminated faces of sandwich panel for local buckling loads.Wang [49] employed the B-Spline Rayleigh-Ritz method based on first order shear deformation theory to study the buckling problem of skew composite laminates.Sundersan et al. [46] have studied the influence of partial edge compression of composite plates by using a finite element method based on FSDT for analyzing isotropic plates.The structural behavior of sandwich laminated composites cannot be assessed satisfactorily by a FSDT, as the core and face sheets deform in different ways due to wide variation of their material properties, which is identified by a kink in the variation of in-plane displacements across the thickness at the interface between the core and stiff face layers.So as to predict the actual parabolic variation of shear stress, this (FSDT) required shear correction factor.This has inspired many researchers to develop a number of refined plate theories.
A higher order variation of in-plane displacement through the thickness is considered to represent the actual warping of the plate cross-section for overcoming the need of shear correction factor so it is called as higher order shear deformation theory (HSDT).The warping of the cross section also allows higher order variation of transverse shear stresses/strains across the depth [38].By choosing appropriate shape functions, four shear deformation theories presented by Aydogdu [10], are named as: parabolic shear deformation beam theory, hyperbolic shear deformation beam theory, first order shear deformation beam theory and exponential shear deformation beam theory, for the buckling analysis of composite beam.The theories developed for buckling analysis by Reddy and Phan [39], Kant and Kommineni [26], Khedir and Reddy [28] , Moita et.al [33,34], Vuksanovic [48] and many other falls under this category.Song and Waas [44] presented a higher order theory for the buckling and vibration analysis of composite beams and the accuracy of HSDT was demonstrated compared to 1-D Euler-Bernoulli, 2-D classical elasticity theory and Timoshenko beam theory.Frostig [23] obtained the general as well as local buckling loads for sandwich panels consisting of two faces and a soft orthotropic core using HSDT.Karama et.al [27] presented a multilayered laminated composite model by introducing a sine function as a transverse shear stress function where a FE based software package Abaqus was used to check the efficiency of the model.Babu and Kant [11] presented two C 0 isoparametric finite element formulations, one based on FSDT and other based on HSDT to investigate the effect of skew angle on buckling coefficient.Matsunaga [29,30] developed one dimensional global higher order theory, in which the fundamental equations were derived based on the power series expansions of continuous displacement components to analyze the vibration and buckling problems.A new triangular element was developed by Chakrabarti and Sheikh [14,15], for the buckling analysis of composite plate using Reddy's higher order theory.The effect of partial edge compression was investigated [15].Nayak et.al [4] developed a shear deformable plate bending element based on a third order shear deformable theory (i.e.HSDT).
Zhen and Wanji [51] presented global theories with higher order shear deformation and zigzag theories satisfying continuity of transverse shear stresses at interfaces for analyzing the global response of sandwich laminated beams.The effects of the number of higher order terms in the shear deformation as well as inter-laminar continuity of shear stress on global response of laminated beams and soft core sandwich beams were studied in this paper [51].
The actual behavior of a composite laminate is that the transverse strain may be discontinuous but the corresponding shear stress must be continuous at the layer interface [41], while HSDT shows discontinuity in the transverse shear stress distribution and continuity in the variations of the corresponding transverse shear strain across the thickness.This is mainly due to the different values of shear rigidity at the adjacent layers.
The above disparity leads to the development of layer-wise theories, which started with discrete layer theories.In discrete layer theories the unknown displacement components are taken at all the layer interfaces including top and bottom surfaces of the structure.Discrete layer theories proposed by Srinivas [45], Toledeno and Murakami [47], Robbins and Reddy [25] , Cetkovic [13] and many others assume unique displacement field in each layer and displacement continuity across the layers.In these theories, the number of unknowns increases directly with the increase in the number of layers due to which it required huge computational involvement.
The improvement of layer-wise theories comes in the form of zigzag theories (ZZT), where the unknowns at different interfaces are defined in terms of those at the reference plane.The number of unknowns are made independent of the number of layers by introducing the transverse shear stress continuity condition at the layer interfaces of the laminate and the in plane displacements have piece-wise variation across the plate thickness in this theory (ZZT).Averill [6] developed a C 0 finite element based on first order zigzag theory and overcome the C 1 continuity requirement by incorporating the concepts of independent interpolations and penalty functions.Hermitian functions were used by Di Sciuva [20,21] to approximate the transverse displacement in the formulations.Carrera [22] used two different fields along the laminate thickness direction for displacement and transverse shear stress respectively for his formulation.Averill and Yip [7] developed a C 0 finite element based on cubic zigzag theory, using interdependent interpolations for transverse displacement and rotations and penalty function concepts.
In some improved version of these theories, the condition of zero transverse shear stresses at the plate/beam top and bottom was also satisfied.Kapuria et.al [40] developed 1D zigzag theory for the analysis of simply supported beams.The effect of the laminate layup and the thickness to span ratio was investigated.Zigzag models for laminated composite beams were developed by using trigonometric terms to represent the linear displacement field, transverse shear strains and stresses [42,43].However, the zigzag theory has a problem in its finite element implementation as it requires C 1 continuity of the transverse displacement at the nodes.
To combine the benefits of the discrete layer wise and higher order zigzag theories, Cheung et.al [50], Yip and Averill [24], Icardi [8,9] and many other authors developed theories which are known as sub-laminated models.Aitharaju and Averill [3] developed a new C 0 FE based on a quadratic zigzag layer-wise theory.For eliminating shear locking phenomenon, the shear strain field is also made field consistent.The transverse normal stress was assumed to be constant through the thickness of the laminate.The new FE was applied to model the beam as combination of different sub-laminates.
Ramtekekar et.al [36], Bambole and Desai [12] and many more developed a mixed FE approach for accurate analysis of stresses, where the stress components are assumed as the field variables at interface nodes along with displacement field variables.Dafedar et.al [17] and Dafedar and Desai [16] have presented a theory based on mixed higher order theory for the buckling analysis of laminated composite structures.Two sets of mixed models HYF1 and HYF2 have been presented by selectively incorporating non-linear components of Green's strain tensor.Individual layer theory (ILT) based HYF1 and HYF2 models have been developed by considering a local Cartesian co-ordinate system at the mid-surface of each individual layer and at mid-surface of the entire structure respectively.
The transverse deformation is very significant in case of a laminated sandwich structure having a soft core as there is abrupt change in the values of transverse shear rigidity and thickness of face sheet and the core.As such to achieve sufficient accuracy, unknown transverse displacement fields across the depth in addition to that in the reference plane are essential to represent the variation of transverse deflection.This can be done by using sub-laminate plate theories but the number of unknowns will increase with the increase in the number of sub-laminates.On the other hand, introduction of additional unknowns in the transverse displacement fields invites additional C 1 continuity requirements in its finite element implementation by using the zigzag theory as mentioned earlier.However, the application of a C 1 continuous finite element is not encouraged in a practical analysis.
Recently Pandit et.al [31,32] proposed a higher order zigzag theory for the static and buckling analysis of sandwich plates with soft compressible core.In this model [31,32] a ninenode isoparametric element with 11 field variables per node was employed.To overcome the problem of C 1 continuity the authors [31,32] have used separate shape functions to define the derivatives of transverse displacements such C 0 finite element was used for the implementation of the theory.It has imposed some constrains, which are enforced variationally through penalty approach.However, the problem of choosing suitable value for the penalty stiffness multiplier is a well known problem in the finite element method.
In this paper an attempt is made to do the stability analysis of laminated sandwich beam having soft compressible core based on higher order zigzag theory by using a C 0 beam finite element model recently developed by the authors [1].The model satisfies the transverse shear stress continuity conditions at the layer interfaces and the conditions of zero transverse shear stress at the top and bottom of the beam.The isoparametric quadratic beam element has three nodes with seven displacement field variables (i.e., in-plane displacements and transverse displacements at the reference mid surface, at the top and at the bottom of the beam along with rotational field variable only at the reference mid surface) at each node.The displacement fields are chosen in such a manner that there is no need to impose any penalty stiffness in the formulation.The element may also be matched quite conveniently with other C 0 elements.The present beam model is used to solve stability problems of laminated sandwich beams having different geometry boundary conditions.

MATHEMATICAL FORMULATIONS
The in-plane displacement field (Figure 1) is chosen as follows: where,u 0 denotes the in-plane displacement of any point on the mid surface, θ x is the rotation of the normal to the middle plane about the z -axis,n u and n l are number of upper and lower layers respectively, β x and θ x are the higher order unknown, α i xu and α i xl are the slopes of i-th layer corresponding to upper and lower layers respectively and H (z − z u i ) and H (−z + z l i ) are the unit step functions.The transverse displacement is assumed to vary quadratically through the core thickness and constant over the face sheets (as shown in Fig. 2)and it may be expressed as, W = l 1 w u + l 2 w 0 + l 3 w l for core = w u for upper face layers = w u for lower face layers where w u , w 0 and w l are the values of the transverse displacement at the top layer, middle layer and bottom layer of the core, respectively, and l 1 , l 2 and l 3 are Lagrangian interpolation functions in the thickness co-ordinate.The stress-strain relationship based on a plane stress condition of an orthotropic layer/ lamina (say k -th layer) having any fiber orientation with respect to structural axes system (x -z ) may be expressed as where {σ}, {ε} and [ QK ] are the stress vector, the strain vector and the transformed rigidity matrix of k -th lamina, respectively.Utilizing the conditions of zero transverse shear stress at the top and bottom surfaces of the beam and imposing the conditions of the transverse shear stress continuity at the interfaces between the layers along with the conditions, u = u u at the top and u =u l at the bottom of the beam, β x ,η x ,α i xu ,α i xl , (∂w u /∂ x )and (∂w l /∂ x ) may be expressed in terms of the displacement u 0, θ x , u u and u l as where, T and the elements of [A] are dependent on material properties.It is to be noted that last two entries of the vector {B} helps to define the derivatives of transverse displacement at the top and bottom faces of the beam in terms of the displacements u 0 , θ x , u u , and u l to overcome the problem of C 1 continuity as mentioned before.
Using the above equations, the in-plane displacement field as given in equation ( 1) may be expressed as where, the coefficients b ′ i s are function of thickness coordinates, unit step functions and material properties.
The generalized displacement vector {δ} for the present beam model can now be written with the help of equations ( 2) and ( 5 Using strain-displacement relation and equations ( 1)-( 4), the strain field may be expressed in terms of unknowns (for the structural deformation) as where, The linear strain part is and non-linear strain part can be written as where, and the elements of [H ], [A G ] and [H G ] are functions of z and unit step functions.The values of l 1 , l 2 , l 3 , b i 's and the elements of [H ] are used as reported in the appendix (A,B and C) by Chakrabarti et al. [1].With the quantities found in the above equations, the total potential energy of the system under the action of transverse load may be expressed as where U s is the strain energy and U ext is the energy due to the externally applied in-plane load.
Using equations ( 3) and ( 6) , the strain energy (U s ) is given by where, and the energy due to application of external in-plane load can be calculated as where, Latin American Journal of Solids and Structures 9(2012) 367 -381 and the stress matrix [S i ] may be expressed in terms of in-plane stress components of the i-th layer as In the present problem, a three-node quadratic element with seven field variables (u 0 , w 0 , θ x , u u , w u , u l and w l ) per node is employed.Figure 3 shows the node numbering and natural coordinate of the element.Using finite element method the generalized displacement vector δ at any point may be expressed as where, δ = u 0 , w 0 , θ x , u u , w u , u l , w l T as defined earlier, δ i is the displacement vector corresponding to node i, N i is the shape function associated with the node i and n is the number of nodes per element, which is three in the present study.
With the help of equation ( 12), the strain vector {ε} that appeared in equation ( 6) may be expressed in terms of field variables as where is the strain-displacement matrix in the Cartesian coordinate system.The elemental potential energy as given in equation ( 7) may be rewritten with the help of equations ( 8)- (13) as where, where, [K e ] and [K G ] are stiffness matrix and geometrical stiffness matrix.
Latin American Journal of Solids and Structures 9(2012) 367 -381 The equilibrium equation can be obtained by minimizing Π e as given in equation ( 14) with respect to {δ} as where λ is a buckling load factor.The displacements, stresses (in case of static analysis) and buckling loads in laminated sandwich beams are calculated by developing a numerical code to implement the above mentioned operations involved in the proposed FE model.The skyline technique has been used to store the global stiffness matrix in a single array and Gaussian decomposition scheme is adopted for the static solution, while simultaneous iteration technique is used for solving the buckling equation (17).

NUMERICAL RESULTS
In this section, some examples on the buckling analysis of laminated sandwich beam having a soft core is analyzed by using the proposed FE model based on higher order zigzag theory.The accuracy and applicability of the proposed FE model based are demonstrated by solving a number of problems having different features.The results obtained are compared with the published results in some cases and finally many new results are also generated.The following different boundary conditions are used 1.Simply supported boundary condition (H): The field variables u 0 , w 0 , w u and w l are restrained while θ x , u u and u l are unrestrained in one boundary.In the other boundary, the field variables w 0 , w u and w l are restrained while u 0 , θ x , u u and u l are unrestrained.

Clamped boundary condition (C):
All the nodal field variables at the boundary are fully restrained.
3. Free boundary condition (F): All the nodal field variables at the boundary are unrestrained.
The non-dimensional quantities used, to show different results are as follows where h and l are depth and span length of the beam respectively and E T f is the transverse modulus of elasticity of face layer.

Three layered laminated composite beam
For the convergence studies of critical buckling load, a three layer beam is analyzed taking mesh divisions 2, 5, 10, 16, 20 and 30.The layers are of equal thickness.The material properties [28] used in this example is as shown in Table 1 (Example 1).The beam is analyzed by considering four different boundary conditions and different values of thickness ratio (l/h).The results obtained by using the proposed FE model are presented in Table 2.It may be observed in Table 2 that the buckling loads are converged at mesh division 20.A mesh division 20 is taken for all subsequent analysis to get accurate results.The present results corresponding to arbitrary boundary conditions are compared with those published; FSDT [10] and TSDT [28].It may be observed in Table 2 that the present results are in good agreement with the other results.

Multilayered composite beam
In this example the effect of number of layers is investigated by analyzing a simply supported laminated composite beam.Different layups are considered for the analysis by increasing the layer numbers from 3 to 10 i.e. symmetric layup (2,5,7,9) and non-symmetric layup (4,6,8,10) to investigate the effect of bending-extension coupling.The thickness of each layer is assumed to be equal.The properties of materials [30] used in this problem are shown in the Latin American Journal of Solids and Structures 9(2012) 367 -381 Table 1(Example 2).The critical buckling load has been calculated by considering different thickness ratios and varying the number of layers.The effect of bending-extension coupling may be seen from Table 3, which reduces the buckling load.The results are shown in Table 3, along with the results reported in the literature based on global local higher order theory (GLHT) and ZZT [51] and GLHT [30].It may be observed from Table 3 that the present results based on the higher order zigzag theory for the critical buckling load are on the lower side as compared to the other results based on ZZT [51] and GLHT [30], while it shows a good agreement with the results based on GLHT [51] .

Three layer soft core sandwich beam
A three layer simply supported sandwich beam is analyzed in this example using the present FE model.The ratio of the thickness of core to thickness of face sheets (t c /t f ) is considered to be 5, 25 and 50.The material properties [16] for the face sheets and core are given in Table 1 (Example 3).In Table 4, the non dimensional critical buckling load is reported for different thickness ratio (l/h).Results reported by Dafedar and Desai [16] based on higher order mixed formulation and the results reported by Zhen and Wanji [51] based on GLHT and ZZT along with the results obtained by using Abaqus (Ver.6.8)software package are considered for the comparison.
To model the laminated sandwich beam in Abaqus, 8000 eight node 3D shell elements (S8R) are used.From Table 5, it can be concluded that the present model performance is quite satisfactory as compared those presented by Zhen and Wanji [51] and Dafedar and Desai [16].The present 1D results based on the refined theory are expected to be closer to the 3D analysis results.This is observed in Table 5 as the present results are quite close to the 3D analysis results obtained by using Abaqus.

Multilayered laminated soft-core sandwich beam
A multilayered laminated sandwich beam (0 ○ /90 ○ /0 ○ /90 ○ /C/90 ○ /0 ○ /90 ○ /0 ○ ) is analyzed in this example which is a new problem.Material properties [37] are as shown in Table 1 (Example 4).The values of non-dimensional buckling load are presented in Table 5 for different thickness   The ratio of thickness of core to thickness of the face sheets is considered to be 5 and 50 for the plot.The plot shows expected variation.

CONCLUSION
The stability behavior of laminated sandwich beam having a soft core is studied in this paper by using a C 0 efficient finite element (FE) model developed by the authors.A refined higher order zigzag shear deformation theory is used to model the in-plane displacement field and a quadratic displacement field is used to define the transverse displacement in the proposed FE model.The nodal field variables are chosen in an efficient manner to overcome the problem of continuity requirement of the derivatives of transverse displacements; and there is no need to use penalty functions in the formulation as used by many previous researchers.problems.The results obtained by using the present FE model are successfully compared with many published results and also with the results obtained by using FE based software package Abaqus.Many new results are also presented which should be useful for future research.

Figure 1
Figure 1 General lamination scheme and displacement configuration.

Figure 2
Figure 2 Variation of transverse displacement (w ) through the thickness of laminated sandwich beam.
) as {δ} = {u 0 w 0 θ x u u w u u l w l } T Latin American Journal of Solids and Structures 9(2012) 367 -381

Figure 3
Figure 3 Node numbering and natural co-ordinate of the beam element.

t c /t f l /
critical buckling load for different boundary conditions mainly; Hinge-Hinge (H − H), Clamped-Hinge (C − H), Clamped-Clamped (C − C) and Clamped-Free (C − F ) over different thickness ratio (l/h) is shown in the Figure 4 (a − b).

Figure 4
Figure 4 Variation of critical buckling load.

Table 1
Material properties for laminated beams

Table 3
Comparison of buckling stresses of multi-layered laminated beams.

Table 4
Comparison of critical loads for three layered soft-core sandwich beam.ranging from 5 to 100 by considering the ratio of the thickness of core to thickness of face sheets as 25 for a simply supported boundary condition.To check the accuracy of present FE model to this new problem of multilayered sandwich beam, present results are compared with those obtained by using finite element software package Abaqus (Ver.6.8).

Table 5
Comparison of critical loads for multilayered soft-core sandwich beam.