## 3 Solution Formulation Based on Energy-Momentum Conservation Algorithm
\begin{array}{ccc}\frac{ \partial p_{i,X}}{\partial p_{i,X}}&\frac{\partial p_{i,X}}{\partial p_{i,m}}& \frac{\partial p_{i,Z}}{\partial p_{i,m}}\\ \frac{\partial p_{i,X}}{\partial p_{i,m}}&\frac{\partial p_{i,Y}}{\partial p_{i, m}}&\frac{\partial p_{i,Z}}{\partial p_{i,m}}\\ \mathbf{A}_{2i}=\left[\begin{array}{ccc}\frac{\partial^{2}\mathbf{p}_{i}}{ \partial p_{i,2}^{2}}&\frac{\partial^{2}\mathbf{p}_{i}}{\partial p_{i,m}^{2} \partial p_{i,m}}&\frac{\partial^{2}\mathbf{p}_{i}}{\partial p_{i,m}^{2}}\\ \frac{\partial^{2}\mathbf{p}_{i}}{\partial p_{i,m}^{2}\partial p_{i,m}}&\frac{ \partial^{2}\mathbf{p}_{i}}{\partial p_{i,m}^{2}}\end{array}\right]\left\{ \begin{array}{l}\dot{p}_{i,n}\\ \dot{p}_{i,n}\end{array}\right\}\otimes\left\{\begin{array}{l}\dot{p}_{i,n} \\ \dot{p}_{i,m}\end{array}\right\}\end{array} \right. \tag{23}\]  
Substituting Equations (16)-(18) into Equation (15), extracting terms containing \(\mathbf{A}_{2i}\) to define the element equivalent load matrix \(\mathbf{F}_{m}^{*}\), and combining terms related to the global nodal variables \(\mathbf{u}_{G}\) to assemble the element equivalent mass matrix \(\mathbf{M}^{e}\), we obtain the following:  
\[\frac{dK^{e}}{dt}=\frac{d\mathbf{u}_{G}^{T}}{dt}\bigg{(}\mathbf{M}^{e}\frac{d^ {2}\mathbf{u}_{G}}{dt^{2}}+\mathbf{F}_{m}^{*}\bigg{)} \tag{24}\]  
where  
\[\mathbf{F}_{m}^{*}=\frac{\rho h^{3}}{12} \begin{Bmatrix}\int_{-1}^{1}\int_{-1}^{ 1}N_{1}\sum\limits_{j=1}^{9}N_{j}\mathbf{A}_{11}^{T}\mathbf{A}_{2j}Id\xi d \eta\\ \vdots\\ \mathbf{0}\\ \int_{-1}^{1}\int_{-1}^{1}N_{0}\sum\limits_{j=1}^{9}N_{j}\mathbf{A}_{19}^{T} \mathbf{A}_{2j}Id\xi d\eta\end{Bmatrix}  \tag{25}\]  
\[\mathbf{M}^{e}= \begin{bmatrix}\mathbf{M}_{11}^{u}&0&\dots&\mathbf{M}_{1j}^{u} &0&\dots&\mathbf{M}_{19}^{u}&0\\ 0&\mathbf{M}_{11}^{u}&\dots&0&\mathbf{M}_{1j}^{u}&\dots&0&\mathbf{M}_{19}^{u}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\ \mathbf{M}_{11}^{u}&0&\dots&\mathbf{M}_{1j}^{u}&0&\dots&\mathbf{M}_{19}^{u}&0\\ 0&\mathbf{M}_{11}^{u}&\dots&0&\mathbf{M}_{1j}^{u}&\dots&0&\mathbf{M}_{19}^{u}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots\\ \mathbf{M}_{91}^{u}&0&\dots&\mathbf{M}_{1j}^{u}&0&\dots&\mathbf{M}_{19}^{u}&0\\ 0&\mathbf{M}_{91}^{u}&\dots&0&\mathbf{M}_{1j}^{u}&\dots&0&\mathbf{M}_{19}^{u}\\ \end{bmatrix}  \tag{26}\]  
\[\mathbf{M}_{ij}^{u}=\rho h \begin{bmatrix}\int\limits_{A}N_{i}N_{j}dA&0&0\\ 0&\int\limits_{A}N_{i}N_{j}dA&0\\ 0&0&\int\limits_{A}N_{i}N_{j}dA\end{bmatrix}  \tag{27}\]  
\[\mathbf{M}_{ij}^{u}=\frac{\rho h^{3}}{12}\mathbf{A}_{1i}^{T}\int\limits_{A}N_{i} N_{j}dA\mathbf{A}_{1j} \tag{28}\]The first-order derivatives of the element strain energy and external forces with respect to time are expressed as follows:  
\[\frac{dV_{i}^{e}}{dt}=\frac{dV_{i}^{e}}{d\mathbf{u}_{G}^{T}}\frac{d\mathbf{u}_{G }}{dt}=\frac{d\mathbf{u}_{G}^{T}}{dt}\mathbf{F}_{\text{int}}^{e} \tag{29}\]  
\[\frac{dV_{ext}}{dt}=-\frac{d\mathbf{u}_{G}^{T}}{dt}\mathbf{F}_{ext} \tag{30}\]  
Substituting the kinetic energy, strain energy, and external forces of all elements in the system into Equation (12) yields the following:  
\[\frac{d\mathbf{u}_{G}^{T}}{dt}\Big{(}\mathbf{M}\frac{d^{2}\mathbf{u}_{G}}{dt^ {2}}+\mathbf{F}_{m}\Big{)}+\frac{d\mathbf{u}_{G}^{T}}{dt}\mathbf{F}_{\text{int }}=\frac{d\mathbf{u}_{G}^{T}}{dt}\mathbf{F}_{ext} \tag{31}\]  
where \(\mathbf{M}\), \(\mathbf{F}_{m}\) and \(\mathbf{F}_{\text{int}}\) are assembled from the equivalent mass matrix \(\mathbf{M}^{e}\), equivalent load matrix \(\mathbf{F}_{m}^{e}\), and internal forces \(\mathbf{F}_{\text{int}}^{e}\) of all system elements, respectively. Within time step \([t_{n},t_{n+1}]\), two generalized midpoints are formulated using the generalized energy-momentum-conserving algorithm [14]:  
\[\left\{ \begin{array}{l}t_{n+1-\alpha_{m}}=(1-\alpha_{m})t_{n+1}+\alpha_{m}t_ {n}\\ t_{n+1-\alpha_{f}}=\big{(}1-\alpha_{f}\big{)}t_{n+1}+\alpha_{f}t_{n}\end{array} \right. \tag{32}\]  
Referring to the Newmark-\(\beta\) method, generalized midpoint interpolation formulas for nodal velocities and accelerations are developed:  
\[\left\{