沥青路面轴对称动力响应问题中解析解对比

Comparative of analytical solutions for axial symmetric dynamic response of asphalt pavements

摘要

当前有关路面力学的研究多集中于解析解推导以及影响因素探究上，缺乏不同解析解方法的差异性研究。为此，本文分别以传递矩阵法、刚度矩阵法、波传递法以及反射、透射矩阵法推导出了考虑材料横观各向同性以及层间接触状态的轴对称路面结构动力响应解析解。在此基础上，结合4种路面结构研究了上述4种解析解的计算效率以及适用范围。结果表明：①在矩阵规模上，由小到大依次为反射、透射矩阵法、传递矩阵法、刚度矩阵法以及波传递法；②在求解速度上，传递矩阵法耗时显著少于其余3种解析解方法，且反射、透射矩阵法耗时最长；③在数值稳定性上，传递矩阵法在进行有限厚度的动力响应计算时极易出现数值溢出问题，而其余解析解方法均具有较好的数值稳定性；④在实际应用上，传递矩阵法尤其适用于求解半空间体表面弯沉；刚度矩阵法仅可计算各结构层表面的动力响应；而波传递法和反射、透射矩阵法可求解路面内部任意位置的动力响应。

Abstract

The current researches on pavement mechanics usually focus on the derivation of analytical solutions and the exploration of the factors affecting them， there are few studies on the difference of analytical solutions. In order to solve this problem， transfer matrix method， stiffness matrix method， wave transfer method and reflection and transmission matrix method were applied in the deduction process of the analytical solutions considering the transverse isotropy of materials and the contact state between layers. On this basis， the calculation efficiency and the range of adaptability of the above four analytical solutions are studied by combining four kinds of pavement structures. The results show that： ①for the matrix scale， reflection and transmission matrix method， transfer matrix method， stiffness matrix method and wave transfer method increase in turn； ② for the solution speed， the transfer matrix method takes much less time than the other three methods， while the reflection and transmission matrix method consumes the most time； ③in the numerical stability， the transfer matrix method is very easy to overflow data when calculating the dynamic response of finite thickness， while the other analytical solution methods have good numerical stability； ④in practical application， the transfer matrix method is particularly suitable for the surface of semi space problems， and the stiffness matrix method can only calculate the dynamic response of the surface of each structural layer， however， the wave transfer method and the reflection and transmission matrix method can solve the dynamic response of any position in the pavement.

Graphical abstract

关键词

层状体系力学 / 传递矩阵 / 刚度矩阵 / 波传递 / 反射、透射矩阵

Key words

mechanics of the layered system / transfer matrix method / stiffness matrix method / wave transfer method / reflection and transmission matrix method

层状体系力学是路面力学的主要研究方向，其出现对路面结构的设计、施工及检测具有重要意义^［1］。 随着对路面力学认识的不断加深，路面结构中所存在的黏弹性^［2-5］、横观各向同性^［5-8］甚至复杂层间条件^［9-11］等因素被相继考虑。陆建飞等^［12］在考虑材料横观各向同性及饱和土等因素下，在频域内对埋置动荷载展开了分析。马宪永等^［13］在考虑路面结构黏弹性、横观各向同性及复杂层间接触的基础上，开展了表面移动荷载的研究。张军辉等^［14］在考虑横观各向同性及层间非完全连续条件的基础上，结合BP神经网络与实数编码的遗传算法开展了路面结构参数反演可行性研究。Li等^［15］在考虑环境温度、材料黏弹性及非线性特征的有限元模型基础上，应用遗传算法——神经网络优化算法开展了路面结构参数反演。

尽管在路面力学中考虑的因素越来越多，但在解析解的求解方法上，只有系数递推法^［16］、传递矩阵法^{［14，17，18］}、刚度矩阵法^{［5，11，19，20］}、精细积分法^{［1，21-23］}、波传递法^［13］以及反射与透射矩阵法^{［12，24］}。这些方法中，由Thomson^［25］以及Haskell^［26］提出的传递矩阵法是最早用于计算路面力学响应的方法，并广泛应用于层状静力体系^［17］。但由于其矩阵元素中同时含有正指数和负指数部分，因而易发生数据溢出现象，这在动力计算时尤为突出。因此，后续提出的各种解析解方法均将消除数据溢出置于核心位置。其中，刚度矩阵法通过构建单层刚度矩阵避免了数据溢出问题，并在此基础上组装整体刚度矩阵以求解力学响应^［27］。但随着结构层层数的增加，刚度矩阵规模不断增大，最终影响计算效率。精细积分法则通过将路面结构划分为一系列厚度极小的亚层，并用级数近似表示正指数来消除数据溢出问题，但过多的亚层往往会显著降低计算效率^［28］。波传递法与反射、透射矩阵法相似，均在传递矩阵法的基础上建立了各结构层层顶下行波及层底上行波表达式，进而建立了波传播矩阵方程以求解力学响应^［29］。由于这两种方法在波传递矩阵计算时仅包含负指数部分，因而均具有良好的数值稳定性。

总体上看，目前的研究主要集中在解析解推导及影响因素探究方面，而鲜有不同解析解间差异性对比的研究。针对这一问题，本文以考虑材料横观各向同性及层间接触状态的轴对称路面结构为例，开展了不同解析解方法的对比研究。首先，基于Hankel-Laplace积分变换及现代控制理论，得到了考虑横观各向同性的动力响应通解表达式。进而，在考虑层间条件的基础上，依次应用传递矩阵法、刚度矩阵法、波传递法以及反射、透射矩阵法推导了Hankel-Laplace域的动力解析解，并结合Durbin法以及高斯插值积分法获得了时域动力响应。最后，以4种路面结构为例对上述不同解析解方法进行对比研究。研究结果可为路面力学解析解的选择与改进提供参考。

１ 考虑横观各向同性的弹性层状体系

如图1所示，将路面结构简化为一轴对称N层层状体系结构，其每一层层底的坐标分别为z₁，z₂，…，z_N，考虑材料横观各向同性以及各结构层层间的不完全连续，表面作用轴对称荷载p（r，t）。本节在轴对称动力平衡方程、几何方程以及考虑横观各向同性的物理方程基础上，结合Hankel-Laplace变换得到了考虑横观各向同性的轴对称动力体系常微分方程组，为后续传递矩阵法、刚度矩阵法、波传递法以及反射与透射矩阵法求解奠定基础。

1.1　动力基本方程

柱坐标系下动力平衡方程^［16］可以表示为：

