# A Mixed Stress/Displacement Approach Model of Homogeneous Shells for Elastodynamic Problems

Authors: Axel Fernando Domínguez Alvarado; Alberto Díaz Díaz
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper presents the development of a model of homogeneous, moderately thick shells for elastodynamic problems. The model is obtained by adapting and modifying SAM-H model (stress approach model of homogeneous shells) developed by Domínguez Alvarado and Díaz in (2018) for static problems. In the dynamic version of SAM-H presented herein, displacements and stresses are approximated by polynomials of the out-of-plane coordinate. The stress approximation coincides with the static version of SAM-H when dynamic effects are neglected. The generalized forces and displacements appearing in the approximations are the same as those involved in a classical, moderately thick shell model (CS model) but the stress approximation adopted herein is more complex: the 3D motion equations and the stress boundary conditions at the faces of the shell are verified. The generalized motion and constitutive equations of dynamic SAM-H model are obtained by applying a variant of Euler–Lagrange equation which includes pertinently Hellinger–Reissner functional. In the constitutive equations, Poisson’s effect of out-of-plane normal stresses on in-plane strains is not ignored; this is one important feature of SAM-H. To test the accuracy of dynamic SAM-H model, the following structures were considered: a hollow sphere and a catenoid. In each case, eigenfrequencies are first calculated and then a frequency analysis is performed applying a harmonic load. The results are compared to those of a CS model, MITC6 (mixed interpolation of tensorial components with 6 nodes per element) shell element calculations, and solid finite element computations. In the two problems, CS, MITC6, and dynamic SAM-H models yield accurate eigenfrequencies and eigenmodes. Nevertheless, the frequency analysis performed in each case showed that dynamic SAM-H provides much more accurate amplitudes of stresses and displacements than the CS model and the MITC6 shell finite element technique.


---

## Body

## 1. Introduction

Shell-like structures are a common solution for light and innovative designs; their structural analysis may be performed by making use of finite element calculations. The application of 3D solid elements requires high computational resources. Surface elements (2D elements) are an alternative to obtain accurate results for a reasonable computing cost [1].Surface elements are based on shell models, which can be classified according to their ability to study thin [2, 3], moderately thick [4, 5], or thick shells [6–8]. The modelling in some shell models starts by proposing a displacement and/or stress field form in the thickness direction; in [9], Carrera classifies these models as axiomatic. In thin shell models, the displacement approach is similar to that in Kirchhoff–Love plate theory: 3 generalized displacements depending only on the in-plane coordinates are considered and the rotations are deduced from the generalized displacements and their spatial derivatives. In moderately thick shell models, the displacement approach is similar to that in Mindlin plate theory: 3 generalized displacements and 2 generalized rotations depending only on the in-plane coordinates are considered. Some examples of displacement approach models are Sanders model [10], the model built by Reddy in [5], and the Donnell–Mushtari–Vlasov (DMV) model [3]. Sanders shell theory derives six generalized equations of motion by means of the principle of virtual work. The model obtained in [5], hereinafter referred as to the classical shell model (CS model), is based on Hamilton’s principle and has 5 generalized motion equations. In DMV theory, the generalized motion equations are reduced from 5 to 3 by neglecting some terms of transverse generalized forces and assuming symmetric in-plane generalized forces. In the three theories above, the applied loads at the inner and outer faces of the shell do not appear explicitly in the motion equations. Other inconvenience of these models is the neglection of Poisson’s effect of out-of-plane normal stresses on in-plane strains; this can lead to an inaccurate displacement evaluation, particularly for thick and moderately thick shells [1]. Besides, the stress field deduced from the displacement approach models does not evaluate out-of-plane stresses accurately. In [11, 12], Tornabene and Viola used a displacement approach model and proposed a postprocessing technique (the differential quadrature method) for recovering and correcting the stress field in static problems. Another alternative to enhance the stress evaluation is the application of stress approach models; some examples of such models are those developed by Trefftz in [13] and Voyiadjis in [7, 8]. The Trefftz model does not verify the stress boundary conditions at the faces of the shell whilst the models proposed by Voyiadjis do verify these conditions and also the 3D equilibrium equations. To the best of our knowledge, in the literature, there are no pure stress approach models which are applied to dynamic problems and verify the 3D motion equations. Mixed approach models approximate both displacement and stress fields simultaneously [6, 14, 15]. Fang et al. developed in [6] the equations of their static mixed model by making use of the Hellinger–Reissner functional proposed in [16]. The stress field in their model fulfilled the stress boundary conditions at the faces of the shell but did not verify the 3D equilibrium equations.In order to evaluate the accuracy of dynamic shell models, researchers usually calculate eigenfrequencies and eigenmodes [17–21]. The eigenmodes and eigenfrequencies are compared to those obtained by other models and solid finite elements. Convergence speed is often discussed. In [22], Leissa analyzed and compared several theories of thin shells such as those developed by Love [2], Sanders [10], Flügge [23], Naghdi [4], and Donnell–Mushtari–Vlasov (DMV) [3]. Several loading conditions, geometries, and thickness-to-radius of curvature ratios were considered. In some cases, a 3D exact solution was available, and the accuracy of the eigenvalues obtained by the different models was analyzed. The DMV model was the least accurate. In dynamic problems where a harmonic load is imposed, it is important to carry out a frequency domain analysis and to determine the amplitudes of displacements and stresses. However, it is less usual to find these analysis with shell models in the literature. In [24], Poultangari and Nikkhah-Bahrami calculated the displacement amplitudes in thin cylinders given by the Flügge shell model [23] combined with an extended wave method and validated these results by comparing them with those obtained by solid finite elements.The resolution of the equations of shell models can be performed by making use of the finite element method. Many commercial finite element software packages include surface elements (shell elements). A well-known family of shell elements using a displacement approach model is the MITCn (mixed interpolation of tensorial components, where n refers to the number of nodes in the element) [25]. The development of these elements is based on pertinent displacement and strain interpolations [26]; the stress-strain constitutive equation of the shell model yields the element rigidity matrix. These elements do not exhibit shear and membrane locking numerical problems which commonly appear in low-order elements [27].In the present paper, a mixed stress/displacement approach model for elastodynamic problems of moderately thick, homogeneous shells is developed. The model is inspired by the stress approach model SAM-H for static shells proposed by Domínguez Alvarado and Díaz in [1]. This model proved to be very accurate in the evaluation of stresses since 3D equilibrium equations and boundary conditions at the faces of the shell are satisfied. The static SAM-H model equations were obtained by making use of the Hellinger–Reissner functional [16] of elastostatic problems. The generalized constitutive equations of the model take into account Poisson’s effect of out-of-plane normal stresses on in-plane strains, and this consideration provided higher accuracy to the model, particularly for thicker shells. The dynamic model developed in this paper coincides with the static SAM-H model when accelerations are neglected. For this reason, the name SAM-H is conserved and the new model is referred to as the dynamic SAM-H model. In the dynamic SAM-H model, 3D displacement and stress approximations are proposed. The latter approximation satisfies the stress boundary conditions at the faces of the shell and the 3D motion equation: this is the first originality of this study. The main goal is to provide an accurate estimation of the stress field. In order to obtain dynamic SAM-H equations, the Euler–Lagrange equation is applied. In this equation, the Hellinger–Reissner functional is used as a potential energy of the Lagrangian. The generalized motion and constitutive equations dynamic SAM-H involve explicitly the applied normal and shear stresses at the faces of the shell: this is a second originality of the model. As in the static version, another important feature of the dynamic model is the inclusion of Poisson’s effect of out-of-plane normal stresses on in-plane strains. The model equations are then applied to the dynamic analysis of two pressurized structures: a hollow sphere and a catenoid. For each structure, eigenvalues/eigenfrequencies and amplitudes of stresses/displacements in case of a forced vibration are compared to those obtained with solid finite elements, CS model, and MITC6 surface elements.This article is organized as follows. In Section2, the Hellinger–Reissner functional and Euler–Lagrange equation are recalled and combined to establish the method to derive the equations of the shell model. In Section 3, the Euler–Lagrange equation is applied to obtain dynamic SAM-H equations. In Section 4, two application examples are shown and discussed to validate the model equations. Section 5 contains the main conclusions of this study.


## 2. Hellinger–Reissner Functional Combined with Euler–Lagrange Equation for Shell Modelling

