# An SPH formulation for general plate and shell structures with finite
deformation and large rotation

Dong Wu<EMAIL_ADDRESS>Chi Zhang<EMAIL_ADDRESS>Xiangyu Hu
<EMAIL_ADDRESS>TUM School of Engineering and Design, Technical University
of Munich, 85748 Garching, Germany  Huawei Technologies Munich Research
Center, 80992 Munich, Germany

###### Abstract

In this paper, we propose a reduced-dimensional smoothed particle
hydrodynamics (SPH) formulation for quasi-static and dynamic analyses of plate
and shell structures undergoing finite deformation and large rotation. By
exploiting Uflyand–Mindlin plate theory, the present surface-particle
formulation is able to resolve the thin structures by using only one layer of
particles at the mid-surface. To resolve the geometric non-linearity and
capture finite deformation and large rotation, two reduced-dimensional linear-
reproducing correction matrices are introduced, and weighted non-singularity
conversions between the rotation angle and pseudo normal are formulated. A new
non-isotropic Kelvin-Voigt damping is proposed especially for the both thin
and moderately thick plate and shell structures to increase the numerical
stability. In addition, a shear-scaled momentum-conserving hourglass control
algorithm with an adaptive limiter is introduced to suppress the mismatches
between the particle position and pseudo normal and those estimated with the
deformation gradient. A comprehensive set of test problems, for which the
analytical or numerical results from literature or those of the volume-
particle SPH model are available for quantitative and qualitative comparison,
are examined to demonstrate the accuracy and stability of the present method.

###### keywords:

SPH , Uflyand–Mindlin plate theory , Finite deformations , Geometric non-
linearity , Thin and moderately thick plate/shell structures , Quasi-static
and dynamic analyses , Reduced-dimensional linear-reproducing correction
matrices , Non-singularity , Non-isotropic damping , Hourglass modes

††journal: Elsevier Journal

##  1 Introduction

For computational continuum dynamics, as alternatives to conventional mesh-
based methods, e.g. finite element method (FEM) and finite volume method
(FVM), meshless methods have flourished in the past decades [1, 2, 3, 4].
Smoothed particle hydrodynamics (SPH), initially developed by Lucy [5] and
Gingold and Monaghan [6] for astrophysical simulations, is one typical
example. In SPH, the continuum is modeled by particles associated with
physical properties such as mass and velocity, and the governing equations are
discretized in the form of particle interactions using a Gaussian-like kernel
function [7, 3, 8]. Since a significant number of physical system abstractions
can be realized through particle interactions, SPH has been used to model
multi-physical systems within a unified computational framework [9], which is
able to achieve seamless monolithic, strong and conservative coupling [10,
11].

To achieve such a unified computational framework, it is crucial to discretize
all relevant physics equations using effective and efficient SPH methods. In
the case of plate and shell structures which are omnipresent thin structures
in scientific and engineering fields such as shipbuilding [12, 13], aerospace
[14], and medical treatment [15], etc., the traditional full-dimensional or
volume-particle SPH method, is not computationally efficient [16]. Since there
are well-developed and matured reduced-dimensional theories, such as
Kirchhoff-Love [17] and Uflyand-Mindlin (or Mindlin-Reissner) [18, 19, 20,
21], for plate and shell structures based on mid-surface reconstruction, it is
expected to develop the computationally much more efficient reduced-
dimensional or surface-particle SPH method with a single-layer of particles
only.

The early meshless methods for plates and shells were based on Petrov or
element-free Galerkin formulation [22, 23, 24, 25], or the reproducing kernel
particle method [26, 27, 28]. As for SPH, Maurel and Combescure [29] first
developed a surface-particle SPH method for total Lagrangian quasi-static and
dynamic analyses of moderately thick plates and shells based on the Uflyand-
Mindlin theory and the assumption of small deformation. In their work, besides
an artificial viscosity term to alleviate numerical instability issues, a
stress point method is applied to temper hourglass or zero-energy modes which
exhibit in the traditional SPH method using collocated particles for both
deformation and stress. While being effective on preventing zero-energy modes,
using stress points may faces several issues, such as how to locate or
generate these points for complex geometries, complicated numerical algorithms
and the compensation of computational efficiency [30, 31]. Nevertheless, this
method was later applied in large deformation analyses by Ming et al . [32]
and dynamic damage-fracture analyses by Caleyron et al. [33]. Lin et al. [34]
developed a similar method for quasi-static analyses, but applied an
artificial viscosity term based on membrane and shearing decomposition. Ming
et al. [35] first considered finite deformation by taking all strain terms
into account with the help of Gauss-Legendre quadrature for more accurately
capturing of non-linear stress. In all of the surface-particle SPH methods
mentioned above, the rotation angles of mid-surface are directly obtained from
the pseudo normal in governing equations under the assumption of small
rotation.

In this work, we propose a collocated surface-particle SPH formulation for
total Lagrangian quasi-static and dynamic analyses of general plate or shell
structures, which may be thin or have moderate thickness, involving finite
deformation or/and large rotation. First, to better resolve the geometric non-
linearity induced by finite deformation and large rotation, two new reduced-
dimensional correction matrices for linearly reproducing position and normal
direction are introduced, and a weighted conversion algorithm, which achieves
non-singularity under large rotation, is proposed. Second, a new non-isotropic
Kelvin-Voigt damping base on Ref. [36] is proposed for achieve good numerical
stability for both thin and moderately thick plate or shell structures. Third,
in order to address hourglass modes using collocated particles only other than
introducing extra stress points, drawing the inspiration from Refs. [37, 30],
a shear-scaled momentum-conserving formulation with an adaptive limiter is
developed by mitigating the discrepancy between the actual particle position
and pseudo normal and those estimated by the deformation gradient. A set of
numerical examples involving quasi-static and dynamic analyses for both thin
and moderately thick plate or shell structures are given. The results are
compared with analytical, numerical solutions in literature or/and those
obtained by the volume-particle SPH method to demonstrate the numerical
accuracy and stability of the present method.

The remainder of this manuscript is organized as follows. Section 2 introduces
the theoretical model of plates and shells, including the kinematics,
constitutive relation, stress correction and conservation equations. The
proposed surface-particle SPH formulation, including the reduced-dimensional
linear-producing correction matrices, weighted conversion algorithm, non-
isotropic damping and momentum-conserving hourglass control, is described in
Section 3. Numerical examples are presented and discussed in Section 4 and
then concluding remarks are given in Section 5. For a better comparison and
future opening for in-depth studies, all the computational codes of this work
are released in the open-source repository of SPHinXsys [38, 9] at
https://github.com/Xiangyu-Hu/SPHinXsys.

##  2 Theoretical models

We first introduce the theoretical mode of 3D plate, and then that of 3D shell
in which material points may possess different initial normal directions
leading to different initial local coordinate systems. After that, we briefly
describe the 2D plate and shell models, which resolve the plane strain
problem, as a simplification of the 3D counterparts.

###  2.1 3D plate model

We consider the Uflyand–Mindlin plate theory [18, 19] to account for
transverse shear stress which is significant for moderately thick plates. The
theory implies that the plate behavior can be represented by one layer of
material points at its mid-surface, as shown in Figure 1.

####  2.1.1 Kinematics

We introduce $\bm{X}=\left(X,Y,Z\right)$ to represent the global coordinate
system, and $\bm{\xi}=\left(\xi,\eta,\zeta\right)$ and
$\bm{x}=\left(x,y,z\right)$, associated with so-called pseudo-normal vector
$\bm{n}$, to denote the initial and current local coordinate systems,
respectively. Note that the initial local coordinate system is same with the
global one for plate. Each material point possesses five degrees of freedom,
viz., three translations $\bm{u}=\left\\{u,v,w\right\\}^{\operatorname{T}}$
and two rotations
$\bm{\theta}=\left\\{\theta,\varphi\right\\}^{\operatorname{T}}$ expressed in
the global coordinates. Positive values of $\theta$ and $\varphi$ indicate
that the plate is rotated anticlockwise around the coordinate axis when the
axis points toward the observer and the coordinate system is right-handed. The
two rotations are used to update the pseudo-normal
$\bm{n}=\left\\{n_{1},n_{2},n_{3}\right\\}^{\operatorname{T}}$ which is also
defined in the global coordinate system and remains straight but is not
necessarily perpendicular to the mid-surface, i.e., the pseudo normal may be
different with the real normal $\bm{n}_{r}$, as shown in Figure 1. Note that
$\bm{n}^{0}=\left\\{0,0,1\right\\}^{\operatorname{T}}$ denotes the pseudo-
normal in the initial configuration with the superscript
$\left(\bullet\right)^{0}$ denoting the initial configuration.

Figure 1: Schematic of a 3D plate model.

For a 3D plate, the position $\bm{r}$ of a material point at a distance $\chi$
away from the mid-surface along the pseudo normal $\bm{n}$ can be expressed as

$\bm{r}\left(\xi,\eta,\chi,t\right)=\bm{r}_{m}\left(\xi,\eta,t\right)+\chi\bm{n}\left(\xi,\eta,t\right),\quad\chi\in\left[-d/2,d/2\right],$
(1)

where $d$ is the thickness, $\bm{r}_{m}$ the position of the material point at
the mid-surface with the subscript $\left(\bullet\right)_{m}$ denoting the
mid-surface. Note that since the thickness is assumed to be constant during
deformation and the pseudo normal $\bm{n}$ represents the plate thickness
direction, the distance $\chi$ is always between $-d/2$ and $d/2$. Therefore,
the displacement $\bm{u}$ of the material point can be determined by

$\bm{u}\left(\xi,\eta,\chi,t\right)=\bm{u}_{m}\left(\xi,\eta,t\right)+\chi\Delta\bm{n}\left(\xi,\eta,t\right),$
(2)

where $\Delta\bm{n}=\bm{n}-\bm{n}^{0}$. Then we can define the deformation
gradient tensor as

$\mathbb{F}=\nabla^{0}\bm{r}=\nabla^{0}\bm{u}+\mathbb{I}=\left(\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}\right),$
(3)

where $\nabla^{0}\equiv\partial/\partial\bm{\xi}$ is the gradient operator
with respect to the initial configuration, $\mathbb{I}$ the identity matrix,
and $\bm{a}_{1}$, $\bm{a}_{2}$, $\bm{a}_{3}$ are specified by

$\begin{cases}\bm{a}_{1}=\bm{r}_{m,\xi}+\chi\bm{n}_{\xi}\\\
\bm{a}_{2}=\bm{r}_{m,\eta}+\chi\bm{n}_{\eta}\\\ \bm{a}_{3}=\bm{n}\end{cases}$
(4)

with
$\nabla^{0}\bm{r}_{m}\equiv(\bm{r}_{m,\xi},\bm{r}_{m,\eta})^{\operatorname{T}}$
and $\nabla^{0}\bm{n}\equiv(\bm{n}_{\xi},\bm{r}_{\eta})^{\operatorname{T}}$.
The deformation gradient tensor can be decomposed into two components as

$\mathbb{F}=\mathbb{F}_{m}+\chi\mathbb{F}_{n},$ (5)

where
$\mathbb{F}_{m}=\left(\bm{r}_{m,\xi}^{\operatorname{T}},\bm{r}_{m,\eta}^{\operatorname{T}},\bm{n}^{\operatorname{T}}\right)$
and
$\mathbb{F}_{n}=\left(\bm{n}_{\xi}^{\operatorname{T}},\bm{n}_{\eta}^{\operatorname{T}},0\right)$.
Furthermore, the real normal $\bm{n}_{r}$ is given as

$\bm{n}_{r}=\frac{\bm{r}_{m,\xi}\times\bm{r}_{m,\eta}}{\left|\bm{r}_{m,\xi}\times\bm{r}_{m,\eta}\right|}.$
(6)

####  2.1.2 Constitutive relation

With the deformation gradient tensor $\mathbb{F}$, the Green-Lagrangian strain
tensor $\mathbb{E}$ can be obtained as

$\mathbb{E}=\frac{1}{2}\left(\mathbb{F}^{\operatorname{T}}\mathbb{F}-\mathbb{I}\right)=\frac{1}{2}\left(\mathbb{C}-\mathbb{I}\right),$
(7)

where $\mathbb{C}$ is the right Cauchy deformation gradient tensor. The
Eulerian Almansi strain black$\epsilon$ can be converted from $\mathbb{E}$ as

${\color[rgb]{1,1,1}\contour{black}{\epsilon}}=\mathbb{F}^{-\operatorname{T}}\cdot\mathbb{E}\cdot\mathbb{F}^{-1}=\frac{1}{2}\left(\mathbb{I}-\mathbb{F}^{-\operatorname{T}}\mathbb{F}^{-1}\right).$
(8)

When the material is linear and isotropic, the Cauchy stress black$\sigma$
reads

$\begin{split}{\color[rgb]{1,1,1}\contour{black}{$\sigma$}}&=K\operatorname{tr}\left({\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}\right)\mathbb{I}+2G\left({\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}-\frac{1}{3}\operatorname{tr}\left({\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}\right)\mathbb{I}\right)\\\
&=\lambda\operatorname{tr}\left({\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}\right)\mathbb{I}+2\mu{\color[rgb]{1,1,1}\contour{black}{$\epsilon$}},\\\
\end{split}$ (9)