\left\{\begin{array}{l}\frac{\partial {\sigma }_r}{\partial r}+\frac{\partial {\tau }_{zr}}{\partial z}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=\rho \frac{{\partial }^{\mathrm{2}}{u}_r}{\partial t^{\mathrm{2}}}\\ \frac{\partial {\sigma }_z}{\partial z}+\frac{\partial {\tau }_{zr}}{\partial r}+\frac{{\tau }_{zr}}{r}=\rho \frac{{\partial }^{\mathrm{2}}{u}_z}{\partial t^{\mathrm{2}}}\end{array}\right.（1）

式中：σ_r 、{\sigma }_{\theta }、σ_z 、{\tau }_{zr}分别为径向、环向、竖向和切向应力分量；r和z分别为距离坐标原点的水平和垂直距离；ρ为密度；u_r 和u_z 分别为水平以及垂直方向位移；t为时间。

在柱坐标系下考虑横观各向同性的物理方程可表示为：

\left[\begin{array}{c}{\sigma }_r\\ {\sigma }_{\theta }\\ {\sigma }_z\\ {\tau }_{zr}\end{array}\right]=\left[\begin{array}{cccc}{C}_{\mathrm{11}}& C_{\mathrm{12}}& C_{\mathrm{13}}& \mathrm{0}\\ C_{\mathrm{12}}& C_{\mathrm{11}}& C_{\mathrm{13}}& \mathrm{0}\\ C_{\mathrm{13}}& C_{\mathrm{13}}& C_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& C_{\mathrm{44}}\end{array}\right]\left[\begin{array}{c}{\epsilon }_r\\ {\epsilon }_{\theta }\\ {\epsilon }_z\\ {\gamma }_{zr}\end{array}\right]（2）

式中：C_{\mathrm{11}}=kn(\mathrm{1}-n{\mu }_{\mathrm{v}}^{\mathrm{2}})；C_{\mathrm{12}}=kn({\mu }_{\mathrm{h}}+n{\mu }_{\mathrm{v}}^{\mathrm{2}})；C_{\mathrm{13}}=kn{\mu }_{\mathrm{v}}(\mathrm{1}+{\mu }_{\mathrm{v}})；C_{\mathrm{33}}=k(\mathrm{1}-{\mu }_{\mathrm{h}}^{\mathrm{2}})；C_{\mathrm{44}}=\mathrm{0.5}{E}_{\mathrm{v}}(\mathrm{1}+{\mu }_{\mathrm{v}}{)}^{-\mathrm{1}}；n为模量比，n=E_{\mathrm{h}}{E}_{\mathrm{v}}^{-\mathrm{1}}；k=E_{\mathrm{v}}\cdot [(\mathrm{1}+{\mu }_{\mathrm{h}})(\mathrm{1}-{\mu }_{\mathrm{h}}-\mathrm{2}n{\mu }_v^{\mathrm{2}}{)]}^{-\mathrm{1}}；ε_r 、{\epsilon }_{\theta }、ε_z、γ_zr分别为径向、环向、竖向和切向应变分量；μ_v和μ_h分别为垂直和水平方向的泊松比；E_h和E_v分别为水平和垂直方向上的模量。柱坐标下的几何方程为：

\left[\begin{array}{c}{\epsilon }_r\\ {\epsilon }_{\theta }\\ {\epsilon }_z\\ {\gamma }_{zr}\end{array}\right]={\left[\begin{array}{cccc}\frac{\partial }{\partial r}& \frac{\mathrm{1}}{r}& \mathrm{0}& \frac{\partial }{\partial z}\\ \mathrm{0}& \mathrm{0}& \frac{\partial }{\partial z}& \frac{\partial }{\partial r}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}{u}_r\\ u_z\end{array}\right]（3）

1.2　求解问题构建

基于状态空间法，选取u_r 、u_z 、σ_z 、τ_zr 为状态变量。将物理方程（式（2））以及几何方程（式（3））代入动力平衡方程（式（1）），消去变量{\sigma }_{\theta }、{\sigma }_r，可得：

\left\{\begin{array}{l}{C}_{\mathrm{11}}\left(\frac{{\partial }^{\mathrm{2}}}{\partial r^{\mathrm{2}}}+\frac{\partial }{r\partial r}-\frac{\mathrm{1}}{r^{\mathrm{2}}}-\rho \frac{{\partial }^{\mathrm{2}}}{\partial t^{\mathrm{2}}}\right)u_r+\\ C_{\mathrm{13}}\frac{{\partial }^{\mathrm{2}}{u}_z}{\partial r\partial z}+\frac{\partial {\tau }_{zr}}{\partial z}=\mathrm{0}\\ \rho \frac{{\partial }^{\mathrm{2}}{u}_z}{\partial t^{\mathrm{2}}}-\left(\frac{\partial }{\partial r}+\frac{\mathrm{1}}{r}\right){\tau }_{zr}-\frac{\partial {\sigma }_z}{\partial z}=\mathrm{0}\\ C_{\mathrm{44}}\frac{\partial u_r}{\partial z}+C_{\mathrm{44}}\frac{\partial u_z}{\partial r}-{\tau }_{zr}=\mathrm{0}\\ C_{\mathrm{13}}\left(\frac{\partial }{\partial r}+\frac{\mathrm{1}}{r}\right)u_r+C_{\mathrm{33}}\frac{\partial u_z}{\partial z}-{\sigma }_z=\mathrm{0}\end{array}\right.（4）

对t执行Laplace变换，

\overline{f}(s)={\int }_{\mathrm{0}}^{+\infty }f(t){\mathrm{e}}^{-st}\mathrm{d}t

对r执行Hankel变换^［30］，

{\tilde{f}}_{\nu }(\xi )={\int }_{\mathrm{0}}^{+\infty }rf(r)J_{\nu }(\xi r)\mathrm{d}r

J_{\nu }(\cdot )为第ν阶的Bessel函数。由以上可得：

\frac{\mathrm{d}\widetilde{\overline{\mathit{S}}}}{\mathrm{d}z}={\mathit{A}}_j(\xi ,s)\widetilde{\overline{\mathit{S}}}（5）

式中：\widetilde{\overline{\mathit{S}}}为Hankel-Laplace域内状态向量，\widetilde{\overline{\mathit{S}}}={\left[\begin{array}{cccc}{\tilde{\overline{u}}}_{r\mathrm{1}}(z)& {\tilde{\overline{u}}}_{z\mathrm{0}}(z)& {\tilde{\overline{\tau }}}_{zr\mathrm{1}}(z)& {\tilde{\overline{\sigma }}}_{z\mathrm{0}}(z)\end{array}\right]}^{\mathrm{T}}，下标0和1分别代表第0阶和第1阶的Hankel变换。{\mathit{A}}_j(\xi ,s)=\left[\begin{array}{cccc}\mathrm{0}& \xi & C_{\mathrm{44}}^{-\mathrm{1}}& \mathrm{0}\\ -C_{\mathrm{13}}{C}_{\mathrm{33}}^{-\mathrm{1}}\xi & \mathrm{0}& \mathrm{0}& C_{\mathrm{33}}^{-\mathrm{1}}\\ C_{\mathrm{11}}{\xi }^{\mathrm{2}}-{\xi }^{\mathrm{2}}{C}_{\mathrm{13}}^{\mathrm{2}}{C}_{\mathrm{33}}^{-\mathrm{1}}+\rho s^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \xi C_{\mathrm{13}}{C}_{\mathrm{33}}^{-\mathrm{1}}\\ \mathrm{0}& \rho s^{\mathrm{2}}& -\xi & \mathrm{0}\end{array}\right]（6）

1.3　边界条件

如图1所示，在Hankel-Laplace域内，荷载边界条件为{\tilde{\overline{\sigma }}}_{z\mathrm{0}}{|}_{z=\mathrm{0}}=-\tilde{\overline{p}}(\xi ,s)，{\tilde{\overline{\tau }}}_{zr\mathrm{1}}{|}_{z=\mathrm{0}}=\mathrm{0}；其中，\tilde{\overline{p}}(\xi ,s)为Hankel-Laplace变换后的表面荷载条件，\tilde{\overline{p}}(\xi ,s)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{\int }_{\mathrm{0}}^{+\mathrm{\infty }}p(r,t){\mathrm{e}}^{-st}{J}_{\mathrm{0}}(\xi r)\mathrm{d}r\mathrm{d}t；位移边界条件为{\tilde{\overline{u}}}_{z\mathrm{0}}\left|{}_{z=z_N}\right.=\mathrm{0}，{\tilde{\overline{u}}}_{r\mathrm{1}}\left|{}_{z=z_N}\right.=\mathrm{0}，当z_N\to +\mathrm{\infty }时便退化为半空间问题，即\widetilde{\overline{\mathit{S}}}(+\mathrm{\infty })\left|{}_N\right.={\mathrm{0}}_{\mathrm{4}\times \mathrm{1}}。此外，由于实际道路各结构层材料间所存在差异性，不可避免会出现层间部分甚至完全滑动的情况。因而，采用Matsui等^［31］提出的剪切弹簧模拟路面层间接触状态。

(\mathrm{1}-{\alpha }_x^n)\left[{\tilde{\overline{u}}}_{r\mathrm{1}}^{n+\mathrm{1}}(z_n)-{\tilde{\overline{u}}}_{r\mathrm{1}}^n(z_n)\right]={\alpha }_x^n{\beta }_x^n{\tilde{\overline{\tau }}}_{zr\mathrm{1}}^n(z_n)（7）

式中：上标n表示结构层层位；{\alpha }_x^n为第n层与第（n+1）层的层间滑移参数，；其中，当{\alpha }_x^n=\mathrm{0}时，为完全连续状态，当{\alpha }_x^n=\mathrm{0.999}时，表示完全滑动状态；{\beta }_x^n为与n层与（n+1）层材料有关的系数，与模量和泊松比有关，按{\beta }_x^n=b[(\mathrm{1}+{\mu }_{\mathrm{v}n}^{})E_{\mathrm{v}n}^{-\mathrm{1}}+(\mathrm{1}+{\mu }_{\mathrm{v}(n+\mathrm{1})}^{})E_{\mathrm{v}(n+\mathrm{1})}^{-\mathrm{1}})]计算。

进一步的，可获得层间边界条件：

{\mathit{F}}_n\cdot \widetilde{\overline{\mathit{S}}}\left|{}_n\right.(z_n)=\widetilde{\overline{\mathit{S}}}\left|{}_{n+\mathrm{1}}\right.(z_n)

式中：{\mathit{F}}_n=\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}& {\tilde{\mathit{F}}}_n\\ {\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}& {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]；{\tilde{\mathit{F}}}_n=\left[\begin{array}{cc}\frac{{\alpha }_x^n{\beta }_x^n}{\mathrm{1}-{\alpha }_x^n}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right]，{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}为2阶单位矩阵，{\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}为2阶零矩阵；假设第N层与刚性基岩处于完全连续状态，即F_N={\mathit{I}}_{\mathrm{4}\times \mathrm{4}}。

2 解析解推导

2.1　传递矩阵法

在式（5）所示的常微分方程组中，由于 A_j （ξ，s）特征方程为λ⁴+κ₁λ²+κ₂形式，因而必定存在两对互为相反数的特征值^［32］。根据现代控制理论，\widetilde{\overline{\mathit{S}}}(z)存在如下通解：

\widetilde{\overline{\mathit{S}}}(z)=\mathit{P}{\mathrm{e}}^{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(-a,-b,a,b)z}{\mathit{P}}^{-\mathrm{1}}\widetilde{\overline{\mathit{S}}}(\mathrm{0})（8）

式中：-a、-b、a、b为矩阵 A_j （ξ，s）特征值； P 为矩阵 A_j （ξ，s）特征向量。代入层间接触条件，可获得相邻层间的传递矩阵：

\widetilde{\overline{\mathit{S}}}(z_n)\left|{}_{n+\mathrm{1}}\right.={\mathit{F}}_n\cdot {\mathit{\Xi }}_n(z_n-z_{n-\mathrm{1}})\widetilde{\overline{\mathit{S}}}(z_{n-\mathrm{1}})\left|{}_n\right.（9）

在传递矩阵（式（9））的基础上，底部状态向量和顶部状态向量间可建立如下联系：