A homogeneous, doubly curved shell lying in the volumeΩ and delimited by the boundary ∂Ω is considered. Its thickness h is uniform, and ω denotes its middle surface. The principal curvatures are κ1 and κ2. Orthogonal curvilinear coordinates ξ1,ξ2,ξ3 are adopted; the unit vectors in the directions of the lines of curvature are e1 and e2 whilst the unit vector e3 normal to the middle surface is obtained by the cross product e1×e2. The surface metric tensor components are gαβ (α,β=1,2), and the following coefficients are defined: a1=g11 and a2=g22. The material is orthotropic and one of its orthotropy directions is direction 3.In what follows, the following assumptions and notations are adopted:(i)
3D vectors and tensors are written with bold face characters(ii)
First- and second-order tensors of the 2D space (defined by the basise1,e2) are noted in bold underlined and double underlined characters, respectively(iii)
r˙ and r¨ denote first and second time derivatives of field r(iv)
Greek subscripts take the values 1 and 2(v)
The outer (ξ3=h/2) and inner (ξ3=−h/2) faces of the shell are subjected to the stress vectors s+=τ¯++σ+e3 and s−=τ¯−+σ−e3, respectively; τ¯+ and τ¯− are the applied shear stress vectors; and σ+ and σ− are the applied normal stresses(vi)
A 3D body load vectorf=f1e1+f2e2+f3e3, uniform through the thickness direction, is applied; it is worth defining the in-plane body load vector f¯=f1e1+f2e2(vii)
The curvaturesκ1 and κ2 vary smoothly, in such a manner that ∂κβ/∂ξα≃0(viii)
Any term multiplied byhκαn is neglected if n≥2


### 2.1. Hellinger–Reissner Functional of Elastostatic Problems

On the boundaries∂Ωu and ∂Ωs of Ω, the body is subjected to the displacement vector ug and stress vector sg, respectively. The Hellinger–Reissner functional for elastic problems [16] applied to the shell is(1)H.R.u∗,σ∗=∫Ωσ∗:εu∗−f⋅u∗−weσ∗dΩ−∫∂Ωuσ∗⋅n⋅u∗−ugdS−∫∂Ωssg⋅u∗dS,where ⋅ and : stand for simple and double tensor contractions, respectively; n is the 3D outward-pointing normal vector to the boundary of Ω; u∗ is a piecewise C1 first-order tensor field; σ∗ is a piecewise C1 second-order symmetric tensor field; ε is Cauchy’s infinitesimal strain tensor; weσ∗=1/2σ∗:S:σ∗ is the elastic energy density; and S is the 3D, fourth-order compliance tensor. In what follows, a field with an uppercase ∗ denotes a field that is not necessarily the solution of the problem and varies in its corresponding virtual space. An integration by parts for the term containing σ∗:εu∗ yields an alternative form of the Hellinger–Reissner functional:(2)H.R.u∗,σ∗=−∫Ωdivσ∗⋅u∗+f⋅u∗+weσ∗dΩ+∫∂Ωuσ∗⋅n⋅ugdS+∫∂Ωsσ∗⋅n−sg⋅u∗dS,where div is the divergence operator in the 3D space. Reissner’s theorem [16] states that(i)
The solution of the mechanical problem is the coupleu,σ which makes the H.R. functional stationary (notice that no uppercase ∗ appears for the solution)(ii)
The stationarity of H.R. with respect tou∗ yields the equilibrium equations and the boundary conditions on the stress vector whilst its stationarity with respect to σ∗ provides the constitutive equations and the boundary conditions on the displacement field


### 2.2. Euler–Lagrange Equation

In Lagrangian mechanics, the evolution of a dynamic system is obtained by solving the Euler–Lagrange equation:(3)ddt∂ℒ∂q˙i−∂ℒ∂qi=0.ℒ is called the Lagrangian and is defined by(4)ℒqi,q˙i=Kq˙i−Πqi,where qi are the generalized coordinates (1≤i≤n, n is the number of generalized coordinates) and K and Π are the kinetic and potential energies, respectively. In this work, Π is considered as the Hellinger–Reissner functional defined in equation (1). In linear elastic problems, the generalized coordinates qi are the components of the stress and displacement fields: σ∗ and u∗. The expression of the kinetic energy is(5)Ku˙∗=12∫Ωρu˙∗⋅u˙∗dΩ,where ρ is the volumetric mass density of the material. Then, the Euler–Lagrange equation (3) yields(6)ddt∂K∂u˙∗+∂HR∂u∗=0,(7)∂HR∂σ∗=0.The coupleu,σ that satisfies equations (6) and (7) is the solution of the problem. From equation (6), the 3D motion equations and stress boundary conditions are obtained. From equation (7), one deduces the 3D constitutive equations and the displacement boundary conditions.


### 2.3. Case of a Shell Model Based on a Mixed Approach

When dealing with plates or shells, mixed-approach models provide an approximation of the 3D stress and displacement fields. In this approximation, polynomials of the coordinate across the thickness are usually applied. The coefficients of these polynomials involve the generalized forces (FG) and displacements (UG). These generalized fields depend only on the in-plane coordinates ξ1,ξ2. The stress approximation may be chosen in order to verify the 3D stress boundary conditions at the faces of the shell and/or to obtain a good estimation of the 3D motion equations.The Euler–Lagrange equation can be applied to derive the generalized equations of mixed-approach models. The generalized coordinates areFG and UG. The Euler–Lagrange equation yields(8)ddt∂K∂U˙G∗+∂HR∂UG∗=0,(9)∂HR∂FG∗=0.The generalized motion equations and force edge conditions are obtained from equation (8). The generalized constitutive equations and displacement edge conditions are given by equation (9).


## 2.1. Hellinger–Reissner Functional of Elastostatic Problems

On the boundaries∂Ωu and ∂Ωs of Ω, the body is subjected to the displacement vector ug and stress vector sg, respectively. The Hellinger–Reissner functional for elastic problems [16] applied to the shell is(1)H.R.u∗,σ∗=∫Ωσ∗:εu∗−f⋅u∗−weσ∗dΩ−∫∂Ωuσ∗⋅n⋅u∗−ugdS−∫∂Ωssg⋅u∗dS,where ⋅ and : stand for simple and double tensor contractions, respectively; n is the 3D outward-pointing normal vector to the boundary of Ω; u∗ is a piecewise C1 first-order tensor field; σ∗ is a piecewise C1 second-order symmetric tensor field; ε is Cauchy’s infinitesimal strain tensor; weσ∗=1/2σ∗:S:σ∗ is the elastic energy density; and S is the 3D, fourth-order compliance tensor. In what follows, a field with an uppercase ∗ denotes a field that is not necessarily the solution of the problem and varies in its corresponding virtual space. An integration by parts for the term containing σ∗:εu∗ yields an alternative form of the Hellinger–Reissner functional:(2)H.R.u∗,σ∗=−∫Ωdivσ∗⋅u∗+f⋅u∗+weσ∗dΩ+∫∂Ωuσ∗⋅n⋅ugdS+∫∂Ωsσ∗⋅n−sg⋅u∗dS,where div is the divergence operator in the 3D space. Reissner’s theorem [16] states that(i)
The solution of the mechanical problem is the coupleu,σ which makes the H.R. functional stationary (notice that no uppercase ∗ appears for the solution)(ii)
The stationarity of H.R. with respect tou∗ yields the equilibrium equations and the boundary conditions on the stress vector whilst its stationarity with respect to σ∗ provides the constitutive equations and the boundary conditions on the displacement field


## 2.2. Euler–Lagrange Equation

In Lagrangian mechanics, the evolution of a dynamic system is obtained by solving the Euler–Lagrange equation:(3)ddt∂ℒ∂q˙i−∂ℒ∂qi=0.ℒ is called the Lagrangian and is defined by(4)ℒqi,q˙i=Kq˙i−Πqi,where qi are the generalized coordinates (1≤i≤n, n is the number of generalized coordinates) and K and Π are the kinetic and potential energies, respectively. In this work, Π is considered as the Hellinger–Reissner functional defined in equation (1). In linear elastic problems, the generalized coordinates qi are the components of the stress and displacement fields: σ∗ and u∗. The expression of the kinetic energy is(5)Ku˙∗=12∫Ωρu˙∗⋅u˙∗dΩ,where ρ is the volumetric mass density of the material. Then, the Euler–Lagrange equation (3) yields(6)ddt∂K∂u˙∗+∂HR∂u∗=0,(7)∂HR∂σ∗=0.The coupleu,σ that satisfies equations (6) and (7) is the solution of the problem. From equation (6), the 3D motion equations and stress boundary conditions are obtained. From equation (7), one deduces the 3D constitutive equations and the displacement boundary conditions.


