Isogeometric analysis of thin Reissner–Mindlin shells: locking phenomena and B-bar method

Abstract

We propose a local type of B-bar formulation, addressing locking in degenerated Reissner–Mindlin shell formulation in the context of isogeometric analysis. Parasitic strain components are projected onto the physical space locally, i.e. at the element level, using a least-squares approach. The formulation allows the flexible utilization of basis functions of different orders as the projection bases. The introduced formulation is much cheaper computationally than the classical $\bar{B}$ method. We show the numerical consistency of the scheme through numerical examples, moreover they show that the proposed formulation alleviates locking and yields good accuracy even for slenderness ratios of $10^5$, and has the ability to capture deformations of thin shells using relatively coarse meshes. In addition it can be opined that the proposed method is less sensitive to locking with irregular meshes.

Appendices

Appendix A: Using an approximated normal vector field

In Eqs. 8 and 9, the normal vector field $\varvec{n}(\xi ,\eta )$ should be commonly built as $\sum R_A \varvec{n}_A$, where $\varvec{n}_A$ needs to be specified for each control point: it denotes the normal vector of a mid-surface point associated with the control point $\varvec{x}_A$. Although no approximation error occurs, it is time-consuming to solve for the mid-surface projecting points of $\varvec{x}_A$. In this paper the normal vectors at the Greville abscissae $(\tilde{\xi },\tilde{\eta })$ are adopted to construct a normal vector field approximately [21], the field is re-written as $\sum R_A(\xi ,\eta )\tilde{\varvec{n}}_A$. Following Eq. 2, these normal vectors $\tilde{\varvec{n}}_A$ at the Greville abscissae points can be easily calculated. Thus we have

$$\begin{aligned} \varvec{x}_P(\xi ,\eta ,\zeta )=\sum _{A=1}^{nm}R_A\left( \varvec{x}_A+\zeta \frac{h}{2} \tilde{\varvec{n}}_A\right) , \end{aligned}$$

(A.1)

and

$$\begin{aligned} \varvec{u}_P(\xi ,\eta ,\zeta )=\sum _{A=1}^{nm}R_A\left( \varvec{u}_A+\zeta \frac{h}{2} \varvec{\theta }_A \times \tilde{\varvec{n}}_A\right) . \end{aligned}$$

(A.2)

Appendix B: Detailed formulation of degenerated Reissner–Mindlin shell element

In the following we introduce some necessary aspects to obtain the element stiffness matrix:

• (a) Using Voigt notation, the relation between the strain vector and the stress vector is expressed as

  $$\begin{aligned} \varvec{\sigma } =\varvec{D}_g \pmb {\varepsilon }, \end{aligned}$$

  (B.1)

where $\varvec{D}_g$ is the global constitutive matrix, and

$$\begin{aligned} \varvec{D}_g =\varvec{T}^\mathrm{T} \varvec{D}_l \varvec{T}, \end{aligned}$$

(B.2)

here $\varvec{D}_l$ is a given local constitutive matrix. To fulfill the plane stress state $\sigma _{33}=0$, the local constitutive matrix is given by

$$\begin{aligned} \varvec{D}_l= \frac{E}{1-\nu ^2} \begin{bmatrix} 1 &{} \nu &{} 0 &{} 0 &{} 0 &{} 0 \\ \nu &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \frac{1-\nu }{2} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \kappa \frac{1-\nu }{2} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \kappa \frac{1-\nu }{2} \end{bmatrix}, \end{aligned}$$

(B.3)

in which $\kappa$ is the shear correction factor and it is set to be 5/6 in this research. To fit the physical shape of the mid-surface, a transformation is employed for $\varvec{D}_l$, specifically

