**2. Preliminaries**

Two issues are addressed in this section. The first one deals with the model problems considered in this research. The second one is presentation of the nature of the locking phenomena.

## *2.1. Model Problems*

Let us consider a wide range of problems of linear elasticity covered by the standard local (strong) formulation

$$\begin{aligned} \sigma\_{\stackrel{ij}{j}}^{ij} + f^i &= 0 \\ \sigma^{ij} &= D^{ijkl} \varepsilon\_{kl} \\ \varepsilon\_{kl} &= \frac{1}{2} (\boldsymbol{u}\_{k,l} + \boldsymbol{u}\_{l,k}) \end{aligned} \quad \text{or} \quad \mathbf{x} \in V \tag{1}$$

composed of the equilibrium, constitutive, and geometrical equations, respectively. Above, *i*, *j*, *k*, *l* = 1, 2, 3, while the smooth components *f i* ∈ *L*<sup>2</sup>(*V*) stand for the known body loading vector, with *V* representing volume of the body. The tensor components *σij* and *εkl* stand for stresses and strains, while the unknown vector components *uk* denote displacements. The terms *Dijkl* are components of the fourth-order tensor of elasticities. The elasticities are symmetric *Dijkl* = *Djikl* = *Dijlk* = *Dklij* and satisfy the strong and uniform ellipticity condition: *<sup>D</sup>ijklξijξkl* ≥ *αξijξkl*, *ξij* = *ξji*, where *ξij* is any real-valued tensor of the second rank and *α* is a positive constant. Note that such elasticities may correspond to three-dimensional stress and strain states or to the constrained, plane stress or plane strain, states. Note also that, for the specific case of isotropy, the elasticities can be expressed by the Young's modulus *E* and Poisson's ratio *ν*.

The set (1) has to be completed with the displacement and stress boundary conditions of the form

$$\begin{aligned} \sigma^{ij} u\_{\rangle} &= \ \ p^{i}, \quad \mathbf{x} \in \mathbb{S} \mathbf{p} \\ u\_{i} &= \ \ w\_{i}, \quad \mathbf{x} \in \mathbb{S}\_{W} \end{aligned} \tag{2}$$

where *i*, *j* = 1, 2, 3, while the smooth components *pi* ∈ *L*<sup>2</sup>(*S*) and the terms *wi* stand for the known values of stress vector components and displacement vector components on the parts *SP* and *SW* of the surface *S* ≡ *∂V* of the body *V*, with *S* = *SP* ∪ *SW* and *SP* ∩ *SW* = ∅. The vector components *nj* represent unit outward normal to the surface part *SP* = *ST* ∪ *SB* composed of the top *ST* and bottom *SB* surfaces of the body. The way they are determined is explained below.

Equation (1) is valid within the body volume *V*. Such a volume may represent thick- or thin-walled structures or complex structures containing such parts. Due to the model character of the considerations of this section, we limit ourselves to the first case. The corresponding thin-walled domain of the volume *V* is sufficiently smooth (Lipschitzian or smoother), open, and bounded region. As we apply the 3D Cartesian description of the problem, i.e., we employ Cartesian coordinates *x*, the following explicit (curvilinear) and implicit (Cartesian) definitions of the volume *V* are appropriate:

$$\begin{aligned} \mathcal{V} &= \left\{ \boldsymbol{\eta} \in \mathbb{R}^3 \colon (\boldsymbol{\eta}^1, \boldsymbol{\eta}^2) \in \mathbb{S}\_{M\prime}, \boldsymbol{\eta}^3 \in (-\operatorname{t}/2, \operatorname{t}/2) \right\} \\ &\equiv \left\{ \mathbf{x}(\boldsymbol{\eta}) \in \mathbb{R}^3 \colon \mathbf{x} = \boldsymbol{\mathcal{F}}(\boldsymbol{\eta}), \ (\boldsymbol{\eta}^1, \boldsymbol{\eta}^2) \in \mathbb{S}\_{M\prime}, \boldsymbol{\eta}^3 \in (-\operatorname{t}/2, \operatorname{t}/2) \right\} \end{aligned} \tag{3}$$