## 2.3. Case of a Shell Model Based on a Mixed Approach

When dealing with plates or shells, mixed-approach models provide an approximation of the 3D stress and displacement fields. In this approximation, polynomials of the coordinate across the thickness are usually applied. The coefficients of these polynomials involve the generalized forces (FG) and displacements (UG). These generalized fields depend only on the in-plane coordinates ξ1,ξ2. The stress approximation may be chosen in order to verify the 3D stress boundary conditions at the faces of the shell and/or to obtain a good estimation of the 3D motion equations.The Euler–Lagrange equation can be applied to derive the generalized equations of mixed-approach models. The generalized coordinates areFG and UG. The Euler–Lagrange equation yields(8)ddt∂K∂U˙G∗+∂HR∂UG∗=0,(9)∂HR∂FG∗=0.The generalized motion equations and force edge conditions are obtained from equation (8). The generalized constitutive equations and displacement edge conditions are given by equation (9).


## 3. Dynamic SAM-H Model

SAM-H (stress approach model of homogeneous shells) is a model originally developed for static problems of moderately thick shells [1]. Herein, this model is adapted and modified so as to enable the resolution of dynamic problems. Basically, the modification consists in adding an approximation of the 3D displacement field in order to use the methodology shown in Section 2.3. A new model including accelerations in the motion equations is obtained. When dynamic effects are neglected, the new model equations are the same as those of the static SAM-H equations developed in [1]. It is for this reason that we decided to preserve the name SAM-H for the new model even if it is based on a mixed approach: the new model will be referred to as the dynamic SAM-H model.


### 3.1. Approximation of Stress and Displacement Fields

For this approximation, the following polynomial basis of 4-th degreeξ3-polynomials is defined:(10)P0ξ3=1,P1ξ3=ξ3h,P2ξ3=−6ξ3h2+12,P3ξ3=−2ξ3h3+310ξ3hP4ξ3=−143ξ3h4+ξ3h2−140.The displacement approximation is chosen as follows:(11)u∗=U1∗ξ1,ξ2,tP0ξ3+hϕ1ξ1,ξ2,tP1ξ3U2∗ξ1,ξ2,tP0ξ3+hϕ2ξ1,ξ2,tP1ξ3U3∗ξ1,ξ2,tP0ξ3,where U1, U2, U3, ϕ1, and ϕ2 are the generalized displacements.The generalized forces of the model are membrane forcesN¯¯∗, shear forces Q¯∗, and moments M¯¯∗; their components are defined as follows:(12)Nαβ∗=∫−h/2h/2σαβ∗1+ξ3R3−βdξ3,Qα∗=∫−h/2h/2σα3∗1+ξ3R3−αdξ3,Mαβ∗=∫−h/2h/2ξ3σαβ∗1+ξ3R3−βdξ3,where R1 and R2 are the principal radii of curvature; the subscript 3−α is 2 if α=1, and 1 if α=2. The same can be said for the subscript 3−β. The stress tensor components are approximated by(13)σαβ∗ξ1,ξ2,ξ3,t=∑n=01σαβn∗ξ1,ξ2,tPnξ3,(14)σα3∗ξ1,ξ2,ξ3,t=∑n=03σα3n∗ξ1,ξ2,tPnξ3,(15)σ33∗ξ1,ξ2,ξ3,t=∑n=04σ33n∗ξ1,ξ2,tPnξ3,where σijn∗ (1≤i,j≤3) are stress coefficients which are expressed as linear combinations of the generalized forces. The polynomial degree of each stress component and the coefficients σijn∗ were chosen so as to verify the 3D stress boundary conditions at the faces of the shell and the 3D motion equation:(16)divσ∗+f=ρu¨∗.In Appendix A, the expressions of the coefficientsσijn∗ are shown.It is important to note that the displacement approximation is the same as that in the CS model. The generalized forces are the same in SAM-H and CS models. Nevertheless, the CS model obtains the 3D stress field from the 3D constitutive equations which require the calculation of the 3D strains by making use of the displacement approximation. It is for this reason that the 3D stress field in the CS model does not verify the 3D motion equations and the boundary conditions at the faces of the shell.


### 3.2. Generalized Motion Equations and Boundary Conditions on Generalized Forces

After introducing the 3D stress and displacement approximations in equation (2) and applying equation (8), one identifies the following generalized motion equations:(17)div¯N¯¯+κ¯¯⋅Q¯+hκ1+κ22τ¯+−τ¯−+τ¯++τ¯−+hf¯=ρhU¯¨+ρh3κ1+κ212ϕ¯¨,(18)divQ¯+σ+−σ−hκ1+κ22+σ++σ−−N11κ1−N22κ2+hf3=ρhU¨3,(19)div¯M¯¯−Q¯+h2κ1+κ24τ¯++τ¯−+h2τ¯+−τ¯−+h3κ1+κ212f¯=ρh3κ1+κ212U¯¨+ρh312ϕ¯¨,where κ¯¯=κ1e1⊗e1+κ2e2⊗e2, in which ⊗ stands for the tensor product operation. In the generalized motion equations above, div¯ and div stand for the divergences of second order tensor and vector fields in the 2D space of coordinates (ξ1, ξ2), respectively. If the accelerations were neglected in these equations, one would obtain the equilibrium equations of static SAM-H model developed in [1]. Let us mention that these generalized motion equations are similar to those in the CS model except for the applied stresses at the faces of the shell which do not appear in the equations of the CS model.Equation (8) also provides the edge conditions on the generalized forces:(20)N¯¯⋅n¯0=F¯g,Q¯⋅n¯0=F3g,M¯¯⋅n¯0=C¯g,where n¯0=n10e1+n20e2 is the outward-pointing normal vector to the edge of the middle surface ω and F¯g, F3g, and C¯g are the given membrane force vector, shear force, and bending moment vector at the shell edges, respectively. The edge forces and bending moments shown in equation (20) are deduced from the applied 3D stress vector sg at the boundary of the solid as follows:(21)Fig=∫−h/2h/2P0+hn102κ1+n202κ2P1sigdξ3,Cαg=∫−h/2h/2h212n102κ1+n202κ2P0ξ3sαgdξ3+∫−h/2h/2hP1ξ31−h6n102κ1+n202κ2sαgdξ3.


### 3.3. Generalized Constitutive Equations and Boundary Conditions on Generalized Displacements

By introducing the stress and displacement approximation in the Hellinger–Reissner functional in equation (1), one identifies the following generalized strain tensors:(22)ε¯¯∗=grad¯¯U¯∗+U3∗κ¯¯,χ¯¯∗=grad¯¯ϕ¯∗,d¯∗=grad¯¯U3∗+ϕ¯∗−κ¯¯⋅U¯∗,where U¯∗=U1∗e1+U2∗e2, ϕ¯∗=ϕ1∗e1+ϕ2∗e2, and κ¯¯′=κ2e1⊗e1+κ1e2⊗e2 in which ⊗ stands for the tensor product operation. The operators grad¯¯ and grad¯ are 2D gradients of vector and scalar fields, respectively.After introducing the 3D stress and displacement approximations in equation (2) and applying equation (9), one obtains the generalized constitutive equations. So as to ease the derivation of these equations, the expression of the elastic energy density weσ∗=1/2σ∗:S:σ∗ is simplified as explained in Appendix B. In order to condense the obtained constitutive equations, the following vectors of generalized forces and strains are defined:(23)q=N^M^,ε=ε^χ^,d¯=d1d2,where(24)N^=N11N22N12N21,M^=M11M22M12M21,ε^=ε11ε22ε12ε21,χ^=χ11χ22χ12χ21.The constitutive equations are the same as those obtained in [1] for the static version of the model; they are written as follows:(25)q=Kε+σ+k++σ−k−+f3k3,Q¯=Ld¯+L+τ¯++L−τ¯−,where K is an 8×8 stiffness matrix, k+, k−, and k3 are 8-component stiffness vectors, and L, L+, and L− are 2×2 stiffness matrices. The expressions of these matrices are shown in Appendix C. As already mentioned in [1], for the static SAM-H model, the consideration of k+, k−, k3, L+, and L− represents an important originality which provides significant accuracy to SAM-H model, particularly when dealing with moderately thick shells. The terms k+ and k− are related to the consideration of Poisson’s effect of out-of-plane normal stresses on in-plane strains. Let us note that the form of these generalized constitutive equations is similar to that in the CS model. The expressions of the terms in the stiffness matrices are not the same in SAM-H and CS models. Another important difference is that the applied stresses at the shell faces and the volumetric force f3 are not taken into account in the constitutive equations of the CS model. In fact, Poisson’s effect of out-of-plane normal stresses on in-plane strains is neglected in the CS model as in most plate and shell models.Equation (9) also provides the following edge conditions on the generalized displacement fields:(26)U¯=U¯g,U3=U3g,ϕ¯=ϕ¯g,where U¯g, U3g, andϕ¯gare the given in-plane displacement vector, out-of-plane displacement, and rotation vector at the shell edges, respectively.


