## 3 Numerical Solution
### Time Integration  
The time integration of the equations of motion of the whole system is performed as proposed in [32], where the \(\theta\) method [43] was used for integrating the equations of motion. First, the equations of motion for the block (8)-(10) are rewritten in matrix form:  
\[\mathbf{M}_{k}\vec{q}_{k}(t)=\mathbf{F}_{k}(t)+\boldsymbol{r}_{k}(t), \tag{16}\]  
where  
\[\mathbf{M}_{k}= \begin{bmatrix}m_{k}&0&0\\ 0&m_{k}&0\\ 0&0&I_{k}\end{bmatrix} ,\quad\vec{\boldsymbol{q}}_{k}(t)= \begin{bmatrix} \tilde{X}_{k}(t)\\ \tilde{Y}_{k}(t)\\ \tilde{\varphi}_{k}(t)\end{bmatrix} ,\quad\mathbf{F}_{k}(t)= \begin{bmatrix} F_{x}(t)\\ F_{y}(t)\\ M_{k}(t)\end{bmatrix} ,\quad\mathbf{r}_{k}(t)= \begin{bmatrix}r_{x,k}(t)\\ r_{y,k}(t)\\ M_{r}(t)\end{bmatrix} . \tag{17}\]  
The motion of the block is described by the position vector \(q_{k}(t)=\left[X_{k}(t),Y_{k}(t),\varphi_{k}(t)\right]^{T}\) and its derivatives with respect to time are denoted by the dots over the vector. In the equations above, \(\mathbf{M}_{k}\) is the block mass matrix, \(\vec{\boldsymbol{q}}_{k}(t)\) is the acceleration vector, \(\mathbf{F}_{k}(t)\) is the vector of external loads, and \(\boldsymbol{r}_{k}(t)\) is the vector of the contact force. The index \(k\) indicates a reference to the block. We consider some arbitrary time interval \((t_{i},t_{i+1}]\) of length \(\Delta t\). Integrating (16) over \((t_{i},t_{i+1}]\), the velocity at the end of the time interval follows as  
\[\dot{\mathbf{q}}_{k}(t_{i+1})=\dot{\mathbf{q}}_{k}(t_{i})+\mathbf{M}_{k}^{-1}\int_{t_{i }}^{t_{i+1}}\mathbf{F}_{k}(t)\,dt+\mathbf{M}_{k}^{-1}\int_{t_{i}}^{t_{i+1}}\mathbf{r}_{ k}(t)\,dt. \tag{18}\]  
The position at \(t_{i+1}\) is evaluated as  
\[\mathbf{q}_{k}(t_{i+1})=\mathbf{q}_{k}(t_{i})+\int_{t_{i}}^{t_{i+1}}\dot{\mathbf{q}}_{k}(t )\,dt. \tag{19}\]  
Hereafter, approximations of all quantities at time instants \(t_{i}\) and \(t_{i+1}\) are referred to through the indices \(i\), and \(i+1\), e.g., \(\mathbf{q}_{k}(t_{i})\approx\mathbf{q}_{k,i}\), \(\mathbf{q}_{k}(t_{i+1})\approx\mathbf{q}_{k,i+1}\), and so on. The contact force variation with time is unknown; therefore, only the impulse of the force is calculated. More precisely, in contact dynamics, the mean value of the impulse is sought:  
\[\mathbf{r}_{k,i+1}=\frac{1}{\Delta t}\int_{t_{i}}^{t_{i+1}}\mathbf{r}_{k}(t)\,dt, \tag{20}\]  
which basically corresponds to a constant contact force over the interval \((t_{i},t_{i+1}]\), and therefore, we refer to it as the contact force. Applying the \(\theta\) method in (18) and (19) results in discretized expressions for the velocity and position at \(t_{i+1}\)  
\[\dot{\mathbf{q}}_{k,i+1}=\dot{\mathbf{q}}_{k,i}+\mathbf{M}_{k}^{-1}[\Delta t\theta\mathbf{ F}_{k,i+1}+\Delta t(1-\theta)\mathbf{F}_{k,i}]+\mathbf{M}_{k}^{-1}\Delta t\mathbf{r}_{k,i+1}, \tag{21}\]  
\[\mathbf{q}_{k,i+1}=\mathbf{q}_{k,i}+\Delta t\theta\dot{\mathbf{q}}_{k,i+1}+\Delta t(1- \theta)\dot{\mathbf{q}}_{k,i}, \tag{22}\]  
where \(\theta\) is some parameter between \(0\) and \(1\). With (21) and (22), an approximate discrete solution for the block position is obtained. For the beam, we repeat the procedure. We integrate (15) and approximate the integrals with the \(\theta\) method, which gives  
\[ \begin{split}\mathbf{M}_{e}\dot{\mathbf{q}}_{e,i+1}=\mathbf{M}_{e} \dot{\mathbf{q}}_{e,i}&+\Delta t\theta(-\mathbf{C}\dot{\mathbf{q}}_{e,i +1}-\mathbf{K}\mathbf{q}_{e,i+1}+\mathbf{p}_{i+1})\\ &+\Delta t(1-\theta)(-\mathbf{C}\dot{\mathbf{q}}_{e,i}-\mathbf{K}\bm {q}_{e,i}+\mathbf{p}_{i})+\Delta t\mathbf{r}_{e,i+1}.\end{split}  \tag{23}\]  
Now, after approximating the positions at \(t_{i+1}\), and grouping all the terms with \(\dot{\mathbf{q}}_{e,i+1}\) on the left side of the equation, the approximate velocities and positions at the end of the time-step are  
\[