Accurate and locking-free analysis of beams, plates and shells using solid elements

Abstract

This paper investigates the capacity of solid finite elements with independent interpolations for displacements and strains to address shear, membrane and volumetric locking in the analysis of beam, plate and shell structures. The performance of the proposed strain/displacement formulation is compared to the standard one through a set of eleven benchmark problems. In addition to the relative performance of both finite element formulations, the paper studies the effect of discretization and material characteristics. The first refers to different solid element typologies (hexahedra, prisms) and shapes (regular, skewed, warped configurations). The second refers to isotropic, orthotropic and layered materials, and nearly incompressible states. For the analysis of nearly incompressible cases, the B-bar method is employed in both standard and strain/displacement formulations. Numerical results show the enhanced accuracy of the proposed strain/displacement formulation in predicting stresses and displacements, as well as producing locking-free discrete solutions, which converge asymptotically to the corresponding continuous problems.

Appendices

Vectors and Matrices

Displacements $\varvec{u}$, strains $\varvec{\varepsilon }$, stresses $\varvec{\sigma }$ and forces $\varvec{f}$ are represented following Voigt’s notation as vectors

$$\begin{aligned} \varvec{u}&=\left( u_x,u_y,u_z\right) ^T \end{aligned}$$

(35)

$$\begin{aligned} \varvec{\varepsilon }&= \left( \varepsilon _x,\varepsilon _y, \varepsilon _z,\gamma _{xy},\gamma _{yz},\gamma _{xz}\right) ^T \end{aligned}$$

(36)

$$\begin{aligned} \varvec{\sigma }&= \left( \sigma _x,\sigma _y,\sigma _z,\tau _{xy}, \tau _{yz},\tau _{xz}\right) ^T \end{aligned}$$

(37)

$$\begin{aligned} \varvec{f}&= \left( f_x,f_y,f_z\right) ^T \end{aligned}$$

(38)

The differential symmetric gradient operator relating the displacements with the strains has the following form

$$\begin{aligned} \varvec{S}^T = \begin{bmatrix} \partial _x &{} 0 &{} 0 &{} \partial _y &{} 0 &{} \partial _z\\ 0 &{}\partial _y &{} 0 &{} \partial _x &{} \partial _z &{} 0 \\ 0 &{} 0 &{} \partial _z &{} 0 &{} \partial _y &{} \partial _x \\ \end{bmatrix} \end{aligned}$$

(39)

The projection matrix, introduced in Eq. (7), is

$$\begin{aligned} \bar{\varvec{G}}^T = \begin{bmatrix} n_x &{} \quad 0 &{} \quad 0 &{} \quad n_y &{} \quad 0 &{} n_z \\ 0 &{}\quad n_y &{} \quad 0 &{} \quad n_x &{} \quad n_z &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad n_z &{} \quad 0 &{} \quad n_y &{} \quad n_z \\ \end{bmatrix} \end{aligned}$$

(40)

where $\varvec{n}=\left( n_x,n_y, n_z\right) ^T$ is the outward normal vector at the boundary of the analysed domain $\varGamma _t$.

The discrete strain-displacement matrix (or discrete symmetric gradient operator) is expressed as

$$\begin{aligned} \varvec{B}_u = \left[ \varvec{B}_{u_{1}}, \dots , \varvec{B}_{u_{i}}, \dots , \varvec{B}_{u_{n}}\right] \end{aligned}$$

(41)

for $1 \le i \le n_n$, with $n_n$ being the number of nodes in the element. The submatrix $\varvec{B}_{u_{i}}$ and its volumetric part $\varvec{B}_{u_{i}}^{vol}$ are expressed in Voigt’s notation as

$$\begin{aligned} \varvec{B}_{u_{i}}= & {} \begin{bmatrix} \partial N_{i,1} &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad \partial N_{i,2} &{} \quad 0 \\ 0 &{} \quad 0 \quad &{} \quad \partial N_{i,3} \\ \partial N_{i,2} &{} \quad \partial N_{i,1} &{}\quad 0 \\ \partial N_{i,3} &{} \quad 0 &{} \quad \partial N_{i,1}\\ 0 &{} \quad \partial N_{i,3} &{} \quad \partial N_{i,2} \end{bmatrix} \end{aligned}$$

(42)

$$\begin{aligned} \varvec{B}_{u_{i}}^{vol}= & {} \frac{1}{3} \begin{bmatrix} \partial N_{i,1} &{} \quad \quad \partial N_{i,2} &{} \quad \partial N_{i,3} \\ \partial N_{i,1} &{} \quad \partial N_{i,2} &{} \quad \partial N_{i,3} \\ \partial N_{i,1} &{} \quad \partial N_{i,2} &{} \quad \partial N_{i,3} \\ 0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0 \end{bmatrix} \end{aligned}$$