## 3.1. Approximation of Stress and Displacement Fields

For this approximation, the following polynomial basis of 4-th degreeξ3-polynomials is defined:(10)P0ξ3=1,P1ξ3=ξ3h,P2ξ3=−6ξ3h2+12,P3ξ3=−2ξ3h3+310ξ3hP4ξ3=−143ξ3h4+ξ3h2−140.The displacement approximation is chosen as follows:(11)u∗=U1∗ξ1,ξ2,tP0ξ3+hϕ1ξ1,ξ2,tP1ξ3U2∗ξ1,ξ2,tP0ξ3+hϕ2ξ1,ξ2,tP1ξ3U3∗ξ1,ξ2,tP0ξ3,where U1, U2, U3, ϕ1, and ϕ2 are the generalized displacements.The generalized forces of the model are membrane forcesN¯¯∗, shear forces Q¯∗, and moments M¯¯∗; their components are defined as follows:(12)Nαβ∗=∫−h/2h/2σαβ∗1+ξ3R3−βdξ3,Qα∗=∫−h/2h/2σα3∗1+ξ3R3−αdξ3,Mαβ∗=∫−h/2h/2ξ3σαβ∗1+ξ3R3−βdξ3,where R1 and R2 are the principal radii of curvature; the subscript 3−α is 2 if α=1, and 1 if α=2. The same can be said for the subscript 3−β. The stress tensor components are approximated by(13)σαβ∗ξ1,ξ2,ξ3,t=∑n=01σαβn∗ξ1,ξ2,tPnξ3,(14)σα3∗ξ1,ξ2,ξ3,t=∑n=03σα3n∗ξ1,ξ2,tPnξ3,(15)σ33∗ξ1,ξ2,ξ3,t=∑n=04σ33n∗ξ1,ξ2,tPnξ3,where σijn∗ (1≤i,j≤3) are stress coefficients which are expressed as linear combinations of the generalized forces. The polynomial degree of each stress component and the coefficients σijn∗ were chosen so as to verify the 3D stress boundary conditions at the faces of the shell and the 3D motion equation:(16)divσ∗+f=ρu¨∗.In Appendix A, the expressions of the coefficientsσijn∗ are shown.It is important to note that the displacement approximation is the same as that in the CS model. The generalized forces are the same in SAM-H and CS models. Nevertheless, the CS model obtains the 3D stress field from the 3D constitutive equations which require the calculation of the 3D strains by making use of the displacement approximation. It is for this reason that the 3D stress field in the CS model does not verify the 3D motion equations and the boundary conditions at the faces of the shell.


## 3.2. Generalized Motion Equations and Boundary Conditions on Generalized Forces

After introducing the 3D stress and displacement approximations in equation (2) and applying equation (8), one identifies the following generalized motion equations:(17)div¯N¯¯+κ¯¯⋅Q¯+hκ1+κ22τ¯+−τ¯−+τ¯++τ¯−+hf¯=ρhU¯¨+ρh3κ1+κ212ϕ¯¨,(18)divQ¯+σ+−σ−hκ1+κ22+σ++σ−−N11κ1−N22κ2+hf3=ρhU¨3,(19)div¯M¯¯−Q¯+h2κ1+κ24τ¯++τ¯−+h2τ¯+−τ¯−+h3κ1+κ212f¯=ρh3κ1+κ212U¯¨+ρh312ϕ¯¨,where κ¯¯=κ1e1⊗e1+κ2e2⊗e2, in which ⊗ stands for the tensor product operation. In the generalized motion equations above, div¯ and div stand for the divergences of second order tensor and vector fields in the 2D space of coordinates (ξ1, ξ2), respectively. If the accelerations were neglected in these equations, one would obtain the equilibrium equations of static SAM-H model developed in [1]. Let us mention that these generalized motion equations are similar to those in the CS model except for the applied stresses at the faces of the shell which do not appear in the equations of the CS model.Equation (8) also provides the edge conditions on the generalized forces:(20)N¯¯⋅n¯0=F¯g,Q¯⋅n¯0=F3g,M¯¯⋅n¯0=C¯g,where n¯0=n10e1+n20e2 is the outward-pointing normal vector to the edge of the middle surface ω and F¯g, F3g, and C¯g are the given membrane force vector, shear force, and bending moment vector at the shell edges, respectively. The edge forces and bending moments shown in equation (20) are deduced from the applied 3D stress vector sg at the boundary of the solid as follows:(21)Fig=∫−h/2h/2P0+hn102κ1+n202κ2P1sigdξ3,Cαg=∫−h/2h/2h212n102κ1+n202κ2P0ξ3sαgdξ3+∫−h/2h/2hP1ξ31−h6n102κ1+n202κ2sαgdξ3.


## 3.3. Generalized Constitutive Equations and Boundary Conditions on Generalized Displacements

By introducing the stress and displacement approximation in the Hellinger–Reissner functional in equation (1), one identifies the following generalized strain tensors:(22)ε¯¯∗=grad¯¯U¯∗+U3∗κ¯¯,χ¯¯∗=grad¯¯ϕ¯∗,d¯∗=grad¯¯U3∗+ϕ¯∗−κ¯¯⋅U¯∗,where U¯∗=U1∗e1+U2∗e2, ϕ¯∗=ϕ1∗e1+ϕ2∗e2, and κ¯¯′=κ2e1⊗e1+κ1e2⊗e2 in which ⊗ stands for the tensor product operation. The operators grad¯¯ and grad¯ are 2D gradients of vector and scalar fields, respectively.After introducing the 3D stress and displacement approximations in equation (2) and applying equation (9), one obtains the generalized constitutive equations. So as to ease the derivation of these equations, the expression of the elastic energy density weσ∗=1/2σ∗:S:σ∗ is simplified as explained in Appendix B. In order to condense the obtained constitutive equations, the following vectors of generalized forces and strains are defined:(23)q=N^M^,ε=ε^χ^,d¯=d1d2,where(24)N^=N11N22N12N21,M^=M11M22M12M21,ε^=ε11ε22ε12ε21,χ^=χ11χ22χ12χ21.The constitutive equations are the same as those obtained in [1] for the static version of the model; they are written as follows:(25)q=Kε+σ+k++σ−k−+f3k3,Q¯=Ld¯+L+τ¯++L−τ¯−,where K is an 8×8 stiffness matrix, k+, k−, and k3 are 8-component stiffness vectors, and L, L+, and L− are 2×2 stiffness matrices. The expressions of these matrices are shown in Appendix C. As already mentioned in [1], for the static SAM-H model, the consideration of k+, k−, k3, L+, and L− represents an important originality which provides significant accuracy to SAM-H model, particularly when dealing with moderately thick shells. The terms k+ and k− are related to the consideration of Poisson’s effect of out-of-plane normal stresses on in-plane strains. Let us note that the form of these generalized constitutive equations is similar to that in the CS model. The expressions of the terms in the stiffness matrices are not the same in SAM-H and CS models. Another important difference is that the applied stresses at the shell faces and the volumetric force f3 are not taken into account in the constitutive equations of the CS model. In fact, Poisson’s effect of out-of-plane normal stresses on in-plane strains is neglected in the CS model as in most plate and shell models.Equation (9) also provides the following edge conditions on the generalized displacement fields:(26)U¯=U¯g,U3=U3g,ϕ¯=ϕ¯g,where U¯g, U3g, andϕ¯gare the given in-plane displacement vector, out-of-plane displacement, and rotation vector at the shell edges, respectively.