\widetilde{\overline{\mathit{S}}}(z_N)\left|{}_N\right.={\mathrm{\Xi }}_N(z_N-z_{N-\mathrm{1}}){\displaystyle \coprod _{n=\mathrm{1}}^{N-\mathrm{1}}}{\mathit{F}}_n{\mathrm{\Xi }}_n(z_n-z_{n-\mathrm{1}})\widetilde{\overline{\mathit{S}}}(\mathrm{0})\left|{}_{\mathrm{1}}\right.（10）

令\overline{\mathit{\Xi }}={\mathrm{\Xi }}_N(z_N-z_{N-\mathrm{1}}){\displaystyle \coprod _{n=\mathrm{1}}^{N-\mathrm{1}}}{\mathit{F}}_n{\mathrm{\Xi }}_n(z_{n+\mathrm{1}}-z_n)，代入前述应力和位移边界条件，可得：

\left\{\begin{array}{l}{\tilde{\overline{u}}}_{z\mathrm{0}}^{\mathrm{1}}(\mathrm{0})=\frac{{\overline{\mathrm{\Xi }}}_{\mathrm{14}}/{\overline{\mathrm{\Xi }}}_{\mathrm{11}}-{\overline{\mathrm{\Xi }}}_{\mathrm{24}}/{\overline{\mathrm{\Xi }}}_{\mathrm{21}}}{{\overline{\mathrm{\Xi }}}_{\mathrm{12}}/{\overline{\mathrm{\Xi }}}_{\mathrm{11}}-{\overline{\mathrm{\Xi }}}_{\mathrm{22}}/{\overline{\mathrm{\Xi }}}_{\mathrm{21}}}\tilde{\overline{p}}(\xi ,s)\\ {\tilde{\overline{u}}}_{r\mathrm{1}}^{\mathrm{1}}(\mathrm{0})=\frac{{\overline{\mathrm{\Xi }}}_{\mathrm{14}}/{\overline{\mathrm{\Xi }}}_{\mathrm{12}}-{\overline{\mathrm{\Xi }}}_{\mathrm{24}}/{\overline{\mathrm{\Xi }}}_{\mathrm{22}}}{{\overline{\mathrm{\Xi }}}_{\mathrm{11}}/{\overline{\mathrm{\Xi }}}_{\mathrm{12}}-{\overline{\mathrm{\Xi }}}_{\mathrm{21}}/{\overline{\mathrm{\Xi }}}_{\mathrm{22}}}\tilde{\overline{p}}(\xi ,s)\end{array}\right.（11）

式中：{\overline{\mathrm{\Xi }}}_{ij}为矩阵\overline{\mathit{\Xi }}第i行第j列的矩阵元素。

若z_N\to +\mathrm{\infty }时，有\widetilde{\overline{\mathit{S}}}(+\mathrm{\infty })\left|{}_N\right.={\mathrm{0}}_{\mathrm{4}\times \mathrm{1}}，必定有{\mathit{P}}_N^{-\mathrm{1}}\widetilde{\overline{\mathit{S}}}(z_{N-\mathrm{1}})\left|{}_N\right.=[x_{\mathrm{1}N},{x_{\mathrm{2}}}_N{,\mathrm{0,0}]}^{\mathrm{T}}，x_{\mathrm{1}N}和{x_{\mathrm{2}}}_N为与边界条件相关的待定系数。结合式（10），有\widetilde{\overline{\mathit{S}}}(z_{N-\mathrm{1}})\left|{}_N\right.={\displaystyle \coprod _{n=\mathrm{1}}^{N-\mathrm{1}}}{\mathit{F}}_n{\Xi }_n(z_n-z_{n-\mathrm{1}})\widetilde{\overline{\mathit{S}}}(\mathrm{0})\left|{}_{\mathrm{1}}\right.。为便于表示，令\tilde{\mathit{\Xi }}={\mathit{P}}_N^{-\mathrm{1}}{\displaystyle \coprod _{n=\mathrm{1}}^{N-\mathrm{1}}}{\mathit{F}}_n{\Xi }_n(z_{n+\mathrm{1}}-z_n)，类似有：

\left\{\begin{array}{l}{\tilde{\overline{u}}}_{z\mathrm{0}}^{\mathrm{1}}(\mathrm{0})=\frac{{\tilde{\mathrm{\Xi }}}_{\mathrm{34}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{31}}-{\tilde{\mathrm{\Xi }}}_{\mathrm{44}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{41}}}{{\tilde{\mathrm{\Xi }}}_{\mathrm{32}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{31}}-{\tilde{\mathrm{\Xi }}}_{\mathrm{42}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{41}}}\tilde{\overline{p}}(\xi ,s)\\ {\tilde{\overline{u}}}_{r\mathrm{1}}^{\mathrm{1}}(\mathrm{0})=\frac{{\tilde{\mathrm{\Xi }}}_{\mathrm{34}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{32}}-{\tilde{\mathrm{\Xi }}}_{\mathrm{44}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{42}}}{{\tilde{\mathrm{\Xi }}}_{\mathrm{31}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{32}}-{\tilde{\mathrm{\Xi }}}_{\mathrm{41}}/{\tilde{\mathrm{\Xi }}}_{\mathrm{42}}}\tilde{\overline{p}}(\xi ,s)\end{array}\right.（12）

传递矩阵法思路清晰，矩阵规模小，计算速度快。但在实际计算中，由于传递矩阵\mathit{\Xi }(z)同时包含了正指数和负指数，反复相乘很容易出现数据溢出问题，因此仅能较为稳定地计算表面力学响应。为提高其数值稳定性，在每次计算传递矩阵后，按本次计算的最小元素对传递矩阵执行整体缩小处理^［33］。

2.2　刚度矩阵法

令{\mathrm{e}}^{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(-az,-bz,az,bz)}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{\Lambda }}^{-\mathrm{1}}(z),\mathit{\Lambda }(z))，\widetilde{\overline{\mathit{S}}}=[{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}},{\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}{]}^{\mathrm{T}}，其中：\mathit{\Lambda }(z)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathrm{e}}^{az},{\mathrm{e}}^{bz})，{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}=[{\tilde{\overline{u}}}_{r\mathrm{1}}^{}{\tilde{\overline{u}}}_{z\mathrm{0}}^{}{]}^{\mathrm{T}},{\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}=[{\tilde{\overline{\tau }}}_{zr\mathrm{1}}{\tilde{\overline{\sigma }}}_{z\mathrm{0}}{]}^{\mathrm{T}}。在式（8）基础上，考虑层间条件有：

\left[\begin{array}{c}{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_n)\left|{}_{n+\mathrm{1}}\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}(z_n)\left|{}_{n+\mathrm{1}}\right.\end{array}\right]={\mathit{F}}_n\left[\begin{array}{cc}{\mathit{P}}_{\mathrm{1}n}& {\mathit{P}}_{\mathrm{2}n}\\ {\mathit{P}}_{\mathrm{3}n}& {\mathit{P}}_{\mathrm{4}n}\end{array}\right]\times

\left[\begin{array}{cc}{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}})& \\ & {\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}})\end{array}\right]\times

\left[\begin{array}{cc}{\mathit{P}}_{\mathrm{1}n}^{-\mathrm{1}}& {\mathit{P}}_{\mathrm{2}n}^{-\mathrm{1}}\\ {\mathit{P}}_{\mathrm{3}n}^{-\mathrm{1}}& {\mathit{P}}_{\mathrm{4}n}^{-\mathrm{1}}\end{array}\right]\left[\begin{array}{c}{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{n-\mathrm{1}})\left|{}_n\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}(z_{n-\mathrm{1}})\left|{}_n\right.\end{array}\right]（13）

式中：{\mathit{P}}_{\mathrm{1}n}~{\mathit{P}}_{\mathrm{4}n}，{\mathit{P}}_{\mathrm{1}n}^{-\mathrm{1}}~{\mathit{P}}_{\mathrm{4}n}^{-\mathrm{1}}分别为{\mathit{P}}_n以及{\mathit{P}}_n^{-\mathrm{1}}的2阶子矩阵。

构建如式（14）所示的单元刚度矩阵。

\left[\begin{array}{c}{\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}(z_{n-\mathrm{1}})\left|{}_n\right.\\ -{\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}(z_n)\left|{}_{n+\mathrm{1}}\right.\end{array}\right]=\left[\begin{array}{cc}{\overline{\mathit{K}}}_{\mathrm{1}}^n& {\overline{\mathit{K}}}_{\mathrm{2}}^n\\ -{\overline{\mathit{K}}}_{\mathrm{3}}^n& -{\overline{\mathit{K}}}_{\mathrm{4}}^n\end{array}\right]\left[\begin{array}{c}{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{n-\mathrm{1}})\left|{}_n\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_n)\left|{}_{n+\mathrm{1}}\right.\end{array}\right]（14）

进而，通过式（13）和式（14），可得：