$$\begin{aligned} \varvec{T}= \begin{bmatrix} t_{1x}t_{1x} &{} t_{1y}t_{1y} &{} t_{1z}t_{1z} &{} t_{1x}t_{1y} &{} t_{1y}t_{1z} &{} t_{1z}t_{1x} \\ t_{2x}t_{2x} &{} t_{2y}t_{2y} &{} t_{2z}t_{2z} &{} t_{2x}t_{2y} &{} t_{2y}t_{2z} &{} t_{2z}t_{2x} \\ n_{x}n_{x} &{} n_{y}n_{y} &{} n_{z}n_{z} &{} n_{x}n_{y} &{} n_{y}n_{z} &{} n_{z}n_{x} \\ 2t_{1x}t_{2x} &{} 2t_{1y}t_{2y} &{} 2t_{1z}t_{2z} &{} t_{1x}t_{2y}+t_{2x}t_{1y} &{} t_{1y}t_{2z}+t_{2y}t_{1z} &{} t_{1z}t_{2x}+t_{2z}t_{1x} \\ 2t_{2x}n_{x} &{} 2t_{2y}n_{y} &{} 2t_{2z}n_{z} &{} t_{2x}n_{y}+n_{x}t_{2y} &{} t_{2y}n_{z}+n_{y}t_{2z} &{} t_{2z}n_{x}+n_{z}t_{2x} \\ 2t_{1x}n_{x} &{} 2t_{1y}n_{y} &{} 2t_{1z}n_{z} &{} t_{1x}n_{y}+n_{x}t_{1y} &{} t_{1y}n_{z}+n_{y}t_{1z} &{} t_{1z}n_{x}+n_{z}t_{1x} \end{bmatrix}. \end{aligned}$$

(B.4)

At each Gauss quadrature point, vectors in Eqs. 2 and 3 are calculated, then the transformation in Eq. B.2 is performed.

• (b) Based on Eq. A.1, the Jacobi matrix is calculated by

$$\begin{aligned} \varvec{J}&= \begin{bmatrix} \frac{\partial \varvec{x}_P^x}{\partial \xi } &{} \frac{\partial \varvec{x}_P^y}{\partial \xi } &{} \frac{\partial \varvec{x}_P^z}{\partial \xi } \\ \frac{\partial \varvec{x}_P^x}{\partial \eta } &{} \frac{\partial \varvec{x}_P^y}{\partial \eta } &{} \frac{\partial \varvec{x}_P^z}{\partial \eta } \\ \frac{\partial \varvec{x}_P^x}{\partial \zeta } &{} \frac{\partial \varvec{x}_P^y}{\partial \zeta } &{} \frac{\partial \varvec{x}_P^z}{\partial \zeta } \end{bmatrix} =\begin{bmatrix} \sum \frac{\partial R_A}{\partial \xi } (\varvec{x}_A^x+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^x) &{} \sum \frac{\partial R_A}{\partial \xi } (\varvec{x}_A^y+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^y) &{} \sum \frac{\partial R_A}{\partial \xi } (\varvec{x}_A^z+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^z) \\ \sum \frac{\partial R_A}{\partial \eta } (\varvec{x}_A^x+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^x) &{} \sum \frac{\partial R_A}{\partial \eta } (\varvec{x}_A^y+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^y) &{} \sum \frac{\partial R_A}{\partial \eta } (\varvec{x}_A^z+\zeta \frac{h}{2} \tilde{\varvec{n}}_A^z) \\ \sum R_A (\varvec{x}_A^x+\frac{h}{2} \tilde{\varvec{n}}_A^x) &{} \sum R_A (\varvec{x}_A^y+\frac{h}{2} \tilde{\varvec{n}}_A^y) &{} \sum R_A (\varvec{x}_A^z+\frac{h}{2} \tilde{\varvec{n}}_A^z) \end{bmatrix}, \end{aligned}$$

(B.5)

then the physical gradient of the shape functions are

