横观各向同性混凝土路面层状结构动力响应求解与分析*

武鸿滨，杨林青，张宇

（1. 河南省收费还贷高速公路管理中心航空港管理处，河南郑州450019；2. 中山大学土木工程学院，广东珠海519082）

路基是路面结构的支撑体，经常出现的路面损坏现象大都是路基稳定性差、强度不足，在外荷载作用下导致过量变形所致。而路基的施工质量是获得坚实而又稳定的路基、保证路基路面整体具有良好使用性能的关键。根据交通运输部《交通运输行业发展统计公报》，我国公路总里程由2010 年的400 万km 增加到2018 年的470 万km，高速公路总里程达到了13.6 万km。随着国家政策的不断务实、路网的逐步完善，以及交通运输一体化和服务质量需求的不断提升，公路交通管理将转变为建养并重、管养结合的模式，公路管理部门面临的主要任务将从建设公路网转向养护、改善和改建现有公路网。

由于路面的荷载主要来自于汽车的运动荷载，对道路工程中路面结构的动力反应问题的求解越来越受到研究学者的重视。随着近年来落锤式弯沉仪（Falling Weight Deflectometer，FWD）的广泛应用，为进行路面结构参数的反演计算，出现了很多路面结构弯沉的计算方法。丰晓等［1］基于弹性层状体系理论建立了沥青路面各层应变计算简式。任瑞波等［2］基于传递矩阵方法，提出了FWD动荷载作用下多层粘弹性结构动力响应的解析解，其中方程的逆变换采用了Durbin 的Laplace 逆变换的数值方法［3］。卢正等［4］将车辆、路面和路基结构作为一个整体系统建模，提出了该结构整体耦合动力学分析模型。Picoux 等［5］和Hamim 等［6］采用有限单元法建立了一种脉冲荷载作用下柔性路面动力响应的数值模型。王建章等［7］同样利用ABAQUS软件建立半刚性基层沥青路面结构模型进行半刚性基层受力特点和功能分析。

Li 等［9］采用一种软计算方法求解得到了FWD动荷载作用下路面弯沉响应；Grenier 等［10］和Al-Khoury等［11］利用谱单元法建立了FWD动荷载作用下柔性路面结构响应模型；吴春颖等［12］提出一类用于三维层状体系的谱单元，用于求解三维层状体系受冲击荷载的动力问题。虽然路面结构动力响应的算法多种多样，但这些计算方法或者以静荷载作用下弹性层状理论为基础，或者假设为各向同性属性。而试验表明，由于地基形成原因和路面的碾压施工工艺，路面结构层和地基均存在不同程度的各向异性特性，尤其是水平和竖直向材料参数的不同，因此利用各向同性假设不能得出真实的结构设计参数。

近年来越来越多的学者研究横观各向同性层状地基动力响应问题，较早针对横观各向同性半无限地基的是Stoneley［13］，随后Synge［14］和Buch⁃wald［15］得到了横观各向同性半无限地基在简谐荷载作用下的动力响应，但是这些研究都是求解均质半无限地基的响应，对于层状地基将不适用。Rajapakse和Wang［16-17］得到了横观各向同性半无限地基在内部时间-简谐荷载作用下的格林函数，同时Shodja 和Eskandari［18］得到了横观各向同性半无限地基在轴对称简谐荷载作用下Substrate-coating系统的动力响应，应用这个方法，Khojasteh 等［19］研究横观各向同性半无限地基在内部简谐荷载作用下的动力响应问题。国内艾志勇等［20］提出了一种解析层元法求解3D横观各向同性地基动力响应；Lin 等［21］和韩泽军等［22-23］采用精细积分算法求解得到了横观各向同性层状地基在任意荷载作用下的动力响应，并针对地基的横观各向同性属性进行了参数分析，结果表明：横观各向同性属性对地表动力响应有显著的影响，随着水平向和竖直向材料参数差异的增大，动力响应均表现出与各向同性情况较大的差异，为得到更精确的计算结果，在实际工程中均应该予以考虑。