\left\{\begin{array}{l}{\overline{\mathit{K}}}_{\mathrm{1}}^n=-\left[\left({\mathit{P}}_{\mathrm{1}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{3}n}\right){\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{1}n}^{-\mathrm{1}}\right.+\\ \left.\left({\mathit{P}}_{\mathrm{2}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{4}n}\right){\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{3}n}^{-\mathrm{1}}\right]{\overline{\mathit{K}}}_{\mathrm{2}}^n\\ {\overline{\mathit{K}}}_{\mathrm{2}}^n=\left[\left({\mathit{P}}_{\mathrm{1}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{3}n}\right){\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{2}n}^{-\mathrm{1}}\right.+\\ {\left.\left({\mathit{P}}_{\mathrm{2}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{4}n}\right){\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{4}n}^{-\mathrm{1}}\right]}^{-\mathrm{1}}\\ {\overline{\mathit{K}}}_{\mathrm{3}}^n={\mathit{P}}_{\mathrm{3}n}{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{2}n}^{-\mathrm{1}}+\\ {\mathit{P}}_{\mathrm{4}n}{\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{4}n}^{-\mathrm{1}}-\\ {\overline{\mathit{K}}}_{\mathrm{4}}^n\left[\left({\mathit{P}}_{\mathrm{1}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{3}n}\right){\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{1}n}^{-\mathrm{1}}\right.+\\ \left.\left({\mathit{P}}_{\mathrm{2}n}+{\tilde{\mathit{F}}}_n{\mathit{P}}_{\mathrm{4}n}\right){\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{3}n}^{-\mathrm{1}}\right]\\ {\overline{\mathit{K}}}_{\mathrm{4}}^n=\left[{\mathit{P}}_{\mathrm{3}n}{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}){\mathit{P}}_{\mathrm{2}n}^{-\mathrm{1}}+\right.\\ {\mathit{P}}_{\mathrm{4}n}{\mathit{\Lambda }}_n(z_n-z_{n-\mathrm{1}})\left.{\mathit{P}}_{\mathrm{4}n}^{-\mathrm{1}}\right]{\overline{\mathit{K}}}_{\mathrm{2}}^n\end{array}\right.（15）

可见，在计算{\overline{\mathit{K}}}_{\mathrm{1}}^n~{\overline{\mathit{K}}}_{\mathrm{4}}^n时，通过相减、相除运算大大避免了数据溢出问题。代入位移边界条件，可得：{\overline{\mathit{K}}}_{\mathrm{1}}^N{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{N-\mathrm{1}})\left|{}_N\right.={\widetilde{\overline{\mathit{S}}}}_{\mathrm{2}}(z_{N-\mathrm{1}})\left|{}_N\right.。当z_N\to +\mathrm{\infty }时，无法由式（15）直接计算{\overline{\mathit{K}}}_{\mathrm{1}}^N；由\widetilde{\overline{\mathit{S}}}(z_{N-\mathrm{1}})\left|{}_N\right.={\mathit{P}}_N^{}[x_{\mathrm{1}N},{x_{\mathrm{2}}}_N{,\mathrm{0,0}]}^{\mathrm{T}}，进而可得：{\overline{\mathit{K}}}_{\mathrm{1}}^N={\mathit{P}}_{\mathrm{3}N}^{}/{\mathit{P}}_{\mathrm{1}N}^{}。结合单元刚度矩阵，可构建下列整体刚度矩阵方程：

{\left[\begin{array}{ccccc}{\overline{\mathit{K}}}_{\mathrm{1}}^{\mathrm{1}}& {\overline{\mathit{K}}}_{\mathrm{2}}^{\mathrm{1}}& & & \\ -{\overline{\mathit{K}}}_{\mathrm{3}}^{\mathrm{1}}& -{\overline{\mathit{K}}}_{\mathrm{4}}^{\mathrm{1}}+{\overline{\mathit{K}}}_{\mathrm{1}}^{\mathrm{2}}& {\overline{K}}_{\mathrm{2}}^{\mathrm{2}}& & \\ & -{\overline{\mathit{K}}}_{\mathrm{3}}^n& -{\overline{\mathit{K}}}_{\mathrm{4}}^n+{\overline{\mathit{K}}}_{\mathrm{1}}^{n+\mathrm{1}}& {\overline{\mathit{K}}}_{\mathrm{2}}^{n+\mathrm{1}}& \\ & & -{\overline{\mathit{K}}}_{\mathrm{3}}^{N-\mathrm{2}}& -{\overline{\mathit{K}}}_{\mathrm{4}}^{N-\mathrm{2}}+{\overline{K}}_{\mathrm{1}}^{N-\mathrm{1}}& {\overline{\mathit{K}}}_{\mathrm{2}}^{N-\mathrm{1}}\\ & & & -{\overline{\mathit{K}}}_{\mathrm{3}}^{N-\mathrm{1}}& -{\overline{\mathit{K}}}_{\mathrm{4}}^{N-\mathrm{1}}+{\overline{\mathit{K}}}_{\mathrm{1}}^N\end{array}\right]}_{\mathrm{2}N}\times \left[\begin{array}{c}{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(\mathrm{0})\left|{}_{\mathrm{1}}\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{\mathrm{1}})\left|{}_{\mathrm{2}}\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{n-\mathrm{1}})\left|{}_n\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{N-\mathrm{2}})\left|{}_{N-\mathrm{1}}\right.\\ {\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{N-\mathrm{1}})\left|{}_N\right.\end{array}\right]=\left[\begin{array}{c}\mathrm{0}\\ -\tilde{\overline{p}}(\xi ,r)\\ {\mathrm{0}}_{(\mathrm{2}N-\mathrm{2})\times \mathrm{1}}\end{array}\right]（16）

求解矩阵方程（式（16）），即可得到各结构层表面的状态向量{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_{n-\mathrm{1}})\left|{}_n\right.。刚度矩阵法由于需要构建2N阶的矩阵方程，相比传递矩阵法需更多计算时间，但基本不存在数据溢出问题。

2.3　波传递法

令{\left[\begin{array}{cc}{W}_d(z)& W_u(z)\end{array}\right]}^{\mathrm{T}}={\mathit{P}}^{-\mathrm{1}}\widetilde{\overline{\mathit{S}}}(z)，称W_d为下行波，W_u为上行波。在式（8）的基础上，可构建如式（17）所示的波传递方程，即第n层结构层的某一深度z(z_{n-\mathrm{1}}\le z\le z_n)处波向量由z_{n-\mathrm{1}}处的下行波W_d(z_{n-\mathrm{1}})\left|{}_n\right.以及z_n处的上行波W_u(z_n)\left|{}_n\right.确定。

{\left[\begin{array}{cc}{W}_d(z)\left|{}_n\right.& W_u(z)\left|{}_n\right.\end{array}\right]}^{\mathrm{T}}=

\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[{\mathit{\Lambda }}^{-\mathrm{1}}(z-z_{n-\mathrm{1}}),{\mathit{\Lambda }}^{-\mathrm{1}}(z_n-z_{})\right]\times

{\left[\begin{array}{cc}{W}_d(z_{n-\mathrm{1}})\left|{}_n\right.& W_u(z_n)\left|{}_n\right.\end{array}\right]}^{\mathrm{T}}（17）

代入层间条件，任意相邻结构层应满足下式：

\begin{array}{l}\begin{array}{l}{\mathit{F}}_n{\mathit{P}}_n\left[\begin{array}{cc}{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}})& \\ & {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot \\ \left[\begin{array}{c}{W}_d(z_{n-\mathrm{1}})\left|{}_n\right.\\ W_u(z_n)\left|{}_n\right.\end{array}\right]-\end{array}\\ {\mathit{P}}_{n+\mathrm{1}}\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}& \\ & {\mathit{\Lambda }}_{n+\mathrm{1}}^{-\mathrm{1}}(z_{n+\mathrm{1}}-z_n)\end{array}\right]\cdot \end{array}

\left[\begin{array}{c}{W}_d(z_n)\left|{}_{n+\mathrm{1}}\right.\\ W_u(z_{n+\mathrm{1}})\left|{}_{n+\mathrm{1}}\right.\end{array}\right]={\mathrm{0}}_{\mathrm{4}\times \mathrm{1}}（18）

代入应力边界条件，可得：

\left[\begin{array}{cc}{\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}& {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]{\mathit{P}}_{\mathrm{1}}\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}& \\ & {\mathit{\Lambda }}_{\mathrm{1}}^{-\mathrm{1}}(z_{\mathrm{1}})\end{array}\right]\cdot \left[\begin{array}{c}{W}_d(z_{\mathrm{0}})\left|{}_{\mathrm{1}}\right.\\ W_u(z_{\mathrm{1}})\left|{}_{\mathrm{1}}\right.\end{array}\right]=

\left[\begin{array}{c}\mathrm{0}\\ -\tilde{\overline{p}}(\xi ,r)\end{array}\right]（19）

代入位移边界条件，可得：

