一种用于PFC的砂土剪切带宽度自动识别方法

陈曦; 刘宗祺; 崔柳生; 唐建彬

doi:10.11779/CJGE2022S2039

摘要: 从细观尺度研究颗粒材料的应变局部化有助于深刻理解其物理力学机制，而剪切带的出现和演化是应变局部化的典型特征。目前在离散元法（DEM）模拟中对砂土应变局部化特征（如剪切带宽度、倾角、分布类型等性质）的观测方法还主要依赖手动、目测和估计等方式。针对这一问题，采用PFC^2D开展了大量密砂双轴压缩模拟试验，提出了一种基于“α-shape”和“Delaunay三角化”算法的剪切带及其宽度的自动识别方法。进而，基于多组砂土试样开展双轴压缩试验，验证了所提出的剪切带及其宽度自动识别方法的有效性。

关键词:

Abstract: Studying the strain localization of granular materials from the meso-scale is helpful for deeply understanding the physical and mechanical natures of strain localization, and the occurrence and evolution of shear bands are the typical characteristics of strain localization. In the numerical simulations of discrete element method (DEM), observing and measuring the evolution of strain localization of sand (such as the width, inclination angle and distribution type of shear bands, etc.) greatly rely on the manual, visual and estimation ways. Towards this problem, many numerical experiments on the dense sands in biaxial compression are carried out by using PFC^2D. Based on the numerical experiments, an identification method based on "α-shape" and "Delaunay triangulation" algorithms is developed to automatically identify the shear bands and their widths. Furthermore, based on some groups of biaxial compression tests on the dense sands, the effectiveness of the proposed automatic identification method is validated.

Keywords:

0. 引言

目前，落锤式弯沉仪（FWD）是应用最为广泛的路面结构快速无损检测手段之一，而层状弹性静力体系仍是FWD测试中最为常用的力学模型。在该模型中，路面结构被简化为各向同性且完全连续的弹性层状体系^[1]。现有研究表明，由于现场碾压成型的施工方式以及各结构层材料的差异，路面结构存在明显的横观各向同性以及非连续层间接触状态^[2-3]。除此之外，由于路基土具有较为明显的应力、湿度依赖特性，路基模量还存在着显著的空间非均匀分布特征。由此可见，原有层状弹性体系理论难以真实反映沥青路面的力学行为，有必要在FWD检测中综合考虑沥青路面的黏弹性、横观各向同性、层间连续条件，以及路基模量的非均匀分布特征。

路面材料方面，Wang等^[4]结合沥青混合料三轴试验，得出沥青层水平模量仅为竖向模量的20%~50%。类似的，Tutumluer等^[5]通过试验得出粒料类材料的水平与竖向的模量比为0.03~0.21。进而，不少学者针对这种特征开展了理论研究。Al-Qadi等^[6]研究了碎石基层材料横观各向同性特性对路面结构动力响应的影响，结果表明考虑粒料类基层的横观各向同性特性可以更准确地计算路面力学响应。董泽蛟等^[7]对比了移动荷载作用下各向同性以及横观各向同性路面结构应变场，结果显示基于各向同性的路面设计会低估永久变形和应变。Yan等^[8]结合解析层元法推导了考虑横观各向同性的路面结构动力响应解析解，表明横观各向同性对路表弯沉盆具有明显的影响。Ai等^[9-10]推导了简谐荷载作用下横观各向同性多层介质动力响应解析解，计算结果表明在移动或原位荷载作用下横观各向同性均能对竖向位移产生显著影响。另外，路面结构层间接触条件也是不可忽略的因素。Mousa等^[11]在对路面现场芯样进行分析后发现，现有服役阶段路面结构层间大多处于半接触状态。冉武平等^[12]指出，层间黏接力不足可引起沥青面层的寿命衰减40%~80%。颜可珍等^[13]通过对比不同层间接触状态下的路面结构应力响应，发现层间接触状况对路表弯沉、沥青层层底拉应变以及基层层底竖向应变均有较大影响。Jiang等^[14]研究了面层–基层、基层–底基层层间接触状态对路面结构寿命的影响，指出层间条件对柔性基层沥青路面寿命的影响要大于荷载分布。

路基材料方面，许多研究表明路基材料具有明显的非线性特征^[15-16]。Wang等^[17]和Li等^[18]在NCHRP 1-28A三参数模型基础上，开发了考虑基层和路基材料非线性的路面结构有限元计算方法，并结合ANN-GA算法开展了路面非线性参数反演研究。Zhang等^[19-20]和Peng等^[21]通过大量华南黏土的动三轴试验以及土水特征曲线结果，建立了考虑湿度和应力依赖的路基土动回弹模量预估模型，开发了路基非线性有限元计算方法，并结合南昌实际气候开展数值计算，指出路基结构动回弹模量具有显著的时空非均匀分布特征。Muho等^[22]在平面应变假设基础上，推导了考虑路基剪切模量沿深度方向非均匀分布的动力响应解析解，并就不同的路基剪切模量分布开展计算，表明剪切模量的非均匀分布对路表弯沉具有较大影响。

当前，路面结构力学响应计算主要分为有限元法和解析解法。有限元法对复杂边界条件和本构关系具有良好的适应性，被广泛应用于材料非线性、多物理场耦合计算^[23-24]。相比有限元法，解析解是一种计算精度高、速度快的求解方法。随着路面力学解析解的长足发展，路面结构中存在的黏弹性^[25]、横观各向同性^[26-27]甚至复杂层间条件^[28]等因素被相继考虑。但当前尚无全面考虑材料黏弹性、横观各向同性、层间接触条件以及路基模量沿深度非均匀分布等因素的沥青路面结构动力响应解析计算方法。为此，本文主要就FWD荷载作用下的路面结构动力响应解析解开展研究。首先，基于状态空间法以及刚度矩阵法获得了考虑以上特征的沥青路面动力响应解析解；然后，应用ABAQUS以及文献对比验证了解析解的准确性；最后，分析了路面及路基结构参数对路表弯沉的影响。

1. 多层动力体系构建

如图 1（a）所示，将FWD作用下的路面结构简化为一轴对称N层层状半空间体系，其每一层层底坐标分别为z1，z2，…。其中第1层—第N-1层为考虑材料横观各向同性以及层间条件的路面结构，第N层为考虑模量非均匀分布的半空间路基。如图 1（b）所示，参考Vrettos^[29-30]对关于各类颗粒土质模量非均匀分布的研究，采用下式表征路基竖向模量沿竖向的分布规律：

${E_{{\text{v}}N}}(z) = {E_{{\text{v}}0}} + ({E_{{\text{v}}\infty }} - {E_{{\text{v}}0}})(1 - {{\text{e}}^{ - \alpha z}}) ,$

(1)

图 1 沥青路面力学模型

Figure 1. Mechanical model for asphalt pavement

下载: 全尺寸图片幻灯片

式中， ${E_{{\text{v}}N}}(z)$ 为路基竖向模量沿深度方向的表达式， ${E_{{\text{v}}0}}$ 为路基顶面模量， ${E_{{\text{v}}\infty }}$ 为无限远处模量， ${E_{{\text{v}}\infty }}$ > ${E_{{\text{v}}0}}$ ， $\alpha$ 为描述模量沿深度渐变速度的无量纲系数，e为自然常数，e≈2.718，z为距离路基顶面的垂直距离。

1.1 基本方程

柱坐标系下动力平衡方程为

$\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 }^{2}{u}_{r}}{\partial {t}^{2}}, \\ \frac{\partial {\sigma }_{z}}{\partial z}+\frac{\partial {\tau }_{zr}}{\partial r}+\frac{{\tau }_{zr}}{r}=\rho \frac{{\partial }^{2}{u}_{z}}{\partial {t}^{2}}。 \end{array} \right\}$

(2)