尽管以往的算法主要针对横观各向同性层状地基动力响应进行求解分析，而应用到道路动力响应求解的较少。Masad 等［24］建立了车辆荷载作用下各向异性路面结构有限元模型；刘能源等［25］考虑了路面结构层的横观各向同性属性，同样采用有限元方法建立沥青加铺层路面三维有限元模型。基于有限元法的数值计算不仅在计算边界采用人工吸收边界近似模拟，降低精度，而且计算效率较低。

本文考虑路面结构层的横观各向同性属性，首先应用傅里叶变换建立频率-波数域内的状态方程，然后利用精细积分方法求解频率-波数域内的动力柔度矩阵，通过傅里叶逆变换得到路面点频率-空间域内的动力刚度矩阵，最后通过有理函数对离散动力刚度矩阵进行拟合，建立混合变量的时域动力方程，求解得到路面弯沉盆位移响应。本文算法在水平向可以得到解析表达，竖直向通过精细积分算法可以达到计算机意义上的精度。数值算例验证了本文算法的准确性，在此基础上针对路面结构和层状地基厚度以及材料属性进行了广泛的参数分析，以期为工程实际提供可靠的理论依据。

1 层状地基波动方程

1.1 问题描述

FWD 通过将标准质量的重锤从一定高度自由落下，形成脉冲式荷载作用于承载板并通过承载板传到沥青路面，导致路面产生瞬时变形，布设在离荷载作用中心点处不同距离的高精度弯沉盆检测计记录下路面形变，这便得到了弯沉盆在动态荷载作用下的弯沉结果，其原理如图1所示。本文采用弹性层状理论体系，路面结构层自上而下可分为面层（沥青层）、基层和地基层，并且层与层之间的交界面水平，厚度为hi= zi- zi-1(i =1，2，…，l)，考虑各自的横观各向同性属性，每层有5个弹性控制参数，分别为水平向和竖直向杨氏模量Ehhi和Ehvi、水平向和竖直向泊松比μhhi和μhvi以及竖直面的剪切模量Ghvi，每层的质量密度为ρi。建立如图所示的坐标系，水平向为x方向，竖直向为z 方向。假设有n 个弯沉盆检测计设置在地表，其中心点位置距离落锤中心的距离分别为li(i =1，2，…，n)。在落锤处(l1= 0)施加荷载幅值为P 的动态荷载，计算所有弯沉盆的位移时程曲线。

图1 FWD原理图Fig. 1 Schematic diagram of the FWD system

1.2 状态方程建立

首先选取任意路面结构层或者地基层，假设其任一点位移、应力和应变向量为：

对于横观各向同性属性，其应力-应变关系为：

其中，D 为弹性矩阵，其与弹性参数之间的关系为：

本文采用Fourier 变换对波动方程进行空间域和波数域的变换，因此定义与位移幅值相关的傅里叶变换为：

2 精细积分方法

很明显方程（14）是一阶常微分方程，它的解是指数形式，直接求解容易导致指数溢出。本文采用精细积分算法进行求解，不仅数值计算稳定，而且能够达到计算机意义上的精度。下面对精细积分算法求解该一阶常微分方程进行简单介绍，详细求解过程可参见文献［27］。

1）选取任意材料属性相同的结构层或者地基层k（k = 1，2，3，…，l），其厚度为hk，将第k 层沿竖直向划分为2N1（N1为正整数）个相同厚度的子层（如图2所示），每个子层厚为η = hk/2N1。接着，采取同样的方法，将每个子层进一步划分为2N2（N2为正整数）个微细薄层，则每个微细薄层厚度为τ = η/2N2。至此，每个细微薄层的厚度已经足够小，其对应的系数矩阵Q(τ)、G(τ)、F(τ)和E(τ)可利用Taylor级数展开求解并作为初始矩阵。

图2 精细积分算法层区间划分Fig. 2 The layer interval division of the precise integra⁃tion method

3 频率-空间域内路面弯沉点刚度矩阵

在路面落锤处施加半径为Δr 的竖直向均布荷载，由于荷载为轴对称情况，因此后续计算在柱坐标系下进行，荷载在柱坐标系下的表达为：

图3 笛卡尔坐标系和柱坐标系的转换Fig. 3 The transformation between the Cartesian and cylindrical coordinate

方程（34）的积分虽然是个双重积分，但是对于ϕ的积分可以解析求解，因此，该方程可以转化为关于波数κ的无限积分。数值算例表明，被积函数关于波数κ有良好的收敛性，无限积分由有限积分代替依然可以得到较高的精度。