## 4. Results and Discussion

In this section, SAM-H is applied to the resolution of two structural dynamic problems. In the first problem, an analytical resolution of the model equations is performed. In the other problem, SAM-H equations are solved numerically by making use of a commercial finite element software: COMSOL Multiphysics (5.3a version). In this software, the PDE physics interface is applied and the equations are implemented manually in a coefficient form. So as to validate the SAM-H model, its results are compared to those obtained with solid finite element calculations, the CS model developed in [5], and MITC6 shell finite elements [25] included in COMSOL Multiphysics.


### 4.1. Analytical Results for a Pressurized Hollow Sphere

For an analytical resolution of the SAM-H dynamic equations, a hollow sphere subjected to a harmonic internal pressure is considered; the external face of the sphere is stress free. The radius of the middle surface and the thickness-to-radius ratio are denoted byR and η, respectively. The material is isotropic and characterized by its Young’s modulus E, Poisson’s ratio ν, and density ρ. The internal pressure p is given by pt=p0sin2πf¯t, where p0 and f¯ are the amplitude of the internal pressure and the excitation frequency, respectively. SAM-H results are solved analytically and compared to those obtained with the classical shell model (CS) and solid finite elements. Owing to a spherical symmetry (only radial displacements are allowed), the 3D solid problem can be reduced to a 1D problem (the spatial coordinate is the radial position) and the shell models are reduced to a 0D problem. The solid problem is solved by making use of COMSOL Multiphysics software.Let us start by analyzing the eigenfrequencies of the vibration modes that are excited by such a load. SAM-H and CS models yield only one vibration mode; the eigenfrequencies are identical for both models:(27)fSAM‐H=fCS=12π2EρR21−ν.The number of eigenfrequencies in the solid problem is infinite. The first eigenfrequencyf1 is related to a vibration mode analogue to those obtained by the shell models. The other eigenfrequencies fn (n≥2) are related to vibrations which are not reproduced by the shell models due to the simplicity of the adopted expression for the kinematic field u3 (see equation (11)). For the particular problem studied herein and considering a moderately thick shell (η≤0.3), the determination of the fundamental eigenfrequency and its corresponding vibration mode is enough to study the dynamic response of the structure. The relative error in the calculation of the fundamental eigenfrequency with SAM-H and CS models is defined by(28)δf1=fSAM‐H−f1f1.In Figure1(a), this relative error is plotted against the thickness to radius ratio η for 5 values of Poisson’s ratio (−0.4, −0.2, 0, 0.2, and 0.4). It can be seen that for moderately thick shells, the predicted eigenfrequency by both shell models (SAM-H and CS models) is very accurate: the relative error δf1 is less than 1.5%. It is worth mentioning that neither Young’s modulus nor the density affects the curves obtained in Figure 1(a).Figure 1
Relative errorδf1 vs. thickness to radius ratio (a) and correction factor vs. thickness to radius ratio (b).
(a)
(b)Let us now make a frequency analysis applying a harmonic pressure. Transient vibrations are not calculated, and only the forced, steady-state vibrations are analyzed. The radial displacement amplitude given by the CS model is(29)U3CS=p0Rη1−ν2E−ρω¯2R21−ν,where ω¯=2πf¯ is the excitation angular frequency. SAM-H model yields the following displacement amplitude:(30)U3SAM=kSAM×U3CS,where kSAM can be considered as a correction factor appearing in SAM-H model to improve the determination of the displacement amplitude; the expression of kSAM is given by(31)kSAM=1−η−q−,q−=2η45η−1050ν+865ην−1414ην2+1050ν221001−ν2.It is worth mentioning that the second term (−η) appearing in the expression of kSAM comes from a contribution of the normal stress σ− in the motion equation (18) whilst the term q− comes from the constitutive equations owing to Poisson’s effect of out-of-plane normal stresses on in-plane strains. In Figure 1(b), the correction factor kSAM is plotted against the thickness to radius ratio η for 5 values of Poisson’s ratio (−0.4, −0.2, 0, 0.2, and 0.4). One can observe that for very thin shells (η≤0.01), no amendment is performed by SAM-H (kSAM≃1), but for thicker shells, the correction factor has to be taken into account, particularly for a negative Poisson’s ratio.Let us consider the following parameters and properties:R=1m, η=0.3, ρ=1000kg/m3, E=1GPa, ν=−0.4, and p0=1kPa. As shown in Figure 1(b), it corresponds to the case where the largest differences between SAM-H and CS models are expected. In Figure 2(a), the amplitudes of the displacements at the middle surface provided by these models and solid finite elements (SFE) are plotted against the excitation frequency. The selected range of the excitation frequency includes only the first eigenfrequency of the sphere. At 130 Hz, as compared to SFE, the relative errors in the displacement amplitude are −5.1% and 49.5% for SAM-H and CS models, respectively. It can be seen that SAM-H makes a better approximation of the displacement amplitude than the CS model. The displacement amplitude peak of the shell models occurs at a 190 Hz frequency. In structural design, for safety reasons, it is a common practice to avoid resonance effects by ensuring that the fundamental eigenfrequency is greater enough than the excitation frequency. Let us assume that the excitation frequencies are at least 20% less than the fundamental eigenfrequency calculated by SAM-H. In Figure 2(b), the amplitude of the in-plane stresses σ11 calculated at the middle surface by SAM-H, CS, and SFE techniques is plotted against the excitation frequency f¯. In this figure, the curves of the out-of-plane normal stresses σ33 evaluated at the middle surface by SAM-H and SFE are also included. Let us recall that the CS model is unable to evaluate this stress. It can be seen that within the range of excitation frequencies, SAM-H results of σ11 are more accurate than those of the CS model and the evaluation of σ33 by SAM-H at the middle surface is also accurate. Let us plot in Figure 3 the normal stress σ33 as a function of the position ξ3 across the thickness for two excitation frequencies: 50 and 150 Hz.Figure 2
Displacement amplitude vs. excitation frequency (a);σ11 and σ33 stress amplitudes vs. excitation frequency (b). CS, SAM-H, and SFE results.
(a)
(b)Figure 3
σ
33 amplitude vs. ξ3 (SAM-H and SFE results).The case of a pressurized sphere has shown that SAM-H and CS provide the same prediction of the fundamental eigenfrequency, and the value is very accurate. Nevertheless, SAM-H provides a much better approximation of displacement and stress field amplitudes in a frequency analysis.


### 4.2. Numerical Results for a Catenoid

In this section, a 0.9 m thick catenoid is considered; its middle surface is shown in Figure4(a) and is defined by the following Cartesian coordinates:(32)x=acoshscosϕ,xy=−acoshssinϕ,z=as,where a=3m, s∈−1,0 and ϕ∈0,2π. The curvilinear coordinates ξ1 and ξ2 are selected as follows:(33)ξ1=asinhs+sinh1,ξ2=aϕcoshs.Figure 4
Catenoid middle surface (a) and example of a 2D meshing (b).
(a)
(b)The principal curvatures have opposite signs:κ1=−κ2. The average thickness-to-radius ratio η¯ is(34)η¯=∫ωηdω∫ωdω≃0.21,and the maximum thickness-to-radius ratio is ηmax=0.3. Consequently, the catenoid is virtually a thick shell. The material is isotropic and its Young’s modulus, Poisson’s ratio, and density are E=25GPa, ν=0.2, and ρ=2300kg/m3, respectively.For all the shell models (SAM-H, CS, and MITC6), triangular elements with 6 nodes using quadratic polynomials were applied; a quasiuniform element size was selected for the meshing. In Figure4(b), an example of the meshing of the middle surface is shown. For the solid finite element calculations, tetrahedral quadratic elements were used. All the numerical results of solid finite element computations shown in this section converge.


#### 4.2.1. Eigenfrequencies and Convergence Analysis

