桐油对固化土搅拌和易性及强度的影响研究

刘行; 邓婷婷; 邓永锋; 詹良通; 刘松玉

doi:10.11779/CJGE20221508

桐油对固化土搅拌和易性及强度的影响研究 English Version

刘行^{1, 2,},
邓婷婷²,
邓永锋^2, ,,
詹良通³,
刘松玉²

1.
中国交通建设股份有限公司, 北京 100088
2.
东南大学交通学院, 江苏南京 211189
3.
浙江大学建筑工程学院, 浙江杭州 310058

基金项目:

国家重点研发计划课题项目 2019YFC1806004

国家自然科学基金项目 42272322

详细信息

作者简介:
刘行(1998—)男，硕士，主要从事环境岩土工程方面研究工作。E-mail: liuhang21@ccccltd.cn

通讯作者:
邓永锋, E-mail: noden@seu.edu.cn

中图分类号: TU43
计量
- 文章访问数: 379
- HTML全文浏览量: 39
- PDF下载量: 119
出版历程
- 收稿日期: 2022-12-04
- 网络出版日期: 2024-04-09
- 刊出日期: 2024-03-31

Mixing workability and strength enhancement of cement-stabilized clay with tung oil

1.
China Communications Construction Co., Ltd., Beijing 100088, China
2.
School of Transportation, Southeast University, Nanjing 211189, China
3.
School of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China

摘要

摘要: 搅拌法在穿越黏土地层（塑性指数大于25）施工时，黏土容易黏附在搅拌桨和转杆上，形成团聚包裹，极大地影响搅拌的和易性和固化土工程性质。桐油是一种干性天然植物油，具有较好疏水和润滑性，在传统生土材料中得到了广泛应用。针对水泥固化土搅拌和易性和强度提升问题，系统地研究了桐油改性黏土的液塑限、固化土的经验黏附比和无侧限抗压强度的变化规律，开展了电子显微镜扫描（SEM）、压汞试验（MIP）和傅里叶红外光谱分析（FT-IR）试验，探究了桐油对固化土搅拌和易性和强度提升的微观机理。结果表明，桐油可显著提高固化土的搅拌和易性，最优掺量约为干土质量的4.0%左右；对均匀搅拌的水泥固化土而言，掺入桐油后28 d无侧限抗压强度可提升约1/3。综合对搅拌和易性与强度的提升效果，桐油在黏土地层搅拌施工中具有良好的应用前景，成果可为理论研究和工程实践提供参考。
- 桐油 /
- 固化土 /
- 搅拌和易性 /
- 无侧限抗压强度 /
- 微观机理
Abstract: When passing through the clay stratum with the plastic index more than 25 and the water content less than 50%, the clay is easy to adhere to the mixing blade and rotating rod, indicating significant influences on the mixing workability and engineering performance of the cement-stabilized soil by using the deep mixing method. Tung oil, a kind of dry natural vegetable oil with excellent hydrophobic and lubricating properties, has been widely used to treat the raw soils as the construction materials. To investigate the influences of the tung oil on the mixing workability and strength of cement-stabilized clays, the Atterberg limit, empirical stickiness ratio and unconfined compressive strength are measured. Furthermore, the scanning electron microscopy, mercury intrusion porosimetry and Fourier infrared spectroscopy are involved to explore the microscopic mechanism of the tung oil on the cement-stabilized clays. The results reveal that the incorporation of the tung oil can significantly improve the workability during the mixing process, and the optimal addition of the tung oil accounts for 4.0% of the dry clay. For the homogeneous stabilized clay with the tung oil, the strength can be enhanced for about 30%. The tung oil has a good application prospect for the in-situ stabilization by the deep mixing method, which is an important reference to overcome the mixing difficulty in this kind of clay strata.
- tung oil /
- cement-stabilized clay /
- mixing workability /
- unconfined compressive strength /
- micro-mechanism

HTML全文

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 水泥土搅拌桩法施工中的“糊钻”现象

Figure 1. Clay adhesion to mixing blade during deep mixing

下载: 全尺寸图片幻灯片

图 2 改性黏土界限含水率与桐油掺量的关系

Figure 2. Atterberg limit vs content of tung oil

下载: 全尺寸图片幻灯片

图 3 不同桐油掺量下固化土的搅拌时的黏附情况

Figure 3. Adhesion during agitation

下载: 全尺寸图片幻灯片

图 4 固化土的经验黏附比与桐油（水）的掺量的关系

Figure 4. Empirical viscosity ratio vs content of tung oil (water)

下载: 全尺寸图片幻灯片

图 5 固化土强度和桐油（水）掺量的关系

Figure 5. UCS vs content of tung oil (water)

下载: 全尺寸图片幻灯片

图 6 养护28 d后固化土的SEM扫描电子显微镜图

