An efficient general curvilinear coordinates finite element method for the linear dynamic study of thickness-independent shells

Abstract

yyTo date, a large number of finite element methods have been developed to study the dynamics of shell structures. Most of them are generally based on the degenerated solid approach and other less in shell theories, but introducing, in this last case, some assumptions to analyze this problem: some of them refer to shallow shells (slightly curved shells), or consider thin shells neglecting shear deformation, or dispense some terms in their stress-strain developments like the off-diagonal components of the curvature or metric tensors (orthogonal coordinates). In the present work, we present an improved finite element method for the linear dynamic analysis of shells, from thin to moderately thick and thick shells, developed in general curvilinear coordinates, based on a refined shear deformation shell theory and free of the well-known shear locking effect. Exact constitutive equations, including higher order moments-strains relations, are also deduced for the adequate analysis of thick shells. To circumvent the shear locking problems, the mixed interpolation of the tensorial components (MITC) of the linear strain tensor is used. An exhaustive study of different surfaces is performed, especially in doubly curved shells, and interesting conclusions of the higher order modes of vibration and the strain energy of the element are derived. Other desirable features like a low computational effort, a straightforward extension to nonlinear formulation and applications for composite shells are found in this novel and general formulation. Very good results in the proposed practical cases have been found.

Keywords:
Vibrations; finite elements; moderately thick and thick shells

1 INTRODUCTION

The vibrations of shells are a topic of undoubted interest in the field of civil, mechanical and aeronautical engineering. Today, we find many examples where this type of phenomenon occurs, for example, the fuselage of supersonic jets, or the dynamic phenomena to which structures in Civil Engineering are subjected such as wind, earthquakes, impacts, etc., which may induce oscillatory motions in structures. There are structural elements that we can study, such as plates or shells, and their dynamics must be calculated taking this into account, Aksu (1997).

The first studies of free vibrations in structural elements date back to the 1800s. The pioneering works of Germain (1821) and later Kirchhoff (1850), introduced the famous biharmonic equation relating the deflections of a plate with the transverse loads applied. However, this theory applies only to thin plates: the displacements obtained are lower than in reality, and it overestimates the stresses, natural frequencies, and buckling loads.

It was Reissner (1945) who introduced the shear deformation in calculations for moderately thick plates and shells and even proposed the dynamic formulation of the problem including the rotary inertia in the kinetic energy term. Since then, there have emerged numerous theories, the so-called higher-order shear deformation plate theories (HOSDPT), with infinite variants that analyze both dynamic and static problems. A good review of these theories can be found in Mantari et al (2011) and Meiche et al. (2011). Analytical solutions for these equations are rarely possible, so it has to resort to approximate methods.

With respect to shell structures, the book of Wempner and Talaslidis (2002) offers a comprehensive and rigorous treatment of theories of deformable bodies and shells using general tensor analysis. The compendium of Kienzler et al. (2004) treats different shell theories, reviewing the current state of the art and higher-order shear deformation shell theories including zigzag theories and layerwise laminated shell theories.

In a similar way to the calculations for plates, Alhazza and Alhazza (2004) and Leissa (1993) have shown that closed-form solutions for general shell structures under dynamic actions are not obtainable. Bhimaraddi (1991) studies the free vibration of doubly curved shallow shells with zero twist on a rectangular planform using the three-dimensional elasticity equations. He studies the equilibrium equations of the shells using double Fourier series under the assumption that the ratio of the shell's thickness to its middle surface radius is negligible compared to unity. The difficulty, in the current state of knowledge, in obtaining general analytical solutions to the proposed problem, including shear deformation, requires the use of numerical methods for the modelling of the problem.

The finite element method offers a wide range of possibilities to address this problem, see for example the book of Zienkiewicz and Taylor (2000). Yang et al. (2004) present a good review of the finite element methods for shell structures that have been developed since its inception back in the 1960's, presenting different techniques adopted by many authors to formulate efficient and versatile shell finite elements free of the well known problem of the different kinds of locking of the stiffness matrix. Membrane and shear lockings are distinct problems associated with shear deformation shell elements with thin thicknesses, and the trapezoidal (curvature thickness locking) and volumetric lockings are related to “solid-shell” elements and 3D elements. Enhanced assumed strain (EAS) elements, assumed natural strain (ANS) elements, elements based on multifield principles, like the Hellinger-Reissner and Hu-Washizu principles, are some of the many methods to address these problems and are presented in that review.

Within the wide range of finite element methods proposed to date, most are called degenerated solid-shell elements, as originally proposed by Ahmad et al. (1970), which discretize the equations of the 3D-continuum concurrently introducing physical assumptions at discrete points.

With this approach, Lee and Han (2001) have studied the free vibration of isotropic plates and shells using a modified nine-node degenerated shell element with the ANS method.

Others have concentrated their efforts on specific geometries, for example, Kang and Leissa (2005), who study paraboloid shells of revolution with variable thickness and solid paraboloids.

Finite element formulations based on different refined shell theories have been used profusely as well. Although they are apparently different approaches, they share a common physical meaning, as shown in various works even in the nonlinear geometric regime, Büchter and Ramm (1992).

There are not many publications studying shell structures in the dynamic regime including shear deformation and large curvatures. For example, Liew and Lim (1995, 1996), Stavridis and Armenakas (1988) and Stavridis (1998) propose shallow shell elements also using the degenerated approach. Others, like Aksu (1997) and Kumbasar and Aksu (1995), formulate shell elements valid for moderately thick and really thick shells but assume orthogonal curvilinear coordinates along the lines of curvature of the reference surface of the shell. As shown by Qatu (1999), shear deformation higher-order theories are also available for the dynamic analysis of shells.

The pioneering work of Reddy and Liu (1985), who developed a 2-D higher-order theory for laminated elastic shells accounting for a parabolic distribution of the transverse shear strains through the thickness of the shell and tangential stress-free boundary conditions on the boundary surface of the shell, was the starting point for a series of theories which try to reproduce the behavior of thick shells, especially in the case of high frequencies. In this regard, Fard et al. (2014) presents a good review of different variants of the finite element method which study this problem and proposes an improved higher order double curved sandwich shell theory to analyze the free vibrations of double curved thick composite sandwich panels. In addition, they considered the terms $\left(1+\frac{z}{R_{\alpha }}\right)$, corresponding to the components of the mixed curvature tensor expressed in orthogonal coordinates, and integrated them exactly through the thickness. Leissa and Chang (1996) demonstrated that by considering the terms $\left(1+\frac{z}{R_{\alpha }}\right)$ and using first-order shear deformation shell theories, better results are obtained than with higher-order theories neglecting these terms.

As far as the present authors' knowledge extends, there is little work in the literature dealing with shell elements using non-orthogonal reference systems on the shell, however, there are the excellent contributions of Chapelle and Bathe (1998), Wempner and Talaslidis (2002).

It is important to note that other variants of the classical MEF have been developed to study the vibrations of shells: the element free Galerkin methods and other meshless methods, like the generalized differential quadrature and the radial basis functions, for example, have been successfully implemented by Liu et al. (2006) and Tornabene et al. (2014), respectively.

A different way to tackle the vibrations of shells is with the use of 3D elements. Such elements have a number of advantages, such as the possibility of using three-dimensional constitutive laws and not making any kinematic assumptions in deriving the strain--displacement relations, which is an implicit hypothesis in the shell elements. As for disadvantages, the major computational effort required, as well as the different kinds of locking of the stiffness matrix, particularly pronounced for low order elements, are the most significant. As with shell elements, different techniques have been developed to circumvent these problems. An interesting review of the state of art for these elements, proposing a cure for the locking problem, is Schwarze and Reese (2009).

In the present paper, we formulate a general shell element in general curvilinear coordinates, not necessarily orthogonal coordinates, based on a consistent shell theory and valid for any kind of surface. Despite an initial major knowledge of the geometry of the shell, including the metric tensor, the Christoffel symbols, and the curvature tensor, the construction of the stiffness matrix and the mass matrix of the element is completely general. In order to ease the shear locking of the stiffness matrix for thin shells and bending-dominated problems Chapelle and Suarez (2008), we have opted for the method of mixed interpolation of tensorial components (MITC). Originally introduced by Bathe and Dvorkin (1986) for low order shell elements, it fits comfortably into the proposed formulation.

The proposed shell element is versatile, lacks of zero energy modes (hourglass modes) and is free of the locking effect for thin shells. It also has some other desirable features, such as a uniform convergence to the exact solution independently of the shell geometry and thickness.

In the presented examples, we compare the obtained solution with the results obtained by other competitive finite elements, Liew and Lim (1995, 1996), and add new results for doubly curved shells of moderate thickness.

2 SHELLS GEOMETRY

For further developments, we will recall some important concepts related to the geometry of shells and we will use the concepts of tensor calculus in profusion. An excellent book that the reader can consult for more specific aspects of tensor analysis is Nayak (2012).

In order to work with shells, we admit that the reference surface of the shell can be parameterized by the curvilinear coordinates ${\vartheta }_{\alpha }(\alpha =\mathrm{1,2})$, not necessarily orthogonal.

Figure 1:
Curvilinear local reference frame at the medium surface of the shell.

As we know, in the undeformed configuration of the shell, the position vector $r$ of a point can be written as,

$r=r_0+{\alpha }^3{a}_3$ (1)

where $r_0$ is the position vector of a point in the middle surface, ${\alpha }^3$ is the material coordinate measuring the length along the normal to the to the mid-surface, and $a_3$ is the normal vector to the mid-surface, see Figure 1.

If we differentiate this expression, we obtain,

$g_{\alpha }=r_{,\alpha }={\left(r_0+{\alpha }^3{a}_3\right)}_{,\alpha }=a_{\alpha }+{\alpha }^3{a}_{\mathrm{3,}\alpha }=\left({\delta }_{\alpha }^{\rho }-{\alpha }^3{b}_{\alpha }^{\rho }\right)a_{\rho }$ (2)

where $g_{\alpha }$are the covariant base vectors which define the plane tangent to an arbitrary surface of the shell and $b_{\alpha }^{\rho }$are the components of the mixed curvature tensor.

The relation between the metric tensors of the shell space $g_{\alpha \beta }$ and the middle surface $a_{\alpha \beta }$ takes the form,

$g_{\alpha \beta }=a_{\alpha \beta }-2{\alpha }^3{b}_{\alpha \beta }+{\left({\alpha }^3\right)}^2{c}_{\alpha \beta }$ (3)

where $b_{\alpha \beta }$ is the curvature tensor, the second fundamental form of the reference surface. and $c_{\alpha \beta }$ is the third fundamental form, a symmetrical tensor whose expression is,

$c_{\alpha \beta }=b_{\alpha }^{\gamma }{b}_{\gamma \beta }$ (4)