The basis of the catenoid is fixed whilst its upper boundary is stress free. This means that for SAM-H and CS models:U1=U2=U3=0 and ϕ1=ϕ2=0 at ξ1=0; N11=N21=Q1=0 and M11=M21=0 at ξ1=a. No displacement boundary conditions are imposed at the inner and outer faces of the shell. For these shell models, in order to analyze the convergence of results, different meshes were tested; each mesh has a quasiuniform element size and in mesh number i, the element size is smaller than in mesh number i−1. The first ten eigenfrequencies were calculated, and the approximate error in the calculation of eigenfrequency number j is analyzed:(35)δfij=fij−fi−1jfij100%,where fij is the estimation of eigenfrequency j using mesh i. The slowest convergence was observed for the last eigenfrequency (j = 10) for all shell computations. In Figure 5, the approximate error δf10 is plotted against the number of degrees of freedom (DoF) for SAM-H, CS, and MITC6 models. The convergence of both SAM-H and CS calculations is faster than that of MITC6 computations. In what follows, the mesh is selected so as to have values of eigenfrequencies with an approximate error less than 0.01%. In Table 1, the first ten eigenfrequencies calculated by SFE are listed. The eigenfrequencies calculated by SAM-H, CS, and MITC6 models and associated to the eigenmodes related to the SFE eigenfrequencies are also listed in this table. As compared to SFE results, a relative error less than 4% is obtained for all shell models considered herein. Since the structure is virtually a thick shell, all models yield good quality results for these eigenfrequencies. In Figure 6, the normalized displacement amplitude of the first mode is shown considering SFE, SAM-H, CS, and MITC6 calculations. Once again, the shell models provide similar eigenmodes to that obtained with SFE.Figure 5
Convergence of SAM-H, CS, and MITC6 models.Table 1
First 10 eigenfrequencies (in Hz) computed by SFE, SAM-H, CS, and MITC6.
j
SFE
SAM-H
CS
MITC6
1, 2
119.60
118.48
119.77
118.71
3, 4
122.98
121.77
122.82
122.82
5
141.84
140.09
142.38
142.67
6, 7
154.23
152.43
154.47
151.26
8
198.17
196.62
196.06
196.68
9, 10
198.32
196.36
196.60
191.72Figure 6
First eigenmode predicted by (a) SAM-H, (b) CS, (c) MITC6, and (d) SFE models.


#### 4.2.2. Frequency Domain Analysis for the Case of Internal Pressure

The basis of the catenoid remains fixed and the upper edge is stress free. A harmonic pressurep is applied at the inner face of the shell (internal pressure); its amplitude p0 is 1 MPa, and its frequency f¯ is varied so as to perform a frequency domain analysis. The outer face of the shell is stress free. To reduce the computational cost of SFE calculations, the 3D problem is reduced to a 2D axisymmetric problem. In Figure 7, the amplitudes of U3 at the top of the catenoid calculated by shell models and SFE (at the middle surface) are shown. It can be seen that SAM-H is slightly more accurate than CS and MITC6 for this case. For all the models, a displacement peak is expected at a frequency close to 83 Hz.Figure 7
Displacement amplitude vs. excitation frequency curve predicted by SFE, SAM-H, CS, and MITC6 models.For safety reasons, the excitation frequencyf¯ should be less than 83 Hz. Let us consider the case when f¯ = 70 Hz and analyze the amplitude of stresses along the blue curve parallel to direction 1 shown in Figure 4(a). In Figures 8(a) and 8(b), the amplitudes of the in-plane normal stresses σ11 and σ22 at the internal face of the shell (ξ3=−h/2) are plotted against the ξ1 coordinate, respectively. Globally, the shape and maximum values of the SFE curves in Figures 8(a) and 8(b) are approximated better by SAM-H calculations. It is worth mentioning that σ22 exhibits the highest value among the other stress components in the curvilinear coordinate system. As compared to SFE, the relative errors in the maximum value of σ22 predicted by SAM-H, MITC6 and CS models are 2.8%, −9.7%, and −10.5%, respectively. SAM-H in-plane normal stresses are then more accurate than those of CS and MITC6 techniques.Figure 8
σ
11 vs. ξ1 (a) and σ22 vs. ξ1 (b) curves at the external face of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)In Figure9(a) the amplitude of the shear stress σ13 at the middle surface is plotted against the curvilinear coordinate ξ1. The evolution of the SFE σ13 stress along the line of curvature is approximated more accurately by SAM-H than MITC6 and CS models. The maximum value of the SFE shear stress is better approximated by SAM-H and MITC6 models than by the CS model. Let us now analyze the accuracy of SAM-H in the calculation of the out-of-plane normal stress σ33. It is worth recalling that CS and MITC6 techniques do not calculate this stress. At the inner and outer faces of the shell, since SAM-H verifies the 3D boundary conditions, SAM-H yields the amplitudes σ33=1MPa and σ33=0MPa, respectively. At the middle surface, the amplitudes of the stresses σ33 calculated by SFE and SAM-H are plotted against ξ1 in Figure 9(b). As in the case of the pressurized hollow sphere in Section 4.1, SAM-H yields good quality predictions of σ33.Figure 9
σ
13 vs. ξ1 (a) and σ33 vs. ξ1 (b) curves at the middle surface of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)


## 4.1. Analytical Results for a Pressurized Hollow Sphere

For an analytical resolution of the SAM-H dynamic equations, a hollow sphere subjected to a harmonic internal pressure is considered; the external face of the sphere is stress free. The radius of the middle surface and the thickness-to-radius ratio are denoted byR and η, respectively. The material is isotropic and characterized by its Young’s modulus E, Poisson’s ratio ν, and density ρ. The internal pressure p is given by pt=p0sin2πf¯t, where p0 and f¯ are the amplitude of the internal pressure and the excitation frequency, respectively. SAM-H results are solved analytically and compared to those obtained with the classical shell model (CS) and solid finite elements. Owing to a spherical symmetry (only radial displacements are allowed), the 3D solid problem can be reduced to a 1D problem (the spatial coordinate is the radial position) and the shell models are reduced to a 0D problem. The solid problem is solved by making use of COMSOL Multiphysics software.Let us start by analyzing the eigenfrequencies of the vibration modes that are excited by such a load. SAM-H and CS models yield only one vibration mode; the eigenfrequencies are identical for both models:(27)fSAM‐H=fCS=12π2EρR21−ν.The number of eigenfrequencies in the solid problem is infinite. The first eigenfrequencyf1 is related to a vibration mode analogue to those obtained by the shell models. The other eigenfrequencies fn (n≥2) are related to vibrations which are not reproduced by the shell models due to the simplicity of the adopted expression for the kinematic field u3 (see equation (11)). For the particular problem studied herein and considering a moderately thick shell (η≤0.3), the determination of the fundamental eigenfrequency and its corresponding vibration mode is enough to study the dynamic response of the structure. The relative error in the calculation of the fundamental eigenfrequency with SAM-H and CS models is defined by(28)δf1=fSAM‐H−f1f1.In Figure1(a), this relative error is plotted against the thickness to radius ratio η for 5 values of Poisson’s ratio (−0.4, −0.2, 0, 0.2, and 0.4). It can be seen that for moderately thick shells, the predicted eigenfrequency by both shell models (SAM-H and CS models) is very accurate: the relative error δf1 is less than 1.5%. It is worth mentioning that neither Young’s modulus nor the density affects the curves obtained in Figure 1(a).Figure 1
Relative errorδf1 vs. thickness to radius ratio (a) and correction factor vs. thickness to radius ratio (b).
(a)
(b)Let us now make a frequency analysis applying a harmonic pressure. Transient vibrations are not calculated, and only the forced, steady-state vibrations are analyzed. The radial displacement amplitude given by the CS model is(29)U3CS=p0Rη1−ν2E−ρω¯2R21−ν,where ω¯=2πf¯ is the excitation angular frequency. SAM-H model yields the following displacement amplitude:(30)U3SAM=kSAM×U3CS,where kSAM can be considered as a correction factor appearing in SAM-H model to improve the determination of the displacement amplitude; the expression of kSAM is given by(31)kSAM=1−η−q−,q−=2η45η−1050ν+865ην−1414ην2+1050ν221001−ν2.It is worth mentioning that the second term (−η) appearing in the expression of kSAM comes from a contribution of the normal stress σ− in the motion equation (18) whilst the term q− comes from the constitutive equations owing to Poisson’s effect of out-of-plane normal stresses on in-plane strains. In Figure 1(b), the correction factor kSAM is plotted against the thickness to radius ratio η for 5 values of Poisson’s ratio (−0.4, −0.2, 0, 0.2, and 0.4). One can observe that for very thin shells (η≤0.01), no amendment is performed by SAM-H (kSAM≃1), but for thicker shells, the correction factor has to be taken into account, particularly for a negative Poisson’s ratio.Let us consider the following parameters and properties:R=1m, η=0.3, ρ=1000kg/m3, E=1GPa, ν=−0.4, and p0=1kPa. As shown in Figure 1(b), it corresponds to the case where the largest differences between SAM-H and CS models are expected. In Figure 2(a), the amplitudes of the displacements at the middle surface provided by these models and solid finite elements (SFE) are plotted against the excitation frequency. The selected range of the excitation frequency includes only the first eigenfrequency of the sphere. At 130 Hz, as compared to SFE, the relative errors in the displacement amplitude are −5.1% and 49.5% for SAM-H and CS models, respectively. It can be seen that SAM-H makes a better approximation of the displacement amplitude than the CS model. The displacement amplitude peak of the shell models occurs at a 190 Hz frequency. In structural design, for safety reasons, it is a common practice to avoid resonance effects by ensuring that the fundamental eigenfrequency is greater enough than the excitation frequency. Let us assume that the excitation frequencies are at least 20% less than the fundamental eigenfrequency calculated by SAM-H. In Figure 2(b), the amplitude of the in-plane stresses σ11 calculated at the middle surface by SAM-H, CS, and SFE techniques is plotted against the excitation frequency f¯. In this figure, the curves of the out-of-plane normal stresses σ33 evaluated at the middle surface by SAM-H and SFE are also included. Let us recall that the CS model is unable to evaluate this stress. It can be seen that within the range of excitation frequencies, SAM-H results of σ11 are more accurate than those of the CS model and the evaluation of σ33 by SAM-H at the middle surface is also accurate. Let us plot in Figure 3 the normal stress σ33 as a function of the position ξ3 across the thickness for two excitation frequencies: 50 and 150 Hz.Figure 2
Displacement amplitude vs. excitation frequency (a);σ11 and σ33 stress amplitudes vs. excitation frequency (b). CS, SAM-H, and SFE results.
(a)
(b)Figure 3
σ
33 amplitude vs. ξ3 (SAM-H and SFE results).The case of a pressurized sphere has shown that SAM-H and CS provide the same prediction of the fundamental eigenfrequency, and the value is very accurate. Nevertheless, SAM-H provides a much better approximation of displacement and stress field amplitudes in a frequency analysis.


