## 2 Formulation of Novel Co-Rotational Quadrilateral Elements
\tag{6}\]  
For a quadrilateral shell element, the normal vector at node \(i\) in the initial configuration is computed as follows:  
\[\overline{\mathbf{p}}_{i0}=\frac{\partial\left[\sum\limits_{i=1}^{9}N_{i}( \xi,\eta)\mathbf{X}_{i0}\right]}{\partial\xi}\times\left.\frac{\partial\left[ \sum\limits_{i=1}^{9}N_{i}(\xi,\eta)\mathbf{X}_{i0}\right]}{\partial\eta} \right|_{(\xi_{i},\eta_{i})} \tag{7}\]  
where \(N_{i}(\xi,\eta)\) denotes the Lagrangian interpolation functions of node \(i\). The normal vectors of adjacent elements at shared nodes are averaged to ensure geometric continuity across element boundaries. Figure 1: Co-rotational framework for quadrilateral element. Coordinates of any point on the element neutral surface can be interpolated by means of shape functions from the element node coordinates:  
\[\mathbf{x}=\left\{ \begin{array}{l}x\\ y\\ z\end{array} \right\}=\sum\limits_{i=1}^{m}N_{i}\mathbf{x}_{i0} \tag{8}\]  
Similarly, the dynamic variables of any point in the element under the local coordinate system can also be obtained by interpolation:  
\[\mathbf{t}=\left\{ \begin{array}{l}u\\ v\\ w\end{array} \right\}=\left\{ \begin{array}{l}\sum\limits_{i=1}^{9}N_{i}(\xi, \eta)u_{i}\\ \sum\limits_{i=1}^{9}N_{i}(\xi,\eta)v_{i}\\ \sum\limits_{i=1}^{9}N_{i}(\xi,\eta)w_{i}\end{array} \right\} \tag{9}\]  
\[\left\{ \begin{array}{l}r_{x}\\ r_{y}\end{array} \right\}=\left\{ \begin{array}{l}\sum\limits_{i=1}^{9}N_{i}( \xi,\eta)r_{ix}\\ \sum\limits_{i=1}^{9}N_{i}(\xi,\eta)r_{iy}\end{array} \right\} \tag{10}\]  
Therefore, the displacement and coordinate functions of any point P in the element can be obtained:  
\[\left\{ \begin{array}{l}\mathbf{t}_{p}=\left\{\begin{array}{l}u_{P}\\ v_{P}\\ w_{P}\end{array}\right\}=\sum\limits_{i=1}^{9}N_{i}(\xi,\eta)\left\{\begin{array} []{l}u_{i}\\ v_{i}\\ w_{i}\\ x_{i0}\\ y_{i0}\\ z_{i0}\end{array}\right\}+\frac{1}{2}\zeta h\sum\limits_{i=1}^{9}N_{i}(\xi, \eta)\left\{\begin{array}{l}r_{ix}-r_{ix0}\\ r_{iy}-r_{iy0}\\ r_{iz}-r_{iz0}\\ r_{iz0}\end{array}\right\}\end{array} \right. \tag{11}\]  
where \(\left\{x_{i0}\quad y_{i0}\quad z_{i0}\right\}\) is the local coordinate value of node i in the initial configuration; and \(\left\{r_{ix0}\quad r_{iy0}\quad r_{iz0}\right\}\) represents the three components of the normal vector of node i in the local coordinate system under the initial configuration. The Green-Lagrange strains are adopted to describe the elastic deformation of the shell element, partitioned into three components: membrane strain vector \(\varepsilon_{m}\), shear strain vector \(\gamma\), and bending strain vector \(z_{l}\chi\). \[\varepsilon=\varepsilon_{m}+z_{l}\chi=\varepsilon_{m}+\frac{1}{2}\zeta h\chi \tag{12}\]  
where \(\varepsilon\) represents the membrane strain components of the element; \(h\) denotes the thickness of the shell element; and \(\zeta\) is the coordinate along the thickness direction in the natural coordinate system. The potential energy functional of the elemental is expressed as follows:  
\[\Pi=\frac{1}{2}\underset{y}{\int}\left(\varepsilon_{m}+z_{l}\chi\right)^{T} \mathbf{D}_{1}(\varepsilon_{m}+z_{l}\chi)dV+\frac{1}{2}\underset{V}{\int} \mathbf{\gamma}^{T}\mathbf{D}_{2}\gamma dV-W_{ext} \tag{13}\]  
where \(\mathbf{D}_{1}\) and \(\mathbf{D}_{2}\) are the elastic-moduli matrices, and \(W_{ext}\) represents the external forces.