where *F* is the reversible map converting the curvilinear coordinates *η* into the Cartesian ones *x*. The function *t* = *<sup>t</sup>*(*η*1, *η*<sup>2</sup>) measures the symmetric thickness in the direction *η*3 normal to the mid-surface *SM*. The mid-surface *SM*, where *η*3 = 0, of the thin-walled body can be any sufficiently smooth (Lipschitzian or smoother), two-dimensional, open and bounded region. Formal definition of the thin-walled body geometry also needs introduction of the lateral part *SL* of the body boundary *S* as

well as the upper (top) *ST* and lower (bottom) *SB* parts of this boundary. More details on the applied definition of the thin-walled geometry can be found in [48].

The considered model problem, described by the local (strong) formulation (1), completed with the boundary conditions (2), can also be presented in the weak variational form:

$$B(\mathfrak{u}, \mathfrak{v}) = L(\mathfrak{v}) \tag{4}$$

where the admissible displacement vector is *v* ∈ V, while the corresponding space reads V = {*v* ∈ (*H*<sup>1</sup>(*V*))<sup>3</sup> : *v* = **0** on *SW*}. The solution function for displacements is *u* ∈ *w* + V, with *w* standing for the lift of the Dirichlet data (see [41]). This lift is consistent with the second Equation (2). The bilinear form *<sup>B</sup>*(*<sup>v</sup>*, *u*) and linear form *<sup>L</sup>*(*v*) from the above functional represent the virtual strain energy and the virtual work of the external body and surface loadings, *f* and *p*, respectively, i.e.,

$$\begin{aligned} \label{eq:SDAC-1} B(\mathfrak{u}, \mathfrak{v}) &= \int\_{V} \mathfrak{e}^{T}(\mathfrak{v}) \mathsf{D} \, \mathfrak{e}(\mathfrak{u}) \, dV \\ L(\mathfrak{v}) &= \int\_{V} \mathfrak{v}^{T} f \, dV + \int\_{S\_{\mathfrak{P}}} \mathfrak{v}^{T} \mathfrak{p} \, dS \end{aligned} \tag{5}$$

Above, *D* stands for the matrix representation of the elasticity constant tensor present in the second Equation (1), while *ε* is the six-component vectorial representation of the strain tensor defined with the third Equation (1).

Let us introduce now the finite element approximation of the variational functional (4). The approximation is based on the general rules of *hp* approximations [41], applied to 3D-based hierarchical shell models of order *q* presented in [38,45]. The approximated variational functional reads

$$B(\mathfrak{u}^{hpq}, \mathfrak{v}^{hpq}) = L(\mathfrak{v}^{hpq}) \tag{6}$$

where the approximated admissible displacements *vhpq* ∈ V*hpq* belong to the space V*hpq* = {*vhpq* ∈ (*H*<sup>1</sup>(*V*))<sup>3</sup> : *vhpq* = **0** on *SW*}. Additionally, *uhpq* ∈ *whpq* + V*hpq*, with the approximated values *whpq* of the lift *w*.

Note that the approximated variational formulation (6) can also be expressed in the language of finite elements

$$
\mathbb{K}q^{hpq} = \mathbb{F} \tag{7}
$$

where *K* and *F* are the global stiffness matrix and the global forces vector, while *qhpq* stands for the displacement dofs vector corresponding to the solution field *uhpq* of displacements.

## *2.2. Locking Phenomena*

It was demonstrated in a numerous works (see [49,50], for example) that the three-dimensional elasticity description of the strain energy *U* ≡ 12*B*(*<sup>u</sup>*, *u*) of the thin-walled body *V* tends to the following limit value when the thickness of the body tends to zero, *t* → 0:

$$\begin{split} \mathcal{U} \mathcal{U} &= \frac{1}{2} \int\_{V} \sigma^{ij} \varepsilon\_{ij} \, dV \equiv \frac{1}{2} \int\_{V} \left[ \sigma^{a\beta} \varepsilon\_{a\beta} + \sigma^{3} \varepsilon\_{33} + 2 \sigma^{3\beta} \varepsilon\_{3\beta} \right] dV \\ &\to \frac{1}{2} t \int\_{S\_{M}} \frac{E}{1 - \nu^{2}} [(1 - \nu) \gamma\_{a\beta} \gamma\_{a\beta} + \nu \gamma\_{a a} \gamma\_{\beta\beta}] \, dS\_{M} + \frac{1}{2} t \int\_{S\_{M}} \frac{E}{1 + \nu} (\gamma\_{a\beta} \gamma\_{a\beta} + \gamma\_{\beta\beta} \gamma\_{\beta\beta}) \, dS\_{M} \\ &\quad + \frac{1}{2} t^{3} \int\_{S\_{M}} \frac{E}{12 (1 - \nu^{2})} [(1 - \nu) \kappa\_{a\beta} \kappa\_{a\beta} + \nu \kappa\_{a a} \kappa\_{\beta\beta}] \, dS\_{M} \\ &=: \quad \mathcal{U}\_{W} + \mathcal{U}\_{\mathcal{S}} + \mathcal{U}\_{\mathcal{b}} \end{split} \tag{8}$$

Above, the three-dimensional components of the stress and strain tensors, *σij* and *<sup>ε</sup>ij*, *i*, *j* = 1, 2, 3, of the three-dimensional theory of elasticity can be expressed in the thin limit with the two-dimensional longitudinal strain components of the mid-surface *SM*, *γαβ*, *α*, *β* = 1, 2; transverse strain components, *γα*3 and *γ*3*β*, *α*, *β* = 1, 2; and the components of the tensor of curvature variation, *καβ*, *α*, *β* = 1, 2. Additionally, in the thin limit, the three-dimensional isotropic elasticity constants can be replaced with the isotropic plane stress constants of the first-order shell theory, expressed by the Young's modulus *E* and the Poisson's ratio *ν* and resulting from the plane stress assumption *σ*33 = 0. In addition, the kinematic condition of no elongation of the normals to the mid-surface, *γ*33 = 0, comes into play. The thin limit strain energy can then be approximated with the limit values of the membrane, shear, and bending parts, *Um*, *Us*, *Ub* , of the strain energy *U*. The quantity *k* is the shear strain correction factor, equal to 5/6 and resulting from the applied first-order model which leads to false (constant) transverse-shear stresses.

Depending on the geometry, loading and kinematic boundary conditions, one may distinguish between the bending-dominated problems, when

$$
\Box I\_b \gg \Box I\_m + \Box I\_s \tag{9}
$$

and membrane-dominated or shear–membrane-dominated problems, where

*Um Ub*, *Us* ∼= 0 *Um* + *Us Ub*(10)

respectively. Then, it results from (9) that, in the analytical solution of the bending-dominated problems, the shear and membrane strains should be equal or very close to zero, i.e.,

$$
\gamma\_{13} = \gamma\_{31} \cong 0, \quad \gamma\_{23} = \gamma\_{32} \cong 0 \tag{11}
$$

for the plate and shell structures, and

$$
\gamma\_{11} \cong 0, \quad \gamma\_{12} \cong 0, \quad \gamma\_{22} \cong 0 \tag{12}
$$

for the shell structures where the coupling of the membrane and bending strains exists. Note that the decoupling of the membrane and bending strains shown by the limit expression of the relation (8) exists

only in the thin limit (*t* → 0) for the first-order shells, while for the first-order plates it holds for any thickness *t*.

The shear and shear–membrane lockings are purely numerical phenomena which happen for poor discretizations and result from insufficiently accurate approximation of the analytical constraints (11) and/or (12). If such inaccurate approximation occurs then the shear and/or membrane strains are too large and the shear and membrane strain energies grow. These energies dominate over the bending strain energy, leading to numerical over-stiffening of the structure, which in turn gives too small (locked) values of displacements of the numerical solution. Remembering that *Uhpq* ≡ 12*B*(*uhpq*, *uhpq*) ≡ 12 (*qhpq*)*<sup>T</sup>Kqhpq*, one has