式中 ${\sigma _r}$ ， ${\sigma _\theta }$ ， ${\sigma _z}$ ， ${\tau _{zr}}$ 为径向、环向、竖向和切向应力分量；r，z为距离坐标原点的水平和垂直距离；ρ为密度； ${u_r}$ ， ${u_z}$ 为水平以及垂直方向位移；t为时间。考虑到路基竖向模量E_v为关于z的函数，对柱坐标系下物理方程进行如下改写：

$\left[ {\begin{array}{*{20}{c}} {{\sigma _r}} \\ {{\sigma _\theta }} \\ {{\sigma _z}} \\ {{\tau _{zr}}} \end{array}} \right] = {E_{\text{v}}}\left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}}&0 \\ {{C_{12}}}&{{C_{11}}}&{{C_{13}}}&0 \\ {{C_{13}}}&{{C_{13}}}&{{C_{33}}}&0 \\ 0&0&0&{{C_{44}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\varepsilon _r}} \\ {{\varepsilon _\theta }} \\ {{\varepsilon _z}} \\ {{\gamma _{zr}}} \end{array}} \right] ,$

(3)

式中， ${C_{11}}, {\text{ }}{C_{12}}, {\text{ }}{C_{13}}, {\text{ }}{C_{33}}, {\text{ }}{C_{44}}$ 为表征材料横观各向同性的材料参数， ${C_{11}} = k{n_{\text{e}}}(1 - {n_{\text{e}}}\mu _{\text{v}}^2)$ ， ${C_{12}} = k{n_{\text{e}}}({\mu _{\text{h}}} + {n_{\text{e}}}\mu _{\text{v}}^{\text{2}})$ ， ${C_{13}} = k{n_{\text{e}}}{\mu _{\text{v}}}(1 + {\mu _{\text{v}}})$ ， ${C_{33}} = k(1 - \mu _{\text{h}}^2)$ ， ${C_{44}} = {G_{\text{v}}} = 0.5{E_{\text{v}}} \cdot$ ${(1 + {\mu _{\text{v}}})^{ - 1}}$ ， $k = {[(1 + {\mu _{\text{h}}})(1 - {\mu _{\text{h}}} - 2{n_{\text{e}}}\mu _{\text{v}}^{\text{2}})]^{ - {\text{1}}}}$ ，n_e为模量比， ${n_{\text{e}}} = {E_{\text{h}}}E_{\text{v}}^{ - 1}$ ； ${\varepsilon _r}$ ， ${\varepsilon _\theta }$ ， ${\varepsilon _z}$ ， ${\gamma _{zr}}$ 为径向、环向、竖向和切向应变分量； ${E_{\text{h}}}$ ， $\mu _{\text{v}}^{}$ 分别为垂直方向的模量以及泊松比； ${E_{\text{v}}}$ ， $\mu _{\text{h}}^{}$ 分别为水平方向上的模量以及泊松比。柱坐标下的几何方程为

$\left[ {\begin{array}{*{20}{c}} {{\varepsilon _r}} \\ {{\varepsilon _\theta }} \\ {{\varepsilon _z}} \\ {{\gamma _{zr}}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\frac{\partial }{{\partial r}}}&{\frac{{\text{1}}}{r}}&{\text{0}}&{\frac{\partial }{{\partial z}}} \\ {\text{0}}&{\text{0}}&{\frac{\partial }{{\partial z}}}&{\frac{\partial }{{\partial r}}} \end{array}} \right]^{\text{T}}}\left[ {\begin{array}{*{20}{c}} {{u_r}} \\ {{u_z}} \end{array}} \right] 。$

(4)

1.2 边界条件

通常，FWD荷载可简化为半正弦脉冲荷载。进而，图 1（a）所示的荷载以及位移边界条件可写为

$\left. \begin{array}{l} {\sigma }_{z}(0)=-{p}_{\mathrm{max}}\mathrm{sin}\left(\frac{{\rm{ \mathsf{ π} }} }{{T}_{0}}t\right)H(t)H({T}_{0}-t)H({r}_{0}-r)\text{ }, \\ {\tau }_{zr}(0)={\sigma }_{r}(0)=0\text{ }, \\ {\sigma }_{z}(\infty )={\tau }_{zr}(\infty )={\sigma }_{r}(\infty )=0\text{ }, \\ {u}_{r}(\infty )={u}_{z}(\infty )=0\text{ }, \end{array} \right\}$

(5)

式中，p_max为荷载峰值，T₀为冲击时间，r₀为作用半径，通常为0.15 m，H(∙)为阶跃函数。此外，由于实际不同路面结构层材料的差异性，不可避免的会出现层间部分甚至完全滑动的情况，这里采用Matsui等^[31]所提出的剪切弹簧来模拟路面层间接触状态：

$(1 - \alpha _x^n)\left[ {u_{r(n + 1)}^{}({z_n}) - u_{rn}^{}({z_n})} \right] = \alpha _x^n\beta _x^n\tau _{zrn}^{}({z_n}) 。$

(6)

式中下标n为路面结构层层位(1≤n≤N-1)； $\alpha _x^n$ 为第n层与第(n+1)层的层间滑移系数， $0 \leqslant \alpha _x^n < 1$ ；其中，当 $\alpha _x^n = 0$ 时为完全连续状态，当 $\alpha _x^n = 0.99{\text{9}}$ 时表示完全滑动状态； $\beta _x^n$ 为与第n层与第(n+1)层模量与泊松比有关的系数， $\beta _x^n = b\left[ {(1 + \mu _{{\text{v}}n}^{})E_{{\text{v}}n}^{ - 1} + (1 + \mu _{{\text{v}}(n + 1)}^{})E_{{\text{v}}(n + 1)}^{ - 1}} \right]$ 。除切向位移 ${u_r}$ 外，其余应力响应均为完全连续。选取 ${u_r}$ ， ${u_z}$ ， ${\tau _{zr}}$ ， ${\sigma _z}$ 为状态变量，并设S(z)为深度z处的状态向量，即S(z)=[u_r(z)，u_z(z)，τ_zr(z)， ${\sigma _z}$ (z)]。可获得如下所示的相邻层间状态向量关系：

${F_n} \cdot S\left| {_n} \right.({z_n}) = S\left| {_{n + 1}} \right.({z_n}) ,$

(7)

式中，F_n为第n层与第(n+1)层的层间接触矩阵， ${F_n} = \left[ {\begin{array}{*{20}{c}} {{I_{2 \times 2}}}&{{{\tilde F}_n}} \\ {{{\text{0}}_{2 \times 2}}}&{{I_{2 \times 2}}} \end{array}} \right]$ ， ${\tilde F_n} = \left[ {\begin{array}{*{20}{c}} {\alpha _x^n\beta _x^n{{(1 - \alpha _x^n)}^{ - 1}}}&0 \\ 0&0 \end{array}} \right]$ ，I_2×2为2阶单位矩阵，0_2×2为2阶零矩阵， $S\left| {_n} \right.(z)$ 表示第n层深度z处的状态向量。

1.3 偏微分方程组构建

对于路面结构，联立式（2）~（4），消去式（2）中变量 ${\sigma _r}$ ， ${\sigma _\theta }$ ， ${\tau _{zr}}$ ， ${\sigma _z}$ ，可构建仅含 ${u_r}$ 和 ${u_z}$ 的偏微分方程组；类似的，对于路基结构，同样可构建仅有 ${u_r}$ 和 ${u_z}$ 的偏微分方程组。

