Buckling analysis of laminated composite elliptical shells using the spline finite strip procedure

According to the first order shear deformation plate assumption, current paper numerically investigates the buckling problem of cylindrical shells with oval cross-section under simply supported boundary condition subjecting to uniform and non-uniform loads via spline finite strip method. Equilibrium governing equations are generated based upon the virtual work's principle. Shell displacement function is assumed to be as dual development of degree 3 spline functions and Lagrange polynomial functions. The analysis method used at any time with any type of loading can calculate deformation and buckling. The step-by-step critical load is calculated using the determinant of the stiffness matrix. In the present paper, it is supposed that the applied linear elastic materials consist of isotropic and non-isotropic materials. Numerical studies are performed to indicate the effects of shell geometry, materials, and layered materials on the buckling load. For the validation purpose, the results of current research are compared with those of previous researches and ABAQUS software.


Introduction
Elastic instability (buckling) phenomena, caused by applying the axial compressive forces and shearing loads to the structures, is one of the important mechanical characteristics that has always been interested in the design of engineering structures (e.g., beams, plates, tanks, columns, shells) by many researchers [1][2][3][4][5][6][7].Elastic instability (buckling) phenomena could be a shape of flimsiness that happens in a flexible framework, such as buckling of pillars and pieces subjected to compressive loads.Due to the high application demand, shell-type structures, particularly, cylindrical shells having superior mechanical characteristics and excellent high strength performance are applied widely in many engineering problems such as aerospace [8], nuclear reactors [9], the marine industry [10], civil engineering [11], and mechanical engineering [12].So far, numerous researches investigated the buckling behavior of shell-type structures.For example [13], experimentally and numerically reported the buckling problem of cylindrical shells with longitudinal joint subjected to the lateral external pressure.On the basis of extended Kantorovich method [14], predicated buckling and post-buckling behavior of cylindrical shells under localized external pressure [15].Numerically studied the effect of periodical imperfections on the changes of the buckling load of closed elastic isotropic shallow conical shells [16].Used an effective technique for analytically achieving the buckling load of a cylindrical shell under nonuniform external pressure [17].Carried out the first-order shear deformation in conjunction with Rayleigh-Ritz technique to find the buckling load of the graphene-reinforced porous cylindrical shell.
In recent years, the composite material has been widely used for their superior performance in reinforcing materials and improving specific stiffness, specific strength, resisting and environmental forcing conditions in structural parts.The development of applications using composite materials is because of their desirable features including great specific strength to stiffness ratio, long fatigue life, and magnetic transparency.It should be noted that the realization the relation between fiber orientation, advanced damage, and the corresponding deformation areas is necessary when utilizing such materials in safety-critical applications.The composite cylindrical shells are abundantly employed in the various engineering fields such as aerospace, marine, automobile.So far, many research projects have been carried out focusing on the investigation of the application of composite materials in engineering structures like composite cylindrical and conical shells, [18][19][20][21][22][23][24][25].Furthermore, recently, buckling problem of composite cylindrical shells are reported by Refs.[25][26][27][28][29][30][31][32].Fathollahi et al. [33] examined the Investigation of the impact of radiation on 2D MHD Al 2 O 3 unsteady water stream through parallel compression plate AGM and HPM.The innovation of this study, the dimensionless values of different fluids, such as Prandtel number and coefficient of friction, influence the study of velocity and heat transfer parameters of oxide nanofluids.Shadman et al. [34] studied the combined fins of baffle and chamfer on threatened tension surfaces under the influence of nanofluids and magnetic parameters for rotary joints in computer hardware.The innovation of this paper is to study the parameters of Mwcnt and TiO2 nanofluids transmitted from different fins and blades on the elongated surface and compare their results with each other.Pasha et al. [35] consider the application of numerical techniques to micropolar fluid flow and heat transfer in permeable plates.The novelty of this research is to study heat transfer and fluid flow through two equal sheets with Adomian Decomposition Method (ADM) and Variable Iterative Method (VIM).Low-speed non-linear impact response of a cylindrical porous metal shell reinforced with graphene platelets under axial motion with geometrical imperfections reviewed by Zhang et al. [36,37].Nonlinear primary resonance analysis of metallic foams reinforced with initial load-bearing graphene platelets, double-curved shells with imperfect geometries examined by Gui-Lin et al. [38].Katiyar et al. [39] investigated the vibrational response of a functionally classified two-dimensional geometrically discontinuous surface lying on an elastic foundation in a thermal environment with initial imperfections.The preceding review explicitly denoted that there are an increasing number of investigations on the buckling problem of thin-walled composite cylindrical shells with circular cross-section [40][41][42][43].Composite materials are created by combining materials into an overall structure that has distinctive properties than the person components.Composite structures are used in a variety of industries from aerospace, marine, aviation, transportation, sports/recreation to construction [44][45][46][47].However, there seems to be no systematic, analytical or numerical studies on utilizing the FSDT and spline finite strip method (SFSM) to investigate the buckling problem of the laminated composite elliptical cylindrical shell under simply supported boundary condition subjected to the uniform and non-uniform axial loadings.The displacement components are defined as the multiplication expansion of the 3rd degree spline functions in the longitudinal direction and the Lagrange interpolation functions in the crosswise direction.The governing buckling equations are inferred using the vibrational rule and solved with consideration of Sanders nonlinear theory.Analysis was conducted in two steps; in the first step, the shell subjecting to the load and displacements are calculated by static analysis and then the geometric hardness matrix is calculated with respect to Sanders nonlinear theory and FSDT.In the second step, combining the linear and geometric hardness matrices, the critical shell load is calculated by analyzing the special value.At the end, the results inferred from the current analysis are compared with those of previous studies and the results of limited version of the ABAQUS software.The purpose of writing this article and its application in this way is a systematic numerical and analytical study with two new methods, FSDT and SFSM, to investigate the buckling of composite cylindrical shell under a special condition, which is mentioned for the first time in this article.The innovation of this article includes 2 parts related to each other.In the first stage of the new ideation, all displacements and the shell under load were checked with static analysis and the geometric matrix was calculated according to the FSDT theory.In the second stage of innovation in this article, the critical load is calculated by eigenvalue analysis of geometric matrices.

