# Finite Element Model for Linear Elastic Thick Shells Using Gradient Recovery Method

Authors: Achille Germain Feumo; Robert Nzengwa; Joseph Nkongho Anyi
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This research purposed a new family of finite elements for spherical thick shell based on Nzengwa-Tagne’s model proposed in 1999. The model referred to hereafter as N-T model contains the classical Kirchhoff-Love (K-L) kinematic with additional terms related to the third fundamental form governing strain energy. Transverse shear stresses are computed andC0 finite element is proposed for numerical implementation. However, using straight line triangular elements does not guarantee a correct computation of stress across common edges of adjacent elements because of gradient jumps. The gradient recovery method known as Polynomial Preserving Recovery (PPR) is used for local interpolation and applied on a hemisphere under diametrically opposite charges. A good agreement of convergence results is observed; numerical results are compared to other results obtained with the classical K-L thin shell theory. Moreover, simulation on increasing values of the ratio of the shell shows impact of the N-T model especially on transverse stresses because of the significant energy contribution due to the third fundamental form tensor present in the kinematics of this model. The analysis of the thickness ratio shows difference between the classical K-L theory and N-T model when the ratio is greater than 0.099.


---

## Body

## 1. Introduction

N-T’s theoretical approach which was mathematically and rigorously deduced from three-dimensional linear elastic curvilinear media through multiple scaling and limit analysis is a more general Kirchhoff-Love (K-L) model. The displacement here is a two-degree polynomial of the thickness parameterz while the strain tensor which is planar contains the change in the third fundamental form in addition to the change of the first and second.For more than thirty years, numerous articles, books, and theses have addressed the problem of shell. Plates and more generally thin shells represent over 70% of industrial calculations [1, 2]. The mechanical models have been validated by some well-known benchmarks. Some locally stiffened thin shells or more generally thick shells have not received the same development probably because some few existing models do not account for transverse stresses and have not been mathematically established. Moreover it is well known that transverse stresses can not be neglected as the shell becomes thicker. Without any ad hoc geometrical hypothesis, it was deduced that the strain tensor is planar; that is, ∈i3=0 and, as a consequence, the limit displacement reads u=uαgα+u3g3, with uα(x,z)=ξα(x)-zθα(x)+z2ψα(x), u3=ξ3(x) (where {g1,g2,g3} is the three-dimensional contravariant shell basis); the in-plane strain ∈αβ=eαβ-zkαβ+z2Qαβ and the stresses σαβ depend on the in-plane strains and σi3 are computed from appropriate differential equations; see [3].The form of the displacement clearly shows thatC0 finite element is discontinuous across adjacent elements and usually provides inaccurate results at elements boundaries. Some authors indeed proposed different numerical methods: the finite difference method (FDM); finite volume method (FVM); finite element method for various simulations like magnetohydrodynamic (MHD) [4–6]; Curved Triangular Finite Elements (CTFE) of [7]; a weighted average method to calculate stresses, which provided good results for both interior and boundary elements; and L2-projection to also calculate stresses, by considering “discrete smoothing” and least squares fitting at the gauss points. But this last method was limited to one element; consequently, the smoothen stress is still discontinuous across element boundaries. This problem was completely solved 20 years ago, when Sheikholeslami et al. introduced their Superconvergent Patch Recovery (SPR) [8–10], where the discrete least square fitting was performed on an element patch, a set of elements having the same vertex. This SPR method produces a continuous stress field, which is superconvergent under uniform mesh. Soon after, Wiberg et al. incorporated equilibrium and boundary conditions to enhance SPR [11, 12] and discussed strategies to improve the finite element solution uh itself (other than the stress, which is essentially the gradient of uh) [12]. Recently, Naga and Zhang and their colleagues proposed an alternative strategy, called Polynomial Preserving Recovery (PPR) [13–17], to recover the gradient. Theoretical analyses revealed that PPR has better superconvergence (over any mesh) properties than SPR and the numerical tests indicated that the a posteriori error estimator based on PPR is as good as or better than that of SPR ([18, 19], remarks: page 323).In order to validate the N-T model, the program is first tested successfully on the widely known hemisphere under diametrically opposite charges of thin shells. Then thin shell theory of K-L and thick shell theory of N-T are implemented in order to evaluate the impact of thickness ratio. Next, simulations with various values of the characteristic shell parameter (thickness ratio) are implemented in order to reveal the contribution of Gauss curvature (change of the third fundamental form) in the stiffness energy.In addition to that, Section2 presents materials and methods: a brief description of the N-T model is presented without all the mathematical development, including gradient recovery method. A variational formulation of the shell equation and the resolution of transverse stresses equations are done. Section 3 is devoted to the discretization of the N-T model. Finite element spaces are next described and also the discretization scheme is layout. All numerical integrations are performed on a reference triangular element using the gradient recovery PPR method. Section 4 is completely devoted to the validation of the finite element model. Finally, results are discussed and concluded.


## 2. Materials and Methods

### 2.1. The N-T Model of Thick Shells

LetΩ=M∈Ω,OM→=om→+za3;-h/2≺z≺h/2,m∈S (S is midsurface, h≻0 is the thickness, and x is the coordinate of m in S) denote a shell. We assume the surface S is bounded and sufficiently smooth for all subsequent computations. Let {a1,a2,a3} and {a1,a2,a3} denote the covariant and contravariant basis of the midsurface and {g1,g2,g3} and {g1,g2,g3}, respectively, the covariant and contravariant basis of the shell. Then(1)gα=δατ-zbατaτ=ματaτ,g3=a3gα=μτα-1aτ,g3=a3dΩ=1-2zbαα+z2det⁡bβαdzdS=ϕx,zdzdS,where bαρ=aργbγα and bγα denote curvature tensor components and aργ is the contravariant component of the metric of the midsurface S. The repeated index convention is adopted. Values of Greek indices α,β take range in the set 1,2 while Latin indices i,j take their values in the set 1,2,3. A vector field can be expressed component wise indifferently in gi-basis or ai-basis as follows:(2)v=vix,zgi=v¯ix,zai,vα=ματv¯τ,v¯α=ματ-1vτ.Then the strain tensor (see [3]) is given by(3)ϵαβu=12uα/β+uβ/α=12ματ∇βu¯τ-bτβu¯3+μβτ∇αu¯τ-bατu¯3ϵ3αu=12uα/3+u3/α=12ματu¯τ,3+u¯3,α+bατu¯τϵ33u=u¯3,3(/ and ∇ indicate covariant derivation in Ω and S, resp., while f,α=∂f/∂xα). The equation ∈i3(u)=0 yields the following results:(4)u¯α=ματξτ-z∂αξ3u¯3=ξ3forξα,ξ3function  ofx=x1,x2,uα=ξα-z∂αξ3+2bατξτ+z2bντbανξτ+bατ∂αξ3u3=ξ3.The strain tensor now reads(5)ϵαβu=eαβξ-zkαβξ+z2Qαβξeαβξ=12∇αξβ+∇βξα-2bαβξ3kαβξ=∇αbβνξν+bαν∇βξν+bβν∇αξν+∇α∇βξ3-bαλbλβξ3Qαβξ=12bαν∇βbνφu˙φ+bανbνφ∇βu˙φ+bβνbνφ∇αu˙φ+bβν∇νbαφu˙φ+bαν∇β∇νu˙3+bβν∇α∇νu˙3.Let us recall thateαβ, kαβ, and Qαβ are, respectively, changes in the first, second, and third fundamental forms tensors, while the displacement (6)u=uαx,zgα+ξ3xa3=u¯αx,zaα+ξ3xa3is a more generalized Kirchhoff-Love displacement [3] and can also be written in Reissner-Mindlin format as follows:(7)u=ξαx-zθαx+z2ψαxgα+ξ3xa3=ξαx+zβαxaα+ξ3xa3,where βα=ξ3,α+2bανξν=θα+zbανξν and θα are the rotations angles. Strain tensors are miscalculated in some literatures because wrong basis vectors were used. The condition ∈i3(u)=0 modifies the constitutive law which, expressed with the Lamé constants, now reads(8)σαβu=λ¯gγδgαβ+μgαγgδβ+gαδgγβ∈γδu,where λ¯=2λμ/λ+2μ, or equivalently with Young’s modulus E and Poisson’s coefficient ν, reads(9)σαβu=E1-ν2νgγδgαβ+1-νgαγgδβ+gαδgγβ∈γδuConsider a shell of thicknessh, clamped on a part of its border Γ0, subject to volume forces fα and f3 and to surface forces h-α and h-3 on the rest of its border Γ1. Suppose the forces are sufficiently smooth; then the transverse stresses are solutions to the differential equations:(10)∂σα3∂z+2Γλ3ασλ3+Γλ3λσα3=-σ,βαβ+Γβλασλβ-fασα3-h2=-h¯-ασα3+h2=h¯+α,where fα∈L2(Ω); σα3∈H1(-h/2,h/2;H-1(S)) and(11)∂σ33∂z+Γα3ασ33=-σ,α3α+Γαλ3σλα+Γβλβσλ3-f3σ33-h2=-h¯-3σ33+h2=h¯+3,where f3∈H1(Ω); σ33∈H2(-h/2,h/2;H-2S).The generic point of the sphere is described by(12)OM→X=Rsin⁡ycos⁡xY=Rsin⁡xsin⁡yZ=Rcos⁡y,where x and y are curvilinear coordinates, 0≤x≤2Π, and 0≤y≤Π; R is the radius of the sphere, and X,Y,Z are the global coordinates.The covariant and contravariant metric tensors on the middle surfaceS are defined by(13)aαβ=R2sin2⁡y00R2aαβ=1R2sin2⁡y001R2.The covariant and contravariant curvature tensors onS (second fundamental form) are given by(14)bαβ=Rsin2⁡y00R;bαβ=1R3sin2⁡y001R3;bβα=1R001R.Leth be the thickness of the sphere; the Christoffel symbols are defined by(15)Γ111=Γ122=Γ212=Γ222=Γ221=0;Γ121=Γ211=cotany;Γ112=-12sin⁡2y;Γ231=Γ132=Γ123=Γ213=0;Γ131=Γ232=1h-R;Γ113=R-hsin2⁡y;Γ223=R-h.


### 2.2. Variational Equations

