## 2 Methodology
### Finite Element Modeling
\begin{array}{cc}\dot{u}_{\Then the corresponding matrix form of the beam motion equation can be written as follows:  
\[\mathbf{M}_{b}\ddot{\mathbf{u}}_{b}(t)+\mathbf{C}_{b}\dot{\mathbf{u}}_{b}(t)+ \mathbf{K}_{b}\mathbf{u}_{b}(t)=\mathbf{P}_{b}(t)=\mathbf{H}\mathbf{P}(t) \tag{26}\]  
By using Hermitian cubic interpolation as the shape function, the equivalent nodal load vector \(\mathbf{P}_{b}(t)\) could be given. The superscript \(T\) means the matrix transpose operation.  
\[\mathbf{P}_{b}(t)=\begin{bmatrix}0&\dots&N_{2}&\dots&0&\dots&0\\ 0&\dots&0&\dots&N_{1}&\dots&0\end{bmatrix}^{T}\begin{bmatrix}p_{1}(t)\\ p_{2}(t)\end{bmatrix}=\mathbf{H}\mathbf{P}(t) \tag{27}\]  
When the front axle or the rear axle is not on the beam, the corresponding shape function \(N_{1}\) or \(N_{2}\) is a null vector. Otherwise, this shape function and its derivative will take the following form:  
\[\begin{array}{l}N_{j}=\left[\begin{array}{c}\frac{2s^{3}}{L_{c}^{3}}- \frac{3\sigma^{2}}{L_{c}^{2}}+1&\frac{s^{3}}{L_{c}^{2}}-\frac{2s^{2}}{L_{c}}+s &\frac{3\sigma^{2}}{L_{c}^{2}}-\frac{2s^{3}}{L_{c}^{3}}&\frac{s^{3}}{L_{c}^{2} }-\frac{s^{2}}{L_{c}}\end{array}\right]\\ N^{\prime}_{j}=\left[\begin{array}{c}\frac{6s^{3}}{L_{c}^{3}}-\frac{6s}{L_{ c}^{2}}&\frac{3\sigma^{2}}{L_{c}^{2}}-\frac{4s}{L_{c}}+1&\frac{6s}{L_{c}^{2}}- \frac{6s^{2}}{L_{c}^{3}}&\frac{3\sigma^{2}}{L_{c}^{2}}-\frac{2s}{L_{c}}\\ s=x_{c}(t)-(n-1)L_{c},(n-1)L_{c}\leq x_{c}(t)\leq nL_{c},j=1,2\end{array}\right]\end{array} \tag{28}\]  
where \(x_{c}(t)\) is the contact point location with the front or rear axle during the nth beam element and \(L_{c}\) is the element length.  
Therefore, the combination of Equations (21) and (26) could describe this interaction system with the vehicle and the beam. The Newmark-Beta integration [34] is used to calculate the responses of the vehicle as well as the beam.  
\[\begin{array}{l}\left[\begin{array}{cc}\mathbf{M}_{c}&0\\ 0&\mathbf{M}_{b}\end{array}\right]\left[\begin{array}{c}\ddot{\mathbf{u}}_{ b}(t)\\ \ddot{\mathbf{u}}_{b}(t)\end{array}\right]+\left[\begin{array}{cc}\mathbf{C}_{c }&-\mathbf{C}_{t1}\mathbf{H}^{T}\\ -\mathbf{H}\mathbf{C}_{t2}&\mathbf{C}_{b}+\mathbf{H}\mathbf{C}_{t3}\mathbf{H} ^{T}\end{array}\right]\left[\begin{array}{c}\dot{\mathbf{u}}_{b}(t)\\ \dot{\mathbf{u}}_{b}(t)\end{array}\right]+\left[\begin{array}{cc}\mathbf{K}_ {c}&-\mathbf{C}_{t1}\mathbf{V}_{e}\mathbf{H}^{T}-\mathbf{K}_{t1}\mathbf{H}^{T }\\ -\mathbf{H}\mathbf{K}_{t2}&\mathbf{K}_{b}+\mathbf{H}\mathbf{C}_{t3}V_{e}\mathbf{H }^{T}+\mathbf{H}\mathbf{K}_{t3}\mathbf{H}^{T}\end{array}\right]\left[ \begin{array}{c}\mathbf{u}_{b}(t)\\ \mathbf{u}_{b}(t)\end{array}\right]\\ -\mathbf{H}\mathbf{P}_{b}-\mathbf{H}\mathbf{C}_{t3}\dot{\mathbf{u}}_{b}(t)- \mathbf{H}\mathbf{K}_{t3}\dot{\mathbf{u}}_{b}(t)\end{array}\right]\end{array}  \tag{29}\]  
The contact points of the vehicles are assumed to be always touching on the surface of the beam and the road. Their responses could be derived from the corresponding time-varied-location beam response \(u_{bc}(x(t),\,t)\) and the surface roughness \(u_{r}(x(t))\). The responses (displacements, velocities, and accelerations) of the three trails' contact points can be expressed as follows:  
\[\begin{array}{l}u_{c,k}(t)=\left[ \begin{array}{ccccc}&0&\dots&N_{k}&\dots& 0\\ \dot{u}_{c,k}(t)=\left[\begin{array}{ccccc}&0&\dots&N_{k}&\dots&0\\ 0&\dots&N_{k}&\dots&0\\ \end{array}\right]\dot{\mathbf{u}}_{b}(t)+\dot{u}_{r}(x_{c,k}(t))=\dot{u}_{bc} (V_{e}(t-t_{k1}),t)+\dot{u}_{r}(V_{e}(t-t_{k1}))\\ \ddot{u}_{c,k}(t)=\left[\begin{array}{ccccc}&0&\dots&N_{k}&\dots&0\\ \end{array}\right]\ddot{\mathbf{u}}_{b}(t)+\ddot{u}_{r}(x_{c,k}(t))=\ddot{u}_{ bc}(V_{e}(t-t_{k1}),t)+\ddot{u}_{r}(V_{e}(t-t_{k1}))\\ k=1,2,3\end{array}  \tag{30}\]