## 4.2. Numerical Results for a Catenoid

In this section, a 0.9 m thick catenoid is considered; its middle surface is shown in Figure4(a) and is defined by the following Cartesian coordinates:(32)x=acoshscosϕ,xy=−acoshssinϕ,z=as,where a=3m, s∈−1,0 and ϕ∈0,2π. The curvilinear coordinates ξ1 and ξ2 are selected as follows:(33)ξ1=asinhs+sinh1,ξ2=aϕcoshs.Figure 4
Catenoid middle surface (a) and example of a 2D meshing (b).
(a)
(b)The principal curvatures have opposite signs:κ1=−κ2. The average thickness-to-radius ratio η¯ is(34)η¯=∫ωηdω∫ωdω≃0.21,and the maximum thickness-to-radius ratio is ηmax=0.3. Consequently, the catenoid is virtually a thick shell. The material is isotropic and its Young’s modulus, Poisson’s ratio, and density are E=25GPa, ν=0.2, and ρ=2300kg/m3, respectively.For all the shell models (SAM-H, CS, and MITC6), triangular elements with 6 nodes using quadratic polynomials were applied; a quasiuniform element size was selected for the meshing. In Figure4(b), an example of the meshing of the middle surface is shown. For the solid finite element calculations, tetrahedral quadratic elements were used. All the numerical results of solid finite element computations shown in this section converge.


### 4.2.1. Eigenfrequencies and Convergence Analysis

The basis of the catenoid is fixed whilst its upper boundary is stress free. This means that for SAM-H and CS models:U1=U2=U3=0 and ϕ1=ϕ2=0 at ξ1=0; N11=N21=Q1=0 and M11=M21=0 at ξ1=a. No displacement boundary conditions are imposed at the inner and outer faces of the shell. For these shell models, in order to analyze the convergence of results, different meshes were tested; each mesh has a quasiuniform element size and in mesh number i, the element size is smaller than in mesh number i−1. The first ten eigenfrequencies were calculated, and the approximate error in the calculation of eigenfrequency number j is analyzed:(35)δfij=fij−fi−1jfij100%,where fij is the estimation of eigenfrequency j using mesh i. The slowest convergence was observed for the last eigenfrequency (j = 10) for all shell computations. In Figure 5, the approximate error δf10 is plotted against the number of degrees of freedom (DoF) for SAM-H, CS, and MITC6 models. The convergence of both SAM-H and CS calculations is faster than that of MITC6 computations. In what follows, the mesh is selected so as to have values of eigenfrequencies with an approximate error less than 0.01%. In Table 1, the first ten eigenfrequencies calculated by SFE are listed. The eigenfrequencies calculated by SAM-H, CS, and MITC6 models and associated to the eigenmodes related to the SFE eigenfrequencies are also listed in this table. As compared to SFE results, a relative error less than 4% is obtained for all shell models considered herein. Since the structure is virtually a thick shell, all models yield good quality results for these eigenfrequencies. In Figure 6, the normalized displacement amplitude of the first mode is shown considering SFE, SAM-H, CS, and MITC6 calculations. Once again, the shell models provide similar eigenmodes to that obtained with SFE.Figure 5
Convergence of SAM-H, CS, and MITC6 models.Table 1
First 10 eigenfrequencies (in Hz) computed by SFE, SAM-H, CS, and MITC6.
j
SFE
SAM-H
CS
MITC6
1, 2
119.60
118.48
119.77
118.71
3, 4
122.98
121.77
122.82
122.82
5
141.84
140.09
142.38
142.67
6, 7
154.23
152.43
154.47
151.26
8
198.17
196.62
196.06
196.68
9, 10
198.32
196.36
196.60
191.72Figure 6
First eigenmode predicted by (a) SAM-H, (b) CS, (c) MITC6, and (d) SFE models.


### 4.2.2. Frequency Domain Analysis for the Case of Internal Pressure

The basis of the catenoid remains fixed and the upper edge is stress free. A harmonic pressurep is applied at the inner face of the shell (internal pressure); its amplitude p0 is 1 MPa, and its frequency f¯ is varied so as to perform a frequency domain analysis. The outer face of the shell is stress free. To reduce the computational cost of SFE calculations, the 3D problem is reduced to a 2D axisymmetric problem. In Figure 7, the amplitudes of U3 at the top of the catenoid calculated by shell models and SFE (at the middle surface) are shown. It can be seen that SAM-H is slightly more accurate than CS and MITC6 for this case. For all the models, a displacement peak is expected at a frequency close to 83 Hz.Figure 7
Displacement amplitude vs. excitation frequency curve predicted by SFE, SAM-H, CS, and MITC6 models.For safety reasons, the excitation frequencyf¯ should be less than 83 Hz. Let us consider the case when f¯ = 70 Hz and analyze the amplitude of stresses along the blue curve parallel to direction 1 shown in Figure 4(a). In Figures 8(a) and 8(b), the amplitudes of the in-plane normal stresses σ11 and σ22 at the internal face of the shell (ξ3=−h/2) are plotted against the ξ1 coordinate, respectively. Globally, the shape and maximum values of the SFE curves in Figures 8(a) and 8(b) are approximated better by SAM-H calculations. It is worth mentioning that σ22 exhibits the highest value among the other stress components in the curvilinear coordinate system. As compared to SFE, the relative errors in the maximum value of σ22 predicted by SAM-H, MITC6 and CS models are 2.8%, −9.7%, and −10.5%, respectively. SAM-H in-plane normal stresses are then more accurate than those of CS and MITC6 techniques.Figure 8
σ
11 vs. ξ1 (a) and σ22 vs. ξ1 (b) curves at the external face of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)In Figure9(a) the amplitude of the shear stress σ13 at the middle surface is plotted against the curvilinear coordinate ξ1. The evolution of the SFE σ13 stress along the line of curvature is approximated more accurately by SAM-H than MITC6 and CS models. The maximum value of the SFE shear stress is better approximated by SAM-H and MITC6 models than by the CS model. Let us now analyze the accuracy of SAM-H in the calculation of the out-of-plane normal stress σ33. It is worth recalling that CS and MITC6 techniques do not calculate this stress. At the inner and outer faces of the shell, since SAM-H verifies the 3D boundary conditions, SAM-H yields the amplitudes σ33=1MPa and σ33=0MPa, respectively. At the middle surface, the amplitudes of the stresses σ33 calculated by SFE and SAM-H are plotted against ξ1 in Figure 9(b). As in the case of the pressurized hollow sphere in Section 4.1, SAM-H yields good quality predictions of σ33.Figure 9
σ
13 vs. ξ1 (a) and σ33 vs. ξ1 (b) curves at the middle surface of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)