Let the border ofS, ∂S=γ0∪γ1 be partitioned in two parts and the border of the shell ∂Ω=Γ0∪Γ1 with Γ0=γ0×-h/2,h/2, and Γ1=γ1×-h/2,h/2∪Γ-∪Γ+. We denote Γ-=S×-h/2 and Γ+=S×h/2. Suppose the shell is clamped on Γ0 and subject to volume and surface forces as stated above; then the three-dimensional variational equation related to the equilibrium equation reads(16)findu∈IHΓ01such  that:∫Ωσu:ϵvdΩ=∫Ωf·vdΩ+∫Γ1h¯·vdS=Lv∀v∈IHΓ01,where IHΓ01Ω=uiΩ→R,ui∈L2Ω;∇jui∈L2Ωandui=0onΓ0 is the Sobolev space and: and · denote, respectively, tensors and vectors scalar products. The variational formulation under Kirchhoff-Love’s approach which is given by(17)A0u,v=Eh1-ν2∫SEut1-νDeβαtDeαβ+νDeααtDeββEv+h212Eut1-νDkβαtDkαβ+νDkααtDkββEvdS=Lv,(18)A1u,v=Eh1-ν2∫SEut1-νDeβαtDeαβ+νDeααtDeββEv+h212Eut1-νDkβαtDkαβ+νDkααtDkββEv+h212Eut1-νDeβαtDQαβ+νDeααtDQββEv+h212Eut1-νDQβαtDeαβ+νDQααtDeββEv+h480Eut1-νDQβαtDQαβ+νDQααtDQββEvdS=Lv˙is a truncated thick shell or the best first-order thick shell variational equation under N-T’s model; see [3]. It can be observed that this equation has similarities with familiar equations in engineering literature.Let(19)H1S=φ∈L2S,∇αφ∈L2SIH1S=H1S2IHγ01S=ηα∈IH1S,ηα=0onγ0H2S=φ∈H1S,∇αφ∈H1SHγ02S=φ∈H2S,φ=∂αφ=0onγ0Uad=IHγ01S×Hγ02S.Nzengwa and Tagne Simo established existence and unicity of solutions of this truncated problem inUad. We should remember that the displacement calculated after ξ is(20)uα=ξα-z∂αξ3+2bατξτ+z2bντbανξτ+bατ∂αξ3u3=ξ3or in the Nzengwa-Tagne format(21)uα=ξαx-zβαx+z2ψαxu3=ξ3xwhich suggests C0 finite element implementation.The linearized membrane strain tensor is expressed as follows:(22)eβαu=aαλeβλu=12aαλuβ,λ+uλ,β-aαλΓ¯λβρuρ-aαλbλβu3.Using the formulas(23)Tβαu=aαλTβλu=TβαEu,we can put (22) in the following form:(24)eβαu=DeβαEu,where(25)EuT=u1;∂1u1;∂2u1;u2;∂1u2;∂2u2;u3;-2b11u1-2b12u2-∂1u3;-2b11∂1u1-2b12∂1u2-2b11∂2u1-2b12∂2u2;-2b21u1-2b22u2-∂2u3;-2b21∂1u1-2b22∂1u2;-2b21∂2u1-2b22∂2u2;b¯1u1+b¯¯1u2+b11∂1u3+b12∂2u3;b¯1∂1u1+b¯¯1∂1u2;b¯1∂2u1+b¯¯1∂2u2;b¯2u1+b¯¯2u2+b21∂1u3+b22∂2u3;b¯2u1+b¯¯2u2+b21∂1u3+b22∂2u3;b¯2∂1u1+b¯¯2∂1u2;b¯2∂2u1+b¯¯2∂2u2,where(26)b¯1=b11b11+b21b12,b¯¯1=b12b11+b22b12,b¯2=b11b21+b21b22,b¯¯2=b12b21+b22b22;(27)Deβα=-aανΓβν1;aα1Iβ,1;12aαλJλ,β;-aαλΓβλ2;12aαλJλ,β;aα2Jβ,2;-bβα;0;0;0;0;0;0;0;0;0;0;0;0;I=1001,J=0110.The linearized part of the change of the curvature tensor of the middle surface is expressed as follows:(28)Dkβα=12aαμbμ,β1+Γ¯λβ1bμλ-Γ¯μβλbλ10;0;12aαμbμ,β2+Γ¯λβ2bμλ-Γ¯μβλbλ2;0;0;bαμbμβ;-aαμΓ¯μβ1;-Iα,1aα1;12Jα,βaαμ;-aαμΓ¯μβ2;12Jα,βaαμ;Iα,2aα2;0,0,0,0,0,0.The linearized part of the change of the third fundamental form tensor can be expressed as follows:(29)DQβα=aαμbν1bβ,μν+Γ¯γβνbμγ-Γ¯μβγbγν;0;0;aαμbν2bβ,μν+Γ¯γβνbμγ-Γ¯μβγbγν;0;0;0;bβ,μ1+Γ¯γβ1bμγ-Γ¯μβγbγ1;0;0;bβ,μ2+Γ¯γβ2bμγ-Γ¯μβγbγ2;0;0;-aαμΓ¯μβ1;Iα,1aα1;12Jα,βaαμ;-aαμΓ¯μβ2;12Jα,βaαμ;Iα,2aα2.


### 2.3. Resolution of the Transverse Stress Ordinary Differential Equation

Using expressions (27), (28), and (29), let τ=-4/h-R, Γjki be the three dimension Christoffel symbols, and(30)σ¯013=-E1-ν2a11e11+νe22,x1+1-νa11e12,x2+1-νcotanya11e12+1-νcotanya22e22.σ¯113=-E1-ν2a11K11+νk22,x1+1-νa11k12,x2+1-νcotanya11k12+1-νcotanya22k22σ¯213=-E1-ν2a11Q11+νQ22,x1+1-νa11Q1,x2+1-νcotanya11Q12+1-νcotanya22Q22Then the first shear stress of (10) is given as follows:(31)σ13x,y,z=-h--1exp⁡-τz+h2+1-exp⁡-τz+h2σ¯013τ+σ¯113τ2+2σ¯213τ3+z-h2-σ¯113τ-2σ¯213τ2+z2+h24σ¯213τ3Eu.In the same way, let(32)σ¯023=-E1-ν21-νa22e2,x12+a11e12,x2+12sin⁡2ya11e11+νe22σ¯123=-E1-ν21-νa22k2,x12+a11k12,x2+12sin⁡2ya11k11+νk22σ¯223=-E1-ν21-νa22Q2,x12+a11Q12,x2+12sin⁡2ya11Q11+νQ22.The second shear stress (10) is defined as follows:(33)σ23x,y,z=-h--2exp⁡-τz+h2+1-exp⁡-τz+h2σ¯023τ+σ¯123τ2+2σ¯223τ3+z-h2-σ¯123τ-2σ¯223τ2+z2+h24σ¯223τ3Eu.Let(34)γ1=σ¯0,x113+σ¯0,x223τ+σ¯1,x113+σ¯1,x223τ2+σ¯2,x113+σ¯2,x223τ3γ2=σ¯1,x113+σ¯1,x223τ+2σ¯2,x113+σ¯2,x223τ2γ3=2σ¯2,x113+σ¯2,x223τ3;ς=2h-R=Γ131+Γ232γ0=γ1+h2γ2-h44γ3γ¯1=Eh-Rsin2⁡y1-ν2a11e11+νe22-Ecotany1+νa11e12+Eh-R1-ν2a22e22+νe11γ¯2=Eh-Rsin2⁡y1-ν2a11k11+νk22-Ecotany1+νa11k12+Eh-R1-ν2a22k22+νk11γ¯3=Eh-Rsin2⁡y1-ν2a11Q11+νQ22-Ecotany1+νa11Q12+Eh-R1-ν2a22Q22+νQ11.Then the third expression of (11) is defined as follows:(35)σ33x,y,z=γ0ς-τ-γ2zς-τ-1ς-τ2+γ3z2ς-τ-2zς-τ2+2ς-τ3exp⁡-τz-h2+γ¯1-γ0ς+γ2-γ¯2zς-1ς2+γ¯3-γ3z2ς-2zς2+2ς3+exp⁡-τz-h2h-+3+γ0-γ¯1ς+γ¯2-γ2h2ς-1ς2+γ3-γ¯3h24ς-hς2+2ς3+γ0ς-τ+γ2h2ς-τ-1ς-τ2-γ3h24ς-τ-hς-τ2+2ς-τ3Eu.


### 2.4. Polynomial Preserving Recovery Method