where $\lambda$ and $\mu$ are the Lamé constants, $K=\lambda+2\mu/3$ the bulk
modulus and $G=\mu$ the shear modulus. The relationship between the two moduli
is given by

$E=2G\left(1+\nu\right)=3K\left(1-2\nu\right),$ (10)

where $E$ denotes the Young’s modulus and $\nu$ the Poisson’s ratio.

####  2.1.3 Stress correction

As the thickness is significantly less than the length and width of plate, the
following boundary conditions hold when the plate is free from external forces
on its surfaces where $\chi=\pm\frac{d}{2}$ or $z=\pm\frac{d}{2}$

$\left.\sigma^{l}_{xz}\right|_{z=\pm\frac{d}{2}}=0,\quad\left.\sigma^{l}_{yz}\right|_{z=\pm\frac{d}{2}}=0,$
(11)
$\left.\sigma^{l}_{zz}\right|_{z\in\left[-\frac{d}{2},\frac{d}{2}\right]}=0,$
(12)

with the superscript $\left(\bullet\right)^{l}$ denoting the current local
coordinates. Taking the boundary condition Eq. (12) and constitutive Eq. (9)
into account, the following relation of strains holds [26]

$\bar{\epsilon}^{l}_{zz}=\frac{-\nu\left(\epsilon^{l}_{xx}+\epsilon^{l}_{yy}\right)}{1-\nu},$
(13)

where the current local strain
${\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}^{l}$ is obtained by

${\color[rgb]{1,1,1}\contour{black}{\epsilon}}^{l}=\mathbb{Q}{\color[rgb]{1,1,1}\contour{black}{\epsilon}}\mathbb{Q}^{\operatorname{T}}.$
(14)

Here, $\mathbb{Q}$ is the orthogonal transformation matrix from the global to
current local coordinates. Following Batoz and Dhatt [39], $\mathbb{Q}$ can be
given as

$\mathbb{Q}=\begin{bmatrix}n_{3}+\frac{(n_{2})^{2}}{1+n_{3}}&-\frac{n_{1}n_{2}}{1+n_{3}}&-n_{1}\\\
-\frac{n_{1}n_{2}}{1+n_{3}}&n_{3}+\frac{(n_{1})^{2}}{1+n_{3}}&-n_{2}\\\
n_{1}&n_{2}&n_{3}\\\ \end{bmatrix}.$ (15)

To satisfy the boundary conditions of Eq. (11), the transverse shear stress
should be corrected as [40]

$\bar{\sigma}^{l}_{xz}=\bar{\sigma}^{l}_{zx}=\kappa\sigma^{l}_{xz},\quad\bar{\sigma}^{l}_{yz}=\bar{\sigma}^{l}_{zy}=\kappa\sigma^{l}_{yz},$
(16)

where $\kappa$ denotes the shear correction factor which is typically set to
$5/6$ for the rectangular section of the isotropic plate. Taking the corrected
strain $\bar{{\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}}^{l}$ into
constitutive Eq. (9) and then applying Eq. (16), the corrected current local
Cauchy stress $\bar{{\color[rgb]{1,1,1}\contour{black}{$\sigma$}}}^{l}$ is
obtained.

####  2.1.4 Conservation equations

The mass conservation equation can be written as

$\rho=J_{m}^{-1}\rho^{0},$ (17)

where $J_{m}=\det(\mathbb{F}_{m})$, $\rho^{0}$ and $\rho$ represent the
initial and current densities, respectively. The momentum conservation
equation is

$\rho\ddot{\bm{u}}^{l}=\nabla\cdot\left(\bar{{\color[rgb]{1,1,1}\contour{black}{\sigma}}}^{l}\right)^{\operatorname{T}}$
(18)

or

$\rho\begin{bmatrix}\ddot{u}^{l}\\\ \ddot{v}^{l}\\\
\ddot{w}^{l}\end{bmatrix}=\begin{bmatrix}\frac{\partial\bar{\sigma}_{xx}^{l}}{\partial
x}+\frac{\partial\bar{\sigma}_{xy}^{l}}{\partial
y}+\frac{\partial\bar{\sigma}_{xz}^{l}}{\partial z}\\\
\frac{\partial\bar{\sigma}_{yx}^{l}}{\partial
x}+\frac{\partial\bar{\sigma}_{yy}^{l}}{\partial
y}+\frac{\partial\bar{\sigma}_{yz}^{l}}{\partial z}\\\
\frac{\partial\bar{\sigma}_{zx}^{l}}{\partial
x}+\frac{\partial\bar{\sigma}_{zy}^{l}}{\partial
y}+\frac{\partial\bar{\sigma}_{zz}^{l}}{\partial z}\end{bmatrix}.$ (19)

With Eqs. (11) and (12), we can integrate Eq. (19) along $\chi$ or
$z\in\left[-d/2,d/2\right]$ as

$d\rho\begin{bmatrix}\ddot{u}_{m}^{l}\\\ \ddot{v}_{m}^{l}\\\
\ddot{w}_{m}^{l}\end{bmatrix}=\begin{bmatrix}\frac{\partial
N_{xx}^{l}}{\partial x}+\frac{\partial N_{xy}^{l}}{\partial y}\\\
\frac{\partial N_{yx}^{l}}{\partial x}+\frac{\partial N_{yy}^{l}}{\partial
y}\\\ \frac{\partial N_{zx}^{l}}{\partial x}+\frac{\partial
N_{zy}^{l}}{\partial y}\end{bmatrix},$ (20)

where the stress resultant $\mathbb{N}^{l}$ is calculated by the
Gauss–Legendre quadrature rule as

$\mathbb{N}^{l}=\int_{-d/2}^{d/2}\bar{{\color[rgb]{1,1,1}\contour{black}{\sigma}}}^{l}\left(z\right)dz=\sum\limits_{ip=1}^{N}\bar{{\color[rgb]{1,1,1}\contour{black}{\sigma}}}^{l}\left(z_{ip}\right)A_{ip}.$
(21)

Here, $z_{ip}$ is the integral point, $A_{ip}$ the weight, and $N$ the number
of the integral point. Since the quadrature rule is conducted to yield an
exact result for polynomials of degree $2N-1$ or lower [41], $N$ is determined
by the applied constitutive relation.

By multiplying both sides of Eq. (18) by $z$ and integrating along
$z\in\left[-d/2,d/2\right]$, the angular momentum conservation equation can be
obtained as

$\frac{d^{3}}{12}\rho\begin{bmatrix}\ddot{n}_{1}^{l}\\\ \ddot{n}_{2}^{l}\\\
\ddot{n}_{3}^{l}\end{bmatrix}=\begin{bmatrix}\frac{\partial
M_{xx}^{l}}{\partial x}+\frac{\partial M_{xy}^{l}}{\partial y}\\\
\frac{\partial M_{yx}^{l}}{\partial x}+\frac{\partial M_{yy}^{l}}{\partial
y}\\\ \frac{\partial M_{zx}^{l}}{\partial x}+\frac{\partial
M_{zy}^{l}}{\partial y}\\\ \end{bmatrix}+\begin{bmatrix}-N_{xz}^{l}\\\
-N_{yz}^{l}\\\ 0\end{bmatrix},$ (22)

where the moment resultant $\mathbb{M}^{l}$ is calculated as

$\mathbb{M}^{l}=\int_{-d/2}^{d/2}z\bar{{\color[rgb]{1,1,1}\contour{black}{\sigma}}}^{l}\left(z\right)dz=\sum\limits_{ip=1}^{N}z_{ip}\bar{{\color[rgb]{1,1,1}\contour{black}{\sigma}}}^{l}\left(z_{ip}\right)A_{ip}.$
(23)

Note that

$\int_{-d/2}^{d/2}z\frac{\partial\bar{\sigma}_{xz}^{l}}{\partial
z}dz={\Big{[}}z\bar{\sigma}_{xz}^{l}{\Big{]}}_{-d/2}^{d/2}-\int_{-d/2}^{d/2}\bar{\sigma}_{xz}^{l}dz=-N_{xz}^{l}.$
(24)

Therefore, the two governing equations, including the evolution of mid-surface
displacement and pseudo normal, respectively, for the 3D plate can be
described as

$\begin{cases}d\rho\bm{\ddot{u}}_{m}^{l}=\nabla^{l}\cdot\left(\mathbb{N}^{l}\right)^{\operatorname{T}}\\\
\frac{d^{3}}{12}\rho\bm{\ddot{n}}^{l}=\nabla^{l}\cdot\left(\mathbb{M}^{l}\right)^{\operatorname{T}}+{\bf{Q}}^{l},\end{cases}$
(25)

where

$\mathbb{N}^{l}=\begin{bmatrix}N_{xx}^{l}&N_{xy}^{l}&0\\\
N_{yx}^{l}&N_{yy}^{l}&0\\\ N_{zx}^{l}&N_{zy}^{l}&0\\\
\end{bmatrix},\mathbb{M}^{l}=\begin{bmatrix}M_{xx}^{l}&M_{xy}^{l}&0\\\
M_{yx}^{l}&M_{yy}^{l}&0\\\ M_{zx}^{l}&M_{zy}^{l}&0\\\
\end{bmatrix},{\bf{Q}}^{l}=\begin{bmatrix}-N_{xz}^{l}\\\ -N_{yz}^{l}\\\ 0\\\
\end{bmatrix}.$ (26)

In total Lagrangian formulation, the conservation equations above are
converted into

$\begin{cases}d\rho^{0}\bm{\ddot{u}}_{m}=\left(\mathbb{F}_{m}\right)^{\operatorname{-T}}\nabla^{0}\cdot\left(J_{m}\mathbb{N}^{\operatorname{T}}\right)\\\
\frac{d^{3}}{12}\rho^{0}\bm{\ddot{n}}=\left(\mathbb{F}_{m}\right)^{\operatorname{-T}}\nabla^{0}\cdot\left(J_{m}\mathbb{M}^{\operatorname{T}}\right)+J_{m}\mathbb{Q}^{\operatorname{T}}{\bf{Q}}^{l},\end{cases}$
(27)

where $\mathbb{N}=\mathbb{Q}^{\operatorname{T}}\mathbb{N}^{l}\mathbb{Q}$ and
$\mathbb{M}=\mathbb{Q}^{\operatorname{T}}\mathbb{M}^{l}\mathbb{Q}$ are the
stress and moment resultants, respectively, in global coordinates.

###  2.2 3D shell model

Based on the 3D plate model, the 3D shell model is obtained by introducing the
initial local coordinate system and the transformation matrix from the global
to initial local coordinate system.

####  2.2.1 Kinematics

The kinematics of shell can be constructed in the initial local coordinates
denoted with the superscript $\left(\bullet\right)^{L}$. Each material point
possesses five degrees of freedom, viz., three translations
$\bm{u}^{L}=\left\\{u^{L},v^{L},w^{L}\right\\}^{\operatorname{T}}$ and two
rotations
$\bm{\theta}^{L}=\left\\{\theta^{L},\varphi^{L}\right\\}^{\operatorname{T}}$
as shown in Figure 2. The pseudo-normal vector is also presented in initial
local coordinates by
$\bm{n}^{L}=\left\\{n^{L}_{1},n^{L}_{2},n^{L}_{3}\right\\}^{\operatorname{T}}$,
especially denoted by $\bm{n}^{0,L}=\left\\{0,0,1\right\\}^{\operatorname{T}}$
in the initial local configuration.

Figure 2: Schematic of a 3D shell model.

The local position $\bm{r}^{L}$ of a material point can be expressed as

$\bm{r}^{L}\left(\xi,\eta,\chi,t\right)=\bm{r}^{L}_{m}\left(\xi,\eta,t\right)+\chi\bm{n}^{L}\left(\xi,\eta,t\right),\quad\chi\in\left[-d/2,d/2\right].$
(28)

The local displacement $\bm{u}^{L}$ can thus be obtained by

$\bm{u}^{L}\left(\xi,\eta,\chi,t\right)=\bm{u}^{L}_{m}\left(\xi,\eta,t\right)+\chi\Delta\bm{n}^{L}\left(\xi,\eta,t\right),$
(29)

where $\Delta\bm{n}^{L}=\bm{n}^{L}-\bm{n}^{0,L}$. Similar to 3D plates, the
local deformation gradient tensor of 3D shells can be defined as

$\mathbb{F}^{L}=\nabla^{0,L}\bm{r}^{L}+\nabla^{0,L}\bm{n}^{L}-\nabla^{0,L}\bm{n}^{0,L}=\left({\bm{a}^{L}_{1}},{\bm{a}^{L}_{2}},{\bm{a}^{L}_{3}}\right),$
(30)

where $\nabla^{0,L}\equiv\partial/\partial\bm{\xi}$ is the gradient operators
defined in the initial local configuration, and $\bm{a}^{L}_{1}$,
$\bm{a}^{L}_{2}$, $\bm{a}^{L}_{3}$ are detailed by

$\begin{cases}\bm{a}^{L}_{1}=\bm{r}^{L}_{m,\xi}+\chi\bm{n}^{L}_{\xi}-\chi\bm{n}^{0,L}_{\xi}\\\
\bm{a}^{L}_{2}=\bm{r}^{L}_{m,\eta}+\chi\bm{n}^{L}_{\eta}-\chi\bm{n}^{0,L}_{\eta}\\\
\bm{a}^{L}_{3}=\bm{n}^{L}.\end{cases}$ (31)

