Dual-mixed hp-version axisymmetric shell finite element using NURBS mid-surface interpolation

A new, general hp-version axisymmetric finite element is derived for the boundary value problems of thin linearly elastic shells of revolution, applying a complementary strain energy-based three-field dual-mixed variational principle. For the interpolation of the mid-surface geometry, non-uniform rational B-splines—NURBS—is used. The independent field variables of the weak formulation are the a priori non-symmetric stress tensor, the displacement vector, and the infinitesimal skew-symmetric rotation tensor. The theoretical model of the shell formulation is based on a consistent dimensional reduction process and a systematic variable-number reduction procedure. The inverse of the unvaried three-dimensional constitutive equation is employed since neither the classical kinematical assumptions nor the stress hypotheses are built in the mathematical model; namely, both the through-the-thickness variation and the normal stress to the shell mid-surface are not excluded. The new hp axisymmetric shell finite element is tested by a representative model problem for extremely thin and moderately thick, singly and doubly curved shells of negative and positive Gaussian curvature. Following from the numerical experiments, the constructed hp-shell finite element gives locking-free results not only for the displacement but also for the stresses.


Introduction
The most commonly used numerical technique for solving the mathematical problems of solid mechanics is the finite element method (FEM). This discretization method built in many commercial software packages is typically based on the principle of virtual work relaxing the translational equilibrium equation and the traction boundary condition (BC). The displacement BC and the rotational equilibrium equation (symmetry condition for the stress tensor) are fulfilled in strong form. In this case, the basic variable is the displacement field [23].
Nevertheless, these standard, widely used, displacement-based finite element (FE) software products are not free from numerical difficulties occurring even in the context of linear elasticity. Most of these FEs implemented in the commercial softwares aim to achieve the required accuracy with the use of mesh refinement technique, i.e., h-version scheme [23]. Inaccuracies and instabilities can be experienced at modeling the deformation and stress state of incompressible materials, yielding the volumetric locking effect [58]. This numerical problem appears in such a way that the convergence for the relative errors of both the stress-and the displacement variables is either very slow or nonexistent.
However, a faster convergence is provided by the p-version approximation technique. In this case, a coarse FE mesh is kept fixed, and the polynomial degree p of the displacement approximation function is increased on each element [14,58]. Although the p-and the combined hp-type displacement-based FEs give sufficiently accurate solutions for the displacement field, slow convergence and low accuracy can be provided for the trace of the stress tensor having more significance in the engineering practice than the displacement vector [14,58].
The classical shell theories are usually derived by the use of the dimensional reduction process. During this procedure, the original three-dimensional (3D) shell equation is replaced with an approximate two-dimensional (2D) one in such a way that all mechanical quantities of interest, i.e., displacements, strains, and stresses, are expanded into series in the thickness direction with respect to an appropriate function basis, which can be monomials, scaled Legendre-type polynomials or trigonometric functions. Then, the series expansions are truncated at various orders, generating a hierarchy of families of shell models. The well-known Koiter-and Naghdi-type shell theories are the first-and second-order members of the displacement-based, approximate shell models [30,36], which are originally elaborated on plates by Kirchhoff and Love, as well as Reissner and Mindlin, respectively, in [33,52]. The consistency of these shell models and many other ones is investigated in [24,25,53,54]. Additional modified, displacement-based shell theories are developed in [1,10,19,34,35,49].
The conventional shell FEs exhibit several, further numerical difficulties. These numerical problems arise from the very small value of one of the characteristic structural lengths, i.e., the shell becomes extremely thin. This phenomenon is known as a numerical locking effect. The standard, low-order, h-version shell FEs with Naghdi-type kinematics have the most serious locking phenomena such as the shear-and membrane locking, as well as the boundary layer effect [6,37,48]. Beside the use of the 2D shell FEs, the application of the solid-shell and degenerated FEs is also widespread. During the assembling procedure of these latter mentioned 3D shell FEs, the integrations with respect to the thickness coordinate are carried out numerically. However, neither the dimensionally reduced (2D) nor the 3D shell FEs are free from the locking effects; moreover, the development of the degenerated and solid-shell FEs can be accompanied by further numerical problems such as the curvature-, dilatational-, thickness-and trapezoidal locking; see the thorough study in [65].
To avoid these kinds of locking phenomena, the three main research directions are: (i) the mixed FEMs, (ii) the application of higher-order FEs and finite cell method (FCM), and (iii) the reduced integration accompanied by stabilization techniques. The numerical performance of the latter mentioned FEs is not satisfactory up to now, and additional techniques are required for the stabilization of their stiffness matrix [13,21,37]. The mixed interpolation of the strain tensor components (MITC) is another way to remove these types of lockings. The theoretical basis of the MITC FEs can be read in detail, for instance, in [17]. The numerical efficiency of these FEs was demonstrated by several test problems, among others, in [5,32]. Later, several additional improvements of the MITC shell FEs have been published, for example, in [28,66,67]. Besides, numerous h-version primal-mixed shell FEs and their improvements were elaborated, for example, in [22,68].
Some of these convergence problems can be avoided by the use of the p-version approximation technique [26,51,55,58]. However, a very complicated task is to establish automatically the FE mesh of a physical domain especially in the case of complex, real engineering structures. To circumvent this trouble, the FCM was developed, for example, in [44]. The FCM operates on an initial rectangular grid such that an adaptive integration technique is applied, which partitions the bounding cell into integration subcells during the mesh refinement. The FCM, or in other words fictitious domain method, is essentially an immersed boundary method combined with p-version FEs. A comprehensive overview about the beneficial behavior of the FCM can be found in [44].
An alternative strategy to circumvent the numerical difficulties appeared in shell FE modeling, is to use complementary strain energy-based, dual-mixed (DM) variational principles, which allow us to approximate directly the stress tensor as primary variable beside the displacement vector, thereby providing higher convergence rate and accuracy for the stresses. In contrast to the primal-mixed variational theorems, in the DM principles the divergence operator is shifted from the displacement vector to the stress tensor. Although a higher-order convergence can be achieved in the FE computation of stresses with the application of DM FEs than which is possible with the use of the primal-mixed ones, the simultaneous exact satisfaction of the symmetry condition for the stress tensor and C 0 continuity requirement on the stress vector causes additional difficulties during the development of DM FEs. A huge challenge is to construct these kinds of stable and effective FEs; see, for example, in [2,45,57].
One of the possibilities, to overcome the numerical difficulties mentioned in connection with the two-field DM FEs, is the fulfillment of the symmetry condition of the stress tensor in a weak form via the rotation tensor as Lagrangian multiplier. The structural FEs based on this approach can give more accurate results and faster convergence not only for the stress-but also for the displacement field. However, the development of these FEs is a much more difficult task as compared to primal and primal-mixed (displacement-based) formulations. The difficulties arise from the complexity of the theoretical background and the selection of the appropriate (stable) polynomial approximation space for the independent variables. The aforementioned procedure results in a three-field DM variational formulation; see the detailed descriptions, for instance, in [50,61,62]. This idea was originally suggested in [18]. The structural FEs based on this approach are elaborated on 2D and 3D problems of linear elasticity in [3,8,12,27,46,50,56] and on cylindrical shell problems in [59,63], as well as [29] using first-order stress functions for the a priori satisfaction of the translational equilibrium equation.
Up to now, each of the aforementioned DM variational principles uses symmetric functional settings. It means that the trial-and test function spaces are the same, yielding the Bubnov-Galerkin-type (BG) FE discretization schemes. The other strategy is to apply the Petrov-Galerkin-type (PG) FE method in order to overcome the locking effects. This methodology was first proposed in [7] and is based on variational principles in which all partial derivatives are shifted from the trial functions to the test functions. For this reason, these are called ultraweak variational principles as well. In contrast to the mixed variational principles, these methods can use non-symmetric functional settings, i.e., the trial and test basis functions can be different; therefore, the Petrov-Galerkin scheme might be needed for the FE discretization. Nevertheless, in the work of [7], the discontinuity property characterizes both trial and test functions, while in another version of PG, the word "discontinuous" refers only to the test functions. Practically, both types of PG methods use discontinuous or, in other words, broken test spaces allowing the approximation of the test functions at the element level, thereby introducing additional unknowns defined at all interelement boundaries. These FE discretization schemes are called discontinuous PG (DPG) methods. Detailed descriptions and several computational tests of the DPG method for the 2D, 3D, and axisymmetric shell problems of linear elasticity can be read in [16,20,39].
In accordance with the above particularized computational difficulties and research directions for their avoidance, the main aim of this paper is the construction of a locking-free and robust (uniformly stable) dual-and mixed axisymmetric hp-shell FE, which later can be extended to more complicated, coupled initialboundary value problems of elastic shells (see, for example, in [61,62]). Throughout the paper, the indicial notation of tensors and the usual summation convention is used. The range of the Latin indices appearing in the right subscript is 1, 2, 3, and that of the Greek indices is 1 and 2, except when special operations are introduced.
The text of the paper is divided into nine main Sections. After the overview of the relevant scientific background (Sect. 1), the applied three-field dual and mixed variational formulation is presented for the 3D problems of linear elasticity in Sect. 2. Our focus is turned to a comprehensive differential geometric description of thin axisymmetric shells in Sect. 3, applying NURBS interpolation to the generation of the shell mid-surface. In Sect. 4, a consistent dimensional reduction procedure is derived, and then, our discussion is proceeded with a systematic variable-number reduction process in Sect. 5, which is the final step in the derivation of the theoretical background. In Sects. 6 and 7, the resulting mathematical basis is followed by the development of the general DM axisymmetric hp-shell FE associated with appropriately selected polynomial stress-and displacement interpolations for the bending-shearing problems coupled with extension-compression. In Sect. 8, the computational performance of the demonstrated hp-shell FE is tested by a representative model problem for three different shell mid-surface geometries, focusing on the boundary layer effect in shell deformation. The paper closes with a brief summary and some comments and milestones on the outline of the future research in Sect. 9.

