The Criterion

In the criterion for detection of the dominance of bending strains, the solutions of two local problems from the equilibrated residual method are applied. The first solution corresponds to the following local discretization parameters: *H* = *h*, *P*1 = *p*, and *Q* = *q*. In the second problem, set to be free of the locking phenomena, one has: *H* = *h*, *P*2 = *pmax* = 8, and *Q* = *q*. The value of *pmax* = 8 is taken from the research of Pitkaranta [4] who demonstrated that for such a value the membrane locking disappears in shell structures regardless of the thickness. For this value also the shear locking in plate structures

disappears. It is well known that this type of locking is weaker than the membrane one in shells (compare remarks in [5,32], for example). The proposed criterion reads

$$\begin{cases} \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{U}} (\boldsymbol{\mu}^{\mathrm{HP}\_{1}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) + \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{S}} (\boldsymbol{\mu}^{\mathrm{HP}\_{1}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) & \geq \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{U}} (\boldsymbol{\mu}^{\mathrm{HP}\_{1}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) \\\ \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{U}} (\boldsymbol{\mu}^{\mathrm{HP}\_{2}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) + \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{S}} (\boldsymbol{\mu}^{\mathrm{HP}\_{2}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) & < \sum\_{c=1}^{E} \overset{\varepsilon}{\mathcal{U}}\_{\mathcal{U}} (\boldsymbol{\mu}^{\mathrm{HP}\_{2}\mathbf{Q}} \boldsymbol{\mathcal{Q}}) \end{cases} \tag{29}$$

The first condition detects dominance of the membrane and shear strains over the bending strains in the local problems corresponding to the global problem under consideration, while the second condition reflects the change of the dominance in the local problems. The form of the criterion corresponds to the shell structures, where coupling between shear, membrane and bending strains exists. In the case of plate structures, the coupling between the membrane strains and the shear and bending strains disappears and because of this the membrane strains have to be removed from the above criterion.

Fulfillment of the above criterion means that the true nature of the problem is bending-dominated and that the locking phenomenon appears in the assessed global problem. As a consequence, the locking has to be removed through increase of the approximation order *p*. The above criteria are necessary to confirm the results from the detection and assessment tools presented in the next sections, as those tools are not capable of distinguishing between the membrane- (and/or) shear-dominated problems and the bending-dominated ones.

#### 3.1.3. Sensitivity Analysis of the Local Solutions

The main purpose of the sensitivity analysis performed in this section is detection of the situation that the assessed global problem solution differs to the solution of the problem potentially free of locking to such an extent that the difference may sugges<sup>t</sup> the presence of the locking phenomena. Note that the theoretical ratio of the true bending-dominated solution and locked membrane-dominated (and/or shear-dominated) solution is *C* · *t* 2, with *C* being a constant. This ratio results from the proportionality of the membrane (and/or shear) part of the strain energy and the bending part of this energy to *t* and *t* 3, respectively, as shown in (8).

To compare the above two solutions, the following two local problems for each element of the plate or shell structure, or such parts of complex structures, has to be solved:

$$\stackrel{\varepsilon}{B}(\stackrel{\varepsilon}{\mathfrak{u}^{HP\_1Q}}, \stackrel{\varepsilon}{\mathfrak{v}^{HP\_1Q}}) - \stackrel{\varepsilon}{L}(\stackrel{\varepsilon}{\mathfrak{v}^{HP\_1Q}}) - \int\_{S\_{\varepsilon}\backslash\{P\cup\mathcal{W}\}} \stackrel{\varepsilon}{\langle \mathfrak{v}^{HP\_1Q}\rangle} \, ^T \langle \stackrel{\varepsilon}{r}(\stackrel{\varepsilon}{\mathfrak{u}^{hop}})\rangle \, dS\_{\varepsilon} = \, ^0 \tag{30}$$

where *P*1 = *p*, *H* = *h* and *Q* = *q* correspond to the global problem under consideration, and:

$$\stackrel{\text{e}}{B}(\stackrel{\text{e}}{\mathbf{u}}^{\text{H}}, \stackrel{\text{e}}{\mathbf{v}}^{\text{H}}, \stackrel{\text{e}}{\mathbf{v}}^{\text{H}\text{P}\_{2}\text{Q}}) - \stackrel{\text{e}}{\text{L}}(\stackrel{\text{e}}{\mathbf{v}}^{\text{H}\text{P}\_{2}\text{Q}}) - \int\_{S\_{\varepsilon}\backslash\langle\mathbb{P}\cup\mathcal{W}\rangle} \left(\stackrel{\text{e}}{\mathbf{v}}^{\text{H}\text{P}\_{2}\text{Q}}\right)^{T} \left<\stackrel{\text{e}}{\mathbf{v}}(\stackrel{\text{h}\text{p}\text{q}}{\mathbf{v}})\right> dS\_{\varepsilon} = \,\,\,\tag{31}$$

with *P*2 = *pmax* = 8, *H* = *h*, and *Q* = *q* corresponding to the problem potentially free of locking.

The following criterion is set to assess sensitivity of the solution to the change of the longitudinal approximation order from *P*1 to *P*2:

$$\sum\_{c=1}^{E} \frac{1}{2} \stackrel{\varepsilon}{B} (\stackrel{\varepsilon}{\mathfrak{u}^{H P\_1 Q}}, \stackrel{\varepsilon}{\mathfrak{u}^{H P\_1 Q}}) < a \sum\_{c=1}^{E} \frac{1}{2} \stackrel{\varepsilon}{B} (\stackrel{\varepsilon}{\mathfrak{u}^{H P\_2 Q}}, \stackrel{\varepsilon}{\mathfrak{u}^{H P\_2 Q}}) \tag{32}$$

where *a* is the coefficient which determines the ability of the adaptive method to overcome the locking phenomena with the standard *hp*-adaptive procedure. In the case of the applied procedure based on Texas three-step strategy [40] and error estimation based on the equilibrated residual method [43,44], the reasonable (verified numerically) value is suggested to be equal to *a* = 0.1. Fulfillment of the above criterion suggests the possibility of locking appearance. This possibility has to be confirmed by the criterion (29) of bending-dominance.

### 3.1.4. The Detection Algorithm

The algorithm of detection of the locking phenomena introduced into the standard three-step adaptive strategy proposed by Oden [40] is presented in Figure 1. The original part of the standard adaptation includes three steps: initial (*i* = 1), intermediate (*i* = 3) called also *h*-step, and final (*i* = 4) called also the *p*-step. Either local *h*-refinement or local *p*-enrichment is performed within each element of the latter two steps, based on the estimated approximation error values obtained from the equilibrated residual method of Ainsworth and Oden [43,44]. The modification of the initial mesh (*i* = 2) is performed as the second step when necessary. This additional step is composed of two stages: the detection of the phenomenon, and the changes in the mesh which consist in the global *p*-enrichment performed within all elements of the thin- or thick-walled structure or such a part of the complex structure.

The details of the additional step modifying the initial mesh are presented in Figure 2. At first, we search for the thick- or thin-walled parts of the generally complex structure. Such a structure can be composed of thick- and thin-walled parts, solid parts, and transition parts joining the previous two types of geometries. Of course, simple geometries composed of a single thin- or thick-walled part are also possible. Then, for each element of such a thin- or thick-walled part, the interelement stresses on the element boundary have to be calculated. Subsequently, for two local problems described with Equations (30) and (31), the terms of these equations are generated and the equations are solved. Next, one has to sum the element energies from (32) and compare two sums for two types of the local problems. Finally, if the criterion for locking detection is met, one has to confirm this result with the criterion for bending dominance detection (29). If it is fulfilled, then the modification of the initial mesh is performed, which is based on *p*-enrichment of all elements of the thin- or thick-walled part of the structure. Such elements of the modified structure are all equipped with the longitudinal order of approximation equal to *p* = *pmax* = 8. Next, the standard steps, *i* = 3 and *i* = 4, of *hp*-adaptation can be performed.

It should be underlined that calculation of the interelement stresses, performed in the above algorithm, is skipped in this paper as it is performed in the standard way presented in [43,44,46,54,55].

We also apply the standard *p*-enrichment algorithm performed at the modification stage of the additional adaptation step (*i* = 2). Such modification does not require any additional actions in comparison to the standard *p*-adaptation [41,42].

**Figure 1.** Three-step adaptive procedure extended with the fourth modification step.

**Figure 2.** Flow diagram of the overall algorithm for the locking detection.

#### *3.2. Calculation of the Optimized Value of p*

The idea presented in this section lies in making an attempt to assess the locking phenomena strength. This strength is measured with the *p* value necessary to remove the locking phenomenon. The criterion (32) assumes that the value of *p* = 8 is necessary for the removal.

This value corresponds to the coarsest mesh possible, composed of one rectangle (hexahedra) or two triangles (prisms). However, depending on the structure thickness *t* and the initial global discretization characterized by the value of *h*, the removal may require much lower (more optimal) value of the longitudinal order *p*. Moreover, in such circumstances, *p* = *pmax* = 8 leads to too rich discretizations. Note also that setting *p* = 8 in the modified mesh reduces further *hp*-adaptation to *h*-adaptation only as the maximum value of the longitudinal order of approximation *p* = *pmax* = 8 is applied in the modified mesh. We address both situations in numerical tests presented in the next sections of the paper.

## 3.2.1. The Idea

The idea lies in solution of the sequence of the local problems for each element *e* of a thin- or thick-walled structure or such a part of a complex structure. The sequence can be characterized with

$$\stackrel{\text{e}}{\text{B}}(\stackrel{\text{e}}{\text{u}}^{\text{HP}/\text{Q}}, \stackrel{\text{e}}{\text{v}}^{\text{HP}/\text{Q}}) - \stackrel{\text{e}}{\text{L}}(\stackrel{\text{e}}{\text{v}}^{\text{HP}/\text{Q}}) - \int\_{S\_{\text{\textell}} \backslash \langle P \cup \mathcal{W} \rangle} (\stackrel{\text{e}}{\text{v}}^{\text{HP}/\text{Q}})^{T} \left\langle \stackrel{\text{e}}{\text{r}} (\stackrel{\text{h}}{\text{r}}^{\text{hp}}) \right\rangle dS\_{\text{c}} = \,\tag{33}$$

where *Pj* = *P*1, ... , *Pmax* − 1, with *P*1 = *p* and *Pmax* = 8. The solutions of the above sequence are compared with the solution to the problem potentially free of locking phenomena, i.e.,

$$\stackrel{\text{def}}{\text{B}}(\stackrel{\text{e}}{\text{u}}^{\text{f}}, \stackrel{\text{f}}{\text{v}}^{\text{e}}, \stackrel{\text{e}}{\text{v}}^{\text{H}}{\text{P}}\_{2}, \stackrel{\text{f}}{\text{v}}) - \stackrel{\text{e}}{\text{L}}(\stackrel{\text{e}}{\text{v}}^{\text{H}}{\text{P}}\_{2}, \stackrel{\text{g}}{\text{Q}}) - \int\_{S\_{\text{c}}\backslash\{\text{P}\cup\mathcal{W}\}} \stackrel{\text{e}}{\text{v}}(\stackrel{\text{H}\prime p}{\text{v}}^{\text{Q}})^{T} \left< \stackrel{\text{e}}{\text{r}}(\stackrel{\text{h}\prime p\text{q}}{\text{s}}) \right> dS\_{\mathcal{E}} = \;=\;\mathbf{0} \tag{34}$$

where *P*2 = *Pmax* = 8. The comparisons are based on the summary criterion

$$\sum\_{\varepsilon=1}^{E} \frac{1}{2} \stackrel{\varepsilon}{B} (\stackrel{\varepsilon}{u}^{HP\_kQ}, \stackrel{\varepsilon}{u}^{HP\_kQ}) \ge a \sum\_{\varepsilon=1}^{E} \frac{1}{2} \stackrel{\varepsilon}{B} (\stackrel{\varepsilon}{u}^{HP\_2Q}, \stackrel{\varepsilon}{u}^{HP\_2Q}) \tag{35}$$

where, as above, one may set *a* = 0.1. For all particular orders of approximation *Pk*, fulfilling the above criterion, the adaptive algorithm can obtain the global adapted solution with the standard *hp*-adaptive algorithm. We choose the solution *Pl* fulfilling (35), such that:

$$\sum\_{\varepsilon=1}^{E} \frac{1}{2} \overset{\varepsilon}{\mathcal{B}} (\boldsymbol{\mu}^{HP\_{\mathbb{P}}Q}, \boldsymbol{\mu}^{HP\_{\mathbb{P}}Q}) = \min\_{P\_k} \left\{ \dots, \sum\_{\varepsilon=1}^{E} \frac{1}{2} \overset{\varepsilon}{\mathcal{B}} (\boldsymbol{\mu}^{HP\_kQ}, \boldsymbol{\mu}^{HP\_kQ}), \dots \right\} \tag{36}$$

The solution corresponding to *Pl* is enough and the most optimal for removal of the phenomenon. Note that for *Pl* = *P*1 removal of the locking is not necessary, as the standard *hp*-adaptation can be performed without changing the longitudinal order of approximation. This order is equal to its value from the assessed global problem, i.e., *P*1 = *p*. For *Pl* = *Pmax*, one obtains the same result as for the detection criterion (32). Note that the search for the optimized value of *p* should be performed if and only if the bending dominance criterion (29) is fulfilled.

## 3.2.2. The Algorithm

The structure of the algorithm remains the same as for the locking detection algorithm described in Section 3.1.4. In the adaptive algorithm presented in Figure 1, the only difference is that the modified mesh is generated with the optimized value of the longitudinal order of approximation within each element of the thin- or thick-walled part of the structure. This optimized value replaces the maximum value *pmax* = 8 applied in the case of the detection without optimization. In the algorithm presented in Figure 2, apart from the problem (31) or (34), the sequence of the local problems (33) is solved instead of the problem (30). After obtaining the solutions, the summary criteria (35) and (36) are checked, and the optimized value of the longitudinal order of approximation is established. If the optimized value is greater than the current value from the initial mesh, the confirmation from the criterion for the detection of bending dominance (29) has to be checked. If it is met, then the mesh is modified by increasing the longitudinal order in any element of the thick- or thin-walled part of the structure. Finally, the standard *hp*-adaptation (*i* = 3 and *i* = 4) can be performed.

#### **4. Verification and Utilization of the Proposed Tools**

Our numerical verification of the proposed numerical tools for detection and optimization of the longitudinal order of approximation is based on comparison of the results obtained from the detection and optimization tools and the results from the corresponding global solution of the model problems under consideration. The model problems concern all possible situations of locking existence or not, i.e., shear locking, membrane-shear locking, and the lack of locking. The performed comparisons should confirm that the detection and optimization based on the local, element solutions can replace the corresponding global analysis.

## *4.1. Model Problems*

It is well known from the literature presented in the state-of-the-art section that in the case of the bending-dominated plates, the shear locking may appear, depending on the plate thickness and the applied discretization. In the case of the bending-dominated shells, the shear–membrane locking is possible. On the other hand, in the case of the membrane-dominated shells, the locking phenomena do not appear. This is the reason for our choice of three model problems introduced below.

The first problem concerns a bending-dominated square plate. The plate longitudinal dimensions are equal to *l* = 3.1415 × 10−<sup>2</sup> m, while its basic thickness is *t* = 0.03 × 10−<sup>2</sup> m. A symmetric quarter of the plate can be seen in the Section 5.2.1. The plate is clamped along its edges. The vertical surface traction of the value equal to *p* = 4.0 × 10<sup>6</sup> N/m<sup>2</sup> is applied to the upper surface of the plate. The traction acts downwards. The plate, as well as the next two shell structures, is made of steel. For this material the Young's modulus is *E* = 2.1 × 10<sup>11</sup> N/m2, while the Poisson's ratio equals *ν* = 0.3.

The second model problem is a bending-dominated half-cylindrical shell. Its length is equal to *l* = 3.1415 × 10−<sup>2</sup> m, while its semi-circle circumference is *π R* ≈ 3.1415 × 10−<sup>2</sup> m, where *R* = 1.0 × 10−<sup>2</sup> m is the radius of the shell middle surface. The shell thickness is *t* = 0.03 × 10−<sup>2</sup> m. A symmetric quarter of the shell is displayed in Section 5.2.3. The shell is clamped along its straight edges, and the curved edges are free. The shell is loaded with the vertical traction of the value *p* = 4.0 × 10<sup>6</sup> N/m2, directed downwards. The third example is a symmetric half of a membrane-dominated cylindrical shell. The shell dimensions are analogous as for the bending-dominated shell, i.e., *l* = 3.1415 × 10−<sup>2</sup> m, *π R* ≈ 3.1415 × 10−<sup>2</sup> m, with *R* = 1.0 × 10−<sup>2</sup> m, and *t* = 0.03 × 10−<sup>2</sup> m. A symmetric octant of the shell is presented in Section 5.2.4. The symmetry boundary conditions are applied along the straight edges and one (left) curved edge of the octant, while there are no rotations along the other (right) curved edge. The shell is loaded with the internal pressure *p* = 4.0 × 10<sup>6</sup> N/m<sup>2</sup> acting outwards.

#### *4.2. Local Problems Solutions Versus Global Solutions*

In our tests, we applied the results from the algorithms for the detection of the bending dominance. We calculated sums of elemental strain energy components present in the criterion (29). In the case of the plate problem, the shear and bending energy components were computed, while, in the shell examples, the sum of shear and membrane components and the bending component of the energy were determined. The assessed global plate problem was characterized with *p* = 1, while in the shell examples the assessed global problems were characterized with *p* = 2. Then, the energy components were calculated for the sequence of problems (30) characterized with *P*1 = *p*, *p*+1, ... , 8. The averaged values of the interelement stresses were used in the local problems definition. The transverse order of approximation corresponded to the second-order hierarchical shell model (*Q* = *q* = 2) for all examples. The used mesh was characterized with *H* = *h* = *l*/2.

The above results from the sequence of the local problems were compared with the results from the global problems, where *q* = 2 and *h* = *l*/2. In the case of the plate, the following longitudinal orders of approximation were applied *p* = 1, 2, ... , 8 in the global problems, while in the case of two shell examples, the values of *p* = 2, 3, ... , 8 were taken into account. The results of the comparisons are presented in Figures 3–8 for the plate and bending- and membrane-dominated shells, respectively. The figures present the ratios of the shear (the plate) or the sum of shear and membrane (the shells) energies to the entire strain energy and the ratios of the bending energy to the entire strain energy obtained from the local problems (top) and the global problems (bottom). Values of these ratios are presented versus the longitudinal order of approximation.

Comparison of the top and bottom figures for the model problems leads to the conclusion that the detection tools based on solution of the element local problems give practically the same results as the solutions to the global problems, both qualitatively and quantitatively. In both bending-dominated problems, the change of the membrane dominance resulting from the locking phenomena to bending dominance due to removal of the locking is perfectly indicated. In addition, the membrane-dominated problem has been perfectly identified, as shown by the third couple of figures—the membrane energy component dominates over the bending one for all values of the longitudinal order of approximation.

**Figure 3.** The sum of local energies—the bending-dominated plate.

**Figure 4.** The global energy—the bending-dominated plate.

**Figure 6.** The global energy—the bending-dominated shell.

**Figure 7.** The sum of local energies—the membrane-dominated shell.

**Figure 8.** The global energy—the membrane-dominated shell.

#### **5. Effectivity of the Method in Model Problems**

In this section, we present examples of application of the algorithms for detection and/or assessment of the locking phenomena, and the standard mesh modification algorithms as well, in the adaptive analysis of model problems. The adaptive procedure is derived from the three-step strategy [40] composed of the global solution and error estimation for the initial, intermediate (or *h*-) and final (or *p*-) steps. The global problems are based on 3D-based hierarchical modelling and approximations for complex structures [45]. Our models and approximations are related, based or derived from [30,31,41,42], respectively. The error estimation is performed using the equilibrated residual method [43,54], adjusted to the 3D-based formulation of the hierarchical models and approximations [46]. The three-step strategy is completed with the modification step performed after the initial one. The algorithms from the previous sections are used within this additional step. The averaged stresses are used instead of the equilibrated ones in the local problems.

### *5.1. Problems and Methodology*

The three model problems in Section 4.1 are applied again. Two model bending-dominated plate examples are considered. In the first one, the data are exactly the same as in the mentioned section. In the second plate problem, we solve the plate of the same dimensions *l*/2 and *t* as for a quarter of the plate from the first plate. In the second problem, all four lateral sides of the plate are clamped. In the first problem, two of them correspond to symmetry boundary conditions, and only two sides are clamped. In the second problem, one deals with the most coarse mesh possible and the lowest longitudinal order of approximation possible. As a result, the strongest possible locking appears for the assumed plate thickness. The third example is the bending-dominated shell from Section 4.1. The only difference is the shell thickness which is now equal to *t* = 0.003 × 10−<sup>2</sup> m. The last example concerns the membrane-dominated shell. In this example, all data are exactly the same as in Section 4.1.

## *5.2. Numerical Examples*

In the analysis of the locking phenomena, the following data were treated as independent: the problem type (the type of the strain dominance), and the structure length *l* and thickness *t* resulting in the thinness ratio *l*/*t*. In the numerical analysis, also the initial mesh data, the relative value of the target admissible error *γT*, and the ratio *γI*/*γT* of the intermediate (after *h*-step) error to the target (after *p*-step) error were treated as independent quantities.

As results of the analysis, we present the meshes corresponding to three performed courses of adaptation: standard *hp*-adaptation and such adaptations preceded by the modification of the initial mesh with the increased longitudinal order of approximation *p*, equal either to the maximum or optimized value. The results are completed with the adaptive convergence curves corresponding to these three courses of adaptation. The convergence curves present the approximation error as a function of the number *N* of degrees of freedom (dofs). The absolute error is defined as a negative difference of the total strain energy *U* corresponding to the global solution under consideration and the reference energy *Ur* replacing the unknown exact value. The energies are calculated in accordance with the strain energy definition from Section 2. Due to the exponential character of *hp*-convergence, the curves are plotted as log(*Ur* − *U*) versus log*<sup>N</sup>*. The reference energy *Ur* is obtained numerically from calculations performed on over-killed meshes with the global discretization parameters equal to: *p* = 9, *q* = 2, *m* = 9, where *m* = *l*/2*h*. Apart from the absolute error values, the relative error values (*Ur* − *U*)/*Ur* are also presented and discussed.

#### 5.2.1. A Quarter of a Bending Dominated Plate