Let us recall that the linear stress-strain relation does not guarantee smoothness across elements because of gradient jumps since the displacement is onlyC0. The shell is meshed with straight lines triangular elements obtained by linear transformation of a reference 2D triangle. Jumps across elements will be prevented by implementing the gradient recovery methods that we briefly present hereafter. Let uh be C0 finite element approximation of the solution u, and let(36)E1=I,J,KE2=I,J,Lbe two triangular elements with a common node I. We denote(37)U1=UhI,UhJ,UhK,U2=UhI,UhJ,UhLare the vector of nodal values in E1 and E2, respectively. Let Ni be the shape functions matrix; then restriction of uh in each element reads(38)uh∣E1=N1U1t,uh∣E2=N2U2t,uh∣E1′=N1′U1t,uh∣E2′=N2′U2t,where Ni′ stands for the derivative according to x or y of Ni. At the common node I, the derivative of uh, uh′(I), is not necessarily the same in both elements. This means that on a patch including a node I common to different elements the gradient of uh, ∇uh, cannot be calculated because of the jump on each part of the patch. Consequently ∇uh is not a good approximation of ∇u across elements and stress and strain calculated using ∇u cannot be continuous across edges. How to define a unique ∇uh at a node I common to different elements has been addressed by gradient recovery methods. The idea is to define a local operator Gh such that Ghuh(I) is unique through any patch and Ghuh-∇u has better approximation than ∇uh-∇u. Consequently convergence issues of these methods should be considered also. In this analysis we consider the 2D PPR gradient recovery methods which also guarantees a superconvergence property of Ghuh to ∇u independently of mesh size. The method is described as follows.LetZi be a node where ∇u is to be determined. Let Zij be the nodes of all triangles Ej having Zi as common vertex. Suppose uh|Ej∈Pk+1, the set of polynomials of order k+1; then the gradient recovery PPR operator consists in defining p(x)∈Pk+1 such that(39)∑j=1npx-uh2Zij=minqx∈Pk+1⁡∑j=1nqx-uh2ZijGhuhZi=∇pZi.This value also depends on the sampling points chosen. In the 2D case it is proved (cf. [14]) that if n≧m=(k+2)(k+3)/2 and if the sum of two adjacent angles in the mesh is not more than π, then Ghuh(Zi) is unique for any Zi. In this work the angle condition will be implemented in the triangularization. The number m of sampling nodes will be respected as follows:(40)GhuhZi=∇pxi,yi=∇p~0,0≃∇uxi,yi.For the 2D case, three types of nodes can be distinguished: internal nodes, boundary nodes at a corner, and boundary nodes out of a corner.(i)
For the in-plane displacementuα(α=1;2), the existence and uniqueness solution are possible if (u1,u2)∈H2(S)2; we fit linear polynômial using the same regular pattern:(41)P1x,y=1,x,ya0,a1,a3t=1,ξ,ςa^0,a^1,a^2t.We scale by a factorh with x=hξ and y=hς, (ξ,ς) with respect to the six derivative values at the barycentric center of each element on the patch. Now we define e=1,1,1,1,1,1, a^=a^0,a^1,a^2=a0,ha1,ha3, A=e,ξ→,ς→, H=diag⁡1,h,h.Let uαi be the values of uα at the nodal points ZIii=1,2,3,4,5,6; we denote u¯=uαii=16; we hereby determine a^ so that Aa^=b¯h.Then, we obtain(42)P1x,y=1,x,yH-1AtA-1Atb¯h;notice that(43)∂uα∂xx,y=∂P1∂xx,y=0,1,0H-1AtA-1Atb¯h=P1Xα1b¯h∂uα∂yx,y=∂P1∂yx,y=0,0,1H-1AtA-1Atb¯h=P1Xα2b¯h.For the transversal displacement,u3, the existence and uniqueness of a gradient recovery solution are possible if (u3)∈H2(S); we fit quadratic polynômial using the same regular pattern:(44)P2x,y=1,x,y,x2,xy,y2a1,a2,a3,a4,a5,a6t=1,ξ,ς,ξ2,ξς,ς2a^1,a^2,a^3,a^4,a^5,a^6t;then, we scale by a factor h with x=hξ and y=hς, (ξ,ς) with respect to the six derivative values at the coordinates nodes of each element on the patch. Let e=(1,1,1,1,1,1,1), a^=a^1,a^2,a^3,a^4,a^5,a^6=a1,ha2,ha3,h2a4,h2a5,h2a6 be the vectors of the matrix B=(e,ξ→,ς→,ξ2→,ξς→,ς2→) and consider the diagonal matrix H1=diag⁡(1,h,h,h2,h2,h2) and also b¯h=(u30,u31,u32,u33,u34,u35,u36) the approximation vector of u3 in a nodal point ZI. Then we fit P2(x,y)=1,x,y,x2,xy,y2H1-1BtB-1Btb¯h and obtain the recovered gradient at the patch center point as follows:(45)∂u3∂xx,y=∂P2∂xx,y=0,1,0,2x,y,0H1-1BtB-1Btb¯h=P2X31b¯h∂u3∂yx,y=∂P2∂yx,y=0,0,1,0,x,2yH1-1BtB-1Btb¯h=P2X32b¯h∂2u3∂x2x,y=∂2P2∂x2x,y=0,1,0,2,0,0H1-1BtB-1Btb¯h=P2X311b¯h∂2u3∂x∂yx,y=∂2P2∂x∂yx,y=0,0,0,0,1,0H1-1BtB-1Btb¯h=P2X312b¯h∂2u3∂y2x,y=∂2P2∂y2x,y=0,0,0,0,0,2H1-1BtB-1Btb¯h=P2X322b¯h.WithGhui given at each vertex by the same processes in (43) and (45), we are able to form a recovered gradient field by using the finite element basis functions. Recovering the gradient at a boundary vertex is more delicate. Using efficient strategy computational experiment indicated in [14], to recover the gradient at a vertex zI∈∂S, we look for the nearest layer of vertices around zI that contain at least one internal vertex. Let this layer be the nth one and denote the internal vertices in this layer by z1,z2,…,zm, where m≥1. The union of the sampling points used in recovering the gradient at z1,z2,…,zm and the mesh nodes in the first n layers around zI constitute the set of sampling point for recovering the gradient at zI; see [15].


## 2.1. The N-T Model of Thick Shells

LetΩ=M∈Ω,OM→=om→+za3;-h/2≺z≺h/2,m∈S (S is midsurface, h≻0 is the thickness, and x is the coordinate of m in S) denote a shell. We assume the surface S is bounded and sufficiently smooth for all subsequent computations. Let {a1,a2,a3} and {a1,a2,a3} denote the covariant and contravariant basis of the midsurface and {g1,g2,g3} and {g1,g2,g3}, respectively, the covariant and contravariant basis of the shell. Then(1)gα=δατ-zbατaτ=ματaτ,g3=a3gα=μτα-1aτ,g3=a3dΩ=1-2zbαα+z2det⁡bβαdzdS=ϕx,zdzdS,where bαρ=aργbγα and bγα denote curvature tensor components and aργ is the contravariant component of the metric of the midsurface S. The repeated index convention is adopted. Values of Greek indices α,β take range in the set 1,2 while Latin indices i,j take their values in the set 1,2,3. A vector field can be expressed component wise indifferently in gi-basis or ai-basis as follows:(2)v=vix,zgi=v¯ix,zai,vα=ματv¯τ,v¯α=ματ-1vτ.Then the strain tensor (see [3]) is given by(3)ϵαβu=12uα/β+uβ/α=12ματ∇βu¯τ-bτβu¯3+μβτ∇αu¯τ-bατu¯3ϵ3αu=12uα/3+u3/α=12ματu¯τ,3+u¯3,α+bατu¯τϵ33u=u¯3,3(/ and ∇ indicate covariant derivation in Ω and S, resp., while f,α=∂f/∂xα). The equation ∈i3(u)=0 yields the following results:(4)u¯α=ματξτ-z∂αξ3u¯3=ξ3forξα,ξ3function  ofx=x1,x2,uα=ξα-z∂αξ3+2bατξτ+z2bντbανξτ+bατ∂αξ3u3=ξ3.The strain tensor now reads(5)ϵαβu=eαβξ-zkαβξ+z2Qαβξeαβξ=12∇αξβ+∇βξα-2bαβξ3kαβξ=∇αbβνξν+bαν∇βξν+bβν∇αξν+∇α∇βξ3-bαλbλβξ3Qαβξ=12bαν∇βbνφu˙φ+bανbνφ∇βu˙φ+bβνbνφ∇αu˙φ+bβν∇νbαφu˙φ+bαν∇β∇νu˙3+bβν∇α∇νu˙3.Let us recall thateαβ, kαβ, and Qαβ are, respectively, changes in the first, second, and third fundamental forms tensors, while the displacement (6)u=uαx,zgα+ξ3xa3=u¯αx,zaα+ξ3xa3is a more generalized Kirchhoff-Love displacement [3] and can also be written in Reissner-Mindlin format as follows:(7)u=ξαx-zθαx+z2ψαxgα+ξ3xa3=ξαx+zβαxaα+ξ3xa3,where βα=ξ3,α+2bανξν=θα+zbανξν and θα are the rotations angles. Strain tensors are miscalculated in some literatures because wrong basis vectors were used. The condition ∈i3(u)=0 modifies the constitutive law which, expressed with the Lamé constants, now reads(8)σαβu=λ¯gγδgαβ+μgαγgδβ+gαδgγβ∈γδu,where λ¯=2λμ/λ+2μ, or equivalently with Young’s modulus E and Poisson’s coefficient ν, reads(9)σαβu=E1-ν2νgγδgαβ+1-νgαγgδβ+gαδgγβ∈γδuConsider a shell of thicknessh, clamped on a part of its border Γ0, subject to volume forces fα and f3 and to surface forces h-α and h-3 on the rest of its border Γ1. Suppose the forces are sufficiently smooth; then the transverse stresses are solutions to the differential equations:(10)∂σα3∂z+2Γλ3ασλ3+Γλ3λσα3=-σ,βαβ+Γβλασλβ-fασα3-h2=-h¯-ασα3+h2=h¯+α,where fα∈L2(Ω); σα3∈H1(-h/2,h/2;H-1(S)) and(11)∂σ33∂z+Γα3ασ33=-σ,α3α+Γαλ3σλα+Γβλβσλ3-f3σ33-h2=-h¯-3σ33+h2=h¯+3,where f3∈H1(Ω); σ33∈H2(-h/2,h/2;H-2S).The generic point of the sphere is described by(12)OM→X=Rsin⁡ycos⁡xY=Rsin⁡xsin⁡yZ=Rcos⁡y,where x and y are curvilinear coordinates, 0≤x≤2Π, and 0≤y≤Π; R is the radius of the sphere, and X,Y,Z are the global coordinates.The covariant and contravariant metric tensors on the middle surfaceS are defined by(13)aαβ=R2sin2⁡y00R2aαβ=1R2sin2⁡y001R2.The covariant and contravariant curvature tensors onS (second fundamental form) are given by(14)bαβ=Rsin2⁡y00R;bαβ=1R3sin2⁡y001R3;bβα=1R001R.Leth be the thickness of the sphere; the Christoffel symbols are defined by(15)Γ111=Γ122=Γ212=Γ222=Γ221=0;Γ121=Γ211=cotany;Γ112=-12sin⁡2y;Γ231=Γ132=Γ123=Γ213=0;Γ131=Γ232=1h-R;Γ113=R-hsin2⁡y;Γ223=R-h.


## 2.2. Variational Equations