在所有弯沉盆处施加荷载求解得到所有监测点的频率域位移响应便可构建系统动力方程为

4 时域动力方程的建立

在得到频率域内离散动力刚度矩阵的基础上，采用Ruge等［28］提出的有理函数进行拟合，即

5 数值算例

5.1 算法验证

对本文的算法进行校验，验证其是否合理。而现有文献较少针对横观各向同性路面结构进行求解，所以本文采用退化后的程序进行求解，与已有关于各向同性路面结构的结果进行对比。假设整个路面结构为各向同性属性，路面结构的材料参数如表1所示。选取包括荷载作用点在内的10个弯沉检测点，用与荷载作用点距离r 表示，分别为0、0.25、0.375、0.50、0.625、0.75、1.00、1.25、1.5 和2.00 m。将本文计算得到的弯沉点最大位移结果与文献［11］采用的谱单元法计算结果进行对比，如图4所示。两者基本完全拟合，从而验证了本文算法的准确性，而常规的谱单元法只能处理各向同性路面结构，对于横观各向同性路面结构将不再适用。

图4 本文算法结果与文献［11］结果比较Fig. 4 The comparison between present results and reference［11］

5.2 边界条件对路面结构动力响应的影响

一般路面结构动力计算问题中均假定地基为均质半无限空间，而忽略地基的层状特性。而有些情况，地基在一定深度存在岩性介质，其对路面结构的动力响应有显著的影响，为分析其对弯沉盆位移响应的影响，选择地基底部为弹性半无限体与刚性基岩两种情况进行对比计算。地基底部为弹性半无限体情况与算例5.1相同；底部为刚性基岩的情况选择如表1所示的路面层状结构材料参数，其中地基层厚为5.0 m，底部为刚性基岩。

表1 层状路面结构材料参数Table 1 The material properties of the layered pavement structure

利用本文算法，求得路面结构施加荷载点2.0 m 范围内的动力响应的最大竖向位移，如图5 所示。由图可见，地基为刚性基岩的情况对路面结构动力响应有显著的影响，虽然最大位移沿距离落锤点的变化趋势基本相同，但是曲线基本发生整体平移。实际工程中，道路设计者应当勘测清楚地基底部工况条件，不同边界条件对路面结构动力响应的影响极为不同。

图5 不同边界条件对路面结构动力响应的影响Fig. 5 The influence of the different boundary conditions on dynamic response of pavement structure

5.3 横观各向同性属性对路面结构动力响应的影响

由于路面碾压和其他因素，实际工程中的道路路面结构并非简单的各向同性属性，而是具有一定程度的各向异性，本文考虑各个结构层的横观各向同性属性，选择不同的横观各向同性系数n=Ehh/Ehv，研究其对路面弯沉值的影响规律。

5.3.1 面层横观各向同性属性为了研究面层的横观各向同性属性对路面结构动力响应的影响，选择如表1 所示的路面层状结构模型，假设基层、底基层和地基的材料参数不变，而面层假定为横观各向同性属性，横观各向同性系数n=Ehh/Ehv分别取为0.5、1.0、1.5 和2.0 四种情况，即保持水平向杨氏模量Ehh不变，根据系数n 选取不同的竖直向杨氏模量。利用本文提出的算法，针对四种情况分别进行计算得到路面层状结构表面弯沉盆的位移响应，将落锤处（r=0.0 m）的位移时程响应以及在r=2.0 m 范围内的最大位移响应绘制成曲线，如图6所示。

图6 面层横观各向同性属性对路面动力响应的影响Fig. 6 The influence of the transversely isotropic property of the surface layer on the dynamic response of the pavement

由图6（a）可知，面层横观各向同性属性对于路面弯沉时程曲线影响较小，对于位移上升段和下降段的影响均可忽略不计，只是对最大位移峰值有略微显著的影响，随着n的增大，其最大位移有增大的趋势，但增大幅值较小。由图6（b）可见，面层横观各向同性属性对于路面最大竖向位移的影响范围也较小，只在0.2 m 范围内有略微明显的影响，随着r 的增大，影响逐渐减小，甚至可以忽略不记。总体而言，面层的横观各向同性属性由于其厚度较小，其对于路面弯沉的影响整体较小。

