## 2 Theory of Energy Conserving Integration  
### Integration Method  
The equation of motion for a nonlinear structure can be written as  
\[\mathbf{M}\mathbf{a}+\mathbf{C}\mathbf{v}+\mathbf{q}_{\text{in}}(\mathbf{u}) =\mathbf{q}_{\text{ex}} \tag{1}\]  
where \(\mathbf{M}\) and \(\mathbf{C}\) are the mass and damping matrixes, respectively; \(\mathbf{u}\), \(\mathbf{v}\), \(\mathbf{a}\), \(\mathbf{q}_{\text{in}}\) and \(\mathbf{q}_{\text{ex}}\) represent the vectors of displacement, velocity, acceleration, internal force and external load, respectively. With implicit midpoint rule [21], the discretized solution to Equation (1) at the \(n+1\) step can be expressed as  
\[\mathbf{v}_{n+1}=\mathbf{v}_{n}+\Delta t\mathbf{M}^{-1}[\mathbf{q}_{\text{ex, m}}-\mathbf{C}\mathbf{v}_{\text{m}}-\mathbf{q}_{\text{in,m}}(\mathbf{u}_{\text{m}})] \tag{2}\]  
\[\mathbf{u}_{n+1}=\mathbf{u}_{n}+\Delta t\frac{\mathbf{v}_{n+1}+\mathbf{v}_{n}} {2} \tag{3}\]  
in which \(\mathbf{v}_{\text{m}}\), \(\mathbf{u}_{\text{m}}\), and \(\mathbf{q}_{\text{ex,m}}\) are the average velocity, displacement, and external load, respectively, _i.e._,  
\[\mathbf{v}_{\text{m}}=\frac{\mathbf{v}_{n+1}+\mathbf{v}_{n}}{2} \tag{4}\]  
\[\mathbf{u}_{\text{m}}=\frac{\mathbf{u}_{n+1}+\mathbf{u}_{n}}{2} \tag{5}\]  
\[\mathbf{q}_{\text{ex,m}}=\frac{\mathbf{q}_{\text{ex,m}+1}+\mathbf{q}_{\text{ ex,m}}}{2} \tag{6}\]  
Letting  
\[\mathbf{a}_{\text{m}}=\frac{(\mathbf{v}_{n+1}+\mathbf{v}_{n})}{\Delta t} \tag{7}\]  
then Equation (2) becomes  
\[\mathbf{M}\mathbf{a}_{\text{m}}+\mathbf{C}\mathbf{v}_{\text{m}}+\mathbf{q}_{ \text{in,m}}(\mathbf{u}_{\text{m}})=\mathbf{q}_{\text{ex,m}} \tag{8}\]  
As a modified format of the implicit midpoint method, the ECM was firstly proposed for nonlinear dynamics by [PERSON] and [PERSON] [16]. [PERSON] and [PERSON] presented an implementation of energy conserving method for structures with beam elements [17]. The difference of the ECM and implicit midpoint algorithm is the calculation of internal force. With [PERSON]'s approach, the average internal force of a plane beam element is expressed as  
\[\mathbf{q}_{\text{in,m}}=\theta_{\text{a}}N_{\text{av}}\mathbf{b}_{\text{m}} +\theta_{\text{b}}\mathbf{B}_{\text{m}}^{T}\mathbf{M}_{\text{av}}+\mathbf{B} _{\text{0}}^{T}\mathbf{M}_{\text{av}} \tag{9}\]in which, \(N_{\rm av}\) and \(\mathbf{M}_{\rm av}\) are average value of axial force and end bending moments in a basic (simply supported) system, respectively; \(\mathbf{b}_{\rm m}\), \(\mathbf{B}_{\rm m}\) and \(\mathbf{B}_{0}\) are nonlinear transformation matrices of midpoint internal forces from basic to global systems (see Figure 1), and they are expressed as  
\[\mathbf{b}_{\rm m}=\frac{\mathbf{b}_{n}+\mathbf{b}_{n+1}}{2} \tag{10}\]  
\[\mathbf{b}_{n}^{T}=\left[\begin{array}{cccc}-c_{n}&-s_{n}&0&c_{n}&s_{n}&1\\ \end{array}\right] \tag{11}\]  
\[\mathbf{B}_{\rm m}=\frac{(\mathbf{B}_{n}+\mathbf{B}_{n+1})}{2} \tag{12}\]  
\[\mathbf{B}_{n}^{T}=\left[\begin{array}{cccc}-s_{n}/l_{n}&c_{n}/l_{n}&1&s_{n} /l_{n}&-c_{n}/l_{n}&0\\ -s_{n}/l_{n}&c_{n}/l_{n}&0&s_{n}/l_{n}&-c_{n}/l_{n}&1\\ \end{array}\right] \tag{13}\]  
\[\mathbf{B}_{0}^{T}=\left[\begin{array}{cccc}0&0&1&0&0&0\\ 0&0&0&0&0&1\\ \end{array}\right] \tag{14}\]  
where \(l\) and \(\Delta\alpha\) are the length and rigid-body rotation angle of the beam element, respectively; \(\rm s=sine(\Delta\alpha)\), and \(\rm c=cosine(\Delta\alpha)\). The parameters featuring energy conserving in Equation (9) are \(\rm\theta_{a}\) and \(\rm\theta_{b}\), which can be expressed by  
\[\rm\theta_{a}=\frac{2}{1+\cos\Delta\alpha} \tag{15}\]  
\[\rm\theta_{b}=\frac{\Delta\alpha}{\sin\Delta\alpha} \tag{16}\]