Modeling of the problem and theoretical formulation
As schematically demonstrated in Fig. 1, an elastic composite laminate elliptic cylinder is composed of N number of laminas, where the position of each point of the cylindrical shell in the local coordinate systems is indicated with {x,θ,z}.It is also assumed θε[0, 2π] and the lengths of semi-major and semi-minor axes of the elliptic meridian to be, respectively, A and B. To consider the effects of deformation at each shell cross section, the terms u, v, w are supposed to represent, respectively, longitudinal, tangential, and radial displacements.
It must be noted that for each point of the elliptical cross section, the mean radius r(α), curvature χ(α), and central angel α can be defined as follows [28]:

Introduction of SFSM
The spline function is a mathematical polynomial function named in terms of the degree of function.Spline function can be expressed in the form of B1-B2-B3-B4 spline and B5 spline.Here, cubic B3 spline function is applied to demonstrate the displace-N.Korkeai et al. ments.The B3 spline function (grade 3) is one of the best choices to solve the shell problems, because this function makes it easy to link these problems.Fig. 2 shows a vertical view of a striped element of length L and width b.On the other hand, Fig. 3 reveals how the nodal lines are changed along a strip element.The length of the sheet is divided into n sections (q represents the length of each section), and these divisions are from the points φ − 1 to φ q+1 , with the point φ 0 as the first support and the point φ q as the end support of the strip.The two points outside of the longitudinal range of the strip and on the nodal line are used to apply boundary conditions and to modify the spline functions of the boundary areas.In the longitudinal direction, on each nodal line, q+1 nodes numbering from node − 1 to node q+1 are used.
In the crosswise direction of the strip, the Lagrangian functions of degree two with three nodes are used as indicated in Fig. 4, which can be presented as [28]: where L j (η) is the term of the shape function relating to the j th nodal line of the strip.Also, η refers to the transverse curvilinear coordinate.However, the general displacement function can be estimated in components of the curvilinear coordinate system using B3spline and Lagrangian functions as follows [28]: in which F expresses the generic generalized displacement function and φ denotes the general term of the B3-spline functions.
According to the FSDT [29], the displacement field (u, v, w) for an optional point can be written as [28]: where u, v and w refer to the middle surface displacements of shell domain along x, y, and z directions, respectively.Also, ψ x and ψ θ refer to the normal transverse rotations at the edge of the sheet around the x and θ axes respectively.On the basis of Sanders strain-Fig.2. A strip with its knots and nodal lines.displacement relations, the linear and nonlinear components of strain-displacement can be defined as [29]: It should be noted that in Eq. ( 5) by setting C 1 = 0, C 2 = C 3 = 1, Sander's strain-displacement relationships can be obtained.By taking into account the plane stress assumption, the constitutive stress-strain relationships for a single lamina can be denoted as [29]: ) refer to the transformed reduced stress -stiffness values and Q ij (i, j = 4, 5) represent the shear-stiffness terms [30].By integrating stresses through the laminate thickness, the stress resultants can be obtained as [30]: where N x , N θ , N xθ , Q x and Q θ , respectively, are the normal membrane forces, shear membrane forces and shear forces along the thickness per unit length.Also, M x , M θ and M xθ represent bending moments and twisting moment per unit length.Furthermore, A ij , B ij and D ij refer to the membrane rigidity, bending-membrane rigidity and bending rigidity, which can be denoted as [30]:  ( in which k i and k j are shear correction coefficients.

