An advanced shell model for the analysis of geometrical and material nonlinear shells

In this paper layered shells subjected to static loading are considered. The displacements of the Reissner–Mindlin theory are enriched by a an additional part. These so-called fluctuation displacements include warping displacements and thickness changes. They lead to additional terms for the material deformation gradient and the Green–Lagrangian strain tensor. Within a nonlinear multi-field variational formulation the weak form of the boundary value problem accounts for the equilibrium of stress resultants and couple resultants, the local equilibrium of stresses, the geometrical field equations and the constitutive equations. For the independent shell strains an ansatz with quadratic shape functions is chosen. This leads to a significant improved convergence behaviour especially for distorted meshes. Elimination of a set of parameters on element level by static condensation yields an element stiffness matrix and residual vector of a quadrilateral shell element with the usual 5 or 6 nodal degrees of freedom. The developed model yields the complicated three-dimensional stress state in layered shells for elasticity and elasto-plasticity considering geometrical nonlinearity. In comparison with fully 3D solutions present 2D shell model requires only a fractional amount of computing time.


Introduction
Shell elements which account for the layer sequence of a laminated structure are able to predict the deformation behaviour of the reference surface in an accurate way. Also the assumption of a linear shape of the in-plane strains through the thickness is sufficient, if the shell is not too thick. In contrast to that only constant transverse shear strains through the thickness are obtained within the Reissner-Mindlin theory. As a consequence only the average of the transverse shear stresses is accurate. Neither the shape of the stresses is correct nor the boundary conditions at the outer surfaces are ful- filled. Within the Kirchhoff theory the transverse shear strains are set to zero. By assumption the thickness normal stresses are neglected in a standard shell theory. This is necessary to avoid unphysical stresses due to inextensibility assumptions in thickness direction.
In several publications the equilibrium equations are exploited within post-processing procedures to obtain the interlaminar stresses, e.g. [1,2] for the transverse shear stresses and e.g. [3] for the thickness normal stresses. In a standard version first and second derivatives of the in-plane stresses require bi-quadratic or bi-cubic shape functions. The essential restriction of the approach is the fact that these stresses are not embedded in the variational formulation and an immediate extension to geometrical and physical nonlinearity is not possible.
The standard displacement field is enhanced by layerwise (zig-zag) functions through the thickness in e.g. [4][5][6]. An enhanced modeling approach for multilayer anisotropic plates based on the dimension reduction method is presented in [7]. An assessment of the variable separation method is discussed in [8].
The use of brick elements or so-called solid shell elements represents a computationally expensive approach, e.g. [18,19]. For a sufficient accurate evaluation of the interlaminar stresses each layer must be discretized with several elements (≈4-10) in thickness direction. Especially for large scale industrial problems with a multiplicity of load steps and several iterations in each load step this is not a feasible approach.
The purpose of this paper is to present a shell model which is able to compute the load-deflection behaviour and the complicated three-dimensional stress-state of geometrical nonlinear, elasto-plastic, layered shells. It is based on Refs. [20][21][22]. Present work delivers essential new developments.
(i) The displacements of the Reissner-Mindlin kinematics are enriched by warping displacements and thickness changes. This yields additional terms for the material deformation gradient, the Green-Lagrangian strain tensor and associated variation. The thickness integration of the local equilibrium equations are extended to elastoplasticity. (ii) The proposed nonlinear variational formulation leads to Euler-Lagrange equations which include besides the usual shell equations in terms of stress resultants, the local equilibrium in terms of stresses, a constraint which enforces the correct shape of the displacement fluctuations through the thickness, the geometric field equations and the constitutive equations. (iii) A new ansatz for the independent shell strains with quadratic shape functions is proposed. It leads to a significant improved convergence behaviour especially for distorted meshes. Furthermore, quadratic functions allow the computation of second derivatives. These are necessary when rewriting the local equilibrium equation for the thickness direction with the derivatives of the transverse shear strains. (iv) Elimination of a certain set of parameters by static condensation from the resulting system of equations yields the element stiffness matrix and the element residual vector. The derived four-node element possesses the usual 5 or 6 nodal degrees of freedom. This is an essential feature since standard geometrical boundary conditions can be applied and the element is applicable also to shell intersection problems. In comparison with fully 3D discretizations present 2D shell model requires only a fractional amount of computing time.
The paper is organized as follows. In Sect. 2 the shell theory is derived. The associated finite element formulation is presented in Sect. 3. Several geometrical and physical nonlinear examples are presented in Sect. 4. The computed displacements and stresses of the developed model are compared with costly 3D computations. Furthermore, we compare with the displacements and stresses computed with Reissner-Mindlin shell elements.