####  2.2.2 Stress correction and conservation equation

With the local deformation gradient tensor $\mathbb{F}^{L}$, the local
Eulerian Almansi strain ${\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}^{L}$
can be calculated by the Eq. (8). After that, the current local
${\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}^{l}$ is obtained according to
the coordinate transformation as

${\color[rgb]{1,1,1}\contour{black}{\epsilon}}^{l}=\mathbb{Q}\left(\mathbb{Q}^{0}\right)^{\operatorname{T}}{\color[rgb]{1,1,1}\contour{black}{\epsilon}}^{L}\mathbb{Q}^{0}\mathbb{Q}^{\operatorname{T}},$
(32)

where $\mathbb{Q}^{0}$, the orthogonal transformation matrix from the global
to initial local coordinates, is calculated from Eq. (15) while the current
pseudo normal $\bm{n}$ is replaced by the initial one $\bm{n}^{0}$. And then
the corrected strain
$\bar{{\color[rgb]{1,1,1}\contour{black}{$\epsilon$}}}^{l}$ is estimated by
applying Eq. (13). After getting the current local Cauchy stress
${\color[rgb]{1,1,1}\contour{black}{$\sigma$}}^{l}$ by Eq. (9), the corrected
one $\bar{{\color[rgb]{1,1,1}\contour{black}{$\sigma$}}}^{l}$ is obtained by
Eq. (16).

Note that the total Lagrangian conservation equations of a 3D shell has the
same form as Eqs. (27) with
$\mathbb{F}_{m}=\left(\mathbb{Q}^{0}\right)^{\operatorname{T}}\mathbb{F}_{m}^{L}\mathbb{Q}^{0}$.

###  2.3 2D plate/shell model

If a plate/shell is assumed to be a strip that is very long and has a finite
width, and the transverse load is assumed to be uniform along the length, the
analysis can be simplified at any cross section as a plane strain problem
[42].

The kinematics of 2D plate and shell can all be built in initial local
coordinates, as the transformation matrix $\mathbb{Q}^{0}$ from the global to
initial local coordinates for the plate is the unit matrix. The 2D model is in
the global $X$-$Z$ plane, and each material point possesses three degrees of
freedom, viz., two translations
$\bm{u}^{L}=\left\\{u^{L},w^{L}\right\\}^{\operatorname{T}}$ and one rotation
$\bm{\theta}^{L}=\left\\{\varphi^{L}\right\\}^{\operatorname{T}}$ expressed in
the initial local coordinates. The pseudo-normal vector is presented in the
initial local coordinates by
$\bm{n}^{L}=\left\\{n^{L}_{1},n^{L}_{3}\right\\}^{\operatorname{T}}$,
especially denoted by $\bm{n}^{0,L}=\left\\{0,1\right\\}^{\operatorname{T}}$
in the initial local configuration. Similar to 3D model, the local position
$\bm{r}^{L}$ of a material point can be expressed as

$\bm{r}^{L}\left(\xi,\chi,t\right)=\bm{r}^{L}_{m}\left(\xi,t\right)+\chi\bm{n}^{L}\left(\xi,t\right),\quad\chi\in\left[-d/2,d/2\right],$
(33)

the local displacement $\bm{u}^{L}$ can be evaluated as

$\bm{u}^{L}\left(\xi,\chi,t\right)=\bm{u}^{L}_{m}\left(\xi,t\right)+\chi\Delta\bm{n}^{L}\left(\xi,t\right).$
(34)

and the local deformation gradient tensor is written as

$\mathbb{F}^{L}=\nabla^{0,L}\bm{r}^{L}+\nabla^{0,L}\bm{n}^{L}-\nabla^{0,L}\bm{n}^{0,L}=\left({\bm{a}^{L}_{1}},{\bm{a}^{L}_{3}}\right),$
(35)

where $\bm{a}^{L}_{1}$ and $\bm{a}^{L}_{3}$ are given by

$\begin{cases}\bm{a}^{L}_{1}=\bm{r}^{L}_{m,\xi}+\zeta\bm{n}^{L}_{\xi}-\zeta\bm{n}^{0,L}_{\xi}\\\
\bm{a}^{L}_{3}=\bm{n}^{L}.\end{cases}$ (36)

The coordinate transformation matrix $\mathbb{Q}$ from global to current local
coordinates is simplified from Eqs. (15) as

$\mathbb{Q}=\begin{bmatrix}n_{3}&-n_{1}\\\ n_{1}&n_{3}\\\ \end{bmatrix},$ (37)

and the 2D transformation matrix $\mathbb{Q}^{0}$ from global to initial local
coordinates can also calculated by Eq. (37) while the current pseudo normal
$\bm{n}$ is replaced by the initial one $\bm{n}^{0}$. The corrected relation
of strains is simplified from Eq. (13) as

$\bar{\epsilon}^{l}_{zz}=\frac{-\nu\epsilon^{l}_{xx}}{1-\nu}.$ (38)

Finally, the 2D conservation equation is identical to 3D Eq. (27) with

$\mathbb{N}^{l}=\begin{bmatrix}N^{l}_{xx}&0\\\ N^{l}_{zx}&0\\\
\end{bmatrix},\mathbb{M}^{l}=\begin{bmatrix}M^{l}_{xx}&0\\\ M^{l}_{zx}&0\\\
\end{bmatrix},{\bf{Q}}^{l}=\begin{bmatrix}-N^{l}_{xz}\\\ 0\\\ \end{bmatrix}.$
(39)

##  3 SPH method for plate and shell structures

In this section, we first introduce the reduced-dimensional SPH method, and
detail the proposed formulations for plate and shell structures, including the
discretization of conservation equations, non-singular conversion algorithm
for the kinematics between rotation angles and pseudo normal, and the
algorithms to increase numerical stability and alleviate hourglass modes.
After that, the time-integration schemes are presented.

###  3.1 Reduced-dimensional SPH method

In full-dimensional SPH method, the smoothed field $f(\bm{r})$ is obtained as

$f(\bm{r})=\int_{\Omega}f(\bm{r}^{\prime})W(\bm{r}-\bm{r}^{\prime},h)d\bm{r}^{\prime},$
(40)

where $f(\bm{r}^{\prime})$ is the original continuous field before smoothing,
$\Omega$ the entire space and $W(\bm{r}-\bm{r}^{\prime},h)$ a Gaussian-like
kernel function with smoothing length $h$ denoting the compact support. By
carrying out the integration of Eq. (40) along the thickness of the
plate/shell structure, we can obtained the reduced-dimensional smoothed field
by

$f\left(\mathbf{r}\right)\approx\int_{\widehat{\Omega}}f\left(\mathbf{r}^{\prime}\right)\widehat{W}\left(\mathbf{r}-\mathbf{r}^{\prime},h\right)d\mathbf{r}^{\prime},$
(41)

where $\widehat{\Omega}$ denotes the reduced space and
$\widehat{W}(\bm{r}-\bm{r}^{\prime},h)$ the reduced kernel function. Note that
Eqs. (40) and (41) have identical forms of formulation. A reduced-dimensional
fifth-order Wendland kernel [43] reads

$\widehat{W}(q,h)=\alpha\begin{cases}\left(1+2q\right)\left(1-q/2\right)^{4}&\text{if}\quad
0\leq q\leq 2\\\ 0&\text{otherwise}\end{cases},$ (42)

where $q=\left|\mathbf{r}-\mathbf{r}^{\prime}\right|/h$ and the constant
$\alpha$ is equal to $\frac{3}{4h}$ and $\frac{7}{4\pi h^{2}}$ for 2D and 3D
problems, respectively. Also note that the reduced kernel function has
identical form with the full-dimensional counterpart except different
dimensional normalizing constant parameter, allowing the integration of unit
can be satisfied in the reduced space. Due to the almost identical forms, in
present work from here, we do not identify the full- and reduced-dimensional
formulations unless explicitly mentioned.

In the reduced-dimensional SPH method, similarly to the full-dimensional
counterpart [7], the gradient of the variable field $f(\bm{r})$ at a surface
particle $i$ can be approximated as

$\begin{split}\nabla f_{i}&=\int_{\Omega}\nabla
f(\bm{r})W(\bm{r}_{i}-\bm{r},h)d\bm{r}\\\ &=-\int_{\Omega}f(\bm{r})\nabla
W(\bm{r}_{i}-\bm{r},h)d\bm{r}\approx-\sum\limits_{j}f_{j}\nabla
W_{ij}V_{j},\end{split}$ (43)

where $V$ is the reduced particle volume, i.e. length and area for 2D and 3D
problems, respectively. Here, the summation is conducted over all the
neighboring particles $j$ located at the support domain of the particle $i$,
and $\nabla W_{ij}=-\frac{\partial W\left(\bm{r}_{ij},h\right)}{\partial
r_{ij}}\bm{e}_{ij}$ is the gradient of the kernel function with
$\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ and
$\bm{e}_{ij}=\bm{r}_{ij}/|\bm{r}_{ij}|$ denoting the unit vector pointing from
particle $j$ to $i$. Equation (43) can be modified into a strong form as

$\nabla f_{i}=\nabla f_{i}-f_{i}\nabla 1\approx\sum\limits_{j}f_{ij}\nabla
W_{ij}V_{j},$ (44)

where $f_{ij}=f_{i}-f_{j}$ is the interparticle difference value. This strong-
form derivative operator can be used to determine the local structure of a
field, such as the deformation gradient tensor. And Eq. (43) can also be
modified into a weak form as

$\nabla f_{i}=f_{i}\nabla 1+\nabla
f_{i}\approx-\sum\limits_{j}\left(f_{i}+f_{j}\right)\nabla W_{ij}V_{j}.$ (45)

This weak-form derivative operator is applied here for solving the
conservation equations. Thanks to its anti-symmetric feature, i.e., $\nabla
W_{ij}=-\nabla W_{ji}$, the momentum conservation of the particle system is
ensured [7].

Also note that, in the present work, the reduced-dimensional SPH method is
employed for total Lagrangian formulation [44], such as Eq. (27). Therefore,
the smoothing kernel function and its derivatives are only evaluated once,
also denoted with superscript $\left(\bullet\right)^{0}$ at the initial
configuration, and kept unchanged during the simulation.

###  3.2 First-order consistency corrections

For the full-dimensional SPH in total Lagrangian formulation, in order to
remedy the 1st-order inconsistency which is caused by incomplete kernel
support at domain boundary or with irregular particle distribution, the
symmetric correction matrix $\mathbb{B}^{0}_{i}$ for each particle [45, 3] is
introduced for each particle to satisfy the linear-reproducing condition

$\left(\sum\limits_{j}{\bm{r}^{0}_{ij}\otimes\nabla^{0}W_{ij}V^{0}_{j}}\right)\mathbb{B}^{0}_{i}=\mathbb{I}.$
(46)

Then the strong-form approximations of gradient Eq. (44) is modified as

$\nabla^{0}f_{i}\approx\left(\sum\limits_{j}f_{ij}\nabla^{0}W_{ij}V^{0}_{j}\right)\mathbb{B}^{0}_{i},$
(47)

and the weak-form approximations of divergence Eq. (45) as

$\nabla^{0}\cdot
f_{i}\approx-\sum\limits_{j}\left(f_{i}\mathbb{B}^{0}_{i}+f_{j}\mathbb{B}^{0}_{j}\right)\nabla^{0}W_{ij}V^{0}_{j}.$
(48)

In the reduced-dimensional SPH, we generalize the linear-reproducing condition
as

$\left[\mathbb{G}^{\operatorname{T}}\mathbb{Q}^{0}_{i}\left(\sum\limits_{j}{\bm{q}^{0}_{ij}\otimes\nabla^{0}W_{ij}V^{0}_{j}}\right)\left(\mathbb{Q}^{0}_{i}\right)^{\operatorname{T}}\mathbb{G}\right]\mathbb{B}^{0,L}_{i}=\mathbb{K}_{i},$
(49)

where $\bm{q}^{0}_{ij}$ is the initial inter-particle difference of a linear
vector, $\mathbb{Q}^{0}_{i}$ is the transformation matrix from the global to
initial local coordinates, and $\mathbb{G}$ is a reducing matrix, i.e.,

$\mathbb{G}=\begin{bmatrix}1\\\ 0\\\
\end{bmatrix}~{}\text{and}~{}\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\
\end{bmatrix}$ (50)

for 2D and 3D problems, respectively. It ensures that the corrections are
carried out within the local reduced space. Similarly, the strong-form
approximations of gradient Eq. (44) is modified as

$\nabla^{0}f_{i}\approx\left(\sum\limits_{j}f_{ij}\nabla^{0}W_{ij}V^{0}_{j}\right)\widetilde{\mathbb{B}}^{0}_{i},$
(51)

where
$\widetilde{\mathbb{B}}^{0}_{i}=\left(\mathbb{Q}^{0}_{i}\right)^{\operatorname{T}}\mathbb{G}\mathbb{B}^{0,L}_{i}\mathbb{G}^{\operatorname{T}}\mathbb{Q}^{0}_{i}$
and the weak-form approximations of divergence Eq. (45) as

$\nabla^{0}\cdot
f_{i}\approx-\sum\limits_{j}\left(f_{i}\widetilde{\mathbb{B}}^{0}_{i}+f_{j}\widetilde{\mathbb{B}}^{0}_{j}\right)\nabla^{0}W_{ij}V^{0}_{j}.$
(52)

