A high precision shear flexible triangular element for vibration of composite shells

A high precision triangular shallow shell element is proposed and it is applied to free vibration analysis of composite and isotropic shells. The Mindlin’s hypothesis is followed to include the effect of shear deformation. The formulation is made in an efficient manner to make the element free from shear locking problem. The element has some internal nodes, which are eliminated through static condensation technique to improve the computational elegance of the element. In the present vibration problem, the implementation of the static condensation became possible with the help of an efficient mass lumping scheme. It is quite interesting that the effect of rotary inertia can be included in the recommended scheme for lumped mass matrix. Numerical examples covering a wide range of problems are solved and the results obtained are compared with the published results in many cases, which show the precision and range of applicability of the proposed element. The performance of the proposed technique for rotary inertia is found to be excellent. Some new results are produced, which may be useful in future research.


Introduction
The finite element method [2,15,28,29] is regarded as the most powerful tool specifically in structural analysis problems.The analysis of plate and shell is one of the first problems where finite element was applied in early sixties of last century.The initial attempts were made with Kirchhoff's hypothesis where a number of problems were faced.The major problem was concerned with the satisfaction of normal slope continuity at the element edges, which could not be solved yet.In the subsequent study, the Mindlin's hypothesis was followed to include the effect of shear deformation where the above mentioned continuity problem could be avoided.This is achieved by considering the transverse displacement (w) and rotations of the normal (θ x and θ y ) as the independent displacement components.In this group the most popular elements are based on isoparametric formulation.Though these elements are quite elegant but they suffer from certain problems like shear locking, stress extrapolation, spurious modes and something else.
Keeping these aspects in view, some investigations have been carried out to find out some amendments, alternative formulations or some other techniques (e.g.[3,10,13,14,20,27,30]) to achieve an element, which may be applied with improved accuracy without facing the above problems.This is clearly reflected in the recent review papers on shell finite elements [28].The requirement for the development of such elements has been increased with the popularity of fibre reinforced laminated composites as structural material where the effect of shear deformation is quite significant.Actually composite is weak in shear due to its low shear modulus compared to extensional rigidity.Moreover, the layered configuration in these structures causes further complications.The transverse strains (normal as well as shear) are found to be discontinuous at the layer interfaces, which can be represented by zigzag through Fig. 1.A typical shell element projected in the base plane.thickness variation of displacement components [7].On the other hand, transverse stresses become continuous at the layer interfaces [5,6].The effect of these refinements in the evaluation of global response parameters like natural frequencies of vibration is not so prominent as that found in the evaluation of stresses specifically transverse stresses [8].All these have inspired the researchers (e.g.[1,4,5,8,16,21,[23][24][25][26]) to make this an active area of research in recent days.However, the present study is based on Mindlin's hypothesis in order to have a simple representation of the structural deformation.
In this context, the triangular Mindlin plate element developed by one of the authors (Sengupta [20]) is quite interesting.This trouble free element based on a simple idea is found to give very good results even with a coarse mesh in the static analysis of isotropic plates [20].It has inspired the authors to develop the proposed shallow shell element using the basic concept of the high precision triangular Mindlin plate element [20].The performance of the element has been found to be very good in the prediction of structural response under static load [21] while that is tested for the prediction of vibration characteristics in the present study.In the shell element the transverse displacement w is expressed by a complete fourth order polynomial while complete cubic polynomials are used to express the in-plane displacements (u and v) and the rotations (θ x and θ y ) of the normal.Thus the interpolation function of w is one order higher than that of θ x and θ y , which has helped to make this element free from locking and other relevant problems.
The fifty-five unknowns in these five polynomials can be easily expressed in terms of fifty-five nodal displacements of the element as shown in Fig. 1.With this the stiffness and mass matrices of an order of fifty-five can be derived following the usual steps of finite element technique.Though these matrices can be used directly in these forms but their implementation and use will be quite cumbersome due to the presence of the internal nodes.This can be avoided by eliminating the degrees of freedom of the internal nodes.Actually Sengupta [20] has done it using the technique of static condensation, which is quite easy to implement in a static analysis [20].This is somewhat difficult in the present vibration problem due to the presence of the mass matrix in the governing equation.The problem will be much more severe if a consistent mass matrix is used where the matrix is more populated and having coupling of the degrees of freedom for the internal and external nodes.On the other hand the effect of rotary inertia can easily be incorporated in the formulation with consistent mass matrix.In the present element, the condensation technique of Guyan [11] and something else cannot be applied, as it contributes a significant amount of mass at the internal nodes.In this situation, the only alternative is to have a mass matrix having no mass contribution at the internal nodes.This is achieved with the proposed mass lumping scheme where the mass of an element is distributed at the external nodes only.Interestingly, some recommendation has also been made for the incorporation of rotary inertia in lumped mass matrix.The detail of the lumping scheme is presented in the formulation.
Numerical examples of isotropic and composite shells having different curvature, aspect ratio, boundary condition, number of layers, fibre orientation and thickness ratio (h/a) are solved by the proposed element.A large number of these results are compared with the results available in literature, which clearly reflects the performance of the proposed element.

Formulation
The formulation is based on shallow shell theory with the usual assumptions.The effect of shear deformation is taken into account following the Mindlin's hypothesis.The area co-ordinate system [29] is used to do the formulation.
According to earlier discussion the independent displacement components may be expressed as follows and where Table 1 Frequency parameters (ωa 2 (ρh)/D) of an isotropic doubly curved shell panel (Fig. 2) clamped at y = 0 and free at the other sides (a/b = 1, a/Rx = 0.2, ν = 0. The unknowns ({α}, {β}, {γ}, {µ} and {λ}) in the above Eq.( 1) can be expressed in terms of nodal displacements by substitution of these equations at the different nodes appropriately, which leads to and the matrix [A] having an order of 55 × 55 can be obtained with the co-ordinates of the different nodes.
The generalised stress strain relationship of a laminated shell may be written as Table 2 Frequency parameters (100ωa ρ/E) of an isotropic cylindrical shell panel (Fig. 3) simply supported at the four sides (a/b = 1, a/R = 0.5, ν = 0. where ) Table 3 Frequency parameters (ωa 2 ρ/E 2 h 2 ) of a cross-ply spherical shell panel (Fig. 2) simply supported at the four sides The rigidity matrix [D] constitutes of the contributions of its individual layers.Using the material properties and fiber orientation of the individual layers it can be easily obtained following the steps available in any standard text on mechanics of fiber reinforced laminated composites.For the evaluation of rigidity parameters for transverse shear (A 44 , A 45 , A 54 and A 55 ), a shear correction factor of 5/6 is used.Now the field variables as expressed in Eq. ( 1) may be substituted in Eq. ( 5) to express the strain vector {ε} as With the help of Eq. ( 2) the strain vector {ε} in Eq. ( 7) may be expressed in terms of nodal displacement vector {X} as Once the matrices [B] and [D] are obtained, the element stiffness matrix [K e ] can be easily derived using the help of Virtual work technique and it may be expressed as In a similar manner the consistent mass matrix of an element can be derived with the help of the Eqs (1) and ( 2) and it may be expressed as Table 4 Frequency parameters (ωa 2 ρ/E 2 h 2 ) of a cross-ply (0/90/0) cylindrical shell panel (Fig. 3) simply supported at the four sides (a/b = 1, where In the above equation [N 1 ] and [N 2 ] are null matrices of order 1 × 10 and 1 × 15 respectively.The first three terms of the mass matrix in Eq. ( 10) are associated with the movement of mass along u, ν and w respectively, which are found to contribute the major inertia.The last two terms are associated with rotary inertia and their contribution is expected to be significant in thick shell only.The integration found in Eqs ( 9) and ( 10) is carried out numerically [30] to evaluate the element stiffness and mass matrices.
Though the consistent mass matrix presented in Eq. ( 10) includes all the contributions including rotary inertia but the difficulty of its use in this form is discussed earlier.It has also been mentioned that the difficulty can be overcome by using a lumped mass matrix based on the mass lumping scheme presented below.In this context two similar mass lumping schemes are proposed where the effect of rotary inertia is taken into account in one case.
In the first lumping scheme (MLS-1), the mass of an element m e is distributed at w of its external nodes where the ratio of distribution is dependent on the corresponding diagonal elements of the consistent mass matrix [M e ] presented in Eq. (10).This is similarly done for u, v and it is as follows Table 5 Frequency parameters (ωa 2 ρ/E 2 h 2 ) of an angle-ply (30/-30/...) cylindrical shell panel (Fig. 3) simply supported at the four sides and where m ul i,i , m vl i,i and m wl i,i are the ith diagonal elements corresponding to u, v and w of the proposed lumped mass matrix [M l ] and m i,i is the ith diagonal element of the consistent mass matrix [M e ].
In the second lumping scheme (MLS-2), the effect of rotary inertia is considered in addition to those (m ul i,i , m vl i,i , and m wl i,i ) taken in the previous case (MLS-1).These additional mass contributions are taken in the following manner Table 6 Frequency parameters (ωa 2 ρ/E 2 h 2 ) of an angle-ply (45/-45/45) spherical shell panel (Fig. 2) having different boundary conditions at the four sides (Rx/a = 3.0, a/b = 1, results, the structure is analysed with an existing shear deformable element to compare the results obtained by the proposed element.In this context a separate computer program (ISSE9) based on nine noded isoparametric element is written where the assumptions are identical to those used in the proposed element.

Isotropic doubly curved shell
A cantilevered doubly shell panel as shown in Fig. 2 (clamped at y = 0) is analysed for three different values of R y /R x (-2.0, -1.0 and 2.0) taking h/a = 0.01 and 0.1.The analysis is carried out by the proposed element with three different mesh divisions taking MLS-2 as the mass lumping scheme.The first six natural frequencies obtained by the proposed element are presented in Table 1 with the analytical solution of Leissa et al. [17].As the effect of shear deformation and rotary inertia is not considered by Leissa et al. [17], the structure is also analysed by the isoparametric element (ISSE9) and the results obtained are included in Table 1.The results indicate the performance of the proposed element in doubly curved shell panel having complex geometry.In one case (h/a = 0.01 and R y /R x = −2) the analysis with the isoparametric element is carried out with four different mesh sizes (see Table 1).It shows that the isoparametric element requires higher mesh size compared to the proposed element to attain the convergence or to get the desirable accuracy.This infers that the proposed element possesses improved performance compared to isoparametric element.This feature has been highlighted in details by Sengupta [20] with his plate element based on similar concept.For the isoparametric element, a mesh size of 20x20 appears to be sufficient and it is used in the other cases.This is also followed in the subsequent examples where it is used.

Isotropic cylindrical shell
A cylindrical shell panel as shown in Fig. 3, is analysed by the proposed element with five different mesh divisions taking simply supported boundary condition at the four sides.Using both the mass lumping schemes (MLS-1 and MLS-2) the analysis is carried out for thickness ratio (h/a) of 0.1 and 0.2.The first six natural frequencies obtained in the present analysis are presented in Table 2 with those of Lim and Liew [18] who have solved the problem analytically using higher order shear deformation theory (HSDT).They [18] have also presented some results based on Mindlin's hypothesis i.e., first order shear deformation theory (FSDT).Table 2 shows the performance of the proposed element in terms of convergance rate and solution accuracy.It also indicates that lumping scheme MLS-2 has performed better particularly in the case of higher thickness ratio and it is expected as discussed earlier.This has been throughly studied by Haldar [12] taking different types of plate and shell problems.Based on these observations, MLS-2 may be recommended for the analysis of shell of any thickness while MLS-1 may be used when thickness ratio is reasonable small.

Cross-ply spherical shell
A cross-ply spherical shell panel on square base simply supported at the four sides is analysed with the proposed element (mesh size: 8 × 8, 10 × 10 and 12) taking MLS-2 as mass lumping scheme.The study is made for unsymmetrical (0/90) and symmetrical (0/90/90/0) ply arrangements taking curvature ratio R x /a = R y /b = 3, 5 and 10 (Fig. 2).In all the cases the thickness ratio (h/a) is taken as 0.1.The first six natural frequencies obtained by the proposed element are presented in Table 3 with the fundamental frequency of Reddy [19] obtained analytically.To compare the results for higher modes the problem is also solved by the isoparametric element and the frequencies obtained are included in Table 3.The table shows that the results obtained from different sources agreed well.

Cross-ply cylindrical shell
A cross-ply (0/90/0) cylindrical shell panel simply supported at the four sides is analysed for R x /a = 4, 5 and 10 (Fig. 3) taking h/a = 0.01 and 0.1.Similar to the previous example the analyis is carried out with the proposed element (mesh size: of 8 × 8, 10 × 10 and 12 × 12) using mass lumping scheme MLS-2 and the isoparametric element (mesh size: 20 × 20).This analysis scheme is followed in the solution of the problems in the subsequest examples.The first six natural frequencies obtained by both the elements are presented with the fundamental frequency of Reddy [19] in Table 4, which shows that the results agreed well.The material properties and other data are given in Table 4, which is done in the other tables also.

Angle-ply cylindrical shell
An angle-ply (30/-30/...) cylindrical shell panel simply supported at the four sides is analysed for the subtended angle φ = 20 • , 30 • and 45 • (Fig. 3) taking h/a = 0.05.The analysis is done for two layer (30/-30) and layer (30/-30/30/-30) ply arrangements using the proposed element and the isoparametric element.The first six natural frequencies obtained by these elements are presented with the fundamental frequency of Soldatos [22] in Table 5.It shows that the results obtained by the proposed element are close to those obtained from other means.

Angle-ply spherical shell
An angle-ply (45/-45/45) spherical shell panel having a square base is analysed by the proposed element and the isoparametric element taking R x /a = R y /b = 3 (Fig. 2) and h/a = 0.1 and 0.2.The study is made for different combination of boundary conditions at the four sides, which include simply supported, clamped and free edges.The first six natural frequencies obtained in the present analysis are presented in Table 6, which shows that the agreement between the results obtained by the two elements is good.

Conclusions
A high precision triangular shallow shell element is proposed.The effect of shear deformation is incorporated in the formulation and it is done in such a manner, which has made the element free from locking in shear.This is based on the concept used in a shear deformable plate element developed by one of the authors of this paper.An effective use of the element in vibration analysis is due the mass lumping scheme proposed in this study.This may be considered as one of the contributions of this paper, as it can be used in any dynamic analysis with an element having internal nodes.The element is applied to the analysis of isotropic and composite shells having different geometry, thickness ratio, stacking sequence, boundary condition and some similar features.The results obtained by the proposed element are compared with the results obtained from other sources, which shows the performance of the element in terms of accuracy and range of applicability.A number of new results particularly for the composite shells are generated.It is expected that these new results will be useful in future research.

Fig. 2 .
Fig. 2. Geometry and axis system of a doubly curved shell panel.

Fig. 3 .
Fig. 3. Geometry and axis system of a cylindrical shell panel.
bResults based on higher order shear deformation theory (HSDT).c Results based on first order shear deformation theory (FSDT).