Let the border ofS, ∂S=γ0∪γ1 be partitioned in two parts and the border of the shell ∂Ω=Γ0∪Γ1 with Γ0=γ0×-h/2,h/2, and Γ1=γ1×-h/2,h/2∪Γ-∪Γ+. We denote Γ-=S×-h/2 and Γ+=S×h/2. Suppose the shell is clamped on Γ0 and subject to volume and surface forces as stated above; then the three-dimensional variational equation related to the equilibrium equation reads(16)findu∈IHΓ01such  that:∫Ωσu:ϵvdΩ=∫Ωf·vdΩ+∫Γ1h¯·vdS=Lv∀v∈IHΓ01,where IHΓ01Ω=uiΩ→R,ui∈L2Ω;∇jui∈L2Ωandui=0onΓ0 is the Sobolev space and: and · denote, respectively, tensors and vectors scalar products. The variational formulation under Kirchhoff-Love’s approach which is given by(17)A0u,v=Eh1-ν2∫SEut1-νDeβαtDeαβ+νDeααtDeββEv+h212Eut1-νDkβαtDkαβ+νDkααtDkββEvdS=Lv,(18)A1u,v=Eh1-ν2∫SEut1-νDeβαtDeαβ+νDeααtDeββEv+h212Eut1-νDkβαtDkαβ+νDkααtDkββEv+h212Eut1-νDeβαtDQαβ+νDeααtDQββEv+h212Eut1-νDQβαtDeαβ+νDQααtDeββEv+h480Eut1-νDQβαtDQαβ+νDQααtDQββEvdS=Lv˙is a truncated thick shell or the best first-order thick shell variational equation under N-T’s model; see [3]. It can be observed that this equation has similarities with familiar equations in engineering literature.Let(19)H1S=φ∈L2S,∇αφ∈L2SIH1S=H1S2IHγ01S=ηα∈IH1S,ηα=0onγ0H2S=φ∈H1S,∇αφ∈H1SHγ02S=φ∈H2S,φ=∂αφ=0onγ0Uad=IHγ01S×Hγ02S.Nzengwa and Tagne Simo established existence and unicity of solutions of this truncated problem inUad. We should remember that the displacement calculated after ξ is(20)uα=ξα-z∂αξ3+2bατξτ+z2bντbανξτ+bατ∂αξ3u3=ξ3or in the Nzengwa-Tagne format(21)uα=ξαx-zβαx+z2ψαxu3=ξ3xwhich suggests C0 finite element implementation.The linearized membrane strain tensor is expressed as follows:(22)eβαu=aαλeβλu=12aαλuβ,λ+uλ,β-aαλΓ¯λβρuρ-aαλbλβu3.Using the formulas(23)Tβαu=aαλTβλu=TβαEu,we can put (22) in the following form:(24)eβαu=DeβαEu,where(25)EuT=u1;∂1u1;∂2u1;u2;∂1u2;∂2u2;u3;-2b11u1-2b12u2-∂1u3;-2b11∂1u1-2b12∂1u2-2b11∂2u1-2b12∂2u2;-2b21u1-2b22u2-∂2u3;-2b21∂1u1-2b22∂1u2;-2b21∂2u1-2b22∂2u2;b¯1u1+b¯¯1u2+b11∂1u3+b12∂2u3;b¯1∂1u1+b¯¯1∂1u2;b¯1∂2u1+b¯¯1∂2u2;b¯2u1+b¯¯2u2+b21∂1u3+b22∂2u3;b¯2u1+b¯¯2u2+b21∂1u3+b22∂2u3;b¯2∂1u1+b¯¯2∂1u2;b¯2∂2u1+b¯¯2∂2u2,where(26)b¯1=b11b11+b21b12,b¯¯1=b12b11+b22b12,b¯2=b11b21+b21b22,b¯¯2=b12b21+b22b22;(27)Deβα=-aανΓβν1;aα1Iβ,1;12aαλJλ,β;-aαλΓβλ2;12aαλJλ,β;aα2Jβ,2;-bβα;0;0;0;0;0;0;0;0;0;0;0;0;I=1001,J=0110.The linearized part of the change of the curvature tensor of the middle surface is expressed as follows:(28)Dkβα=12aαμbμ,β1+Γ¯λβ1bμλ-Γ¯μβλbλ10;0;12aαμbμ,β2+Γ¯λβ2bμλ-Γ¯μβλbλ2;0;0;bαμbμβ;-aαμΓ¯μβ1;-Iα,1aα1;12Jα,βaαμ;-aαμΓ¯μβ2;12Jα,βaαμ;Iα,2aα2;0,0,0,0,0,0.The linearized part of the change of the third fundamental form tensor can be expressed as follows:(29)DQβα=aαμbν1bβ,μν+Γ¯γβνbμγ-Γ¯μβγbγν;0;0;aαμbν2bβ,μν+Γ¯γβνbμγ-Γ¯μβγbγν;0;0;0;bβ,μ1+Γ¯γβ1bμγ-Γ¯μβγbγ1;0;0;bβ,μ2+Γ¯γβ2bμγ-Γ¯μβγbγ2;0;0;-aαμΓ¯μβ1;Iα,1aα1;12Jα,βaαμ;-aαμΓ¯μβ2;12Jα,βaαμ;Iα,2aα2.


## 2.3. Resolution of the Transverse Stress Ordinary Differential Equation

Using expressions (27), (28), and (29), let τ=-4/h-R, Γjki be the three dimension Christoffel symbols, and(30)σ¯013=-E1-ν2a11e11+νe22,x1+1-νa11e12,x2+1-νcotanya11e12+1-νcotanya22e22.σ¯113=-E1-ν2a11K11+νk22,x1+1-νa11k12,x2+1-νcotanya11k12+1-νcotanya22k22σ¯213=-E1-ν2a11Q11+νQ22,x1+1-νa11Q1,x2+1-νcotanya11Q12+1-νcotanya22Q22Then the first shear stress of (10) is given as follows:(31)σ13x,y,z=-h--1exp⁡-τz+h2+1-exp⁡-τz+h2σ¯013τ+σ¯113τ2+2σ¯213τ3+z-h2-σ¯113τ-2σ¯213τ2+z2+h24σ¯213τ3Eu.In the same way, let(32)σ¯023=-E1-ν21-νa22e2,x12+a11e12,x2+12sin⁡2ya11e11+νe22σ¯123=-E1-ν21-νa22k2,x12+a11k12,x2+12sin⁡2ya11k11+νk22σ¯223=-E1-ν21-νa22Q2,x12+a11Q12,x2+12sin⁡2ya11Q11+νQ22.The second shear stress (10) is defined as follows:(33)σ23x,y,z=-h--2exp⁡-τz+h2+1-exp⁡-τz+h2σ¯023τ+σ¯123τ2+2σ¯223τ3+z-h2-σ¯123τ-2σ¯223τ2+z2+h24σ¯223τ3Eu.Let(34)γ1=σ¯0,x113+σ¯0,x223τ+σ¯1,x113+σ¯1,x223τ2+σ¯2,x113+σ¯2,x223τ3γ2=σ¯1,x113+σ¯1,x223τ+2σ¯2,x113+σ¯2,x223τ2γ3=2σ¯2,x113+σ¯2,x223τ3;ς=2h-R=Γ131+Γ232γ0=γ1+h2γ2-h44γ3γ¯1=Eh-Rsin2⁡y1-ν2a11e11+νe22-Ecotany1+νa11e12+Eh-R1-ν2a22e22+νe11γ¯2=Eh-Rsin2⁡y1-ν2a11k11+νk22-Ecotany1+νa11k12+Eh-R1-ν2a22k22+νk11γ¯3=Eh-Rsin2⁡y1-ν2a11Q11+νQ22-Ecotany1+νa11Q12+Eh-R1-ν2a22Q22+νQ11.Then the third expression of (11) is defined as follows:(35)σ33x,y,z=γ0ς-τ-γ2zς-τ-1ς-τ2+γ3z2ς-τ-2zς-τ2+2ς-τ3exp⁡-τz-h2+γ¯1-γ0ς+γ2-γ¯2zς-1ς2+γ¯3-γ3z2ς-2zς2+2ς3+exp⁡-τz-h2h-+3+γ0-γ¯1ς+γ¯2-γ2h2ς-1ς2+γ3-γ¯3h24ς-hς2+2ς3+γ0ς-τ+γ2h2ς-τ-1ς-τ2-γ3h24ς-τ-hς-τ2+2ς-τ3Eu.


## 2.4. Polynomial Preserving Recovery Method