Here, we introduce the correction matrix
$\widetilde{\mathbb{B}}^{0}_{i}=\widetilde{\mathbb{B}}^{0,\bm{r}}_{i}$,
$\bm{q}^{0}_{ij}=\bm{r}_{ij}^{0}$ and $\mathbb{K}_{i}$ is the reduced identity
matrix denoted as

$\mathbb{K}_{i}=\mathbb{K}^{\bm{r}}=\begin{bmatrix}1\end{bmatrix}~{}\text{and}~{}\begin{bmatrix}1&0\\\
0&1\end{bmatrix}$ (53)

for 2D and 3D problems, respectively, to correct the position-based
quantities. Similarly, we introduce the correction matrix
$\widetilde{\mathbb{B}}^{0}_{i}=\widetilde{\mathbb{B}}^{0,\bm{n}}_{i}$,
$\bm{q}^{0}_{ij}=\bm{n}_{ij}^{0}$ and

$\mathbb{K}_{i}=\mathbb{K}^{\bm{n}}_{i}=\begin{bmatrix}1/R_{i}^{L}\end{bmatrix}~{}\text{and}~{}\begin{bmatrix}1/R_{1,i}^{L}&0\\\
0&1/R_{2,i}^{L}\end{bmatrix},$ (54)

where $R_{i}^{L}$, $R_{1,i}^{L}$ and $R_{2,i}^{L}$ are the curvature radii of
particle $i$ for 2D and 3D problems, respectively, to correct rotation-based
quantities.

###  3.3 Discreization of conservation equations

With two correction matrices obtained from Eq. (49) and following Eq. (52),
the momentum equations (27) are discretized as

$d\rho_{i}^{0}\bm{\ddot{u}}_{m,i}=\sum\limits_{j}\left(J_{m,i}\mathbb{N}_{i}\left(\mathbb{F}_{m,i}\right)^{\operatorname{-T}}\widetilde{\mathbb{B}}^{0,\bm{r}}_{i}+J_{m,j}\mathbb{N}_{j}\left(\mathbb{F}_{m,j}\right)^{\operatorname{-T}}\widetilde{\mathbb{B}}^{0,\bm{r}}_{j}\right)\nabla^{0}W_{ij}V_{j}^{0},$
(55)

and

$\begin{split}\frac{d^{3}}{12}\rho_{i}^{0}\bm{\ddot{n}}_{i}&=\sum\limits_{j}\left(J_{m,i}\mathbb{M}_{i}\left(\mathbb{F}_{m,i}\right)^{\operatorname{-T}}\widetilde{\mathbb{B}}^{0,\bm{n}}_{i}+J_{m,j}\mathbb{M}_{j}\left(\mathbb{F}_{m,j}\right)^{\operatorname{-T}}\widetilde{\mathbb{B}}^{0,\bm{n}}_{j}\right)\nabla^{0}W_{ij}V_{j}^{0}\\\
&+{J_{m,i}\left(\mathbb{Q}_{i}^{0}\right)}^{\operatorname{T}}\bm{Q}^{l}_{i}.\end{split}$
(56)

###  3.4 Kelvin–Voigt type damping

Following Ref. [36], when calculating the current local Cauchy stress by using
the constitutive Eq. (9), an artificial damping stress
${\color[rgb]{1,1,1}\contour{black}{$\sigma$}}_{d}^{l}$ based on the Kelvin-
Voigt type damper is introduced here as

${\color[rgb]{1,1,1}\contour{black}{\sigma}}_{d}^{l}=J_{m}^{-1}\mathbb{Q}\left(\mathbb{Q}^{0}\right)^{\operatorname{T}}\mathbb{F}^{L}\dot{\mathbb{E}}^{L}{\color[rgb]{1,1,1}\contour{black}{\gamma}}\left(\mathbb{F}^{L}\right)^{\operatorname{T}}\mathbb{Q}^{0}\mathbb{Q}^{\operatorname{T}},$
(57)

where the numerical viscosity matrix

${\color[rgb]{1,1,1}\contour{black}{\gamma}}=\begin{bmatrix}\rho ch/2&0\\\
0&\rho cs/2\\\ \end{bmatrix}~{}\text{and}~{}\begin{bmatrix}\rho ch/2&0&0\\\
0&\rho ch/2&0\\\ 0&0&\rho cs/2\\\ \end{bmatrix}$ (58)

where $c=\sqrt{K/\rho}$ and $s=\min(h,d)$, for 2D and 3D problems,
respectively. Note that, different from Ref. [36], where an isotropic
numerical damping is applied, the present damping leads to a smaller out-of-
plane contribution when $d\leq h$, which makes it suitable for both thin and
moderately thick plate and shell structures. The change rate of the Green-
Lagrangian strain tensor is given as

$\dot{\mathbb{E}}^{L}=\frac{1}{2}\left[\left(\dot{\mathbb{F}}^{L}\right)^{\operatorname{T}}\mathbb{F}^{L}+\left(\mathbb{F}^{L}\right)^{\operatorname{T}}\dot{\mathbb{F}}^{L}\right].$
(59)

Here, the change rate of the deformation gradient tensor of particle $i$ is

${\mathbb{\dot{F}}}^{L}_{i}=\nabla^{0,L}\bm{\dot{u}}_{i}^{L}=\nabla^{0}\bm{\dot{u}}^{L}_{m,i}+\chi\nabla^{0}\bm{\dot{n}}^{L}_{i},$
(60)

where

$\begin{cases}\nabla^{0}\bm{\dot{u}}^{L}_{m,i}=\mathbb{Q}_{i}^{0}\left(\sum\limits_{j}\bm{\dot{u}}_{m,ij}\otimes\nabla^{0}W_{ij}V_{j}^{0}\right)\widetilde{\mathbb{B}}_{i}^{0,\bm{r}}\left(\mathbb{Q}_{i}^{0}\right)^{\operatorname{T}}\\\
\nabla^{0}\bm{\dot{n}}^{L}_{i}=\mathbb{Q}_{i}^{0}\left(\sum\limits_{j}\bm{\dot{n}}_{ij}\otimes\nabla^{0}W_{ij}V_{j}^{0}\right)\widetilde{\mathbb{B}}_{i}^{0,\bm{n}}\left(\mathbb{Q}_{i}^{0}\right)^{\operatorname{T}}\end{cases}$
(61)

are obtained following the consistency condition Eq. (49) and the strong-form
correction Eq. (51).

###  3.5 Hourglass control

Inspired from Refs. [37, 30] in full-dimensional SPH for total Lagrangian
solid dynamics, we introduce a hourglass control algorithm here to alleviate
the hourglass modes in the dynamics of plate and shell structures.

First, we estimate the position of the inter-particle middle point linearly
using the deformation gradient tensor from particles $i$ and $j$,
respectively, as

$\bm{r}_{i+\frac{1}{2}}=\bm{r}_{m,i}-\frac{1}{2}\mathbb{F}_{m,i}\bm{r}_{m,ij}^{0},\quad\bm{r}_{j-\frac{1}{2}}=\bm{r}_{m,j}+\frac{1}{2}\mathbb{F}_{m,j}\bm{r}_{m,ij}^{0}.$
(62)

One can find that the inconsistency beyond linear estimation
$\hat{\bm{r}}_{ij}=\bm{r}_{i+\frac{1}{2}}-\bm{r}_{j-\frac{1}{2}}$ is

$\hat{\bm{r}}_{ij}=\bm{r}_{m,ij}-\frac{1}{2}\left(\mathbb{F}_{m,i}+\mathbb{F}_{m,j}\right)\bm{r}_{m,ij}^{0}.$
(63)

To suppress the position inconsistency $\hat{\bm{r}}_{ij}$, we introduce an
extra correction term to the discrete momentum conservation Eq. (55) as

$d\rho_{i}^{0}\bm{\ddot{u}}^{cor}_{m,i}=\sum\limits_{j}\alpha
G\beta_{ij}\gamma_{ij}^{\bm{r}}D\hat{\bm{r}}_{ij}\frac{\partial
W\left(\bm{r}_{ij}^{0},h\right)}{\partial r_{ij}^{0}}V_{j}^{0}\\\ $ (64)

where $\beta_{ij}=W^{0}_{ij}/W_{0}$ leads to less correction to further
neighbors,
$\gamma_{ij}^{\bm{r}}=\min\left(2\left|\hat{\bm{r}}_{ij}\right|/\left|\bm{r}_{m,ij}\right|,1\right)$
is an adaptive limiter for less correction on the domain where the
inconsistency is less pronounced, D the dimension, and parameter
$\alpha=0.002$ according to the numerical experiment and remains constant
throughout this work. Note that, since the inconsistency decreases with
decreasing particle spacing, different from Refs. [30], the present correction
is purely numerical and vanishes with increasing resolution. Similarly, for
the predicted pseudo normal, the difference of the intermediate point can be
described as

$\hat{\bm{n}}_{ij}=\bm{n}_{ij}-\bm{n}_{ij}^{0}-\frac{1}{2}\left(\mathbb{F}_{\bm{n},i}+\mathbb{F}_{\bm{n},j}\right)\bm{r}_{ij}^{0}.$
(65)

Similar with Eq. (64), the extra correction term added to the discrete angular
momentum conservation Eq. (56) is

$\frac{d^{3}}{12}\rho_{i}^{0}\bm{\ddot{n}}^{cor}_{i}=\sum\limits_{j}\alpha
Gd^{2}\beta_{ij}\gamma_{ij}^{\bm{n}}D\hat{\bm{n}}_{ij}\frac{\partial
W\left(\bm{r}_{ij}^{0},h\right)}{\partial r_{ij}^{0}}V_{j}^{0},$ (66)

where the adaptive limiter is
$\gamma_{ij}^{\bm{n}}=\min\left(2\left|\hat{\bm{n}}_{ij}\right|/\left|\bm{n}_{ij}-\bm{n}_{ij}^{0}\right|,1\right)$.
Note that, different with Refs. [37, 30], the present correction force is
introduced in particle pairwise pattern, implying momentum conservation [7].
Also note that, the correction force is scaled to the shear, rather than
Young’s, modulus, due to the fact that the hourglass modes are characterized
by shear deformation [31].

###  3.6 Conversion between rotations and pseudo normal

Different from the mid-surface displacement, which can be numerically
integrated directly from its evolution equation, the pseudo normal is not
suitable for direct numerical integration since its unit magnitude may not be
maintained strictly. Under the assumption of small rotation, one may have the
simplified relation between the pseudo normal and rotations, i.e.
$\bm{\ddot{\theta}}=(-\ddot{n}_{2},\ddot{n}_{1})$, so that one can obtain the
rotation increment, and update rotation matrix $\mathbb{R}$ using Rodrigues
formula [46, 34] and finally the integrated pseudo normal by
$\bm{n}=\mathbb{R}\bm{n}^{0}$.

In the present work, the numerical integration of pseudo normal is carried out
without the assumption of small rotation, that is, we strictly identify the
rotations and pseudo normal by using the original evolution equation and
obtain their conversion relations
$\bm{\ddot{\theta}}=\bm{\ddot{\theta}}(\bm{\ddot{n}},\bm{\dot{\theta}},\bm{\theta})$.
Different from using the rotation matrix $\mathbb{R}$ based on Rodrigues
formula, we update the pseudo normal $\bm{n}^{L}$ with [46, 40]

$\bm{n}^{L}=\mathbb{R}^{L}_{\eta}\mathbb{R}^{L}_{\xi}\bm{n}^{0,L},$ (67)

where $\mathbb{R}_{\xi}^{L}\equiv\mathbb{R}_{\xi}(\theta^{L})$ and
$\mathbb{R}^{L}_{\eta}\equiv\mathbb{R}_{\eta}(\varphi^{L})$ are the local
rotation matrices respected to the axes $\xi$ and $\eta$, respectively, or,
equivalently, with the change rate

$\bm{\dot{n}}^{L}=\dot{\mathbb{R}}^{L}_{\eta}\dot{\mathbb{R}}^{L}_{\xi},$ (68)

where
$\dot{\mathbb{R}}^{L}_{\xi}\equiv\mathbb{R}_{\xi}(\theta^{L},\dot{\theta}^{L})$
and
$\dot{\mathbb{R}}^{L}_{\eta}\equiv\mathbb{R}_{\eta}(\varphi^{L},\dot{\varphi}^{L})$.
Here, the rotations and their change rates are numerically integrated directly
with the help of conversion relations.

Specifically, for a 2D problem, $\mathbb{R}^{L}_{\xi}$ is a unit matrix, and
$\mathbb{R}^{L}_{\eta}$ can be described as

$\mathbb{R}^{L}_{\eta}=\begin{bmatrix}\cos\varphi^{L}&\sin\varphi^{L}\\\
-\sin\varphi^{L}&\cos\varphi^{L}\\\ \end{bmatrix}.$ (69)

Then, one has the relation as

$\bm{n}^{L}=(\sin\varphi^{L},\cos\varphi^{L})^{\operatorname{T}},$ (70)

its 1st-order time derivative corresponding Eq. (68)

$\dot{\bm{n}}^{L}=(\cos\varphi^{L}\cdot\dot{\varphi}^{L},-\sin\varphi^{L}\cdot\dot{\varphi}^{L})^{\operatorname{T}},$
(71)