Virtual work
The principle of virtual work is employed to obtain equilibrium equations.On the basis of the principle of virtual work, the resultant virtual work exerted by internal and external forces for virtual displacements should be zero as [30]: where δw int and δw ext express, respectively, the internal and external virtual works.According to the stress vector σ and virtual strain vector δε, internal virtual work can be presented as By integration thorough thickness of the strip in Eq. ( 10), the internal virtual work can be rewritten in terms of generalized stresses strains as By substituting linear and nonlinear strains in Eq. ( 11), the linear hardness matrix [K E ] and the geometric hardness matrix [K G ] can be expressed as [30]: where Δ is displacements and rotations vector of the shell and δ{ Δ} is virtual displacements and rotations.
The external virtual work exerted by concentrated and distributed out of plane loadings (radial loadings) can be defined as [30]: where P i refer to the concentrated force at point i and δw i indicates the virtual displacement conjugate.Furthermore, q(x,y) expresses the radial distributed exerted loading.

Implementing boundary conditions
In the SFSM, as in the finite element method (FEM), it is necessary to correct the functions used in boundary conditions (BCs) to satisfy the displacement and geometric boundary conditions.Correction of the geometric BCs on the edges of the shell is possible by eliminating the related degrees of freedom, and in order to correct the displacement BCs, it is necessary to correct the spline functions at the beginning and end of the shell, that is, the supporting conditions.Assuming that the degree of freedom f (ξ, η) is in the first support (the beginning of the shell), the nodal line numbers 1 and zero is considered.At the first supporting point, the nodal line number 1 is ξ = 0 and η = − 1.

