## 2 The Vehicle-Pavement-Foundation System: Modeling and Solution
### Represent the Pavent with the E-B Beam Model  
The differential equation for the uniform motion of an infinitely long E-B beam on a viscoelastic Winkler foundation under the action of tire force is given by the following equation:  
\[EI\frac{\partial^{4}w}{\partial x^{4}}+m\frac{\partial^{2}w}{\partial t^{2}}+c \frac{\partial w}{\partial t}+kw+k_{N}w^{3}=F(t)\,\delta(x-vt), \tag{10}\]  
where \(w\) represents the vertical displacement of the beam, \(E\) is the elastic modulus of the beam, \(I\) is the moment of inertia of the beam's cross-section, \(m\) is the mass per unit length of the beam, and \(c\) is the damping coefficient of the foundation; \(k\) and \(k_{N}\) are the linear stiffness and nonlinear stiffness of the foundation, respectively. \(v\) signifies the velocity of the vehicle, \(t\) indicates time and the Dirac function is utilized to define the vehicle's location. The initial conditions and boundary conditions for pavement vibration are set as  
\[\begin{split} w(x,t)\Big{|}_{t=0}=\frac{\partial w(x,t)}{ \partial t}\Big{|}_{t=0}=0,\\ \lim_{x\to\pm\infty}\frac{\partial^{n}w(x,t)}{\partial x_{n}}=0 \quad(n=0,1,2,3).\end{split} \tag{11}\]  
Next, the Adomian decomposition method is employed to solve Equation (10). This method has been proven to provide a solution in a rapidly convergent series often with less computational work than traditional approaches [8; 38]. Assume that the solution of Equation (10) can be expressed as  
\[w=\sum_{j=0}^{\infty}w_{j}. \tag{12}\]  
The nonlinear term in Equation (10) is expanded into an infinite series form  
\[k_{N}w^{3}=k_{N}\sum_{n=0}^{\infty}B_{n}, \tag{13}\]  
where the sequence \(B_{n}\) represents the Adomian polynomials, which are obtained by the following formula  
\[B_{n}=\frac{1}{n!}\Bigg{[}\frac{d^{n}}{d\lambda^{n}}\Bigg{(}\sum_{k=0}^{ \infty}\lambda^{k}w_{k}\Bigg{)}^{3}\Bigg{]}_{\downarrow_{\lambda=0}}. \tag{14}\]  
The first several Adomian polynomials are explicitly expressed as  
\[\begin{split} B_{0}&=w_{0}^{3},\quad B_{1}=3w_{0}^{ 2}w_{1},\\ B_{2}&=3(w_{0}w_{1}^{2}+w_{0}^{2}w_{2}),\quad B_{3} =w_{1}^{3}+6w_{0}w_{1}w_{2}+w_{0}^{2}w_{3},\\ B_{4}&=3(w_{1}^{2}+w_{2}+w_{0}w_{2}^{2}+2w_{0}w_{1} w_{3}+w_{0}^{2}w_{4}),\\ &\ldots\end{split} \tag{15}\]  
Define the following differential operator  
\[L[\cdot]=EI\frac{\partial^{4}}{\partial x^{4}}+m\frac{\partial^{2}}{\partial t ^{2}}+c\frac{\partial}{\partial t}+k. \tag{16}\]Substituting Equations (12), (13) and (15) into Equation (10), yields  
\[\begin{split} L[w_{0}]&=F(t)\delta(x-v\,t),\\ L[w_{1}]&=-k_{N}B_{0},\\ L[w_{2}]&=-k_{N}B_{1},\\ &\cdots\end{split} \tag{17}\]  
The terms \(w_{j}\) can be recursively calculated from Equation (17). Define  
\[\phi_{j}=\left\{\begin{array}{ll}\frac{||w_{j+1}||}{||w_{j}||},&||w_{j}||\
eq 0,\\ 0,&||w_{j}||=0,\end{array}\right. \tag{18}\]  
where \(||w_{j}||=max\{|w_{j}|\}\). It has been proven [44] that decomposition series (12) rapidly converges to the exact solution if \(0\leq\phi_{j}<1\), \(j=0,1,2,\cdots\). The approximate solution obtained by superimposing the first several terms of \(w_{j}\) can be written as  
\[w\approx w_{0}+w_{1}+w_{2}+\cdots, \tag{19}\]  
where  
\[\begin{split} w_{0}&=u_{1}(e^{i\omega t}+e^{-i\omega t })+u_{2},\\ w_{1}&=-k_{N}[u_{3}+u_{4}(e^{i\omega t}+e^{-i\omega t })+u_{5}(e^{i2\omega t}+e^{-i2\omega t})+u_{6}(e^{i3\omega t}+e^{-i3\omega t}) ],\\ w_{2}&=3k_{N}^{2}[u_{7}+u_{8}(e^{i\omega t}+e^{-i \omega t})+u_{9}(e^{i2\omega t}+e^{-i2\omega t})+u_{10}(e^{i3\omega t}+e^{-i3 \omega t})\\ &\quad+u_{11}(e^{i4\omega t}+e^{-i4\omega t})+u_{12}(e^{i5\omega t }+e^{-i5\omega t})]\end{split}\]  
The expressions for \(u_{i}\), \(i=1,2,\cdots,12\), and the corresponding solution process can be found in Appendix A. Based on the form of the solution, it is clear to see that, compared to the linear stiffness case, the presence of nonlinear stiffness of the foundation introduces higher harmonic terms in the solution.