Let us recall that the linear stress-strain relation does not guarantee smoothness across elements because of gradient jumps since the displacement is onlyC0. The shell is meshed with straight lines triangular elements obtained by linear transformation of a reference 2D triangle. Jumps across elements will be prevented by implementing the gradient recovery methods that we briefly present hereafter. Let uh be C0 finite element approximation of the solution u, and let(36)E1=I,J,KE2=I,J,Lbe two triangular elements with a common node I. We denote(37)U1=UhI,UhJ,UhK,U2=UhI,UhJ,UhLare the vector of nodal values in E1 and E2, respectively. Let Ni be the shape functions matrix; then restriction of uh in each element reads(38)uh∣E1=N1U1t,uh∣E2=N2U2t,uh∣E1′=N1′U1t,uh∣E2′=N2′U2t,where Ni′ stands for the derivative according to x or y of Ni. At the common node I, the derivative of uh, uh′(I), is not necessarily the same in both elements. This means that on a patch including a node I common to different elements the gradient of uh, ∇uh, cannot be calculated because of the jump on each part of the patch. Consequently ∇uh is not a good approximation of ∇u across elements and stress and strain calculated using ∇u cannot be continuous across edges. How to define a unique ∇uh at a node I common to different elements has been addressed by gradient recovery methods. The idea is to define a local operator Gh such that Ghuh(I) is unique through any patch and Ghuh-∇u has better approximation than ∇uh-∇u. Consequently convergence issues of these methods should be considered also. In this analysis we consider the 2D PPR gradient recovery methods which also guarantees a superconvergence property of Ghuh to ∇u independently of mesh size. The method is described as follows.LetZi be a node where ∇u is to be determined. Let Zij be the nodes of all triangles Ej having Zi as common vertex. Suppose uh|Ej∈Pk+1, the set of polynomials of order k+1; then the gradient recovery PPR operator consists in defining p(x)∈Pk+1 such that(39)∑j=1npx-uh2Zij=minqx∈Pk+1⁡∑j=1nqx-uh2ZijGhuhZi=∇pZi.This value also depends on the sampling points chosen. In the 2D case it is proved (cf. [14]) that if n≧m=(k+2)(k+3)/2 and if the sum of two adjacent angles in the mesh is not more than π, then Ghuh(Zi) is unique for any Zi. In this work the angle condition will be implemented in the triangularization. The number m of sampling nodes will be respected as follows:(40)GhuhZi=∇pxi,yi=∇p~0,0≃∇uxi,yi.For the 2D case, three types of nodes can be distinguished: internal nodes, boundary nodes at a corner, and boundary nodes out of a corner.(i)
For the in-plane displacementuα(α=1;2), the existence and uniqueness solution are possible if (u1,u2)∈H2(S)2; we fit linear polynômial using the same regular pattern:(41)P1x,y=1,x,ya0,a1,a3t=1,ξ,ςa^0,a^1,a^2t.We scale by a factorh with x=hξ and y=hς, (ξ,ς) with respect to the six derivative values at the barycentric center of each element on the patch. Now we define e=1,1,1,1,1,1, a^=a^0,a^1,a^2=a0,ha1,ha3, A=e,ξ→,ς→, H=diag⁡1,h,h.Let uαi be the values of uα at the nodal points ZIii=1,2,3,4,5,6; we denote u¯=uαii=16; we hereby determine a^ so that Aa^=b¯h.Then, we obtain(42)P1x,y=1,x,yH-1AtA-1Atb¯h;notice that(43)∂uα∂xx,y=∂P1∂xx,y=0,1,0H-1AtA-1Atb¯h=P1Xα1b¯h∂uα∂yx,y=∂P1∂yx,y=0,0,1H-1AtA-1Atb¯h=P1Xα2b¯h.For the transversal displacement,u3, the existence and uniqueness of a gradient recovery solution are possible if (u3)∈H2(S); we fit quadratic polynômial using the same regular pattern:(44)P2x,y=1,x,y,x2,xy,y2a1,a2,a3,a4,a5,a6t=1,ξ,ς,ξ2,ξς,ς2a^1,a^2,a^3,a^4,a^5,a^6t;then, we scale by a factor h with x=hξ and y=hς, (ξ,ς) with respect to the six derivative values at the coordinates nodes of each element on the patch. Let e=(1,1,1,1,1,1,1), a^=a^1,a^2,a^3,a^4,a^5,a^6=a1,ha2,ha3,h2a4,h2a5,h2a6 be the vectors of the matrix B=(e,ξ→,ς→,ξ2→,ξς→,ς2→) and consider the diagonal matrix H1=diag⁡(1,h,h,h2,h2,h2) and also b¯h=(u30,u31,u32,u33,u34,u35,u36) the approximation vector of u3 in a nodal point ZI. Then we fit P2(x,y)=1,x,y,x2,xy,y2H1-1BtB-1Btb¯h and obtain the recovered gradient at the patch center point as follows:(45)∂u3∂xx,y=∂P2∂xx,y=0,1,0,2x,y,0H1-1BtB-1Btb¯h=P2X31b¯h∂u3∂yx,y=∂P2∂yx,y=0,0,1,0,x,2yH1-1BtB-1Btb¯h=P2X32b¯h∂2u3∂x2x,y=∂2P2∂x2x,y=0,1,0,2,0,0H1-1BtB-1Btb¯h=P2X311b¯h∂2u3∂x∂yx,y=∂2P2∂x∂yx,y=0,0,0,0,1,0H1-1BtB-1Btb¯h=P2X312b¯h∂2u3∂y2x,y=∂2P2∂y2x,y=0,0,0,0,0,2H1-1BtB-1Btb¯h=P2X322b¯h.WithGhui given at each vertex by the same processes in (43) and (45), we are able to form a recovered gradient field by using the finite element basis functions. Recovering the gradient at a boundary vertex is more delicate. Using efficient strategy computational experiment indicated in [14], to recover the gradient at a vertex zI∈∂S, we look for the nearest layer of vertices around zI that contain at least one internal vertex. Let this layer be the nth one and denote the internal vertices in this layer by z1,z2,…,zm, where m≥1. The union of the sampling points used in recovering the gradient at z1,z2,…,zm and the mesh nodes in the first n layers around zI constitute the set of sampling point for recovering the gradient at zI; see [15].


## 3. Discretization of N-T Model

### 3.1. Finite Element Space

The continuous truncated variational equation (18) is defined in the space Uad. In order to define a finite element space, we begin by recalling the following results. Let Tj be a triangle with nodes a1=(x1,y1), a2=(x2,y2), and a3=(x3,y3). Let P1 and P2 be the sets of the first- and second-order polynomials with basis 1,x,y and 1,x,y,x2,y2,xy, respectively, then λ1,λ2,λ3 and λ1,λ2,λ3,4λ1λ2,4λ1λ3,4λ2λ3 also generate P1 and P2 where the barycentric basis functions of the triangle are defined by(46)λ1x,y=12∆y3-y2x2-x-x3-x2y2-yλ2x,y=12∆y1-y3x3-x-x1-x3y3-yλ3x,y=12∆y2-y1x1-x-x2-x1y1-y2∆=2∫Tjdxdy=2det⁡111x1x2x3y1y2y3=x3-x2y1-y2-x1-x2y3-y2.LetTh be a triangularization of the midsurface of the shell and nh the number of triangles. We denote by(47)Xh1=vh∈C0S¯,vh∣Tj∈P1Tj,∀Tj,j=1,…,nhthe subspace of dimension Nh (which is the number of nodes) and(48)Xh2=vh∈C0S¯,vh∣Tj∈P2Tj,∀Tj,j=1,…,nh,∇vh=Ghvh=∇P2;S¯ is the closure of S, the subspace of dimension Nh2 which is the number of unisolvent points (nodes and edge mid points).The spaceUad=HΓ01S2×HΓ02(S) is approximated with the finite element space of dimension 3N=2Nh+Nh2(49)℧h=Xh1×Xh1×Xh2.


### 3.2. The Discrete Scheme

LetTh be a triangularization of the shell’s midsurface S and Tj∈Th be a triangle of vertices a1,a2,a3 and midedges a4,a5,a6. Let ηh∈℧h as defined above. Then(50)ηαhx,yTj∈P1,α=1,2;1≤j≤nh,η3hx,yTj∈P2,1≤j≤nh.By using the barycentric polynomials and writingλ4=λ1, 1≤l≤2; 1≤α, β≤2, and 1≤j≤nh,(51)ηαhx,yTj=∑k=13ηαhakλkx,y=∑k=13η-αkλkx,yη3hx,yTj=∑k=13η3hakλkx,y+∑k=46η3hakλk-3x,yλk-2x,y=∑k=13η-3kλkx,y+∑k=46η-3kλk-3x,yλk-2x,y∂αηlhx,yTj=∑k=13∂αη3hakλkx,y,∂αηlhak=Ghηlhak∂αβ2η3hx,yTj=∑k=13∂αβ2η3hakλkx,y,∂αβ2η3hak=Gh∂αη3hak.In the above formula,Gh is the gradient operator and η-34, η-35, and η-36 are, respectively, the unknown values of η3 in the midedges a1a2, a1a3, and a2a3. Let(i)
PνIαβ=PνXαβ and P2Iααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is internal of the mesh,(ii)
PνCαβ=PνXαβ and P2Cααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is in the boundary at the corner of the mesh,(iii)
PνMαβ=PνXαβ and P2Mααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is in the boundary not at the corner of the mesh.Exploiting the above, let us take linear element on uniform triangular mesh of a regular pattern; we hereby investigate the case of the element where all three nodes are internal of the mesh; the matrices deduced from the discretization of the gradient are given as follows:(52)A1=λ1x1,y10000000λ2x1+h,y100000000000λ3x1,y1-htA2=0000λ1x2-h,y200λ2x2,y200000000000λ3x2-h,y2-h0tA3=000λ1x3,y3+h00000λ2x3+h,y3+h0000λ3x3,y3000000t.When we return to energy function withm=0,1, we have(53)Amu¯,v¯=∫SEu¯TCgmEv¯dS=∑j=1nh∑k=13∫TjEu¯kTjtCgmEv¯kTjdxdy.Using (43), (45), and (51) in which all three nodes of the triangular element are internal,  let  i=1,2,3; α,β=1,2 and we consider (54)αi=λix,yρk=4λk-3x,yλk-2x,y,k=4,5,6βiα=P1IiαA1:,i+P1IiαA2:,i+P1IiαA3:,iγiα=P2I3αA1:,i+P2I3αA2:,i+P2I3αA3:,iτiαβ=P2I3αβA1:,i+P2I3αβA2:,i+P2I3αβA3:,i.Then the deformation vector can be written asEuTj=BTjduTj where(55)BTjt=α1α2α3000000000β11β21β31000000000β12β12β12000000000000α1α2α3000000000β11β21β31000000000β12β22β32000000000000α1α2α3ρ4ρ5ρ6-2b11α1-2b11α2-2b11α3-2b12α1-2b12α2-2b12α3-γ11-γ21-γ31000-2b11β11-2b11β21-2b11β31-2b12β11-2b12β21-2b12β31-τ111-τ211-τ311000-2b11β12-2b11β22-2b11β32-2b12β12-2b12β22-2b12β32-τ112-τ212-τ312000-2b21α1-2b21α2-2b21α3-2b22α1-2b22α2-2b22α3-γ12-γ22-γ32000-2b21β11-2b21β21-2b21β31-2b22β11-2b22β21-2b22β31-τ112-τ212-τ312000-2b21β12-2b21β22-2b21β32-2b22β12-2b22β22-2b22β32-τ122-τ222-τ322000b¯1α1b¯1α2b¯1α3b¯¯1α1b¯¯1α2b¯¯1α3b11γ11+b12γ12b11γ21+b12γ22b11γ31+b12γ32000b¯1β11b¯1β21b¯1β21b¯¯1β11b¯¯1β21b¯¯1β31b11τ111+b12τ112b11τ211+b12τ212b11τ311+b12τ312000b¯1β12b¯1β22b¯1β22b¯¯1β12b¯¯1β22b¯¯1β32b11τ112+b12τ122b11τ212+b12τ222b11τ312+b12τ322000b¯2α1b¯2α2b¯2α3b¯¯2α1b¯¯2α2b¯¯2α3b21γ11+b22γ12b21γ21+b22γ22b21γ31+b22γ32000b¯2β11b¯2β21b¯2β31b¯¯2β11b¯¯2β21b¯¯2β31b21τ111+b22τ112b21τ211+b22τ212b21τ311+b22τ312000b¯2β12b¯2β22b¯2β32b¯¯2β12b¯¯2β22b¯¯2β32b21τ112+b22τ122b21τ212+b22τ222b21τ312+b22τ322000,(56)du¯Tj=u¯1j1,u¯1j2,u¯1j3,u¯2j1,u¯2j2,u¯2j3,u¯3j1,u¯3j2,u¯3j3,u¯3j4,u¯3j5,u¯3j6t.


### 3.3. Stiffness Matrix

