## 3 Theoretical Formulation
### Discrete Dynamic Equation  
In this section, according to the kinetic energy and potential energy of the system in Section 3.1 and the Chebyshev spectral approximation in Section 3.2, the discrete dynamic equation of the AFG beam with a variable cross-section on an inhomogeneous foundation subjected to a non-uniformly distributed moving load is derived. The displacement fields \(w(x,t)\) and \(\alpha(x,t)\) can be approximated using Chebyshev expansions:  
\[w(x,t)=\sum_{k=0}^{N-1}a_{1k}(t)\mathcal{T}_{k}(x),\quad\alpha(x,t)=\sum_{k=0}^{ N-1}a_{2k}(t)\mathcal{T}_{k}(x) \tag{17}\]  
and be discretized by sampling points  
\[x_{k}=\frac{L}{2}-\frac{L}{2}\cos[k\pi/(N-1)],\quad w_{k}=w(x_{k},t),\quad \alpha_{k}=\alpha(x_{k},t). \tag{18}\]  
where \(k=0,1,2,\ldots,N-1\). Define  
\[ \begin{split}& h_{1}(x)=\rho(x)A(x),\quad h_{2}(x)=\rho(x)I(x), \quad h_{3}(x)=E(x)I(x),\\ & h_{4}(x)=\kappa G(x)A(x),\quad h_{5}(x)=k_{f}(x),\quad h_{6}=c_ {f}(x)\end{split}  \tag{19}\]  
Denote the vector \(\{w_{k}\}\) as \(\mathbf{w},\{\alpha_{k}\}\) as \(\mathbf{\pi}\). According to Equations (2) and (16), the kinetic energy of the system can be computed as  
\[T=\frac{1}{2}\dot{\mathbf{w}}^{\mathsf{T}}\mathbf{V}_{h_{1}}\dot{\mathbf{w}}+\frac{1}{2} \dot{\mathbf{\pi}}^{\mathsf{T}}\mathbf{V}_{h_{2}}\dot{\mathbf{\kappa}} \tag{20}\]  
where \(\dot{\mathbf{w}}\) and \(\dot{\mathbf{\kappa}}\) are the sampled vectors of \(\dot{w}\) and \(\dot{\mathbf{\kappa}}\) at the Gauss-Lobatto points, respectively. According to Equations (3) and (16), the potential energy of the beam can be computed as  
\[U_{b}=\frac{1}{2}\mathbf{\pi}^{\mathsf{T}}\Big{(}\mathbf{Q}_{1}^{\mathsf{T}}\mathbf{V}_{h_ {3}}\mathbf{Q}_{1}\Big{)}\mathbf{\alpha}+\frac{1}{2}\mathbf{w}^{\mathsf{T}}\Big{(}\mathbf{Q}_ {1}^{\mathsf{T}}\mathbf{V}_{h_{4}}\mathbf{Q}_{1}\Big{)}\mathbf{w}-\mathbf{w}^{\mathsf{T}}\Big{(} \mathbf{Q}_{1}^{\mathsf{T}}\mathbf{V}_{h_{4}}\mathbf{\alpha}+\frac{1}{2}\mathbf{\pi}^{\mathsf{ T}}\mathbf{V}_{h_{4}}\mathbf{\alpha} \tag{21}\]  
Similarly, the formula for calculating the elastic potential energy of the foundation and the Rayleigh dissipation function of the viscoelastic foundation reaction can be obtained:  
\[U_{f}=\frac{1}{2}\mathbf{w}^{\mathsf{T}}\mathbf{V}_{h_{6}}\mathbf{w},\quad D=\frac{1}{2} \dot{\mathbf{w}}^{\mathsf{T}}\mathbf{V}_{h_{6}}\dot{\mathbf{w}} \tag{22}\]  
where \(\mathbf{V}_{h_{i}}\) is the weighted inner product matrix about the weight function \(h_{i},i=1,2,\ldots,6\), \(\mathbf{Q}_{1}\) is the first derivative matrix, as shown in Equation (14). To calculate the work done by the non-uniformly distributed moving load \(q\), define a moving frame \(\partial-\dot{x}\dot{y}\), as shown in Figure 2. The range of the moving load is indicated by the red dotted line, which is denoted as \(\bar{w}(\dot{x},t)\). Appling the Chebyshev spectral approximation and the Gauss-Lobatto sampling on the interval \([0,\bar{l}]\), the work can be computed by  
\[W_{q}=\int_{0}^{\bar{l}}q(\dot{x},t)\bar{w}(\dot{x},t)\mathrm{d}\bar{x}=\bar{ \mathbf{q}}^{\mathsf{T}}\bar{\mathbf{V}}\bar{\mathbf{w}} \tag{23}\]  
where \(\bar{\mathbf{V}}\) is the Chebyshev inner product matrix defined in Equation (15) and where the superscript \(\tilde{\ }\) is used to indicate its integral domain \([0,\bar{l}]\). \(\bar{\mathbf{q}}\) and \(\bar{\mathbf{w}}\), respectively, are the Gauss-Lobatto sampling vectors of functions \(q(\dot{x},t)\) and \(\bar{w}(\dot{x},t)\) on the interval \([0,\bar{l}]\) in the coordinate system \(\partial-\dot{x}\dot{y}\). Denote the coordinates of these sampling points in the coordinate system \(o-xy\) as  
\[\{x_{0},x_{1},x_{2},\ldots,x_{N-1}\} \tag{24}\]  
The value of elements of \(\bar{\mathbf{w}}\) can be calculated using Equation (17)  
\[\bar{\mathbf{w}}= \begin{bmatrix}w(x_{0},t)\\ w(x_{1},t)\\ w(x_{2},t)\\ \vdots\\ w(x_{N-1},t)\end{bmatrix}