Figure 6. SEM images after 28 days' curing

下载: 全尺寸图片幻灯片

图 7 养护28 d后固化土的孔径分布密度图

Figure 7. Pore-size distribution after 28 days' curing

下载: 全尺寸图片幻灯片

图 8 养护28 d后固化土的红外光谱图

Figure 8. FT-IR images after 28 days' curing

下载: 全尺寸图片幻灯片

图 9 桐油的固化作用示意图

Figure 9. Stabilization mechanism of tung oil incorporation

下载: 全尺寸图片幻灯片

表 1 土体基本参数

Table 1 Physical parameters of soils

液限 ${w_{\text{L}}}$ /%	塑限 ${w_{\text{P}}}$ /%	天然含水率 $w$ /%	细砂含量/% (75~425 μ ${\rm{m}}$ )	粉粒含量/ % (2~75 μ ${\rm{m}}$ )	黏粒含量/ % (< 2 μ ${\rm{m}}$ )
43.3	25.5	45.4	1.9	70.1	28.0

下载: 导出CSV

表 2 试验配比方案

Table 2 Test schemes

试验组别	水泥掺量/%	桐油掺量/%	总含水率/%	总水灰比
1	8	0	50.0	6.25
2	8	2.0	50.0	6.25
3	8	4.0	50.0	6.25
4	8	6.0	50.0	6.25
5	8	8.0	50.0	6.25
6	12	0	75.0	6.25

下载: 导出CSV

参考文献(26)