Using local element, stiffness matrix can be expressed as follows:(57)K0j=∫TjBTjtCg0BTjdxdy;1≤j≤NhK1j=∫TjBTjtCg1BTjdxdy;1≤j≤Nh,where(58)Cg0=Eh1-ν21-νDeβαtDeαβ+Eh1-ν2νDeααtDeλλ+Eh3121-ν21-νDkβαtDkαβ+Eh3121-ν2νDkααtDkλλCg1=Eh1-ν21-νDeβαtDeαβ+Eh1-ν2νDeααtDeλλ+Eh3121-ν21-νDkβαtDkαβ+Eh3121-ν2νDkααtDkλλ+Eh3121-ν21-νDQβαtDeαβ+Eh3121-ν2νDQααtDkλλ+Eh3121-ν21-νDeβαtDQαβ+Eh3121-ν2νDeααtDQλλ+Eh5801-ν21-νDQβαtDQαβ+Eh5801-ν2νDQααtDQλλis the generalized 19×19 behavior matrix.LetF0(x,y)=(BTjtCg0[BTj]) and F1(x,y)=(BTjtCg1[BTj]); then we have(59)K0j=∫TjF0x,ydxdy;1≤j≤Nh=13F023x1+16x2+16x3,23y1+16y2+16y3+13F016x1+23x2+16x3,16y1+23y2+16y3+13F016x1+16x2+23x3,16y1+16y2+23y3K1j=∫TjF1x,ydxdy;1≤j≤Nh=13F123x1+16x2+16x3,23y1+16y2+16y3+13F116x1+23x2+16x3,16y1+23y2+16y3+13F116x1+16x2+23x3,16y1+16y2+23y3.Using (18) the second member of the variational equation is written as follows:(60)Lv˙=∑j=1N∫Tjdv¯TjtBTjtPJTjdS+∮ϖ1∗ϖ2∗dv¯TjtBTjtGJ∗Tjdζ=∑j=1nhdv¯Tjfj,where fj=[F1j1,F1j2,F1j3,F2j1,F2j2,F2j3,F3j1,F3j2,F3j3,F3j4, F3j5,F3j6]T and dv-Tj is defined as in (56); here J∗ is the determinant of the Jacobian which toggles between the linear element dζ of curve ϕ(γ1) to real element dγ. These formulas take into account the load spread over the entire surface and average load at the edges.


## 3.1. Finite Element Space

The continuous truncated variational equation (18) is defined in the space Uad. In order to define a finite element space, we begin by recalling the following results. Let Tj be a triangle with nodes a1=(x1,y1), a2=(x2,y2), and a3=(x3,y3). Let P1 and P2 be the sets of the first- and second-order polynomials with basis 1,x,y and 1,x,y,x2,y2,xy, respectively, then λ1,λ2,λ3 and λ1,λ2,λ3,4λ1λ2,4λ1λ3,4λ2λ3 also generate P1 and P2 where the barycentric basis functions of the triangle are defined by(46)λ1x,y=12∆y3-y2x2-x-x3-x2y2-yλ2x,y=12∆y1-y3x3-x-x1-x3y3-yλ3x,y=12∆y2-y1x1-x-x2-x1y1-y2∆=2∫Tjdxdy=2det⁡111x1x2x3y1y2y3=x3-x2y1-y2-x1-x2y3-y2.LetTh be a triangularization of the midsurface of the shell and nh the number of triangles. We denote by(47)Xh1=vh∈C0S¯,vh∣Tj∈P1Tj,∀Tj,j=1,…,nhthe subspace of dimension Nh (which is the number of nodes) and(48)Xh2=vh∈C0S¯,vh∣Tj∈P2Tj,∀Tj,j=1,…,nh,∇vh=Ghvh=∇P2;S¯ is the closure of S, the subspace of dimension Nh2 which is the number of unisolvent points (nodes and edge mid points).The spaceUad=HΓ01S2×HΓ02(S) is approximated with the finite element space of dimension 3N=2Nh+Nh2(49)℧h=Xh1×Xh1×Xh2.


## 3.2. The Discrete Scheme

LetTh be a triangularization of the shell’s midsurface S and Tj∈Th be a triangle of vertices a1,a2,a3 and midedges a4,a5,a6. Let ηh∈℧h as defined above. Then(50)ηαhx,yTj∈P1,α=1,2;1≤j≤nh,η3hx,yTj∈P2,1≤j≤nh.By using the barycentric polynomials and writingλ4=λ1, 1≤l≤2; 1≤α, β≤2, and 1≤j≤nh,(51)ηαhx,yTj=∑k=13ηαhakλkx,y=∑k=13η-αkλkx,yη3hx,yTj=∑k=13η3hakλkx,y+∑k=46η3hakλk-3x,yλk-2x,y=∑k=13η-3kλkx,y+∑k=46η-3kλk-3x,yλk-2x,y∂αηlhx,yTj=∑k=13∂αη3hakλkx,y,∂αηlhak=Ghηlhak∂αβ2η3hx,yTj=∑k=13∂αβ2η3hakλkx,y,∂αβ2η3hak=Gh∂αη3hak.In the above formula,Gh is the gradient operator and η-34, η-35, and η-36 are, respectively, the unknown values of η3 in the midedges a1a2, a1a3, and a2a3. Let(i)
PνIαβ=PνXαβ and P2Iααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is internal of the mesh,(ii)
PνCαβ=PνXαβ and P2Cααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is in the boundary at the corner of the mesh,(iii)
PνMαβ=PνXαβ and P2Mααβ=P2Xααβ the coefficient of the approximation of the gradient when the vertex is in the boundary not at the corner of the mesh.Exploiting the above, let us take linear element on uniform triangular mesh of a regular pattern; we hereby investigate the case of the element where all three nodes are internal of the mesh; the matrices deduced from the discretization of the gradient are given as follows:(52)A1=λ1x1,y10000000λ2x1+h,y100000000000λ3x1,y1-htA2=0000λ1x2-h,y200λ2x2,y200000000000λ3x2-h,y2-h0tA3=000λ1x3,y3+h00000λ2x3+h,y3+h0000λ3x3,y3000000t.When we return to energy function withm=0,1, we have(53)Amu¯,v¯=∫SEu¯TCgmEv¯dS=∑j=1nh∑k=13∫TjEu¯kTjtCgmEv¯kTjdxdy.Using (43), (45), and (51) in which all three nodes of the triangular element are internal,  let  i=1,2,3; α,β=1,2 and we consider (54)αi=λix,yρk=4λk-3x,yλk-2x,y,k=4,5,6βiα=P1IiαA1:,i+P1IiαA2:,i+P1IiαA3:,iγiα=P2I3αA1:,i+P2I3αA2:,i+P2I3αA3:,iτiαβ=P2I3αβA1:,i+P2I3αβA2:,i+P2I3αβA3:,i.Then the deformation vector can be written asEuTj=BTjduTj where(55)BTjt=α1α2α3000000000β11β21β31000000000β12β12β12000000000000α1α2α3000000000β11β21β31000000000β12β22β32000000000000α1α2α3ρ4ρ5ρ6-2b11α1-2b11α2-2b11α3-2b12α1-2b12α2-2b12α3-γ11-γ21-γ31000-2b11β11-2b11β21-2b11β31-2b12β11-2b12β21-2b12β31-τ111-τ211-τ311000-2b11β12-2b11β22-2b11β32-2b12β12-2b12β22-2b12β32-τ112-τ212-τ312000-2b21α1-2b21α2-2b21α3-2b22α1-2b22α2-2b22α3-γ12-γ22-γ32000-2b21β11-2b21β21-2b21β31-2b22β11-2b22β21-2b22β31-τ112-τ212-τ312000-2b21β12-2b21β22-2b21β32-2b22β12-2b22β22-2b22β32-τ122-τ222-τ322000b¯1α1b¯1α2b¯1α3b¯¯1α1b¯¯1α2b¯¯1α3b11γ11+b12γ12b11γ21+b12γ22b11γ31+b12γ32000b¯1β11b¯1β21b¯1β21b¯¯1β11b¯¯1β21b¯¯1β31b11τ111+b12τ112b11τ211+b12τ212b11τ311+b12τ312000b¯1β12b¯1β22b¯1β22b¯¯1β12b¯¯1β22b¯¯1β32b11τ112+b12τ122b11τ212+b12τ222b11τ312+b12τ322000b¯2α1b¯2α2b¯2α3b¯¯2α1b¯¯2α2b¯¯2α3b21γ11+b22γ12b21γ21+b22γ22b21γ31+b22γ32000b¯2β11b¯2β21b¯2β31b¯¯2β11b¯¯2β21b¯¯2β31b21τ111+b22τ112b21τ211+b22τ212b21τ311+b22τ312000b¯2β12b¯2β22b¯2β32b¯¯2β12b¯¯2β22b¯¯2β32b21τ112+b22τ122b21τ212+b22τ222b21τ312+b22τ322000,(56)du¯Tj=u¯1j1,u¯1j2,u¯1j3,u¯2j1,u¯2j2,u¯2j3,u¯3j1,u¯3j2,u¯3j3,u¯3j4,u¯3j5,u¯3j6t.


## 3.3. Stiffness Matrix

Using local element, stiffness matrix can be expressed as follows:(57)K0j=∫TjBTjtCg0BTjdxdy;1≤j≤NhK1j=∫TjBTjtCg1BTjdxdy;1≤j≤Nh,where(58)Cg0=Eh1-ν21-νDeβαtDeαβ+Eh1-ν2νDeααtDeλλ+Eh3121-ν21-νDkβαtDkαβ+Eh3121-ν2νDkααtDkλλCg1=Eh1-ν21-νDeβαtDeαβ+Eh1-ν2νDeααtDeλλ+Eh3121-ν21-νDkβαtDkαβ+Eh3121-ν2νDkααtDkλλ+Eh3121-ν21-νDQβαtDeαβ+Eh3121-ν2νDQααtDkλλ+Eh3121-ν21-νDeβαtDQαβ+Eh3121-ν2νDeααtDQλλ+Eh5801-ν21-νDQβαtDQαβ+Eh5801-ν2νDQααtDQλλis the generalized 19×19 behavior matrix.LetF0(x,y)=(BTjtCg0[BTj]) and F1(x,y)=(BTjtCg1[BTj]); then we have(59)K0j=∫TjF0x,ydxdy;1≤j≤Nh=13F023x1+16x2+16x3,23y1+16y2+16y3+13F016x1+23x2+16x3,16y1+23y2+16y3+13F016x1+16x2+23x3,16y1+16y2+23y3K1j=∫TjF1x,ydxdy;1≤j≤Nh=13F123x1+16x2+16x3,23y1+16y2+16y3+13F116x1+23x2+16x3,16y1+23y2+16y3+13F116x1+16x2+23x3,16y1+16y2+23y3.Using (18) the second member of the variational equation is written as follows:(60)Lv˙=∑j=1N∫Tjdv¯TjtBTjtPJTjdS+∮ϖ1∗ϖ2∗dv¯TjtBTjtGJ∗Tjdζ=∑j=1nhdv¯Tjfj,where fj=[F1j1,F1j2,F1j3,F2j1,F2j2,F2j3,F3j1,F3j2,F3j3,F3j4, F3j5,F3j6]T and dv-Tj is defined as in (56); here J∗ is the determinant of the Jacobian which toggles between the linear element dζ of curve ϕ(γ1) to real element dγ. These formulas take into account the load spread over the entire surface and average load at the edges.