and 2nd-order derivative

$\ddot{\bm{n}}^{L}=(-\sin\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}+\cos\varphi^{L}\cdot\ddot{\varphi}^{L},-\cos\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}-\sin\varphi^{L}\cdot\ddot{\varphi}^{L})^{\operatorname{T}}.$
(72)

Note that Eq. (72) suggests two theoretically equivalent conversion relations

$\ddot{\varphi}^{L}=\frac{\ddot{n}_{1}^{L}+\sin\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}}{\cos\varphi^{L}}\quad\text{and}\quad\ddot{\varphi}^{L}=\frac{\ddot{n}_{2}^{L}+\cos\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}}{-\sin\varphi^{L}}.$
(73)

Although each of them can be used to obtain the rotation angle $\varphi^{L}$
and its change rate with direct numerical integration and hence the pseudo
normal with Eq. (67), there are singularities at large rotation angles
$\varphi^{L}=0.5\pi+k\pi$ (1st relation) or $\varphi^{L}=k\pi$ (2nd relation)
with $k=0,1,2,3,...$ [47, 46, 48]. In order to eliminate the singularities, we
propose to uses both relations with a weighted average as

$\begin{split}\ddot{\varphi}^{L}&=\left(\cos\varphi^{L}\right)^{2}\frac{\ddot{n}_{1}^{L}+\sin\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}}{\cos\varphi^{L}}+\left(\sin\varphi^{L}\right)^{2}\frac{\ddot{n}_{2}^{L}+\cos\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}}{-\sin\varphi^{L}}\\\
&=\cos\varphi^{L}\left(\ddot{n}_{1}^{L}+\sin\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}\right)-\sin\varphi^{L}\left(\ddot{n}_{2}^{L}+\cos\varphi^{L}\cdot\left(\dot{\varphi}^{L}\right)^{2}\right),\end{split}$
(74)

which cancels both denominators. Note that the present formulation recovers
$\ddot{\varphi}^{L}=\ddot{n}_{1}^{L}$ under the assumption of small rotation.

As for 3D problems, the rotation matrices $\mathbb{R}^{L}_{\xi}$ and
$\mathbb{R}^{L}_{\eta}$ are

$\mathbb{R}^{L}_{\xi}=\begin{bmatrix}1&0&0\\\
0&\cos\theta^{L}&-\sin\theta^{L}\\\ 0&\sin\theta^{L}&\cos\theta^{L}\\\
\end{bmatrix},$ (75)

and

$\mathbb{R}^{L}_{\eta}=\begin{bmatrix}\cos\varphi^{L}&0&\sin\varphi^{L}\\\
0&1&0\\\ -\sin\varphi^{L}&0&\cos\varphi^{L}\\\ \end{bmatrix}.$ (76)

Similarly, one has the relation between rotations and pseudo normal [49]

$\bm{n}^{L}=(\cos\theta^{L}\sin\varphi^{L},-\sin\theta^{L},\cos\theta^{L}\cos\varphi^{L})^{\operatorname{T}},$
(77)

its 1st-order time derivatives corresponding Eq. (68)

$\begin{cases}\dot{n}^{L}_{1}=-\sin\theta^{L}\sin\varphi^{L}\dot{\theta}^{L}+\cos\theta^{L}\cos\varphi^{L}\dot{\varphi}^{L}\\\
\dot{n}^{L}_{2}=-\cos\theta^{L}\dot{\theta}^{L}\\\
\dot{n}^{L}_{3}=-\sin\theta^{L}\cos\varphi^{L}\dot{\theta}^{L}-\cos\theta^{L}\sin\varphi^{L}\dot{\varphi}^{L},\end{cases}$
(78)

and 2nd-order derivatives

$\begin{cases}\begin{split}\ddot{n}^{L}_{1}=&-\sin\theta^{L}\sin\varphi^{L}\ddot{\theta}^{L}-\cos\theta^{L}\sin\varphi^{L}({\dot{\theta}}^{L})^{2}-2\sin\theta^{L}\cos\varphi^{L}\dot{\theta}^{L}\dot{\varphi}^{L}\\\
&-\cos\theta^{L}\sin\varphi^{L}(\dot{\varphi}^{L})^{2}+\cos\theta^{L}\cos\varphi^{L}\ddot{\varphi}^{L}\end{split}\\\
\ddot{n}^{L}_{2}=\sin\theta^{L}(\dot{\theta}^{L})^{2}-\cos\theta^{L}\ddot{\theta}^{L}\\\
\begin{split}\ddot{n}^{L}_{3}&=-\sin\theta^{L}\cos\varphi^{L}\ddot{\theta}^{L}-\cos\theta^{L}\cos\varphi^{L}({\dot{\theta}}^{L})^{2}+2\sin\theta^{L}\cos\varphi^{L}\dot{\theta}^{L}\dot{\varphi}^{L}\\\
&-\cos\theta^{L}\cos\varphi^{L}(\dot{\varphi}^{L})^{2}-\cos\theta^{L}\sin\varphi^{L}\ddot{\varphi}^{L}.\end{split}\end{cases}$
(79)

Note that, one can obtained 3 theoretically equivalent conversion relations,
respectively, by 1st and 3rd expressions of Eq. (79) as

$\begin{cases}\ddot{\theta}^{L}=-\left(\ddot{n}^{L}_{3}\cos\varphi^{L}+\ddot{n}^{L}_{1}\sin\varphi^{L}+\left(\dot{\varphi}^{L}\right)^{2}\cos\theta^{L}+\left(\dot{\theta}^{L}\right)^{2}\cos\theta^{L}\right)/\sin\theta^{L}\\\
\ddot{\varphi}^{L}=\left(\ddot{n}^{L}_{1}\cos\varphi^{L}-\ddot{n}^{L}_{3}\sin\varphi^{L}+2\dot{\varphi}^{L}\dot{\theta}^{L}\sin\theta^{L}\right)/\cos\theta^{L},\end{cases}$
(80)

1st and 2nd expressions

$\begin{cases}\ddot{\theta}^{L}=\left(\sin\theta^{L}\left(\dot{\theta}^{L}\right)^{2}-\ddot{n}^{L}_{2}\right)/\cos\theta^{L}\\\
\begin{split}\ddot{\varphi}^{L}&=(\ddot{n}^{L}_{1}\cos\theta^{L}+\left(\dot{\varphi}^{L}\right)^{2}\cos^{2}\theta^{L}\sin\varphi^{L}+\left(\dot{\theta}^{L}\right)^{2}\sin\varphi^{L}-\ddot{n}^{L}_{2}\sin\varphi^{L}\sin\theta^{L}\\\
&+2\dot{\varphi}^{L}\dot{\theta}^{L}\cos\varphi^{L}\cos\theta^{L}\sin\theta^{L})/\cos\varphi^{L}\cos^{2}\theta^{L},\end{split}\end{cases}$
(81)

and 2nd and 3rd expressions

$\begin{cases}\ddot{\theta}^{L}=\left(\sin\theta^{L}\left(\dot{\theta}^{L}\right)^{2}-\ddot{n}^{L}_{2}\right)/\cos\theta^{L}\\\
\begin{split}\ddot{\varphi}^{L}&=-(\ddot{n}^{L}_{3}\cos\theta^{L}+\left(\dot{\varphi}^{L}\right)^{2}\cos\varphi^{L}\cos^{2}\theta^{L}+\left(\dot{\theta}^{L}\right)^{2}\cos\varphi^{L}-\ddot{n}^{L}_{2}\cos\varphi^{L}\sin\theta^{L}\\\
&-2\dot{\varphi}^{L}\dot{\theta}^{L}\cos\theta^{L}\sin\varphi^{L}\sin\theta^{L})/\sin\varphi^{L}\cos^{2}\theta^{L}.\end{split}\end{cases}$
(82)

Again, each of these conversion relations suffers singularities at large
rotations similar to that of 2D formulations. To eliminate the singularities,
we first apply the weighted average to the conversion between between
$\ddot{\theta}^{L}$ and $\bm{\ddot{n}}^{L}$ with Eqs. (80) and (81) as

$\begin{split}\ddot{\theta}^{L}&=-\left(\ddot{n}^{L}_{3}\cos\varphi^{L}+\ddot{n}^{L}_{1}\sin\varphi^{L}+\left(\dot{\varphi}^{L}\right)^{2}\cos\theta^{L}+\left(\dot{\theta}^{L}\right)^{2}\cos\theta^{L}\right)\sin\theta^{L}\\\
&+\left(\sin\theta^{L}\left(\dot{\theta}^{L}\right)^{2}-\ddot{n}^{L}_{2}\right)\cos\theta^{L}.\end{split}$
(83)

Then, for the conversion relation between $\ddot{\varphi}^{L}$ and
$\bm{\ddot{n}}^{L}$, according to Eq. (80), we can rewrite the relation as

$\cos\theta^{L}=B/\ddot{\varphi}^{L},$ (84)

where $B$ denotes the numerator of the 2nd expression in Eq. (80). We further
denote the numerators of the 2nd expressions in Eqs. (81) and (82),
respectively, as $B_{1}$ and $B_{2}$. Inserting Eq. (84) into Eqs. (81) and
(82), we have

$\begin{cases}\ddot{\varphi}^{L}=\left(B^{2}\cos\varphi^{L}\right)/B_{1}\\\
\ddot{\varphi}^{L}=\left(B^{2}\sin\varphi^{L}\right)/B_{2},\end{cases}$ (85)

and obtain the weighted average of the conversion relation as

$\ddot{\varphi}^{L}=\frac{B_{1}B^{2}\cos\varphi^{L}+B_{2}B^{2}\sin\varphi^{L}}{B_{1}^{2}+B_{2}^{2}}.$
(86)

Again, one can easily find that the present relations recover
$\bm{\ddot{\theta}}^{L}=(\ddot{\theta}^{L},\ddot{\varphi}^{L})=(-\ddot{n}_{2}^{L},\ddot{n}_{1}^{L})$
for small rotations.

###  3.7 Time stepping

For the time integration of plate and shell dynamics, we use the position-
based Verlet scheme [50]. At the beginning of each time step, the deformation
gradient, particle position, rotation angles and pseudo normal are updated to
the midpoint of the $n$-th time step as

$\begin{cases}\mathbb{F}^{L,n+\frac{1}{2}}=\mathbb{F}^{L,n}+\frac{1}{2}\Delta
t\dot{\mathbb{F}}^{L,n}\\\
\bm{r}_{m}^{n+\frac{1}{2}}=\bm{r}_{m}^{n}+\frac{1}{2}\Delta
t\bm{\dot{u}}_{m}^{n}\\\
\bm{\theta}^{L,n+\frac{1}{2}}=\bm{\theta}^{L,n}+\frac{1}{2}\Delta
t\bm{\dot{\theta}}^{L,n}\\\
\bm{n}^{L,n+\frac{1}{2}}=\bm{n}^{L,n}+\frac{1}{2}\Delta
t\bm{\dot{n}}^{L,n}.\end{cases}$ (87)

After the stress correction and Gauss-Legendre quadrature Eqs. (21) and (23),
the conservation equations are solved to obtain the $\bm{\ddot{u}}_{m}^{n+1}$
and $\bm{\ddot{n}}^{n+1}$. After transforming $\bm{\ddot{n}}^{n+1}$ to
$\bm{\ddot{n}}^{L,n+1}$, $\bm{\ddot{\theta}}^{L,n+1}$ is obtained through the
conversion relation between the pseudo normal and rotation angle, i.e., Eq.
(74) for 2D problems and Eqs. (83) and (86) for 3D problems. At this point,
the velocity and angular velocity are updated by

$\begin{cases}\bm{\dot{u}}_{m}^{n+1}=\bm{\dot{u}}_{m}^{n}+\Delta
t\bm{\ddot{u}}_{m}^{n+1}\\\
\bm{\dot{\theta}}^{L,n+1}=\bm{\dot{\theta}}^{L,n}+\Delta
t\bm{\ddot{\theta}}^{L,n+1},\end{cases}$ (88)

and the change rate of pseudo normal $\bm{\dot{n}}$ is updated by Eq. (71) or
(78). Finally, the change rate of the deformation gradient
$\dot{\mathbb{F}}^{L,n+1}$ is estimated by Eq. (60), and then the deformation
gradient, density, particle position, rotation angles and pseudo normal are
updated to the new time step with

$\begin{cases}\mathbb{F}^{L,n+1}=\mathbb{F}^{L,n+\frac{1}{2}}+\frac{1}{2}\Delta
t\dot{\mathbb{F}}^{L,n+1}\\\
\rho^{n+1}=\left(J_{m}^{n+1}\right)^{-1}\rho^{0}\\\
\bm{r}_{m}^{n+1}=\bm{r}_{m}^{n+\frac{1}{2}}+\frac{1}{2}\Delta
t\bm{\dot{u}}_{m}^{n+1}\\\
\bm{\theta}^{L,n+1}=\bm{\theta}^{L,n+\frac{1}{2}}+\frac{1}{2}\Delta
t\bm{\dot{\theta}}^{L,n+1}\\\
\bm{n}^{L,n+1}=\bm{\theta}^{L,n+\frac{1}{2}}+\frac{1}{2}\Delta
t\bm{\dot{n}}^{L,n+1}.\end{cases}$ (89)