$\left. \begin{array}{l} \text{ }[{C}_{11}\left(\frac{{\partial }^{2}{u}_{r}}{\partial {r}^{2}}+\frac{\partial {u}_{r}}{\partial r}\frac{1}{r}-\frac{{u}_{r}}{{r}^{2}}\right)+\text{ }\\ \text{ }({C}_{13}+{C}_{44})\frac{{\partial }^{2}{u}_{z}}{\partial z\partial r}+{C}_{44}\frac{{\partial }^{2}{u}_{r}}{\partial {z}^{2}}]{E}_{\text{v}}=\rho \frac{{\partial }^{2}{u}_{r}}{\partial {t}^{2}}, \\ [({C}_{13}+{C}_{44})\left(\frac{{\partial }^{2}{u}_{r}}{\partial r\partial z}+\frac{1}{r}\frac{\partial {u}_{r}}{\partial z}\right)+\\ \text{ }{C}_{33}\frac{{\partial }^{2}{u}_{z}}{\partial {z}^{2}}+{C}_{44}\left(\frac{{\partial }^{2}{u}_{z}}{\partial {r}^{2}}+\frac{\partial {u}_{z}}{\partial r}\frac{1}{r}\right)]{E}_{\text{v}}=\rho \frac{{\partial }^{2}{u}_{z}}{\partial {t}^{2}}。 \end{array} \right\}$

(8)

$\left. \begin{array}{l} \text{ }[{C}_{11}\left(\frac{{\partial }^{2}{u}_{r}}{\partial {r}^{2}}+\frac{\partial {u}_{r}}{\partial r}\frac{1}{r}-\frac{{u}_{r}}{{r}^{2}}\right)+\\ ({C}_{13}+{C}_{44})\frac{{\partial }^{2}{u}_{z}}{\partial z\partial r}\text{+}{C}_{44}\frac{{\partial }^{2}{u}_{r}}{\partial {z}^{2}}]{E}_{\text{v}}(z)+\\ {C}_{44}\left(\frac{\partial {u}_{r}}{\partial z}+\frac{\partial {u}_{z}}{\partial r}\right)\frac{\partial {E}_{v}(z)}{\partial z}=\rho \frac{{\partial }^{2}{u}_{r}}{\partial {t}^{2}}, \\ [({C}_{13}+{C}_{44})\left(\frac{{\partial }^{2}{u}_{r}}{\partial r\partial z}+\frac{1}{r}\frac{\partial {u}_{r}}{\partial z}\right)+\\ {C}_{33}\frac{{\partial }^{2}{u}_{z}}{\partial {z}^{2}}+{C}_{44}\left(\frac{{\partial }^{2}{u}_{z}}{\partial {r}^{2}}+\frac{\partial {u}_{z}}{\partial r}\frac{1}{r}\right)]{E}_{\text{v}}(z)\text{+}\\ \left({C}_{13}\frac{\partial {u}_{r}}{\partial r}+{C}_{13}\frac{{u}_{r}}{r}+{C}_{33}\frac{\partial {u}_{z}}{\partial z}\right)\frac{\partial {E}_{v}(z)}{\partial z}=\rho \frac{{\partial }^{2}{u}_{z}}{\partial {t}^{2}}。 \end{array} \right\}$

(9)

2. 解析解推导

2.1 路面部分响应通解推导

对式（8）中t进行Laplace变换， ${\bar f}(s) = \int_0^{ + \infty } {f(t) \cdot }$ ${{\text{e}}^{ - st}}{\text{d}}t$ ；对r执行Hankel变换， ${\tilde f_k}(\xi ) = \int_0^{ + \infty } {rf(r){{\text{J}}_k} \cdot }$ $(\xi r){\text{d}}r$ ，J_k(∙)为k阶Bessel函数。可得关于z的常系数常微分方程组：

$\left. \begin{array}{l} \frac{{\text{d}}^{2}{\tilde{\overline{u}}}_{r\text{1}}}{\text{d}{z}^{2}}-{\varphi }_{2}\frac{\text{d}{\tilde{\overline{u}}}_{z0}}{\text{d}z}-{\varphi }_{1}{\tilde{\overline{u}}}_{r\text{1}}\text{=0 }, \\ \frac{{\text{d}}^{2}{\tilde{\overline{u}}}_{z0}}{\text{d}{z}^{2}}\text{+}{\varphi }_{4}\frac{\text{d}{\tilde{\overline{u}}}_{r\text{1}}}{\text{d}z}-{\varphi }_{3}{\tilde{\overline{u}}}_{z0}\text{=0}。 \end{array} \right\}$

(10)

式中 ${\varphi _1} = ({C_{11}}{\xi ^2}{\text{ + }}\rho {s^2}E_{\text{v}}^{ - 1})C_{44}^{ - 1}$ ； ${\varphi _2} = \xi ({C_{13}} + {C_{44}})C_{44}^{ - 1}$ ； ${\varphi _3} = {\text{(}}{C_{44}}{\xi ^2}{\text{ + }}\rho {s^2}E_{\text{v}}^{ - 1}{\text{)}}C_{33}^{ - 1}$ ； ${\varphi _4} = \xi ({C_{13}} + {C_{44}})C_{33}^{ - 1}$ ； ${\tilde {{\bar u}_{z0}}}$ 表示对变量u_z执行Laplce变换以及0阶Hankel变换的结果； ${\tilde {{\bar u}_{r{\text{1}}}}}$ 表示对变量u_r执行Laplce变换以及1阶Hankel变换的结果；s和ξ分别为Laplace变换以及Hankel变换积分变量。求解式（10）可得

$\left. \begin{array}{l} {\tilde{\overline{u}}}_{r\text{1}}={A}_{1}{\text{e}}^{-(h-a)z}+{A}_{2}{\text{e}}^{-(h-b)z}+{A}_{3}{e}^{-az}+{A}_{4}{\text{e}}^{-bz}, \\ {\tilde{\overline{u}}}_{z0}=\frac{{a}^{2}-{\varphi }_{1}}{a{\varphi }_{2}}{A}_{1}{\text{e}}^{-(h-z)a}+\frac{{b}^{2}-{\varphi }_{1}}{b{\varphi }_{2}}{A}_{2}{\text{e}}^{-(h-z)b}-\\ \frac{{a}^{2}-{\varphi }_{1}}{a{\varphi }_{2}}{A}_{3}{\text{e}}^{-az}-\frac{{b}^{2}-{\varphi }_{1}}{b{\varphi }_{2}}{A}_{4}{\text{e}}^{-bz}, \\ a, b=[0.5({\varphi }_{1}+{\varphi }_{3}-{\varphi }_{2}{\varphi }_{4})\pm \text{ }\\ {\sqrt{{({\varphi }_{2}{\varphi }_{4})}^{2}-2{\varphi }_{2}{\varphi }_{4}({\varphi }_{1}+{\varphi }_{3})+{({\varphi }_{1}-{\varphi }_{3})}^{2}}]}^{0.5}, \end{array} \right\}$

(11)

式中，参数A₁，A₂，A₃，A₄为与边界条件相关的待定常数，h为结构层厚度。结合式（3），（4），可获得 ${\tilde {{\bar \tau} _{zr1}}}$ 以及 ${\tilde {{\bar \sigma} _{z0}}}$ 的通解表达式。因而，可获得路面结构状态向量 $\tilde {\bar S}(z) = M(\xi , s) \cdot {\text{diag}}\left[ {{{\text{e}}^{ - (h - z)a}}, {{\text{e}}^{ - (h - z)b}}, {{\text{e}}^{ - az}}, {{\text{e}}^{ - bz}}} \right] \cdot w$ 。其中，w=[A₁, A₂ A₃ A₄]^T，系数矩阵M(ξ, s)为

$M(\xi , s) = \left[ {\begin{array}{*{20}{c}} {{D_{\text{d}}}}&{{D_{\text{u}}}} \\ {{S_{\text{d}}}}&{{S_{\text{u}}}} \end{array}} \right],$

(12)