Shell kinematics
Let B 0 be the three-dimensional Euclidean space occupied by the shell with thickness h in the undeformed configuration. The position vector X of any point P ∈ B 0 is parametrized with convected coordinates ξ i where X 0 and N denote the position vector and normal vector of the reference surface , respectively. The thickness coordinate ξ 3 is defined in the range h − ≤ ξ 3 ≤ h + , where h − and h + denote the coordinates of the outer surfaces. The coordinate on the boundary = u ∪ σ of is s. The usual summation convention is used, where Latin indices range from 1 to 3 and Greek indices range from 1 to 2. Commas denote partial differentiation with respect to the coordinates ξ i . The kinematic assumption for the position vector of the deformed shell space B reads Here, x 0 describes the position vector of the current reference surface t . The director vector d(ϕ) of the current configuration is not perpendicular to the deformed reference surface, thus shear deformations are accounted for. We introduce the vector field v(ξ 1 , ξ 2 ) = [u, ϕ] T , where u = x 0 − X 0 and ϕ contains rotational parameters. In Eq. (2) the assumptions of the Reissner-Mindlin theory are extended by the displacement fluctuation field For this purpose the shell is subdivided in thickness direction in M numerical layers. For laminated shells M usually corresponds to the number of physical layers. The vector α contains displacements at nodes in thickness direction and is element-wise constant. The matrix contains layer-wise cubic hierarchic functions where a i is an assembly matrix, which relates the 12 degrees of freedom of layer i to the K components of α. For M layers this leads to K = 9 · M + 3. Furthermore, ζ is a normalized coordinate of layer i defined in the range −1 ≤ ζ ≤ 1 and 1 n denotes a unit matrix of order n.
With the parametrization of the shell (1) and kinematic assumption (2) one can express the material deformation gradient as with the contravariant metric coefficients G i j and The vectorũ, 3 follows fromũ, 3 where h i denotes the thickness of layer i. Inserting the deformation gradient into the Green- where the index g refers to geometrical strains. The components read with (6) considering d · d = N · N = 1 as well as N, α ·N = d, α ·d = 0 with The higher order curvatures ρ αβ are neglected. This is admissible for sufficiently thin structures with L/h 1, see e.g. [23,24]. Here, L is a characteristic length of a plate or lowest curvature radius of a shell. Hence, the in-plane strains {E 11 , E 22 , 2E 12 } are linear functions of ξ 3 . Our numerical investigations show that for L/h ≥ 10 this is a good approximation. Using Voigt notation the Green-Lagrangian strains of a point in shell space with coordinate ξ 3 are obtained with The components of F and E g are transformed to the elementwise constant Cartesian element coordinate system t i , see Sect. 3. The transformations are standard and therefore are not displayed here. For the below described variational formulation the virtual Green-Lagrangian strain tensor has to be specified Using Voigt notation it holds with δε g = [δε 11 , δε 22 , 2δε 12 , δκ 11 , δκ 22 , 2δκ 12 , δγ 1 , The below presented finite element equations are iteratively solved. To maintain quadratic convergence in the Newton-Raphson scheme the vectors [d, g 1 , g 2 ] of the last converged load increment are used inÃ 2 and A 2 . In order to assess the effect of this approximation a variation of the load step size has to be performed. This is done by means of the numerical examples in Sect. 4. In this context also the second variation of the Green-Lagrangian strains has to be derived. One obtains The linearized virtual shell strains δε g are specified e.g. in Ref. [25].

Equilibrium equations and a constraint
This part follows Ref. [22] considering the new relations of Sect. 2.1. Furthermore, the extension to elasto-plastic material behaviour is derived. In this section we use Voigt notation. The components of the stresses, strains and material laws refer to the element-wise constant Cartesian coordinate system t i , see Sect. 3.

Definition of stress resultants
The transformation of the Second Piola-Kirchhoff stress tensor S = S T to the First Piola-Kirchhoff stress tensor P is obtained in a standard way by P = F S. Thus, with Figure 1 depicts static loading p − and p + acting at the outer surfaces of the shell with coordinates ξ 3 = h − and ξ 3 = h + , respectively. The stress boundary conditions at the outer surfaces read The components of the Second Piola-Kirchhoff stress tensor are organized in a vector using Voigt notation In an analogous way the physical strains E and associated variations δE are introduced. The components are organized in vectors as the geometrical counterparts in Eqs. (12) and (14). Now, the relation of S to the stress resultants is defined by thickness integration of the internal virtual work per unit area. This yields with δE = A 1 δε + A 2 δα, where δε is the vector of the virtual independent shell strains with the determinant of the shifter tensorμ and The components of ∂ ε W = [n 11 , n 22 , n 12 , m 11 , m 22 , m 12 , q 1 , q 2 ] T are membrane forces, bending moments and shear forces, whereas the components of ∂ α W are higher order stress resultants. Remark: In contrast to the independent stress resultants denoted by σ , we use for the stress resultants (21) computed with S via the constitutive law the notation ∂ ε W . Thus, the notation ∂ ε W and ∂ α W has for nonlinear material laws not the meaning of derivatives of a strain energy density with respect to ε and α.

Thickness integration of the equilibrium equations
Neglecting body forces the component representation of the equilibrium equation Div P = P kl , l t k = 0 reads Subsequently, the integral form of f = 0 is derived. This is done in two steps. At first the derivatives of the stresses P iα are reformulated. In a second step the terms P i3 , 3 are treated. The derivatives of the stresses P iα with respect to ξ α using (17) yield In-plane stresses S αβ as well as transverse shear stresses S 3α enter in this equation.
In the following the term S 31 , 1 +S 32 , 2 is reformulated. For this purpose the first two equations of (22) are rewritten with (17) and (23) as followŝ For arbitraryF = 0 and p α p * 3 the transverse shear stresses follow from (24) as (25) and the sum of the derivatives yields In a first step we assume elastic orthotropic material behaviour. Hence, the strain energy density (E) = 1 2 E T C E can be written in terms of the constant elasticity matrix C and the physical strains E. They are obtained with the physical shell strains ε via E = A 1 ε +Ã 2 α . Hence, the constitutive equation may be written in the following standard form Due to the varying fibre orientation in laminated shells the constants C i j = C ji differ for each individual layer. For the in-plane stresses holds with E = A 1 ε +Ã 2 α and Eq. (27) We introduce the matrices For a symmetric shape of C i j through the thickness and With layerwise constant C i j the integration in (29) 2 and (30) can be done by summation over layers and analytical integration in each layer.
The derivatives of the stresses P iα (23) yield with (26)- where The derivatives of E 33 follows from (9) 3 and yield E 33 , β = g 3 , β ·ũ 3 = d, β ·ũ 3 . However, the numerical tests show that the second part in (32) leads to negligible contributions. The vector λ contains derivatives of the membrane strains and where the reformulation of ∂ α W withμ = 1 according to (21) 2 is obtained with integration by parts. Here, one needs the relation with A 2 and P i3 from (14) and (17), respectively. Considering the stress boundary conditions (18) one obtains for the boundary term where p − = p − N and p + = p + N. Furthermore, 0 n denotes a zero vector with n components. The weighted integral of p * is constant as is taken from the last load increment. Therefore it can be summarized with q Now the integral form of f = 0 (22) is formulated with δũ = δα and Eqs. (31)- (37) With δα = 0 one obtains the equilibrium of higher order stress resultants

Thickness integration of a constraint
The displacement fieldũ must fulfill a constraint. To specify this equation we introduce the equilibrium of virtual stresses δ P iα . As C 23 and p * are constant in Eq. (31) one obtains The integral form of With Eq. (41) reads −δλ T D 32 α = 0 where δλ = 0. One obtains the constraint which enforces the correct shape of the displacementsũ through the thickness. It has the descriptive meaning that the superposed warping displacements and thickness changes must not lead to additional stress resultants.

Extension to elasto-plastic material behaviour
The extension to elasto-plastic material behaviour is possible. In the following we consider small strain plasticity with an additive decomposition of the total strains E = E el + E pl in elastic and plastic part and Eq. (27) is replaced by The T leads to further contributions in p * according to (32). It has to be replaced by The vector p * depends on the deformation. To maintain quadratic convergence in the Newton iteration process, the constant values of the last load increment are used. As already written above this requires a variation of the load step size to assess this approximation.

Linearized weak form of the boundary value problem
In the following the weak form of the boundary value problem and associated linearization is derived. For a compact representation the vector θ and admissible variation δθ are introduced Furthermore, summarizes the constitutive equation, the equilibrium of higher order stress resultants (39) and the constraint (43). Hence, the weak form of the boundary value problem can be written as The shell is loaded statically by surface loads p + = p + N and p − = p − N acting at the outer surfaces, see Fig. 1 and by boundary forcest on the boundary σ . Hence, the virtual work of the external forces reads With application of integration by parts to the integral of δε T g σ in Eq. (48) and use of standard arguments of variational calculus one obtains as Euler-Lagrange equations the equilibrium of stress resultants and couple resultants, the geometric field equation ε g − ε = 0 and ψ = 0 in along with the static boundary conditions t =t on σ .
The associated nonlinear finite element equations are iteratively solved using Newton's method. For this purpose variational equation (48) is linearized. With A 2 = A 2 (α) one has to consider Eq. (16). The following matrices are introduced and D 23 according to (34). The matrix C T = ∂ E S denotes the algorithmic tangent matrix. The singular matrix D = D T is of order 29 + K . The integration of the submatrices D i j is performed by summation over M layers and three point Gauss integration for each layer. With dξ 3 = h i 2 dζ one obtains where C i j denotes the respective integrand of D i j . Furthermore, ζ j and W j are the normalized thickness coordinate and weighting factor of the particular Gauss point, respectively. For the determinant of the shifter tensor the first approximationμ = 1 is taken. In case of linear elasticity C T is constant and the numerical integration is exact.
With displacement independent loadsp,m,t and consideration of (50) one obtains the linearized weak form The second variation δd of the current director vector has been derived e.g. in Ref. [25].

Finite element formulation
The approximation of initial and current geometry of the shell reference surface applying the isoparametric concept for 4node elements is specified in detail in Refs. [25,26] Hence, the Jacobian matrix J reads The superscript h refers to the finite element approximation of the particular quantity, and commas denote the partial derivative with respect to ξ or η.
The shape of the quantities δθ h is chosen as follows To avoid shear locking, ansatz functions of the assumed strain method [27] are incorporated in B, see Ref. [25].
The matrix N σ for the interpolation of δσ h = N σ δσ as well as σ h is chosen as follows The coefficient matrices T 0 σ ,T 0 σ describe a transformation of contravariant tensor components to the constant base system t i . The matrices and the constantsξ andη are specified in [25]. The parameter vectorsσ and δσ contain 8 parameters for the constant part and 6 parameters for the varying part of the stress field. The interpolation of the membrane forces and bending moments corresponds to the membrane part in Ref. [28]. The original approach for plane stress problems was published in Ref. [29]. Regarding requirements on the interpolation functions to fulfill the patch test and to ensure stability of the discrete system of equations we refer to the discussion in Ref. [25].
The interpolation for ϑ and equivalent δϑ is chosen as follows whereε ∈ R 8 ,λ 1 ∈ R 6 ,λ 2 ∈ R 4 ,λ 3 ∈ R l , l = n + 9 and α ∈ R K . The submatrices associated with ε h read The coefficient matrices T 0 ε andT 0 cause a transformation of contravariant tensor components to the constant element base system t i . The entries J 0 αβ are the components of J according to (54) evaluated at the element center. Detailed investigations on the use of ansatz functions for contravariant stress and strain components in the framework of a Hu-Washizu functional are contained in Ref. [30].
The interpolation matrix M m n is chosen as where the index n ∈ {0, 2, 4, 6, 7, 9, 11} has the meaning that optionally the first n columns are taken. With n = 0 the matrix is omitted. Due to the factor j 0 / j and the constant coefficient matrix T 0 ε the integrals of the functions (ξ, η, ξ η, ξ 2 η, η 2 ξ) in N 2 ε over the element domain e vanish.
The shape factor c considers the deviation of the element geometry from a square and warping of the particular element. For this purpose the metric coefficients G αβ of the initial reference surface are evaluated at the element center A principal axis analysis leads to with λ 2 min > 0, as G 0 is positive definite. Introducingx = [x,ŷ] T one obtains for the quadratic formx T G 0x := r 2 =x 2 λ 2 min +ŷ 2 λ 2 max and thus the equation of an ellipsê Here, a and b (a ≥ b) denote the semiaxes of the ellipse which can be inscribed in the flat projection h 0 of a distorted element, see Fig. 2a. Furthermore, the element warping d according to Fig. 2b is computed. One obtains with the position vectors X i of the four nodes and the normal vector t 3 . Hence, c is defined as With N α = 1 K the vector α h is constant within each element. Considering Eq. (57) the interpolation for λ h = [λ h ε1 , λ h κ1 , λ h κ2 ] T is performed with λ h = N 2 λλ 2 +N 3 λλ 3 using the matrices where Concerning the predefinition = 100 h we refer to the investigations in Ref. [20]. The FE approximation of the external virtual work ofp,m andt leads to Here, numel denotes the total number of finite shell elements to discretize the problem and f a corresponds to the element load vector of a standard displacement formulation. Furthermore, it holds where for the special casem · δd = 0 the matrix K g is specified in detail in Ref. [25]. For the general case the term m · δd = 0 leads to a symmetric contribution in K g . The second variation of d has been derived in e.g. Ref. [25] such that the consideration of the load term is standard. Hence, inserting eq. (55) and the corresponding equation for θ h into the linearized variational equation (52) yields The vectorψ corresponds to ψ according to (47) without vector σ . The integrals in (73) are computed numerically applying a 3×3 Gauss integration scheme considering dA = |X h 0 , ξ ×X h 0 , η | dξ dη. It is important to note, that although D is singular, H is regular. The matrix L is expressed with (57) The columns with quadratic shape functions in N 3 ε are not orthogonal to columns 9-12 of N σ and thus lead to entries in L 3 . They are consistently omitted when setting L 3 = 0 in L, f e and f s .
We continue with L[g(θ h , δθ h ), θ h ] = 0 , where δθ h = 0 and obtain for each element where r denotes the vector of element nodal forces. Since σ and ϑ are interpolated discontinuously across the element boundaries the parameters σ and θ can be eliminated from the set of equations. This is done applying a standard Gaussian elimination procedure [31] to the system of equations (75). One obtains the tangential element stiffness matrix k e T , the element residual vectorf e and (72) reduces to where J 0 denotes the Jacobian matrix (54) evaluated at the element center. The derivatives of E pl i j (ξ, η) with respect to ξ and η are approximated using a central difference scheme and thus are elementwise constant. Here, ξ = η = √ 0.6 are the coordinates of the applied 3 × 3 Gauss integration scheme.
Remark: The derivatives of the deformation gradient in (23) with respect to the coordinates ξ α require the computation of second derivatives of the displacements. Thus, at least

Examples
The displacements and stresses obtained with present 4-node shell element are compared with 3D reference solutions. For this purpose we use the 8-noded solid shell element [18] with a sufficient number of elements in thickness direction of the plate or shell. We add results computed with shell element [34]. It is based on standard Reissner-Mindlin kinematic assumptions along with thickness strains which have a constant and linear shape through the thickness. This leads to an interface for three-dimensional constitutive laws. The element was designed for homogeneous shell structures. With application to layered shells it leads to restraints especially for the thickness normal stresses. Furthermore, we compare with literature results. These were obtained with the 4-node shells elements [35] and [36]. The used notation of the different element formulations is summarized in Table 1.

Comparison of relative computing times
In Table 2 relative computing times (stiffness computation and solution of the system of equations) using present 2D shell element and the 3D solid shell element are displayed. The 3D meshes are generated with 4 elements in thickness direction of each layer. This is necessary to map the complicated shape of the stresses through the thickness, see [37]. A fast direct parallel solver is used along with a Windows PC (1 Intel i9-7900X CPU, 10 cores, 3.3 GHz). In each row of the table the 3D times are normalized with respect to the 2D times. The 3D times for the finest meshes are not computable with the used hardware. This shows that fully 3D computations are costly, which restricts the applicability to small problems.

Hemispherical shell
The first problem is the hemispherical shell with an 18 • cutout subjected to alternating radial point loads P at its equator, shown in Fig. 3a. This geometrically non-linear example is often cited as a benchmark problem for shell elements. It is a test for the ability to model rigid body modes and inextensible bending, [33]. Geometrical and material data are R = 10, ϕ = 18 • , thickness h = 0.04 and E = 6.825 · 10 7 , ν = 0.3. Considering symmetry one quarter of the structure corresponding to the region ABCD in Fig.  3a is discretized using regular and distorted meshes. Here, the principal mesh distortion is described in Fig. 3b for a 4 × 4 mesh. Each edge is discretized using the aspect ratios L 1 :  Fig. 3c. We employ the boundary conditions u y = β = 0 on AD, u x = β = 0 on BC and u z = 0 at a point on AB, e.g. at A. Fig. 4 shows load displacement curves for the regular and distorted mesh 8 × 8. The defined converged solution is computed with a 128 × 128 mesh. Results are only presented for P − u x A ; similar output can be obtained for P − u y B . For comparison we add results from Ref. [35] using the MITC4+ element. The convergence behaviour of the displacement u x A at P = 400 versus the number of elements N is depicted in Fig. 5 for both meshes. It is obvious that n = 11 leads to a significant improvement of the element behaviour. The results for n = 9 are close to those for n = 11. In this context we compare the effort to setup the tangential stiffness matrix for different parameters n. The computing time T is normalized to T = 1 for the case n = 0. The comparison with the further parameters leads to T = 1.21 for n = 7, T = 1.28 for n=9 and T = 1.32 for n = 11. The numbers prove that the effort especially for the last two cases is comparable. Within the following examples the parameter n = 9 is not considered. Finally a plot of the shape factor c according to Eq. (66) is shown for a regular and distorted 16 × 16 mesh in Fig. 6.

Twisted beam
The twisted beam problem according to Fig. 7 was originally introduced in [33]. Geometrical and material data are L = 12, b = 1.1, thickness h = 0.0032 and E = 29 · 10 6 , ν = 0.22, respectively. The cantilever beam is clamped at one  Fig. 8. Here, a ratio L max /L min = 2 is chosen, where L max and L min denote the longest and shortest element length in the flat projection, respectively. Figure 9 depicts load deflection curves for the regular mesh and the distorted mesh using different parameters n . Available results using the MITC4+ element [35] are added. The converged solution is obtained employing a 32×192 regular mesh. Again, the quadratic terms in Eq. (61) (n = 11) are necessary to produce accurate results.
For the distorted mesh the convergence behaviour of the final displacements u y A and −u z A versus the number of ele-  Figs. 10 and 11, respectively. Again, n = 11 leads to a significant improvement of the element behaviour. This is only obtained with the shape factor c according to Eq. (66), as the comparison with the curves n = 11 (c = 0) shows.

Hook problem
Next, we consider the hook problem shown in Fig. 12, referred to in linear analysis as the Raasch challenge, [38]. For the FE-discretization we use N × 2 N × 3 N elements with N elements in width direction, 2 N elements for the first arch (radius R 1 ) and 3 N elements for the second arch (radius R 2 ), see Fig. 12.
Geometrical and material data are R 1 = 14, ϕ 1 = 60 • , R 2 = 46, ϕ 2 = 150 • , b = 20, thickness h = 0.02 and E = 3.3 · 10 3 , ν = 0.3, respectively. The structure is fully clamped at one end and is loaded by a shear load P  For the solution, we use a regular and a distorted 4 × 8× 12 mesh. The principal distorted mesh is shown in Fig. 13 with respect to a flat projection together with a perspective view of the structure. Following [36], L max /L min = 1.5 is chosen for the first arch and L max /L min = 2.0 for the second arch.  Figure 14 shows the resulting load-displacement curves P − u z A of point A . The displacements using the MITC4+ element [35] are included. Similar results can be found for the +HW element (see Fig. 12b, 13b of [36]). The defined converged solutions are obtained with a 32 × 64 × 96 regular mesh.    Figure 15 shows the convergence behaviour of the final displacement u z A of point A versus the number of elements N in width direction. Results for the elements MITC4, MITC4+ and +HW are taken from Fig. 12a and 13a in [36]. The supe-

Elasto-plastic analysis of a square plate
A simply supported square plate (soft support) subjected to a constant load is considered. The length and thickness of the plate are L = 200 cm and h = 4 cm, respectively. The origin of the (x, y, z)− coordinate system is chosen at the center of the plate. The material parameters are the elastic constants E and ν as well as yield stress y 0 and linear hardening parameter ξ for J 2 -plasticity. They are chosen as follows E = 21000 kN/cm 2 ν = 0.3 y 0 = 36 kN/cm 2 ξ = 420 kN/cm 2 .
(80) The plate is loaded by a constant load λp withp = 1 N/cm 2 which is applied one half each via the lower and upper surface. With respect to symmetry one quarter of the plate is discretized with a mesh of 20×20 quadrilateral elements. Present element and the Reissner-Mindlin shell element [34] are used with N = 5 layers of equal thickness. For the 3D reference solution of the stresses the mesh has to be refined to obtain sufficient accurate results. We choose a 40 × 40 × 20 mesh. Besides the symmetry conditions the boundary conditions are u z = 0 at (L/2, y, z) and (x, L/2, z).
Geometrical linear and material nonlinear calculations are performed. Load deflection curves for the center displacement w = −u z (0, 0) are shown in Fig. 16. Both shell results agree with the 3D reference solution. The yield line theory, which assumes elasto-plasticity without hardening, predicts a bearing loadp ylt = 6 · y 0 · (h/L) 2 which leads to a load factorλ = 86.4.
For the maximum load the normal stresses S x x and S zz are evaluated near the plate center at the coordinates (x p1 , y p ) =

Elasto-plastic analysis of a sandwich plate
The next example is a square sandwich plate with a layup according to Fig. 23. The origin of the (x, y, z)-coordinate system lies in the center of the plate. The thickness of the core is t c and the thickness of one face layer is t f . All edges are simply supported (soft support). A constant load λp with p = 10 −3 N/mm 2 is applied at the upper surface. The material properties for isotropy are E c , ν c for the elastic core and E f , ν f for the elasto-plastic face layers with yield stress y 0 and linear hardening parameter ξ for a J 2 -plasticity model. One quarter is discretized with a mesh of N × N elements considering symmetry. Two layers are used for the core and one for each face layer. For the 3D solution 6 elements are used in thickness direction of the core and one element for each face layer. All data are summarized as follows. With shear correction factor k = 0.03, computed with element formulation [37], there is good agreement in the elastic range. For inelasticity the factor is applied to the algorithmic tangent matrix. This leads in the inelastic range to a response which is to soft. The stresses are evaluated for load factor λ = 12. To obtain sufficient accurate results especially for the thickness normal stresses the 3D mesh is refined with N = 80. The normal stresses S x x and S zz are computed near the plate center at x p1 = y p = L/80. They are plotted as function of z in Figs. 25 and 26, respectively. The transverse shear stresses are computed at x p2 = 61/80 · L, y p = L/80 and are shown in Fig. 27. In all diagrams there is good agreement between present 2D and the 3D results. This holds for the displace- ments and the computed stresses. The Reissner-Mindlin results for the stresses are not included in the diagrams. As already shown with the previous example only for the inplane stresses and elasticity one obtains usable results. Here, even the load deflection behaviour is not correct in the inelastic range.  The geometrical nonlinear computation is performed applying the load in one single step. The maximum displacement computed with the different models amounts to u z = −145.7 mm. The stresses S x x , S zz , S xz and S yz are evaluated at the coordinates x p = 13/48 · L, y p = 5/48 · L. They are plotted as function of the thickness coordinate z in Figs. 29, 30, 31, and 32. The transverse shear stress S yz according to Fig. 32 shows a very complicate shape. Even the sign of the stresses changes several times. There is a good agreement of present 2D solution with the 3D reference stresses. For comparison we add the results based on the Reissner-Mindlin theory. As already written in the previous examples only the in-plane stresses agree with the reference solution. The interlaminar stresses deviate considerable form the reference solution. The integrals through the thickness lead to correct stress resultants within the chosen mesh density. Finally, the deformed configuration is depicted in Fig. 33.

Stiffened cylindrical shell
As last example the stiffened cylindrical shell according to Fig. 34 is investigated. The figure shows a cross-section of the shell and a coarse finite element mesh of half the structure   Comparative results are computed with a 3D discretization of the skin. In thickness direction of the skin each physical layer is discretized with 2 elements. To obtain sufficient accurate results for the thickness normal stresses an in-plane refined mesh with l × m × n = 144 × 96 × 4 is used. For the discretization of the stiffeners present shell element is used.
The geometrical nonlinear computations are performed with displacement control and a step size w = 20 mm. In Fig. 35 load F is plotted versus the prescribed displacement w.
For two deformed configurations w = 200 mm and w = 400 mm the stresses S 11 , S 33 and S 13 at a point P of the reference surface with coordinates ξ 1 p = (13/96 · π/2) · R, ξ 2 p = x 2 p = 7/64 · L are displayed in Figs. 36, 37 and 38 in dependence of z = ξ 3 . In all diagrams there is good agreement of present element with the reference solution.
For comparison results of the Reissner-Mindlin model are added. Again, one obtains good results for the in-plane stress S 11 , whereas this is not the case for S 33 and S 13 . Figure

Conclusions
An advanced 2D shell model for layered shell structures is presented. With the developed finite element formulation the load deflection behaviour of geometrical and material nonlinear shell problems and the complicated three-dimensional stress state in layered shells can be computed. The matrices A 2 and A 2 (Eqs. (12) and (14) ) are evaluated with the displacements of the last converged load step. This also holds for D 23 and p * (Eqs. (34) 2 and (45) ). All computed examples show by a variation of the load step size that no visible differences in the load deflection curves occur, which means that the applied approximation is admissible. The proposed ansatz for the independent shell strains with quadratic functions leads to a significant improved convergence behaviour especially when meshes with distorted elements are used. The computed displacements and stresses are in good agreement with the results of 3D comparative models. For structures with many layers present 2D shell model requires only a fractional amount of the computing time in comparison with fully 3D discretizations.
The numerical examples furthermore show, that shell elements based on the Reissner-Mindlin theory are able to predict the correct shape of the in-plane stresses only for the case of elasticity. For elasto-plasticity the stresses fulfill the yield condition but not the equilibrium equations in terms of stresses. This holds for all components. As the weak form of equilibrium in terms of stress resultants is basis of the finite element formulation, the integrals of the stresses lead to correct stress resultants and thus to correct load deflection curves for sufficient fine meshes. For sandwich plates and shells the load deflection behaviour can only be predicted with a shear correction factor in the elastic range.
Funding Open Access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.