(43)

where $N_i$ is the shape function of node i and $\partial N_{i,j}$ is its derivative with respect to the jth Cartesian coordinate ($j=[1:3]$). The deviatoric part is obtained by

$$\begin{aligned} \varvec{B}_{u_{i}}^{dev} = \varvec{B}_{u_{i}} - \varvec{B}_{u_{i}}^{vol}. \end{aligned}$$

(44)

Principle of virtual Work

This Appendix presents the derivation of equation (7) from equation (5) in two steps. First, Eq. (5) is premultiplied by an arbitrary virtual displacement $\delta \varvec{u}$ and integrated over the spatial domain $\varOmega$

$$\begin{aligned} \int _{\varOmega } \delta \varvec{u}^T \left[ \varvec{\mathcal {S}}^T \left( \varvec{D_s}\varvec{\varepsilon }\right) \right] \, d\varOmega + \int _{\varOmega } \delta \varvec{u}^T \varvec{f}\, d\varOmega = 0 \quad \forall \delta \varvec{u}\end{aligned}$$

(45)

Then, the Divergence Theorem is applied on the first term of the above equation yielding

$$\begin{aligned}&\int _{\varOmega } \delta \varvec{u}^T \left[ \varvec{\mathcal {S}}^T \left( \varvec{D_s}\varvec{\varepsilon }\right) \right] \, d\varOmega \nonumber \\&\quad =- \int _{\varOmega } \left( \varvec{\mathcal {S}}\delta \varvec{u}\right) ^T \, \left( \varvec{D_s}\varvec{\varepsilon }\right) d\varOmega + \int _{\varGamma } \delta \varvec{u}^T \left( \bar{\varvec{G}}^T\varvec{D_s}\varvec{\varepsilon }\right) \, d\varGamma \nonumber \\&\quad = - \int _{\varOmega } \left( \varvec{\mathcal {S}}\delta \varvec{u}\right) ^T \, \left( \varvec{D_s}\varvec{\varepsilon }\right) d\varOmega + \underbrace{\int _{\varGamma _u} \delta \varvec{u}^T \left( \bar{\varvec{G}}^T\varvec{D_s}\varvec{\varepsilon }\right) \, d\varGamma }_{=0} \nonumber \\&\qquad + \int _{\varGamma _t} \delta \varvec{u}^T \left( \bar{\varvec{G}}^T\varvec{D_s}\varvec{\varepsilon }\right) \, d\varGamma \nonumber \\&\quad = - \int _{\varOmega } \left( \varvec{\mathcal {S}}\delta \varvec{u}\right) ^T \, \left( \varvec{D_s}\varvec{\varepsilon }\right) d\varOmega + \int _{\varGamma _t} \delta \varvec{u}^T \left( \bar{\varvec{G}}^T\varvec{D_s}\varvec{\varepsilon }\right) \, d\varGamma \end{aligned}$$

(46)

In the previous derivation, Eq. (3) is used on the integral over $\varGamma$ and adopted the assumption that the prescribed displacements vanish on the boundary $\varGamma _u$. Finally, substituting Eq. (46) into Eq. (45) the final version of the Principle of Virtual Work for the mixed $\varepsilon /u$ formulation is obtained

$$\begin{aligned}&\int _{\varOmega } \left( \varvec{\mathcal {S}}\delta \varvec{u}\right) ^T \left( \varvec{D_s}\varvec{\varepsilon }\right) d\varOmega \nonumber \\&\quad = \int _{\varOmega } \delta \varvec{u}^T \varvec{f}\, d\varOmega + \int _{\varGamma _t} \delta \varvec{u}^T \left( \bar{\varvec{G}}^T\varvec{D_s}\varvec{\varepsilon }\right) \, d\varGamma \end{aligned}$$

(47)

presented in equation (7).

Variational multiscale stabilization method

This section presents the stabilization procedure leading to the final system of Eq. (24) of the $\varepsilon /u$ indepedent interpolation formulation. The stabilisation procedure adopted herein consists in the modification of the discrete variational form using the Orthogonal Subscales Method, introduced in [61] within the framework of the Variational Multiscale Stabilization methods [62, 63].

The stabilization of the problem is achieved by substituting the approximated strains in Eq. (9) with the following form

$$\begin{aligned} \varvec{\varepsilon }&\cong \varvec{\hat{\varvec{\varepsilon }}}+ \tau _\varepsilon \left( \tilde{\varepsilon } - \hat{\varepsilon } \right) \nonumber \\&= \varvec{N}_\varepsilon \varvec{E} + \tau _\varepsilon (\varvec{B}_u\varvec{U}-\varvec{N}_\varepsilon \varvec{E})\nonumber \\&= (1-\tau _\varepsilon )\varvec{N}_\varepsilon \varvec{E} + \tau _\varepsilon \varvec{B}_u \varvec{U} \end{aligned}$$

