Shell elements

Strain Calculation

The first step to calculating stresses is to calculate the strains in the shell element. All the computations are performed in the natural coordinate system of the shell ($r_1$, $r_2$ and $r_3$). The geometry of the element at any time t is defined as:

$$^t x_i = \sum_{k=1}^4 h_k ^t {x_i}^k + \frac{r_3}{2} \sum_{k=1}^4 a_k h_k ^tV_{ni}^k$$

Therefore, the incremental displacements of any point within the shell element with natural coordinate $r_i$ at time $t$ in the global cartesian coordinate system ($^ox_1$, $^ox_2$ and $^ox_3$) are:

$$u_i = \sum_{k=1}^4 h_k u_i^k + \frac{r_3}{2} \sum_{k=1}^4 a_k h_k (-^tV_{2i}^k \alpha_k + ^tV_{1i}^k \beta_k)$$

If $^t \mathbf{g_i} = \partial ^t \mathbf{x}/\partial r_i$ are the covariant base vectors, then the Green-Lagrange strain components can be written as:

$$\tilde{\epsilon}_{ij} = \frac{1}{2}(^t \mathbf{g_i} \cdot ^t \mathbf{g_j} - ^0 \mathbf{g_i} \cdot ^0 \mathbf{g_j})$$

The expressions for the small and large strains in the 11, 22, 12, 13 and 23 directions are provided in equations 21 to 24 of Dvorkin and Bathe (1984). The expressions for the shear strains in the 13 and 23 directions also include a correction for shear locking. In MOOSE the small strain increments are computed in ADComputeIncrementalShellStrain and the small + large strain increments are computed in ADComputeFiniteShellStrain). Note that strain in the 33 direction is not computed.

Apart from the strain increment, two other matrices (B and BNL) are computed by the strain class. These two matrices connect the nodal displacements/rotations to the small and large strains, respectively, i.e., $e=B u$ and $\eta = B_{NL} u$. These two matrices are then used to compute the residuals in ADStressDivergenceShell.

Elasticity tensor

The shell element assumes a plane stress condition, i.e., the normal stress in the 33 direction should be zero. To enable the plane stress condition, an edited constitutive tensor is created in the local cartesian coordinate system ($\hat{C}_{mnop}$). This local constitutive tensor is then transformed into contravariant constitutive tensor using:

$$\tilde{C}_{ijkl} = (g^i \cdot \hat{e}_m)(g^j \cdot \hat{e}_n)(g^k \cdot \hat{e}_o)(g^l \cdot \hat{e}_p) \hat{C}_{mnop}$$

where $g^i$ are the contravariant base vectors. The contravariant base vectors $g^1$, $g^2$ and $g^3$ are determined from the covariant base vectors $g_1$, $g_2$ and $g_3$ using the following equations: $$g^1=\dfrac{g_2 \times g_3}{g_1 \cdot (g_2 \times g_3)}$$

$$g^2=\dfrac{g_3 \times g_1}{g_1 \cdot (g_2 \times g_3)}$$

$$g^3=\dfrac{g_1 \times g_2}{g_1 \cdot (g_2 \times g_3)}$$ These contravariant base vectors satisfy the orthogonality condition: $g^i \cdot g_j =\delta_{ij}$, where $\delta_{ij}$ is Kronecker delta. The $\hat{e}_i$ vectors are the local cartesian coordinate system computed using the covariant base vectors, i.e.,

$$\hat{e}_3 = \frac{g_3}{|g_3|}; \hat{e}_1 = \frac{g_2 \times \hat{e}^3}{|g_2 \times \hat{e}^3|}; \hat{e}_2 = \hat{e}_3 \times \hat{e}_1$$

The elasticity tensor at each quadrature point (qp) is computed in ADComputeIsotropicElasticityTensorShell.

Stress calculation

The stress tensor at each qp is computed as:

$$\tilde{\tau}^{ij} = \tilde{C}^{ijkl} \tilde{\epsilon}_{kl}$$

The total stress at each qp is computed in ADComputeShellStress.

Stress divergence

Finally the residuals are computed and assembled in ADStressDivergenceShell for both small and large strain scenarios. For small strain scenarios, the residual at each qp is computed by multiplying the stress tensor at that qp with the corresponding components of the B matrix. Additionally, for large strain scenarios, the old stress tensor is multiplied with corresponding components of the $B_{NL}$ matrix and that is also added to the residual at each qp. These two components for the residual are same as $K_{L} u$ and $K_{NL} u$ in equation 25 of Dvorkin and Bathe (1984).

Note that there are quadrature points both along the plane as well as along the thickness of the element. The order of Gauss quadrature rule along the thickness is provided as input by the user to all the above mentioned objects.

Inertia

This element's generalized inertia forces are obtained directly from the kinematic of the shell by computing the inertial forces component-wise. Inertia forces can be expressed in terms of a mass matrix as described in Bolourchi (1979). The mass matrix takes the following form of a volume integral:

$$\mathbf{M} = \int_{V}^{} \rho \mathbf{H}^{T} \mathbf{H} dV$$

where $\mathbf{H}$ is an element-wise matrix that contains the interpolation functions for displacement and rotational degrees of freedom. In the implementation, the volume integral is simplified, and it is assumed that the element thickness is constant. Similarly, one orientation matrix is used per element, that is, changes in orientation in a single element's geometry are neglected.

References