$$\mathbf{K}q^{hpq} = \mathbf{F}, \quad \mathbf{F}\_I = \text{const}, \; \mathbf{K}\_{I\downarrow} \uparrow \Rightarrow q\_{I}^{hpq} \downarrow \tag{13}$$

where *KI J*, *FI* are components of the global stiffness matrix and the global forces vector, with *I*, *J* = 1, 2, ... , *N* and *N* being the global number of degrees of freedom (dofs). The components *qhpq J* represent the global vector of unknown displacement dofs *qhpq*corresponding to the approximated field of displacements *<sup>u</sup>hpq*.

#### **3. Locking Detection, Assessment, and Resolution**

The proposed ideas of locking detection, assessment, and resolution are based on application of three existing numerical techniques. The first one is the equilibrated residual method of error estimation [43,44], here applied to solution of the local problems of two types. The second idea is sensitivity analysis based on comparison of the solutions to two local problems which differ with the longitudinal orders of approximation. The results of such a comparison can be used for the locking detection or assessment. If the locking is detected, then the mesh can be modified through the increase of the longitudinal order of approximation, by means of the standard *p*-enrichment technique [41,42].

The mentioned three techniques can be combined with the three-step *hp*-adaptive procedure controlled by the approximation error [51,52] or the automatic, iterative *hp*-adaptation driven by the interpolation error [41,42,53].

#### *3.1. The Idea and Algorithm of a Posteriori Phenomenon Detection*

The idea of a posteriori detection of the numerical locking was originally proposed by Zboi ´nski [38]. Here, we develop and verify this idea. It consists in solution of two local problems for each element of the thick- or thin-walled part of the structure. The solutions to these problems are obtained from the equilibrated residual method. These two problems differ with the longitudinal order of approximation *p*. In the first problem, it is equal to its current value from the global problem under consideration. In the second local problem, its value corresponds to the problem potentially free of locking. For two mentioned solutions, the strain energy norms are calculated and compared. This corresponds to sensitivity analysis where the sensitivity of the solutions of the local problems to the change of the longitudinal order of approximation is assessed. Once the locking is detected, one may modify the mesh by increasing the current order of the longitudinal approximation to its value corresponding to the problems free of locking. The similar approach is used also for determination of the optimized value of the longitudinal approximation order which corresponds to the minimum value sufficient for removal of the phenomenon.

#### 3.1.1. Solutions from the Equilibrated Residual Method

Let us recall now the equilibrated residual method of a posteriori error estimation, invented by Ainsworth and Oden [43]. All relations of this section are taken from the cited work. Here, we apply

this method to the estimation of the numerical solutions of the 3D-based hierarchical shell problems [46]. The principal relation of the method is [44]

$$\begin{aligned} B(\mathfrak{e}, \mathfrak{v}) &= \ \mathcal{B}(\mathfrak{u}^{HPQ} - \mathfrak{u}^{hpq}, \mathfrak{v}) = \mathcal{B}(\mathfrak{u}^{HPQ}, \mathfrak{v}) - \mathcal{B}(\mathfrak{u}^{hpq}, \mathfrak{v}) \\ &= \ \mathcal{L}(\mathfrak{v}) - \mathcal{B}(\mathfrak{u}^{hpq}, \mathfrak{v}), \ \forall \mathfrak{v} \in \mathcal{V}^{HPQ}(V) \end{aligned} \tag{14}$$