## 4.2.1. Eigenfrequencies and Convergence Analysis

The basis of the catenoid is fixed whilst its upper boundary is stress free. This means that for SAM-H and CS models:U1=U2=U3=0 and ϕ1=ϕ2=0 at ξ1=0; N11=N21=Q1=0 and M11=M21=0 at ξ1=a. No displacement boundary conditions are imposed at the inner and outer faces of the shell. For these shell models, in order to analyze the convergence of results, different meshes were tested; each mesh has a quasiuniform element size and in mesh number i, the element size is smaller than in mesh number i−1. The first ten eigenfrequencies were calculated, and the approximate error in the calculation of eigenfrequency number j is analyzed:(35)δfij=fij−fi−1jfij100%,where fij is the estimation of eigenfrequency j using mesh i. The slowest convergence was observed for the last eigenfrequency (j = 10) for all shell computations. In Figure 5, the approximate error δf10 is plotted against the number of degrees of freedom (DoF) for SAM-H, CS, and MITC6 models. The convergence of both SAM-H and CS calculations is faster than that of MITC6 computations. In what follows, the mesh is selected so as to have values of eigenfrequencies with an approximate error less than 0.01%. In Table 1, the first ten eigenfrequencies calculated by SFE are listed. The eigenfrequencies calculated by SAM-H, CS, and MITC6 models and associated to the eigenmodes related to the SFE eigenfrequencies are also listed in this table. As compared to SFE results, a relative error less than 4% is obtained for all shell models considered herein. Since the structure is virtually a thick shell, all models yield good quality results for these eigenfrequencies. In Figure 6, the normalized displacement amplitude of the first mode is shown considering SFE, SAM-H, CS, and MITC6 calculations. Once again, the shell models provide similar eigenmodes to that obtained with SFE.Figure 5
Convergence of SAM-H, CS, and MITC6 models.Table 1
First 10 eigenfrequencies (in Hz) computed by SFE, SAM-H, CS, and MITC6.
j
SFE
SAM-H
CS
MITC6
1, 2
119.60
118.48
119.77
118.71
3, 4
122.98
121.77
122.82
122.82
5
141.84
140.09
142.38
142.67
6, 7
154.23
152.43
154.47
151.26
8
198.17
196.62
196.06
196.68
9, 10
198.32
196.36
196.60
191.72Figure 6
First eigenmode predicted by (a) SAM-H, (b) CS, (c) MITC6, and (d) SFE models.


## 4.2.2. Frequency Domain Analysis for the Case of Internal Pressure

The basis of the catenoid remains fixed and the upper edge is stress free. A harmonic pressurep is applied at the inner face of the shell (internal pressure); its amplitude p0 is 1 MPa, and its frequency f¯ is varied so as to perform a frequency domain analysis. The outer face of the shell is stress free. To reduce the computational cost of SFE calculations, the 3D problem is reduced to a 2D axisymmetric problem. In Figure 7, the amplitudes of U3 at the top of the catenoid calculated by shell models and SFE (at the middle surface) are shown. It can be seen that SAM-H is slightly more accurate than CS and MITC6 for this case. For all the models, a displacement peak is expected at a frequency close to 83 Hz.Figure 7
Displacement amplitude vs. excitation frequency curve predicted by SFE, SAM-H, CS, and MITC6 models.For safety reasons, the excitation frequencyf¯ should be less than 83 Hz. Let us consider the case when f¯ = 70 Hz and analyze the amplitude of stresses along the blue curve parallel to direction 1 shown in Figure 4(a). In Figures 8(a) and 8(b), the amplitudes of the in-plane normal stresses σ11 and σ22 at the internal face of the shell (ξ3=−h/2) are plotted against the ξ1 coordinate, respectively. Globally, the shape and maximum values of the SFE curves in Figures 8(a) and 8(b) are approximated better by SAM-H calculations. It is worth mentioning that σ22 exhibits the highest value among the other stress components in the curvilinear coordinate system. As compared to SFE, the relative errors in the maximum value of σ22 predicted by SAM-H, MITC6 and CS models are 2.8%, −9.7%, and −10.5%, respectively. SAM-H in-plane normal stresses are then more accurate than those of CS and MITC6 techniques.Figure 8
σ
11 vs. ξ1 (a) and σ22 vs. ξ1 (b) curves at the external face of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)In Figure9(a) the amplitude of the shear stress σ13 at the middle surface is plotted against the curvilinear coordinate ξ1. The evolution of the SFE σ13 stress along the line of curvature is approximated more accurately by SAM-H than MITC6 and CS models. The maximum value of the SFE shear stress is better approximated by SAM-H and MITC6 models than by the CS model. Let us now analyze the accuracy of SAM-H in the calculation of the out-of-plane normal stress σ33. It is worth recalling that CS and MITC6 techniques do not calculate this stress. At the inner and outer faces of the shell, since SAM-H verifies the 3D boundary conditions, SAM-H yields the amplitudes σ33=1MPa and σ33=0MPa, respectively. At the middle surface, the amplitudes of the stresses σ33 calculated by SFE and SAM-H are plotted against ξ1 in Figure 9(b). As in the case of the pressurized hollow sphere in Section 4.1, SAM-H yields good quality predictions of σ33.Figure 9
σ
13 vs. ξ1 (a) and σ33 vs. ξ1 (b) curves at the middle surface of the shell obtained by SFE, CS, MITC6, and SAM-H models.
(a)
(b)


## 5. Conclusions

To conclude, a mixed stress/displacement approach model for dynamic problems of moderately thick shells was developed. It is an extension of the static SAM-H model developed by Domínguez Alvarado and Díaz in [1]; for this reason, the name SAM-H was preserved for the dynamic model. The dynamic model equations were obtained by introducing the Hellinger–Reissner functional into a Lagrangian and applying the Euler–Lagrange equation. The main goal of SAM-H is to determine accurately the 3D stresses: their approximation is inspired by the static version of SAM-H. The 3D displacement approximation is required to obtain the generalized motion equations: an approximation equivalent to that in first-order displacement approach models was chosen. The stress approximation verifies the 3D stress boundary conditions at the faces of the shell and the 3D motion equations. The constitutive equations take into account Poisson’s effect of out-of-plane normal stresses on in-plane strains; this is an important feature which is ignored in most shell models. Two structures were analyzed to test the accuracy of SAM-H: a hollow sphere and a catenoid. An analytical resolution was proposed for the study of the first structure whilst a finite element method was adopted to solve the model equations for the analysis of the second one. For both structures, eigenfrequencies were determined; after this, a harmonic internal pressure was applied and a frequency domain analysis was performed. SAM-H results were compared to those of MITC6 surface elements, the classic shell (CS) model, and solid finite elements (SFE). Both shell models and MITC6 elements yield accurate eigenfrequencies and eigenmodes even for a thick catenoid; their values are virtually the same and similar to those of SFE. Nevertheless, in a frequency domain analysis, the amplitudes of displacements and stresses of CS and MITC6 are not as accurate as those of SAM-H model which approximates SFE results with great precision.The numerical resolution of SAM-H equations does not imply a higher computational cost than that of the CS model because the same number of generalized displacements are considered. Besides, we have proved that the finite element calculations for the resolution of SAM-H equations converge faster than MITC6 calculations. The polynomial stress approximation in SAM-H may seem complex, but since the stress calculation should be performed in a postprocessing stage, it does not represent a computational drawback. For these reasons, SAM-H can be implemented in commercial finite element software for the dynamic analysis of homogeneous, moderately thick shells: calculation of eigenfrequencies and eigenmodes, frequency domain analysis, transient response, etc.



---

*Source: 5429615-2020-01-13.xml*