式中， ${D_{\text{d}}} = \left[ {\begin{array}{*{20}{c}} 1&1 \\ {({a^2} - {\varphi _1}){{(a{\varphi _2})}^{ - 1}}}&{({b^2} - {\varphi _1}){{{\text{(}}b{\varphi _2}{\text{)}}}^{ - 1}}} \end{array}} \right]$ ， ${D_{\text{u}}} = \left[ {\begin{array}{*{20}{c}} {\text{1}}&{\text{1}} \\ { - {\text{(}}{a^2} - {\varphi _1}{\text{)(}}a{\varphi _2}{{\text{)}}^{ - 1}}}&{ - ({b^2} - {\varphi _1}){{(b{\varphi _2})}^{ - 1}}} \end{array}} \right],$ ${S_{\text{d}}} = \left[\begin{array}{cc}{C}_{44}{E}_{\text{v}}\left(a-\frac{{a}^{2}-{\varphi }_{1}}{a{\varphi }_{2}}\xi \right)& {C}_{44}{E}_{\text{v}}\left(b-\frac{{b}^{2}-{\varphi }_{1}}{b{\varphi }_{2}}\xi \right)\\ {E}_{\text{v}}\left(\xi {C}_{13}+{C}_{33}\frac{{a}^{2}-{\varphi }_{1}}{{\varphi }_{2}}\right)& {E}_{\text{v}}\left(\xi {C}_{13}+{C}_{33}\frac{{b}^{2}-{\varphi }_{1}}{{\varphi }_{2}}\right)\end{array}\right],$ ${S}_{\text{u}}= \left[ {\begin{array}{*{20}{c}} { - {C_{44}}{E_{\text{v}}}\left( {a - \frac{{{a^2} - {\varphi _1}}}{{a{\varphi _2}}}\xi } \right)}&{ - {C_{44}}{E_{\text{v}}}\left( {b - \frac{{{b^2} - {\varphi _1}}}{{b{\varphi _2}}}\xi } \right)} \\ {{E_{\text{v}}}\left( {\xi {C_{13}} + {C_{33}}\frac{{{a^2} - {\varphi _1}}}{{{\varphi _2}}}} \right)}&{{E_{\text{v}}}\left( {\xi {C_{13}} + {C_{33}}\frac{{{b^2} - {\varphi _1}}}{{{\varphi _2}}}} \right)} \end{array}} \right]。$

2.2 路基部分响应通解推导

对式（9）中t执行Laplace变换，变量r分别执行Hankel变换，可得变系数常微分方程组：

$\left. \begin{array}{l} \left[-{C}_{11}{\xi }^{2}{\tilde{\overline{u}}}_{r1}-\xi \left({C}_{13}+{C}_{44}\right)\frac{\text{d}{\tilde{\overline{u}}}_{z0}}{\text{d}z}+{C}_{44}\frac{{\text{d}}^{2}{\tilde{\overline{u}}}_{r1}}{\text{d}{z}^{2}}\right]{E}_{\text{v}N}(z)+\text{ }\\ {C}_{44}\left(\frac{\text{d}{\tilde{\overline{u}}}_{r0}}{\text{d}z}-\xi {\tilde{\overline{u}}}_{z0}\right)\frac{\text{d}{E}_{vN}(z)}{\text{d}z}=\rho {s}^{2}{\tilde{\overline{u}}}_{r1}, \\ \left[\left({C}_{13}+{C}_{44}\right)\xi \frac{\text{d}{\tilde{\overline{u}}}_{r1}}{\text{d}z}+{C}_{33}\frac{{\text{d}}^{2}{\tilde{\overline{u}}}_{z0}}{\text{d}{z}^{2}}-{C}_{44}{\xi }^{2}{\tilde{\overline{u}}}_{z0}\right]{E}_{\text{v}N}(z)+\\ \left(\xi {C}_{13}{\tilde{\overline{u}}}_{r1}+{C}_{33}\frac{\text{d}{\tilde{\overline{u}}}_{z0}}{\text{d}z}\right)\frac{\text{d}{E}_{\text{v}N}(z)}{\text{d}z}=\rho {s}^{2}{\tilde{\overline{u}}}_{z0}\text{ }。 \end{array} \right\}$

(13)

与路面结构所获得的常微分方程组不同，由于路基竖向模量E_vN(z)为与z相关函数，因而不能直接求解。这里结合Vrettos^[32-33]所提出的求解方法，引入 $\psi$ =Ξe^-αz。其中，Ξ为无量纲系数，

$\Xi = 1 - {E_{{\text{v}}0}}E_{{\text{v}}\infty }^{ - 1} 。$

(14)

进而，可得如下式所示的关于 $\psi$ 的变系数常微分方程组：

$\left. \begin{array}{l} -\left[{C}_{11}{\xi }^{2}(1-\psi )+\rho {s}^{2}{E}_{\text{v}\infty }^{-1}\right]{\tilde{\overline{u}}}_{r1}+\\ {C}_{44}{\alpha }^{2}\psi (1-2\psi ){\tilde{\overline{{u}}'}}_{r1}+{\alpha }^{2}{C}_{44}(1-\psi ){\psi }^{2}{\tilde{\overline{{u}}''}}_{r1}-\\ \alpha \psi {C}_{44}\xi {\tilde{\overline{u}}}_{z0}+\alpha \psi \xi (1-\psi )({C}_{13}+{C}_{44}){\tilde{\overline{{u}}'}}_{z0}=0\text{ }, \\ \alpha \psi \xi {C}_{13}{\tilde{\overline{u}}}_{r1}-\alpha \psi (1-\psi )({C}_{13}+{C}_{44})\xi {\tilde{\overline{{u}}'}}_{r1}-\\ \left[{C}_{44}(1-\psi ){\xi }^{2}+\rho {s}^{2}{E}_{\text{v}\infty }^{-1}\right]{\tilde{\overline{u}}}_{z0}+\\ {\alpha }^{2}{C}_{33}\psi (1-2\psi ){\tilde{\overline{{u}}'}}_{z0}+{\alpha }^{2}{C}_{33}(1-\psi ){\psi }^{2}{\tilde{\overline{{u}}''}}_{z0}=0\text{ }。 \end{array} \right\}$

(15)

式中f'，f''分别表示对 $\psi$ 的1阶和2阶导数； $\delta$ = E_v0/E_v∞。

为便于表示，后续表达式中均省略变量 $\psi$ 。采用Frobenius法，假设 ${\tilde {{\bar u}_{r1}}}$ 以及 ${\tilde {{\bar u}_{z0}}}$ 存在如下式所示的级数形式通解：

$\left. \begin{array}{l} {\tilde{\overline{u}}}_{z0}={\displaystyle {\sum }_{n=0}^{\infty }{a}_{n}{\psi }^{n+m}\text{ }}, \\ {\tilde{\overline{u}}}_{r1}={\displaystyle {\sum }_{n=0}^{\infty }{b}_{n}{\psi }^{n+m}\text{ }, } \end{array} \right\}$

(16)

式中，a_n，b_n为级数系数，m为待定常数。将位移通解表达式（式（16））代入变系数常微分方程组（式（15））中，可得下列级数等式关系：

$\left. \begin{array}{l} {\displaystyle {\sum }_{n=0}^{\infty }\left[{C}_{44}{\alpha }^{2}{(n+m)}^{2}-{C}_{11}{\xi }^{2}-\theta \right]{a}_{n}{\psi }^{n+m}+}\\ \alpha \xi ({C}_{13}+{C}_{44}){\displaystyle {\sum }_{n=0}^{\infty }(n+m){b}_{n}{\psi }^{n+m}\text{ }}\\ ={\displaystyle {\sum }_{n=1}^{\infty }\left[{\alpha }^{2}{C}_{44}(n+m-1)(n+m)-{C}_{11}{\xi }^{2}\right]{a}_{n-1}{\psi }^{n+m}+}\\ \alpha \xi {\displaystyle {\sum }_{n=1}^{\infty }\left[{C}_{13}(n+m-1)+{C}_{44}(n+m)\right]{b}_{n-1}{\psi }^{n+m}, }\\ -\alpha \xi ({C}_{13}+{C}_{44}){\displaystyle {\sum }_{n=0}^{\infty }(n+m){a}_{n}{\psi }^{n+m}\text{+}}\\ {\displaystyle {\sum }_{n=0}^{\infty }\left[{\alpha }^{2}{C}_{33}{(n+m)}^{2}-{C}_{44}{\xi }^{2}-\theta \right]{b}_{n}{\psi }^{n+m}}\\ =-\alpha \xi {\displaystyle {\sum }_{n=1}^{\infty }\left[{C}_{13}(n+m)+{C}_{44}(n+m-1)\right]{a}_{n-1}{\psi }^{n+m}}+\\ {\displaystyle {\sum }_{n=1}^{\infty }\left[{\alpha }^{2}{C}_{33}(n+m-1)(n+m)-{C}_{44}{\xi }^{2}\right]{b}_{n-1}{\psi }^{n+m}\text{ }, } \end{array} \right\}$