The derivatives of the tangent vectors and the normal vector to the surface can be obtained from the Weingarten equations,

$a_{\alpha ,\beta }={\Gamma }_{\alpha \beta }^{\gamma }{a}_{\lambda }+b_{\alpha \beta }{a}_3$ (5)

$a_3{}_{,\beta }=-b_{\beta }^{\lambda }{a}_{\lambda }$ (6)

In these equations we have introduced the Christoffel symbols of the second kind ${\Gamma }_{\alpha \beta }^{\gamma }$, which are not tensors.

Equation (2) can be written in a more compact way as,

$g_{\alpha }=\mu a_{\alpha }$ (7)

where

$\mu :=g^i\otimes a_i$ (8)

$\mu$ is the so-called shifter tensor, the two point tensor which transforms the components of a vector in one coordinate system, the shell space, to the components of the same vector in another coordinate system, the reference surface, and has the expression,

$\mu ={\delta }_{\alpha }^{\rho }-{\alpha }^3{b}_{\alpha }^{\rho }.$ (9)

This contra-covariant tensor plays a fundamental role in the development of the theories of shells because it allows us to relate the tangent vectors to the medium surface and those related to any other surface of the shell (shell space).

For second order tensors, like the metric tensor, we have,

$g_{ij}={\mu }_i^a{\mu }_j^b{a}_{ab}$ (10)

If we consider, for example, a fourth order tensor, like the elastic constitutive tensor ${\overline{C}}^{\alpha \beta \lambda \mu }$, we can write

$C^{\alpha \beta \lambda \mu }{\overrightarrow{a}}_{\alpha }\otimes {\overrightarrow{a}}_{\beta }\otimes {\overrightarrow{a}}_{\lambda }\otimes {\overrightarrow{a}}_{\mu }={\overline{C}}^{\alpha \beta \lambda \mu }{\overrightarrow{g}}_{\alpha }\otimes {\overrightarrow{g}}_{\beta }\otimes {\overrightarrow{g}}_{\lambda }\otimes {\overrightarrow{g}}_{\mu }$ (11)

And hence,

${\overline{C}}^{\alpha \beta \lambda \mu }=C^{\rho \sigma \tau \omega }{\left({\mu }^{-1}\right)}_{\rho }^{\alpha }{\left({\mu }^{-1}\right)}_{\sigma }^{\beta }{\left({\mu }^{-1}\right)}_{\tau }^{\lambda }{\left({\mu }^{-1}\right)}_{\omega }^{\mu }$ (12)

Similarly, one can obtain the relation between any geometric features of any surface of the shell and those of the middle surface through the shifter tensor.

The Christoffel symbols of the 2nd kind can be expressed as:

${\overline{\Gamma }}_{\alpha \beta }^{\gamma }=g_{\alpha ,\beta }{g}^{\gamma }={\left({\mu }_{\alpha }^{\delta }{a}_{\delta }\right)}_{,\beta }·{\left({\mu }^{-1}\right)}_{\xi }^{\gamma }{a}^{\xi }$ (13)

Now further developing this expression, we obtain

${\tilde{\Gamma }}_{33}^{\gamma }={\tilde{\Gamma }}_{3\alpha }^3={\tilde{\Gamma }}_{\alpha 3}^3=0$ (14)

${\tilde{\Gamma }}_{\alpha \beta }^{\gamma }={\tilde{\Gamma }}_{\beta \alpha }^{\gamma }={\Gamma }_{\alpha \beta }^{\gamma }-t{\zeta }_{\lambda }^{\gamma }{b}_{\alpha |\beta }^{\lambda }$ (15)

${\tilde{\Gamma }}_{\alpha \beta }^3=b_{\alpha \beta }-t{c}_{\alpha \beta }$ (16)

${\tilde{\Gamma }}_{3\beta }^{\alpha }=-{\zeta }_{\lambda }^{\alpha }{b}_{\beta }^{\lambda }$ (17)

with ${\zeta }_{\lambda }^{\alpha }$as the inverse of the shifter tensor,

${\zeta }_{\lambda }^{\alpha }{\mu }_{\beta }^{\gamma }={\delta }_{\beta }^{\alpha }$ (18)

${\zeta }_{\beta }^{\alpha }=\frac{1}{\mu }{\delta }_{\beta \lambda }^{\alpha \gamma }{\mu }_{\gamma }^{\lambda }$ (19)

These expressions are very interesting for developments concerning thick shells or those concerning high order theories, where it cannot be assumed that the geometrical properties of the shell space are not independent of the thickness.

To finish this section, it is very important for further developments, to know the expression of the differential volume $dV$ in terms of known quantities of the middle surface.

By definition,

$dV=g=\det g_{\alpha \beta }=\left[g_1{g}_2{g}_3\right]$ (20)

and taking into account the relation (7), we have,

$dV=g=\det g_{\alpha \beta }=\left[g_1{g}_2{g}_3\right]=\left({\delta }_1^{\rho }-{\alpha }^3{b}_1^{\rho }\right)a_{\rho }\left({\delta }_2^{\gamma }-{\alpha }^3{b}_2^{\gamma }\right)a_{\gamma }{a}_3=\left({\delta }_1^{\rho }-{\alpha }^3{b}_1^{\rho }\right)\left({\delta }_2^{\gamma }-{\alpha }^3{b}_2^{\gamma }\right)a{e}_{\rho \gamma 3}=a\left|{\delta }_{\beta }^{\alpha }-{\alpha }^3{b}_{\beta }^{\alpha }\right|=a\left[1-2{\alpha }^{3}H+{\left({\alpha }^3\right)}^{2}K\right],$ (21)

where H y K are the mean and Gaussian curvature of the middle surface, invariants of the mixed curvature tensor.

The relation between H and K can be established from,

$K{a}_{\alpha \beta }-2H{b}_{\alpha \beta }+c_{\alpha \beta }=0$ (22)

Hence, we can write,

$dV=\mu dA=\left(1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right)dA$ (23)

which relates the differential of the volume to the area as a function of the mean and Gaussian curvature of the middle surface.

3 SHELLS KINEMATIC

The expression of the components of the linear strain tensor can be written as, Wempner (1982)

${\epsilon }_{\alpha \beta }=\frac{1}{2}\left(V_{\alpha \parallel \beta }+V_{\beta \parallel \alpha }\right)$ (24)

where $V$ are the displacements of any point of the shell space and the double vertical bar denotes the covariant derivative with respect to a general basis of vectors tangent to a general surface of the shell space, not necessarily the middle shell surface. Recall that for a first order tensor, the covariant derivative has a more explicit form, given by,

$V_{\alpha |\beta }=V_{\alpha ,\beta }-{\Gamma }_{\alpha \beta }^{\gamma }{V}_{\gamma }$ (25)

As seen above, the tangent vectors to any surface of the shell $V_{\alpha }$ and those tangents to the middle surface $v_{\beta }$ are related through the shifter tensor by

$V_{\alpha }={\mu }_{\alpha }^{\beta }{v}_{\beta }$ (26)

$V_3=v_3$ (27)

At this point, we can express the covariant derivatives of the three dimensional components of the displacement vectors as a function of the two-dimensional components of the displacement vectors through the shifter tensor by,

$V_{\alpha \Vert \beta }={\mu }_{\alpha }^{\lambda }\left(v_{\lambda |\beta }-b_{\lambda \beta }{v}_3\right)$ (28)

So we can rewrite the expression (24),

${\epsilon }_{\alpha \beta }=\frac{1}{2}\left({\mu }_{\alpha }^{\lambda }\left(v_{\lambda |\beta }-b_{\lambda \beta }{v}_3\right)+{\mu }_{\beta }^{\lambda }\left(v_{\lambda |\alpha }-b_{\lambda \alpha }{v}_3\right)\right)$ (29)

where the single bar denotes the covariant derivative in the two-dimensional shell surface space.

Considering the expression of the shifter tensor, the generalized strains can be written

${\epsilon }_{\alpha \beta }=\left(\frac{1}{2}\left(v_{\alpha |\beta }+v_{\beta |\alpha }\right)-b_{\alpha \beta }{v}_3\right)+\left(c_{\alpha \beta }{v}_3-\frac{1}{2}\left(b_{\alpha }^{\lambda }{v}_{\lambda |\beta }+b_{\beta }^{\lambda }{v}_{\lambda |\alpha }\right)\right)\left({\alpha }^3\right)$ (30)

${\epsilon }_{\alpha 3}=\frac{1}{2}\left(V_{\alpha \Vert 3}+V_{3\Vert \alpha }\right)=\frac{1}{2}\left({\mu }_{\alpha }^{\lambda }{v}_{\lambda \mathrm{,3}}+v_{\mathrm{3,}\alpha }+b_{\alpha }^{\lambda }{v}_{\lambda }\right)$ (31)

${\epsilon }_{33}=v_{\mathrm{3,3}}$ (32)

At this point we must introduce the corresponding kinematic hypothesis, which in our case it is a first shear deformation shell theory.

$v_{\alpha }=u_{\alpha }+{\alpha }^3{\theta }_{\alpha }$ (33)

It physically indicates that after the deformation process, the fibers normal to the middle surface remain straight but not necessarily normal. After a more or less lengthy calculation, we obtain:

${\epsilon }_{\alpha \beta }=\left(\frac{1}{2}\left(u_{\alpha |\beta }+u_{\beta |\alpha }\right)-b_{\alpha \beta }{u}_3\right)+\left(\frac{1}{2}\left({\theta }_{\alpha |\beta }+{\theta }_{\beta |\alpha }-b_{\beta }^{\lambda }{u}_{\lambda |\alpha }-b_{\alpha }^{\lambda }{u}_{\lambda |\beta }\right)+c_{\alpha \beta }{u}_3\right)\left({\alpha }^3\right)+\left(\frac{1}{2}\left(b_{\beta }^{\lambda }{\theta }_{\lambda \mid \alpha }+b_{\alpha }^{\lambda }{\theta }_{\lambda \mid \beta }\right)\right){\left({\alpha }^3\right)}^2$ (34)

${\epsilon }_{\alpha 3}=\frac{1}{2}\left({\theta }_{\alpha }+u_{\mathrm{3,}\alpha }+b_{\alpha }^{\lambda }{u}_{\lambda }\right)$ (35)

So the linear strain tensor,$E_{\alpha \beta }$ can be expressed in matrix form as follows,

$E_{\alpha \beta }=\left(\begin{array}{cc}{\epsilon }_{\alpha \beta }& {\epsilon }_{\alpha 3}\\ {\epsilon }_{3\alpha }& {\epsilon }_{33}\end{array}\right)=\left(\begin{array}{cc}\left(e_{\alpha \beta }+{\alpha }^3{\rho }_{\alpha \beta }+{\left({\alpha }^3\right)}^2{\kappa }_{\alpha \beta }\right)& \frac{1}{2}{\zeta }_{\alpha }\\ \frac{1}{2}{\zeta }_{\alpha }& 0\end{array}\right)$ (36)