$$\begin{aligned} \begin{bmatrix} \frac{\partial R_A}{\partial x}\\ \frac{\partial R_A}{\partial y}\\ \frac{\partial R_A}{\partial z} \end{bmatrix} = \varvec{J}^{-1} \begin{bmatrix} \frac{\partial R_A}{\partial \xi }\\ \frac{\partial R_A}{\partial \eta }\\ \frac{\partial R_A}{\partial \zeta } \end{bmatrix} = \varvec{J}^{-1} \begin{bmatrix} \frac{\partial R_A}{\partial \xi }\\ \frac{\partial R_A}{\partial \eta }\\ 0 \end{bmatrix}. \end{aligned}$$

(B.6)

Moreover for $\zeta \frac{h}{2}R_A:=\tilde{R}_A$, their physical gradients are given by

$$\begin{aligned} \begin{bmatrix} \frac{\partial \tilde{R}_A}{\partial x}\\ \frac{\partial \tilde{R}_A}{\partial y}\\ \frac{\partial \tilde{R}_A}{\partial z} \end{bmatrix} = \varvec{J}^{-1} \begin{bmatrix} \frac{\partial \tilde{R}_A}{\partial \xi }\\ \frac{\partial \tilde{R}_A}{\partial \eta }\\ \frac{\partial \tilde{R}_A}{\partial \zeta } \end{bmatrix} = \varvec{J}^{-1} \begin{bmatrix} \zeta \frac{h}{2}\frac{\partial R_A}{\partial \xi }\\ \zeta \frac{h}{2}\frac{\partial R_A}{\partial \eta }\\ \frac{h}{2}R_A \end{bmatrix}. \end{aligned}$$

(B.7)

• (c) Taking a look at the strain formula Eq. 5, we form the strain matrix as

$$\begin{aligned} \varvec{B}= \begin{bmatrix} R_{A,x} &{} 0 &{} 0 &{} 0 &{} \tilde{R}_{A,x}\tilde{\varvec{n}}_A^z &{} -\tilde{R}_{A,x}\tilde{\varvec{n}}_A^y\\ 0 &{} R_{A,y} &{} 0 &{} -\tilde{R}_{A,y}\tilde{\varvec{n}}_A^z &{} 0 &{} \tilde{R}_{A,y}\tilde{\varvec{n}}_A^x\\ 0&{} 0 &{} R_{A,z} &{} \tilde{R}_{A,z}\tilde{\varvec{n}}_A^y &{} -\tilde{R}_{A,z}\tilde{\varvec{n}}_A^x &{} 0 \\ R_{A,y} &{} R_{A,x} &{} 0 &{} -\tilde{R}_{A,x}\tilde{\varvec{n}}_A^z &{} \tilde{R}_{A,y}\tilde{\varvec{n}}_A^z &{} \tilde{R}_{A,x}\tilde{\varvec{n}}_A^x - \tilde{R}_{A,y}\tilde{\varvec{n}}_A^y\\ 0&{} R_{A,z} &{} R_{A,y} &{} \tilde{R}_{A,y}\tilde{\varvec{n}}_A^y-\tilde{R}_{A,z}\tilde{\varvec{n}}_A^z &{} -\tilde{R}_{A,y}\tilde{\varvec{n}}_A^x &{} \tilde{R}_{A,z}\tilde{\varvec{n}}_A^x\\ R_{A,z} &{} 0 &{} R_{A,x} &{} \tilde{R}_{A,x}\tilde{\varvec{n}}_A^y &{} \tilde{R}_{A,z}\tilde{\varvec{n}}_A^z-\tilde{R}_{A,x}\tilde{\varvec{n}}_A^x &{} -\tilde{R}_{A,z}\tilde{\varvec{n}}_A^y \end{bmatrix}. \end{aligned}$$

(B.8)

Finally the element stiffness matrix is calculated by

$$\begin{aligned} \varvec{K}_e=\int _{\varOmega _e} \varvec{B}^\mathrm{T} \varvec{D}_g \varvec{B} \mathrm{d} \varOmega = \int _{-1}^{1}\int _{\eta _e}\int _{\xi _e} \varvec{B}^\mathrm{T} \varvec{D}_g \varvec{B} |\varvec{J}| \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta . \end{aligned}$$

(B.9)