\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}& {\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]{\mathit{P}}_N\left[\begin{array}{cc}{\mathit{\Lambda }}_N^{-\mathrm{1}}(h_N-h_{N-\mathrm{1}})& \\ & {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot \left[\begin{array}{c}{W}_d(z_{N-\mathrm{1}})\left|{}_N\right.\\ W_u(z_N)\left|{}_N\right.\end{array}\right]={\mathrm{0}}_{\mathrm{4}\times \mathrm{1}}（20）

若z_N\to +\mathrm{\infty }时，由式（18）可得：

{\mathit{F}}_{N-\mathrm{1}}{\mathit{P}}_{N-\mathrm{1}}\left[\begin{array}{cc}{\mathit{\Lambda }}_{N-\mathrm{1}}(h_{N-\mathrm{1}}-h_{N-\mathrm{2}})& \\ & {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot

\left[\begin{array}{c}{W}_d(z_{N-\mathrm{2}})\left|{}_{N-\mathrm{1}}\right.\\ W_u(z_{N-\mathrm{1}})\left|{}_{N-\mathrm{1}}\right.\end{array}\right]-

{\mathit{P}}_N\left[\begin{array}{c}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\\ {\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot W_d(z_{N-\mathrm{1}})\left|{}_N\right.={\mathrm{0}}_{\mathrm{4}\times \mathrm{1}}（21）

结合式（17）~（21），可构造下列矩阵方程。其中，式（22）为有限厚度情况，式（23）为半空间情况。

{\left[\begin{array}{ccc}\left[\begin{array}{cc}{\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}& {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot {\mathit{P}}_{\mathrm{1}}{\tilde{\mathit{\Lambda }}}_{\mathrm{12}}^{}(z_{\mathrm{1}})& & \\ {\mathit{F}}_{\mathrm{1}}{\mathit{P}}_{\mathrm{1}}{\tilde{\mathit{\Lambda }}}_{\mathrm{11}}^{}(z_{\mathrm{1}})& -{\mathit{P}}_{\mathrm{2}}{\tilde{\mathit{\Lambda }}}_{\mathrm{22}}^{}(z_{\mathrm{2}})& \\ & ...& ...\\ & & \left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}& {\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot {\mathit{P}}_N{\tilde{\mathit{\Lambda }}}_{N\mathrm{1}}^{}(z_N)\end{array}\right]}_{\mathrm{4}N}\times \left[\begin{array}{c}{W}_d(z_{\mathrm{0}})\left|{}_{\mathrm{1}}\right.\\ W_u(z_{\mathrm{1}})\left|{}_{\mathrm{1}}\right.\\ ...\\ W_d(z_{N-\mathrm{1}})\left|{}_N\right.\\ W_u(z_N)\left|{}_N\right.\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\mathrm{0}\\ -\tilde{\overline{p}}(\xi ,r)\end{array}\\ {\mathrm{0}}_{(\mathrm{4}N-\mathrm{2})\times \mathrm{1}}\end{array}\right]（22）

{\left[\begin{array}{ccc}\left[\begin{array}{cc}{\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}& {\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\cdot {\mathit{P}}_{\mathrm{1}}{\tilde{\mathit{\Lambda }}}_{\mathrm{12}}^{}(z_{\mathrm{1}})& & \\ {\mathit{F}}_{\mathrm{1}}{\mathit{P}}_{\mathrm{1}}{\tilde{\mathit{\Lambda }}}_{\mathrm{11}}^{}(z_{\mathrm{1}})& -{\mathit{P}}_{\mathrm{2}}{\tilde{\mathit{\Lambda }}}_{\mathrm{22}}^{}(z_{\mathrm{2}})& \\ & ...& \\ & {\mathit{F}}_{N-\mathrm{1}}{\mathit{P}}_{N-\mathrm{1}}{\tilde{\mathit{\Lambda }}}_{(N-\mathrm{1})\mathrm{1}}^{}(z_{N-\mathrm{1}})& -{\mathit{P}}_N\cdot \left[\begin{array}{c}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\\ {\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right]\end{array}\right]}_{\mathrm{4}N-\mathrm{2}}\times \left[\begin{array}{c}{W}_d(z_{\mathrm{0}})\left|{}_{\mathrm{1}}\right.\\ W_u(z_{\mathrm{1}})\left|{}_{\mathrm{1}}\right.\\ ...\\ W_d(z_{N-\mathrm{1}})\left|{}_N\right.\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\mathrm{0}\\ -\tilde{\overline{p}}(\xi ,r)\end{array}\\ {\mathrm{0}}_{(\mathrm{4}N-\mathrm{4})\times \mathrm{1}}\end{array}\right]（23）

式中：

{\tilde{\mathit{\Lambda }}}_{n\mathrm{1}}^{}(z_n-z_{n-\mathrm{1}})=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}}),{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\right]，{\tilde{\mathit{\Lambda }}}_{n\mathrm{2}}^{}(z_n-z_{n-\mathrm{1}})=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[{\mathit{I}}_{\mathrm{2}\times \mathrm{2}},{\mathit{\Lambda }}_n^{-\mathrm{1}}(z_n-z_{n-\mathrm{1}})\right]。

求解式（22）或式（23）所示的矩阵方程，可获得{\begin{array}{cc}[W_d(z_{n-\mathrm{1}})\left|{}_n\right.& W_u(z_n)\left|{}_n\right.\end{array}]}^{\mathrm{T}}，结合式（17）即可获得任意深度处的状态向量。波传递法改写了传递矩阵，使矩阵仅包含负指数，因而具有良好的数值稳定性。相比刚度矩阵法，波传递法的矩阵规模更大，因而在构建矩阵方程以及求解上需消耗更多的时间。

2.4　反射、透射矩阵法

通过式（8）可知\widetilde{\overline{\mathit{S}}}(z)最后的表达式为{\mathrm{e}}^{-az},{\mathrm{e}}^{-bz},e^{az},{\mathrm{e}}^{bz}组合形式。因此，\widetilde{\overline{\mathit{S}}}(z)的通解也可写为\widetilde{\overline{\mathit{S}}}(z)=\mathit{M}(\xi ,s){\mathrm{e}}^{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(-a,-b,a,b)z}\mathit{\Phi }的形式，该形式可通过建立仅含变量u_r^{}和u_z^{}的偏微分方程获得。其中，\mathit{\Phi }=[x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}}{]}^{\mathrm{T}}为由边界条件确定的4×1阶待定系数矩阵，\mathit{M}(\xi ,s)具体表达形式如下。

\mathit{M}(\xi ,s)=\left[\begin{array}{cc}{\mathit{D}}_d& {\mathit{D}}_u\\ {\mathit{S}}_d& {\mathit{S}}_u\end{array}\right]（24）

式中：

{\mathit{D}}_d=\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ -\frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\alpha {\phi }_{\mathrm{2}}}& -\frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\delta {\phi }_{\mathrm{2}}}\end{array}\right]；

{\mathit{D}}_u=\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\alpha {\phi }_{\mathrm{2}}}& \frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\delta {\phi }_{\mathrm{2}}}\end{array}\right]；

{\mathit{S}}_d=\left[\begin{array}{cc}-C_{\mathrm{44}}\left(\alpha -\frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\alpha {\phi }_{\mathrm{2}}}\xi \right)& -C_{\mathrm{44}}\left(\delta -\frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\delta {\phi }_{\mathrm{2}}}\xi \right)\\ \xi C_{\mathrm{13}}+C_{\mathrm{33}}\frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{{\phi }_{\mathrm{2}}}& \xi C_{\mathrm{13}}+C_{\mathrm{33}}\frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{{\phi }_{\mathrm{2}}}\end{array}\right]；

{\mathit{S}}_u=\left[\begin{array}{cc}{C}_{\mathrm{44}}\left(\alpha -\frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\alpha {\phi }_{\mathrm{2}}}\xi \right)& C_{\mathrm{44}}\left(\delta -\frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{\delta {\phi }_{\mathrm{2}}}\xi \right)\\ \xi C_{\mathrm{13}}+C_{\mathrm{33}}\frac{{\alpha }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{{\phi }_{\mathrm{2}}}& \xi C_{\mathrm{13}}+C_{\mathrm{33}}\frac{{\delta }^{\mathrm{2}}-{\phi }_{\mathrm{1}}}{{\phi }_{\mathrm{2}}}\end{array}\right]；{\phi }_{\mathrm{1}}=(C_{\mathrm{11}}{\xi }^{\mathrm{2}}+\rho s^{\mathrm{2}})C_{\mathrm{44}}^{-\mathrm{1}}；{\phi }_{\mathrm{2}}=\xi (C_{\mathrm{13}}+C_{\mathrm{44}})C_{\mathrm{44}}^{-\mathrm{1}}；

{\phi }_{\mathrm{3}}=(C_{\mathrm{44}}{\xi }^{\mathrm{2}}+\rho s^{\mathrm{2}})C_{\mathrm{33}}^{-\mathrm{1}}；{\phi }_{\mathrm{4}}=\xi (C_{\mathrm{13}}+C_{\mathrm{44}})C_{\mathrm{33}}^{-\mathrm{1}}。

类似于波传递法，令\widetilde{\overline{\mathit{S}}}(z)=\mathit{M}(\xi ,s)\mathit{w}(z)，其中\mathit{w}(z)=[{\mathit{w}}_d^{\mathrm{T}}(z),{\mathit{w}}_u^{\mathrm{T}}{(z)]}^{\mathrm{T}}，{\mathit{w}}_d^{}(z)=[x_{\mathrm{1}}{\mathrm{e}}^{-\alpha z},x_{\mathrm{2}}{\mathrm{e}}^{-\delta z}{]}^{\mathrm{T}}称为上行波，{\mathit{w}}_u^{}(z)=[x_{\mathrm{3}}{\mathrm{e}}^{\alpha z},x_{\mathrm{4}}{\mathrm{e}}^{\delta z}{]}^{\mathrm{T}}称为下行波。

假设z、z_R为同一层内的不同深度且满足，设{\mathit{T}}^{(w)}(z,z_R)为波传播矩阵，且满足\mathit{w}(z)={\mathit{T}}^{(w)}(z,z_R)\mathit{w}(z_R)。可得：

{\mathit{T}}^{(w)}(z,z_R)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[{\mathit{\Lambda }}^{-\mathrm{1}}(z-z_R),\mathit{\Lambda }(z-z_R)\right]

由此传递矩阵\mathit{\Xi }(z,z_R)也可表示为：

\mathit{\Xi }(z,z_R)=\mathit{M}(\xi ,s){\mathit{T}}^{(w)}(z,z_R){\mathit{M}}^{-\mathrm{1}}(\xi ,s)

假设z_n^{(+)}和z_n^{(-)}从上下两侧无限接近z=z_n的截面坐标，且满足。设{\mathit{T}}^{(w)}(z_n^{(+)},z_n^{(-)})为波传播矩阵，即：