For the numerical stability, the time-step size $\Delta t$ is given by

$\Delta t=\text{CFL}\min\left(\Delta t_{1},\Delta t_{2},\Delta t_{3}\right),$
(90)

where

$\begin{cases}\Delta
t_{1}=\min\left(\frac{h}{c_{v}+|\bm{\dot{u}}_{m}|_{max}},\sqrt{\frac{h}{|{\bm{\ddot{u}}_{m}}|_{max}}}\right)\\\
\Delta
t_{2}=\min\left(\frac{1}{c_{v}+|\bm{\dot{\theta}}^{L}|_{max}},\sqrt{\frac{1}{|{\bm{\ddot{\theta}}^{L}}|_{max}}}\right)\\\
\Delta
t_{3}=h\left(\frac{\rho\left(1-\nu^{2}\right)/E}{2+\left(\pi^{2}/12\right)\left(1-\nu\right)\left[1+1.5\left(h/d\right)^{2}\right]}\right)^{1/2}.\\\
\end{cases}$ (91)

Note that the time-step size $\Delta t_{3}$ follows the Refs. [34, 51] and
depends on the thickness and material properties, and the present Courant-
Friedrichs-Lewy (CFL) number is set as 0.6. An overview of the complete
solution procedure is presented in Algorithm 1.

1 Setup parameters and initialize the simulation;

2 Construct the particle-neighbor list and compute the kernel values;

3 Compute the correction matrices $\widetilde{\mathbb{B}}^{0,\bm{r}}$ and
$\widetilde{\mathbb{B}}^{0,\bm{n}}$ for each particle (Section 3.2);

4 while  _simulation is not finished_ do

5 Compute the time-step size $\Delta t$ using Eq. (90);

6 Update the deformation gradient tensor $\mathbb{F}^{L}$, particle position
$\bm{r}_{m}$, rotation angle $\bm{\theta}^{L}$ and pseudo normal $\bm{n}$ for
half time step $\Delta t/2$;

7 Compute and correct the Cauchy stress
${\color[rgb]{1,1,1}\contour{black}{$\sigma$}}^{l}$ (Sections 2.1.2 and
2.1.3);

8 Compute the resultants $\mathbb{N}^{l}$ and $\mathbb{M}^{l}$, and shear
force $\bm{Q}^{l}$ (Eq. (26));

9 Compute the acceleration $\bm{\ddot{u}}_{m}$ (Eqs. (55) and (64)) and
$\bm{\ddot{n}}$ (Eqs. (56) and (66));

10 Compute the angular acceleration $\bm{\ddot{\theta}}^{L}$ (Eq. (74) for 2D
problems, and Eqs. (83) and (86) for 3D problems);

11 Update the velocity $\bm{\dot{u}}_{m}$ and angular velocity
$\bm{\dot{\theta}}^{L}$ for a time step $\Delta t$;

12 Compute the change rate of pseudo normal $\bm{\dot{n}}^{L}$ using Eq. (71)
or (78);

13 Compute the change rate of the deformation gradient tensor
$\partial\mathbb{F}^{L}/\partial t$ (Eq. (60));

14 Update the deformation gradient tensor $\mathbb{F}^{L}$, density $\rho$,
particle position $\bm{r}_{m}$, rotation angle $\bm{\theta}^{L}$ and pseudo
normal $\bm{n}$ for another half time step $\Delta t/2$;

15

16 end while

Terminate the simulation.

Algorithm 1 The present SPH method for plate/shell structures.

##  4 Numerical examples

To demonstrate the accuracy and stability of the proposed surface-particle SPH
method (denoted as shell method), this section investigates a series of
benchmark tests where analytical or numerical reference data from the
literature or/and volume-particle SPH method (denoted as volume method) are
available for qualitative and quantitative comparison. The smoothing length
$h=1.15~{}dp$, where $dp$ denotes the initial particle spacing, is employed in
all the following simulations.

###  4.1 2D oscillating plate strip

The first example involves a plate strip with initial uniform transverse
velocity along the length with one edge fixed and the others free, which has
previously been theoretically [52] and numerically [53, 54, 31] investigated
in the literature. As shown in Figure 3(a), this plate strip is assumed to be
infinitely long along the $y$-axis with a finite width $a=0.2~{}\text{m}$
along the $x$-axis. To demonstrate that both thin and moderately thick
structures can be simulated, this plate strip is modeled with the thicknesses
$d=0.01~{}\text{m}$ and $0.001~{}\text{m}$. The material properties are set as
follows: density $\rho_{0}=1000.0~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=2.0~{}\text{MPa}$, and Poisson’s ratio $\nu$ varies for different cases.
Figure 3(b) shows the discrete model of the chosen cross-section with clamped
edges at $x=0$.

Figure 3: 2D oscillating plate strip: Initial configuration with width
$a=0.2~{}\text{m}$.

The transverse velocity $v_{z}$ is applied to the plate strip as

$v_{z}(x)=v_{f}c\frac{f(x)}{f(a)},$ (92)

where $v_{f}$ is a constant that varies with different cases, and

$\begin{split}f(x)&=\left(\sin(ka)+\sinh(ka)\right)\left(\cos(kx)-\cosh(kx)\right)\\\
&-\left(\cos(ka)+\cosh(ka)\right)\left(\sin(kx)-\sinh(kx)\right)\end{split}$
(93)

with $k$ determined by

$\cos(ka)\cosh(ka)=-1$ (94)

and $ka=1.875$. The frequency $\omega$ of the oscillating plate strip is
theoretically given by

$\omega^{2}=\frac{Ed^{2}k^{4}}{12\rho\left(1-\nu^{2}\right)}.$ (95)

Figure 4 shows the particles with von Mises stress $\bar{\sigma}$ contour for
the case of $d=0.001~{}\text{m}$, $v_{f}=0.01$, $\nu=0.4$, and the initial
particle spacing $dp=a/40=0.005~{}\text{m}$. It should be noted that the
present method predicts smooth deformation and stress fields without
singularities for large rotations (more than $\pi$).

Figure 4: 2D oscillating plate strip: Deformed particle configuration colored
by von Mises stress $\bar{\sigma}$ of the mid-surface at serial time instants
with the width $a=0.2~{}\text{m}$, thickness $d=0.02~{}\text{m}$,
$v_{f}=0.01$, and spatial particle resolution $a/dp=40$. The material is
modeled with density $\rho_{0}=1000.0~{}\text{kg}/\text{m}^{3}$, Young’s
modulus $E=2.0~{}\text{MPa}$, and Poisson’s ratio $\nu=0.4$.

Three different spatial resolutions, $a/dp=40$, $a/dp=80$, and $a/dp=160$, are
tested in the convergence study. Figure 5 shows the time history of vertical
position $z$ of the strip endpoint with $d=0.01~{}\text{m}$, $v_{f}=0.025$ and
$\nu=0.4$.

Figure 5: 2D oscillating plate strip: Time history of the vertical position
$z$ observed at the plate strip endpoint when $d=0.01~{}\text{m}$,
$v_{f}=0.025$ and $\nu=0.4$.

It can be observed that typical 2nd-order convergence has been achieved. In
addition, a long-term simulation is performed herein to demonstrate the
numerical stability of the proposed formulation. For quantitative validation,

Table 1: 2D oscillating plate strip: Quantitative validation of the oscillation period for $a=0.2~{}\text{m}$ and $d=0.01~{}\text{m}$ with various $v_{f}$ and $\nu$. $v_{f}$ | $\nu$ | $T_{\text{Shell model}}$ | $T_{\text{Theoretical}}$ | Error
---|---|---|---|---
0.025 | 0.22 | 0.58137 | 0.54018 | 7.63%
0.05 | 0.22 | 0.57715 | 0.54018 | 6.92%
0.1 | 0.22 | 0.56801 | 0.54018 | 5.15%
0.025 | 0.30 | 0.56804 | 0.52824 | 7.53%
0.05 | 0.30 | 0.56308 | 0.52824 | 6.60%
0.1 | 0.30 | 0.55481 | 0.52824 | 5.03%
0.025 | 0.40 | 0.54447 | 0.50752 | 7.28%
0.05 | 0.40 | 0.53683 | 0.50752 | 5.78%
0.1 | 0.40 | 0.53252 | 0.50752 | 4.93%
Table 2: 2D oscillating plate strip: Quantitative validation of the oscillation period for $a=0.2~{}\text{m}$ and $d=0.001~{}\text{m}$ with various $v_{f}$ and $\nu$. $v_{f}$ | $\nu$ | $T_{\text{Shell model}}$ | $T_{\text{Theoretical}}$ | Error
---|---|---|---|---
0.0025 | 0.22 | 5.80249 | 5.40182 | 7.42%
0.005 | 0.22 | 5.75544 | 5.40182 | 6.55%
0.01 | 0.22 | 5.64181 | 5.40182 | 4.44%
0.0025 | 0.30 | 5.66756 | 5.28243 | 7.29%
0.005 | 0.30 | 5.61006 | 5.28243 | 6.20%
0.01 | 0.30 | 5.49156 | 5.28243 | 3.96%
0.0025 | 0.40 | 5.42826 | 5.07519 | 6.96%
0.005 | 0.40 | 5.34224 | 5.07519 | 5.26%
0.01 | 0.40 | 5.27522 | 5.07519 | 3.94%

Tables 1 and 2 detail the oscillation period $T$ for a wide range of $v_{f}$
and $\nu$, obtained by the present method with the spatial particle resolution
$a/dp=160$, when thickness $d=0.01~{}\text{m}$ and $0.001~{}\text{m}$,
respectively, and the comparison to theoretical solution obtained form small
perturbation analysis. The differences, which are less than 8.00% for
$\nu=0.22$ and decrease to about 5.00% when the Poisson’s ratio is increased
to 0.4, are acceptable.

###  4.2 3D square plate

In this section, a 3D square plate under different types of boundary
conditions is considered for quasi-steady analyses, as shown in Figure 6.

Figure 6: 3D square plate: Problem setup with $a=b=254~{}\text{mm}$ and
thickness $d=25.4~{}\text{mm}$.

With side length $a=b=$ 254 mm and thickness $d=$ 25.4 mm, the plate material
is defined with density $\rho_{0}=1600~{}\text{kg}/\text{m}^{3}$, Young’s
modulus $E=53.7791~{}\text{GPa}$ and Poisson’s ratio $\nu=0.3$. Three types of
boundary conditions denoted as SS0, SS1 and SS3 following Refs. [42, 34] are
implemented as

  * $\blacktriangleright$

SS0: constrained mass center without constrained boundaries;

  * $\blacktriangleright$

SS1: $u=w=\varphi=0$ on edges parallel to $x$-axis and $v=w=\theta=0$ on edges
parallel to $y$-axis;

  * $\blacktriangleright$

SS3: $u=v=w=0$ on all edges.

Note that, for the case of SS0, the outer square ring with width $d$ is
imposed with negative pressure $q_{02}$. The uniformly distributed loads are
parameterized by the loading factors $\bar{P}$ and $\bar{P}_{1}$ as
$q_{0}=\bar{P}E\left({d/a}\right)^{4}$,
$q_{01}=\bar{P}_{1}E\left({d/a}\right)^{4}$ and
$q_{02}(2ad+2bd+4d^{2})=q_{01}ab$, so that the applied negative force along
the $z$-axis prevents the center of mass from moving.

For comprehensive validation, a convergence study of tests with SS0 is
conducted, and the results are compared with those obtained by the volume
method released in the SPHinXsys repository [9]. Figure 7 shows the particle
distribution and stress fileds under the loading factor $\bar{P}_{1}=25$ with
the spatial discretization $d/dp=8$.

Figure 7: 3D square plate: Particles colored by von Mises stress
$\bar{\sigma}$ of tests with SS0 obtained by the present shell (left) and
volume (right) methods under the loading factor $\bar{P}_{1}=25$. Note that
the left panel shows the stress $\bar{\sigma}$ of the plate mid-surface. The
material is modeled with the density $\rho_{0}=1600~{}\text{kg}/\text{m}^{3}$,
Young’s modulus $E=53.7791~{}\text{GPa}$ and Poisson’s ratio $\nu=0.3$. The
spatial particle resolution is set as $d/dp=8$. Figure 8: 3D square plate:
Load-deflection curves of tests with SS0 under three different spatial
resolutions, and their comparison with those of the volume method [9].

Figure 8 shows the non-dimensional deflection $\bar{w}_{C}=w_{C}/d$ and
$\bar{w}_{A}=w_{A}/d$ probed at the central point $C$ and corner point $A$,
respectively, obtained by both SPH shell and volume methods. It should be
emphasized that there are only quite small differences between the results of
the present reduced-dimensional and full-dimensional models.

The particles colored by von Mises stress $\bar{\sigma}$ at the mid-surface
for three spatial discretizations, $a/dp=20$, $a/dp=40$ and $a/dp=80$, with
the SS1 and SS3 boundary conditions under $\bar{P}=200$ are shown in Figure 9.

(a) SS1 boundary conditon (b) SS3 boundary conditon

Figure 9: 3D square plate: Particles colored by von Mises stress
$\bar{\sigma}$ of the mid-surface with particle refinement under the loading
factor $\bar{P}=200$ and SS1 and SS3 boundary conditions. The plate material
has parameters of the density $\rho_{0}=1600~{}\text{kg}/\text{m}^{3}$,
Young’s modulus $E=53.7791~{}\text{GPa}$ and Poisson’s ratio $\nu=0.3$.