where,

$e_{\alpha \beta }=\frac{1}{2}\left(u_{\alpha |\beta }+u_{\beta |\alpha }\right)-b_{\alpha \beta }{u}_3$ (37)

${\rho }_{\alpha \beta }=\frac{1}{2}\left({\theta }_{\alpha |\beta }+{\theta }_{\beta |\alpha }-b_{\beta }^{\lambda }{u}_{\lambda |\alpha }-b_{\alpha }^{\lambda }{u}_{\lambda |\beta }\right)+c_{\alpha \beta }{u}_3$ (38)

${\kappa }_{\alpha \beta }=\frac{1}{2}\left(b_{\beta }^{\lambda }{\theta }_{\lambda \mid \alpha }+b_{\alpha }^{\lambda }{\theta }_{\lambda \mid \beta }\right)$ (39)

${\zeta }_{\alpha }={\theta }_{\alpha }+u_{\mathrm{3,}\alpha }+b_{\alpha }^{\lambda }{u}_{\lambda }$ (40)

Taking into account that for thin shells we can introduce the Kirchhoff-Love assumption, normals to the shell middle surface remain straight and normal to it during deformations and unstretched in length

${\theta }_{\alpha }=-\left(u_{3|\alpha }+b_{\alpha }^{\lambda }{u}_{\lambda }\right)$ (41)

In this case, displacements of a point of the shell space can be related to those in the reference surface by,

$v_{\alpha }=u_{\alpha }-{\alpha }^3\left(u_{3|\alpha }+b_{\alpha }^{\lambda }{u}_{\lambda }\right)a^{\alpha }$ (42)

and the curvature tensor has the expression,

${\rho }_{\alpha \beta }=\frac{1}{2}{u}_{3|\alpha \beta }+b_{\beta |\alpha }^{\lambda }{u}_{\lambda }+b_{\beta }^{\lambda }{u}_{\lambda \mid \alpha }+b_{\alpha }^{\lambda }{u}_{\lambda \mid \beta }-b_{\alpha }^{\lambda }{b}_{\lambda \beta }^{}{u}_3$ (43)

Note that in case we use this expression in finite element formulations, we have to impose first order continuity between shell elements, which is difficult to achieve for arbitrary geometries when using standard polynomials as shape functions.

4 STATEMENT OF THE PROBLEM.

There are different ways to get to the governing equations of the dynamic problem, either for analytical or approximate solutions. In the present work, we choose Hamilton's principle.

The problem to be solved is reduced to minimizing the functional:

$\delta \underset{time}{\overset{}{{\displaystyle \int }}}\underset{A}{\overset{}{{\displaystyle \iint }}}{L}^{\text{*}}\sqrt{a}d\xi d\eta dt=0$ (44)

$\begin{array}{l}\delta \underset{time}{\overset{}{{\displaystyle \int }}}\underset{A}{\overset{}{{\displaystyle \iint }}}{L}^{\text{*}}dt=\delta \underset{time}{\overset{}{{\displaystyle \int }}}\left[{\displaystyle \underset{v}{\iiint }(T-{\pi }^{\mathrm{int}})dV}+{\displaystyle \underset{A}{\iint }\left(\tau {*}^{\alpha \beta }{\nu }_{\beta }{\mu }_{\alpha }^{\gamma }{v}_{\gamma }+\tau {*}^{3\beta }{\nu }_{\beta }{v}_3\right)}dA\right]dt=\\ \delta \underset{time}{\overset{}{{\displaystyle \int }}}\left[{\displaystyle \underset{v}{\iiint }\left[\left(\frac{1}{2}\rho {\dot{V}}^{\alpha }{\dot{V}}_{\alpha }+{\dot{V}}_3{}^2\right)-{\pi }^{\mathrm{int}}\right]dV}+{\displaystyle \underset{A}{\iint }\left(\tau {*}^{\alpha \beta }{\nu }_{\beta }{\mu }_{\alpha }^{\gamma }{u}_{\gamma }+\tau {*}^{3\beta }{\nu }_{\beta }{u}_3\right)}dA\right]dt=0\end{array}$ (45)

where $L^{*}$ is the called Lagrangian density function, which is defined as the difference between the kinetic energy T and the sum of the potential energy of deformation (strain energy density), ${\pi }^{\mathrm{int}}$, plus the work of the external forces, defined in the previous equation. δ is the variation symbol, ρ is the material density and dots indicate differentiation with respect to time. $\tau {*}^{\alpha \beta }$are the stress components corresponding to the specified edge tractions and $\nu_\beta$ are the components of the unit normal vector with respect to the mid-surface.

By differentiating the internal strain energy and the kinetic energy and discretizing, we reach the classical eigenvalue problem depending on the stiffness matrix and the mass matrix.

4.1 Strain energy of the shell element and equilibrium equations.

In order to formulate the stiffness matrix of the shell element, we develop an expression of the internal strain energy. In general, we can write,

${\pi }^{\mathrm{int}}=\frac{1}{2}\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }^{ij}{E}_{ij}dV$ (46)

where ${\tau }^{ij}$is the symmetric stress tensor and $E_{ij}$ is the linear strain tensor.

Considering the expression (7) it can be written,

