## 3 Solution Formulation Based on Energy-Momentum Conservation Algorithm
\begin{array}{l}\ddot{\mathbf{u}}\Big{(}t_{n+1-\alpha_{f}}\Big{)}= \frac{1-\alpha_{f}}{\beta\Delta t}[\mathbf{u}(t_{n+1})-\mathbf{u}(t_{n})]- \frac{(1-\alpha_{f})\gamma-\beta}{\beta}\ddot{\mathbf{u}}(t_{n})-\frac{( \gamma-2\beta)(1-\alpha_{f})}{2\beta}\Delta t\ddot{\mathbf{u}}(t_{n})\\ \ddot{\mathbf{u}}(t_{n+1-\alpha_{m}})=\frac{1-\alpha_{m}}{\beta\Delta t^{2}}[ \mathbf{u}(t_{n+1})-\mathbf{u}(t_{n})]-\frac{1-\alpha_{m}}{\beta\Delta t} \ddot{\mathbf{u}}(t_{n})-\frac{1-\alpha_{m}-2\beta}{2\beta}\ddot{\mathbf{u}}(t _{n})\end{array} \right. \tag{33}\]  
where \(\alpha_{m}\), \(\alpha_{f}\), \(\beta\) and \(\gamma\) are integration parameters governed by the spectral radius \(\rho_{\infty}\). Applying the variational operations to both sides of the motion equilibrium equation (Equation (26)) and substituting the formulations of velocities and accelerations at generalized midpoint instant (Equation (28)), all known quantities at time \(t_{n}\) are extracted to yield the equivalent force vector \(\mathbf{F}_{a}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\):  
\[\mathbf{F}_{a}\Big{(}t_{n+1-\alpha_{f}}\Big{)}=\mathbf{M}\Big{(}t_{n+1-\alpha _{f}}\Big{)}\bigg{[}\frac{(1-\alpha_{m})}{\beta\Delta t^{2}}\mathbf{u}_{G}(t _{n})+\frac{(1-\alpha_{m})}{\beta\Delta t}\dot{\mathbf{u}}_{G}(t_{n})+\frac{( 1-\alpha_{m}-2\beta)}{2\beta}\ddot{\mathbf{u}}_{G}(t_{n})\bigg{]} \tag{34}\]  
This yields a system of governing differential equations for motion equilibrium, with \(\mathbf{u}_{G}(t_{n+1})\) as the primary unknown:  
\[\mathbf{M}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\frac{(1-\alpha_{m})}{\beta\Delta t ^{2}}\mathbf{u}_{G}(t_{n+1})+\mathbf{F}_{\text{int}}\Big{(}t_{n+1-\alpha_{f}} \Big{)}=\mathbf{F}_{ext}\Big{(}t_{n+1-\alpha_{f}}\Big{)}+\mathbf{F}_{a}\Big{(} t_{n+1-\alpha_{f}}\Big{)}-\mathbf{F}_{m}\Big{(}t_{n+1-\alpha_{f}}\Big{)} \tag{35}\]  
where the internal force vector \(\mathbf{F}_{\text{int}}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\) and tangent stiffness matrix \(\mathbf{k}_{G}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\) at the generalized midpoint configuration in the global coordinate system are interpolated from quantities at time steps \(t_{n}\) and \(t_{n+1}\):  
\[ \begin{array}{l}\mathbf{F}_{\text{int}}\Big{(}t_{n+1-\alpha_{f}}\Big{)}= \mathbf{T}^{T}\mathbf{f}_{\text{int}}\Big{(}t_{n+1-\alpha_{f}}\Big{)}=\mathbf{T} ^{T}h\int_{\mathcal{V}}\Big{[}\mathbf{B}_{m}^{T}\Big{(}t_{n+1-\alpha_{f}} \Big{)}\mathbf{D}_{1}\mathbf{\varepsilon}_{m}\Big{(}t_{n+1-\alpha_{f}}\Big{)} \\ \quad+\frac{h^{2}}{12}\mathbf{B}_{b}^{T}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\mathbf{D }_{1}\mathbf{x}\Big{(}t_{n+1-\alpha_{f}}\Big{)}+\mathbf{B}_{\gamma}^{T}\Big{(} t_{n+1-\alpha_{f}}\Big{)}\mathbf{D}_{2}\gamma\Big{(}t_{n+1-\alpha_{f}}\Big{)}\Big{]}dV \end{array}  \tag{36}\]  
Then, at the (k + 1)-th iteration, the equilibrium equation is updated to  
\[\overline{\textbf{K}}_{T}^{k}\Big{(}t_{n+1-\alpha_{f}}\Big{)}\Delta\textbf{u}^{k +1}(t_{n+1})=\overline{\textbf{P}}^{k}\Big{(}t_{n+1-\alpha_{f}}\Big{)} \tag{44}\]  
After each iteration is completed, the equivalent stiffness matrix and equivalent load vector need to be updated based on the displacement changes. After the completion of the (K + 1)-th iteration step, the work done by the residual load is \(E^{k}\), and the ratio of \(E^{k}\) to the residual load work \(E\) in the first iteration step are calculated as follows:  
\[cc=\frac{E^{k}}{E}=\frac{\Delta\textbf{u}^{k+1}(t_{n+1})\overline{\textbf{P}} ^{k}\Big{(}t_{n+1-\alpha_{f}}\Big{)}}{\Delta\textbf{u}^{1}(t_{n+1})\overline{ \textbf{P}}^{1}\Big{(}t_{n+1-\alpha_{f}}\Big{)}} \tag{45}\]  
When \(cc\leq 1\times 10^{-6}\) and \(k\geq 3\), the calculation converges, terminates the iteration, and proceeds to the iteration process of the next time step; otherwise, we set \(k=k+1\) to continue the iteration. For clarity, the complete iterative calculation flow chart is presented in Figure 2.  
Figure 2: Calculation flow chart.