(48)

where $\tau _\varepsilon =\left[ 0,1\right]$ is a stabilization parameter. Observe that for $\tau _\varepsilon = 0$ the stabilization effect is lost, while for $\tau _\varepsilon = 1$ the strain interpolation of the standard irreducible formulation is recovered

$$\begin{aligned} \varvec{\varepsilon }\cong \tilde{\varepsilon } = \varvec{B}_u\varvec{U} \end{aligned}$$

(49)

The use of equation (48) in equations (6)-(7) gives the final stabilized set of equations for the mixed $\varepsilon /u$ FE formulation

$$\begin{aligned}&-(1-\tau _\varepsilon )\int _{\varOmega } \delta \varvec{E}^T \varvec{N}_\varepsilon ^T \varvec{D_s}\varvec{N}_\varepsilon \varvec{E} \, d\varOmega \nonumber \\&\quad + (1-\tau _\varepsilon ) \int _{\varOmega } \delta \varvec{E}^T \varvec{N}_\varepsilon ^T \varvec{D_s}\varvec{B}_u \varvec{U} \, d\varOmega = 0 \quad \forall \delta \varvec{E} \end{aligned}$$

(50)

$$\begin{aligned}&(1-\tau _\varepsilon )\int _{\varOmega } \delta \varvec{U}^T \varvec{B}_u^T \left( \varvec{D_s}\varvec{N}_\varepsilon \varvec{E}\right) d\varOmega \nonumber \\&\quad + \tau _\varepsilon \int _{\varOmega } \delta \varvec{U}^T \varvec{B}_u^T \varvec{D_s}\varvec{B}_u \varvec{U} d\varOmega = \int _{\varOmega } \delta \varvec{U}^T \varvec{N}_u^T \varvec{f}\, d\varOmega \nonumber \\&\quad + \int _{\varGamma _t} \delta \varvec{U}^T \varvec{N}_u^T \varvec{\bar{t}}\, d\varGamma \quad \forall \delta \varvec{U} \end{aligned}$$

(51)

Residual-based stabilisation procedures, like the one in (48) used herein, do not introduce any additional approximation nor any consistency error. For this, the stabilisation technique is variationally consistent, meaning that converging values of the unknowns $\varvec{\varepsilon }$ ad $\varvec{u}$ satisfying the Galerkin system (16)–(17) also satisfy the stabilized form (50)–(51). In particular, considering a converged solution, when the size of the element h tends to zero, $h \rightarrow 0$, $\varvec{\varepsilon }\rightarrow \varvec{N}_\varepsilon \varvec{E} = \varvec{B}_u \varvec{U}$ and the stabilization term vanishes. Considering a non-converged situation, the added terms $\tau _\varepsilon (\varvec{B}_u\varvec{U}-\varvec{N}_\varepsilon \varvec{E})$ are small, as they depend on the difference between two approximations of different order to the same quantity. This means that for a given FE mesh, using different values of the stabilization procedure yields slightly different results (see Appendix D). Nevertheless, the consistency of the residual-based stabilization guarantees that the discrete problem converges to the unique solution. The use of different stabilization parameters on the same mesh is analogous to the use of different FE interpolations of the same order of convergence with the same nodal arrangement.

As shown in [64, 65], the optimal convergence rate in linear problems is obtained reducing the stabilization on mesh refinement, such that

$$\begin{aligned} \tau _\varepsilon = c_\varepsilon \frac{h}{L_0} \end{aligned}$$

(52)

where $c_\varepsilon$ stands for a positive number of the order $c_\varepsilon =O(1)$, h for the finite element size and $L_0$ is the characteristic size of the problem.

Following the above, the stabilized system of equations becomes

$$\begin{aligned} \begin{bmatrix} -\varvec{M}_\tau &{} \varvec{G}_\tau \\ \varvec{G}_\tau ^T &{} \varvec{K}_\tau \end{bmatrix} \begin{bmatrix} \varvec{E}\\ \varvec{U} \end{bmatrix} = \begin{bmatrix} 0 \\ \varvec{F} \end{bmatrix} \end{aligned}$$

(53)

with

$$\begin{aligned} \varvec{M}_\tau&= (1-\tau _\varepsilon ) \varvec{M} \end{aligned}$$

(54)

$$\begin{aligned} \varvec{G}_\tau ^T&= (1-\tau _\varepsilon ) \varvec{G} \end{aligned}$$

(55)

$$\begin{aligned} \varvec{K}_\tau&= \tau _\varepsilon \underbrace{\int _\varOmega \varvec{B}_u^T \varvec{D_s}\varvec{B}_u d\varOmega }_{\varvec{K}} = \tau _\varepsilon \varvec{K} \end{aligned}$$