${\pi }^{\mathrm{int}}=\underset{A}{\overset{}{{\displaystyle \int }}}\underset{\raisebox{1ex}{$-h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overset{\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{{\displaystyle \int }}}\left(\frac{1}{2}{\tau }^{\alpha \beta }{E}_{\alpha \beta }+{\tau }^{\alpha 3}{E}_{\alpha 3}\right)\sqrt{a}\left(1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right)d{\alpha }^{1}d{\alpha }^2$ (47)

Introducing the constitutive equations relating strains and stresses, for a linear elastic material, we have

${\tau }^{\alpha \beta }=C^{\alpha \beta \gamma \delta }{E}_{\gamma \delta },{\tau }^{\alpha 3}=C^{\alpha 3\beta 3}{E}_{\beta 3}$ (48)

If this last expression is developed some integrals appear,

$D_0^{\alpha \beta \gamma \delta }={\displaystyle {\int }_{-\frac{t}{2}}^{\frac{t}{2}}{C}^{\alpha \beta \gamma \delta }}\mu d{\alpha }^3$ (49)

$D_1^{\alpha \beta \gamma \delta }={\displaystyle {\int }_{-\frac{t}{2}}^{\frac{t}{2}}{\alpha }^3{C}^{\alpha \beta \gamma \delta }}\mu d{\alpha }^3$ (50)

$D_2^{\alpha \beta \gamma \delta }={\displaystyle {\int }_{-\frac{t}{2}}^{\frac{t}{2}}{\left({\alpha }^3\right)}^2{C}^{\alpha \beta \gamma \delta }}\mu d{\alpha }^3$ (51)

which are generalized constitutive equations for the shells appearing in the strain energy equation.

$C^{\alpha \beta \gamma \delta }$is the known elastic constant tensor (elasticity tensor) defined as:

$C^{\alpha \beta \gamma \delta }=\frac{1-\nu }{2}\left(a^{\alpha \delta }{a}^{\beta \gamma }+a^{\alpha \gamma }{a}^{\beta \delta }+\frac{2\nu }{1-\nu }{a}^{\alpha \beta }{a}^{\gamma \delta }\right)$ (52)

where υ is the Poisson coefficient and $a^{\alpha \delta }$is the metric tensor of the reference surface.

For clarity, we can develop this last tensor equation by components (Voight notation),

$C=\left(\begin{array}{ccc}{\left(a^{11}\right)}^2& \left(1-\upsilon \right){\left(a^{12}\right)}^2+\upsilon a^{11}{a}^{22}& a^{11}{a}^{22}\\ & {\left(a^{22}\right)}^2& a^{12}{a}^{22}\\ Sym& & \frac{\left(1-\upsilon \right)}{2}{a}^{11}\cdot a^{22}+\frac{\left(1+\upsilon \right)}{2}{\left(a^{12}\right)}^2\end{array}\right)$ (53)

Full expressions for the tensors $D_1^{\alpha \beta \gamma \delta }$ and $D_2^{\alpha \beta \gamma \delta }$can be found integrating over the thickness,

$D_1^{\alpha \beta \gamma \delta }={\displaystyle {\int }_{-\frac{t}{2}}^{\frac{t}{2}}{\alpha }^3{C}^{\alpha \beta \gamma \delta }}\mu d{\alpha }^3=\frac{\text{E}}{2(1+\nu )}\left[\begin{array}{l}2\left(a^{\alpha \gamma }{b}^{\beta \delta }+b^{\alpha \gamma }{a}^{\beta \delta }\right)+2\left(a^{\alpha \delta }{b}^{\beta \gamma }+b^{\alpha \delta }{a}^{\beta \gamma }\right)+\\ \frac{4\nu }{1-\nu }\left(a^{\alpha \beta }{b}^{\gamma \delta }+b^{\alpha \beta }{a}^{\gamma \delta }\right)\end{array}\right]$ (54)

$D_2^{\alpha \beta \gamma \delta }={\displaystyle {\int }_{-\frac{t}{2}}^{\frac{t}{2}}{\left({\alpha }^3\right)}^2{C}^{\alpha \beta \gamma \delta }}\mu d{\alpha }^3=\frac{\text{E}}{2(1+\nu )}\left[\begin{array}{l}\left(3a^{\alpha \gamma }{b}_{\tau }^{\beta }{b}^{\tau \delta }+4b^{\alpha \gamma }{b}^{\beta \delta }+3a^{\beta \delta }{b}_{\tau }^{\alpha }{b}^{\tau \gamma }\right)+\\ \left(3a^{\alpha \delta }{b}_{\tau }^{\beta }{b}^{\tau \gamma }+4b^{\alpha \delta }{b}^{\beta \gamma }+3a^{\beta \gamma }{b}_{\tau }^{\alpha }{b}^{\tau \delta }\right)\\ +\frac{2\nu }{1-\nu }\left(3a^{\alpha \beta }{b}_{\tau }^{\gamma }{b}^{\tau \delta }+4b^{\alpha \beta }{b}^{\gamma \delta }+3a^{\gamma \delta }{b}_{\tau }^{\alpha }{b}^{\tau \beta }\right)\end{array}\right]$ (55)

Assuming that the shell is moderately thick, the above equation can be rewritten as,

${\pi }^{\mathrm{int}}\cong \underset{A}{\overset{}{{\displaystyle \int }}}\underset{\raisebox{1ex}{$-h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\overset{\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{{\displaystyle \int }}}\left(\frac{1}{2}{\tau }^{\alpha \beta }{E}_{\alpha \beta }+{\tau }^{\alpha 3}{E}_{\alpha 3}\right)\sqrt{a}\mu d{\alpha }^{1}d{\alpha }^2$ (56)

If we define generalized stresses in the form,

$\begin{array}{l}{n}^{\alpha \beta }={\displaystyle \int \mu {\mu }_{\alpha }^{\beta }{\tau }^{\alpha \beta }d{\alpha }^3={\displaystyle \int \left[1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right]}\left[{\delta }_{\alpha }^{\beta }-{\alpha }^3{b}_{\alpha }^{\beta }\right]}{\tau }^{\alpha \beta }d{\alpha }^3=D_0^{\alpha \beta \gamma \delta }{e}_{\alpha \beta }+D_1^{\alpha \beta \gamma \delta }{\rho }_{\alpha \beta }\\ +D_2^{\alpha \beta \gamma \delta }{\kappa }_{\alpha \beta }\approx h{C}^{\alpha \beta \gamma \delta }{\epsilon }_{\gamma \delta }\end{array}$ (57)

$\begin{array}{l}{}^1{m}^{\alpha \beta }={\displaystyle \int \mu {\mu }_{\alpha }^{\beta }{\tau }^{\alpha \beta }{\alpha }^{3}d{\alpha }^3={\displaystyle \int \left[1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right]}\left[{\delta }_{\alpha }^{\beta }-{\alpha }^3{b}_{\alpha }^{\beta }\right]}{\tau }^{\alpha \beta }{\alpha }^{3}d{\alpha }^3=D_1^{\alpha \beta \gamma \delta }{e}_{\alpha \beta }+\\ D_2^{\alpha \beta \gamma \delta }{\rho }_{\alpha \beta }\approx \frac{h^3}{12}{C}^{\alpha \beta \gamma \delta }{\rho }_{\gamma \delta }\end{array}$ (58)

$\begin{array}{l}{}^2{m}^{\alpha \beta }={\displaystyle \int \mu {\mu }_{\alpha }^{\beta }{\tau }^{\alpha \beta }{\alpha }^{3}d{\alpha }^3={\displaystyle \int \left[1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right]}\left[{\delta }_{\alpha }^{\beta }-{\alpha }^3{b}_{\alpha }^{\beta }\right]}{\tau }^{\alpha \beta }{\left({\alpha }^3\right)}^{2}d{\alpha }^3\\ =D_2^{\alpha \beta \gamma \delta }{e}_{\alpha \beta }\approx \frac{h^5}{80}{C}^{\alpha \beta \gamma \delta }{\kappa }_{\gamma \delta }\end{array}$ (59)

$q^{\alpha }={\displaystyle \int \mu {\tau }^{\alpha 3}d{\alpha }^3={\displaystyle \int \left[1-2H{\alpha }^3+K{\left({\alpha }^3\right)}^2\right]}}{\tau }^{\alpha 3}d{\alpha }^3\approx h{C}^{\alpha 3\beta 3}{E}_{\beta 3}$ (60)

the elastic strain energy of the shell element can be written

${\pi }^{\mathrm{int}}=\frac{1}{2}\underset{A}{\overset{}{{\displaystyle \int }}}\left(C^{\alpha \beta \gamma \delta }h{e}_{\alpha \beta }{e}_{\gamma \delta }+\frac{h^3}{12}{C}^{\alpha \beta \gamma \delta }{\rho }_{\alpha \beta }{\rho }_{\gamma \delta }+\frac{h^5}{80}{C}^{\alpha \beta \gamma \delta }{\kappa }_{\alpha \beta }{\kappa }_{\gamma \delta }+C^{\alpha 3\beta 3}h{E}_{\alpha 3}{E}_{\beta 3}\right)\sqrt{a}d{\alpha }^{1}d{\alpha }^2$ (61)

where the expressions of the various tensors were derived in the previous section and are given by equations (34) and (36).

In the theory of moderately thick shells, the tensor ${\kappa }_{\alpha \beta }$is usually not taken into account in the derivation of the stiffness matrix. Note that the contribution in the strain energy of this term is multiplied by $\frac{h^5}{80}$and so we understand is much less than the other terms. Bearing this in mind and neglecting this term, we can differentiate this last equation,

$\begin{array}{l}\delta {\pi }^{\mathrm{int}}={\displaystyle \int {\displaystyle {\int }_{\Omega }\left(n^{\alpha \beta }\delta {\epsilon }_{\alpha \beta }+m^{\alpha \beta }\delta {\rho }_{\alpha \beta }+q^{\alpha }\delta {\varsigma }_{\alpha }\right)}}\sqrt{a}d{\alpha }^{1}d{\alpha }^2={\displaystyle \int {\displaystyle {\int }_{\Omega }\left(n^{\alpha \beta }\left(\frac{1}{2}\left(\delta u_{\alpha |\beta }+\delta u_{\beta |\alpha }\right)-b_{\alpha \beta }\delta u_3\right)\right)}}\\ {+}^1{m}^{\alpha \beta }\left(\frac{1}{2}\left(\delta {\theta }_{\alpha |\beta }+\delta {\theta }_{\beta |\alpha }-b_{\beta }^{\lambda }\delta u_{\lambda |\alpha }-b_{\alpha }^{\lambda }\delta u_{\lambda |\beta }\right)+c_{\alpha \beta }\delta u_3\right)+q^{\alpha }\left(\frac{1}{2}\left(\delta {\theta }_{\alpha }+\delta u_{\mathrm{3,}\alpha }+b_{\alpha }^{\lambda }\delta u_{\lambda }\right)\right)\sqrt{a}d{\alpha }^{1}d{\alpha }^2\end{array}$ (62)

The External virtual work by distributed load intensity $\delta W_{\Omega }$ and the external virtual work by shell edge resultant loads $\delta W_{\Gamma }$ are given as:

$\delta W_{\Omega }={\displaystyle {\int }_{\Omega }\left(p^{\alpha }\delta u_{\alpha }+p^3\delta u_3\right)\sqrt{a}d{\alpha }^{1}d{\alpha }^2}$ (63)

$\delta W_{\Gamma }={\displaystyle {\int }_{\Gamma }\left(n^{\alpha \beta }{n}_{\alpha }\delta u_{\beta }+m^{\alpha \beta }{n}_{\alpha }\delta {\theta }_{\beta }+q^{\alpha }{n}_{\alpha }\delta u_3\right)d\Gamma }$ (64)

The finite element discretization of these functionals will be explained later, employing the MITC technique. Continuing with the usual process, the equilibrium equations can be deduced through the variational process:

$n^{\alpha \beta }|{}_{\beta }-\left(b_{\gamma }^{\alpha }{m}^{\gamma \beta }\right)|{}_{\alpha }-b_{\gamma }^{\beta }{q}^{\gamma }+p^{\beta }=\rho u_{\alpha ,tt}$ (65)

$-b_{\gamma \sigma }{n}^{\sigma \gamma }+b_{\gamma \sigma }{b}_{\delta }^{\sigma }{m}^{\delta \gamma }+q^{\alpha }|{}_{\alpha }+p^3=\rho w_{,tt}$ (66)

$m^{\alpha \beta }|{}_{\alpha }-q^{\beta }+c^{\rho }=J{\theta }_{\alpha ,tt}$ (67)

where $J=\frac{\rho h^3}{12}$is the rotational inertia. With the consistent boundary conditions,

$n^{*}{}^{\alpha \beta }{\nu }_{\beta }+b_{\gamma }^{\alpha }{m}^{*\delta \gamma }{\nu }_{\alpha }=N_{st}{}^{\alpha }$ (68)

$q^{*}{}^{\alpha }{\nu }_{\alpha }=Q_{sh}$ (69)

$m^{*\delta \gamma }{\nu }_{\gamma }{\nu }_{\beta }=M_{be}$ (70)

$m^{*\delta \gamma }{t}_{\lambda }{\nu }_{\beta }=H_{tw}$ (71)

Where $M_{be}$is the bending moment resultant and $H_{tw}$ is the twisting moment resultants. $N_{st}$ and $Q_{sh}$ are the in-surface stretching resultant and shear resultant, respectively. The terms ${\nu }_{\gamma },{\nu }_{\beta }$ are the components of the unit normal vector with respect to the mid-surface, and $t_{\lambda }$are the components of the unit tangent vector with respect to the same mid-surface.

4.2 Kinetic energy of the shell finite element in general curvilinear coordinates. Consistent and Lumped mass matrices.

The mass matrix of the finite element is necessary to study the dynamic problem. The lumped and consistent mass matrices are alternative formulations which can be used.

In this case, we propose to develop the consistent mass matrix, using the same shape functions for its discretization than those used for the geometry. As low order elements will be used for the interpolation of the variables, the solutions with the lumped mass matrix are very poor, and that is why we have opted for the consistent mass matrix that most authors agree in emphasizing for its elders benefits, Zienkiewicz and Taylor (2000).

The development of the mass matrix for a shell finite element in general curvilinear coordinates is more complex than for other elements such as three-dimensional elements, but conceptually it is developed in the same way, through the variation of the kinetic energy.

Our starting point is the definition of the kinetic energy,

$dT=\underset{A}{\overset{}{{\displaystyle \int }}}\underset{{\alpha }^3}{\overset{}{{\displaystyle \int }}}\ddot{\overline{R}}\rho ·\delta R\mu d{\alpha }^{3}dA$ (72)

where R is the position vector of an arbitrary shell point in the deformed configuration. If we consider that the underlying shell theory allows for the transverse shear deformation effects, we can write,

$dR=d{R}_0+{\alpha }^{3}dN$ (73)

$\ddot{R}={\ddot{R}}_0+{\alpha }^{3}N$ (74)

Therefore the kinetic energy functional is:

$\delta T=\underset{\text{A},\text{el}}{\overset{}{{\displaystyle \int }}}\left({\ddot{R}}_0\delta R_0\right)\underset{\begin{array}{c}\text{ConstantTerm}\\ \text{C}\text{.T}\text{.}\end{array}}{\underset{{\alpha }^3}{\overset{}{{\displaystyle \int }}}\rho \mu \text{d}{\alpha }^3}+\left(\partial \ddot{N}{R}_0+\partial {\ddot{R}}_{0}N\right)\underset{\begin{array}{c}\text{Linear Term}\\ \text{L}\text{.T}\text{.}\end{array}}{\underset{}{\underset{{\alpha }^3}{\overset{}{{\displaystyle \int }}}}\rho {\alpha }^3\mu \text{d}{\alpha }^3}+\left(\partial \ddot{N}N\right)\underset{\begin{array}{c}\text{CuadraticTerm}\\ \text{C}\text{.T}\text{.}\end{array}}{\underset{}{\underset{{\alpha }^3}{\overset{}{{\displaystyle \int }}}}\rho {\left({\alpha }^3\right)}^2\mu \text{d}{\alpha }^3}$ (75)

where C.T., L.T., and Q.T. stand for constant term, linear term, and quadratic term, respectively. We now carry out the interpolation of the variables involved:

$\delta {\ddot{R}}_0={\displaystyle \sum }_{i=1}^4{N}_i\delta R_0,\delta {\ddot{R}}_0={\displaystyle \sum }_{i=1}^4{N}_i\delta {\ddot{R}}_0,\delta N={\displaystyle \sum }_{i=1}^4{N}_i\delta N_{}^{},\delta {\ddot{N}}_{}={\displaystyle \sum }_{i=1}^4{N}_i\delta {\ddot{N}}_{}^{}$ (76)

Introducing these expressions in (75),

$\begin{array}{l}\delta T={\displaystyle \sum _{n=1}^{N^{el}}{\displaystyle \sum _{m=1}^{N^{el}}\left(\delta {\ddot{R}}_0\cdot \delta R_0\right)}}\underset{\begin{array}{c}\text{ConstantTerm}\\ \text{C}\text{.T}\text{.}\end{array}}{\underset{A_{elem}}{\overset{}{{\displaystyle \int }}}(C.T.)N_n{N}_m\text{d}{A}_{elem}}+{\displaystyle \sum _{n=1}^{N^{el}}{\displaystyle \sum _{m=1}^{N^{el}}\left(\delta {\ddot{R}}_0\cdot \delta N\right)}}\underset{\begin{array}{c}\text{Linear Term}\\ \text{L}\text{.T}\text{.}\end{array}}{\underset{}{\underset{A_{elem}}{\overset{}{{\displaystyle \int }}}(L.T.)}{N}_n{N}_m\text{d}{A}_{elem}}\\ +{\displaystyle \sum _{n=1}^{N^{el}}{\displaystyle \sum _{m=1}^{N^{el}}\left(\delta \ddot{N}\cdot \delta R_0\right)}}\underset{\begin{array}{c}\text{Linear Term}\\ \text{C}\text{.T}\text{.}\end{array}}{\underset{}{\underset{A_{elem}}{\overset{}{{\displaystyle \int }}}}(L.T.)N_n{N}_m\text{d}{A}_{elem}}+{\displaystyle \sum _{n=1}^{N^{el}}{\displaystyle \sum _{m=1}^{N^{el}}\left(\delta \ddot{N}\cdot \delta N\right)}}\underset{\begin{array}{c}\text{CuadraticTerm}\\ \text{C}\text{.T}\text{.}\end{array}}{\underset{}{\underset{A_{elem}}{\overset{}{{\displaystyle \int }}}}(CuadraticT.)N_n{N}_m\text{d}{A}_{elem}}\end{array}$ (77)

The terms in equation (75) can be integrated without difficulty in the light of equation (7). Doing this, we find a matrix expression for the kinetic energy depending on the mass matrix, the nodal variables of the element and its second derivatives. The explicit expression of the mass matrix, M, is shown below.

$M=\left(\begin{array}{ccccc}{a}^{11}\left(h+\frac{1}{12}{h}^{3}K\right)& a^{12}\left(h+\frac{1}{12}{h}^{3}K\right)& 0& a^{11}\left(-\frac{1}{6}H{h}^3\right)& a^{12}\left(-\frac{1}{6}H{h}^3\right)\\ & a^{22}\left(h+\frac{1}{12}{h}^{3}K\right)& 0& a^{12}\left(-\frac{1}{6}H{h}^3\right)& a^{22}\left(-\frac{1}{6}H{h}^3\right)\\ & & \left(h+\frac{1}{12}{h}^{3}K\right)& 0& 0\\ & & & \left(\frac{h^3}{12}+\frac{1}{80}{h}^{5}K\right)a^{11}& \left(\frac{h^3}{12}+\frac{1}{80}{h}^{5}K\right)a^{12}\\ Sym& & & & \left(\frac{h^3}{12}+\frac{1}{80}{h}^{5}K\right)a^{22}\end{array}\right)$ (78)

Note that the terms of the mass matrix depend on the mean and Gaussian curvature (H and K) and on the thickness of the shell, in contrast, for example, to a plate element, where all these quantities except the thickness are 0. If we introduce the Kirchhoff-Love assumption given by equation (43), a more refined expression for the mass matrix for thin shells can be found. At this point, we can solve equation (44).

5 LOCKING OF THE STIFFNESS MATRIX. MIXED INTERPOLATION OF TENSORIAL COMPONENTS

Once an expression for the internal energy of the shell element has been deduced, with the inclusion of the shear strain effects and the kinetic energy necessary for the dynamic formulation of the problem, we need to discretize the middle surface of the shell by the isoparametric formulation, to subsequently integrate over the thickness.

For the interpolation of the displacements u within an element, bilinear shape functions have been used,

$u={\displaystyle \sum }_{i=1}^4{N}_i\left({\xi }_1,{\xi }_2\right){\delta }_i$ (79)

where ${\delta }_i$ is the generalized nodal displacement vector, i.e. ${\delta }_i=(u_1,u_2,u_3,{\theta }_1,{\theta }_2)$and ${\xi }_1,{\xi }_2$are the isoparametric coordinates.

One of the advantages of the proposed formulation is that the basis vectors and other geometric parameters of the surface are known exactly.

As mentioned in the introduction, when the thickness of the shell becomes thinner, the shell elements suffer from the so-called shear locking (and membrane locking if the element are curved), causing spurious shear and membrane strains in bending dominated problems.

To circumvent this drawback, we turn to a well known technique originally proposed by Dvorkin and Bathe (1984); Bathe and Dvorkin (1986) for elements of low order, and Bucalem and Bathe (1993) for higher order elements, called mixed interpolation of tensorial components ($MIT{C}_n$, where n stands for the number of nodes in the element), which can be encompassed within assumed strain methods. In general, the covariant strains are individually interpolated in $MIT{C}_n$finite elements, where a set of $n_{ij}$tying points associated to each strain component $E_{ij}$ is adopted.

Hence, the shear components of the linear strain tensor (in covariant form) are interpolated so that

${\tilde{E}}_{\alpha \beta }^{AS}={\displaystyle \sum }_{k=1}^p{\Lambda }_k^{\alpha \beta }\left({\xi }_1,{\xi }_2\right)E_{\alpha \beta }^{DI}\left({\xi }_{1k},{\xi }_{2k},t\right)$ (80)

In this equation, the assumed strains ${\tilde{E}}_{\alpha \beta }^{AS}$ can be computed from the strains which have been interpolated through the shape functions ${\Lambda }_k^{\alpha \beta }\left({\xi }_1,{\xi }_2\right)$ . These shape functions depend on the isoparametric coordinates $\left({\xi }_1,{\xi }_2\right)$ and correspond to the strain component (α,β) evaluated at point k. The so-called tying points defined on the shell midsurface with coordinates $\left({\xi }_P,{\xi }_P\right)$define the assumed covariant strain components. p is the number of tying points for the covariant strain component ${\tilde{E}}_{\alpha \beta }^{AS}$and ${\Lambda }_k^{\alpha \beta }\left({\xi }_1,{\xi }_2\right)$are the assumed interpolation functions, satisfying,

${{\Lambda }_k^{\alpha \beta }|}_l={\Lambda }_k^{\alpha \beta }\left({\xi }_l,{\xi }_l\right)={\delta }_{kl}$ (81)

We must note that this procedure is carried out at the elemental level for each individual element.

For more details regarding this formulation, we refer to the above-mentioned papers. In our case, the components of the linear strain tensor are given by equations (34) and (36). If we develop these tensor equations, the relation between strains and displacements can be written as,

$\begin{array}{l}\\ B=\left(\begin{array}{ccccc}{\partial }_1-{\Gamma }_{11}^1& -{\Gamma }_{11}^2& -b_{11}& 0& 0\\ -{\Gamma }_{22}^1& {\partial }_2-{\Gamma }_{22}^2& -b_{22}& 0& 0\\ {\partial }_2-2{\Gamma }_{12}^1& {\partial }_1-2{\Gamma }_{12}^2& -2b_{12}& 0& 0\\ -b_1^1{\partial }_1+b_1^{\alpha }{\Gamma }_{\alpha 1}^1& -b_1^2{\partial }_1+b_1^{\alpha }{\Gamma }_{\alpha 1}^2& b_1^{\alpha }{b}_{\alpha 1}& {\partial }_1-{\Gamma }_{11}^1& -{\Gamma }_{11}^2\\ -b_2^1{\partial }_2+b_2^{\alpha }{\Gamma }_{\alpha 2}^1& -b_2^2{\partial }_2+b_2^{\alpha }{\Gamma }_{\alpha 2}^2& b_2^{\alpha }{b}_{\alpha 2}& -{\Gamma }_{22}^1& {\partial }_2-{\Gamma }_{22}^2\\ \left(\begin{array}{l}-b_2^1{\partial }_1+b_2^{\alpha }{\Gamma }_{\alpha 1}^1\\ -b_1^1{\partial }_2+b_1^{\alpha }{\Gamma }_{\alpha 2}^2\end{array}\right)& \left(\begin{array}{l}-b_2^2{\partial }_1+b_2^{\alpha }{\Gamma }_{\alpha 1}^2\\ -b_1^2{\partial }_2+b_1^{\alpha }{\Gamma }_{\alpha 2}^2\end{array}\right)& 2b_1^{\alpha }{b}_{\alpha 2}& {\partial }_2-2{\Gamma }_{12}^1& {\partial }_1-2{\Gamma }_{12}^2\\ b_2^1& b_2^2& {\partial }_2& 0& 1\\ b_1^1& b_1^2& {\partial }_1& 1& 0\end{array}\right)\end{array}$ (82)

with

$\left(\begin{array}{c}{e}_{11}\\ e_{22}\\ e_{12}\\ {\rho }_{11}\\ {\rho }_{22}\\ {\rho }_{12}\\ {\zeta }_{23}\\ {\zeta }_{13}\end{array}\right)=B\left(\begin{array}{c}{u}_1\\ u_2\\ w\\ {\theta }_1\\ {\theta }_2\end{array}\right)$ (83)

which would be properly interpolated by the chosen shape functions. For the construction of the stiffness matrix, note that we use ${\tilde{E}}_{\alpha \beta }^{AS}$instead of $E_{\alpha \beta }^{DI}$defined by the last equation. Also, it is important to remark that if we had considered higher order moments defined by the tensor ${\kappa }_{\alpha \beta }$, this last matrix would have to contain its components adequately ordered.

$\left(\begin{array}{c}{\kappa }_{11}\\ {\kappa }_{22}\\ {\kappa }_{12}\end{array}\right)=\left(\begin{array}{cc}-b_1^1{\partial }_1+b_1^{\alpha }{\Gamma }_{\alpha 1}^1& -b_1^2{\partial }_1+b_1^{\alpha }{\Gamma }_{\alpha 1}^2\\ -b_2^1{N}_{\mathrm{,2}}^I+b_2^{\alpha }{\Gamma }_{\alpha 2}^1& -b_2^2{\partial }_2+b_2^{\alpha }{\Gamma }_{\alpha 2}^2\\ \left(\begin{array}{c}-b_2^1{\partial }_1+b_2^{\alpha }{\Gamma }_{\alpha 1}^1\\ -b_1^1{\partial }_2+b_1^{\alpha }{\Gamma }_{\alpha 2}^2\end{array}\right)& \begin{array}{c}-b_2^2{\partial }_1+b_2^{\alpha }{\Gamma }_{\alpha 1}^2\\ -b_1^2{\partial }_2+b_1^{\alpha }{\Gamma }_{\alpha 2}^2\end{array}\end{array}\right)\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right)$ (84)

It should be noted that the geometric characteristics of the surface are defined only at the nodes of the finite element, and that due to the invariance of the internal energy we can operate directly on the tensorial components without having to resort to the physical components of the variables.

6 RESULTS AND DISCUSSION

To test the validity of these results, we compared the results from using the new shell element with the results obtained by other authors with several types of shell structures: hyperbolic paraboloid shells as well as cylindrical, spherical and elliptic parabolioid shells. In general, doubly-curved shells pose very stringent tests to validate solutions, whether they are static or dynamic. A considerable number of examples, including clamped and simply supported isotropic shells with a wide range of ratios h/R and L/R, will be investigated.

The comparison of the present formulation is made with regard to the original test case, and also with the four-node shell model based on the degenerated solid approach developed by Dvorkin and Bathe (1984). In all cases, the consistent mass matrices have been used, and finite element models have been discretized with 50 by 50 elements in the x and y directions respectively. Throughout these examples, we show the efficiency and versatility of the present formulation.

6.1 Example 1: Corner point supported thin hyperbolic paraboloid

To begin with, Narita and Leissa (1984) and Chakravorty and Bandyopadhyay (1995) and Chakravorty et al. (1995) study the free vibration of several types of surfaces incluiding the hyperbolic paraboloid. The hyperbolic paraboloid test, a surface with negative Gaussian curvature, is quite demanding both surface type (doubly curved shell) and the boundary conditions.

The corner point supported thin hyperbolic paraboloid test of Figure 2 involves a surface with a negative Gaussian curvature, and is quite demanding both as to the type of surface involved (doubly curved shell) and its boundary conditions. The data of the problem is: square planform 1 x 1 m side, constant thickness of h=0.01 m., equal radii of curvature and opposite in sign $R_y=R_x=2m.$, Young modulus, $E=10.92·{10}^{6}N/m^2$ Poisson's ratio $\upsilon =\mathrm{0,3}$ and density $\rho =\text{1}00\text{kg}/{\text{m}}^3$.

Figure 2:
Corner point supported thin hyperbolic paraboloid with curved edges (Example 1, Table 1).

Table 1 presents the frequency associated with the first mode of vibration of the hyperbolic paraboloid according to results obtained by the above-named authors.

The frequency associated with the first mode of vibration of the hyperbolic paraboloid according the new formulation is 17, 21rad / s. The results obtained by the above-named authors are:

Table 1:
Frequency associated with the first mode of vibration of the corner point supported thin hyperbolic paraboloid with curved edges. Comparison with Narita and Leissa (1984), Chakravorty and Bandyopadhyay (1995), and the degenerated solid approach of Dvorkin and Bathe (1984).

Narita and Leissa (1984 )	17.16 rad/s
Chakravorty and Bandyopadhyay (1995)	17.25 rad/s
Deg. Sol. Shell	17.31 rad/s
Proposed	17.21 rad/s

The first three modes of vibration for the corner point supported thin hyperbolic paraboloid with curved edges are depicted in Figure 3:

Figure 3:
First (a), second (b), and third (c) mode of vibration for the thin hyperbolic paraboloid simply supported at the corners (Example 1, Table 1).

6.2 Example 2: Simply supported thin hyperbolic paraboloid shell with square planform

In order to compare the higher order modes of vibration, we compare the results of the present formulation with Liew and Lim (1996) and the degenerated solid approach of Dvorkin and Bathe (1984).

Figure 4:
Simply supported thin hyperbolic paraboloid shell with square planform (Example 2,Table 2).

For a simply supported thin hyperbolic paraboloid shell with square planform (Figure 4), the results are compared for the first, fourth and eight vibration modes given in this work. In this case we take ν=0.3, a⁄b=1, and b⁄h=100. The results are shown in Table 2 for the non dimensional frequency parameter λ given by

$\lambda =\omega ab\sqrt{\frac{\rho h}{D}},D=\frac{E{h}^3}{12\left(1-{\upsilon }^2\right)}$ (85)

Table 2:
Frequency parameter λ for simply supported thin hyperbolic paraboloid shells and various values of $\frac{b}{R_y}$and $\frac{R_y}{R_x}$.

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Liew and Lim [1996]	42.688	83.136	130.86	-1	0.1
Deg. Sol. Shell	43.326	84.379	135.40	-1
Proposed	42.681	84.062	133.80	-1
Liew and Lim [1996]	36.713	81.902	130.71	-0.5	0.1
Deg. Sol. Shell	36.012	83.942	131.46	-0.5
Proposed	36.715	82.887	133.52	-0.5
Liew and Lim [1996]	104.02	110.78	150.76	-1	0.3
Deg. Sol. Shell	103.41	113.06	154.12	-1
Proposed	103.42	110.79	151.22	-1
Liew and Lim [1996]	75.152	102.10	148.35	-0.5	0.3
Deg. Sol. Shell	74.965	100.89	147.13	-0.5
Proposed	75.082	102.25	148.99	-0.5
Liew and Lim [1996]	148.74	159.39	198.54	-1	0.5
Deg. Sol. Shell	146.39	161.26	198.25	-1
Proposed	145.62	160.34	196.79	-1
Liew and Lim [1996]	105.48	136.48	177.97	-0.5	O.5
Deg. Sol. Shell	105.71	135.93	178.74
Proposed	104.66	137.11	178.99

• (ν=0.3,, a⁄b=1, and b⁄h=100). Comparison with Liew and Lim (1996) and the degenerated solid approach of Dvorkin and Bathe (1984).

6.2.1 Mesh convergence analysis for thin hyperbolic paraboloid shell

A mesh convergence analysis has been performed to evaluate the performance of the proposed formulation. Table 3 presents the frequency parameter λ for simply supported thin hyperbolic paraboloid shells.

Table 3:
Mesh convergence analysis of the proposed formulation. Parameter for simply supported thin hyperbolic paraboloid shells and various values of $\frac{b}{R_y}$and $\frac{R_y}{R_x}$. (ν=0.3, a⁄b=1, and b⁄h=100)

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Present (50*50)	42.681	84.062	133.80	-1	0.1
(40*40)	42.681	84.071	133.81
(30*30)	42.683	84.075	133.84
(20*20)	42.686	84.082	133.89
(10*10)	42.697	84.106	134.59
(5*5)	43.041	85.127	136.24
Present	36.715	82.887	133.52	-0.5	0.1
(40*40)	36.715	82.892	133.54
(30*30)	36.718	82.898	133.55
(20*20)	36.720	82.901	133.56
(10*10)	36.735	83.047	134.72
(5*5)	37.182	83.966	135.93
Present	103.42	110.79	151.22	-1	0.3
(40*40)	103.42	110.80	151.23
(30*30)	103.43	110.81	151.25
(20*20)	103.43	110.82	151.30
(10*10)	103.78	111.24	152.09
(5*5)	104.59	112.36	154.00
Present	75.082	102.25	148.99	-0.5	0.3
(40*40)	75.082	102.28	149.02
(30*30)	75.084	102.29	149.03
(20*20)	75.087	102.31	149.03
(10*10)	75.580	102.82	149.96
(5*5)	75.894	103.56	151.86
Present	145.62	160.34	196.79	-1	0.5
(40*40)	145.62	160.36	196.80
(30*30)	145.63	160.36	196.81
(20*20)	145.63	160.38	196.82
(10*10)	145.99	160.42	197.59
(5*5)	146.46	161.68	199.33
Present	104.66	137.11	178.99	-0.5	0.5
(40*40)	104.66	137.12	179.02
(30*30)	104.67	137.12	179.02
(20*20)	104.68	137.14	179.04
(10*10)	104.71	137.20	180.47
(5*5)	105.34	138.09	182.90

6.2.2 Vibration analysis for thin hyperbolic paraboloid shell with distorted meshes

To evaluate the proposed formulation for distorted meshes, the structured meshes (consisting in elements of size L x L) where distorted by different fractions of the element length L, for values 0.1L, 0.2L, and 0.4L, see Figure 5. Table 4 presents the frequency parameter λ for simply supported thin hyperbolic paraboloid shells with distorted mesh.

Figure 5:
Test with distorted meshes. Simply supported thin hyperbolic paraboloid shells, with mesh distorted by different fractions of the element length: (a) 0.1L, (b) 0.2L, and (c) 0.4L. (Table 4).

Table 4:
Vibration analysis for distorted meshes, for the simply supported thin hyperbolic paraboloid shells. Parameter λ for the first frequency mode (ν=0.3, a⁄b=1, and b⁄h=100). 50 x 50 element mesh distorted with $dis{t}_{frac}$ of 0.1L, 0.2L, and 0.4L. Comparison with the deg. solid approach of Dvorkin and Bathe (1984).

	0.1L	0.2L	0.4L	Undistorted mesh	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Proposed finite element	42.681	42.681	42.682	42.681	-1	0.1
Deg. Solid Shell	43.329	43.704	43.981	43.326	-1
Proposed finite element	36.715	36.715	36.717	36.715	-0.5	0.1
Deg. Solid Shell	36.014	36.126	36.532	36.012	-0.5
Proposed finite element	103.42	103.42	103.43	103.42	-1	0.3
Deg. Solid Shell	103.42	103.49	103.69	103.41	-1

6.3 Example 3: Completely clamped thin hyperbolic paraboloid shell with square planform

Table 5 present the results for other boundary conditions, i.e., a completely clamped thin hyperbolic paraboloid (ν=0.3,, a⁄b=1, and b⁄h=100), see Figure 6.

Figure 6:
Clamped thin hyperbolic paraboloid shell with square planform (Example 3, Table 5). 50 by 50 finite element discretization.

The first, fourth and eighth modes of vibration for the case of completely clamped thin hyperbolic paraboloid with $\frac{R_y}{R_x}=-1$ and $\frac{b}{R_y}=0.5$are depicted in Figure 7. Given these results, the efficiency of the propose shell element is of the order of other shell elements (less than 1% error) in all cases.

Table 5:
Frequency parameter $\lambda$ for completely clamped thin hyperbolic paraboloid shells and various values of $\frac{b}{R_y}$and $\frac{R_y}{R_x}$. (ν=0.3, a⁄b=1, b⁄h=100). Comparison with the degenerated solid approach of Dvorkin and Bathe (1984).

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Proposed	50.92	112.99	173.32	-1	0.1
Deg. Sol. Shell	50.73	111.57	170.84	-1
Proposed	46.49	112.74	173.24	-0.5	0.1
Deg. Sol. Shell	46.72	109.71	169.90	-0.5
Proposed	110.93	128.60	184.18	-1	0.3
Deg. Sol. Shell	110.89	127.61	181.31	-1
Proposed	90.88	124.62	184.90	-0.5	0.3
Deg. Sol. Shell	89.02	127.57	188.13	-0.5
Proposed	154.83	167.23	207.95	-1	0.5
Deg. Sol. Shell	154.22	167.07	205.86	-1
Proposed	113.20	153.80	204.83	-0.5	0.5
Deg. Sol. Shell	110.79	154.53	199.97	-0.5
Proposed	197.86	222.42	272.69	-1	0.8
Deg. Sol. Shell	198.62	221.81	272.84	-1
Proposed	150.13	184.31	248.16	-0.5	0.8
Deg. Sol. Shell	151.18	181.75	244.12	-0.5
Proposed	222.84	242.29	312.08	-1	1
Deg. Sol. Shell	223.52	241.52	311.42	-1
Proposed	173.65	198.07	279.37	-0.5	1
Deg. Sol. Shell	172.89	194.83	274.63	-0.5

Figure 7:
First (a), fourth (b), and eighth (c) modes of vibration for the case of completely clamped thin hyperbolic paraboloid with $\frac{R_y}{R_x}=-1$, and $\frac{b}{R_y}=0.5$ (Example 3, Table 5).

6.4 Example 4: Completely clamped thin spherical, cylindrical, hyperbolic paraboloid and elliptic paraboloid shells with circular planform.

This example compares the first, fourth and eight vibration modes for clamped thin shells with circular planform. Several shell geometries are evaluated $\frac{R_y}{R_x}=1$ is a spherical shell, $\frac{R_y}{R_x}\prec 0$ is a hyperbolic paraboloid shell, $\frac{R_y}{R_x}=0$is a cylindrical shell, and $0\prec \frac{R_y}{R_x}\prec 1$is an elliptic paraboloid shell.

By taking ν=0.3, a⁄b=1, b⁄h=100, the results are shown in Table 6 for the non-dimensional frequency parameter $\lambda$defined previously.

Table 6:
Parameter $\lambda$ for completely clamped thin shells with circular planform (Example 4): $R_y=R_x=1$ is a spherical, $R_y=R_x\prec 0$ is hyperbolic paraboloid, $R_y/R_x=0$is cylindrical, and $0\prec R_y/R_x\prec 1$ is elliptic paraboloid. Comparison with Dvorkin and Bathe (1984).

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Proposed	54.04	140.77	204.93	-1	0.1
Deg. Sol. Shell	54.25	139.96	205.31	-1
Proposed	50.06	140.28	204.30	-0.5	0.1
Deg. Sol. Shell	51.03	140.01	205.96	-0.5
Present	50.16	140.54	204.25	0	0.1
Deg. Sol. Shell	48.97	141.65	204.78	0
Present	54.34	141.55	204.78	0.5	0.1
Deg. Sol. Shell	53.77	140.86	204.21	0.5
Present	61.73	143.30	205.93	1	0.1
Deg. Sol. Shell	60.43	144.21	205.05	1
Present	111.80	153.72	217.54	-1	0.3
Deg. Sol. Shell	111.71	152.69	216.88	-1
Present	94.025	149.69	212.98	-0.5	0.3
Deg. Sol. Shell	93.71	150.11	213.84	-0.5
Present	93.48	151.71	212.52	0	0.3
Deg. Sol. Shell	92.91	152.25	213.07	0
Present	108.59	159.70	217.36	0.5	0.3
Deg. Sol. Shell	107.82	158.01	216.39	0.5
Present	135.54	172.99	226.53	1	0.3
Deg. Sol. Shell	133.85	174.29	225.77	1
Present	164.31	175.38	238.89	-1	0.5
Deg. Sol. Shell	162.02	177.13	239.52	-1
Present	121.86	165.70	229.12	-0.5	0.5
Deg. Sol. Shell	120.37	163.94	228.63	-0.5
Present	108.74	170.65	226.87	0	0.5
Deg. Sol. Shell	106.39	169.20	225.74	0
Present	139.49	188.52	238.80	0.5	0.5
Deg. Sol. Shell	137.30	190.80	240.07	0.5
Present	192.79	218.40	262.66	1	0.5
Deg. Sol. Shell	190.51	217.65	262.04	1

6.5 Example 5: Completely clamped moderately thick hyperbolic paraboloid shells

In order to compare the results with moderately thick shells, we use the parameter λ´, given by

$\lambda ´=\omega ab\sqrt{\frac{\rho h}{E}}$ (86)

If we analyze again the completely clamped hyperbolic paraboloid but with a relation b⁄h=10, then within the range of moderately thick shells, (ν=0.3, a⁄b=1, b⁄h=10), we obtain the results presented in Table 7.

Table 7:
Frequency parameter $\lambda ´$ for completely clamped moderately thick hyperbolic paraboloid shells and various values of $\frac{b}{Ry}$and $\frac{Ry}{Rx}$ (Example 5), (ν=0.3, a⁄b=1, b⁄h=10). Comparison with the degenerated solid approach of Dvorkin and Bathe (1984).

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Present	0.997	2.678	3.732	-1	0.1
Deg. Sol. Shell	0.992	2.648	3.709	-1
Present	0.996	2.679	3.732	-0.5	0.1
Deg. Sol. Shell	0.986	2.649	3.713	-0.5
Present	1.035	2.663	3.728	-1	0.3
Deg. Sol. Shell	1.003	2.697	3.717	-1
Present	1.019	2.666	3.730	-0.5	0.3
Deg. Sol. Shell	1.021	2.687	3.734	-0.5
Present	1.111	2.632	3.720	-1	0.5
Deg. Sol. Shell	1.075	2.606	3.678	-1
Present	1.068	2.651	3.723	-0.5	0.5
Deg. Sol. Shell	1.069	2.712	3.742	-0.5
Present	1.262	2.565	3.681	-1	0.8
Deg. Sol. Shell	1.277	2.523	3.641	-1
Present	1.180	2.609	3.705	-0.5	0.8
Deg. Sol. Shell	1.186	2.615	3.674	-0.5
Present	1.383	2.521	3.620	-1	1
Deg. Sol. Shell	1.384	2.513	3.599	-1
Present	1.260	2.582	3.692	-0.5	1
Deg. Sol. Shell	1.255	2.601	3.722	-0.5

The good agreement between the two result sets becomes worst as the shallowness ratio $\frac{b}{Ry}$ increases. The reason is that the present approach is based on a general shell theory, whereas in in Liew and Lim (1996) a shallow shell formulation is considered. The first, fourth and eighth modes of vibration for the completely clamped moderately thick hyperbolic paraboloid with $\frac{Ry}{Rx}=-0.5$, and $\frac{b}{Ry}=0.8$are depicted in Figure 8.

Figure 8:
First (a), fourth (b), and eighth (c) modes of vibration for the case of completely clamped moderately thick hyperbolic paraboloid with $\frac{Ry}{Rx}=-0.5$, and $\frac{b}{Ry}=0.8$ (Example 5, Table 7).

6.6 Example 6: Simply supported thin cylindrical, spherical, and elliptic paraboloid shells

To conclude this work, we also study the natural frequencies of other kinds of surfaces, such as cylindrical, spherical, and elliptic paraboloid shells. In the following table, the boundary conditions are also those of a simply supported shell in order to compare them with the analysis of Liew and Lim [1996], Table 8, (ν=0.3, a⁄b=1, b⁄h=100).

Table 8:
Frequency parameter λ for simply supported thin cylindrical shells and spherical shells (Example 6) for various values of$\frac{b}{Ry}$and $\frac{Ry}{Rx}$ (ν=0.3, a⁄b=1, b⁄h=100). Comparison with Liew and Lim (1996) and the degenerated solid approach of Dvorkin and Bathe (1984).

Mode sequence number	1	4	8	$\frac{R_y}{R_x}$	$\frac{b}{R_y}$
Liew and Lim [1996]	36.841	82.302	131.11	0	0.1
Deg. Sol. Shell	36.024	81.687	130.85	0
Proposed	36.019	82.411	133.57	0
Liew and Lim [1996]	43.027	84.316	132.05	0.5	0.1
Deg. Sol. Shell	44.350	86.062	134.80	0.5
Proposed	42.208	85.154	134.21	0.5
Liew and Lim [1996]	53.049	87.829	133.51	1	0.1
Deg. Sol. Shell	52.907	87.374	133.89	1
Proposed	52.677	88.906	134.62	1
Liew and Lim [1996]	66.574	104.95	151.49	0	0.3
Deg. Sol. Shell	66.777	105.30	151.12	0
Proposed	66.152	106.02	153.73	0
Liew and Lim [1996]	86.927	118.41	158.69	0.5	0.3
Deg. Sol. Shell	87.437	117.06	161.16	0.5
Proposed	86.000	118.85	160.29	0.5
Liew and Lim [1996]	121.99	139.21	169.68	1	0.3
Deg. Sol. Shell	122.23	143.31	173.78	1
Proposed	119.65	141.52	170.83	1
Liew and Lim [1996]	88.431	140.07	182.95	0	0.5
Deg. Sol. Shell	87.624	141.22	180.97	0
Proposed	87.851	142.74	183.62	0
Liew and Lim [1996]	127.81	165.19	200.99	0.5	0.5
Deg. Sol. Shell	126.69	166.32	199.17	0.5
Proposed	128.63	167.28	202.15	0.5
Liew and Lim [1996]	188.59	201.27	239.87	1	0.5
Deg. Sol. Shell	193.56	199.93	239.52
Proposed	200.01	202.36	241.24

Finally, in Table 9 and Table 10, we present different results for completely clamped spherical and cylindrical shells. The results show no significant differences from those obtained while using the other finite elements. We can see that the results obtained with the classical formulation are closer to the shallow shell theory whereas the new formulation is closer to the deep shell theory. For simply supported shells, it can be seen that higher frequencies vary linearly with the shallowness ratio, however, the lower frequencies are not affected by increasing it. In this sense, the internal energy of the shell remains nearly unaltered with respect to the shell curvature for the fundamental frequencies.

Table 9:
Frequency parameter λ´ (=$\omega ab\sqrt{\frac{\rho h}{E}}$) for clamped (at the corners) thin spherical shells and various values of $\frac{a}{R}$. Comparison with ANSYS (ANSYS, 2011) and the degenerated solid approach of Dvorkin and Bathe (1984).

Mode sequence number	1	4	6	$\frac{a}{R_{}}$
Proposed	0.1093	0.3981	0.6594	0
Deg. Sol. Shell	0.1083	0.4011	0.6397	0
ANSYS	0.1082	0.3959	0.6585	0
Proposed	0.1096	2.669	3.733	0.005
Deg. Sol. Shell	0.1089	2.639	3.772	0.005
ANSYS	0.1086	2.657	3.721	0.005
Proposed	0.1115	2.629	3.7228	0.01
Deg. Sol. Shell	0.1197	2.069	3.7174	0.01
ANSYS	0.1091	2.620	3.7214	0.01

Table 10:
Frequency parameter λ´ (=$\omega ab\sqrt{\frac{\rho h}{E}}$) for clamped (at the corners) thin cylindrical shells and various values of $\frac{a}{R}$. Comparison with ANSYS (ANSYS, 2011) and the degenerated solid approach of Dvorkin and Bathe (1984).

Mode sequence number	1	4	6	$\frac{a}{R_{}}$
Proposed	0.0060	0.0169	0.0243	0
Deg. Sol. Shell	0.0061	0.0166	0.0241	0
ANSYS	0.0060	0.0168	0.0240	0
Proposed	0.0093	0.0173	0.0248	0.05
Deg. Sol. Shell	0.0092	0.0174	0.0246	0.05
ANSYS	0.0093	0.0171	0.0246	0.05
Proposed	0.0824	0.1297	0.1657	0.1
Deg. Sol. Shell	0.0824	0.1295	0.1623	0.1
ANSYS	0.0824	0.1294	0.1652	0.1

7 CONCLUSIONS

An efficient shell element for the linear dynamic analysis of moderately thick shells has been formulated. In spite of a greater initial complexity than other finite element models dedicated to this purpose, generally based on the degenerated solid approach, the new proposal has a series of advantages, such as its versatility, since it is based on a general and consistent shell model (tensor form) which permits an exact knowledge of the geometric parameters of the studied surface, and the reduced computational time calculation (compared to degenerate shell elements).

The tensor form of the stiffness and mass matrices (consistent mass matrix) and the generalized constitutive equations have been presented for this shell element, including the higher order moments deduced in the theory. To avoid the known locking problem of the stiffness matrix for thin shells, we have opted for the mixed interpolation of the tensorial components of the linear strain tensor (MITC technique). The results are in very good agreement with those reported in the literature and, besides, some interesting ideas with respect to higher order vibration modes can be highlighted. In this sense, the differences with other competitive finite elements become greater for thicker shells and for higher values of the ratios $\frac{b}{R_y}$.

Indeed, we have focused on several kinds of surfaces, such as cylindrical, spherical, elliptic paraboloid and doubly curved shells, with a wide range of shallowness ratios and different boundary conditions, constituting strict tests for validating the new finite element, as well as due to the lack of results for non-shallow and thick shells (deep shells). The extension of this general dynamic formulation to other attractive problems in structural mechanics like material and geometric nonlinear shell analysis or applications to composite shell models is direct and will be presented in future works.

8 ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support from the University of Córdoba in Spain (Plan Propio de Investigación de la Universidad de Córdoba y del Programa Operativo de fondos FEDER Andalucía), and CONICET (Argentine Council for Scientific and Technical Research). A. E. Albanesi also acknowledges the National Technological University of Argentina (UTN) for Grants PID ENUTNFE0004405 and PID MAUTIFE0005270TC, and the National Agency of Scientific and Technological Promotion of Argentina (ANPCYT) for the Grant PICT2015-3396.

References

• ANSYS, (2011). Academic Research Mechanical, Release 14, Help System.
• Aksu, T., (1997). A finite element formulation for free vibration analysis of shells of general shape. Computers and Structures 65, 687-694.
• Ahmad, S., Irons, B., and Zienkiewicz, O.C., (1970). Analysis of thick and thin shell structures by curved finite elements. International Journal for Numerical Methods in Engineering 2, 419-451.
• Alhazza, K, and Alhazza, A., (2004). A review of the vibrations of plates and shells. The Shock and Vibration Digest; 36(5):377-395.
• Bathe, K.J., and Dvorkin, E.N., (1986). A formulation of general shell elements - the use of mixed interpolation of tensorial components. International Journal for Numerical methods in Engineering. 22:697-722.
• Bhimaraddi, A., (1991). Free vibration analysis of doubly curved shallow shells on rectangular planform using three-dimensional elasticity theory. International Journal of Solids and Structures. 27, 897-913.
• Bucalem, M.L., y Bathe, K.J., (1993). Higher-order MITC general shells elements. International Journal for Numerical methods in Engineering. 36:3729-3754.
• Büchter, N., Ramm, E., (1992). Shell Theory versus Degeneration - A Comparison in Large Rotation Finite Element Analysis. International Journal for Numerical methods in Engineering 34, 39-59.
• Chakravorty, D., and Bandyopadhyay, J., (1995). On the free vibration of shallow shells. Journal of Sound and Vibration, 185(4):673-684.
• Chakravorty, D., Bandyopadhyay, J.N., and Sinha, P.K., (1995). Free vibration of point-supported laminated doubly curved shells-a finite element approach. Comput Struct, 54(2):191-8.
• Chapelle, D., and Bathe, K.J., (1998). Fundamental considerations for the finite element analysis of shell structures. Computers & Structures. 66 (1) January: 19-36.
• Chapelle, D., Suarez, IP., (2008). Detailed reliability assessment of triangular MITC elements for thin shells. Computers and Structures. 86(23- 24):2192-2202.
• Dvorkin E. N., and Bathe K.-J., (1984). A continuum mechanics based four-node shell element for general nonlinear analysis. Eng. Comput. 1(1):77-88.
• Fard, K. M., Livani M., Faramarz Ashenai Ghasemi., (2014). Improved high order free vibration analysis of thick double curved sandwich panels with transversely flexible cores. Latin American Journal of Solids and Structures. 11: 2284-2307.
• Germain, S. (1821). Recherches sur la the’orie des surfaces e’lastiques Mme. V. Courcier edt., Paris.
• Kang, J.H., Leissa, A.W., (2005). Free vibration analysis of complete paraboloidal shells of revolution with variable thickness and solid paraboloids from a three-dimensional theory Comput. Struct. 83: 2594-2608.
• Kienzler, R., Altenbach, H. and Ott, I., (2004). Theories of Plates and Shells: Critical Review and New Applications, Springer, Berlin.
• Kirchhoff, G.R., (1850). Uber das gleichgewicht und die bewegung einer elastischen Scheibe. Journal für die reine und angewandte Mathematik (Crelle's Journal) 40: 51-88.
• Kumbasar, N., Aksu, T., (1995). A finite element formulation for moderately thick shells of general shape. Comput Struct, 54:49-57(9).
• Lee, S. J., and Han, S. E., (2001). Free-vibration analysis of plates and shells with a nine-node assumed natural degenerated shell element, Journal of Sound and Vibration. 241(4), 605-633.
• Leissa, A. W., (1993). Vibration of Shells. Acoustical Society of America. Sewickly.
• Leissa, AW, Chang, J., (1996). Elastic deformation of thick, laminated composite shallow shells. Composite Structures. 35:153-70.
• Liew, K., and Lim, C., (1995). Vibratory behaviour of doubly curved shallow shells of curvilinear planform. Journal of Engineering Mechanics, 121(12):1277-1283.
• Liew, K., and Lim, C., (1996). Vibration of doubly-curved shallow shells. Acta Mechanica, 114:95-119.
• Liu, L., Chua, L. P., Ghista, D. N. (2006). Element-free Galerkin method for static and dynamic analysis of spatial shell structures, Journal of Sound and Vibration, 295:1-2 (388).
• Mantari, J.L., Oktem, A.S., Soares, C.G. (2011). Static and dynamic analysis of laminated composite and sandwich plates and shells by using a new higher-order shear deformation theory. Composite Structures 94: 37-49.
• Meiche, N.E., Tounsi, A., Zlane, N., Mechab, I., Bedia, E.A.A. (2011). A new hyperbolic shear deformation theory for buckling and vibration of functionally graded sandwich plate. International Journal of Mechanical Sciences 53: 237-47.
• Narita, Y. and Leissa, A.W., (1984). Vibrations of corner point supported shallow shells of rectangular planform. Earthquake Engineering Structures. 12(5):651-61.
• Nayak, P.K. (2012). Textbook of Tensor Calculus and Differential Geometry, PHI Learning. New Delhi
• Reissner, E. (1945). The effect of transverse shear deformation on the bending of elastic plates. Journal of Applied Mechanics 12: 69-76.
• Qatu, M.S. (1999). Accurate equations for laminated composite deep thick shells. International Journal of Solids and Structures, 36: 2917-41.
• Reddy JN, Liu CF., (1985). A higher-order shear deformation theory of laminated elastic shells. International Journal Engineering Science;23:440-7.
• Schwarze M. and Reese S., (2009). A reduced integration solid-shell finite element based on the EAS and the ANS concept-Geometrically linear problems. International Journal for Numerical Methods in Engineering, 80(10):1322-55.
• Stavridis, L., (1998). Dynamic analysis of shallow shells of rectangular base. Journal of Sound and Vibration, 218(5):861-882.
• Stavridis, L. and Armenakas, A., (1988). Analysis of shallow shells with rectangular projection. Applications Journal of Engineering Mechanics, 114(6):943-952.
• Tornabene, N., Fantuzzi, AA and Bacciocchi, M., (2014). Free vibrations of free-form doubly-curved shells made of functionally graded materials using higher-order equivalent single layer theories, Composites Part B: Engineering, 67(1), 490-509.
• Wempner, G.A., (1982). Mechanics of Solids, with Applications to thin Bodies, Springer. The Netherlands.
• Wempner, G.A. and Talaslidis, D., (2002). Mechanics of Solids and Shells: Theories and Approximations, CRC Press. Boca Raton.
• Yang, H., Saigal, S., Masud, A., and Kapania, R., (2004). A survey of recent shell finite elements. International Journal for Numerical Methods in Engineering. 47,101-127.
• Zienkiewicz, O. C, and Taylor, L.R., (2000). El Método de los Elementos Finitos. Mecánica de Sólidos. CIMNE (The International Centre for Numerical Methods in Engineering). Barcelona (España).

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  27 June 2019
• Date of issue
  2019

History

• Received
  30 Oct 2018
• Reviewed
  25 Mar 2019
• Accepted
  29 Mar 2019
• Published
  24 Apr 2019