# A 2D Model Which Accounts for Transverse Strains in a Linear Elastic Thick Shell

Authors: Nzengwa Robert; Djopkop Kouanang Landry; Bodol Momha Merlin; Amba Jean Chills; Nkongho Anyi Joseph; Zoa Ambassa
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper we present a 2D four-parameter model that accounts for the variation of transverse deformation through-the-thickness in thick linear elastic shells. This model is deduced directly from the 3D elasticity equations of the traction-displacement boundary value problem. Transverse shear strains and thickness variation are accountable through additional terms which appear in the final equations besides well-known terms in the classical Kirchhoff-Love models for thin shells. A unique solution of the variational equation is established, and numerical and analytical results are compared with satisfaction.


---

## Body

## 1. Introduction

Shell structures are increasingly being used in modern complex constructions. An outstanding problem is how to account for the through-the-thickness stretching and shear strains in a two-dimensional model when the shell becomes thick. Letx=x1,x2 be the coordinates of the generic point of the mid-surface of the shell, and z the transverse coordinate. In order to derive a 2D shell model from the 3D solid, a displacement of the form as shown in following equation was proposed in [1].(1)x=φ+∑n=1∞zndn,where dn are the shell directors and φ the position vector of the mid-surface in the deformed configuration. Several attempts to approach the kinematics of a point on the shell are presented in the literature to offer 2D models, although there are also 3D models. These latter have the advantage of avoiding complex shell finite elements [2]. But numerically, a locking phenomenon arises when the thickness of the shell is small, and in addition, numerical approximation requires a voluminous mesh in three dimensions to approach the exact solution. The numerical implementation is memory-greedy and needs an extensive calculating time, hence the necessity to search for efficient 2D models. In this sense, several authors ([3–6]) have proposed 2D models based on the works of Love, Sander and Kirchhoff, where the proposed kinematics satisfy the Kirchhoff-Love hypotheses, i.e., normal segments to the mid-surface before deformation remain rigid and normal to the deformed mid-surface εxz=εyz=εzz=0. So, transverse strains do not contribute to the strain energy. Also, transverse shear stresses σα3 and the pinch stress σ33 are neglected. In [7], the nonlocal model of strain gradient small-scale approach and the first-order shear deformation shell theory (FSDST), initially proposed by Hencky, Timoshenko-Mindlin, Naghdi-Berry, and Reissner-Mindlin ([8–10]), which are 2D models were proposed. These models calculate transverse shear stresses but do not address variations in thickness (see ([11–13]). Some authors restricted the series of equation (1) into a linear or quadratic formula in z ([14–16]). In this case, the number of parameters (or unknown functions) of the model is 6 or 7 unlike in thin shell theory where they are only 3. Moreover, the strain tensor which is a change in metrics, is either linear or quadratic in z and does not depend on the change of the fundamental forms of the mid-surface as expressed in ([17, 18]) for linear elastic thick shells. The N-T theory implemented in ([19–22]) proposed by [17] introduced the Gauss deformation strain tensor (change of the third fundamental form) and calculated the transverse stresses σxz,σyz,σzz as reactions to a more general plane strain. The setting is quite different from that of the classical shell theory. However, it does not calculate transverse strains. In this paper we present a 2D model also derived from the 3D elasticity by introducing a particular form of admissible displacements obtained as a solution of a torsion loading problem in 3D solid. Let ϕ be the vector field representing the total transverse strains. We assume rotϕ=0 as widely used in 3D solids for torsion problems. The solution of the resulting ordinary differential equation is in-plane quadratic and depends transversally on a stretching function qx,z=wzq¯x, wz is a transverse distribution function. The displacement appears to be the sum of the plane strain thick shell displacement field [18] with a warping in-plane tensor and a through-the-thickness stretching term. The number of parameters is 4. The resulting 2D model accounts for transverse strains and stresses and appears to be an improvement on the linear elastic thick shell theory. Moreover, the transverse stretching can be described by a nonpolynomial function in z or a polynomial function of any order according to the desired accuracy. In the work of [23], a good conformity is observed between the experimental and numerical predictions for the choice of a linear polynomial distribution stretching function.The rest of the paper is organized as follows: Section2 is devoted to the study of the admissible displacement and its consequences on strains and constitutive laws. In Section 3 we establish the existence and uniqueness of the solution to a model problem, the classical traction-displacement variational problem which is the 2D model. We next present convergence tests and compare numerical and analytical results in Section 4. We conclude in Section 5 after some discussions.The repeated index convention is adopted unless otherwise specified. The derivative of a function∂u/∂xi will also be denoted u,i or ∂iu. Greek and Latin indexes range in 1,2 and 1,2,3 respectively. Covariant derivation will be denoted “|” or ∇ in the 3D shell and mid-surface respectively. The scalar product of matrices A and B will be denoted A:B. Further notations will be specified in the text.


## 2. Admissible Displacements and Constitutive Law

LetΩ=S×−h,h (S is the mid-surface and h > 0 is half the thickness) denote a shell. We assume the surface S is bounded and sufficiently smooth [24] for all subsequent computations. Let a1,a2,a3 and a1,a2,a3 denote the covariant and contravariant bases of the mid-surface and G1,G2,G3,G1,G2,G3 the covariant and contravariant bases of the shell. Then(2)x=x1,x2Gα=δαγ−zbαγaγ=μαγaγ,G3=a3Gα=δγα−zbγα−1aγ=μ−1γαaγ,G3=a3,dΩ=G1,G2,a3dx1dx2dz=1−2zH¯+z2K¯dzdS=ψx,zdzdS;2H¯=bαα,K¯=detb,where bαρ=aργbγα and bγα denote the curvature tensor components and aαβ is the contra variant component of the metric of the mid-surface S. We have implicitly assumed that the characteristic parameter of the mid-surface χ=h/R is less than 1; h is half the thickness and R the minimum absolute value of its radius of curvature. A vector field v:Ω¯⟶ℝ3 can be expressed component wise indifferently in the G-base or the a-base as follows:(3)v=vix,zGi=v¯αx,zaα+v3x,za3and deduce thatvα=μαγv¯γ,v¯α=μ−1αγvγ.Then the strain tensor [17] reads(4)∈αβv=vα|β+vβ|α2=μαν∇βv¯ν−bαβv¯3+μβν∇αv¯ν−bνβv¯32,∈α3v=vα|3+v3|α2=μανv¯ν,3+v¯3,α+bανv¯ν2,∈33v=v3,3=dv3dz.We deduce the form of admissible displacement and strain from the followingLemma 1.
Letϕx1,x2,z=ϕi be such that ϕαx,z=∈α3+∈3α=2∈α3,ϕ3x,z=∈33. We suppose rotϕ=0; rotϕi=ϕi+1,i+2−ϕi+2,i+1 mod [3]
Then there exist functionsξixandqx,z such that(5)vα=ξα−z∂αξ3+2bαρξρ+z2bαρbρτξτ+bαρ∂ρξ3,v3=ξ3+qx,z,∈αβv=eαβξ−zkαβξ+z2Qαβξ+qx,zϒαβx,z,∈α3v=∂αq2,∈33v=∂zq,eαβξ=∇​αξβ+∇βξα−2bαβξ32,kαβξ=∇αbβυξυ+bαυ∇βξυ+bβυ∇αξυ+∇α∇βξ3−bατbτβξ3,Qαβξ=bαυ∇βbυτξτ+∇υξ3+bβυ∇αbυτξτ+∇υξ32,ϒαβx,z=−μανbνβ+μβνbνα2,where e, k, and Q respectively change in the first, second, and third fundamental forms of the mid-surface and ϒ is the section warping tensor, bβα the mixed components of the curvature tensor, bρα the covariant components of the curvature tensor, and δαρ Kronecker’s symbol.Proof.
We shall denotev3,1,v3,2 by ∇v3. Let(6)Aβα=μ−1βα,A−1βα=μβα.
The tensorsA and b commute and we deduce from AA−1=I that dA/dz=A2b. From equation (2) we deduce that there exists qx,z such that ϕi=∂iq which is equivalent to ∈α3=∂αq/2,∈33=∂zq or(7)A−1dv¯dz+bv¯+∇v3=∇q,dv3dz=∂zq.
Equation (7) are also equivalent to(8)dAv¯dz+A2∇v3−q=0,dv3dz=∂zq.
Letw=∇v3−q then equation (8) becomes (9)dAv¯dz+A2w=0,dv3dz=∂zq.
A solution to equation (9) is
v3=η3+q,v¯=A−1η¯−z∇η3,η¯=η1,η2. In fact since dw/dz=0, one has(10)dAA−1η¯−zwdz+A2w=dη¯−zAwdz+A2w=−zA2bw−Aw+A2w=−A2zb+A−1w+A2w=−A2w+A2w=0.=0.
Therefore,(11)v¯=A−1η¯−z∇η3,v3=η3+q,
or(12)vα=ηα−z∂αη3+2bαρηρ+z2bαρbρτητ+bαρ∂ρη3,v3=η3+qx,z.
The solutionv can be decomposed as a sum of two displacements namely a plane strain displacement field vη=vαη,η3 [17] and a stretching displacement vq=0,0,q. Therefore the strain can also be decomposed as(13)∈ijv=∈ijvη+∈ijvq=∈η+∈q,where(14)∈αβvη=eαβη−zkαβη+z2Qαβη=∈αβη,∈i3vη=0=∈i3η,∈αβvq=qx,zϒαβx,z=∈αβq,∈α3vq=∂αq2=∈α3q,∈33vq=∂zq=∈33q,and the proof is complete by taking η=ξ.
It should be noted that ifϕ=0, equation (2) will be satisfied. Therefore, plane strains hypothesis which leads to the N-T displacement is another form of equation (2). The in-plane displacement is quadratic in z with an additional warping strain. The stretching function can be expressed as follows:(15)qx,z=wzq¯x,where w is an arbitrary nonconstant function or a polynomial of any degree. The number of unknown functions is 4.
Let us consider the strain decomposition∈=∈η+∈q=∈η+∈q. For a linear elastic isotropic homogeneous material with Lamé constants λ¯,μ¯, the stress-strain constitutive relations is as follows:(16)σ=λ¯∈llG+2μ¯∈,=σ∈η+σ∈q.
Since∈33η=0, we have(17)σ33=λ¯∈ααη+∈ρρq+λ¯+2μ¯∈33q.
Ifσ33=0, then(18)∈33=−λ¯λ¯+2μ¯∈ααη+∈ααq,σα3=2μ¯∈α3q,and plane stress can be reproduced only if ∂αq=0.


## 3. Presentation of the 2D Model

Let the border ofS be partitioned by ∂S=γ0∪γ such that the border of the shell reads ∂Ω, =​Γ0​∪, Γ1∪, Γ−∪Γ+;Γ0=γ0×−h,h;Γ1=γ1×−h,h;Γ−=S×−h;Γ+=S×h. The model problem is the standard traction-displacement boundary value problem defined by: find the displacement field u:x,z∈Ω¯⟶ℝ3, (Ω¯ is the closure of Ω) which solves the following equations:(19)−divσ=finΩ,u.=0onΓ0,σn⟶=ponΓ=Γ1∪Γ−∪Γ+,σij=λ¯∈lluGij+2μ¯∈iju,where n⟶ is the unit outer normal vector; fε, p volume and surface force densities, respectively; divσ=σjij. The variational equation associated to the model problem reads: find u∈IHΓ01Ω=v∈H1Ω3,v=0onΓ0 such that(20)Eu,v=Lv,for allv∈IHΓ01Ω,Eu,v=∫Ωσiju∈ijvdΩ,Lv,=∫ΩfividΩ+∫ΓpividΓ,HmD for a natural number m≠0 and a domain D in ℝ1,ℝ2orℝ3 designates a Sobolev space; H0D=L2D. Volume and surface forces are supposed to be sufficiently regular.The model problem has a unique solution because the operatorE.,. is coercive in the space of admissible displacements IHΓ01Ω. We shall be interested in the particular case where qx,z=wzq¯x. Let Uad=v∈IHΓ01Ω,Rotϕ=0, then Uad is a closed subset of IHΓ01Ω. Let v be in Uad, then because w is nonconstant in z, we deduce from (7) and the boundary condition on the displacement that η1=η2=η3=q¯=∂αη3=0onγ0. It is easy to check that(21)ηα∈Hγ01S,q¯∈Hγ01S,η3∈Hγ02S,Hγ02S=η3∈H2S,η3=∂ν⟶η3=0onγ0,ν⟶ is the outer unit vector on the border of the mid-surface. We denote the set of admissible displacements and the bilinear form, respectively, by(22)Uadw=η,q¯,ηα∈Hγ01S,q¯∈Hγ01S,η3∈Hγ02S,Eu,v=Ewξ,q¯,η,y¯,where u=uξ,q¯,v=vη,y¯ are constructed as in equation (12).Lemma 2.
The operatorEw is coercive in Uadw.Proof.
LetX=Xη,q¯∈Uadwtheneαβη,kαβη,Qαβη,q¯,∂1q¯,∂2q¯∈L2S. We next define(23)q¯L2S2=max∫Ωw2zq¯2xdΩ,∫Ωw'2zq¯2xdΩ,∂αq¯L2S2=∫Ωw2z∂αq¯2xdΩ,X¯=e,k,Q,q¯,∂1q¯,∂2q¯,Θ=X¯∈L2S3×L2S3×L2S3×L2S×L2S×L2S,equipped with the norm(24)X¯=e2+k2+Q2+q¯2+∂1q¯2+∂2q¯21/2,where we have voluntarily omitted the subscript L2(S) on the norms. We deduce from the expression of the strain tensor that the bilinear form Ew can also be considered as a quadratic form on the space Θ where it is also continuous for the weak topology. Therefore, it is bounded in the unit ball and there exist positive constants m, M such that m≤EwX¯,X¯≤M. We deduce that for all X∈Θ, one has mX¯2≤EwX¯,X¯≤MX¯. More over there exists X¯min the unit ball such thatEwX¯m,X¯m=m. We infer that m≻0. In fact suppose m=0. Then(25)0=EwX¯m,X¯m≥2μ∫Ωem−zkm+z2Qm+wzq¯mϒ:em−zkm+z2Qm+wzq¯mϒ+q¯m2+∂1q¯m2+∂2q¯m2.
We deduce thatq¯m=0 and from [17] we have em=km=Qm=0 which is in contradiction with X¯m=1. Therefore EwX,X=EwX¯,X¯≥mX¯2. Again from [17, 25] we deduce that there exists a strictly positive constant C such that EwX,X≥CX2 and the proof is complete.Theorem 1.
The 2D variational equation for the genuine traction-displacement problem has a unique solutionξxq¯x which satisfies the equations: find ξx,q¯x such that ξ.,q¯.∈Uadw;(26)Ewu,v=Lwv,v∈Uadw,uα=ξαx−z∂αξ3x+2bαρξρx+z2bαρbρτξτx+bαρ∂ρξ3x,u3=ξ3x+wzq¯x,vα=ηα−z∂αη3+2bαρηρ+z2bαρbρτητ+bαρ∂ρη3,v3=η3x+wzy¯x,v3=η3x+wzy¯x,Ewu,v=∫S∫−hhσiju∈ijvψx,zdzdS,Lwv=∫S∫−hhfix,z,tviψx,zdzdS+∫Γpix,z,tvidΓ;=∫SFixηi+F4y¯dS+∫γ1Ηixηi+Η4y¯dγ1+∫γ1mαθαdγ1.
Here the stressσ and strain ∈ fields are defined through equations (8)–(13) Fi,F4,Hi,H4,mα,θα are, respectively, resultant surface force, shear force, moment density, and the opposite angle of rotation of a section at the border.Proof.
Resultant forces are also sufficiently smooth and the bilinear form is coercive. Existence and uniqueness of a solution are obtained through the Lax-Milgram lemma.


## 4. Validation of the Model

In the numerical validation section we will use a finite element that we name “MT6” or “T6-m”, which is a modified T6 finite element. Figure1 below presents the shape of both elements.Figure 1
(a) Triangular mesh. (b) MT6 Finite Element.Let us consider a shellΩ=S×−h/2,h/2 whose constitutive law is(27)σiju,q=λ¯∈ααu+∈llqGij+2μ¯∈iju+∈ijq=Cijkl∈klu,q.The variational formulation for a fixed transverse distribution functionw,is as follows:(28)ux,q¯x∈Uad,Eu,q¯;v,y¯=Lv,y¯.For everyv,y¯, ∈Uad with (29)Uad=HΥ01SxHΥ01SxHΥ02SxHΥ01S.LetAαβδτ=λ¯GαβGδτ+2μ¯GατGτβThen(30)Eu,q¯;v,y¯=∫S∫h/2h/2σiju,q∈ijv,wy¯ψdzdS=∫S∫h/2h/2Aαβδτϵδτuϵαβv+wAαβδτΥδτϵαβuy¯+λ¯w′Gαβϵαβuy¯+λ¯Gαβw′+wAαβδτΥδτq¯ϵαβvψdzdS+∫S∫h/2h/2w2AαβδτΥδτΥαβ+λ¯ww'GαβΥαβq¯y¯+λ¯+2μ¯w′2q¯y¯+μ¯Gαβw2∂βq¯∂αy¯ψdzdS.Let(31)D0wαβ=∫h/2h/2wAαβδτΥδτψdz;D1wαβ,=∫h/2h/2zwAαβδτΥδτψdz;D2wαβ=∫h/2h/2z2wAαβδτΥδτψdz;,E0wαβ=∫h/2h/2λ¯Gαβw′ψdz;E1wαβ,=∫h/2h/2zλ¯Gαβw′ψdz;E2wαβ=∫h/2h/2z2λ¯Gαβw′ψdz,F0wαβγ=∫h/2h/2wAαβδτΥδτ+λ¯Gαβw′ψdz;F1wαβγ∫h/2h/2zwAαβδτΥδτ+λ¯Gαβw′ψdz;F2wαβγ,=∫h/2h/2z2wAαβδτΥδτ+λ¯Gαβw′ψdz;,Fw033γ=∫h/2h/2w2AαβδτΥδτΥαβ+λ¯ww′Gαβψdz;Fw133γ,=∫h/2h/2λ¯wGαβΥαβw′ψdz;Iwwαβγ=∫h/2h/2μ¯Gαβw2ψdz.Then an equivalent form of the variational equation is as follows:(32)Eu,q¯;v,y¯=∫SNαβueαβv+MαβuKαβv+M∗αβuQαβvdS+∫SD0wαβeαβu−D1wαβKαβu+D2wαβQαβuy¯dS+∫SE0wαβeαβu−E1wαβKαβu+E2wαβQαβuy¯dS+q¯∫SF0wαβγeαβv−F1wαβγKαβv+F2wαβγQαβvdS+∫Sq¯y¯Fw033γ−Fw133γ+Fw233γdS+∫SIwwαβ∂βq¯∂αy¯dS,where by letting ψ≈1,Gαβ≈aαβ in calculating the efforts, E,ν¯ Young’s modulus and Poisson’s ratio, respectively,(33)Nαβ=Eh1−v¯21−v¯eαβu+v¯1−v¯1−2v¯eρρuaαβ+Eh3121−v¯21−v¯Qαβu+v¯1−v¯1−2v¯Qρρuaαβ,Mαβ=Eh3121−v¯21−v¯Kαβu+v¯1−v¯1−2v¯Kρρuaαβ,M∗αβ=Eh5801−v¯21−v¯Qαβu+v¯1−v¯1−2v¯Qρρuaαβ+Eh3121−v¯21−v¯eαβu+v¯1−v¯1−2v¯eρρuaαβ,λ¯=Ev¯1−2v¯1+v¯;μ¯=E21+v¯,Lv,y¯=∫SP¯ivi+P¯4y¯dS+∫γ1qivi+q4y¯dγ+∫γ1mαθαvdγ,where(34)P¯4=∫h/2h/2f3wdz+wh2p¯+3−w−h2p¯−3.We use iso-parametric Lagrangian symplectic type triangle finite elements which are based on complete polynomial bases [23]. The finite element used has 4 (four) degrees of freedom at the vertex nodes (3 displacements and 1 stretching) and 1 (one) additional degree of freedom at mid-points of edges (1 for transverse displacement). By laying the matrix of shape functions we get(35)N=1−ξ−ηξη0000000000000001−ξ−ηξη000000000000000−1−ξ−η1−21−ξ−η−ξ1−2ξ−η1−2η41−ξ−ηξ4ξη41−ξ−ηη0000000000000001−ξ−ηξη,U^te=u^1u^2u^3v^1v^2v^3w^1w^2w^3w^4w^5w^6q^1q^2q^3.


### 4.1. Hemisphere under Diametrically Opposed Loads

This test case makes it possible to verify the behavior of the element in bending and shear [26]. The test of a thin hemisphere (R/h = 250) subjected to a free base with four concentrated charges (see Figure 2) is used to verify the absence of membrane locking [21]. The reference solution presented by [27] gives displacements according to the direction of the load: UA=VB=0.094mm. Due to symmetry considerations, only a quarter of the hemisphere is meshed. The mesh is regular and the number of elements varies from 2 to 15. We compare our results obtained with the MT6 or T6-m to several families of finite elements of shells. In Figure 3, the displacement convergence at UA according to the number of elements is represented in Table 1. This new model of thick shell does not have a locking problem. In addition, it is memory less greedy and fast converging with very few elements, as shown in Figure 3 for a thin shell.Figure 2
Hemispheric shell benchmark.Figure 3
(a) Convergence of displacements at point A(3D Model), (b) convergence of displacements at point A (2D Model).
(a)(b)Table 1
Hemisphere under diametrically opposite loads, displacementUA×103.
Reference:UA×103=94mmFamily of 3D elementsFamily of 2D elementsMesh step 3D modelPRI6SHB6 initialSHB6Mesh step 2D modelT6_M3×(5 × 5 × 1)×20,1040,5891,602 × 2 × 21, 83×(11 × 11 × 1)×20,4571,2094,274 × 4 × 210, 273×(22 × 22 × 1)×22,0155,42418,056 × 6 × 231, 923×(180 × 180 × 1)×250,7264,1884,708 × 8 × 254, 353×(250 × 250 × 1)×264,3574,8288,4010 × 10 × 272, 3312 × 12 × 287, 2615 × 15 × 293, 616 × 16 × 294. 8


### 4.2. Pinched Cylinder

The aim of this test is to check the behavior of the element (cylindrical shell) in bending and shearing. It has been studied by several authors ([21, 27–30]). The two ends of the shell are closed by an infinitely rigid diaphragm. An eighth of the cylinder is meshed thanks to the symmetries of the problem. The shell is subjected to a concentrated load at point A, P⟶=1 (Figure 4). Geometry, loading, and material data for this test case are listed in Table 2. And comparison is provided in Table 3.(36)W¯C=−WC.E.hP=164.24.Figure 4
Geometry and loading of the clamped cylindrical shell with diaphragm.Table 2
Geometry and mechanical data of the pinched cylinder.
LengthL=6mRadiusR=3mThicknessh=0.03mYoung’s modulusE=3×1010PaPoisson’s coefficientv=0.3Limit conditionsU=W=0on ADConditions of symetryW=0on AB;V=0on BC;U=0on  CDSollictation at Cf=−0.25NTable 3
Transverse displacements resultsWC at point C.
Reference:W¯C=−WC.E.h/P=164.24Mesh step 3D modelPRI6SHB6 initialSHB6Mesh step 2D modelT6_M(10 × 10 × 1)×210,13416,33539,4112 × 2 × 29,79(30 × 30 × 1)×235,91957,33107,824 × 4 × 243,32(50 × 50 × 1)×260,69689,973137,796 × 6 × 274,08(70 × 70 × 1)×280,406112,95151,838 × 8 × 296,09(90 × 90 × 1)×295,31129,06159,310 × 10 × 2116,4712 × 12 × 2135,2315 × 15 × 2155,0116 × 16 × 2158,7518 × 18 × 2163,38


### 4.3. Sphere under Uniform Pressure

The geometric and mechanical parameters of the hollow sphere under consideration are defined in Table4.Table 4
Geometry and mechanical data of the sphere under uniform pressure.
Inner radiusRi=10cmOuter radiusRe=20cmThicknessh=Re−RiYoung’s modulusPoisson’s coefficientE=21MPav=0.3Inner pressurePi=0.1MPaOuter pressurePe=0MPaThe solution is obtained by solving the equilibrium equation in the absence of volume force:(37)∂σrr∂r−2rσrr−σθθ=0,where σrr,σθθ are the radial and circumferential components of the stress tensor.It is shown that the solution is in the formur=Arλ1+Brλ2 [31].With:(38)λ1=−12−1+λ¯2μ¯+λ¯+12−1+λ¯2μ¯+λ¯2−4−2−λ¯2μ¯+λ¯,λ2=−12−1+λ¯2μ¯+λ¯−12−1+λ¯2μ¯+λ¯2−4−2−λ¯2μ¯+λ¯,B=−2μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−12μ¯+λ¯λ2Reλ2−1+λ¯Reλ2−1A,A=−Pi2μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−1−2μ¯+λ¯λ2Riλ2−1+λ¯Riλ2−12μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−1B.The superposition curves shown in Figures5(a)–5(f) and Table 5 show the deviations.Figure 5
(a) Evolution of radial deviation as a function of thickness ratio (khi = 0.0198). (b) Evolution of radial deviation as a function of thickness ratio (khi = 0.0769), (c) evolution of radial deviation as a function of thickness ratio (khi = 0.1818). (d) Evolution of radial deviation as a function of thickness ratio (khi = 0.4), (e) evolution of radial deviation as a function of thickness ratio (khi = 0.667). (f) Evolution of the radial deviation of the mid-surface as a function of the thickness ratio khi.
(a)(b)(c)(d)(e)(f)Table 5
Radial deviation of the shell (analytical and numerical model) according to the ratioh/R, h is the thickness.
Ratioh/R0.01980.07690.18180.40.667zAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/A−h/210,974610,8910.762,85542,48113.111,2411,12529.330,61290,59692.600,42220,40124.97−h/410,949710,8820.622,82522,47312.471,2001,12276.500,54890,57685.090,32470,381517.49010,925310,8710.502,79762,46212.001,1671,11254.700,50520,55549.940,27190,360232.47h/410,901610,8650.342,77252,43512.171,1401,10143.390,47520,527410.980,24210,339540.23h/210,878410,8420.332,74962,41912.021,1171,08413.010,4550,502410.420,22540,315439.92


## 4.1. Hemisphere under Diametrically Opposed Loads

This test case makes it possible to verify the behavior of the element in bending and shear [26]. The test of a thin hemisphere (R/h = 250) subjected to a free base with four concentrated charges (see Figure 2) is used to verify the absence of membrane locking [21]. The reference solution presented by [27] gives displacements according to the direction of the load: UA=VB=0.094mm. Due to symmetry considerations, only a quarter of the hemisphere is meshed. The mesh is regular and the number of elements varies from 2 to 15. We compare our results obtained with the MT6 or T6-m to several families of finite elements of shells. In Figure 3, the displacement convergence at UA according to the number of elements is represented in Table 1. This new model of thick shell does not have a locking problem. In addition, it is memory less greedy and fast converging with very few elements, as shown in Figure 3 for a thin shell.Figure 2
Hemispheric shell benchmark.Figure 3
(a) Convergence of displacements at point A(3D Model), (b) convergence of displacements at point A (2D Model).
(a)(b)Table 1
Hemisphere under diametrically opposite loads, displacementUA×103.
Reference:UA×103=94mmFamily of 3D elementsFamily of 2D elementsMesh step 3D modelPRI6SHB6 initialSHB6Mesh step 2D modelT6_M3×(5 × 5 × 1)×20,1040,5891,602 × 2 × 21, 83×(11 × 11 × 1)×20,4571,2094,274 × 4 × 210, 273×(22 × 22 × 1)×22,0155,42418,056 × 6 × 231, 923×(180 × 180 × 1)×250,7264,1884,708 × 8 × 254, 353×(250 × 250 × 1)×264,3574,8288,4010 × 10 × 272, 3312 × 12 × 287, 2615 × 15 × 293, 616 × 16 × 294. 8


## 4.2. Pinched Cylinder

The aim of this test is to check the behavior of the element (cylindrical shell) in bending and shearing. It has been studied by several authors ([21, 27–30]). The two ends of the shell are closed by an infinitely rigid diaphragm. An eighth of the cylinder is meshed thanks to the symmetries of the problem. The shell is subjected to a concentrated load at point A, P⟶=1 (Figure 4). Geometry, loading, and material data for this test case are listed in Table 2. And comparison is provided in Table 3.(36)W¯C=−WC.E.hP=164.24.Figure 4
Geometry and loading of the clamped cylindrical shell with diaphragm.Table 2
Geometry and mechanical data of the pinched cylinder.
LengthL=6mRadiusR=3mThicknessh=0.03mYoung’s modulusE=3×1010PaPoisson’s coefficientv=0.3Limit conditionsU=W=0on ADConditions of symetryW=0on AB;V=0on BC;U=0on  CDSollictation at Cf=−0.25NTable 3
Transverse displacements resultsWC at point C.
Reference:W¯C=−WC.E.h/P=164.24Mesh step 3D modelPRI6SHB6 initialSHB6Mesh step 2D modelT6_M(10 × 10 × 1)×210,13416,33539,4112 × 2 × 29,79(30 × 30 × 1)×235,91957,33107,824 × 4 × 243,32(50 × 50 × 1)×260,69689,973137,796 × 6 × 274,08(70 × 70 × 1)×280,406112,95151,838 × 8 × 296,09(90 × 90 × 1)×295,31129,06159,310 × 10 × 2116,4712 × 12 × 2135,2315 × 15 × 2155,0116 × 16 × 2158,7518 × 18 × 2163,38


## 4.3. Sphere under Uniform Pressure

The geometric and mechanical parameters of the hollow sphere under consideration are defined in Table4.Table 4
Geometry and mechanical data of the sphere under uniform pressure.
Inner radiusRi=10cmOuter radiusRe=20cmThicknessh=Re−RiYoung’s modulusPoisson’s coefficientE=21MPav=0.3Inner pressurePi=0.1MPaOuter pressurePe=0MPaThe solution is obtained by solving the equilibrium equation in the absence of volume force:(37)∂σrr∂r−2rσrr−σθθ=0,where σrr,σθθ are the radial and circumferential components of the stress tensor.It is shown that the solution is in the formur=Arλ1+Brλ2 [31].With:(38)λ1=−12−1+λ¯2μ¯+λ¯+12−1+λ¯2μ¯+λ¯2−4−2−λ¯2μ¯+λ¯,λ2=−12−1+λ¯2μ¯+λ¯−12−1+λ¯2μ¯+λ¯2−4−2−λ¯2μ¯+λ¯,B=−2μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−12μ¯+λ¯λ2Reλ2−1+λ¯Reλ2−1A,A=−Pi2μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−1−2μ¯+λ¯λ2Riλ2−1+λ¯Riλ2−12μ¯+λ¯λ1Riλ1−1+λ¯Riλ1−1B.The superposition curves shown in Figures5(a)–5(f) and Table 5 show the deviations.Figure 5
(a) Evolution of radial deviation as a function of thickness ratio (khi = 0.0198). (b) Evolution of radial deviation as a function of thickness ratio (khi = 0.0769), (c) evolution of radial deviation as a function of thickness ratio (khi = 0.1818). (d) Evolution of radial deviation as a function of thickness ratio (khi = 0.4), (e) evolution of radial deviation as a function of thickness ratio (khi = 0.667). (f) Evolution of the radial deviation of the mid-surface as a function of the thickness ratio khi.
(a)(b)(c)(d)(e)(f)Table 5
Radial deviation of the shell (analytical and numerical model) according to the ratioh/R, h is the thickness.
Ratioh/R0.01980.07690.18180.40.667zAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/AAnalytical (mm)Numerical (mm)Error (%) (A-N)/A−h/210,974610,8910.762,85542,48113.111,2411,12529.330,61290,59692.600,42220,40124.97−h/410,949710,8820.622,82522,47312.471,2001,12276.500,54890,57685.090,32470,381517.49010,925310,8710.502,79762,46212.001,1671,11254.700,50520,55549.940,27190,360232.47h/410,901610,8650.342,77252,43512.171,1401,10143.390,47520,527410.980,24210,339540.23h/210,878410,8420.332,74962,41912.021,1171,08413.010,4550,502410.420,22540,315439.92


## 5. Discussions

The shell model, the “N” model [18], is used for the calculation of the hemisphere (Figure 2 and (Table 1), the pinched cylinder (Figure 4 and Table 2) and the sphere under uniform pressure (Table 4). The finite element of the triangular type uses stretching as a degree of freedom to follow the displacement in the thickness. Several authors have worked on the hemisphere case using triangular, rectangular, and even volume finite elements, where all displacement components are represented by polynomials of order 2 and 3. In [26], it is indicated that the results obtained with the new SHR6 element converge faster than the initial PRI6 and SHR6 does not lock for the same number of elements. However, this 3D element has an excessively high number of degrees of freedom as the thickness of the shell increases. It would therefore require a significant computational effort to perform the analysis using this element. In comparison, the present element uses only four degrees of freedom per vertex node and one for middle nodes. The hemispherical shell described above is analyzed here using the triangular element presented above. The results from the analyses are compared with those obtained from the solution given in [27]. Convergence tests were carried out for the normal deviation at point A (Figure 3) of the shell. Table 1 shows that MT6 with the model “N” indicates a deviation of 0.42% for a mesh of 15 × 15 × 2, which is reached by CSFE3-sh with a mesh of 12 × 12 × 2 to converge to the reference solution. The DKT12 and DKT18 give 2% and 1.06% with a mesh of 12 × 12 × 2 for the reference solution [20]. As for the SHR6 developed in [26] the difference obtained is 4.49% for a mesh of 3×(250 × 250 × 1) × 2 or 375 000 elements. We realize that compared to the “N” model, this volumetric model [26], which also takes into account the variation of thickness, is very memory greedy and calculating time is excessive. Figure 3 shows an excellent agreement between the results obtained with the element MT6 and those found in the literature for the variation of the normal deviation at point A.The pinched cylinder is a reference example for testing a finite element model, especially for thin shells. The loading case here is not symmetrical (due to the concentrated point load) and the solution cannot be managed by the 3D elasticity theory [20]. The “N” model although being a thick shell model, shows good convergence towards the reference solution. The error is 0.523% for a mesh of 18 × 18 × 2, while the SHR6 of [26] indicates an error of 3% for a mesh of (90 × 90 × 1) × 2 or 16200 elements. The convergence curves are illustrated in (Figures 6(a) and 6(b)). It is a difficult problem to solve numerically when the thickness is small because the locking and hourglass problems are recurrent. The finite element of [26] does not lock, but it needs a very fine mesh and fairly important calculating time.Figure 6
(a) Convergence of displacements at point C (3D model) (b) convergence of displacements at point C (2D model).
(a)(b)The sphere case under uniform pressure is a resolvable case by the 3D elasticity theory. The reference solution, which is a function of the radius, is obtained analytically. The model “N” (by Nzengwa [18]) has the ability to also calculate the radial displacement as a function of the radius. In Table 5, we plotted the following values of the characteristic parameter: h/R = 0.0198, 0.0769, 0.1818, 0.4, and 0.667. A gap of less than 5.55% was observed between the inner and outer radii for the radial deviation for thin shells (example h/R = 0.0198). In general, the average gap is around 10%, as shown in (Figures 5(a)–5(e)). The numerical displacements are smaller than the analytical ones, which highlights the contribution of the flectional energy of the double curvature shell. Table 5 shows that the model is closer to the reference solutions according to the ratio khi for a uniform loading. For nonuniform loads, the 3D elasticity theory can no longer be able to provide analytical solutions, while the model “N” still remains efficient for solving common problems in complex shell structures, either thin or thick.


## 6. Conclusion

The 2D model “N” which is a 4-parameter model has proven its efficiency in handling different loadings in thin or thick shells with the calculation of transverse stresses, strains, and also thickness variation. Increasing the degree of the polynomial of the transverse stretch distribution functionw(z) does not increase the number of parameters. The transverse distribution function w(z) is not limited to polynomials only. The model includes terms found in the classical theory of thin or thick shells and seems to be suitable for stiffened shells that are actually locally thin or thick shell structures. The distribution function w(z) in this paper is z, so the quality of the results can be improved by choosing another expression, which makes the model flexible.



---

*Source: 5485130-2022-10-18.xml*