(56)

Influence of parameter $\tau _\varepsilon$

This Appendix investigates the influence of the stabilization parameter $\tau _\varepsilon$ in the numerical results obtained with the $\varepsilon /u$ FE formulation. The parameter $\tau _\varepsilon$ is defined in all the studied cases through the equation (52), in which intervenes the parameter c aside with the parameters h and $L_0$, associated with the finite element size and the characteristic size of the problem, respectively. Here, we investigate the influence of parameter c, with regard to the case of the clamped square plate.

Figure 47 presents the results obtained using three different values of $c = 5; 1; 1/5$ in equation (52). A value of $c=1$ corresponds to the reference value used for this case $\tau _{\varepsilon ,ref}=h/L_0$. The results show that the convergence rate is very similar for all the selected values of $\tau _\varepsilon$, as analytically predicted [64, 65]. The fact that using different values for c (i.e. different $\tau _\varepsilon$) produces different approximate solutions can be seen as similar to getting different approximate solutions by using meshes with different layouts, as already mentioned in Appendix C. Nevertheless, convergence to the solution, and optimal rate of convergence, are independent from the choice of parameter c.

As can be observed, for the same mesh, the use of a higher value of $\tau _\varepsilon$ results in an increase of the estimated error. This is to be expected, as for the limit value of $\tau _\varepsilon = 1$ the standard irreducible formulation is recovered. On the other end, very small values of $\tau _\varepsilon$ fail to effectively stabilize the $\varepsilon /u$ formulation.

Fig. 47

Clamped square plate: Local error in displacement $u_z$ at the center of the square plate versus the Number of Degrees of Freedom (DOF) for different values of $\tau _\varepsilon$

Table 5 Scordelis Lo-Roof: Vertical displacement at the midpoint of the free edge normalized with the reference solution 0.3086 given in [47]

Table 6 Hemispherical shell: Radial displacement at the load points normalized with the reference solution of 0.0940 given in [41]

Table 7 Pinched cylinder: Vertical displacement at the midpoint normalized with the reference solution of $0.18248 \cdot 10^{-4}$ given in [53]

Comparison with solid-shell and EAS FEs

This Appendix presents a comparison between the numerical results of the standard displacement based linear hexahedron (referred in the tables as Q1), the proposed $\varepsilon /u$ FEs (referred in the tables as Q1Q1) and the reported results of several successful solid-shell and EAS elements for three benchmark shell problems: the Scordelis-Lo roof (Table 5), the hemispherical shell (Table 6) and the pinched cylinder (Table 7).

The following solid-shell and EAS elements are considered:

• Wriggers and Koralc QS/E9 [13]: 3D solid-shell enhanced strain element with 9 enhanced modes based on Taylor expansion with exact symbolic integration.

• Wriggers and Koralc QS/E12 [13]: 3D solid-shell enhanced strain element with 12 enhanced modes based on Taylor expansion with exact symbolic integration.

• Reese [51]: EAS solid-shell based on reduced integration with hourglass stabilization (QISPs).

• Kim et al. [21]: ANS solid-shell with plane stress assumption (XSolid85).

• Alves de Sousa et al. [66] : EAS solid-shell with reduced (in-plane) integration (RESS).

• Areias et al. [25]: EAS solid element with penalty stabilization.

• Kasper and Taylor [20]: Mixed-enhanced strain element with nine enhanced modes (H1/ME9).

• Schwarze and Reese [22]: Reduced integration solid-shell based on the EAS and the ANS concepts.

• Huang et al. [23]: unsymmetric 8-node hexahedral solid-shell (US-ATFHS8).

• Sze et al. [16]: hybrid stress ANS solid-shell.

It is observed that:

• 1 The standard general purpose FEs lock in the tested curved thick shell situations, while the proposed $\varepsilon /u$ FEs and the solid-shell elements do not.

• 2 The special purpose solid-shell elements, enhanced with higher order bending modes, are more accurate than the general purpose $\varepsilon /u$ finite elements. However, their corresponding displacement convergence rate is the same.

• 3 Even if the stable solid-shell elements are notoriously more accurate than the corresponding underlying linear element, the asymptotic rate of convergence of displacements is the same as they do not interpolate with the full second order polynomial needed to achieve higher order convergence.

• 4 Only displacement results are reported in the literature for the solid-shell elements. The mixed $\varepsilon /u$ FEs are devised to yield enhanced strain and stress order of convergence.

• 5 All the reported tests are performed on hexahedral elements, as this is the shape of the solid-shell elements. Mixed $\varepsilon /u$ FEs can be equally shaped as prisms or tetras, without loss of convergence rate.

• 6 All the reported tests are performed in regular meshes. EAS elements often underperform in distored meshes. Huang et al. [23] solve this quaint at the expense of using an unsymmetrical element.