(a) SS1 boundary conditon (b) SS3 boundary conditon

Figure 10: 3D square plate: Load-deflection curves of tests with SS1 and SS3
under three different spatial resolutions, and their comparison with that of
Reddy [42].

It can be observed that the regular particle distribution and smooth stress
field are obtained. Also, both the deformation and von Mises stress
$\bar{\sigma}$ exhibit good convergence properties with particle refinement.
In order to demonstrate the accuracy of the present method, the non-
dimensional deflections $\bar{w}_{C}$ for tests with SS1 and SS3 under various
spatial resolutions are compared to those of the Ref. [42]. As shown in Figs.
10, the numerical results quickly converge to the reference solutions obtained
by the finite element method (FEM) with increasing resolution.

###  4.3 Dynamic response of a 3D square plate

Following Ref. [55], the 3D square plate studied in Section 4.2 is considered
herein with the thickness $d=12.7~{}\text{mm}$ and Young’s modulus
$E=68.94~{}\text{GPa}$. The SS0 and SS3 boundary conditions are applied for
dynamic analyses under a step loading of uniform normal pressure
$q_{01}=q_{0}=2.068427~{}\text{MPa}$. For convergence study, three different
spatial discretizations, i.e., $d/dp=2$, $d/dp=4$ and $d/dp=8$, are
considered.

For quantitative validation, Figure 11 shows the time history of the
deflections $w_{C}$ probed at the central point $C$ and $w_{A}$ at the corner
point $A$ with SS0 boundary condition and its comparison to the results
obtained by the volume method. Also, Figure 12 shows the time history of the
deflection $w_{C}$ with the SS3 boundary condition and its comparison with
that of Ref. [55]. In general, the present results are in good agreements with
those obtained by the volume method and of Ref. [55].

Figure 11: 3D square plate with dynamic response: Time history of the
deflection $w_{C}$ and $w_{A}$ probed at the central point $C$ and corner
point $A$, respectively, with SS0 boundary condition. The material is modeled
with the density $\rho_{0}=1600~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=68.94~{}\text{GPa}$ and Poisson’s ratio $\nu=0.3$. Figure 12: 3D square
plate with dynamic response: Time history of the deflection $w_{C}$ observed
at the central point $C$ with SS3 boundary condition.

###  4.4 3D cantilevered plate

Following Refs [56, 57, 58], the static response of a 3D cantilevered plate
subjected to a distributed end shear load $f_{0}$ is considered.

Figure 13: 3D cantilevered plate: Initial configuration with the length
$a=10~{}\text{m}$, width $b=1~{}\text{m}$ and thickness $d=0.1~{}\text{m}$.

As shown in Figure 13, the plate with length $a=10~{}\text{m}$, width
$b=1~{}\text{m}$ and thickness $d=0.1~{}\text{m}$ is clamped at $y=0$, and has
material parameters of density $\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$,
Young’s modulus $E=1.2~{}\text{MPa}$ and Poisson’s ratio $\nu=0.0$. The shear
load is parameterized by a loading factor $\bar{F}$ as $f_{0}=\bar{F}EI/a^{2}$
with the inertia moment $I=bd^{3}/12$. Three different resolutions, i.e.,
$b/dp=5$, $b/dp=7$ and $b/dp=9$, are considered for convergence study.

Figure 14 shows the particles colored by the vertical displacement under
different loading factor $\bar{F}$ at the spatial resolution $b/dp=9$. A
regular particle distribution and smooth vertical displacement field are
noted.

Figure 14: 3D cantilevered plate: Particles colored by the vertical
displacement $w_{c}$ under the various loading factor $\bar{F}$ at spatial
resolution $b/dp=9$. The material is set as the density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus $E=1.2~{}\text{MPa}$
and Poisson’s ratio $\nu=0.0$. Figure 15:  3D cantilevered plate: Load-
deflection curves with three various spatial discretizations, and their
comparison with that of Payette et al. [58].

Figure 15 gives the displacement $u_{c}$ and $w_{c}$ of the point $C$, defined
in Figure 13, as a function of the loading factor $\bar{F}$ and the initial
particle spacing $dp$, and their comparison with those in Ref. [58]. It can be
noted that the displacement is converging rapidly, again at about 2nd-order,
with increasing resolution, demonstrating the accuracy of the present method.

###  4.5 Scordelis-Lo roof

As shown in Figure 16, the Scordelis-Lo roof with length $a=50~{}\text{m}$,
radius $r=25~{}\text{m}$, thickness $d=0.25~{}\text{m}$ and $\beta=40^{\circ}$
is considered herein, and the material properties are density
$\rho_{0}=36~{}\text{kg}/\text{m}^{3}$, Young’s modulus $E=432~{}\text{MPa}$
and zero Poisson’s ratio.

Figure 16: Scordelis-Lo roof: Initial configuration with the length
$a=50~{}\text{m}$, radius $r=25~{}\text{m}$, thickness $d=0.25~{}\text{m}$ and
$\beta=40^{\circ}$.

The roof is supported at its ends by fixed diaphragms, i.e. the translations
in $x$ and $z$ directions are constrained, and subjected to a gravity loading
of $g=10~{}\text{m}/\text{s}^{2}$.

The FEM solution of the vertical displacement $w$ at the midpoint of the side
edge converges to 0.3024 m as reported in Refs. [59, 60]. A sequentially
refined resolutions of $b/dp=15,20,25,30~{}\text{and}~{}40$ with $b=2r\beta$
denoting the arc length of the roof end are considered to assess the
convergence property of the present method.

Figure 17: Scordelis-Lo roof: Particles colored by the von Mises stress
$\bar{\sigma}$ of the mid-surface obtained by the present method with particle
refinement. The material is set as the density
$\rho_{0}=36~{}\text{kg}/\text{m}^{3}$, Young’s modulus $E=432~{}\text{MPa}$
and Poisson’s ratio $\nu=0.0$. Figure 18: Scordelis-Lo roof: Convergence study
of the displacement obtained by the present method with particle refinement.

Figure 17 shows the particles colored with the von Mises stress $\bar{\sigma}$
of the mid-surface obtained at different resolutions. The regular particle
distribution and smooth stress fields are noted. With increasing resolution, a
clear convergence is exhibited. The profile of displacement $w$ with varying
spatial resolution obtained by the present method is depicted in Figure 18. It
can be noted that the result converges rapidly to $w=0.2991~{}\text{m}$ when
$b/dp=40$ with 1.09% error compared to the solution of Refs. [59, 60].

###  4.6 Pinched hemispherical shell

We now consider a pinched hemispherical shell with an $18^{\circ}$ circular
cutout at its pole following Refs. [61, 62, 63, 57, 58]. As shown in Figure
19(a), the hemispherical shell with the radius $r=10.0~{}\text{m}$ and
thickness $d=0.04~{}\text{m}$ is loaded by four alternating radial point
forces $\bm{F}$, prescribed along the equator at $90^{\circ}$ intervals.

Figure 19: Pinched hemispherical shell: (a) Initial configuration with the
radius of the mid-surface $r=10.0~{}\text{m}$ and thickness
$d=0.04~{}\text{m}$, (b-d) particles colored by the von Mises stress
$\bar{\sigma}$ of the mid-surface under 3 point force magnitudes at spatial
discretication $2\pi r/dp=160$. The material parameters are set as the density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=68.25~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3$.

A linear elastic material with the density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=68.25~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3$ is applied.

Figure 19(b-d) shows the distribution of von Mises stress $\bar{\sigma}$ at
the mid-surface under varying magnitude of the point force $\bm{F}$ . The
regular particle distribution is observed, although slight stress fluctuation
is exhibited near the place where the point force is applied. For quantitative
analysis and convergence study, the radial deflections $w_{A}$ and $w_{B}$ of
monitoring points $A$ and $B$ as a function of the point force magnitude and
resolution are compared with those of Ref. [57]. Three different spatial
discretizations, i.e., $2\pi r/dp=80,160~{}\text{and}~{}240$, are considered
for convergence study. As shown in Figure 20, the results of present SPH shell
model is quickly converging to those of Ref. [57].

Figure 20: Pinched hemispherical shell: Curves of radical displacements of
points $A$ and $B$ as a function of the point force magnitude and spatial
resolution, and their comparison with those of Sze et al. [57].

###  4.7 Pulled-out cylindrical shell

A more challenging benchmark test with large displacements is considered in
this section following Refs. [29, 64]. As shown in Figure 21, a cylindrical
shell with the radius $r=5.0~{}\text{m}$, length $a=10.35~{}\text{m}$ and
thickness $d=0.094~{}\text{m}$ is subjected to a pair of point forces $\bm{F}$
which are equal in magnitude and opposite in direction.

Figure 21: Pulled-out cylindrical shell: Initial configuration with the radius
of the mid-surface $r=5.0~{}\text{m}$, length $a=10.35~{}\text{m}$ and
thickness $d=0.094~{}\text{m}$.

A linear elastic material with the density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=10.5~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3125$ is applied.

Figure 22 shows the distribution of von Mises stress $\bar{\sigma}$ at the
mid-surface under varying magnitude of the point force $\bm{F}$. The regular
particle distribution and smooth stress fields, even close to the place where
the point force is applied, are observed.

Figure 22: Pulled-out cylindrical shell: Particles colored by the von Mises
stress $\bar{\sigma}$ of the mid-surface under 3 point force magnitudes at
spatial resolution of $b/dp=240$. The material parameters are set as the
density $\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=10.5~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3125$. Figure 23: Pulled-out
cylindrical shell: Curves of radical displacements of points $A$, $B$ and $C$
as a function of the point force magnitude and spatial resolution, and their
comparison with those of Maurel and Combescure [29] and Jiang et al. [64].

For quantitative analysis and convergence study, the radial displacements
$w_{A}$, $w_{B}$ and $w_{C}$ of monitoring points $A$, $B$ and $C$ as a
function of the point force magnitude and resolution are compared with those
of Ref. [29, 64]. Three different spatial discretizations, i.e.,
$b/dp=80,160~{}\text{and}~{}240$ with $b=2\pi r$ denoting the circumference
length of the end, are considered for convergence study. As shown in Figure
23, the bifurcation point of the curve is accurately predicted, suggesting
good accuracy and robustness of the present method.

###  4.8 Pinched semi-cylindrical shell

We further consider a pinched semi-cylindrical shell with finite deformation
and rotation following Refs. [65, 66, 57, 67]. As shown in Figure 24(a), the
semi-cylindrical shell with the radius $r=1.016~{}\text{m}$, length
$a=3.048~{}\text{m}$ and thickness $d=0.03~{}\text{m}$ is completely clamped
at a circumferential periphery and experiences a pinching force at the center
of free-hanging periphery.

Figure 24: Pinched semi-cylindrical shell: (a) Initial configuration with the
radius $r=1.016~{}\text{m}$, length $a=3.048~{}\text{m}$ and thickness
$d=0.03~{}\text{m}$, (b-d) particles colored by the von Mises stress
$\bar{\sigma}$ of the mid-surface under 3 point force magnitudes at spatial
resolution $b/dp=80$. The material parameters are set as the density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=20.685~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3$.

Along its longitudinal edges, the vertical direction and the rotation about
the $y$-axis are constrained. The elastic material properties are density
$\rho_{0}=1100~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=20.685~{}\text{MPa}$ and Poisson’s ratio $\nu=0.3$.

Figure 24(b-d) shows the distribution of von Mises stress $\bar{\sigma}$ at
the mid-surface under varying magnitude of the point force $\bm{F}$ . Noted
that the present method features regular particle distribution and smooth
stress fields, even close to the constrained edges and place where the point
force is applied, without singularities for finite rotations (more than
$0.5\pi$). For quantitative analysis and convergence study, the downward
deflection $w_{A}$ of monitoring point $A$ as a function of the point force
magnitude and resolution is compared with that of Ref. [57]. Three different
spatial discretizations, i.e., $\pi r/dp=20,40~{}\text{and}~{}80$, are
considered for convergence study. As shown in Figure 25, the result difference
obtained by the present SPH shell method between different resolution rapidly
decreases as the spatial refinement, and the results agree well with those of
Ref. [57].

Figure 25: Pinched semi-cylindrical shell: Curves of radical displacements of
point $A$ as a function of the point force magnitude and spatial resolution,
and their comparison with those of Sze et al. [57].

##  5 Concluding remarks

In this paper, we present a reduced-dimensional SPH method for quasi-static
and dynamic analyses of both thin and moderately thick plate and shell
structures. By introducing two reduced-dimensional linear-reproducing
correction matrices, the method reproduces linear gradients of the position
and pseudo-normal. The finite deformation is taken into account by considering
all terms of strain with the help of Gauss-Legendre quadrature along the
thickness. To cope with large rotations, the method introduces weighted non-
singularity conversion relation between the rotation angles and pseudo normal.
A non-isotropic Kelvin-Voigt damping and a momentum-conserving hourglass
control algorithm with a limiter are also proposed to increase numerical
stability and to suppress hourglass modes. An extensive set of numerical
examples have been investigated to demonstrate the accuracy and robustness of
the present method. Note that, while the plate and shell structure considered
here have moderate and high modulus, one extension of the present method is
for soft thin structures, such as membranes. Another outlook, along with the
multi-physical modeling within unified computational framework, is to develop
the SPH method for the interaction between fluid and thin structures.