(17)

式中，θ=ρs²/E_v∞。若式（17）成立，则等式两侧ψ相同次幂的系数应对应相等。令n=0，若a₀与b₀存在唯一关系，参数m应满足下列4次方程。

${\alpha ^2}{C_{33}}{C_{44}}{m^4} + \left[ {{\xi ^2}(C_{13}^2 + 2{C_{13}}{C_{44}} - {C_{11}}{C_{33}})} \right. - \\ \theta ({C}_{33}+{C}_{44})]{m}^{2}+({C}_{44}{\xi }^{2}+\theta )({C}_{11}{\xi }^{2}+\theta )=0\text{ }。$

(18)

复数域内，求解式（18）所示的4次方程，m必定存在4个解，假设 $a_0^{(1)} = a_0^{(2)} = a_0^{(3)} = a_0^{(4)} = 1$ ，进而可得

$b_0^{(i)} = \frac{{\alpha \xi ({C_{13}} + {C_{44}}){m_i}}}{{{\alpha ^2}{C_{33}}m_i^2 - {C_{44}}{\xi ^2} - \theta }} ,$

(19)

式中，i=1，2，3，4。

通过式（17），可得a_n，b_n递推关系：

$\left. \begin{array}{l} {a}_{n}^{(i)}=\frac{{F}_{22}^{(i)}(n){\overline{F}}_{1}^{(i)}(n)-{F}_{12}^{(i)}(n){\overline{F}}_{2}^{(i)}(n)}{{F}_{11}^{(i)}(n){F}_{22}^{(i)}(n)-{F}_{12}^{(i)}(n){F}_{21}^{(i)}(n)}, \\ {b}_{n}^{(i)}=\frac{{F}_{11}^{(i)}(n){\overline{F}}_{2}^{(i)}(n)-{F}_{21}^{(i)}(n){\overline{F}}_{1}^{(i)}(n)}{{F}_{11}^{(i)}(n){F}_{22}^{(i)}(n)-{F}_{12}^{(i)}(n){F}_{21}^{(i)}(n)}, \end{array} \right\}$

(20)

式中，n=1, 2, 3, …，

$\left. \begin{array}{l} {F}_{11}^{(i)}(n)={C}_{44}{\alpha }^{2}{(}^{n}-{C}_{11}{\xi }^{2}-\theta \text{ }, \\ {F}_{12}^{(i)}(n)=\alpha \xi ({C}_{13}+{C}_{44})(n+{m}_{i})\text{ }, \\ {F}_{21}^{(i)}(n)=-\alpha \xi ({C}_{13}+{C}_{44})(n+{m}_{i})\text{ }, \\ {F}_{22}^{(i)}(n)={\alpha }^{2}{C}_{33}{(}^{n}-{C}_{44}{\xi }^{2}-\theta \text{ }, \end{array} \right\}$

(21)

$\left. \begin{array}{l} {\overline{F}}_{1}^{(i)}(n)=\left[{\alpha }^{2}{C}_{44}(n+{m}_{i}-1)(n+{m}_{i})-{C}_{11}{\xi }^{2}\right]{a}_{n-1}^{(i)}+\\ \alpha \xi \left[{C}_{13}(n+{m}_{i}-1)+{C}_{44}(n+{m}_{i})\right]{b}_{n-1}^{(i)}\text{ }, \\ {\overline{F}}_{2}^{(i)}(n)=-\alpha \xi \left[{C}_{13}(n+{m}_{i})+{C}_{44}(n+{m}_{i}-1)\right]{a}_{n-1}^{(i)}+\\ \left[{\alpha }^{2}{C}_{33}(n+{m}_{i}-1)(n+{m}_{i})-{C}_{44}{\xi }^{2}\right]{b}_{n-1}^{(i)}\text{ }, \end{array} \right\}$

(22)

由此， ${\tilde {{\bar u}_{r1}}}$ ， ${\tilde {{\bar u}_{z0}}}$ 的通解公式可以表示为

$\left. \begin{array}{l} {\tilde{\overline{u}}}_{r1}={\displaystyle {\sum }_{i=1}^{4}{A}_{i}{\displaystyle {\sum }_{n=0}^{\infty }{a}_{n}^{(i)}{\psi }^{n+{m}_{i}}, }}\\ {\tilde{\overline{u}}}_{z0}={\displaystyle {\sum }_{i=1}^{4}A{\displaystyle {\sum }_{n=0}^{\infty }{b}_{n}^{(i)}{\psi }^{n+{m}_{i}}, }} \end{array} \right\}$

(23)

式中，参数A₁，A₂，A₃，A₄为与边界条件相关的待定常数。

将式（23）所示的通解代入式（3），（4），可进一步获得路基动力响应通解表达式： $\tilde {\bar S}(\psi ) =$ $M(\xi , s, \psi )w$ 。其中，w=[A₁, A₂, A₃, A₄]^T，矩阵 $M(\xi , s, \psi )$ 元素如下所示：

$\left. \begin{array}{l} {M}_{1i}(\psi )={\displaystyle {\sum }_{n=0}^{\infty }{a}_{n}^{(i)}{\psi }^{n+{m}_{i}}\text{ }, }\\ {M}_{2i}(\psi )={\displaystyle {\sum }_{n=0}^{\infty }{b}_{n}^{(i)}{\psi }^{n+{m}_{i}}\text{ }, }\\ {M}_{3i}(\psi )=-{\displaystyle {\sum }_{n=0}^{\infty }{C}_{44}{E}_{\text{v}N}(\psi )\times }\left[\alpha (n+{m}_{i}){a}_{n}^{(i)}\text{+}\xi {b}_{n}^{(i)}\right]{\psi }^{n+{m}_{i}}\text{ }, \\ {M}_{4i}(\psi )={\displaystyle {\sum }_{n=0}^{\infty }{E}_{\text{v}N}(\psi )\times }\left[{C}_{13}\xi {a}_{n}^{(i)}-\alpha {C}_{33}(n+{m}_{i}){b}_{n}^{(i)}\right]{\psi }^{n+{m}_{i}}\text{ }。 \end{array} \right\}$

(24)

2.3 刚度矩阵构建

为便于后续推导，对路面结构通解表达式中的系数矩阵M(ξ, s)、状态向量 $\tilde {\bar S}(z)$ 以及系数矩阵w进行分块。 $M(\xi , s) = \left[ {\begin{array}{*{20}{c}} {{D_{\text{d}}}}&{{D_{\text{u}}}} \\ {{S_{\text{d}}}}&{{S_{\text{u}}}} \end{array}} \right]$ ； $\tilde {\bar S}(z) = {\left[ {\begin{array}{*{20}{c}} {{{\tilde {\bar S}}_1}(z)}&{{{\tilde {\bar S}}_2}(z)} \end{array}} \right]^{\text{T}}}$ ， ${\tilde {\bar S}_1}(z) = {\begin{array}{*{20}{c}} {[{{\tilde {\bar u}}_{r{\text{1}}}}(z)}&{{{\tilde {\bar u}}_{z0}}(z)]} \end{array}^{\text{T}}}$ ， ${\tilde {\bar S}_2} = {\begin{array}{*{20}{c}} {[{{\tilde {\bar \tau} }_{zr1}}(z)}&{{{\tilde {\bar \sigma} }_{z0}}(z)]} \end{array}^{\text{T}}}$ ；w=[w₁, w₂]^T，w₁=[A₁, A₂]^T，w₂=[A₃, A₄]^T。取第n(1 < n≤N-1)层路面结构，并令 ${h_n} = {z_n} - {z_{n - 1}}$ ，可得