## 4. Validation Tests of Our Finite Element

The aim of this study is to investigate the accuracy of N-T thick shells theory for linear elastic shell by using Spherical Shell Finite Element (SSFE). We lay our investigation on a well-known benchmark as hemisphere under diametrically opposite charges given in Figure1 used to evaluate the performance of a shell element. The computed deformed limit surface of a quarter of the hemisphere is shown in Figure 2 plotted with the Matlab R2015a tools and the displacement convergence results are shown in Figure 3 and Table 1. We monitor the displacement in load point A presented in Figure 1.Table 1
Table of results of displacement at pointA for the hemisphere.
Thickness = 0.04 and force = 0.5 Nsol_ref = reference solution (UA∗103)
N
DKT12
DKT18
SKLFE
SFE3 (cmc)
SSFE
sol_ref
2
121
87
—
4.1
—
94
3
—
—
1.6
—
1.61
94
4
108
94
10.2
28
10.21
94
6
102
94
52.3
63
52.3
94
8
99
93
73.8
77
73.8
94
10
98
93
81.56
82
81.55
94
12
96
93
86.3
84
86.32
94
15
88.4
86
88.4
94
20
86.5
87
86.52
94Figure 1
Hemisphere shell benchmark.Figure 2
Deformed configuration of one quarter of the hemisphere.Figure 3
Convergence in the hemisphere.Here, the case study is that of a thin shell of hemisphere subjected to four opposite diametrically concentrated loads at the base proposed by Macneal and Harder [20] which is a standard test. This benchmark is usually used to verify the absence of membrane locking and good representation of rigid bodies motion. The hemisphere undergoes large rotations around the normal of the middle surface. Deformations of inextensible bending membrane are also important and this problem is therefore an excellent test to examine the ability of a shell element to represent the rigid and inextensible modes. The geometrical and mechanical characteristics are indicated for h/R=0.004, the radius R=10m, the thickness h=0.04m, the angle subtended by the north pole of the hemisphere is θ=180, Young’s modulus is E=6.825e+7Pa, Poisson’s ratio is ν=0.3, and the diametrically opposite loads at points A and B are 0.5N. The limits conditions W=0 in the slib bound E and the symmetric conditions are given by V=0 on the edge bound AC and U=0 on the edge bound BD.


### 4.1. Convergence

A reference solution presented in [20] provides for displacement in the direction of the loads as follows: UA=VB=0.094m. Only 1/4 of the hemisphere is discretized because of the symmetry of loads and the geometry. Both Kirchhoff-Love shell theory and N-T shell theory are computed using SSFE model and the respective results are analysed, compared with DKT12 and DKT18 proposed in [21] and SFE3 [22] then commented. Table 1 shows the displacement results at point A and Figure 3 perfectly describes their variation and the rate of convergence. The convergence properties of the method are clearly shown from Figure 3 and Table 1. Then SSFE converge as well as both the semifinite element (SFE) and Discrete Kirchhoff Triangle (DKT) elements for the membrane displacement at load point A.


### 4.2. Scaling and Deviation

Scaling of the ratio2χ=h/R given in Table 2 is proceeded on the range of the following values 0.006, 0.099, 0.12, 0.15, 0.175, and 0.2 of the spherical shell. Notice that the radius R is constant while the thickness varies with the ratio. We observe in Tables 1 and 2 and Figure 3 that, for the thickness ratio 0≺2χ≺0.099, the membrane displacement in load point A is the same for both K-L and N-T models. When the ratio 2χ=h/R is greater or equal to 0.099, the displacement computed for inextensible bending membrane from spherical equation of K-L and N-T is not the same. This means that, above 2χ≃1/10, both K-L and N-T approaches are different for all values of the thickness ratio.Table 2
Deviation of both Reissner-Mindlin (R-M) and N-T thick shell according to the scaling of the ratioh/R.
Ratioh/R
Theory
Mesh
5
×
5
7
×
7
9
×
9
11
×
11
13
×
13
15
×
15
0.006
R-M
1.77
∗
1
0
-
5
3.91
∗
1
0
-
4
0.0035
0.088
0.0954
0.0950
N-T
1.77
∗
1
0
-
5
3.91
∗
1
0
-
4
0.0035
0.0877
0.0954
0.0950
0.099
R-M
1.14
∗
1
0
-
5
2.64
∗
1
0
-
4
0.0018
0.0109
0.0209
0.0204
N-T
1.14
∗
1
0
-
5
2.64
∗
1
0
-
4
0.0018
0.0110
0.0209
0.0204
0.12
R-M
2.91
∗
1
0
-
6
8.76
∗
1
0
-
6
4.52
∗
1
0
-
3
1.07
∗
1
0
-
2
2.76
∗
1
0
-
2
2.47
∗
1
0
-
2
N-T
4.67
∗
1
0
-
6
9.72
∗
1
0
-
6
3.52
∗
1
0
-
3
9.90
∗
1
0
-
3
2.67
∗
1
0
-
2
2.47
∗
1
0
-
2
0.15
R-M
1.07
∗
1
0
-
6
1.20
∗
1
0
-
5
4.67
∗
1
0
-
4
5.02
∗
1
0
-
3
1.57
∗
1
0
-
2
1.39
∗
1
0
-
2
N-T
3.03
∗
1
0
-
6
4.21
∗
1
0
-
6
1.79
∗
1
0
-
4
2.35
∗
1
0
-
3
1.36
∗
1
0
-
2
1.17
∗
1
0
-
2
0.175
R-M
1.23
∗
1
0
-
6
2.87
∗
1
0
-
6
3.20
∗
1
0
-
4
3.22
∗
1
0
-
3
9.41
∗
1
0
-
3
8.96
∗
1
0
-
3
N-T
8.20
∗
1
0
-
7
1.02
∗
1
0
-
6
2.30
∗
1
0
-
4
3.41
∗
1
0
-
4
7.12
∗
1
0
-
3
6.89
∗
1
0
-
3
0.2
R-M
5.36
∗
1
0
-
7
3.51
∗
1
0
-
6
4.89
∗
1
0
-
4
2.89
∗
1
0
-
3
5.88
∗
1
0
-
3
5.52
∗
1
0
-
3
N-T
4.01
∗
1
0
-
7
3.10
∗
1
0
-
6
3.71
∗
1
0
-
4
2.49
∗
1
0
-
4
5.57
∗
1
0
-
3
5.23
∗
1
0
-
3We investigate now the deviation between K-L and N-T displacements in load pointA. With the variation of thickness ratio 2χ=0.06,0.099,0.12,0.15, 0.175, and 0.2, the results plotted in Table 2 and Figures 4–7 clearly show that the deviation of displacement is encountered at the specified values of 2χ above. This deviation increases with the number of meshes at the load point A. We also observe that the deviation increases with thickness ratio.Figure 4
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.006.Figure 5
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.099.Figure 6
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.12.Figure 7
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.15.


### 4.3. Impact of the Transverse Stresses

Surveys of various shear stresses in linear elastic thick shells can be found in the works of [3] where they are mathematically substantiated. We have obtained good convergence results for thin shells with the half side element N=12 and a pressure in the load point PA=0.53Pa for a spherical shell. For the thickness coordinate z∈-h/2,+h/2,whereh=1.2m, the numerical results of the transverse shear stresses through the thickness are computed (see Table 3).Table 3
Transverse shear stresses of a thick spherical shell under load pointA convergence analysis.
z
−0.6
−0.466
−0.2
0.066
0.2
0.467
0.6
σ
13
(
P
a
)
0
0.005
0.001
0.019
0.021
0.025
0.259
σ
23
(
P
a
)
0
0.007
0.018
0.026
0.028
0.034
0.035
σ
33
(
P
a
)
1.6
×
1
0
-
14
−0.132
0.309
−0.417
−0.455
−0.514
−0.533After shear stresses investigation, we observe that they satisfy tangential stress-free boundary conditions at bottom to top surfaces of the panel. These results show predicting ability of SSFE finite element based on N-T’s model. It also appears that the shear stresses vary through the thickness as compared to 3D thick shell equations. We can also appreciate that transverse shear stressesσ13 and σ23 are not negligible in thick shells.


### 4.4. Discussion

