## 1 Introduction
### Contact Constraint Method for Three-Dimensional Pavement Roughness  
Pavement roughness has long been regarded as one of the primary excitations influencing vehicle-bridge coupling vibration systems [33; 34]. Currently, pavement roughness simulation is primarily based on measured road roughness samples [35] and the power spectral density function of pavement roughness [36]. Measured roughness samples are highly specific and accurately represent the roughness of the surveyed road but are less  
Figure 1: Simplified model of the vehicle–bridge coupling system. Figure 2: Mechanical model of the vehicle–bridge coupling system. applicable to other roads. Therefore, this study simulates pavement roughness using the power spectral density function. According to the fitting expression of road surface power spectral density recommended by the Chinese standard GB/T 7031-2005 [37]:  
\[G_{q}(n)=G_{q}(n_{0})\left|\frac{n}{n_{0}}\right|^{-\omega} \tag{3}\]  
where the G\({}_{q}(n)\) is the pavement displacement power spectral density, in units of m\({}^{2}\)/m\({}^{-1}\); \(G_{q}(n_{0})\) represents the pavement roughness coefficient at spatial frequency \(n_{0}\), which depends on the pavement grade; \(n_{0}\) is the spatial reference frequency, typically set to 0.1 m\({}^{-1}\); \(w\) is the frequency exponent, determining the spectral structure of the power spectral density function; and \(n\) denotes a specific frequency within the effective spatial frequency band. This study employs the harmonic superposition method to generate pavement roughness samples, representing road roughness as the sum of numerous cosine functions with random phases. The two-dimensional pavement roughness can be expressed by Equation (4). \[r(x)=\sum_{i=1}^{N}\sqrt{2G_{q}(n_{i})\bigtriangleup n}\cos(2\pi n_{i}x+\varphi _{i}) \tag{4}\]  
where \(r(x)\) represents the two-dimensional pavement roughness; x is the longitudinal coordinate of the road surface; \(N\) is a sufficiently large integer; \(\bigtriangleup n\) is the spatial frequency interval; \(n_{i}\) represents a specific frequency within the effective frequency band of the spatial frequencies; \(\varphi_{i}\) is a random phase between \([0,2\pi]\). To better reflect the real-world pavement conditions of highway bridges, it is necessary to consider three-dimensional pavement roughness. It is assumed that the lateral pavement roughness also follows a zero-mean, ergodic Gaussian distribution [38]. The two-dimensional pavement roughness is extended laterally to obtain the three-dimensional pavement roughness, as shown in Equation (5). \[r(x,y)=\sum_{i=1}^{N}\sqrt{2G_{q}(n_{i})\bigtriangleup n}\cos\left(2\pi n_{i} \sqrt{(x^{2}+y^{2})}+\varphi_{i}(x,y)\right) \tag{5}\]  
where \(r(x,y)\) represents the three-dimensional pavement roughness, x is the longitudinal coordinate of the pavement, y is the transverse coordinate of the pavement, and the other symbols are the same as those in the previous equation. Using the harmonic superposition method, three-dimensional pavement roughness of levels A, B, and C was simulated, as shown in Figure 3. It can be seen that under the same pavement grade, the pavement roughness for the left and right wheels is rarely the same. Moreover, the higher the pavement grade, the greater the fluctuation in pavement roughness amplitude. Due to the considerable randomness of the road surface samples generated by the method detailed above, multiple pavement roughness samples are used under the same pavement grade. The average of the vehicle-bridge coupling system's multiple dynamic responses for the same pavement grade is then taken as the bridge's response value. When studying complex issues such as vehicle-bridge coupling, it is also necessary to start with certain simplified assumptions. This study assumes that the vehicle moves in a straight line and that the wheels remain in vertical contact with the bridge deck. As shown in Figure 4, suppose the vehicle travels along the positive \(X\)-axis.