$\left[ {\begin{array}{*{20}{c}} {{{\tilde {\bar S}}_1}\left| {_n} \right.({z_{n - 1}})} \\ {{{\tilde {\bar S}}_1}\left| {_n} \right.({z_n})} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{D_{{\text{d}}n}}{\varLambda ^{ - 1}}({h_n})}&{{D_{{\text{u}}n}}} \\ {{D_{{\text{d}}n}}}&{{D_{{\text{u}}n}}{\varLambda ^{ - 1}}({h_n})} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{w_1}} \\ {{w_2}} \end{array}} \right] ,$

(25)

$\left[ {\begin{array}{*{20}{c}} {{{\tilde {\bar S}}_2}\left| {_n} \right.({z_{n - 1}})} \\ {{{\tilde {\bar S}}_2}\left| {_n} \right.({z_n})} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{S_{{\text{d}}n}}}&{{S_{{\text{u}}n}}{\varLambda ^{ - 1}}({h_n})} \\ {{S_{{\text{d}}n}}{\varLambda ^{ - 1}}({h_n})}&{{S_{{\text{u}}n}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{w_1}} \\ {{w_2}} \end{array}} \right] ,$

(26)

式中， ${\varLambda ^{ - 1}}(z) = {\text{diag}}({{\text{e}}^{ - az}}, {{\text{e}}^{ - bz}})$ 。引入层间连续条件，并结合式（25），（26），可获得相邻层间的单元刚度矩阵：

$\left[ {\begin{array}{*{20}{c}} {{{\tilde {\bar S}}_2}({z_{n - 1}})\left| {_n} \right.} \\ { - {{\tilde {\bar S}}_2}({z_n})\left| {_{n + 1}} \right.} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\bar K_1^n}&{\bar K_2^n} \\ { - \bar K_3^n}&{ - \bar K_4^n} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde {\bar S}}_1}({z_{n - 1}})\left| {_n} \right.} \\ {{{\tilde {\bar S}}_1}({z_n})\left| {_{n + 1}} \right.} \end{array}} \right] ,$

(27)

式中，

$\left. \begin{array}{l} {\overline{K}}_{1}^{n}=({D}_{\text{d}n}{S}_{\text{u}n}-{D}_{\text{u}n}{S}_{\text{d}n}){\varLambda }^{-1}({h}_{n})\varPhi \text{ }, \\ {\overline{K}}_{2}^{n}=(-{D}_{\text{d}n}{S}_{un}{\left[{\varLambda }^{-1}({h}_{n})\right]}^{2}+{D}_{\text{u}n}{S}_{\text{d}n})\varPhi \text{ }, \\ {\overline{K}}_{3}^{n}=\{{D}_{\text{d}n}{S}_{\text{u}n}\text{ }+{\tilde{F}}_{n}{S}_{\text{d}n}{S}_{\text{u}n}-({D}_{\text{u}n}{S}_{\text{d}n}+{\tilde{F}}_{n}{S}_{\text{d}n}{S}_{\text{u}n}){\left[{\varLambda }^{-1}({h}_{n})\right]}^{2}\}\varPhi \text{ }, \\ {\overline{K}}_{4}^{n}=({D}_{\text{u}n}{S}_{\text{d}n}-{D}_{\text{d}n}{S}_{\text{u}n}){\varLambda }^{-1}({h}_{n})\varPhi \text{ }, \end{array} \right\}$

(28)

式中， $\varPhi = \left\{ {{D_{{\text{d}}n}}{D_{{\text{u}}n}}} \right. + {D_{{\text{u}}n}}{\tilde F_n}{S_{{\text{d}}n}} - ({D_{{\text{d}}n}}{D_{{\text{u}}n}} + {D_{{\text{d}}n}}{\tilde F_n}{S_{{\text{u}}n}}) \times$ ${\left. {{{\left[ {{\varLambda ^{ - 1}}({h_n})} \right]}^2}} \right\}^{ - 1}}$ 。引入边界条件 $\tilde {\bar S}\left| {_N} \right.(z \to \infty ) = {{\text{0}}_{4 \times 1}}$ ，并结合路基通解公式，可获得路基表面状态向量满足下式：

${\tilde {\bar S}_2}({z_{N - 1}})\left| {_N} \right. = \bar K_1^N{\tilde {\bar S}_1}({z_{N - 1}})\left| {_N} \right. ,$

(29)

式中， $\bar K_1^N = {S_{{\text{d}}N}} \cdot D_{{\text{d}}N}^{ - 1}$ ， ${D_{{\text{d}}N}} = \left[ {\begin{array}{*{20}{c}} {{M_{11}}(\Xi )}&{{M_{12}}(\Xi )} \\ {{M_{21}}(\Xi )}&{{M_{22}}(\Xi )} \end{array}} \right]$ ， ${S_{{\text{d}}N}} =$ $\left[ {\begin{array}{*{20}{c}} {{M_{{\text{3}}1}}(\Xi )}&{{M_{{\text{3}}2}}(\Xi )} \\ {{M_{{\text{4}}1}}(\Xi )}&{{M_{{\text{4}}2}}(\Xi )} \end{array}} \right]$ 。

在单元刚度矩阵的基础上，结合式（27），（29），可组建路面结构的总体刚度矩阵方程：

$\left[ {\begin{array}{*{20}{c}} {\bar K_1^1}&{\bar K_2^1}&{}&{} \\ { - \bar K_3^1}&{ - \bar K_4^1 + \bar K_1^2}&{\bar K_2^2}&{} \\ {}&{...}&{...}&{...} \\ {}&{}&{ - \bar K_3^{N - 1}}&{ - \bar K_4^{N - 1} + \bar K_1^N} \end{array}} \right] \times \\ \left[\begin{array}{c}{\tilde{\overline{S}}}_{1}(0)|{}_{1}\\ {\tilde{\overline{S}}}_{1}({z}_{1})|{}_{2}\\ \mathrm{...}\\ {\tilde{\overline{S}}}_{1}({z}_{N-1})|{}_{N}\end{array}\right]=\left[\begin{array}{c}0\\ -\tilde{\overline{p}}(\xi , s)\\ {\text{0}}_{(2N-2)\times 1}\end{array}\right]\text{ },$

(30)

式中， $\tilde {\bar p}(\xi , s)$ 为Hankel-Laplace变换后的表面荷载条件， $\tilde {\bar p}(\xi , s) = \left( {1 + {{\text{e}}^{ - s{T_0}}}} \right){\text{J}_1}({r_0}\xi ){T_0}{{\rm{ \mathsf{ π} }} }{r_0}{p_{\max }}{({s^2}T_0^2 + {{{\rm{ \mathsf{ π} }} }^2})^{ - {\text{1}}}}{\xi ^{ - {\text{1}}}}$ 。求解式（30）所示矩阵方程，即可得到任意层状态向量。

3. 数值计算及验证

3.1 数值计算方法及收敛性分析

前文仅能获得在Hankel-Laplace域内应力响应结果，为得时域结果通常还需进行数值积分逆运算^[23]。考虑到精度要求，采用20节点Gauss插值积分以及Durbin法分别进行Hankel逆变换和Laplace逆变换：

$f(r, {t_l}) = \frac{{2{{\text{e}}^{\gamma l\Delta t}}}}{T}{\rm{R} _e}\left\{ { - \frac{1}{2}{\bar f}(r, z, \gamma )} \right. + \\ {\displaystyle \sum _{k=0}^{{N}_{a}-1}{\displaystyle \sum _{\nu =0}^{L}\left[\overline{f}\left(r, z, \gamma +\text{j}(k+\nu {N}_{\text{a}})\frac{2{\rm{ \mathsf{ π} }} }{T}\right)\cdot {\text{e}}^{j\frac{2{\rm{ \mathsf{ π} }} }{{N}_{\text{a}}}kl}\right]\}\text{ }}},$