C
0 finite element discretization of the first-order N-T model applied the PPR method to solve problems of gradient discontinuities across edges of triangular elements during interpolation. The displacement is two-degree polynomial of the thickness ratio χ while the membrane strain tensor contains the change in the third fundamental form.The convergence of SSFE has been clearly established and this element uses 12 degrees of freedom per triangular element which is robust and less greedy (in terms of memory) than DKT12, DKT18, and SFE that are based on the Kirchhoff-Love (thin shells) and Reissner-Mindlin (thick shells) approaches which neglect the third fundamental form in their shell kinematic equations [23]. The constitutive law through the strain tensor contains the change of third fundamental form Qαβ as shown in the tables; the thickness ratio χ impacts the behavior of the shell because of the significant energy contribution due to Qαβ as the ratio increases unlike the K-L classical model [3].Recall that, in (17), we see that total deformation energy due to K-L and R-M models contains membrane deformation energy Eee and bending deformation energy Ekk; that is, ER-M=Eee+Ekk. Equation (18) shows that the total deformation energy due to N-T model contains additional terms: coupled membrane and Gaussian bending (EeQ,EQe) and Gaussian deformation energy EQQ. Then EN-T=ER-M+EQQ+EeQ+EQe.WhereEee=Eh/Rαee, Ekk=E(h/R)h/R2αkk, EQQ=10-1E(h/R)h/R4αQQ, and EeQ=EQe=10-1E(h/R)h/R2αeQ, which represents the portion of energy contribution, αee,αkk, and αQQ are constants which do not depend on h/R.As we mentioned above, when the thickness ratio2χ is greater than 0.098, this additional energy (EeQ,EQe) influences global deformation energy and shows the difference between N-T and R-M models applied to the spherical thick shells.The investigation of the variation of the thickness ratioh/R for certain values 0≺2χ≺0.099 proves that the Kirchhoff-Love, Reissner-Mindlin, and N-T classical models have the same contribution of total deformation energy. The additional terms containing the change of the third fundamental form do not have here any influence. The energies EeQ and EQQ disappear in the global deformation energy when 2χ≺0.1 or χ2≺0.01 because they are inversely proportional, respectively, to 104 and 108. In this case EN-T=ER-M. But if the thickness ratio is 0.099≤h/R≺1, displacement results encountered for both models are different because EeQ and EQQ do not disappear in the global deformation energy. The consequence is that they enhance global deformation energy and improve the rigidity of the shell structure. Moreover, N-T’s model brings us real facilitation to determine the distribution of transverse shear stresses through the thickness.


## 4.1. Convergence

A reference solution presented in [20] provides for displacement in the direction of the loads as follows: UA=VB=0.094m. Only 1/4 of the hemisphere is discretized because of the symmetry of loads and the geometry. Both Kirchhoff-Love shell theory and N-T shell theory are computed using SSFE model and the respective results are analysed, compared with DKT12 and DKT18 proposed in [21] and SFE3 [22] then commented. Table 1 shows the displacement results at point A and Figure 3 perfectly describes their variation and the rate of convergence. The convergence properties of the method are clearly shown from Figure 3 and Table 1. Then SSFE converge as well as both the semifinite element (SFE) and Discrete Kirchhoff Triangle (DKT) elements for the membrane displacement at load point A.


## 4.2. Scaling and Deviation

Scaling of the ratio2χ=h/R given in Table 2 is proceeded on the range of the following values 0.006, 0.099, 0.12, 0.15, 0.175, and 0.2 of the spherical shell. Notice that the radius R is constant while the thickness varies with the ratio. We observe in Tables 1 and 2 and Figure 3 that, for the thickness ratio 0≺2χ≺0.099, the membrane displacement in load point A is the same for both K-L and N-T models. When the ratio 2χ=h/R is greater or equal to 0.099, the displacement computed for inextensible bending membrane from spherical equation of K-L and N-T is not the same. This means that, above 2χ≃1/10, both K-L and N-T approaches are different for all values of the thickness ratio.Table 2
Deviation of both Reissner-Mindlin (R-M) and N-T thick shell according to the scaling of the ratioh/R.
Ratioh/R
Theory
Mesh
5
×
5
7
×
7
9
×
9
11
×
11
13
×
13
15
×
15
0.006
R-M
1.77
∗
1
0
-
5
3.91
∗
1
0
-
4
0.0035
0.088
0.0954
0.0950
N-T
1.77
∗
1
0
-
5
3.91
∗
1
0
-
4
0.0035
0.0877
0.0954
0.0950
0.099
R-M
1.14
∗
1
0
-
5
2.64
∗
1
0
-
4
0.0018
0.0109
0.0209
0.0204
N-T
1.14
∗
1
0
-
5
2.64
∗
1
0
-
4
0.0018
0.0110
0.0209
0.0204
0.12
R-M
2.91
∗
1
0
-
6
8.76
∗
1
0
-
6
4.52
∗
1
0
-
3
1.07
∗
1
0
-
2
2.76
∗
1
0
-
2
2.47
∗
1
0
-
2
N-T
4.67
∗
1
0
-
6
9.72
∗
1
0
-
6
3.52
∗
1
0
-
3
9.90
∗
1
0
-
3
2.67
∗
1
0
-
2
2.47
∗
1
0
-
2
0.15
R-M
1.07
∗
1
0
-
6
1.20
∗
1
0
-
5
4.67
∗
1
0
-
4
5.02
∗
1
0
-
3
1.57
∗
1
0
-
2
1.39
∗
1
0
-
2
N-T
3.03
∗
1
0
-
6
4.21
∗
1
0
-
6
1.79
∗
1
0
-
4
2.35
∗
1
0
-
3
1.36
∗
1
0
-
2
1.17
∗
1
0
-
2
0.175
R-M
1.23
∗
1
0
-
6
2.87
∗
1
0
-
6
3.20
∗
1
0
-
4
3.22
∗
1
0
-
3
9.41
∗
1
0
-
3
8.96
∗
1
0
-
3
N-T
8.20
∗
1
0
-
7
1.02
∗
1
0
-
6
2.30
∗
1
0
-
4
3.41
∗
1
0
-
4
7.12
∗
1
0
-
3
6.89
∗
1
0
-
3
0.2
R-M
5.36
∗
1
0
-
7
3.51
∗
1
0
-
6
4.89
∗
1
0
-
4
2.89
∗
1
0
-
3
5.88
∗
1
0
-
3
5.52
∗
1
0
-
3
N-T
4.01
∗
1
0
-
7
3.10
∗
1
0
-
6
3.71
∗
1
0
-
4
2.49
∗
1
0
-
4
5.57
∗
1
0
-
3
5.23
∗
1
0
-
3We investigate now the deviation between K-L and N-T displacements in load pointA. With the variation of thickness ratio 2χ=0.06,0.099,0.12,0.15, 0.175, and 0.2, the results plotted in Table 2 and Figures 4–7 clearly show that the deviation of displacement is encountered at the specified values of 2χ above. This deviation increases with the number of meshes at the load point A. We also observe that the deviation increases with thickness ratio.Figure 4
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.006.Figure 5
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.099.Figure 6
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.12.Figure 7
Variation of membrane displacementUA at A for the thickness ratio 2χ=h/R=0.15.


## 4.3. Impact of the Transverse Stresses

Surveys of various shear stresses in linear elastic thick shells can be found in the works of [3] where they are mathematically substantiated. We have obtained good convergence results for thin shells with the half side element N=12 and a pressure in the load point PA=0.53Pa for a spherical shell. For the thickness coordinate z∈-h/2,+h/2,whereh=1.2m, the numerical results of the transverse shear stresses through the thickness are computed (see Table 3).Table 3
Transverse shear stresses of a thick spherical shell under load pointA convergence analysis.
z
−0.6
−0.466
−0.2
0.066
0.2
0.467
0.6
σ
13
(
P
a
)
0
0.005
0.001
0.019
0.021
0.025
0.259
σ
23
(
P
a
)
0
0.007
0.018
0.026
0.028
0.034
0.035
σ
33
(
P
a
)
1.6
×
1
0
-
14
−0.132
0.309
−0.417
−0.455
−0.514
−0.533After shear stresses investigation, we observe that they satisfy tangential stress-free boundary conditions at bottom to top surfaces of the panel. These results show predicting ability of SSFE finite element based on N-T’s model. It also appears that the shear stresses vary through the thickness as compared to 3D thick shell equations. We can also appreciate that transverse shear stressesσ13 and σ23 are not negligible in thick shells.


## 4.4. Discussion

C
0 finite element discretization of the first-order N-T model applied the PPR method to solve problems of gradient discontinuities across edges of triangular elements during interpolation. The displacement is two-degree polynomial of the thickness ratio χ while the membrane strain tensor contains the change in the third fundamental form.The convergence of SSFE has been clearly established and this element uses 12 degrees of freedom per triangular element which is robust and less greedy (in terms of memory) than DKT12, DKT18, and SFE that are based on the Kirchhoff-Love (thin shells) and Reissner-Mindlin (thick shells) approaches which neglect the third fundamental form in their shell kinematic equations [23]. The constitutive law through the strain tensor contains the change of third fundamental form Qαβ as shown in the tables; the thickness ratio χ impacts the behavior of the shell because of the significant energy contribution due to Qαβ as the ratio increases unlike the K-L classical model [3].Recall that, in (17), we see that total deformation energy due to K-L and R-M models contains membrane deformation energy Eee and bending deformation energy Ekk; that is, ER-M=Eee+Ekk. Equation (18) shows that the total deformation energy due to N-T model contains additional terms: coupled membrane and Gaussian bending (EeQ,EQe) and Gaussian deformation energy EQQ. Then EN-T=ER-M+EQQ+EeQ+EQe.WhereEee=Eh/Rαee, Ekk=E(h/R)h/R2αkk, EQQ=10-1E(h/R)h/R4αQQ, and EeQ=EQe=10-1E(h/R)h/R2αeQ, which represents the portion of energy contribution, αee,αkk, and αQQ are constants which do not depend on h/R.As we mentioned above, when the thickness ratio2χ is greater than 0.098, this additional energy (EeQ,EQe) influences global deformation energy and shows the difference between N-T and R-M models applied to the spherical thick shells.The investigation of the variation of the thickness ratioh/R for certain values 0≺2χ≺0.099 proves that the Kirchhoff-Love, Reissner-Mindlin, and N-T classical models have the same contribution of total deformation energy. The additional terms containing the change of the third fundamental form do not have here any influence. The energies EeQ and EQQ disappear in the global deformation energy when 2χ≺0.1 or χ2≺0.01 because they are inversely proportional, respectively, to 104 and 108. In this case EN-T=ER-M. But if the thickness ratio is 0.099≤h/R≺1, displacement results encountered for both models are different because EeQ and EQQ do not disappear in the global deformation energy. The consequence is that they enhance global deformation energy and improve the rigidity of the shell structure. Moreover, N-T’s model brings us real facilitation to determine the distribution of transverse shear stresses through the thickness.


## 5. Conclusion

The finite element SSFE using the PPR method for the spherical shell described in this paper forC0 finite element triangularization has given accurate results and is skilled to design thicker shells. It has proved to be faster and uses less memory than other well-known methods used in the different benchmarks. The N-T model handles spherical thin and thick shell properly because it clearly shows how the change of the third fundamental form enhances the total deformation energy when the ratio χ becomes greater.Transverse stresses which have been predicted by the 2D governing equations of N-T model were calculated numerically and correctly. Structural engineers can therefore design more accurately spherical thick shells.



---

*Source: 5903503-2017-06-14.xml*