[1]	邓婷婷, 刘行, 张峰, 等. 固化/改性土搅拌均匀性测试与评价方法[J]. 东南大学学报(自然科学版), 2022, 52(4): 697-703. doi: 10.3969/j.issn.1001-0505.2022.04.010. DENG Tingting, LIU Hang, ZHANG Feng, et al. Test and evaluation method for mixing uniformity of soil stabilization/modification[J]. Journal of Southeast University (Natural Science Edition), 2022, 52(4): 697-703. (in Chinese) doi: 10.3969/j.issn.1001-0505.2022.04.010.
[2]	陈富, 李海涛. 黄骅港地区深层水泥土搅拌桩施工工艺研究[J]. 岩土工程学报, 2015, 37(增刊1): 156-160. doi: 10.11779/CJGE2015S1030 CHEN Fu, LI Haitao. Construction technology of cement deep mixing piles in Huanghua Port region[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(S1): 156-160. (in Chinese) doi: 10.11779/CJGE2015S1030
[3]	陈晋中, 刘凤翰, 刘松玉. 双向水泥土搅拌桩技术及常见施工问题处理[J]. 建筑技术, 2011, 42(9): 808-810. https://www.cnki.com.cn/Article/CJFDTOTAL-JZJI201109013.htm CHEN Jinzhong, LIU Fenghan, LIU Songyu. Technology of bidirectional deep mixing cement-soil columns and treatment of common construction problems[J]. Architecture Technology, 2011, 42(9): 808-810. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JZJI201109013.htm
[4]	储诚富, 洪振舜, 刘松玉, 等. 用似水灰比对水泥土无侧限抗压强度的预测[J]. 岩土力学, 2005, 26(4): 645-649. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200504030.htm CHU Chengfu, HONG Zhenshun, LIU Songyu, et al. Prediction of unconfined compressive strength of cemented soils with quasi-water-cement ratio[J]. Rock and Soil Mechanics, 2005, 26(4): 645-649. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200504030.htm
[5]	SZYMKIEWICZ F, TAMGA F S, LE KOUBY A, et al. Optimization of strength and homogeneity of deep mixing material by the determination of workability limit and optimum water content[J]. Canadian Geotechnical Journal, 2013, 50(10): 1034-1043. doi: 10.1139/cgj-2012-0327
[6]	SALVATORE E, MODONI G, SPAGNOLI G, et al. Conditioning clayey soils with a dispersant agent for deep soil mixing application: laboratory experiments and artificial neural network interpretation[J]. Acta Geotechnica, 2022, 17(11): 5073-5087. doi: 10.1007/s11440-022-01505-9
[7]	TERASHI M. Theme lecture: deep mixing method-brief state of the art[C]// 14th International Conference on Soil Mechanics and Foundation Engineering. Hamburg, 1997: 2475-2478.
[8]	NOZU M, SAKAKIBARA M, ANH N T. Securing of in-situ cement mixing quality for the expansive soil with the Montmorillonite inclusion[C]// 2015 DFI Deep Mixing Conference. San Francisco, 2015: 845-852.
[9]	SPAGNOLI G, SCHELLER P, DOHERTY P. In situ and laboratory tests on a novel offshore mixed-in-place pile for oil and gas platforms[J]. Journal of Petroleum Science and Engineering, 2016, 145: 502-509. doi: 10.1016/j.petrol.2016.06.027
[10]	LARSSON S. Mixing Processes for Ground Improvement by Deep Mixing[D]. Stockholm: Royal Institute of Technology, 2003.
[11]	CHAN C M. Influence of mix uniformity on the induced solidification of dredged marine clay[J]. Environmental Earth Sciences, 2014, 71(3): 1061-1071. doi: 10.1007/s12665-013-2510-0
[12]	彭小芹, 曹春鹏, 季晓丽, 等. 防水剂对石灰偏高岭土修补砂浆性能的影响[J]. 湖南大学学报(自然科学版), 2016, 43(6): 83-88. https://www.cnki.com.cn/Article/CJFDTOTAL-HNDX201606014.htm PENG Xiaoqin, CAO Chunpeng, JI Xiaoli, et al. Effect of water-repellent admixtures on repair mortars made of lime and metakaolin[J]. Journal of Hunan University (Natural Sciences), 2016, 43(6): 83-88. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HNDX201606014.htm
[13]	唐晓武, 林廷松, 罗雪, 等. 利用桐油和糯米汁改善黏土的强度及环境土工特性[J]. 岩土工程学报, 2007, 29(9): 1324-1329. http://cge.nhri.cn/cn/article/id/12604 TANG Xiaowu, LIN Tingsong, LUO Xue, et al. Strength and geoenvironmental properties of clay improved by tung oil and sticky rice juice[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(9): 1324-1329. (in Chinese) http://cge.nhri.cn/cn/article/id/12604
[14]	LIN H, LIU F Y, LOURENCO S D N, et al. Stabilization of an earthen material with Tung oil: compaction, strength and hydrophobic enhancement[J]. Construction and Building Materials, 2021, 290: 123213. doi: 10.1016/j.conbuildmat.2021.123213
[15]	FANG S, ZHANG H, ZHANG B, et al. A study of Tung-oil-lime putty-A traditional lime based mortar[J]. International Journal of Adhesion and Adhesives, 2014, 48: 224-230. doi: 10.1016/j.ijadhadh.2013.09.034
[16]	TANG X W, YU Y, ZHOU L P, et al. Study on strength and permeability of silt soils improving by tung oil and sticky rice juice[J]. Advances in Civil Engineering, 2020: 8852998.
[17]	ZUMSTEG R, PUZRIN A M. Stickiness and adhesion of conditioned clay pastes[J]. Tunnelling and Underground Space Technology, 2012, 31: 86-96. doi: 10.1016/j.tust.2012.04.010
[18]	BALL R P A, YOUNG D J, ISAACSON J, et al. Research in soil conditioning for EPB tunneling through difficult soils[C]// Rapid Excavation and Tunneling Conference. Las Vegas, 2009: 320-333.
[19]	LIU L, DENG T T, DENG Y F, et al. Stabilization nature and unified strength characterization for cement-based stabilized soils[J]. Construction and Building Materials, 2022, 336: 127544. doi: 10.1016/j.conbuildmat.2022.127544
[20]	唐晓武, 王艳, 林廷松, 等. 桐油和糯米汁改良土体防渗性和耐久性的研究[J]. 岩土工程学报, 2010, 32(3): 351-355. http://cge.nhri.cn/cn/article/id/12414 TANG Xiaowu, WANG Yan, LIN Tingsong, et al. Permeability and durability of soils improved by tung oil and sticky rice juice[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(3): 351-355. (in Chinese) http://cge.nhri.cn/cn/article/id/12414
[21]	ZHANG H Y, ZHU S B, LI M, et al. Water repellency of monument soil treated by tung oil[J]. Geotechnical and Geological Engineering, 2016, 34(1): 205-216. doi: 10.1007/s10706-015-9939-8
[22]	KHAMEHCHIYAN M, CHARKHABI A H, TAJIK M. Effects of crude oil contamination on geotechnical properties of clayey and sandy soils[J]. Engineering Geology, 2007, 89(3/4): 220-229.
[23]	LIU Z B, LIU S Y, CAI Y. Engineering property test of kaolin clay contaminated by diesel oil[J]. Journal of Central South University, 2015, 22(12): 4837-4843.
[24]	LOURENCO S D N, SAULICK Y, ZHENG S, et al. Soil wettability in ground engineering: fundamentals, methods, and applications[J]. Acta Geotechnica, 2018, 13(1): 1-14.
[25]	LIU D, SANDEEP C S, SENETAKIS K, et al. Micromechanical behaviour of a polymer-coated sand[J]. Powder Technology, 2019, 347: 76-84.
[26]	LU Y Q, MILLER J D. Carboxyl stretching vibrations of spontaneously adsorbed and LB-transferred calcium carboxylates as determined by FTIR internal reflection spectroscopy[J]. Journal of Colloid and Interface Science, 2002, 256(1): 41-52.