## 2 Formulation of Novel Co-Rotational Quadrilateral Elements  
A nine-node quadrilateral shell element is presented in this section, with three coordinate systems defined at its mid-surface [13]: the global coordinate system \((X,Y,Z)\), the local coordinate system \((x,y,z)\), and the natural coordinate system \((\xi,\eta,\xi)\). In both the global and local coordinate systems, the element nodal variables \(\mathbf{u}_{G}\) and \(\mathbf{u}_{L}\) incorporate translational displacements and vectorial rotational variables, employing 45 degrees of freedom at each element:  
\[\mathbf{u}_{G}^{T}=\left\langle\mathbf{d}_{1}^{T}\ \mathbf{n}_{g1}^{T}\ \ldots\ \mathbf{d}_{i}^{T}\ \mathbf{n}_{gi}^{T}\ \ldots\ \mathbf{d}_{9}^{T}\ \mathbf{n}_{g0}^{T}\right\rangle \tag{1}\]  
\[\mathbf{u}_{L}^{T}=\left\langle\mathbf{t}_{1}^{T}\ \mathbf{\Theta}_{1}^{T}\ \ldots\ \mathbf{t}_{i}^{T}\ \mathbf{\Theta}_{i}^{T}\ \ldots\ \mathbf{t}_{9}^{T}\ \mathbf{\Theta}_{9}^{T}\right\rangle \tag{2}\]  
where \(\mathbf{d}_{i}^{T}=\left\langle U_{i}\ \ V_{i}\ \ W_{i}\right\rangle\) and \(\mathbf{t}_{i}^{T}=\left\langle u_{i}\ \ v_{i}\ \ w_{i}\right\rangle\) represent the translational displacements vector of node \(i\) in the global and local coordinate systems, respectively; \(\mathbf{n}_{gi}^{T}=\left\langle p_{i,n}\ \ p_{i,m}\right\rangle\) and \(\mathbf{\Theta}_{i}^{T}=\left\langle r_{ix}\ \ \ r_{iy}\right\rangle\) represent the vectorial rotational variables at node \(i\) of the element in the global and local coordinate systems, respectively; \(p_{i,n}\) and \(p_{i,m}\) denote the two minor components of the nodal normal vector \(\mathbf{p}_{i}\) of the shell mid-surface along the global coordinate axis; and \(r_{ix}\) and \(r_{iy}\) denote the components of \(\mathbf{p}_{i}\) along the \(x\)-axis and \(y\)-axis of the local coordinate system, respectively. The specific meanings of all symbols appearing in this paper can be found in the Nomenclature section. Figure 1 illustrates the co-rotational frame of the quadrilateral element. In the initial configuration, the local coordinate system of the element is defined with its origin at the  
diagonal intersection point, where the initial triad vectors are derived from the unit vectors of the diagonals:  
\[\left\{ \begin{array}{l}\mathbf{e}_{x0}=\frac{\mathbf{e}_{130}-\mathbf{e}_{240}} {|\mathbf{e}_{130}-\mathbf{e}_{240}|}\\ \mathbf{e}_{y0}=\frac{\mathbf{e}_{130}+\mathbf{e}_{240}}{|\mathbf{e}_{130}+ \mathbf{e}_{240}|}\\ \mathbf{e}_{z0}=\mathbf{e}_{x0}\times\mathbf{e}_{y0}\end{array} \right. \tag{3}\]  
where \(\mathbf{e}_{130}\) and \(\mathbf{e}_{240}\) are the unit vectors of the quadrilateral diagonals:  
\[\left\{ \begin{array}{l}\mathbf{e}_{130}=\frac{\mathbf{X}_{30}-\mathbf{X}_{1 0}}{|\mathbf{X}_{30}-\mathbf{X}_{10}|}\\ \mathbf{e}_{240}=\frac{\mathbf{X}_{30}-\mathbf{X}_{20}}{|\mathbf{X}_{30}- \mathbf{X}_{20}|}\end{array} \right. \tag{4}\]  
where \(\mathbf{X}_{i0}\) denotes the global coordinates of node \(i\) in the initial configuration. Upon rigid-body motion, the unit vectors of the quadrilateral diagonals in the current configuration are updated as follows:  
\[\left\{ \begin{array}{l}\mathbf{e}_{13}=\frac{\mathbf{X}_{30}-\mathbf{X}_{1 0}+\mathbf{d}_{3}-\mathbf{d}_{1}}{|\mathbf{X}_{30}-\mathbf{X}_{10}+\mathbf{d} _{3}-\mathbf{d}_{1}|}\\ \mathbf{e}_{24}=\frac{\mathbf{X}_{30}-\mathbf{X}_{20}+\mathbf{d}_{3}-\mathbf{ d}_{2}}{|\mathbf{X}_{30}-\mathbf{X}_{20}+\mathbf{d}_{4}-\mathbf{d}_{2}|}\end{array} \right. \tag{5}\]  
The triad vectors \(\left(\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\right)\) of the local coordinate axes in the current configuration are updated as follows:  
\[\left\{ \begin{array}{l}\mathbf{e}_{x}=\frac{\mathbf{e}_{13}-\mathbf{e}_{24}} {|\mathbf{e}_{13}-\mathbf{e}_{24}|}\\ \mathbf{e}_{y}=\frac{\mathbf{e}_{13}+\mathbf{e}_{24}}{|\mathbf{e}_{13}+ \mathbf{e}_{24}|}\\ \mathbf{e}_{z}=\mathbf{e}_{x}\times\mathbf{e}_{y}\end{array} \right.