Analytical solution for dynamic response of asphalt pavement with subgrade modulus varying with depth
-
摘要: 推导了落锤式弯沉仪(FWD)荷载作用下考虑路面结构黏弹性、横观各向同性、层间连续条件以及路基模量沿深度方向非均匀分布的沥青路面动力响应解析解。首先,在考虑横观各向同性轴对称动力基本方程基础上,借助Hankel-Laplace积分变换建立了路面结构部分的常系数常微分方程组以及路基结构部分的变系数常微分方程组;然后,应用Frobenius法以及刚度矩阵法获得了沥青路面动力响应解析解,并编制数值计算程序;最后,结合ABAQUS有限元结果以及文献对比验证了解析解的可靠性,并就路面结构以及路基结构参数进行讨论。结果表明,路面结构中的层间接触条件和横观各向同性能对弯沉产生显著影响,同时路基模量沿深度的非均匀分布对弯沉的影响也不容忽视,建议在实际的FWD道路检测以及路面力学分析中给予充分考虑。Abstract: The analytical solution for the dynamic response of asphalt pavement under falling weight deflectometer (FWD) load is obtained considering the viscoelasticity, transverse isotropy, interlayer contact and non-uniform distribution of subgrade modulus along depth. Firstly, based on the transversely isotropic axisymmetric dynamic equations, the constant coefficient ordinary differential equations for pavement structure and the variable coefficient ordinary differential equations for subgrade are established through the Hankel-Laplace transform. Then, combined with the Frobenius method and the stiffness matrix method, the analytical solution of dynamic response for the asphalt pavement is obtained, and the numerical program is also compiled. Then, the reliability of analytical solution is verified by comparison and the results of finite element method. Finally, the influences of interlayer condition, modulus ratio and distribution of subgrade modulus on the road deflection are analyzed. The results show that the interlayer contact and transverse isotropy in the pavement structure have a significant impact on the road deflection. The influences of the subgrade modulus distribution along depth also cannot be ignored. It is suggested that the above characteristics of pavement structure and subgrade should be fully considered in the FWD tests and mechanical analysis of the pavement.
-
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]对关于各类颗粒土质模量非均匀分布的研究,采用下式表征路基竖向模量沿竖向的分布规律:
EvN(z)=Ev0+(Ev∞−Ev0)(1−e−αz), (1) 式中,EvN(z)为路基竖向模量沿深度方向的表达式,Ev0为路基顶面模量,Ev∞为无限远处模量,Ev∞ > Ev0,α为描述模量沿深度渐变速度的无量纲系数,e为自然常数,e≈2.718,z为距离路基顶面的垂直距离。
1.1 基本方程
柱坐标系下动力平衡方程为
∂σr∂r+∂τzr∂z+σr−σθr=ρ∂2ur∂t2,∂σz∂z+∂τzr∂r+τzrr=ρ∂2uz∂t2。} (2) 式中σr,σθ,σz,τzr为径向、环向、竖向和切向应力分量;r,z为距离坐标原点的水平和垂直距离;ρ为密度;ur,uz为水平以及垂直方向位移;t为时间。考虑到路基竖向模量Ev为关于z的函数,对柱坐标系下物理方程进行如下改写:
[σrσθσzτzr]=Ev[C11C12C130C12C11C130C13C13C330000C44][εrεθεzγzr], (3) 式中,C11, C12, C13, C33, C44为表征材料横观各向同性的材料参数,C11=kne(1−neμ2v),C12=kne(μh+neμ2v),C13=kneμv(1+μv),C33=k(1−μ2h),C44=Gv=0.5Ev⋅ (1+μv)−1,k=[(1+μh)(1−μh−2neμ2v)]−1,ne为模量比,ne=EhE−1v;εr,εθ,εz,γzr为径向、环向、竖向和切向应变分量;Eh,μv分别为垂直方向的模量以及泊松比;Ev,μh分别为水平方向上的模量以及泊松比。柱坐标下的几何方程为
[εrεθεzγzr]=[∂∂r1r0∂∂z00∂∂z∂∂r]T[uruz]。 (4) 1.2 边界条件
通常,FWD荷载可简化为半正弦脉冲荷载。进而,图 1(a)所示的荷载以及位移边界条件可写为
σz(0)=−pmaxsin(πT0t)H(t)H(T0−t)H(r0−r) ,τzr(0)=σr(0)=0 ,σz(∞)=τzr(∞)=σr(∞)=0 ,ur(∞)=uz(∞)=0 ,} (5) 式中,pmax为荷载峰值,T0为冲击时间,r0为作用半径,通常为0.15 m,H(∙)为阶跃函数。此外,由于实际不同路面结构层材料的差异性,不可避免的会出现层间部分甚至完全滑动的情况,这里采用Matsui等[31]所提出的剪切弹簧来模拟路面层间接触状态:
(1−αnx)[ur(n+1)(zn)−urn(zn)]=αnxβnxτzrn(zn)。 (6) 式中下标n为路面结构层层位(1≤n≤N-1);αnx为第n层与第(n+1)层的层间滑移系数,0⩽;其中,当\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)=[ur(z),uz(z),τzr(z), {\sigma _z} (z)]。可获得如下所示的相邻层间状态向量关系:
{F_n} \cdot S\left| {_n} \right.({z_n}) = S\left| {_{n + 1}} \right.({z_n}) , (7) 式中,Fn为第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] ,I2×2为2阶单位矩阵,02×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 ,Jk(∙)为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}}} 表示对变量uz执行Laplce变换以及0阶Hankel变换的结果; {\tilde {{\bar u}_{r{\text{1}}}}} 表示对变量ur执行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) 式中,参数A1,A2,A3,A4为与边界条件相关的待定常数,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=[A1, A2 A3 A4]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) 与路面结构所获得的常微分方程组不同,由于路基竖向模量EvN(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 = Ev0/Ev∞。
为便于表示,后续表达式中均省略变量 \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) 式中,an,bn为级数系数,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) 式中,θ=ρs2/Ev∞。若式(17)成立,则等式两侧ψ相同次幂的系数应对应相等。令n=0,若a0与b0存在唯一关系,参数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),可得an,bn递推关系:
\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) 式中,参数A1,A2,A3,A4为与边界条件相关的待定常数。
将式(23)所示的通解代入式(3),(4),可进一步获得路基动力响应通解表达式:\tilde {\bar S}(\psi ) = M(\xi , s, \psi )w。其中,w=[A1, A2, A3, A4]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=[w1, w2]T,w1=[A1, A2]T,w2=[A3, A4]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域内函数;N0为高斯积分分段总数,N0=ceil(χ/ΔL),χ为高斯插值积分计算上限,ΔL为高斯积分的分段长度,ceil(∙)为向上取整计算;Ak为高斯积分系数;xik为高斯积分第i段第k个积分点,xik=(i-1)∙ΔL+xk,xk为20节点高斯积分点;T是总求解时间,T=60 ms;Na为总求解步,Na=60;Δt为求解步长,Δt=1 ms,T=Na×Δt;tl为求解时刻,tl=l∙T/Na,l为求解序列;L和 \gamma 为计算参数,取L=ceil(200/ Na)=4, \gamma =10/T=166.6;j为单位虚数,j2=-1。
在计算时,需进一步确定式(32)中高斯数值积分上限χ,分段长度ΔL,以及式(24)中无穷级数上限n0。如表 1所示,引入Burgers模型(图 2)来表征沥青面层黏弹性,其Laplace域内表达式如式(33)所示,表中参数取值参考文献[34]。荷载为半正弦波形式,峰值pmax=707 kPa,持续时间T0=30 ms,作用半径0.15m。分三步进行收敛分析;首先,在ΔL=5,n0=100基础上,χ分别取50,100,300,500;其次,在χ=300,n0=100基础上,ΔL分别取2,5,10,15,30;最后,在ΔL=5以及χ=300基础上,n0分别取20,40,60,80,100,200,结果如图 3所示。考虑到精度和效率,本文取χ=300,ΔL=5,n0=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 ne(Eh/Ev) ?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.2z1.0 0.35 1600 ∞ — 式中,E1,E2,η1,η2为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处的弯沉。在计算时,设定路基模量Ev=Ev0+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) 表 2 有限元计算参数取值Table 2. Finite element parameters计算情况 ne1 (Eh/Ev) \alpha _x^1 \alpha _x^2 \alpha _x^3 路基模量Ev /MPa 情况1 1.0 0 0 0 300~200e-0.2z 情况2 1.0 0.999 0.999 0.0 注:表格中未说明参数取值均与表 1相同。 式中, {G_{\text{0}}} 为初始剪切模量, g_k^{\text{p}} 为级数系数, \tau _k^{\text{G}} 为松弛时间。结合图 2,并结合松弛模量等效原则[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} 、路面结构层模量比ne、路基模量非均匀分布系数?和Ξ、路基模量比ne4等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范围内的弯沉盆曲线。
4.2 路面结构层模量比ne影响
设定面层、基层、底基层模量比ne为0.05~1.00,并分别计算其表面弯沉响应。为便于后续分析,在计算时,假定基层模量比ne2与底基层模量比ne3相等,并绘制如图 8所示的模量比–弯沉峰值云图以及弯沉盆曲线。从图 8中可见,路面结构模量比同样对路表弯沉具有显著影响,其影响最大可达70%以上。与滑移系数影响相似,模量比主要影响r < 1.2 m范围内的弯沉盆曲线。
4.3 路基结构非均匀分布系数?以及Ξ影响
值得注意的是,当Ξ < 0.5同时? < 0.3时,矩阵M(\xi , s, \psi )[式(24)]易接近奇异矩阵而导致计算失败。因而,在路基顶面模量Ev0=100 MPa条件下,分别假定Ξ=0.5~0.8,?=0.3~10.0,结果如图 9,10所示。如图所示,相比路面结构参数,路基模量非均匀分布对弯沉盆的影响范围更大且更为均匀,且对荷载中心处弯沉峰值影响可达40%以上。
4.4 路基模量比neN影响
在前述路基非均匀分布参数取值基础上,假定路基模量比neN分别为1.0,1.5,2.0,2.5和3.0。结果如图 11,12所示。由图 11可见,路基模量比neN与路表弯沉峰值间存在较为明显的线性关系,且斜率基本不受路基模量非均匀分布系数?和Ξ影响。由图 12可见,路基模量比对弯沉盆的影响与非均匀分布系数类似,能影响整条弯沉盆曲线,但影响不及其他因素,其对弯沉峰值影响基本在10%以内。
5. 结论
本文推导了FWD作用下轴对称路面结构动力响应解析解,并考虑了路面结构所存在的黏弹性、横观各向同性以及层间接触状态,以及路基模量沿深度方向的非线性分布特征,得到3点结论。
(1)对于路面结构,层间滑移系数αx以及模量比ne均可对路表弯沉峰值具有显著影响,其影响均可达60%以上,但主要影响r < 1.2 m范围内的弯沉盆曲线,在路面力学分析中应充分考虑。
(2)对于路基结构,随着非均匀分布系数?,Ξ以及模量比neN的增大,表面弯沉均呈现下降趋势。其中,路基模量比neN与路表弯沉峰值间存在较为明显的线性关系,且斜率基本不受路基模量非均匀分布的影响。相比路面结构参数,路基结构参数对弯沉盆的影响范围更大。
(3)对于解析解及计算方法,在高斯插值积分上限取χ=300,分段长度ΔL=5,级数项数n0=100时具有较高的计算精度及计算效率,可作为路面结构设计及FWD等原位无损检测设备的高精度力学计算模型。
-
表 1 路面结构层材料参数
Table 1 Material parameters of pavement structure
层位 垂直模量/MPa ne(Eh/Ev) ?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.2z1.0 0.35 1600 ∞ — 表 2 有限元计算参数取值
Table 2 Finite element parameters
计算情况 ne1 (Eh/Ev) 路基模量Ev /MPa 情况1 1.0 0 0 0 300~200e-0.2z 情况2 1.0 0.999 0.999 0.0 注:表格中未说明参数取值均与表 1相同。 -
[1] BILODEAU J P, DORÉ G. Estimation of tensile strains at the bottom of asphalt concrete layers under wheel loading using deflection basins from falling weight deflectometer tests[J]. Canadian Journal of Civil Engineering, 2012, 39(7): 771–778. doi: 10.1139/l2012-063
[2] 马宪永, 全蔚闻, 董泽蛟. 横观各向同性黏弹性沥青路面动力响应解析解[J]. 中国公路学报, 2020, 33(10): 135–145. doi: 10.3969/j.issn.1001-7372.2020.10.008 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 and Transport, 2020, 33(10): 135–145. (in Chinese) doi: 10.3969/j.issn.1001-7372.2020.10.008
[3] 董泽蛟, 全蔚闻, 马宪永, 等. 考虑沥青路面材料参数空间差异性的解析计算及影响分析[J]. 中国公路学报, 2020, 33(10): 91–101. doi: 10.3969/j.issn.1001-7372.2020.10.004 DONG Ze-jiao, QUAN Wei-wen, MA Xian-yong, et al. Analytical solution and effect analysis of asphalt pavement considering the spatial difference of material parameters[J]. China Journal of Highway and Transport, 2020, 33(10): 91–101. (in Chinese) doi: 10.3969/j.issn.1001-7372.2020.10.004
[4] WANG L B, HOYOS L R, WANG J, et al. Anisotropic properties of asphalt concrete: characterization and implications for pavement design and analysis[J]. Journal of Materials in Civil Engineering, 2005, 17(5): 535–543. doi: 10.1061/(ASCE)0899-1561(2005)17:5(535)
[5] TUTUMLUER E, THOMPSON M R. Anisotropic modeling of granular bases in flexible pavements[J]. Transportation Research Record: Journal of the Transportation Research Board, 1997, 1577(1): 18–26. doi: 10.3141/1577-03
[6] AL-QADI I L, WANG H, TUTUMLUER E. Dynamic analysis of thin asphalt pavements by using cross-anisotropic stress-dependent properties for granular layer[J]. Transportation Research Record: Journal of the Transportation Research Board, 2010, 2154(1): 156–163. doi: 10.3141/2154-16
[7] 董泽蛟, 刘美丽, 郑好, 等. 考虑横观各向同性特性的沥青路面动力学分析[J]. 中国公路学报, 2012, 25(5): 18–23. doi: 10.3969/j.issn.1001-7372.2012.05.005 DONG Ze-jiao, LIU Mei-li, ZHENG Hao, et al. Dynamic mechanical analysis of asphalt pavement based on cross-isotropic properties[J]. China Journal of Highway and Transport, 2012, 25(5): 18–23. (in Chinese) doi: 10.3969/j.issn.1001-7372.2012.05.005
[8] 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(3): 275–282. doi: 10.1080/10298436.2014.993187
[9] AI Z Y, LI Z X, CANG N R. Analytical layer-element solution to axisymmetric dynamic response of transversely isotropic multilayered half-space[J]. Soil Dynamics and Earthquake Engineering, 2014, 60: 22–30. doi: 10.1016/j.soildyn.2014.01.010
[10] AI Z Y, YANG J J, LI H T. General solutions of transversely isotropic multilayered media subjected to rectangular time-harmonic or moving loads[J]. Applied Mathematical Modelling, 2019, 75: 865–891. doi: 10.1016/j.apm.2019.07.015
[11] MOUSA M M, ELSEIFI M A, ELBAGALATI O, et al. Evaluation of interface bonding conditions based on non-destructing testing deflection measurements[J]. Road Materials and Pavement Design, 2019, 20(3): 554–571. doi: 10.1080/14680629.2017.1400995
[12] 冉武平, 张玉, 李爽. 沥青路面层间接触状态研究进展[J]. 重庆交通大学学报(自然科学版), 2019, 38(8): 45–52. doi: 10.3969/j.issn.1674-0696.2019.08.08 RAN Wu-ping, ZHANG Yu, LI Shuang. Research progress on interlayer contact state of asphalt pavement[J]. Journal of Chongqing Jiaotong University (Natural Science), 2019, 38(8): 45–52. (in Chinese) doi: 10.3969/j.issn.1674-0696.2019.08.08
[13] 颜可珍, 满建宏, 石挺魏, 等. 考虑层间接触状态的横观各向同性结构动力响应解析解[J]. 湖南大学学报(自然科学版), 2019, 46(11): 97–105. https://www.cnki.com.cn/Article/CJFDTOTAL-HNDX201911011.htm 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. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HNDX201911011.htm
[14] JIANG X, ZENG C, YAO K, et al. Influence of bonding conditions on flexible base asphalt pavement under non-uniform vertical loads[J]. International Journal of Pavement Engineering, 2021, 22(12): 1491–1503. doi: 10.1080/10298436.2019.1697441
[15] GU F, SAHIN H, LUO X, et al. Estimation of resilient modulus of unbound aggregates using performance-related base course properties[J]. Journal of Materials in Civil Engineering, 2015, 27(6): 04014188. doi: 10.1061/(ASCE)MT.1943-5533.0001147
[16] 张军辉, 彭俊辉, 郑健龙. 路基土动态回弹模量预估进展与展望[J]. 中国公路学报, 2020, 33(1): 1–13. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL202001001.htm ZHANG Jun-hui, PENG Jun-hui, ZHENG Jian-long. Progress and prospect of the prediction model of the resilient modulus of subgrade soils[J]. China Journal of Highway and Transport, 2020, 33(1): 1–13. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL202001001.htm
[17] WANG H, LI M Y. Comparative study of asphalt pavement responses under FWD and moving vehicular loading[J]. Journal of Transportation Engineering, 2016, 142(12): 04016069. doi: 10.1061/(ASCE)TE.1943-5436.0000902
[18] 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(4): 490–498. doi: 10.1080/10298436.2017.1309197
[19] ZHANG J H, PENG J H, ZENG L, et al. Rapid estimation of resilient modulus of subgrade soils using performance-related soil properties[J]. International Journal of Pavement Engineering, 2021, 22(6): 732–739. doi: 10.1080/10298436.2019.1643022
[20] ZHANG J H, PENG J H, ZHENG J L, et al. Characterisation of stress and moisture-dependent resilient behaviour for compacted clays in South China[J]. Road Materials and Pavement Design, 2020, 21(1): 262–275.
[21] PENG J H, ZHANG J H, LI J, et al. Modeling humidity and stress-dependent subgrade soils in flexible pavements[J]. Computers and Geotechnics, 2020, 120: 103413. doi: 10.1016/j.compgeo.2019.103413
[22] MUHO E V, BESKOU N D. Dynamic response of an isotropic elastic half-plane with shear modulus varying with depth to a load moving on its surface[J]. Transportation Geotechnics, 2019, 20: 100248.
[23] ZHANG Y Q, GU F, LUO X, et al. Modeling stress-dependent anisotropic elastoplastic unbound granular base in flexible pavements[J]. Transportation Research Record: Journal of the Transportation Research Board, 2018, 2672(52): 46–56.
[24] BEHNKE R, WOLLNY I, HARTUNG F, et al. Thermo-mechanical finite element prediction of the structural long-term response of asphalt pavements subjected to periodic traffic load: tire-pavement interaction and rutting[J]. Computers & Structures, 2019, 218: 9–31.
[25] 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.
[26] 鲁巍巍, 郑健龙. 横观各向同性黏弹性沥青路面的动力响应[J]. 中南大学学报(自然科学版), 2018, 49(4): 964–970. https://www.cnki.com.cn/Article/CJFDTOTAL-ZNGD201804026.htm 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. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZNGD201804026.htm
[27] AI Z Y, LI Z X. Time-harmonic response of transversely isotropic multilayered half-space in a cylindrical coordinate system[J]. Soil Dynamics and Earthquake Engineering, 2014, 66: 69–77.
[28] 张军辉, 范海山, 张石平, 等. 考虑层间接触状态的路面动力响应解析解及参数反演[J]. 中国公路学报, 2021, 34(5): 11–23. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL202105002.htm ZHANG Jun-hui, FAN Hai-shan, ZHANG Shi-ping, et al. Analytical solution for the dynamic responses and parameter inversion of pavement structures considering the condition of interlayer contact[J]. China Journal of Highway and Transport, 2021, 34(5): 11–23. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL202105002.htm
[29] VRETTOS C. Dynamic response of soil deposits to vertical SH waves for different rigidity depth-gradients[J]. Soil Dynamics and Earthquake Engineering, 2013, 47: 41–50.
[30] VRETTOS C. Rectangular footing on soil with depth-degrading stiffness: vertical and rocking impedances under conditional existence of surface waves[J]. Soil Dynamics and Earthquake Engineering, 2014, 65: 294–302.
[31] MATSUI K, MAINA J W, INOUE T. Axisymmetric analysis of elastic multilayer system considering interface slips[C]// International Symposium on Maintenance & Rehabilitation of Pavements & Technological Control Segundo Simposio Sobre Manutencao E Rehabilitacao De Pavimentos E Controle Technologico, 2001, Auburn.
[32] CH V. In-plane vibrations of soil deposits with variable shear modulus: Ⅱ. Line load[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1990, 14(9): 649–662.
[33] CH V. In-plane vibrations of soil deposits with variable shear modulus: Ⅱ. Line load[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1990, 14(9): 649–662.
[34] 任瑞波, 谭忆秋, 张肖宁. FWD动荷载作用下沥青路面层状黏弹体路表弯沉的求解[J]. 中国公路学报, 2001, 14(2): 9–12, 17. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL200102002.htm REN Rui-bo, TAN Yi-qiu, ZHANG Xiao-ning. Solution for solving asphalt pavement multilayered viscoelastic body surface deflection in the FWD dynamic case[J]. China Journal of Highway and Transport, 2001, 14(2): 9–12, 17. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL200102002.htm
[35] 周骜, 谢发祥, 章登精, 等. 基于修正Burgers模型的浇筑式沥青混合料黏弹性参数确定方法[J]. 林业工程学报, 2017, 2(3): 143–149. https://www.cnki.com.cn/Article/CJFDTOTAL-LKKF201703023.htm ZHOU Ao, XIE Fa-xiang, ZHANG Deng-jing, et al. Viscoelastic parameters determination method of pouring type asphalt mixture based on modified Burgers model[J]. Journal of Forestry Engineering, 2017, 2(3): 143–149. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LKKF201703023.htm
-
期刊类型引用(5)
1. 梁智超,李慧文,刘天赐,苏兴茂,侯相琛. 落锤式弯沉仪评估路基模量的试验研究. 公路. 2025(03): 62-69 . 百度学术
2. 刘超超,鲁文龙,郑健龙,吕松涛,张良奇. 考虑路面自重与循环荷载影响的路基结构模量特性研究. 中国公路学报. 2025(05): 11-25 . 百度学术
3. 李梦威,卢正,唐楚轩,胡智,柴少强,刘永. 基于冲击荷载法的碎石土路基力学特性及压实质量评价研究. 岩石力学与工程学报. 2025(06): 1649-1657 . 百度学术
4. 张婷. 基于路基改良设计的土体承载特征及路面结构力学响应研究. 工程技术研究. 2024(02): 22-24 . 百度学术
5. 刘圣洁,赵之仲,李晓波,隋明言,王晓铭,刘桂强,赵瑜隆. 加热碎石封层对沥青路面层间抗剪性能的影响. 山东交通学院学报. 2023(02): 75-80+88 . 百度学术
其他类型引用(5)