## 2 Numerical Simulation
### The Basic Theory of the Vehicle-Bridge Interaction System
Assuming that the average of the weighted residual is zero, the given equation becomes  
\[\int_{0}^{L}\omega\bigg{(}\rho A\frac{\partial^{2}y}{\partial t^{2}}+EI\frac {\partial^{4}y}{\partial x^{4}}-p\bigg{)}\text{d}x=0 \tag{8}\]  
Figure 2: The adopted model and mechanical parameters of the VBI (vehicle–bridge interaction) system: The vehicle is assumed as an RBSM (rigid-body-spring model), while the bridge model is based on the Euler–Bernoulli theorem. where \(\omega\) is the weight function, and \(p\) summarizes the right-hand side terms of Equation (7). The weak form of this equation is expressed as  
\[\int_{0}^{L}\omega\rho A\frac{\partial^{2}y}{\partial t^{2}}\mathrm{d}x+\int_{0} ^{L}\frac{\partial^{2}\omega}{\partial x^{2}}EI\frac{\partial^{2}y}{\partial x^ {2}}\mathrm{d}x=\int_{0}^{L}\omega p\;\mathrm{d}x \tag{9}\]  
Let \(\mathbf{N}(x)\) and \(\mathbf{y}(t)\) be the Hermite interpolation basis functions and numerical solutions of bridge vibration displacements. The approximate solution can be given by \(y(x,t)\approx\mathbf{N}(x)\cdot\mathbf{y}(t)\). Now, assuming that the weight function is also approximated by \(\omega(x)\approx\mathbf{N}(x)\cdot\mathbf{\omega}\), this is called the Galerkin method, which is the most popular scheme for FEM. Substituting both approximated variables into Equation (9), the equation becomes  
\[\mathbf{\omega}^{\mathrm{T}}\big{(}\mathbf{M}_{b}\tilde{\mathbf{y}}(t)+\mathbf{K}_{b} \mathbf{y}(t)-\mathbf{F}_{b}(t)\big{)}=0 \tag{10}\]  
where the obtained matrices are \(\mathbf{M}_{b}=\int_{0}^{L}\mathbf{N}\mathbf{N}^{\mathrm{T}}\mathrm{d}x\) (the global mass matrix) and \(\mathbf{K}_{b}=\int_{0}^{L}\frac{\partial^{2}\mathbf{N}}{\partial x^{2}}\frac{ \partial^{2}\mathbf{N}}{\partial x^{2}}\mathrm{d}x\) (global stiffness matrix), respectively. The external force vector is  
\[ \begin{array}{rl}\mathbf{F}_{b}(t)&=\int_{0}^{L}\mathbf{N}(x)p(x,t)\mathrm{d}x\\ &=\mathbf{N}(x_{1}(t))P_{1}(t)+\mathbf{N}(x_{2}(t))P_{2}(t)+\mathbf{N}(0)R_{\mathrm{A}}(t) +\mathbf{N}(L)R_{\mathrm{B}}(t)\\ &=\mathbf{\mathrm{L}}(t)\bigg{\{}\begin{array}{l}P_{1}(t)\\ P_{2}(t)\end{array}\bigg{\}}+\mathbf{H}(t)=\mathbf{\mathrm{L}}(t)\mathbf{P}(t)+\mathbf{H}(t) \end{array}  \tag{11}\]  
where \(\mathbf{\mathrm{L}}(t)=\big{[}\mathbf{N}(x_{1}(t))\quad\mathbf{N}(x_{2}(t))\big{]}\), and \(\mathbf{H}(t)\) is the reaction force vector. When \(g\) denotes the gravity acceleration (\(g=-9.8\;\mathrm{m/s^{2}}\)), the contact forces acting on the bridge are  
\[P_{1}(t)=\frac{d_{2}m_{s}}{d_{1}+d_{2}}\big{(}g-\ddot{z}_{s1}\big{)}+m_{u1} \big{(}g-\ddot{z}_{u1}\big{)}\]  
\[P_{2}(t)=\frac{d_{1}m_{s}}{d_{1}+d_{2}}\big{(}g-\ddot{z}_{s2}\big{)}+m_{u2} \big{(}g-\ddot{z}_{u2}\big{)} \tag{12}\]  
The position of the contact force changes when the wheel position changes as the vehicle's travel. Then, assuming that Equation (10) is always zero for arbitrary \(\omega(x,t)\), and introducing the damping term, the following equation to be solved is obtained:  
\[\mathbf{M}_{b}\tilde{\mathbf{y}}(t)+\mathbf{C}_{b}\tilde{\mathbf{y}}(t)+\mathbf{K}_{ b}\mathbf{y}(t)=\mathbf{F}_{b}(t) \tag{13}\]  
where the damping is assumed to be \(\mathbf{C}_{b}=\alpha\mathbf{M}_{b}+\beta\mathbf{K}_{b}\). \(\alpha\) and \(\beta\) are called the coefficients of Rayleigh damping. The interactions of this system are represented by the input profiles \(\mathbf{u}(t)\) and the contact forces \(\mathbf{P}(t)\). The input profiles consist of the road profiles \(\mathbf{r}(t)\) and the bridge vibration components \(\widetilde{\mathbf{y}}(t)\):  
\[\mathbf{u}(t)=\mathbf{r}(t)+\widetilde{\mathbf{y}}(t) \tag{14}\]  
The road profile of \(i\)-th wheel \(r_{i}(t)\) is the vector component of \(\mathbf{r}(t)\). This variable is represented as a time function, but it is originally derived from a spatial function \(R(x)\), which denotes road unevenness.