Boundary value problem of linear elasticity
Let V ⊂ R 3 denote a region of the 3D Euclidean space occupied by linearly elastic solid. Let V be bounded by the piecewise smooth surface S = ∂ V = S f ∪ S u with the outward normal vector n , where S f and S u are interpreted as two disjoint, open subsurfaces such that S f ∩ S u = Ø is valid. Hereby the closure of V can be given by V = V ∪ S. We assume that the linearly elastic body is exposed to the volumetric force b k in V and the traction force f k on the surface part S f , as well as the displacement vectorũ k is prescribed on the remaining subsurface S u .
In the above boundary value problem (BVP) of the linearly elastic solid, we seek the displacement vector u k , the infinitesimal strain-and rotation tensor k and ψ k , as well as the stress tensor σ k that satisfy the translational equilibrium equation in which the semicolon followed by an index in the subscript indicates covariant derivative with respect to a general curvilinear coordinate x , the rotational equilibrium equation as the symmetry condition for the stress tensor, where pqr is the covariant permutation tensor, as well as the geometric equation as the additive decomposition of the displacement gradient u k; , the symmetric-and the skew-symmetric part of which are defined, respectively, as where φ s denote the axial vector of ψ k and the constitutive relation (inverse Hooke's law) as well as the displacement BC as Dirichlet BC and the traction BC as Neumann BC. Here the fourth-order tensor C k rs is the elastic compliance tensor with the major symmetry C k rs = C rsk and the minor symmetries C k rs = C k sr = C krs . In the case of isotropic solids, the elastic compliance tensor can be expressed, for example, in terms of the Poisson ratio ν and the elasticity modulus E, as ( g kr g s + g r g ks ) − ν E g k g rs ; therefore, making use of Eq. (2), the constitutive equation (4) can be reduced to the form The complete system of Eqs. (1)-(3) and (4) or, for isotropic body, Eq. (8), along with the BCs (5)-(6), can be considered as the strong formulation for the BVP of linear elastostatics. However, the latter detailed BVP can be cast not only in differential (or strong) form but also in various variational forms. These variational principles restate original BVP of linear elasticity in integral forms. In these integral principles, the finite elements (FEs) and its extended versions can be applied as discretization methods to the approximation of unknown functions.

Three-field dual and mixed variational principle
In accordance with the aims formulated in the Introduction, the h-and p-version axisymmetric shell FEs developed in this paper are based on a three-field DM variational principle. In this variational theorem, the solution of the BVP of linear elasticity can be characterized as the unique stationary point of the complementary strain energy-based functional over the space of all vector fields u k , φ s and all a priori non-symmetric tensor fields σ k satisfying the constitutive relation (4), or for isotropic materials, Eq. (8), and the traction BC (6). After having taken the variation of functional (9) with respect to u k , φ s , and σ k , as well as setting them separately equal to zero and integrating by parts, we obtain the translational-and rotational equilibrium equations (1)-(2) and the geometric equation (3) as Euler-Lagrange equations and the displacement BC (5) as natural BC. Nevertheless, the traction BC (6) is the essential BC, which has to be enforced in the strong form.
The above three-field DM weak formulation of the BVP of linear elastostatics reads: find the trial functions u k , φ s , σ k satisfying Eq. (6) in such a way that for all corresponding test functions δu k , δφ s , δσ k satisfying the homogeneous form of Eq. (6). It is important to emphasize here that this DM variational principle is theoretically (and of course numerically) well suitable for the removal of the volumetric locking effect appearing in nearly incompressible linear elasticity. Namely, the strain tensor computed directly from the stress tensor by Eq. (8) tends to the deviatoric stress tensor as the Poisson ratio is near to one of its limit values, i.e., ν → 0.5, see [15]. It means that the aforementioned variational principle is always robust (uniformly stable) with respect to the Poisson ratio, i.e., volumetric locking-free.
Most of plate and shell theories which appeared in the scientific literature use a modified constitutive law arising from the stress hypothesis as additional assumption or, in other words, restriction. This procedure leads to the system of constitutive equations, which has a same structure as the plane elasticity constitutive equation. Thus, most of the classical shell theories are volumetric locking-free. However, the displacement-based 3D shell or solid FEs (which do not use the mentioned stress hypothesis) are not free from volumetric locking because they use the unmodified 3D Hooke's law.
In spite of the fact that the unmodified 3D constitutive equations are being used, the shell theory presented in this paper will be volumetric locking-free at the theory level because the inverse version (7) or (8) of 3D Hooke's law is applied, i.e., the coefficient 1/(1 − 2ν), which is responsible for the volumetric locking in the classical 3D FE formulations (when ν is near to 0.5), does not appear in the shell formulation developed in this paper. This is the reason why the DM-shell FE based on the variational formulation (10) will be volumetric locking-free or, in other words, any modification of the constitutive law is not needed (to obtain volumetric locking-free shell FE at the theory level).

Differential geometry of shells of revolution
Let a shell of revolution be considered now as a 3D solid body occupying a region of space bounded by the outer/inner axisymmetric surfaces S ± , which are symmetrically situated with respect to the mid-surface S 0 , and the lateral surfaces S 0 and S L , the intersections of which with S 0 are the circles 0 ⊂ S 0 and L ⊂ S L with radii R 0 and R L , respectively, namely the boundary of the shell domain V is given by S = S f ∪ S u = S ± ∪ S 0 ∪ S L (see the meridian section of the shell in Fig. 1).
The shell mid-surface S 0 can be established by the rotation of a curve about the x-axis. In other words, this curve is the meridian section of the mid-surface S 0 with any open plane lying on the x-axis (Fig. 2). The equation of the meridian curve is defined by the twice-differentiable functions where stands for the parameterization of the meridian coordinate ξ 1 and L is the length of the shell mid-surface S 0 in the axial direction x (Fig. 1). The position of any axisymmetric shell point is given by means of the meridian coordinate ξ 1 parameterized in Eq. (12) and the polar angle 0 ≤ ξ 2 < 2π, as well as the arc-length coordinate −d/2 ≤ ξ 3 ≤ +d/2 measured Fig. 1 Meridian section of the axisymmetric shell

Fig. 2
Axisymmetric shell mid-surface-coordinate systems along the normal to S 0 (Figs. 1, 2). The distance d between S + and S − in the direction ξ 3 is the thickness of the shell, which will be assumed to be constant ( Fig. 1). With the knowledge of Eq. (11), the relation between the Cartesian coordinate lines x, y, z and the curvilinear coordinates ξ α is given by: on the basis of which the position vector of any arbitrary mid-surface point P 0 ∈ S 0 can be expressed in the Cartesian coordinate frame as: ( Fig. 2). Introducing the notation A ξ 1 = (x ) 2 + (R ) 2 ≥ 1 and taking definition (11) and parameterization (12) into account, the total arc-lengths, along the meridian curve and the latitude circles lying in the planes perpendicular to the rotation axis x, can be determined, respectively, as where the prime appearing in the superscript means the differentiation with respect to the meridian coordinate ξ 1 . The equations, associated with the BVP of the axisymmetric shell by the related mathematical/mechanical quantities, will be written down in the curvilinear coordinate system attached to the shell mid-surface S 0 (Fig. 2). In order to carry out the indicated differentiations on the mid-surface, its curvature properties have to be determined exactly. For this reason, firstly, the covariant base vectors of the mid-surface coordinate frame are computed from partial differentiation of Eq. (13) in the following forms: Here a Latin index or a Greek index α preceded by a comma in the subscript indicates partial differentiation with respect to the corresponding coordinates ξ or ξ α . In view of the latter results, the covariant and contravariant metric tensors of the mid-surface coordinate system can be represented, respectively, by the matrices with the use of which and Eq. (14) the contravariant base vectors of the mid-surface coordinate frame are obtained as: Carrying out the partial differentiations of the bases a α with respect to the curvilinear coordinates ξ β and using Eqs. (14)- (16), the matrices of the covariant and mixed curvature tensor components, as well as the nonzero mid-surface Christoffel symbols of second kind, can be written, respectively, in the forms [11,60] and as well as For thin shells of revolution, the characteristic lengths measured along the mid-surface S 0 are much larger than the thickness d. In this case, the volume element, the mid-surface element, the outer/inner-and lateral surface elements can be given, respectively, as and where A 0 = A(ξ 1 = 0) and A L = A(ξ 1 = 1). In what follows, the meridian curve will be generated by a NURBS interpolation whose position vector r 0 (ξ 1 ) = x(ξ 1 ) e x + R(ξ 1 ) e R is defined by [47] in which i r = i x e x + i R e R are the position vectors of the control points lying in the plane x-R, e R is the radial unit vector given for an arbitrary but fixed angle coordinate ξ 2 , i w are the weights, n denotes the number of the control points, and i, p N (ξ 1 ) are the p-th degree B-spline basis functions. The control points 0 r and n r coincide with the endpoints of the meridian curve, i.e., with the points belonging to the coordinates ξ 1 0 = 0 and ξ 1 L = 1, respectively, see the definitions in Eqs. (11)- (12). The main advantage of the application of NURBS as meridian curve is its capability of the precise representation of either the conic sections and circular or free-form curves. Using the knot vector the B-spline basis functions are recursively defined by Introducing the rational basis functions Eq. (24) can be rewritten in the form The derivatives of the NURBS curve are obtained by computing the kth derivative of the B-spline basis functions (25)- (26), which can be given by A NURBS curve is infinitely many times continuously differentiable in the interior of a knot span and p − k times differentiable at a knot, where k is the multiplicity of a knot.

Dimensional reduction process
In order to derive a new axisymmetric shell model based on the three-field DM weak formulation (10), a systematic dimensional reduction will be used for the trial and test functions u k , φ s , σ k , and δu k , δφ s , δσ k . As a starting step, the 3D translational equilibrium equation (1) is written down in the mid-surface coordinate frame, which can be given for thin shells of revolution as where the 3D mid-surface covariant derivative σ k | can be put in the form which can be separated into [11,36,60] where and The above five equations and the homogeneous form of Eq. (27) are valid for the test function δσ k as well.
As a next point of the derivation the trial and test functions, as well as the known functions occurring in Eq. (10) are expanded into power series with respect to the thickness coordinate ξ 3 on the mid-surface S 0 : as well as Now all variables will be referred to unit mid-surface, covariant or contravariant base vectors, namely the physical components of the tensorial quantities, denoted by a hat over the variables, will be introduced as follows: in which there is no summation for the uppercase indices. Furthermore, our investigations will be restricted to axisymmetric problems of shells of revolution, i.e., the variables will not depend on the polar coordinate ξ 2 . Thus, upon substitution of the curvature tensor components (17)- (18), the nonzero mid-surface Christoffel symbols of second kind (19)-(21) and the Taylor series expansions (33) and (36) into the translational equilibrium equation (28) by the derivatives (29)-(32), expressing the results in terms of the physical components of tensorial variables (38) and making separation the obtained equations with respect to the powers of the thickness coordinate ξ 3 , we find from which it follows that the stress componentsσ kλ have to be approximated along the thickness coordinate ξ 3 by one degree lower thanσ k3 in order to satisfy any ith scalar expanded translational equilibrium equations. This is admissible because the weak formulation (10) allows the stress tensor to be treated as a priori nonsymmetric. In this paper, the following linear and quadratic stress approximations are applied: as well asb is valid for the components of the volumetric force density vector, namely the scalar, expanded translational equilibrium equations of i = 0 and i = 1 will be retained from Eqs.
As the final step in the dimensional reduction, the relating test functions δσ k and δû k , δφ s are approximated along ξ 3 , respectively, by the same degree polynomials as the trial functionsσ k andû k ,φ s : and According to Eqs. (42), (43) and (45), (46), the number of the unknown trial functions is 33 including the 21 stress coefficients 0σ k , 1σ k , 2σ k3 and the 6 displacement coefficients 0ûk , 0ûk , as well as the 6 rotation coefficients 0φ s , 1φ s , leading to a consistent shell model, see, for example, the consistency property of the displacement-based and stress-based shell theories, in [24,25,31,53,54]. In the next Section, a systematic process will be presented for the reduction in the number of variables.

Traction boundary conditions on the outer/inner surfaces
During the subsequent derivations, the displacement BC will be prescribed only on the lateral surface S u = S 0 , while surface tractions will be prescribed on the remaining surface part S f = S ± ∪ S L of the shell. The traction BCs (6) can be written on the outer/inner surfaces S ± in the form where (f ± ) k are prescribed surface tractions on S ± with outward unit normalsn ± defined byn ± 1 =n ± 2 = 0, n ± 3 = ±1. Upon substitution of the quadratic stress approximation (43) into the above BCs, we get 0σ considering that ξ 3 = ±d/2 on S ± . Introducing the vector-valued load coefficients 0f the combination of Eqs. (51) and (52) yields the transformed forms which can be viewed as field equations for the stress coefficients 0σ k3 , 1σ k3 , and 2σ k3 . As the first step of the variable-number reduction procedure, the number of the 21 stress variables is reduced by 6 by the a priori fulfillment of the constraint equations (53), which essentially leads to the elimination of the stress coefficients 1σ k3 and 2σ k3 from the weak formulation. Since on the basis of the variational theorem the homogeneous form of the BC (6) or the BCs (51)-(52) holds true for the test functions δ 0σ k3 , δ 1σ k3 and δ 2σ k3 , we can write 1σ expressing the symmetry of the stress tensor in an integral average sense along the thickness. By the use of these equations, not only all the rotation coefficients 0φ s , 1φ s , but also the additional six stress coefficients 0σ λ3 , 1σ κ3 , 0σ 12 , and 1σ 12 , are eliminated, as the final point in the variable-number reduction process. Accordingly, the number of the independent variables will be 15 including the 9 stress coefficients 0σ 3λ , 0σ 21 , 1σ 21 and 0σ 11 , 1σ 11 0σ 22 , 1σ 22 , 0σ 33 , as well as the 6 displacement coefficients 0ûk , 1ûk . From the variational formulation (10), it follows that the homogeneous forms of Eqs. (55)- (58) are valid for the test functions δ iσ 21 , δ iσ 12 , and δ iσ 3λ , δ iσ λ3 (i = 0, 1): Substitution of Eqs. (54) and (55) for i = 0 and i = 1, which are divided into two groups of equations: Eqs. (61), (63)- (64), and (66) expressed in terms of the stress coefficients 0σ 11 , 0σ 22 , 0σ 33 , 1σ 11 , 1σ 22 , and 0σ 31 describe the bending-shearing-and the extension-compression problems of the axisymmetric shell, while Eqs. (62) and (65) given in terms of the stress coefficients 0σ 21 , 0σ 32 , 1σ 21 describe its torsional problem. It can easily be proven for axisymmetric problems of shells of revolution that the bending-shearing equations (including tension-compression) are decoupled from the torsional equations. For this reason, our investigation is focused on the bending-shearing problems of axisymmetric shells.

Dual and mixed weak formulation for bending-shearing problems
In this Sections, a DM weak formulation is derived for axisymmetric, bending-shearing problems of a thin shell of revolution. After having vanished the last two integrals from the variational equation (10) by the elimination of all rotations by Eqs. (55)-(58), taking into account Eqs. (11), (12) and (22), (23), the integration with respect to the angle coordinate ξ 2 results in the variational form considering thatn 1 = −1,n 2 =n 3 = 0 on S 0 . Assuming now homogeneous and isotropic thin shells of revolution, inserting successively the elastic compliance tensor (7) expressed in terms of physical components and the linear and quadratic approximations (42), (43), (47), (48), as well as the transformed BCs (53), (54) and (55) (68) for the first integral appearing in the weak form (67). After having used the linear approximations (47) for the test functions δ 0σ k1 , δ 1σ k1 , taking into account Eq. (59) and carrying out the indicated integration, the second, boundary integral occurring in Eq. (67) can be written as: Substituting now the scalar translational equilibrium equations (61), (63), (64), (66) and their homogeneous forms for the corresponding test functions δ iσ 11 , δ iσ 22 , δ 0σ 33 , and δ 0σ 31 (i = 0, 1), as well as the linear approximations (45) and (49) for the trial-and test functions iû1 , iû3 , and δ iû1 , δ iû3 (i = 0, 1) into the last two integrals occurring in Eq. (67), we have Upon substitution of Eqs. (68)-(71) into the weak form (67), we get the final form of the DM variational equation serving as theoretical basis for the DM hp-FE formulation presented in the next Section. The dimensionally reduced shell model developed here for the axisymmetric, bending-shearing-and extension-compression problems of shells of revolution has 10 independent variables including the membrane stresses 0σ 11 , 0σ 22 , 0σ 33 , the bending stresses 1σ 11 , 1σ 22 , and the transverse shear stress 0σ 31 , as well as the mid-surface displacements 0û1 , 0û3 , and the displacements 1û1 , 1û3 being outside the mid-surface S 0 .
Since the traction BCs and their homogeneous forms on the outer/inner surfaces S ± were built in the DM weak form (67) by Eqs. (53), (54), the traction BC on the lateral surface S L has remained as the only subsidiary condition to the variational equation (67), which can be written aŝ σ k n × =σ k1n× 1 =σ k1 = (f × ) k on S L , considering thatn × 1 = 1,n × 2 =n × 3 = 0 are valid for the physical components of the normal vector to S L . Expanding now the prescribed lateral tractions (f × ) k into a truncated power series with respect to the thickness coordinate ξ 3 and using the linear stress approximations (42), we obtain the following traction BC on S L : which can be manipulated into the form 0σ k1 as essential BC, where 0 (f × ) k and 1 (f × ) k are the vector-valued load coefficients prescribed on the lateral surface S L . According to the variational theorem, the homogeneous form of the BC (6) or the BC (72) holds true for the test functions δ iσ k1 (i = 0, 1), i.e.,

hp-version axisymmetric shell finite element
In this Section, the h-and p-version DM, axisymmetric shell FEs will be presented and the relevant polynomial space will be given for the bending-shearing problems (including tension-compression) of homogeneous and isotropic shells of revolution. The most serious locking phenomena can accompany the convergence analyses of the classical shell FEs mainly in the case of the bending-shearing BVPs. Let us consider now the meridian sections of the axisymmetric shell mid-surface, as shown in Figs. 3, 4, and, 5. The master element ω st = {η | −1 ≤ η ≤ 1} is mapped onto the e-th physical element ω el = ξ 1 | ξ 1 e ≤ ξ 1 ≤ ξ 1 e+1 (of the meridian curve) with nodal points e and e + 1 by the transformation   Fig. 6 Two-element mesh Fig. 7 Sixteen-element mesh    Fig. 9 Convergence curves for the relative error in the maximum norm of the membrane stresses 0σ 11 and 0σ 22 : conical shell in which

Fig. 4 Fixed spherical shell with end loads
are the external shape functions, and ξ 1 e and ξ 1 e+1 are the meridian coordinates of the nodal points. Applying Eq. (75), the Jacobi determinant of the coordinate transformation (74) is given by: where e is the length of the e-th axisymmetric element ω el . The applied polynomial function space over the e-th FE is defined by the external shape functions N α for p = 1 and the bubble modes

Fig. 11
Convergence curves for the relative error in maximum norm of the transverse shear stress 0σ 31 and the membrane stress 0σ 33 : conical shell for p 2, where P j−1 (η) and P j−3 (η) are Legendre polynomials [14,58], as shown in Table 1. The DM FE space selected in this way leads to reliable numerical solutions at every p level (Sect. 8). The continuity requirements are imposed directly on the corresponding stress coefficients 0σ 11 , 1σ 11 , and 0σ 31 of the external shape functions introduced in (75).

Numerical example
This Section presents the computational performance of the DM hp axisymmetric shell FE model. The numerical computations were done by an FE code implemented in Octave environment. The hp-FE solutions are obtained directly for the stresses and the displacements.
To show the performance of the DM hp-shell FE model, convergence rates are measured in the relative error of energy norm and maximum norm of stresses and displacements, which are calculated, respectively, by where the indices R E F and F E imply the reference-and FE solutions, respectively. The reference is a solution calculated with high element number n and high polynomial degree p. The robustness and capability of the constructed DM element in the stress and displacement approximations are analyzed for both p-and hextensions. The error measures are computed for two different d/R 0 values (d/R 0 = 0.01 and d/R 0 = 0.0001, respectively) with elasticity modulus E = 2 · 10 5 MPa and Poisson's ratio ν = 0.3.
In our examples, we investigate three axisymmetric shells with different meridian curves, i.e., conical, spherical and hyperbolic mid-surface geometry. The lateral surface S 0 of the shells is fixed, i.e., the prescribed displacements 0ũ1 , 0ũ3 , 1ũ1 , and 1ũ3 at ξ 1 = 0 are set to zero, whereas the lateral surface S L of the shells is subjected to the transverse shear force    Taking into account the boundary layer effect, the meridian curves of the shells are divided into two subdomains with lengths s b and s 1 − s b , where s b = 3.67π/β. The two-element mesh shown in Fig. 6 is used for the p-approximations when the polynomial degree p is varying from 1 to 8 for all the two elements. During the h-extension, the subdomains s 1 − s b and s b as the initial FE mesh are refined uniformly and simultaneously in eight steps from the element number n = 2 to 16 (see the last step in Fig. 7). The polynomial degree of each element is set to p = 1 and kept fixed.

.1 MPa and the bending moment
In the case of conical-, spherical-and hyperboloid shells, the h-and p-convergence of the relative errors obtained for the energy norm computations is shown for moderately thick-and extremely thin shells (d/R 0 = 0.01 and d/R 0 = 0.0001) as the function of the number of degrees of freedom (DOF) on log-log scale in Figs. 8, 14, and 20, respectively. Then, the convergence curves of the relative errors calculated in maximum norm of the membrane stresses 0σ 11 , 0σ 22 , 0σ 33 , the bending stresses 1σ 11 , 1σ 22 , and the transverse shear stress 0σ 31 , as well as the in-plane-and the out-of-plane displacements, 0û1 , 0û3 , and 1û1 , 1û3 , are plotted against the number of DOF for the ratios d/R 0 = 0.01 and d/R 0 = 0.0001 in Figs. 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, and 21, 22, 23, 24, 25, respectively, in the case of conical-, spherical-and hyperboloid shells.
It can be observed that high convergence rates are experienced at both p-and h-approximations, which behave insensitively to the slenderness ratio d/R 0 . This excellent convergence behavior is true for both the

Fig. 19
Convergence curves for the relative error in maximum norm of the tangential displacements 1û1 and 1û3 : spherical shell stress-and the displacement variables, i.e., the membrane stresses 0σ 11 , 0σ 22 , 0σ 33 , the bending stresses 1σ 11 , 1σ 22 , and the transverse shear stress 0σ 31 , as well as the transverse deflection 0û3 , the through-the-thickness variation 1û3 and the tangential displacements 0û1 , 1û1 . As expected, the convergence histories indicate clearly that the new hp-type DM axisymmetric shell FE is free from the locking phenomena.

Summary
In this paper, a new hp-type, dual and mixed axisymmetric shell FE has been constructed for BVPs of linearly elastic shells of revolution. To summarize, the main basic assumptions and the most important restrictions on the shell theory presented are listed briefly as follows: (i) As usual in shell modeling, traction BCs are prescribed on the outer-and inner surfaces S ± ; however, on the lateral surfaces S 0 and S L , either traction-or displacement BC can be prescribed. (ii) For thin shells of revolution, the thickness is assumed to be sufficiently small, i.e., the relations d/L

Fig. 21
Convergence curves for the relative error in maximum norm of the membrane stresses 0σ 11 and 0σ 22 : hyperboloid shell stands for the least principal radius of curvature of the shell mid-surface S 0 . As a consequence, the shifter is equal to the identity tensor, i.e., μ k l ∼ = δ k l , see the definition of the shifter, for example, in [36,59,60]. (iii) Following from the latter assumption, the slenderness ratio d/R min is bounded. The upper bound is approximately 1/10 (for which both 1/15 and 1/25 are also chosen in the scientific literature [64]), while the lower bound belongs to the case when the shell thickness does not lie in the micro-range yet, i.e., d/R min is higher than approximately 10 −6 [4]. (iv) It is assumed that the tensorial quantities are multiple times (at least twice) differentiable functions of the thickness, namely, can be expanded into power series with respect to the thickness coordinate.
The derivation of the theoretical model is based on two key steps, using a three-field DM variational formulation which has the a priori non-symmetric stress tensor, the skew-symmetric, infinitesimal rotation tensor, and the displacement vector as independent field variables, as well as applying NURBS for the interpolation of the meridian curve. The first key step is the consistent dimensional reduction process and the second one is the variable-number reduction procedure. The former one has led to the linear approximation of the displacements, the infinitesimal rotations, and the stress componentsσ kλ , as well as the quadratic approximations of the stress componentsσ k3 along the thickness after truncating consistently the power series expansions of the corresponding equations. As the results of the latter one, the traction BC has been built in the DM variational formulation, and the symmetry condition of the stress tensor has been weakly imposed along the thickness by the elimination of the rotations. This long derivation has resulted in a new shell theory with 15 independent functions including 9 stress ( 0σ 11 , 0σ 22 , 0σ 33 , 1σ 11 , 1σ 22 , 0σ 31 , 0σ 12 , 1σ 12 , 0σ 32 ) and 6 displacement ( 0û1 , 0û2 , 0û3 , 1û1 , 1û2 , 1û3 ) components for the axisymmetric bending-shearing-(coupled with extension-compression) and torsional problems of shells of revolution. Finally, a new hp-version DM axisymmetric shell FE has been developed for the numerical solution.
The numerical efficiency of the developed DM hp-version axisymmetric shell FE has been tested by a representative bending-shearing (including tension-compression) mixed BVP for three different shell mid-surfaces: singly-and doubly curved axisymmetric shells with positive and negative Gaussian curvature, concentrating on the boundary layer effect during the shell deformation. Following from the computational experiments, the demonstrated hp-version shell FE provides robust (uniformly stable) convergence results not only for the stress but also for the displacement variables independently of the thickness value. It means that the constructed hp-version axisymmetric shell FE is free from the locking phenomena concerning the standard shell FEs.
Additionally, a nice feature of the proposed hp-shell FEM is that it can be used for the reliable computation of the deformation state of both extremely thin and moderately thick shells of revolution, obtaining a quadratic

Fig. 25
Convergence curves for the relative error in the maximum norm of the tangential displacements 1û1 and 1û3 : hyperboloid shell distribution for the transverse shear stress along the thickness. Furthermore, a beneficial property of the shell model is that either the through-the-thickness variation or the membrane stress normal to the shell mid-surface is not eliminated as independent unknown, and the derived first-and second-degree approximations for the other, basic variables along the thickness coordinate come from the demonstrated consistent dimensional-and variable-number reduction processes, namely these cannot be treated as any modification of existing shell theories.
Based on the computational experiments and the theoretical background, it is worth further developing the presented hp-version shell FE for (i) more general initial-BVPs (contact and dynamic problems) of shells [61,62], (ii) modeling the mechanical behavior of layered structures (composites) with extremely thin or moderately thick layers, taking into account the slip between them, (iii) automatic hp-adaptivity, applying the error indicators introduced, for example, in [8,20,39], (iv) coupled thermoelasticity and thermoelectroelasticity problems, using the code Artap for FE mesh optimization [41][42][43] and (v) 2D axisymmetric shell problems. This latter extension allows us to analyze more general BCs. These complex BVPs can be solved, for example, by Fourier series expansions of the related variational equations with respect to the polar angle coordinate, applying the DM axisymmetric shell FEs along the meridian, or using 2D DM-shell FEs, see a comparison in [40].
From both mathematical and engineering points of view, it would be worth testing the developed DM hp-shell FEs on the Girkmann benchmark problem being defined as the BVP of a spherical shell stiffened with a foot ring; see the detailed problem description in [38,58].
The presented hp-version axisymmetric shell FEs can be applied to the solution of BVPs occurring in the engineering practice, see, for example, the reservoir-pipe-valve system which is subjected to transient pressure induced by a water hammer [9].