(31)

$\overline{f}(r, s)\approx {\displaystyle \sum _{i=1}^{{N}_{0}}{\displaystyle \sum _{k=1}^{20}{A}_{k}\left[{x}_{ik}\tilde{\overline{f}}\left({x}_{ik}, z, s\right)\text{J}({x}_{ik}r)\right]}}\text{ }。$

(32)

式中 $f(r, t), {\bar f}(r, s), \tilde {\bar f}(\xi , s)$ 为时域、Laplace域以及Hankel-Laplace域内函数；N₀为高斯积分分段总数，N₀=ceil(χ/ΔL)，χ为高斯插值积分计算上限，ΔL为高斯积分的分段长度，ceil(∙)为向上取整计算；A_k为高斯积分系数；x_ik为高斯积分第i段第k个积分点，x_ik=(i-1)∙ΔL+x_k，x_k为20节点高斯积分点；T是总求解时间，T=60 ms；N_a为总求解步，N_a=60；Δt为求解步长，Δt=1 ms，T=N_a×Δt；t_l为求解时刻，t_l=l∙T/N_a，l为求解序列；L和 $\gamma$ 为计算参数，取L=ceil(200/ N_a)=4， $\gamma$ =10/T=166.6；j为单位虚数，j²=-1。

在计算时，需进一步确定式（32）中高斯数值积分上限χ，分段长度ΔL，以及式（24）中无穷级数上限n₀。如表 1所示，引入Burgers模型（图 2）来表征沥青面层黏弹性，其Laplace域内表达式如式（33）所示，表中参数取值参考文献[34]。荷载为半正弦波形式，峰值p_max=707 kPa，持续时间T₀=30 ms，作用半径0.15m。分三步进行收敛分析；首先，在ΔL=5，n₀=100基础上，χ分别取50，100，300，500；其次，在χ=300，n₀=100基础上，ΔL分别取2，5，10，15，30；最后，在ΔL=5以及χ=300基础上，n₀分别取20，40，60，80，100，200，结果如图 3所示。考虑到精度和效率，本文取χ=300，ΔL=5，n₀=100。

${E_{\text{v}}}(s) = \frac{{({E_1}{\eta _2} + {\eta _1}{\eta _2}s){E_2}s}}{{{E_1}{E_2} + ({E_1}{\eta _2} + {E_2}{\eta _1} + {E_2}{\eta _2})s + {\eta _1}{\eta _2}{s^2}}} ,$

(33)

表 1 路面结构层材料参数

Table 1. Material parameters of pavement structure

层位	垂直模量/MPa	n_e(E_h/E_v)	?_v(?_h)	密度/(kg∙m^-3)	厚度/m	层间滑移系数α_x
面层	Burgers 模型	0.2	0.25	2400	0.18	0.8
基层	910	0.4	0.30	2300	0.35	0.7
底基层	350	0.5	0.35	2200	0.20	0.0
路基	300~ 200e^-0.2z	1.0	0.35	1600	∞	—

下载: 导出CSV

| 显示表格

图 2 Burgers模型

Figure 2. Burgers model

下载: 全尺寸图片幻灯片

图 3 数值计算收敛研究

Figure 3. Researches on convergence of numerical computation

下载: 全尺寸图片幻灯片

式中，E₁，E₂，η₁，η₂为Burgers模型参数，s为Laplace变换积分变量，参数取值如表 1所示。

3.2 解析解验证

应用ABAQUS以及已有文献结果进一步验证解析解。如图 4所示，构建4层轴对称路面结构有限元模型，其路面厚度、密度、泊松比取值均与表 1对应参数相同，通过改变接触面上的水平向摩擦系数来实现层间非完全连续；模型单元类型均为CAX4环单元。网格采用渐变划分方式，网格尺寸从左到右为0.03~0.5 m，从上到下为0.06~0.5 m。采用隐式动力求解器，求解时间为60 ms，计算步长1 ms。另外，采用UMAT子程序来实现路基模量沿深度方向的非均匀分布，采用2阶Prony级数（式（34））以模拟Burgers四元黏弹性模型。分别采用ABAQUS以及MATLAB计算如表 2所示的2种情况下r为0.0，0.3 m处的弯沉。在计算时，设定路基模量E_v=E_v0+1-e^-5z来似描述匀质模量场，刚性基岩通过设置6000 MPa的路基来模拟。除有限元验证外，选用文献[25，28]进行对比验证。有限元结果如图 5所示，文献对比如图 6所示。由图 5可见，层间接触条件与路基模量非均匀分布均能产生波反射现象，这更符合FWD实际测试结果；且ABAQUS完成计算需花费2~3 min，而解析解在90 s以内即可完成计算（intel i7 8750H，8G）。由图 6可见，退化后的解析解与文献[25，28]结果高度一致。因而，所推导解析解具有较高的精度以及计算效率。

${G_{\text{R}}}(t) = {G_{\text{0}}}\sum\nolimits_{k = 1}^2 {g_k^{\text{p}}{{\text{e}}^{ - t/\tau _k^G}}} ,$

(34)

图 4 有限元与荷载模型

Figure 4. Finite element model and load model

下载: 全尺寸图片幻灯片

表 2 有限元计算参数取值

Table 2. Finite element parameters

计算情况	n_e1 (E_h/E_v)	$\alpha _x^1$	$\alpha _x^2$	$\alpha _x^3$	路基模量E_v /MPa
情况1	1.0	0	0	0	300~200e^-0.2z
情况2	1.0	0.999	0.999	0.0	300~200e^-0.2z
注：表格中未说明参数取值均与表 1相同。

下载: 导出CSV

| 显示表格

图 5 有限元计算结果对比

Figure 5. Comparison of calculated results

下载: 全尺寸图片幻灯片

图 6 文献对比验证

Figure 6. Comparative verification

下载: 全尺寸图片幻灯片

式中， ${G_{\text{0}}}$ 为初始剪切模量， $g_k^{\text{p}}$ 为级数系数， $\tau _k^{\text{G}}$ 为松弛时间。结合，并结合松弛模量等效原则^[35]，有 ${G}_{0}=920\text{ MPa}，{g}_{1}^{\text{p}}=0.23，{g}_{2}^{\text{p}}=0.77，{\tau }_{1}^{\text{G}}=0.67\text{ s}，{\tau }_{\text{2}}^{\text{G}}=$ $0.09{\text{ s}}$ 。

4. 参数讨论及分析

在参数收敛分析以及可靠性验证基础上，本部分主要对路面结构层层间滑移系数 ${\alpha _x}$ 、路面结构层模量比n_e、路基模量非均匀分布系数?和Ξ、路基模量比n_e4等4个因素开展参数讨论及分析。由于篇幅有限，这里仅对路表弯沉进行讨论，而不涉及剪应力和应变，这并不意味所推导解析解仅能用于计算路表弯沉。计算中未特别说明的参数取值均同表 1。

4.1 路面结构层层间滑移系数 ${\alpha _x}$ 影响

设定面层底部、基层底部滑移系数 ${\alpha _x}$ 分别为0.99，0.9，0.8，0.7，0.6，0.5，0.4，0.3，0.2，0.1，0.0，并分别计算其表面弯沉响应。在此基础上，绘制如图 7所示的层底滑移系数–弯沉峰值云图以及弯沉盆曲线。从图 7中可见，路面结构层间滑移系数对弯沉具有显著影响，其影响可达63%。如图 7（a）所示，面层底部和基层底部滑移系数的影响基本相当，云图等高线基本呈现环形分布。如图 7（b）所示，层间接触条件主要影响r < 1.2 m范围内的弯沉盆曲线。

图 7 滑移系数α_x影响

Figure 7. Influences of slip coefficient on deflection

下载: 全尺寸图片幻灯片

4.2 路面结构层模量比n_e影响