Analytical integration
In some references, for semi-analytic SFSM which uses harmonic functions in the longitudinal direction of the strip, the analytical (Precise) integration has been used [31].But in the SFSM similar to other numerical methods, so far, numerical integration has been used to solve integral problems [32].This paper presents a method in which all integrals are calculated analytically and precisely.In N. Korkeai et al. this method, because of calculating the integrals before the analysis accurately, the speed of the analysis is increased and the problem solutions can be calculated more accurately.Due to the complexity of calculating the integrals of the spline functions compared to Lagrangian functions, this section discusses how to integrate spline functions accurately.The cubic B-splines in four consecutive sections have non zero values, for example, ∅ i between knots i-2 to i+2 involve non-zero values (Fig. 5 (.
The polynomials of a B3-spline, named S1, S2, S3 and S4, are given as on the other hand, all types of integrals involved in forming the tangent stiffness matrix and residual vector can be presented as where φ ′ i = ∂φ i ∂x .By employing areas S 1 to S 4 expressed for each φ i or φ ′ i , it is possible to transform integration over strip length.Fig. 6 shows a range of spline functions ∅ i , ∅ j and ∅ k , which have certain values at the three points i, j and k in the middle range of the shell (Fig. 6a) as well as for the supporting conditions (Fig. 6b).According to Fig. 6, the common areas between the three spline functions are calculated in three points including two middle and one supporting points as follows: It should be noted that in the MATLAB software environment, a subroutine is written to calculate the spline functions integral for estimating the hardness matrices, residuals vector and load vector:

Eigenvalue problem
Finally, the buckling equation can be expressed as eigenvalue problem as [30]: in which λ represents the buckling load of the structure, which can be calculated by achieving the non-trivial solution of relation ( 19) using the subtype repeat method and subroutine which is written in MATLAB software.

Verification of results
Before presenting the main results, the correctness and accuracy of the attained formulations must be established.Firstly, the convergence checking of computations is systematically ensured in a standard trial & error manner, i.e., by calculating the number of elements and longitudinal Sections accumulating the series truncation constants, while looking for stability in the predicted buckling load.Hence, in Table 1, the axial critical pressure is calculated under uniform loading of an isotropic elliptical cylindrical shell with radial ratio, λ = 0.3 and a length of 200 mm, for N three-node element in crosswise direction and q section in longitudinal direction.It is seen that by increasing the number of the elements and longitudinal sections, the buckling load will be more accurate and closer to the buckling factor of the FE.As can be found this table, for 18 or more three-node elements and 20 or more longitudinal sections, the buckling load has a good convergence.Therefore, in order to reduce the analysis time, in the following sections, 18 nodes with 20 longitudinal sections are used.
Table 2 presents the buckling load of the shell against the different amounts of the lengths and variable ratio of λ using FSDT.For validation purposes, the results of the spline finite strip analysis are checked with those of previous studies [33] and ABAQUS FE N. Korkeai et al. software in Table 2.The area of the elliptical cross section is considered constant and equal to 10 4 π.The value of the shell thickness is supposed to be 1 mm and the shell length is considered to be variable.However, the shell is supposed to be isotropic with E = 210 GPa, v = 0.3.In Table 2, the buckling load of the SFSM analysis is calculated for the number of cross-sectional elements and longitudinal sections.The ABAQUS FE analysis has been conducted with 7600 S4R5 elements, which are dedicated to thin wall shells and small strains.
In another comparison study, the results of the analysis by SFSM for a composite elliptical cylindrical shell are demonstrated in Table 3 and are validated with the results of ABAQUS simulation, when: E 11 = 176 GPa, E 22 = 7.04 GPa, G 12 = G 13 = 3.52 GPa ; G 23 = 1.408GPa and ν = 0.25.The effect of layering type and the radial ratio on the critical load of elliptical cylindrical shell has been investigated under a uniform load for the variable length of the shell.
For validation purposes, in different radial ratios and different values of the cylinder length, the buckling force is inferred from the linear analysis of the finite strip eigenvalue and the FE analysis with ABAQUS software and then the results are compared.In ABAQUS analysis, 8000 S4R5 elements are used, while in the finite strip analysis, 18 three-node elements in crosswise direction and 20 longitudinal sections are applied.According to the conducted studies, both analyzes have a good convergence.From the hypothesized layers, the layering mode of (0 ̗ 90 ̗ 90 ̗ 0) s due to the position angle of the layers to the applied load, has more buckling load than singledirectional layering modes of (90) 8 And (0) 8 and the layering mode of (0 ̗ 90) 2s .It is also observed that for any layering mode at different radial ratios, the more elliptical section approaches to the circle, consequently, the more buckling load will be resulted.Also, the buckling load decreases with increasing shell thickness for any material layering     N. Korkeai et al.Fig. 7 gives the buckling modes for the radial ratio of λ = 0.5 for the layering of (90/90/0) s , (0/90) 2s and (0) 8 .According to the data obtained from the three-dimensional contours, the most changes and the buckling state of the cylinders happened in the (0,90) state.Hence, the maximum differences in the buckling occurred in the middle of the layers, and it decreased from this value around.

Benchmark results
Numerical computations are performed in current section to better understanding of influences of various terms (different types of layering, high semi ellipse, wide semi ellipse and full ellipse) on the variations of buckling load of laminated composite elliptical shell under non-uniform loading.It is necessary to mention that specifications of non-isotropic materials used in the next sections are: In Figs. 8 and 9, the effect of the two types of symmetric and nonsymmetrical layering of (θ ̗ ₋θ) s and (θ ̗ ₋θ) 2 under uniform and nonuniform axial loadings against the various values of radius ratio λ on buckling load is investigated.Fig. 8 will have symmetrical states by increasing the position angle of the layers.For both symmetric and non-symmetrical layering modes of (θ ̗ − θ) s and (θ ̗ − θ) 2 , under uniform and non-uniform axial loading, the layers have the same behavior with increasing angle.Both layering modes at angle of 45 • have the lowest buckling load and in the case of − 20-and − 70 • angles have the highest buckling force.Moreover, according to Fig. 9, the layer mode buckling load of (θ ̗ ₋θ) 2 is always higher than the maximum and minimum values of the layer mode buckling force (θ ̗ ₋θ) s .increase.Note that as the λ ratio increases, the curve approaches a horizontal position and the buckling load decreases.Fig. 8. Buckling force of the composite cylindrical shell under the uniform axial force on the elliptical cross section for different radius ratio and the two layering modes of (θ ̗ ₋θ) s and (θ ̗ ₋θ) 2 .

9.
Buckling load of the composite cylindrical shell under non-uniform axial force on the upper half of the ellipse for different radius ratio and the two layering modes of (θ ̗ ₋θ) s and (θ ̗ ₋θ) 2 .
N. Korkeai et al.Fig. 10 reveals that for a constant radius ratio of the cross section λ = 1.5, and for non-symmetrical layering modes of (θ ̗ ₋θ) 2 , when the load is applied to the parameter of high semi ellipse, the buckling load is higher compared to the states that the load is applied on the perimeter of wide semi ellipse and full ellipse.

Conclusions
According to the first order shear deformation plate assumption, current paper numerically investigates the buckling problem of cylindrical shells with oval cross-section under simply supported boundary condition subjecting to uniform and non-uniform loads via spline finite strip method.Equilibrium governing equations are generated based upon the virtual work's principle.Shell displacement function is assumed to be as dual development of degree three spline functions and Lagrange polynomial functions.The analysis method used at any time with any type of loading can calculate deformation and buckling.The step-by-step critical load is calculated using the determinant of the stiffness matrix.In the present paper, it is supposed that the applied linear elastic materials consist of isotropic and non-isotropic materials.Numerical studies are performed to indicate the effects of shell geometry, materials, and layered materials on the buckling load.For the validation purpose, the results of current research are compared with those of previous researches and ABAQUS software.The results are as follows: • In the SFSM considering the very low number of strips, relatively appropriate answers can be achieved in the short time compared to the FEM.This is attributed to the fact that the fewer degrees of freedom than those of the FEM and the separation of integrals in the two directions, which increases the speed of calculating the hardness.• By keeping the shell volume constant, changing the radius ratio (λ) reduces the buckling force of the composite cylindrical shell.In fact, due to the stretching of the cylindrical section, the structure suffers from a geometric defect which leads to its weakening against the applied load.• By reducing the axial load area, a lower tension occurs in the shell resulting in hardening of the shell and reducing buckling load.
Thus, for a shell with a constant cross section under a non-uniform load on wide semi of ellipse perimeter, the buckling occurs faster compared to the state that non-uniform pressure applies on high semi of the cross-section perimeter.• According to the data obtained from the three-dimensional contours, the most changes and the buckling state of the cylinders happened in the (0,90) state.Hence, the maximum differences in the buckling occurred in the middle of the layers, and it decreased from this value around.

Fig. 1 .
Fig. 1.Geometry and coordinate system of an elliptic cylinder.

Fig. 3 .
Fig. 3.A set of local splines in the longitudinal direction.

Table 1
Convergence study of the buckling force analysis of the isotropic elliptical shell for L = 200.

Table 2
Comparison study of the buckling force analysis of the isotropic elliptical cylindrical shell.

Table 3
Critical buckling force of composite cylindrical shell with different layering.