## 4 Solution of the System of Equilibrium Equations
### Numerical Solution of Equations of Motion  
The system of Equations (12) was solved using a numerical method and a discretization with respect to time. For this purpose, the direct differential method with respect to time was applied. In this method, acceleration \(\ddot{\boldsymbol{q}}=\left(\ddot{u}_{i},\ddot{w}_{i},\ddot{\varphi}_{i1} \right)^{T}\) and displacement velocity \(\dot{\boldsymbol{q}}=\left(\dot{u}_{i},\ddot{w}_{i},\dot{\varphi}_{i1}\right) ^{T}\) were approximated in the system of Equations (13) at time instants \(t_{n-1}=t_{n}-\Delta t\), \(t_{n}=n\Delta t\), and \(t_{n+1}=t_{n}+\Delta t\), using Equation (23). \[\ddot{\boldsymbol{q}}^{n}=\frac{\dot{\boldsymbol{q}}^{n}-\dot{\boldsymbol{q}}^ {n-1}}{\Delta t}\text{ and }\dot{\boldsymbol{q}}^{n}=\frac{\Delta\boldsymbol{q}^{n}}{\Delta t} \tag{23}\]  
Using these approximations, the system of Equation (13) becomes the following. \[G(\boldsymbol{q},\lambda)=\Delta\boldsymbol{q}^{n}-\Delta\boldsymbol{q}_{G}^{ n}-\Delta\lambda^{n}\Delta\boldsymbol{q}_{P}^{n}=0 \tag{24}\]  
where \(\Delta\boldsymbol{q}^{n}=\left(\Delta u_{i}^{n},\Delta w_{i}^{n},\Delta \varphi_{i1}^{n}\right)^{T}\) is a vector of searched displacement increments, \(\Delta\lambda^{n}=\lambda^{n}-\lambda^{n-1}\) and is increment of the load parameter in the current time step. The vector of components of the displacement increments for the total load has the following form. \[\Delta\boldsymbol{q}_{P}^{n}=\left\{ \begin{array}{rl}a_{i}P_{x}(s_{i})\\ a_{i}P_{z}(s_{i})\\ a_{fi}[M_{i_{k1}}k(x_{i}-x_{i_{k1}})+M_{i_{k1}}k(x_{i}-x_{i_{k1}})]\end{array} \right. \tag{25}\]  
where \(k(f)=\left\{ \begin{array}{rl}1&if\quad f=0\\ 0&if\quad f\
eq 0\end{array} \right.\) is a selection operator and \(i\) is the node number. The vector of components of the displacement increments for the load parameter \(\lambda^{n-1}\) from the previous time step has a form. \[\Delta\boldsymbol{q}_{G}^{n}=\left\{ \begin{array}{rl}b_{i}\Delta u_{i}^{n-1 }+a_{i}[-N_{i1}\cos\theta_{i1}+N_{i0}\cos\theta_{i0}+Q_{i1}\sin\theta_{i1}-Q_ {i0}\sin\theta_{i0}]\\ b_{i}\Delta w_{i}^{n-1}+a_{i}[N_{i1}\sin\theta_{i1}-N_{i0}\sin\theta_{i0}+Q_{i1 }\cos\theta_{i1}-Q_{i0}\cos\theta_{i0}]\\ b_{fi}\Delta\varphi_{i1}^{n}+a_{fi}(M_{i+2}-M_{i1}-Q_{i1}\Delta s_{i1})\end{array}  \right\}+\lambda^{n-1}\Delta\boldsymbol{q}_{P}^{n} \tag{26}\]  
where \(b_{i}=\frac{m_{i}}{m_{i}+C_{i}\Delta t}\) is a coefficient for linear displacements, \(a_{i}=\frac{\Delta t^{2}}{m_{i}+C_{i}\Delta t}\) is a coefficient for internal forces, \(b_{fi}=\frac{\dot{i}}{\dot{i}+C_{fi}\Delta t}\) is a coefficient for rotations, and \(a_{fi}=\frac{\Delta t^{2}}{\dot{i}+C_{fi}\Delta t}\) is a coefficient for bending moments. The initial conditions for time step \(t_{n=0}=0\) have the following form. \[\Delta\boldsymbol{q}^{n=-1}=0 \tag{27}\]  
The iterative procedure terminates when the solutions of displacements in subsequent time instants have converged, according to the condition below. \[\|\Delta\boldsymbol{q}^{n}\|\leq\varepsilon_{\Delta q} \tag{28}\]  
The time step \(\Delta t\) is determined so that the numerical integration is stable and has the following form.