设定面层、基层、底基层模量比n_e为0.05~1.00，并分别计算其表面弯沉响应。为便于后续分析，在计算时，假定基层模量比n_e2与底基层模量比n_e3相等，并绘制如图 8所示的模量比–弯沉峰值云图以及弯沉盆曲线。从图 8中可见，路面结构模量比同样对路表弯沉具有显著影响，其影响最大可达70%以上。与滑移系数影响相似，模量比主要影响r < 1.2 m范围内的弯沉盆曲线。

图 8 模量比n_e影响

Figure 8. Influences of modulus ratio on deflection

下载: 全尺寸图片幻灯片

4.3 路基结构非均匀分布系数?以及Ξ影响

值得注意的是，当Ξ < 0.5同时? < 0.3时，矩阵 $M(\xi , s, \psi )$ [式(24)]易接近奇异矩阵而导致计算失败。因而，在路基顶面模量E_v0=100 MPa条件下，分别假定Ξ=0.5~0.8，?=0.3~10.0，结果如图 9，10所示。如图所示，相比路面结构参数，路基模量非均匀分布对弯沉盆的影响范围更大且更为均匀，且对荷载中心处弯沉峰值影响可达40%以上。

图 9 路基变模量参数–弯沉峰值云图(r=0.0 m)

Figure 9. Cloud map of non-uniform distribution coefficient and deflection peak of subgrade (r=0.0 m)

下载: 全尺寸图片幻灯片

图 10 弯沉盆曲线(?, Ξ影响)

Figure 10. Curves of deflection basin under influence of ?, Ξ

下载: 全尺寸图片幻灯片

4.4 路基模量比n_eN影响

在前述路基非均匀分布参数取值基础上，假定路基模量比n_eN分别为1.0，1.5，2.0，2.5和3.0。结果如图 11，12所示。由图 11可见，路基模量比n_eN与路表弯沉峰值间存在较为明显的线性关系，且斜率基本不受路基模量非均匀分布系数?和Ξ影响。由图 12可见，路基模量比对弯沉盆的影响与非均匀分布系数类似，能影响整条弯沉盆曲线，但影响不及其他因素，其对弯沉峰值影响基本在10%以内。

图 11 路基模量比

${n_{eN}}$ 影响探究(r=0.0 m)

Figure 11. Researches on influence of subgrade modulus ratio

${n_{eN}}$ (r=0.0 m)

下载: 全尺寸图片幻灯片

图 12 弯沉盆曲线(路基模量比影响)

Figure 12. Curves of deflection basin under influence of subgrade modulus ratio

下载: 全尺寸图片幻灯片

5. 结论

本文推导了FWD作用下轴对称路面结构动力响应解析解，并考虑了路面结构所存在的黏弹性、横观各向同性以及层间接触状态，以及路基模量沿深度方向的非线性分布特征，得到3点结论。

（1）对于路面结构，层间滑移系数α_x以及模量比n_e均可对路表弯沉峰值具有显著影响，其影响均可达60%以上，但主要影响r < 1.2 m范围内的弯沉盆曲线，在路面力学分析中应充分考虑。

（2）对于路基结构，随着非均匀分布系数?，Ξ以及模量比n_eN的增大，表面弯沉均呈现下降趋势。其中，路基模量比n_eN与路表弯沉峰值间存在较为明显的线性关系，且斜率基本不受路基模量非均匀分布的影响。相比路面结构参数，路基结构参数对弯沉盆的影响范围更大。

（3）对于解析解及计算方法，在高斯插值积分上限取χ=300，分段长度ΔL=5，级数项数n₀=100时具有较高的计算精度及计算效率，可作为路面结构设计及FWD等原位无损检测设备的高精度力学计算模型。

图 1 PFC中可用于剪切带观测的几种方式

Figure 1. Some approaches for observing shear band

下载: 全尺寸图片幻灯片

图 2 通过滚球法实现的“α-shape”算法

Figure 2. α-shape algorithm accomplished by rolling ball method

下载: 全尺寸图片幻灯片

图 3 颗粒级配曲线

Figure 3. Grain-size distribution curve of particles

下载: 全尺寸图片幻灯片

图 4 删除转动角较小的颗粒后进行多次“降噪”

Figure 4. Multiple times of "denoising" after deleting soil particles with relatively small rotation

下载: 全尺寸图片幻灯片

图 5 搜索基线与搜索基点的移动

Figure 5. Searching baseline of shear zone and moving of base point

下载: 全尺寸图片幻灯片

图 6 剪切带基线的搜索过程

Figure 6. Searching process for baselines of shear band

下载: 全尺寸图片幻灯片

图 7 基于“α-shape”算法的剪切带边界的识别

Figure 7. Identification of shear band boundary based on 'α-shape' algorithm

下载: 全尺寸图片幻灯片

图 8 3种剪切带类型

Figure 8. Three types of shear bands

下载: 全尺寸图片幻灯片

图 9 5组不同中值粒径试样的剪切带及识别结果

Figure 9. Shear bands and identificated results from five groups of sand samples with different median diameters

下载: 全尺寸图片幻灯片

图 10 5组不同中值粒径试样的单一型剪切带及识别结果

Figure 10. Single shear bands and identificated results from five groups of sand samples with different median diameters

下载: 全尺寸图片幻灯片

表 1 双轴压缩模型的参数取值

Table 1 Parameters for the biaxial compression model

物理意义	数值	物理意义	数值
试样大小	5 cm×10 cm	接触刚度k_n/k_s	2×10⁷ N/m
颗粒密度	2.6×10³ kg/m³	轴向应变速率	5%/min
孔隙率	0.15	目标轴向应变	10%
围压	100 kPa	阻尼系数	0.7
颗粒摩擦系数	0.5	墙壁摩擦系数	0

下载: 导出CSV

表 2 5组试样PFC^2D双轴压缩试验及结果

Table 2 Identificated results of five biaxial loading tests

d_med /mm	识别宽度/mm	相对宽度	网格法宽度/mm	相对误差/%	孔隙率场宽度/mm	相对误差/%
6	8.0	13.3	7.2	-10	11.7	47
7	8.8	12.5	9.1	4	12.6	44
9	9.8	10.9	9.6	-2	13.3	36
12	11.7	9.8	12.6	7	16.2	38
15	13.5	9.0	15.0	11	17.8	32

下载: 导出CSV

[1]	ALSHIBLI K A, STURE S. Shear band formation in plane strain experiments of sand[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(6): 495–503. doi: 10.1061/(ASCE)1090-0241(2000)126:6(495)
[2]	王冬勇. 基于二阶锥规划有限元法的岩土体稳定性及应变局部化研究[D]. 北京: 北京交通大学, 2020. WANG Dong-yong. Stability and Strain Localization Analysis of Geotechnical Problems Based on Second-Order Cone Programming Optimized Finite Element Method[D]. Beijing: Beijing Jiaotong University, 2020. (in Chinese)
[3]	蒋明镜, 李秀梅. 双轴压缩试验中砂土剪切带形成的离散元模拟分析[J]. 山东大学学报(工学版), 2010, 40(2): 52–58. https://www.cnki.com.cn/Article/CJFDTOTAL-SDGY201002011.htm JIANG Ming-jing, LI Xiu-mei. DEM simulation of the shear band of sands in biaxial tests[J]. Journal of Shandong University (Engineering Science), 2010, 40(2): 52–58. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SDGY201002011.htm
[4]	TANG H X, DONG Y, WANG T, et al. Simulation of strain localization with discrete element-Cosserat continuum finite element two scale method for granular materials[J]. Journal of the Mechanics and Physics of Solids, 2019, 122: 450–471. doi: 10.1016/j.jmps.2018.09.029
[5]	TIAN J Q, LIU E L, HE C. Shear band analysis of granular materials considering effects of particle shape[J]. Acta Mechanica, 2020, 231(11): 4445–4461. doi: 10.1007/s00707-020-02771-y