## Theoretical Modeling
### Mid-Point Algorithm  
The dynamic equilibrium equation of the mid-point algorithm is given by:  
\[\mathbf{g}_{m}=\mathbf{q}_{e,m}-\mathbf{q}_{i,m}-\mathbf{q}_{mas,m}=\mathbf{0} \tag{11}\]  
where \(\mathbf{q}_{e,m}\), \(\mathbf{q}_{i,m}\) and \(\mathbf{q}_{mas,m}\) are vectors of mid-point external, internal and inertial nodal forces, respectively.  
The mid-point inertial force vector in present framework is defined as (note it is different from the one given in [37]):  
\[\mathbf{q}_{mas,m}=\frac{1}{\Delta t}(\mathbf{U}_{e,n+1}^{*}\mathbf{M}\mathbf{ U}_{0}^{*}\dot{\mathbf{p}}_{n+1}-\mathbf{U}_{e,n}^{*}\mathbf{M}\mathbf{U}_{0}^{*} \dot{\mathbf{p}}_{n}) \tag{12}\]  
where \(\dot{\mathbf{p}}_{n}\) and \(\dot{\mathbf{p}}_{n+1}\) are, respectively, the vectors of nodal velocities at step \(n\) and \(n+1\), whose translational components are defined in the global reference system \(\mathbf{U}_{g}\) and the rotational ones in the body-attached frame \(\mathbf{U}_{e}\).  
The nodal external and internal forces acting during the time step are represented by their mid-point values:  
\[\mathbf{q}_{e,m} =\frac{\mathbf{q}_{e,n}+\mathbf{q}_{e,n+1}}{2} \tag{13}\] \[\mathbf{q}_{i,m} =(\frac{\mathbf{T}_{n}+\mathbf{T}_{n+1}}{2})^{T}\frac{\mathbf{q} _{il,n}+\mathbf{q}_{il,n+1}}{2} \tag{14}\]  
where \(\mathbf{q}_{il,n}\) and \(\mathbf{q}_{il,n+1}\) are the local nodal internal forces at time \(n\) and \(n+1\), respectively, and the matrices \(\mathbf{T}_{n}\) and \(\mathbf{T}_{n+1}\) are the corresponding tangential transformation matrices.  
The inertial contribution to the effective stiffness matrix is obtained from the variation of Equation (12):  
\[\delta\mathbf{q}_{mas,m}=\mathbf{K}_{trans,m}\delta\mathbf{p}=(\mathbf{K}_{ mas1,m}+\mathbf{K}_{mas2,m})\delta\mathbf{p} \tag{15}\]  
The full expressions of the matrices \(\mathbf{K}_{mas1,m}\) and \(\mathbf{K}_{mas2,m}\) are given as:  
\[\mathbf{K}_{mas1,m}=\left[\begin{array}{cccc}\mathbf{0}&\mathbf{0}&\mathbf{ 0}\\ \mathbf{0}&-\mathbf{S}(\mathbf{U}_{1,n+1}\mathbf{M}_{r}(2\mathbf{w}_{1}+ \mathbf{w}_{2}))&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&-\mathbf{S}(\mathbf{U}_{2,n+1}\mathbf{M}_{r }(\mathbf{w}_{1}+2\mathbf{w}_{2}))\end{array}\right] \tag{16}\]  
\[\mathbf{K}_{mas2,m}=\frac{2}{\Delta t}\mathbf{U}_{e,n+1}^{*}\mathbf{M} \mathbf{D} \tag{17}\]  
where \(\mathbf{M}_{r}=\frac{l}{6}\mathbf{I}_{p}\).  
The mid-point static tangential stiffness matrix \(\mathbf{K}_{{{\it{}}stat,m}}\) is given by:  
\[\delta\mathbf{q}_{i,m}=\mathbf{K}_{{{\it{}}stat,m}}\delta\mathbf{p}=\frac{( \mathbf{T}_{n}+\mathbf{T}_{n+1})^{T}}{2}\frac{\mathbf{K}_{l}}{2}\mathbf{T}_{n+ 1}\delta\mathbf{p}+\frac{1}{2}\mathbf{K}_{t\sigma}(\frac{\mathbf{q}_{il,n}+ \mathbf{q}_{il,n+1}}{2})\delta\mathbf{p} \tag{18}\]  
which comes from the variation of Equation (14). \(\mathbf{K}_{l}\) is the local linear stiffness matrix, and the geometric stiffness matrix \(\mathbf{K}_{t\sigma}\) takes the identical form to that of the statics [36].  
The effective tangential stiffness matrix of the mid-point algorithm is then expressed as the sum of the static term \(\mathbf{K}_{{{\it{}}stat,m}}\) and the dynamic term \(\mathbf{K}_{{{\it{}}tmas,m}}\):  
\[\mathbf{K}_{t,m}=\mathbf{K}_{{{\it{}}stat,m}}+\mathbf{K}_{{{\it{}}tmas,m}} \tag{19}\]In some cases, the approximate conservation of energy of the mid-point algorithm is not sufficient to ensure the numerical stability, which undergoes an uncontrolled growth of energy and finally is unable to reach convergence. A numerical damping correction is introduced as applied to the mid-point algorithm, in which, in place of Equation (14), the mid-point internal force \(\mathbf{q}_{i,m}\) is redefined as:  
\[\mathbf{q}_{i,m}=(\frac{\mathbf{T}_{n}+\mathbf{T}_{n+1}}{2})^{T}(\frac{ \mathbf{q}_{il,n+1}+\mathbf{q}_{il,n}}{2}+\xi\frac{\mathbf{q}_{il,n+1}-\mathbf{ q}_{il,n}}{2}) \tag{20}\]  
where positive parameter \(\xi\) controls the numerical damping.