\mathit{w}(z_n^{(+)})={\mathit{T}}^{(w)}(z_n^{(+)},z_n^{(-)})\mathit{w}(z_n^{(-)})。

结合层间条件：

{\mathit{T}}^{(w)}(z_i^{(+)},z_i^{(-)})={\mathit{M}}_{n+\mathrm{1}}^{-\mathrm{1}}(\xi ,s){\mathit{F}}_n{\mathit{M}}_n^{}(\xi ,s)

假设任意两个截面z_{\alpha }和z_{\beta }（）处的波向量分别为\mathit{w}(z_{\alpha })和\mathit{w}(z_{\beta })，并满足关系：

\left[\begin{array}{c}{\mathit{w}}_d(z_{\beta })\\ {\mathit{w}}_u(z_{\beta })\end{array}\right]=\left[\begin{array}{cc}{\mathit{T}}_{dd}^{(w)}(z_{\beta }^{},z_{\alpha }^{})& {\mathit{T}}_{du}^{(w)}(z_{\beta }^{},z_{\alpha }^{})\\ {\mathit{T}}_{ud}^{(w)}(z_{\beta }^{},z_{\alpha }^{})& {\mathit{T}}_{uu}^{(w)}(z_{\beta }^{},z_{\alpha }^{})\end{array}\right]\left[\begin{array}{c}{\mathit{w}}_d(z_{\alpha })\\ {\mathit{w}}_u(z_{\alpha })\end{array}\right]（25）

式中：{\mathit{T}}_{dd}^{(w)}、{\mathit{T}}_{du}^{(w)}、{\mathit{T}}_{ud}^{(w)}、{\mathit{T}}_{uu}^{(w)}为{\mathit{T}}^{(w)}(z_{\beta }^{},z_{\alpha }^{})的分块矩阵。

引入反射和透射矩阵，改写式（25）为：

\left[\begin{array}{c}{\mathit{w}}_u(z_{\alpha })\\ {\mathit{w}}_d(z_{\beta })\end{array}\right]=\left[\begin{array}{cc}{\mathit{R}}_d^{}(z_{\beta }^{},z_{\alpha }^{})& {\mathit{T}}_u^{}(z_{\beta }^{},z_{\alpha }^{})\\ {\mathit{T}}_d^{}(z_{\beta }^{},z_{\alpha }^{})& {\mathit{R}}_u^{}(z_{\beta }^{},z_{\alpha }^{})\end{array}\right]\left[\begin{array}{c}{\mathit{w}}_d(z_{\alpha })\\ {\mathit{w}}_u(z_{\beta })\end{array}\right]（26）

式中：{\mathit{R}}_d^{}(z_{\beta }^{},z_{\alpha }^{})和{\mathit{R}}_u^{}(z_{\beta }^{},z_{\alpha }^{})分别为下行波、上行波的反射矩阵；{\mathit{T}}_d^{}(z_{\beta }^{},z_{\alpha }^{})和{\mathit{T}}_u^{}(z_{\beta }^{},z_{\alpha }^{})分别为下行波、上行波的透射矩阵。

联立式（25）和式（26），传递矩阵{\mathit{T}}^{(w)}(z_{\beta }^{},z_{\alpha }^{})与反射、透射矩阵有如下对应关系：

\left\{\begin{array}{l}{\mathit{R}}_d^{}(z_{\beta }^{},z_{\alpha }^{})=-{\mathit{T}}_{uu}^{(w)-\mathrm{1}}(z_{\beta }^{},z_{\alpha }^{}){\mathit{T}}_{ud}^{(w)}(z_{\beta }^{},z_{\alpha }^{})\\ {\mathit{R}}_u^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{T}}_{du}^{(w)}(z_{\beta }^{},z_{\alpha }^{}){\mathit{T}}_{uu}^{(w)-\mathrm{1}}(z_{\beta }^{},z_{\alpha }^{})\\ {\mathit{T}}_d^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{T}}_{dd}^{(w)}(z_{\beta }^{},z_{\alpha }^{})-{\mathit{T}}_{du}^{(w)}(z_{\beta }^{},z_{\alpha }^{})\times \\ {\mathit{T}}_{uu}^{(w)-\mathrm{1}}(z_{\beta }^{},z_{\alpha }^{}){\mathit{T}}_{ud}^{(w)}(z_{\beta }^{},z_{\alpha }^{})\\ {\mathit{T}}_u^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{T}}_{uu}^{(w)-\mathrm{1}}(z_{\beta }^{},z_{\alpha }^{})\end{array}\right.（27）

观察式（27），若z_{\alpha }和z_{\beta }位于同一土层时，{\mathit{T}}_d^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{T}}_u^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{\Lambda }}^{-\mathrm{1}}(z_{\beta }^{}-z_{\alpha }^{})，{\mathit{R}}_d^{}(z_{\beta }^{},z_{\alpha }^{})={\mathit{R}}_u^{}(z_{\beta }^{},z_{\alpha }^{})={\mathrm{0}}_{\mathrm{2}\times \mathrm{2}}，可见反射、投射矩阵消除了正指数的影响。假设截面z_{\alpha }、z_{\beta }、z_{\gamma }（）处的波向量分别为\mathit{w}(z_{\alpha })、\mathit{w}(z_{\beta })、\mathit{w}(z_{\gamma })，可得到递推关系：

\mathit{w}(z_{\gamma })={\mathit{T}}^{(W)}(z_{\gamma },z_{\beta }){\mathit{T}}^{(W)}(z_{\beta },z_{\alpha })\mathit{w}(z_{\alpha })（28）

联立式（25）~式（28），有：

\left\{\begin{array}{l}{\mathit{R}}_d(z_{\gamma },z_{\alpha })={\mathit{R}}_d(z_{\beta },z_{\alpha })+{\mathit{T}}_u(z_{\beta },z_{\alpha })\times \\ {\mathit{\Gamma }}_a^{-\mathrm{1}}(z_{\gamma },z_{\beta },z_{\alpha }){\mathit{R}}_d(z_{\gamma },z_{\beta }){\mathit{T}}_d(z_{\beta },z_{\alpha })\\ {\mathit{R}}_u(z_{\gamma },z_{\alpha })={\mathit{R}}_u(z_{\gamma },z_{\beta })+{\mathit{T}}_d(z_{\gamma },z_{\beta }){\mathit{R}}_u(z_{\beta },z_{\alpha })\times \\ {\mathit{\Gamma }}_a^{-\mathrm{1}}(z_{\gamma },z_{\beta },z_{\alpha })T_u(z_{\gamma },z_{\beta })\\ {\mathit{T}}_d(z_{\gamma },z_{\alpha })={\mathit{T}}_d(z_{\gamma },z_{\beta }){\mathit{\Gamma }}_b^{-\mathrm{1}}(z_{\gamma },z_{\beta },z_{\alpha }){\mathit{T}}_d(z_{\beta },z_{\alpha })\\ {\mathit{T}}_u(z_{\gamma },z_{\alpha })={\mathit{T}}_u(z_{\beta },z_{\alpha }){\mathit{\Gamma }}_a^{-\mathrm{1}}(z_{\gamma },z_{\beta },z_{\alpha }){\mathit{T}}_u(z_{\gamma },z_{\beta })\\ {\mathit{\Gamma }}_a^{}(z_{\gamma },z_{\beta },z_{\alpha })={\mathit{I}}_{\mathrm{2}\times \mathrm{2}}-{\mathit{R}}_d(z_{\gamma },z_{\beta }){\mathit{R}}_u(z_{\beta },z_{\alpha })\\ {\mathit{\Gamma }}_b^{}(z_{\gamma },z_{\beta },z_{\alpha })={\mathit{I}}_{\mathrm{2}\times \mathrm{2}}-{\mathit{R}}_u(z_{\beta },z_{\alpha }){\mathit{R}}_d(z_{\gamma },z_{\beta })\end{array}\right.（29）

代入位移边界条件{\widetilde{\overline{\mathit{S}}}}_{\mathrm{1}}(z_N)\left|{}_N\right.={\mathrm{0}}_{\mathrm{2}\times \mathrm{1}}，可得：

{\mathit{w}}_u^{}(z_N)=-{\mathit{D}}_{uN}^{-\mathrm{1}}(\xi ,s){\mathit{D}}_d(\xi ,s){\mathit{w}}_d^{}(z_N)

结合式（26），可知土层L_s(z_N,\mathrm{0})中波向量{\mathit{w}}_u(\mathrm{0})、{\mathit{w}}_d(z_N)满足{\mathit{w}}_u(\mathrm{0})=\tilde{\mathit{K}}{\mathit{w}}_d(\mathrm{0})，其中\tilde{\mathit{K}}如式（30）所示。

\begin{array}{l}\tilde{\mathit{K}}={\mathit{R}}_d^{}(z_N,\mathrm{0})-{\mathit{T}}_u^{}(z_N,\mathrm{0}){\mathit{D}}_u^{-\mathrm{1}}(z_N){\mathit{D}}_d(z_N)\times \\ {\left[{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}+{\mathit{R}}_u^{}(z_N,\mathrm{0}){\mathit{D}}_u^{-\mathrm{1}}(z_N){\mathit{D}}_d(z_N)\right]}^{-\mathrm{1}}{\mathit{T}}_d^{}(z_N,\mathrm{0})\end{array}（30）

当z_N\to +\mathrm{\infty }时，有\tilde{\mathit{K}}={\mathit{R}}_d^{}(\mathrm{\infty },\mathrm{0})。结合荷载边界条件以及式（26）和式（30），可求得：

{\mathit{w}}_d(\mathrm{0})={\left[{\mathit{S}}_d(\mathrm{0})+{\mathit{S}}_u(\mathrm{0})\tilde{\mathit{K}}\right]}^{-\mathrm{1}}\cdot {\left[\begin{array}{cc}\mathrm{0}& -\tilde{\overline{p}}(\xi ,r)\end{array}\right]}^{\mathrm{T}}（31）

式中：\tilde{\mathit{K}}可反复应用式（27）~式（29）获得。

在求得{\mathit{w}}_d(\mathrm{0})以及{\mathit{w}}_u(\mathrm{0})后，结合式（26）即可获得任意位置的状态向量。与波传递法类似，反射、透射矩阵的引入消除了正指数。总体上，反射、透射矩阵法是一种矩阵规模小、数值稳定性高的计算方法。

3 数值积分收敛分析及解析解验证

3.1　数值积分方法及收敛分析

前述解析解方法仅能获得Hankel-Laplace域内的力学响应，为计算在时域内的力学响应，往往需通过数值积分逆变换。考虑到精度要求，本文采用20节点高斯插值积分以及Durbin法^［34］分别进行Hankel逆变换和Laplace逆变换。

f(r,t_l)=\frac{\mathrm{2}{\mathrm{e}}^{\gamma l\mathrm{\Delta }t}}{T}{R}_{\mathrm{e}}\left\{-\frac{\mathrm{1}}{\mathrm{2}}\overline{f}(r,\gamma )+\right.

{\displaystyle \sum _{k=\mathrm{0}}^{N_a-\mathrm{1}}}\left.{\displaystyle \sum _{\nu =\mathrm{0}}^L}\left[\overline{f}\left(r,\gamma +j(k+\nu N_a)\frac{\mathrm{2}\mathrm{\pi }}{T}\right)\cdot e^{j\frac{\mathrm{2}\mathrm{\pi }}{N_a}kl}\right]\right\}（32）

\overline{f}\left(r,s\right)\approx {\int }_{\mathrm{0}}^{\chi }\xi \tilde{\overline{f}}\left(\xi ,s\right)J(\xi r)\mathrm{d}\xi \approx

{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{0}}}}{\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{20}}}{A}_k\left[x_{ik}\tilde{\overline{f}}\left(x_{ik},s\right)J(x_{ik}r)\right]（33）