Above, *uHPQ* stands for either the exact solution of the problem or sufficiently accurate approximation of such a solution. The quantity *e* is the error vector corresponding to the assessed numerical solution *<sup>u</sup>hpq*. The kinematically admissible displacements *v* ≡ *vHPQ* belong to the space V *HPQ* defined in analogy to the spaces introduced for (4) and (6). The discretization parameters *H*, *P*, and *Q* are the counterparts of *h*, *p*, and *q* from the relation (6). In the case of the exact solution (*uHPQ* ≡ *u*), one deals with 1/*H*, *P*, *Q* → <sup>∞</sup>, while in the case of the approximation, one deals with the finite values of the discretization parameters *H*, *P* and *Q*.

After noticing that *e* = *uHPQ* − *uhpq* and introduction of this definition into the error functional (14), the latter simplifies to:

$$B(\mathfrak{u}^{HPQ}, \mathfrak{v}) \;=\; L(\mathfrak{v}), \; \forall \mathfrak{v} \in \mathcal{V}^{HPQ}(\mathcal{V}) \tag{15}$$

This way one searches for the approximation of the exact solution (displacements) of the problem instead of the problem error itself. However, once the displacements are determined, the error can be calculated from the above error definition. The presented approach is applied in the works [30,46], for example.

Let us now follow the work [43] and divide the domain *V* into *E* subdomains *Ve*, corresponding to finite elements *e* = 1, 2, ... , *E*. The global functional (15) and its left- and right-hand side components can now be presented as a sum of the element contributions:

$$\begin{array}{rcl}B(\mathfrak{u}^{HPQ},\mathfrak{v})&=&\sum\_{\mathfrak{e}=1}^{\tilde{E}}\tilde{B}(\mathfrak{u}^{HPQ},\stackrel{\mathfrak{e}}{\mathfrak{v}})\\&\quad\quad\quadL(\mathfrak{v})&=&\sum\_{\mathfrak{e}=1}^{\tilde{E}}\tilde{L}(\tilde{\mathfrak{v}})=\sum\_{\mathfrak{e}=1}^{\tilde{E}}\tilde{L}(\tilde{\mathfrak{v}})-\sum\_{\mathfrak{e}=1}^{\tilde{E}}\tilde{\beta}(\tilde{\mathfrak{v}})\\&\quad&=&\sum\_{\mathfrak{e}=1}^{\tilde{E}}\left[\tilde{L}(\tilde{\mathfrak{v}})+\int\_{S\_{\mathfrak{r}}\backslash\{S\_{\mathfrak{p}}\sqcup S\_{\mathfrak{N}}\}}^{\varepsilon}\tilde{\mathfrak{v}}^{\prime}(\tilde{\mathfrak{r}}^{\mathrm{hpq}})\right]dS\_{\mathfrak{E}}\end{array}\tag{16}$$

where the element admissible displacements are defined as the global admissible displacements projected onto elements *Ve*: *ev* ≡ *<sup>v</sup>*|*Ve*. As demonstrated by Ainsworth and Oden [43], the sum of the auxiliary element functionals *e β*( *ev*) is equal to 0, as the internal load consistency condition must hold.

Above, the vectors *er*(*uhpq*) represent the interelement loading due to equilibrated stresses. These stresses are defined in [43,44], for example. Their definition is as follows

$$<\langle \overset{\varepsilon}{r}(\overset{hps}{u}^{hpq}) \rangle \ = \overset{fc\,\varepsilon}{\underset{a}{\overset{\varepsilon}{r}}} \overset{\varepsilon}{r}(\overset{hps}{u}^{hpq}) + \overset{cf\,f}{\overset{cf\,f}{r}} (\overset{hps}{u}^{hpq}) \tag{17}$$

where

$$\begin{cases} \overset{\varepsilon}{r}(\boldsymbol{u}^{hpq}) = \mathbf{H}(\overset{\varepsilon}{\nu}) \overset{\varepsilon}{\sigma}(\boldsymbol{u}^{hpq}),\\ \overset{\varepsilon}{r}(\boldsymbol{u}^{hpq}) = \mathbf{H}(\overset{\varepsilon}{\nu}) \overset{\varepsilon}{\sigma}(\boldsymbol{u}^{hpq}), \end{cases} \tag{18}$$

and

$$H(\stackrel{\mathcal{C}}{\dot{\boldsymbol{\nu}}}) = \begin{bmatrix} \boldsymbol{\nu}\_1 & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\nu}\_2 & \boldsymbol{0} & \boldsymbol{\nu}\_3 \\ \boldsymbol{0} & \boldsymbol{\nu}\_2 & \boldsymbol{0} & \boldsymbol{\nu}\_1 & \boldsymbol{\nu}\_3 & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\nu}\_3 & \boldsymbol{0} & \boldsymbol{\nu}\_2 & \boldsymbol{\nu}\_1 \end{bmatrix} \tag{19}$$

The components of the unit normal vector *eν* = [*<sup>ν</sup>*1, *ν*2, *<sup>ν</sup>*3]*<sup>T</sup>* are defined on the element surface *Se*. The vectorial six-component representations *eσ* and *fσ* of the stress tensors have to be determined for the element *e* and any of its neighbours *f* . The terms *f eα* and *er* stand for the splitting function diagonal matrices and stress vectors. These matrices are: *f eα* = diag[*<sup>α</sup>*1, *α*2, *<sup>α</sup>*3], with *f eα* = **1** − *e fα*, **1** = diag[1, 1, 1] and directional components *αm*, *m* = 1, 2, 3. The algorithms for calculation of the splitting functions in the case of the internally unconstrained and constrained (e.g., by the Reissner-Mindlin kinematic constraints) are described in [43,44,46], respectively. The alternative is replacement of the equilibrated stresses with their averaged counterparts: *f eα* = *e fα* = diag[ 12 , 12 , 12 ].

Comparing relations (15) and (16) one can notice that the minimization of the global energy functional can be replaced with the minimization of local (element) functionals (see [43,44,46]). The element functionals read

$$\stackrel{\varepsilon}{u}\stackrel{HPQ}{\rightsquigarrow}\in\stackrel{\varepsilon}{\mathcal{V}}^{HPQ}(V\_{\varepsilon})\colon\quad\stackrel{\varepsilon}{\mathcal{B}}(\stackrel{\varepsilon}{u}^{HPQ},\stackrel{\varepsilon}{\mathfrak{v}})=\stackrel{\varepsilon}{L}(\stackrel{\varepsilon}{\mathfrak{v}})+\int\_{S\_{\varepsilon}\backslash\backslash(P\cup Q)}\stackrel{\varepsilon}{\mathfrak{v}}^{T}\langle\stackrel{\varepsilon}{r}(\stackrel{hpq}{u})\rangle\,dS\_{\varepsilon},\quad\forall\ \stackrel{\varepsilon}{\mathfrak{v}}\in\stackrel{\varepsilon}{\mathcal{V}}^{HPQ}(V\_{\varepsilon})\tag{20}$$

where *e uHPQ* is the element solution function and *e* V *HPQ* denotes the space of kinematically admissible element displacements *ev* ≡ *e <sup>v</sup>HPQ*. As shown in [43,44], the solution to the local problems (20) exists and is either unique or unique up to rigid body motions, for Dirichlet and Neumann boundary value problems, respectively. The mixed boundary value problems are also possible. Note that one deals with the Dirichlet problems for clamped elements adjacent to the external boundary of the structure, and with the Neumann problems for elements not adjacent to this boundary, namely for elements from the interior of the structure.

#### 3.1.2. Check on Bending-Dominance of the Solution

Check on bending-dominance of the solution to the global problems of the type (6) was proposed by Zboi ´nski [32]. Such a check is necessary as the locking phenomena are present only in the bending-dominated problems. In the membrane-dominated problems the phenomena do not appear. In addition, in the case of the mixed dominance, the locking may appear.

## Strain Energy Components

In this work, we apply the proposed approach to the local solutions from the equilibrated residual method and to their sum. To apply this approach, one has to take into consideration that the strain energy approximation from (8) corresponds to the structures of the infinitely small thickness, for which the first-order theory can approximate the three-dimensional description. As we intend to check the locking phenomena in the structures of the finite thickness, described with the hierarchical models corresponding to higher-order shell theories, the decomposition of the three-dimensional strain energy into component energies typical for the first-order models can only be done approximately. Let us start with the following decomposition of the strain vector into its longitudinal, transverse and shear parts

$$\mathfrak{e}\_{\mathfrak{e}} = \begin{bmatrix} \mathfrak{e}\_{l}, \mathfrak{e}\_{\mathfrak{u}}, \mathfrak{e}\_{\mathfrak{e}} \end{bmatrix}^{\mathrm{T}} \tag{21}$$

where *εl* = [... ,*εαβ*, ... ]*T*, *εn* = [*<sup>ε</sup>*33]*<sup>T</sup>*, *εs* = [... ,*ε*3*<sup>α</sup>*,*εα*3, ... ]*T* and *α*, *β* = 1, 2 are the local directions tangent to longitudinal directions of the thin-walled structure, and the index 3 corresponds to the local transverse direction. The analogous decomposition of the stress vector *σ* is also possible. Such decompositions lead to the division of the density *υ* of the total strain energy into the following components:

$$\begin{aligned} \upsilon &= \ \sigma^{ij}\varepsilon\_{ij} \\ &= \ \sigma^{a\beta}\varepsilon\_{a\beta} + \sigma^{33}\varepsilon\_{33} + 2\sigma^{3a}\varepsilon\_{3a} \\ &= \ \upsilon\_l + \upsilon\_u + \upsilon\_s \end{aligned} \tag{22}$$

where *i*, *j* = 1, 2, 3 stand for the global Cartesian directions. In the above relation, the symmetries: *σ*3*<sup>α</sup>* = *σ<sup>α</sup>*<sup>3</sup> and *ε*3*α* = *εα*3 were taken into account. The terms *υl*, *υn*, *υs* represent longitudinal, transverse, and shear components, respectively, of the density function. Consequently, the strain energy can also be divided into the corresponding energy components *Ul*, *Un*, *Us*:

$$\begin{aligned} \vert \mathcal{U} \rangle &= \int\_{V} \upsilon \, dV \\ &= \vert \mathcal{U}\_{\text{I}} + \mathcal{U}\_{\text{I}} + \mathcal{U}\_{\text{s}} \end{aligned} \tag{23}$$

Our calculations of the strain components of (21) on the element level correspond to the 3D-based hierarchical shell formulation proposed and developed in [38,45]. In this formulation local strains (defined in two longitudinal and the third transverse directions) and the global displacement dofs are applied. Because of that, the relation defining local strains *ε* (or their components) by the product of the element matrix*e B*itsandtheelementdofsvector*e qhpq* be

strain-displacement (or components) displacement can written as follows *e*

$$\begin{aligned} \varepsilon\_{\varepsilon} &= \stackrel{\varepsilon}{B} \stackrel{\mu\_{\text{I}}}{q}^{\text{h}pq} \\ &= \begin{bmatrix} \stackrel{\varepsilon}{B}\_{\text{I}} & \stackrel{\varepsilon}{B}\_{\text{II}} \stackrel{\varepsilon}{B}\_{\text{S}} \end{bmatrix} \stackrel{\varepsilon}{q}^{\text{h}pq} \\ &= \begin{bmatrix} \stackrel{\varepsilon}{B}\_{\text{I}} \stackrel{\varepsilon}{q}^{\text{h}pq} \stackrel{\varepsilon}{,B}\_{\text{I}} \stackrel{\varepsilon}{q}^{\text{h}pq} \stackrel{\varepsilon}{,B}\_{\text{S}} \stackrel{\varepsilon}{q}^{\text{h}pq} \end{bmatrix} \\ &= \begin{bmatrix} \varepsilon\_{\text{I}} \,\,\varepsilon\_{\text{II}} \,\,\varepsilon\_{\text{S}} \end{bmatrix}^{\text{T}} \end{aligned} \tag{24}$$

where the longitudinal *e Bl*, transverse *e B<sup>n</sup>*, and shear *e Bs* blocks of the strain-displacement matrix include the latter matrix rows corresponding to the in-plane, normal and shear strain components: *εl*, *εn*, *εs*. Note that the vector of displacement dofs is related to the displacement field *e uhpq* with the standard interpolation formula of the element *e*: *e uhpq* = *eN e<sup>q</sup>hpq*, where *eN* stands for the shape function matrix of the element.

Because of the 3D-based shell models applied in our research, further decomposition of the strain energy will be performed approximately. The exact decomposition is easy for the conventional shell models with the degrees of freedom defined on the mid-surface. Then, the dofs corresponding to the even and odd powers of the thickness contribute to the membrane and bending strains, respectively. In the 3D-based formulation, the dofs are defined along (or through) the thickness, and extracting their membrane and bending contributions needs much more complex procedure explained in [32]. The approximate procedure requires replacement of the top and bottom displacements of the symmetric thin-walled structure with

their sums and differences which contribute to the mid-surface displacements and rotations. This allows for distinguishing the following components of the in-plane strains

$$
\varepsilon\_{a\beta} \quad \approx \quad \gamma\_{a\beta} \leftarrow \kappa\_{a\beta} \tag{25}
$$

where *γαβ* stand for the mid-surface strains (membrane strains) and *καβ* are the three-dimensional counterparts of the change of curvature tensor (bending strains). With the above division, the in-plane density *υl* of strain energy can be decomposed into three parts (membrane, bending, and coupling ones):

$$
\boldsymbol{\upsilon}\_{l} \approx \boldsymbol{\upsilon}\_{m} + \boldsymbol{\upsilon}\_{b} + \boldsymbol{\upsilon}\_{c} \tag{26}
$$

that contain products of the membrane strains, bending strains, and their combination, respectively. Consequently, one also has

$$
\mathcal{U}\_I \quad \approx \quad \mathcal{U}\_m + \mathcal{U}\_b + \mathcal{U}\_c \tag{27}
$$

Calculation of the two components of the right-hand side of (25) on the element level needs adequate transformation of the numerical representation of element longitudinal strains *εl*, defined as a matrix product of the corresponding part of the strain-displacement matrix and the vector of displacement dofs. For the 3D-based hierarchical shell formulation, presented in [38,45], this transformation reads

$$\begin{aligned} \varepsilon\_{l} &= \prescript{c}{}{\mathcal{B}}\_{l} \prescript{f}{}{\eta}^{h\eta} \\ &= \ [\mathcal{B}\_{l}, \mathcal{B}\_{l}, \mathcal{B}\_{o}] \times [\boldsymbol{q}\_{l}, \boldsymbol{q}\_{b'}, \boldsymbol{q}\_{o}]^{\intercal} \\ &= \ [\mathcal{B}\_{l} + \mathcal{B}\_{b'}, \mathcal{B}\_{l} - \mathcal{B}\_{b'}, \mathcal{B}\_{o}] \times [\frac{1}{2}(\boldsymbol{q}\_{l} + \boldsymbol{q}\_{b}), \frac{1}{2}(\boldsymbol{q}\_{l} - \boldsymbol{q}\_{b}), \boldsymbol{q}\_{o}]^{\intercal} \\ &= \ [\mathcal{B}\_{ls}, \mathcal{B}\_{\boldsymbol{\theta}}, \mathcal{B}\_{o}] \times [\boldsymbol{q}\_{s'}, \boldsymbol{q}\_{\boldsymbol{q}'}, \boldsymbol{q}\_{o}]^{\intercal} = \mathcal{B}\_{s} \boldsymbol{q}\_{s} + \mathcal{B}\_{\boldsymbol{\theta}} \boldsymbol{q}\_{\boldsymbol{q}} + \mathcal{B}\_{o} \boldsymbol{q}\_{o} \\ &= \ \gamma + \kappa + r \\ &\approx \ \gamma + \kappa \end{aligned} \tag{28}$$

where the sub-blocks *Bt*, *Bb*, and *Bo* correspond to the top *q<sup>t</sup>*, bottom *qb*, and all other *qo* displacement dofs of the element dofs vector *e <sup>q</sup>hpq*. Location of the dofs of three types is on the top and bottom surfaces of the shell element, and apart from these two surfaces, respectively. In turn, the sub-blocks *B<sup>s</sup>*, *<sup>B</sup>φ* correspond to the mid-surface displacement dofs *qs* = 1 2 (*qt*+*qb*) and rotational dofs *qφ* = 1 2 (*qt*<sup>−</sup>*qb*). The approximate character of the performed transformation results from neglecting the mentioned other dofs in the resultant decomposition into the membrane and bending strain contributions.