5.3.2 基层横观各向同性属性为了研究基层的横观各向同性属性对路面结构动力响应的影响，同样选择如表1所示的路面层状结构模型，假设面层、底基层和地基的材料参数不变，而基层假定为横观各向同性属性，横观各向同性系数n=Ehh/Ehv分别取为0.5、1.0、1.5 和2.0 四种情况。利用本文提出的算法，针对四种情况分别进行计算得到路面层状结构表面弯沉盆的位移响应，将落锤处（r=0.0 m）的位移时程响应以及在r=2.0 m 范围内的最大位移响应绘制成曲线，如图7所示。

图7 基层横观各向同性属性对路面动力响应的影响Fig. 7 The influence of the transversely isotropic property of the base course on the dynamic response of the pavement

由图7可知，基层横观各向同性属性对于路面层状结构动力响应的影响比面层显著得多。随着n的增大，路面最大位移逐渐增大，与各向同性属性（n=1.0）相比，横观各向同性情况计算结果均有明显的不同。由图7（b）可见，基层横观各向同性属性的影响范围较面层更大，在r=1.0 m 内最大位移均有差别，且距离落锤位置越近，影响越大。

图8 底基层横观各向同性属性对路面动力响应的影响Fig. 8 The influence of the transversely isotropic property of the subbase on the dynamic response of the pavement

5.3.3 底基层横观各向同性属性底基层作为路基的最底层，其与地基直接相邻，为了研究底基层的横观各向同性属性对路面结构动力响应的影响，假设面层、基层和地基的材料参数不变，而底基层假定为横观各向同性属性，横观各向同性系数n=Ehh/Ehv分别取为0.5、1.0、1.5 和2.0 四种情况。利用本文提出的算法，针对四种情况分别进行计算得到路面层状结构表面弯沉盆的位移响应，将落锤处（r=0.0 m）的位移时程响应以及在r=2.0 m 范围内的最大位移响应绘制成曲线，如图8所示。由图可见，与面层相比，底基层与基层横观各向同性属性对于路面动力位移时程响应的影响基本相同，位移最大值的变化幅度基本相同，但不同的是底基层横观各向同性属性的影响范围比基层更广，由图8 可见，随着n 的变化，最大竖向位移曲线在1.0 m 范围内基本发生了平移，在1.0 ～2.0 m范围内的影响也是很明显的。

5.3.4 同时考虑面层、基层和底基层横观各向同性属性为了研究路面结构层均为横观各向同性属性时对路面结构动力响应的影响，假设面层、基层、底基层以及地基的材料属性均为横观各向同性情况，横观各向同性系数n=Ehh/Ehv分别取为0.5、1.0、1.5 和2.0 四种情况。利用本文提出的算法，针对四种情况分别进行计算得到路面层状结构表面弯沉盆的位移响应，将落锤处（r=0.0 m）的位移时程响应以及在r=2.0 m 范围内的最大位移响应绘制成曲线，如图9所示。由图可见，当同时考虑路面三结构层的横观各向同性属性时，随着n的增大，路面位移响应也相应增大。但是当n小于1 时对位移响应的影响比较明显；当n 大于等于1时，影响比较小。

图9 路面结构层横观各向同性属性对路面动力响应的影响Fig. 9 The influence of the transversely isotropic property of the layered pavement structure on the dynamic response of the pavement

6 结论

本文利用积分变换方法求解得到了频域内弯沉点动力刚度矩阵，进而采用混合变量算法建立了弯沉点动力时域方程，求解得到横观各向同性层状路面结构在动力荷载作用下的竖直向位移响应。在验证本文算法准确性的基础上，针对地基底部边界条件和各个路面结构层的横观各向同性属性进行了参数分析，结果表明：

1）地基底部边界条件对路面结构动力响应有显著的影响，在工程实际中如果地基底部为刚性基岩的情况，为得到更为准确的结果，不能简单的假定地基为均值半无限空间；

2）路面结构层中面层沥青的横观各向同性属性对路面动力响应影响较小，但是基层和底基层的影响比较显著，影响范围也较大，在工程实际中应该予以考虑。