式中：χ为积分上限；N₀为高斯积分分段总数；A_k 为20节点高斯积分系数；x_ik 为积分点，x_ik =(i-1)∙ΔL+x_k，ΔL为高斯积分的分段长度，x_k 为20节点高斯积分点；T为总求解时间；N_a 为总的求解步；Δt为求解时间步增量；t_l 为求解时刻，t_l =l∙T/N_a；L和γ为计算参数，通常需满足L∙N_a =50~5 000，γ∙T=5~10，l为求解序列；取L=ceil(100/N_a )，γ=10/T，ceil(∙)为向上取整计算；ς为单位虚数，ς²=-1。

实际计算时，需确定积分上限χ以及分段长度ΔL^［35］。在表1所示的路面结构基础上，分两步开展收敛研究。首先，在ΔL=5基础上，χ分别取值50，100，150，200，300，500；其次，在χ=300的基础上，ΔL分别取2，3，5，10，20，30，结果如图2~图5所示。由图2可见，4种方法的计算结果较为接近，初步验证了解析解推导的正确性。由图3~图5可知，ΔL取值对计算结果影响较大，且χ≥300，ΔL≤5时计算结果基本收敛。考虑到计效率以及精度要求，本文取χ=300，ΔL=5。

3.2　解析解验证

由于前述4种解析解方法具有高度一致性，因而这里以反射、透射矩阵法为例进行进一步的文献对比验证，所采用的模型参数及对比结果如图6所示。结果表明，所推导的解析解与相关文献结果吻合较好，计算精度较高。

4 解析解方法对比

在验证了所推导解析解可靠性的基础上，采用图7所示的4种路面结构进行测试；其中，路面结构A和C为半空间结构，路面结构B和D为有限厚度结构。本文主要以FWD荷载作用为计算目标，限于篇幅有限，本文仅展示针对路表弯沉的讨论，而不涉及剪应力和应变，这并不意味所推导解析解仅能用于计算路表弯沉。输出为r=0.0和0.3 m位置处的路表弯沉，求解时间T=60 ms，增量Δt=1 ms。相关计算均采用Matlab编程实现，耗时采用Matlab探查器获得，计算所用处理器为Intel i7 8750H （2.2 GHz）。相关计算结果如表2、图8所示。

通过表2和图8的结果可以看到，4种方法计结果基本相同，且从快到慢依次为传递矩阵法、刚度矩阵法、波传递法以及反射透射矩阵法。其中，传递矩阵法的计算耗时仅为其他方法的1/5，在保证数值稳定下可作为首选算法；刚度矩阵法的时间大多消耗在计算单元刚度矩阵上；波传递法在构建矩阵方程以及求解上较刚度矩阵法需花费更多的时间；值得注意的是，尽管反射透射矩阵法的矩阵规模小，但大量矩阵运算使得其计算速度最慢。传递矩阵法在对路面结构B和D展开计算时出现了异常结果NAN，其原因是计算式（11）和式（12）时分母为零造成的。而且，多次试算表明，在有限厚度且厚度大于2 m时较大概率发生计算失败；而其他解析解方法均消除了正指数的影响，基本不存在计算失败的情况。

在适用性方面，对于传递矩阵法，由于在计算路面内部力学响应时传递矩阵需反复相乘，较易出现数据溢出问题，因而仅能计算路面结构的路表弯沉；对于刚度矩阵法，通过构建式（16）所示的刚度矩阵方程即可获得各结构层表面的力学响应，而内部应力通常需通过增设结构层实现，这将进一步增加计算时间；对于波传递法，在获得各结构层表面的状态向量后，其内部力学响应可通过波传递矩阵进一步获得；对于反射透射矩阵法，在获得了表面力学响应的基础上，结合反射、透射矩阵同样也可计算路面结构任意位置的力学响应，有关求解源代码的更多信息可从通信作者处获得。综上所述，4种解析解方法的特点及适用范围如表3所示。

5 结 论

（1）高斯数值积分上限χ取300、分段长度ΔL取5时，各解析解均可较好地完成路面结构动力响应计算。但各个算法在求解思想上有所不同，传递矩阵法和反射、透射矩阵法侧重于寻找不同深度处状态向量间的联系，因而矩阵规模较小；而刚度矩阵法和波传递法更侧重于构建路面结构的整体矩阵方程，因而矩阵规模较大。

（2）在求解速度上，从快到慢依次为传递矩阵法、刚度矩阵法、波传递法以及反射透射矩阵法。其中，传递矩阵法的耗时仅为其他方法的1/5左右；刚度矩阵法较多时间消耗在计算单元刚度矩阵上；波传递法由于矩阵规模大而慢于刚度矩阵法；值得注意的是，尽管反射、透射矩阵法的矩阵规模最小，但大量的矩阵运算使得其计算耗时最大。

（3）在稳定性上，传递矩阵法并未对正指数进行任何处理，因而存在数据溢出问题；刚度矩阵法通过建立单元刚度矩阵消除了正指数；波传递法则通过构建上行波和下行波消除了正指数；而反射透射矩阵法在不改变传递矩阵法基本思路的基础上应用反射矩阵和透射矩阵消除了正指数的影响。

（4）在实际应用上，传递矩阵法非常适合求解半空间的路表弯沉，在满足计算可靠性的条件下可作为首选算法；刚度矩阵法在求解有限厚度力学响应中有着较为出色的计算效率，但仅能计算各结构层表面的力学响应；波传递法和反射透射矩阵法均能计算路面结构内部任意位置的动力响应，但其计算较慢。

参考文献

[1] Liu H, Pan Ernian, Cai Y C. General surface loading over layered transversely isotropic pavements with imperfect interfaces[J]. Advances in Engineering Software, 2018, 115: 268-282.

[2] Lee H S. Viscowave-a new solution for viscoelastic wave propagation of layered structures subjected to an impact load[J]. International Journal of Pavement Engineering, 2014, 15(6): 542-557.

[3] 程永春, 李赫, 李立顶, 等. 基于灰色关联度的矿料对沥青混合料力学性能的影响分析[J]. 吉林大学学报: 工学版, 2021, 51(3): 925-935.

[4] Cheng Yong-chun, Li He, Li Li-ding, et al. Analysis of mechanical properties of asphalt mixture affected by aggregate based on grey relational degree[J]. Journal of Jilin University(Engineering and Technology Edition), 2021, 51(3): 925-935.

[5] 黄宣, 白桃, 张德育, 等. 沥青混合料动静态黏弹性主曲线转化与比较[J]. 公路交通科技, 2024, 41(6): 1-8, 64.

[6] Huang Xuan, Bai Tao, Zhang De-yu, et al. Viscoelastic master curve transformation and comparison of asphalt mixture under dynamic and static loads[J]. Journal of Highway and Transportation Research and Development, 2024, 41(6): 1-8, 64.

[7] 鲁巍巍, 郑健龙. 横观各向同性黏弹性沥青路面的动力响应[J]. 中南大学学报: 自然科学版, 2018, 49(4): 964-970.