## CRediT authorship contribution statement

Dong Wu: Conceptualization, Methodology, Investigation, Visualization,
Validation, Formal analysis, Writing - original draft, Writing - review &
editing; Chi Zhang: Investigation, Writing - review & editing; Xiangyu Hu:
Supervision, Methodology, Investigation, Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgment

D. Wu is partially supported by the China Scholarship Council (No.
201906130189). D. Wu, C. Zhang and X.Y. Hu would like to express their
gratitude to the German Research Foundation (DFG) for its sponsorship of this
research under grant number DFG HU1527/12-4.

## References

##

  * [1] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Computer Methods in Applied Mechanics and Engineering 139 (1-4) (1996) 3–47.
  * [2] G. R. Liu, Y. T. Gu, An introduction to meshfree methods and their programming, Springer Science & Business Media, 2005.
  * [3] M. Liu, G. Liu, Smoothed particle hydrodynamics (SPH): an overview and recent developments, Archives of Computational Methods in Engineering 17 (1) (2010) 25–76.
  * [4] C. Zhang, Y.-j. Zhu, D. Wu, N. A. Adams, X. Hu, Smoothed particle hydrodynamics: Methodology development and recent achievement, Journal of Hydrodynamics 34 (5) (2022) 767–805.
  * [5] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.
  * [6] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Monthly Notices of the Royal Astronomical Society 181 (3) (1977) 375–389.
  * [7] J. J. Monaghan, Smoothed particle hydrodynamics, Reports on Progress in Physics 68 (8) (2005) 1703.
  * [8] J. J. Monaghan, Smoothed particle hydrodynamics and its diverse applications, Annual Review of Fluid Mechanics 44 (2012) 323–346.
  * [9] C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, J. Wang, X. Hu, SPHinXsys: An open-source multi-physics and multi-resolution library based on smoothed particle hydrodynamics, Computer Physics Communications (2021) 108066.
  * [10] H. G. Matthies, J. Steindorf, Partitioned strong coupling algorithms for fluid–structure interaction, Computers & structures 81 (8-11) (2003) 805–812.
  * [11] H. G. Matthies, R. Niekamp, J. Steindorf, Algorithms for strong coupling procedures, Computer Methods in Applied Mechanics and Engineering 195 (17-18) (2006) 2028–2049.
  * [12] B. C. Cerik, K. Lee, S.-J. Park, J. Choung, Simulation of ship collision and grounding damage using Hosford-Coulomb fracture model for shell elements, Ocean Engineering 173 (2019) 415–432.
  * [13] Y.-X. Peng, A.-M. Zhang, F.-R. Ming, S.-P. Wang, A meshfree framework for the numerical simulation of elasto-plasticity deformation of ship structure, Ocean Engineering 192 (2019) 106507.
  * [14] G. Totaro, Z. Gürdal, Optimal design of composite lattice shell structures for aerospace applications, Aerospace Science and Technology 13 (4-5) (2009) 157–164.
  * [15] J. D. Laubrie, J. S. Mousavi, S. Avril, A new finite-element shell model for arterial growth and remodeling after stent implantation, International Journal for Numerical Methods in Biomedical Engineering 36 (1) (2020) e3282.
  * [16] P. W. Cleary, Modelling confined multi-material heat and mass flows using SPH, Applied Mathematical Modelling 22 (12) (1998) 981–993.
  * [17] A. E. H. Love, The Small Free Vibrations and Deformation of a Thin Elastic Shell, Proceedings of the Royal Society of London Series I 43 (1887) 352–353.
  * [18] Y. S. Uflyand, Wave Propagation by Transverse Vibrations of Beams and Plates, PMM: Journal of Applied Mathematics and Mechanics 12 (1948) 287–300.
  * [19] R. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates, Journal of Applied Mechanics 18 (1951) 31–38.
  * [20] I. Elishakoff, F. Hache, N. Challamel, Vibrations of asymptotically and variationally based Uflyand–Mindlin plate models, International Journal of Engineering Science 116 (2017) 58–73.
  * [21] I. Elishakoff, Handbook on Timoshenko-Ehrenfest beam and Uflyand-Mindlin plate theories, World Scientific, 2020.
  * [22] P. Krysl, T. Belytschko, Analysis of thin shells by the element-free Galerkin method, International Journal of Solids and Structures 33 (20-22) (1996) 3057–3080.
  * [23] S. Li, W. Hao, W. K. Liu, Numerical simulations of large deformation of thin shell structures using meshfree methods, Computational Mechanics 25 (2000) 102–116.
  * [24] T. Rabczuk, P. Areias, T. Belytschko, A meshfree thin shell method for non-linear dynamic fracture, International Journal for Numerical Methods in Engineering 72 (5) (2007) 524–548.
  * [25] D. Li, Z. Lin, S. Li, Numerical analysis of Mindlin shell by meshless local Petrov-Galerkin method, Acta Mechanica Solida Sinica 21 (2) (2008) 160–169.
  * [26] B. M. Donning, W. K. Liu, Meshless methods for shear-deformable beams and plates, Computer Methods in Applied Mechanics and Engineering 152 (1-2) (1998) 47–71.
  * [27] J.-S. Chen, D. Wang, A constrained reproducing kernel particle formulation for shear deformable shell in Cartesian coordinates, International Journal for Numerical Methods in Engineering 68 (2) (2006) 151–172.
  * [28] Y. Peng, A. Zhang, F. Ming, A thick shell model based on reproducing kernel particle method and its application in geometrically nonlinear analysis, Computational Mechanics 62 (2018) 309–321.
  * [29] B. Maurel, A. Combescure, An SPH shell formulation for plasticity and fracture analysis in explicit dynamics, International Journal for Numerical Methods in Engineering 76 (7) (2008) 949–971.
  * [30] G. C. Ganzenmüller, An hourglass control algorithm for Lagrangian smooth particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 286 (2015) 87–106.
  * [31] D. Wu, C. Zhang, X. Tang, X. Hu, An essentially non-hourglass formulation for total Lagrangian smoothed particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 407 (2023) 115915.
  * [32] F.-R. Ming, A.-M. Zhang, X.-Y. Cao, A robust shell element in meshfree SPH method, Acta Mechanica Sinica 29 (2) (2013) 241–255.
  * [33] F. Caleyron, A. Combescure, V. Faucher, S. Potapov, Dynamic simulation of damage-fracture transition in smoothed particles hydrodynamics shells, International Journal for Numerical Methods in Engineering 90 (6) (2012) 707–738.
  * [34] J. Lin, H. Naceur, D. Coutellier, A. Laksimi, Efficient meshless SPH method for the numerical modeling of thick shell structures undergoing large deformations, International Journal of Non-Linear Mechanics 65 (2014) 1–13.
  * [35] F. Ming, A. Zhang, S. Wang, Smoothed particle hydrodynamics for the linear and nonlinear analyses of elastoplastic damage and fracture of shell, International Journal of Applied Mechanics 7 (02) (2015) 1550032.
  * [36] C. Zhang, Y. Zhu, Y. Yu, D. Wu, M. Rezavand, S. Shao, X. Hu, An artificial damping method for total Lagrangian SPH method with application in biomechanics, Engineering Analysis with Boundary Elements 143 (2022) 1–13.
  * [37] M. Kondo, Y. Suzuki, S. Koshizuka, Suppressing local particle oscillations in the Hamiltonian particle method for elasticity, International Journal for Numerical Methods in Engineering 81 (12) (2010) 1514–1528.
  * [38] C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, S. Zhang, J. Wang, X. Hu, SPHinXsys: An open-source meshless, multi-resolution and multi-physics library, Software Impacts 6 (2020) 100033.
  * [39] J.-L. Batoz, G. Dhatt, Modélisation des structures par éléments finis: Solides élastiques, Presses Université Laval, 1990.
  * [40] K. Wisniewski, Finite rotation shells: Basic equations and finite elements for Reissner kinematics, CIMNE-Springer, 2010.
  * [41] A. Gil, J. Segura, N. M. Temme, Numerical methods for special functions, SIAM, 2007\.
  * [42] J. N. Reddy, Theory and analysis of elastic plates and shells, CRC press, 2006.
  * [43] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. Comput. Math. 4 (1) (1995) 389–396.
  * [44] P. Randles, L. D. Libersky, Smoothed particle hydrodynamics: some recent improvements and applications, Computer Methods in Applied Mechanics and Engineering 139 (1-4) (1996) 375–408.
  * [45] R. Vignjevic, J. R. Reveles, J. Campbell, SPH in a total Lagrangian formalism, CMC-Tech Science Press- 4 (3) (2006) 181.
  * [46] P. Betsch, A. Menzel, E. Stein, On the parametrization of finite rotations in computational mechanics: A classification of concepts with application to smooth shells, Computer Methods in Applied Mechanics and Engineering 155 (3-4) (1998) 273–305.
  * [47] J. C. Simo, D. D. Fox, M. S. Rifai, On a stress resultant geometrically exact shell model. Part III: Computational aspects of the nonlinear theory, Computer Methods in Applied Mechanics and Engineering 79 (1) (1990) 21–70.
  * [48] P. Singla, D. Mortari, J. L. Junkins, How to avoid singularity for Euler angle set, in: Proceedings of the AAS Space Flight Mechanics Conference, Hawaii, 2004, pp. AAS04–190.
  * [49] T. J. Hughes, K. S. Pister, Consistent linearization in mechanics of solids and structures, Computers & Structures 8 (3-4) (1978) 391–397.
  * [50] C. Zhang, M. Rezavand, X. Hu, A multi-resolution SPH method for fluid-structure interactions, Journal of Computational Physics 429 (2021) 110028\.
  * [51] T. Y. Tsui, P. Tong, Stability of transient solution of moderately thick plate by finite-difference method, AIAA journal 9 (10) (1971) 2062–2063.
  * [52] L. D. Landau, E. M. Lifchits, Course of theoretical physics: Theory of elasticity (1986).
  * [53] J. P. Gray, J. J. Monaghan, R. Swift, SPH elastic dynamics, Computer Methods in Applied Mechanics and Engineering 190 (49-50) (2001) 6641–6662.
  * [54] C. Zhang, X. Y. Hu, N. A. Adams, A generalized transport-velocity formulation for smoothed particle hydrodynamics, Journal of Computational Physics 337 (2017) 216–232.
  * [55] B. Momenan, M. R. Labrosse, A New Continuum-Based Thick Shell Finite Element for Soft Biological Tissues in Dynamics: Part 1-Preliminary Benchmarking Using Classic Verification Experiments, arXiv preprint arXiv:1801.04029 (2018).
  * [56] K. Sze, S.-J. Zheng, A stabilized hybrid-stress solid element for geometrically nonlinear homogeneous and laminated shell analyses, Computer Methods in Applied Mechanics and Engineering 191 (17-18) (2002) 1945–1966.
  * [57] K. Sze, X. Liu, S. Lo, Popular benchmark problems for geometric nonlinear analysis of shells, Finite Elements in Analysis and Design 40 (11) (2004) 1551–1569.
  * [58] G. Payette, J. Reddy, A seven-parameter spectral/hp finite element formulation for isotropic, laminated composite and functionally graded shell structures, Computer Methods in Applied Mechanics and Engineering 278 (2014) 664–704.
  * [59] T. Belytschko, H. Stolarski, W. K. Liu, N. Carpenter, J. S. Ong, Stress projection for membrane and shear locking in shell finite elements, Computer Methods in Applied Mechanics and Engineering 51 (1-3) (1985) 221–258.
  * [60] J. Simo, D. Fox, M. Rifai, On a stress resultant geometrically exact shell model. Part II: The linear theory; computational aspects, Computer Methods in Applied Mechanics and Engineering 73 (1) (1989) 53–92.
  * [61] J. Simo, M. Rifai, D. Fox, On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the-thickness stretching, Computer Methods in Applied Mechanics and Engineering 81 (1) (1990) 91–126.
  * [62] N. Buechter, E. Ramm, Shell theory versus degeneration—a comparison in large rotation finite element analysis, International Journal for Numerical Methods in Engineering 34 (1) (1992) 39–59.
  * [63] L. Jiang, M. Chernuka, A simple four-noded corotational shell element for arbitrarily large rotations, Computers & structures 53 (5) (1994) 1123–1132.
  * [64] L. Jiang, M. W. Chernuka, N. G. Pegg, A co-rotational, updated Lagrangian formulation for geometrically nonlinear finite element analysis of shell structures, Finite Elements in Analysis and design 18 (1-3) (1994) 129–140.
  * [65] N. Stander, A. Matzenmiller, E. Ramm, An assessment of assumed strain methods in finite rotation shell analysis, Engineering Computations 6 (1) (1989) 58–66.
  * [66] B. Brank, F. Damjanić, D. Perić, On implementation of a nonlinear four node shell finite element for thin multilayered elastic shells, Computational Mechanics 16 (5) (1995) 341–359.
  * [67] R. Arciniega, J. Reddy, Tensor-based finite element formulation for geometrically nonlinear analysis of shell structures, Computer Methods in Applied Mechanics and Engineering 196 (4-6) (2007) 1048–1073.