[8] Lu Wei-wei, Zheng Jian-long, Dynamic response of cross-anisotropic viscoelastic asphalt pavement[J]. Journal of Central South University (Science and Technology), 2018, 49(4): 964-970.

[9] Yan K Z, Xu H B, You L Y. Analytical layer-element approach for wave propagation of transversely isotropic pavement[J]. International Journal of Pavement Engineering, 2016, 17(03): 275-282.

[10] 颜可珍, 游凌云, 葛冬冬, 等. 横观各向同性沥青路面结构力学行为分析[J]. 公路交通科技, 2016, 33(4): 1-6.

[11] Yan Ke-zhen, You Ling-yun, Ge Dong-dong, et al. Analysis of structural mechanical behavior of transverse isotropic asphalt pavement[J]. Journal of Highway Transportation Research and Development, 2016, 33(4): 1-6.

[12] 李崛, 张安顺, 张军辉, 等. 级配碎石基层结构动力响应模型测试及数值分析[J]. 吉林大学学报: 工学版, 2023, 53(6): 1782-1789.

[13] Li Jue, Zhang An-shun, Zhang Jun-hui, et al. Model testing and numerical analysis of dynamic response of graded crushed rock base structure[J]. Journal of Jilin University(Engineering and Technology Edition), 2023, 53(06): 1782-1789.

[14] 薛亮, 张维刚, 梁鸿颉. 考虑层间不同状态的沥青路面力学响应分析[J]. 沈阳建筑大学学报: 自然科学版, 2006, 22(04): 575-578.

[15] Xue Liang, Zhang Wei-gang, Liang Hong-jie. The Mechanical analysis of asphalt pavement in different interface condition between layers[J]. Journal of Shenyang University (Natural Science), 2006, 22(4): 575-578.

[16] 李彦伟, 穆柯, 石鑫, 等. 基面层间接触状态对沥青路面力学响应影响[J]. 长安大学学报: 自然科学版, 2014, 34(2): 38-44.

[17] Li Yan-wei, Mu Ke, Shi Xin, et al. Impact of base-surface contact on mechanical response of asphalt pavement[J]. Journal of Chang'an University (Natural Science Edition), 2014, 34(2): 38-44.

[18] 颜可珍, 满建宏, 石挺魏, 等. 考虑层间接触状态的横观各向同性结构动力响应解析解[J]. 湖南大学学报: 自然科学版, 2019, 46(11): 97-105.

[19] Yan Ke-zhen, Man Jian-hong, Shi Ting-wei, et al. Analytical solution for dynamic response of transversely isotropic structures considering the state of interlayer contact state[J]. Journal of Hunan University (Natural Sciences), 2019, 46(11): 97-105.

[20] 陆建飞, 周慧明, 刘洋. 横观各向同性层状饱和土动力问题的反射、透射矩阵方法[J]. 岩土力学, 2018, 39(6): 2219-2226.

[21] Lu Jian-fei, Zhou Hui-ming, Liu Yang. Reflection-transmission matrix method for dynamic response of transversely isotropic multilayered saturated soil[J]. Rock and Soil Mechanics, 2018, 39(6): 2219-2226.

[22] 马宪永, 全蔚闻, 董泽蛟. 横观各向同性黏弹性沥青路面动力响应解析解[J]. 中国公路学报, 2020, 33(10): 135-145.

[23] Ma Xian-yong, Quan Wei-wen, Dong Ze-jiao. Analytical solution for dynamic response of transversely isotropic viscoelastic asphalt pavement[J]. China Journal of Highway Transport, 2020, 33(10): 135-145.

[24] 张军辉, 范海山, 张石平, 等. 考虑层间接触状态的路面动力响应解析解及参数反演[J]. 中国公路学报, 2021, 34(5): 11-23.

[25] Zhang Jun-hui, Fan Hai-shan, Zhang Shi-ping, et al. Analytical solution for dynamic response and parameter inversion of pavement structure considering the condition of interlayer contact[J]. China Journal of Highway and Transport, 2021, 34(5): 11-23.

[26] Li M Y, Wang H. Development of ANN-GA program for backcalculation of pavement moduli under FWD testing with viscoelastic and nonlinear parameters[J]. International Journal of Pavement Engineering, 2019, 20(3): 490-498.

[27] 郭大智, 冯德成. 层状弹性体系力学[M]. 哈尔滨: 哈尔滨工业大学出版社, 2001.

[28] Cai Y C, Sangghaleh A, Pan E. Effect of anisotropic base/interlayer on the mechanistic responses of layered pavements[J]. Computers and Geotechnics, 2015, 65: 250-257.

[29] 艾智勇, 成怡冲. 三维横观各向同性成层地基的传递矩阵解[J]. 岩土力学, 2010, 31(): 25-30.

[30] Ai Zhi-yong, Cheng Yi-chong. Transfer matrix solutions of three-dimensional transversely isotropic multilayered soils[J]. Rock and Soil Mechanics, 2010, 31(Sup.2): 25-30.

[31] Wang L, Rokhlin S. Stable reformulation of transfer matrix method for wave propagation in layered anisotropic media[J]. Ultrasonics, 2001, 39(6): 413-424.

[32] 艾智勇, 董洲, 成怡冲. 多层地基轴对称弹性空间问题的解析层元解[J]. 建筑结构学报, 2012, 33(4): 150-153.

[33] Ai Zhi-yong, Dong Zhou, Cheng Yi-chong. Analytical layer element solution of axisymmetrically elastic problem for multilayered soils[J]. Journal of Building Structures, 2012, 33(4): 150-153.

[34] Chen L. Three-dimensional Green's function for an anisotropic multi-layered half-space[J]. Computational Mechanics, 2015, 56(5): 795-814.

[35] Ai Z Y, Cheng Y C. Extended precise integration method for consolidation of transversely isotropic poroelastic layered media[J]. Computers & Mathematics with Applications, 2014, 68(12): 1806-1818.

[36] 张宇, 杨林青. 基于精细积分算法的混凝土路面层状结构动力响应求解与分析[J]. 中外公路, 2020, 40(5): 29-35.

[37] Zhang Yu, Yang Lin-qing. Analysis of dynamic response for layered structure of concrete pavement based on precise lntegration algorithm[J]. Journal of China & Foreign Highway, 2021, 51(3): 925-935.

[38] 卢正, 姚海林, 胡梦玲, 等. 基于传递-反射矩阵法的层状公路结构动力响应研究[J]. 岩土力学, 2012, 33(12): 3767-3774, 3809.

[39] Lu Zheng, Yao Hai-lin, Hu Meng-ling, et al. Study of dynamic response of multilayered road structures based on transmission-reflection matrix method[J]. Rock and Soil Mechanics, 2012, 33(12): 3767-3774, 3809.

[40] Thomson W T. Transmission of elastic waves through a stratified solid medium[J]. Computational Mechanics, 1950, 21(2): 89-93.

[41] Haskell N A. The dispersion of surface waves on multilayered media[J]. Bull Seismological Soc Am, 1953, 30: 17-34.

[42] Rokhlin S I, Wang L. Stable recursive algorithm for elastic wave propagation in layered anisotropic media: stiffness matrix method[J]. Journal of the Acoustical Society of America, 2002, 112(3): 822-834.

[43] Gao Q, Lin J H, Zhong W X, et al. A precise numerical method for Rayleigh waves in a stratified half space[J]. International Journal for Numerical Methods in Engineering, 2006, 67(6): 771-786.

[44] Zhang Y H, Gao Q. Stability analysis of the mixed variable method and its application in wave reflection and transmission in multilayered anisotropic structures[J]. Archive of Applied Mechanics, 2020, 90(1): 863-867.

[45] 艾智勇, 仓乃瑞, 成怡冲. 解析层元法求解层状横观各向同性地基轴对称问题[J]. 岩土工程学报, 2012, 34(5): 863-867.

[46] Ai Zhi-yong, Cang Nai-rui, Cheng Yi-chong. Analytical layer-element method for axisymmetric problem of transversely isotropic multi-layered soils[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(5): 863-867.

[47] Matsui K, Maina J W, Inoue T. Axi-symmetric analysis of elastic multilayer system considering interface slips[J]. Journal of Japan Society of Civil Engineers, 2000(5): 122-129.

[48] 任瑞波, 钟阳, 殷建华. 路面结构在动荷载作用下路表弯沉的求解[J]. 岩土工程学报, 2000, 22(6): 738-740.

[49] Ren Rui-bo, Zhong Yang, Yin Jian-hua. The solution of road surface deflection in the dynamic case[J]. Chinese Journal of Geotechnical Engineering, 2000, 22(6): 738-740.

[50] Zhang J H, Fan H S, Zhang S P, et al. Time-domain elasto-dynamic model of a transversely isotropic, layered road structure system with rigid substratum under a FWD load[J]. Road Materials and Pavement Design, 2022, 23(12): 2857-2875.

[51] Durbin F. Numerical inversion of laplace transforms: An efficient improvement to dubner and Abate's method[J]. The Computer Journal, 1974, 17(4): 371-376.

[52] 范海山, 张军辉, 郑健龙. 路基模量沿深度非均匀分布沥青路面动力解析解[J]. 岩土工程学报, 2022, 44(6): 1016-1026.

[53] Fan Hai-shan, Zhang Jun-hui, Zheng Jian-long. Analytical solution for dynamic response of asphalt pavement with subgrade modulus varying with depth[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(6): 1016-1026.

基金资助

国家自然科学基金项目(52025085)

国家自然科学基金项目(52338009)