基于半解析有限元的水泥路面荷载响应分析

沈凯仁，陈先华，周 杰，张含宇

(东南大学交通学院，南京210096)

车辆荷载是水泥混凝土路面结构的主要外部作用，荷载引起的力学响应是路面结构设计的重要指标[1]。我国现行规范[2]中，水泥混凝土路面荷载应力的计算，仍主要使用弹性地基上的薄板模型；规范中只给出了临界荷位处的荷载应力计算公式，但未对路面结构荷载应力的分布与变化趋势进行表述。有限元方法通过对路面结构进行单元划分，并逐一定义单元属性，计算各节点的位移和应力，可以有效地反映路面结构的力学响应分布与变化趋势[3]。因此，有限元方法在路面荷载响应分析中不断得到应用，现有研究[4－5]证明了有限元是分析路面力学响应的有效工具。主流的有限元程序(如ABAQUS、ANSYS)分析路面荷载响应时，一般采用三维模型进行模拟，其对专业知识和计算成本要求较高，使其局限于研究领域，难以推广到普遍的工程应用中[6]。因此，有必要结合路面结构的特殊性，对有限元方法进行适当的简化，开发出既能保证计算精度，又能节约计算成本的有限元程序，并进一步将该程序引入实际工程中，以辅佐现行水泥路面设计规范。

将空间问题简化为平面问题是提升有限元运算效率的重要方法之一。例如，挡土墙、隧道、管道等结构，其横截面的几何尺寸、材料属性和外力约束一般沿纵向不发生变化，原本的三维有限元问题可以简化为平面应变问题[7－8]。但对于典型的路面结构，通常各结构层假定为各向同性的均质材料，其横断面的几何尺寸和材料属性沿行车方向不发生变化，但是外部荷载和边界条件沿该方向却有着明显的变化，因此，路面结构问题不能直接简化为平面应变问题[9]。轴对称分析也是简化有限元模型的重要手段之一，但路面荷载的分布情况通常难以满足轴对称模型的要求，因此路面结构问题的轴对称分析方法也受到了诸多限制[10]。基于路面结构的特殊性，可以在横断面上进行二维有限元网格划分，并通过傅里叶级数表示单元节点位移和路面荷载沿行车方向的分量，利用傅里叶级数表示连续函数的能力，节点位移和荷载沿行车方向的变化均得到了良好的表达。同时，路面结构只在横截面进行了二维有限元数值分析，荷载和位移沿行车方向的分量采用傅里叶级数的解析表达，在形成整体平衡方程的过程中，由于傅里叶级数的正交性，整体刚度矩阵变化为对角矩阵，原本的三维有限元问题简化为多个独立的二维平面问题与解析问题的组合，通过叠加二维平面问题的结果，即可得到路面结构真实的受力状态，该方法也被称作半解析有限元方法[11]。而且二维有限元往往只能获得平面内的应力状态，但半解析有限元可以通过选取横断面纵向坐标，获得任意横断面的应力状态，因此它的计算功能与三维有限元更为接近。图1简要对比了半解析有限元与三维有限元的模型建立方法。

已有的工程应用成果[12－16]，均表明半解析有限元在力学分析中具有良好的准确性与计算效率。Liu等[12－13]通过半解析有限元模拟沥青路面弹性层状体系，对比了试验路的路表弯沉值与半解析有限元的数值解，结果具有良好的一致性。Wu等[14]将超声波导波理论与半解析有限元相结合，开发了一种计算列车荷载作用下高铁轨道导波特性的分析程序。Bian等[15－16]通过把Euler-Bernoulli梁模拟的钢轨和半解析有限元模拟的轨下基础相结合，分析了荷载作用下高铁轨下结构的振动响应，分析结果与 ABAQUS的模拟结果吻合良好。基于此，将半解析有限元引入水泥路面荷载响应分析中，开发出准确、高效的有限元分析程序，具有一定的技术可行性和良好的工程应用前景。

图1 半解析有限元与三维有限元模型示意图Fig.1 Diagram of the semi-analytical finite element and 3D finite element model

本文基于半解析有限元原理和MATLAB软件，开发了水泥混凝土路面荷载响应分析程序CP-SAFEM。在下面的章节中，将具体介绍半解析有限元基本理论和CP-SAFEM模型建立的方法，并结合具体工况与现行规范和 ABAQUS进行荷载响应的对比验证，据此证明CP-SAFEM具有良好的计算精度与运行效率；进一步运用CP-SAFEM确定水凝混凝土路面的临界荷位，为CP-SAFEM在具体工程中的应用提供一定的参考。

1 水泥路面半解析有限元的力学公式推导与模型建立

1.1 路面结构的半解析有限元模型

水泥路面半解析有限元模型如图2所示，假定路面结构的几何尺寸和材料属性沿z方向保持不变。在此基础上沿z方向对位移函数和荷载函数进行傅里叶变换，利用傅里叶级数表示连续函数的能力表征位移函数和荷载函数沿z坐标的变化趋势。

图2 水泥路面半解析有限元模型Fig.2 Semi-analytical finite element model of cement pavement

本文采用四节点矩形单元进行单元划分，如图3所示，与二维有限元不同的是，每个节点需要定义三个方向的位移分量。

图3 半解析有限元四节点矩形单元Fig.3 The four-node rectangular unit of semi-analytical finite element

节点位移u的形函数表示成z的傅里叶级数，单元内任一点的位移u定义如下：

式中：k为单元节点编号；uk为节点位移；Nk(x,y,z)为对应的形函数；l为第l项傅里叶级数；L为傅里叶级数的总项数；是xy平面内的节点形函数，与二维平面有限元中使用的形函数相同。同理，位移v和w的表达式与u的位移表达式相似。

荷载函数也可以沿z方向进行相应的傅里叶变化，如式(2)。即若原荷载函数为作用平面xz内的面密度，则变换后的荷载函数为假定z坐标下沿x方向的线密度。

对于路面的边界条件，也可以通过傅里叶变换的奇拓展或偶拓展实现。例如，假定在z=0和z=a的xy平面内仅存在沿z方向的位移，则应对位移u和v进行奇拓展，对位移w进行偶拓展，这样在边界平面内仅存在沿z方向的位移，单元内任一点的位移函数可以表示成下式：

考虑到路面结构三维应力状态下单元应变和位移有如下关系：

将式(3)代入式(4)可得到式(5)。其中，应变-位移矩阵如式(6)。

此外，单元应力与单元应变有如下关系，其中D是各向同性的线弹性材料的本构矩阵。

根据最小势能原理确定单元对整体平衡方程的贡献，单元刚度矩阵和荷载向量表示见式(8)~式(9)：

式中：Nl是表面力作用面的形函数矩阵的傅里叶级数的第l项。同时，式(8)将包含下列积分：

由于积分的正交性，可以得到如下关系：

因此，单元刚度矩阵可表达成式(12)：

式中，Ae为单元积分域。通过叠加法组装整体刚度矩阵和整体荷载向量，根据最小势能原理获得路面结构的平衡方程：

式(13)表明，最终的平衡方程简化为L个独立的方程，可简写为式(14)，其中，K11、F1分别为整体刚度矩阵和荷载向量在傅里叶变换下的第l项。

如果荷载向量的傅里叶变换仅包含特定的某几项，那么只需求解对应的独立方程，这进一步节省了计算成本。在得到节点位移的傅里叶级数后，通过傅里叶逆变换可以得到对应三维空间中的节点位移，进一步根据应变-位移、应力-应变关系可以分析整个路面结构的力学响应。

1.2 层间连接和边界条件的半解析表达

对于普通水泥混凝土路面结构，如图4所示，选取一块水泥混凝土面板作为分析对象，考虑横缝构造和传力杆的传荷作用，传力杆的设置不应妨碍相邻混凝土板的自由伸缩[17]，因此对z=0和z=a的截面仅固定xy方向的位移，z方向自由伸缩。可以通过式(15)的傅里叶变换满足上述边界条件：

图4 普通水泥混凝土路面结构示意图Fig.4 Diagram of cement concrete pavement structure

纵缝边界设置为自由端，同时基层和路基选择与路面板相同的边界条件。此外，在路基底部设置竖向位移为零的法向约束。层间接触面假定为完全连续，即接触面上、下单元共用节点。

1.3 车辆荷载的半解析表达

依照我国路面结构厚度设计规范，本文采用单轴双轮 100 kN作为分析荷载，轮胎接地形状简化为矩形，荷载位置分为板中、纵缝中部和角隅三种情况，如图5所示。对荷载项应进行特定的傅里叶变换以确保在Z域中正确模拟荷载情况，荷载函数的傅里叶变换[18]可以写成：

式中：Pt为荷载压应力；t为轮载个数；Zt1、Zt2分别为轮载起始和结束的z轴坐标。

图5 不同荷载位置示意图Fig.5 Diagram of different load positions

1.4 基于MATLAB的半解析有限元分析程序

本文采用 MATLAB编写力学响应分析程序。通过发挥 MATLAB矩阵运算的优势，快速获得路面力学响应，用于指导层位结构和厚度设计。需要指出，为了避免积分的不可积性，对刚度矩阵和荷载向量分别采用高斯插值求解。同时，在组装整体刚度矩阵的过程中，由于单元形函数需要从局部坐标系过渡到全局坐标系，因此需要借助雅可比矩阵进行坐标变换。

式中：J为雅可比矩阵；|J|为对应的行列式；Wi、Wj为高斯点对应的权系数。

2 CP-SAFEM各模块功能介绍

水泥路面半解析有限元程序由前处理模块、分析模块和后处理模块三个部分组成，软件总体结构如图6所示，各模块又由对应的子程序组成，具体功能如下。

2.1 前处理模块

包括建立路面结构的半解析有限元模型和表达车辆荷载的半解析有限元格式的功能。通过输入模型尺寸、材料参数和荷载参数，定义单元类型和尺寸，实现模型的网格自动划分。进一步自动形成节点编号、单元编号和单元对应节点。不同结构层的层间接触通过相邻单元共用节点模拟为层间完全连续。根据输入的荷载形状与大小，形成在模型空间中的荷载函数。

2.2 分析模块

包括计算分析路面结构的力学响应的功能。依照有限元单元刚度矩阵和单元荷载向量的形成方法，并沿z坐标对位移形函数和荷载形函数进行傅里叶变换，形成单元刚度矩阵和单元荷载向量在局部坐标系下的傅里叶表达。并根据雅各比矩阵和单元在整体中的位置，实现单元刚度矩阵和单元荷载向量的傅里叶表达从局部坐标系下到全局坐标的转换。进一步根据直接刚度法和叠加法逐一装配整体模型中各单元的刚度矩阵和荷载向量。依据最小势能原理求解傅里叶级数各维度下的节点位移，并通过傅里叶逆变换求解出原三维空间中的节点位移。最后根据位移-应变-应力关系求解出整体模型各项力学响应。

图6 CP-SAFEM模块流程图Fig.6 Flow diagram of CP-SAFEM

2.3 后处理模块

包括生成各项响应指标的数据表格、绘制力学响应的结构云图和查询指定节点的力学响应的功能。后处理模块主要是对力学响应进行可视化处理。通过指定横断面的z坐标，可以生成各项力学响应指标与对应单元节点的数据表格；同时也可以通过指定节点坐标查询该节点的各项力学指标。并依据生成的数据表格绘制响应指标沿某一坐标的变化趋势或任意截面的力学响应云图。

3 有限元程序 CP-SAFEM 的合理性验证

结合单层板模型和双层板模型两种常见的水泥路面结构，对比了规范公式、CP-SAFEM 和ABAQUS的荷载应力求解结果；并针对一种典型路面结构，进行了板中荷位作用下 CP-SAFEM 与ABAQUS的各项应力指标对比，进一步比较了两者的精算精度和计算效率。

3.1 单层板模型荷载应力对比结果

对于粒料基层上的单层板水泥路面应力分析，规范中采用弹性地基单层板模型[2]。面层板底面以下部分按弹性地基处理，采用地基当量回弹模量Et表征板下结构层的整体弹性特性，计算公式如式(20)~式(23)：

式中：E0/MPa为路基回弹模量；Ex/MPa为粒料层当量回弹模量；α为回归系数；hx为粒料层总厚度；Ei和hx/m为第i结构层的回弹模量与厚度。

设计轴载在临界荷位处产生的荷载应力σps(路面板层底最大弯拉应力)按式(24)~式(26)计算：

式中：Ps/kN为设计轴载轴重；Ec/MPa、hc/m和vc分别为混凝土面层板的弹性模量、厚度和泊松比；r/m为路面板的相对刚度半径；Dc/(MN·m)为混凝土面层板的截面弯曲刚度。

对于典型的单层板路面结构，结构层材料参数如表1所示，路面单元宽度3.5 m，长度4.5 m，荷载选取100 kN，临界荷位选择路面板纵缝中部[2]，如图7所示。

在CP-SAFEM和ABAQUS中，横截面上采用相同的划分形式，如表2，对于纵向行车方向，ABAQUS划分为50个单元，CP-SAFEM选取40项傅里叶级数。荷载选取单轴双轮100 kN，轮印采用 0.2 m×0.18 m 的矩形模拟，均布压应力为0.694 MPa。通常，对于典型弹性层状体系，傅里叶级数维度选择 40即可保证良好的计算精度[19]，因而本文的傅里叶维度均选择40。

表1 单层板水泥路面各结构层材料参数Table 1 Material parameters of single-layer cement pavement

图7 临界荷位示意图Fig.7 Diagram of the critical load position

表2 路面横断面网格划分规格Table 2 Meshing division of pavement cross section

三种方法的计算结果如图8。对比板底最大弯拉应力，规范计算结果为 1.652 MPa，CP-SAFEM计算结果为 1.643 MPa，ABAQUS计算结果为1.577 MPa。CP-SAFEM 与规范计算结果相差0.54%。对比结果表明，CP-SAFEM 的计算结果与规范、ABAQUS具有良好的一致性。另外，CP-SAFEM不仅求解出最大荷载应力，同时也能绘制不同方向的荷载应力的变化趋势，例如路面板板底的纵向荷载应力远大于横向荷载应力[20]，这对路面板结构和配筋设计具有重要意义。

图8 单层板模型路面层层底拉应力Fig.8 The tensile stress of single-layer plate model pavement

3.2 双层板模型荷载应力对比结果

对于设置了半刚性基层的水泥路面结构，规范中采用弹性地基上的双层板模型，半刚性基层以下部分按弹性地基处理，采用地基当量回弹模量Et表征板下结构层的整体弹性特性，分别求解路面板层底和半刚性基层层底的荷载应力。

设计轴载Ps在临界荷位处产生的路面板层底荷载应力σps(最大弯拉应力)按下式计算：

式中：rg/m为路面板的相对刚度半径；Dc/(MN·m)为混凝土面层板的截面弯曲刚度；Eb/MPa、hb/m和vb分别表示半刚性基层的弹性模量、厚度和泊松比。半刚性基层层底的荷载应力σbps(最大弯拉应力)按式(30)计算：

对于典型的双层板路面结构，结构层材料参数如表3所示，路面单元宽度3.5 m，长度4.5 m，在CP-SAFEM和ABAQUS中，定义相同的横截面网格划分形式，如表4，对于纵向行车方向，ABAQUS划分为50个单元，CP-SAFEM选取40项傅里叶级数。荷载形式与单层板模型中相同。

表3 单层板水泥路面各结构层材料参数Table 3 Material parameters of single-layer cement pavement

表4 路面横断面网格划分规格Table 4 Meshing division of pavement cross section

三种方法的计算结果如图9。对比路面板板底最大弯拉应力，规范计算结果为 1.039 MPa，CP-SAFEM计算结果为0.968 MPa，ABAQUS计算结果为0.947 MPa，CP-SAFEM与规范计算结果相差6.83%。对比半刚性基层板底最大弯拉应力，规范计算结果为0.341 MPa，CP-SAFEM计算结果为0.335 MPa，ABAQUS计算结果为 0.316 MPa，CP-SAFEM与规范计算结果相差1.76%。

图9 双层板模型路面层层底拉应力Fig.9 The tensile stress of double-layer plate model pavement

3.3 板中荷位下的荷载应力对比结果

通过 CP-SAFEM 分析板中荷载作用下水泥混凝土路面力学响应，通过与 ABAQUS分别对比道路表面竖向位移、轮迹下竖向压应力、水泥板底部横向和纵向拉应力四项指标，以此验证CP-SAFEM的正确性和可靠性。各结构层材料参数如表5所示，荷载选择单轮-双轴荷载，轴重100 kN，轮印宽度0.24 m，轮印高度0.2 m，接地压力0.521 MPa。路面宽度取4.8 m，长度取5.0 m。在CP-SAFEM和ABAQUS中，定义相同的结构层，横截面上采用相同的划分形式，如表6，对于纵向行车方向，ABAQUS划分为50个单元，CP-SAFEM选取40项傅里叶级数。同时在CP-SAFEM中选择四节点矩形单元，ABAQUS中选择八节点六面体单元，模型结构如图10、图11所示。

表5 水泥路面各结构层材料参数Table 5 Material parameters of each structural layer of the cement pavement

表6 路面横断面网格划分规格Table 6 Meshing division of pavement cross section

图10 路面结构ABAQUS剖分模型Fig.10 Pavement structure model of ABAQUS

图11 路面结构CP-SAFEM剖分模型Fig.11 Pavement structure model of CP-SAFEM

对于典型的水泥路面结构，最大拉应力发生在轮载下水泥板板底，需对水泥板板底横向和纵向拉应力进行研究；竖向压应力是保证路基稳定性的重要指标。因此本文选择上述四项指标作CP-SAFEM与ABAQUS的对比验证，验证结果如图12所示。从图12中可以发现，CP-SAFEM 的计算结果与ABAQUS均保持良好的一致性，关于各项指标的最大值的具体对比结果见表7。通过对比发现各项指标与 ABAQUS误差均在 5%以内，说明自编程序CP-SAFEM的计算精度满足要求。理论上提高傅里叶变换维度，计算精度会进一步提高。

3.4 CP-SAFEM计算效率的验证

图12 各项力学响应指标对比图Fig.12 Comparison of various mechanical indicators

为了评估CP-SAFEM的计算效率，CP-SAFEM和ABAQUS均在intel core i5-8500 CPU、8G内存上运行，单元、节点个数和计算时间如表8所示。通过对比可知，CP-SAFEM 需要的单元个数仅为ABAQUS的1/50，计算时间仅为ABAQUS的1/6，如果开启 MATLAB多线程优化，计算时间将会进一步缩减[21]。

表7 各项应力指标最大值对比Table 7 Maximum comparison of various stress indicators

表8 CP-SAFEM与ABAQUS单元、节点和运行时间对比Table 8 Comparison of CP-SAFEM and ABAQUS units,nodes and operation time

通过本节分析可以发现，半解析有限元不仅在计算精度上基本达到 ABAQUS水准，同时利用傅里叶变换和正交函数的优势，可以极大地减少计算成本。随着分析对象结构尺寸的增加和材料复杂程度的加深，这一优势将变得更加突出。

3.5 不同荷位下路面板弯拉应力对比分析

在水泥混凝土路面路面荷载应力分析中，通常选择使面层板产生最大应力的荷载位置作为应力计算的临界荷位，临界荷位一般发生在水泥路面板板底。现通过半解析有限元程序CP-SAFEM，采用3.4节水泥混凝土路面结构模型，纵缝边界不考虑拉杆，分别考虑板中、纵缝中部、角隅三种荷位，作路面板层底横向和纵向拉应力的对比分析，据此寻找临界荷位。

在板中、纵缝中部、角隅荷位下的横向弯拉应力分别如图13所示，每一单元大小为0.12 m×0.02 m，为了使云图更为细致地表现应力状态，竖向尺寸放大了 10倍。路面板上部表现为压应力，最大弯拉应力出现在面层层底。其中，对于板中和板边荷位，输出应力截面为z=2.5 m，角隅输出截面为z=4.9 m。类似地，纵向弯拉应力如图14所示。

综合考虑不同荷位下横纵向弯拉应力最大值，结果如表9所示。最大弯拉应力出现在板边荷位下的纵向拉应力，因此选择板边荷位作为临界荷载。对于角隅荷位，由于传力杆的约束作用，弯拉应力相较于其他荷位要小很多。

图13 不同荷位下路面板横向弯拉应力/MPaFig.13 Transverse bending stress of road panel under different loading positions

图14 不同荷位下路面板纵向弯拉应力/MPaFig.14 Longitudinal bending stress of road panel under different loading positions

表9 不同荷位下的最大弯拉应力Table 9 Maximum bending stress at different load levels

依据上述分析，对于横缝设置传力杆的水泥混凝土路面，临界荷位应选择板边荷位，这与文献[22]中水泥路面荷载应力的临界荷位分析结果一致。最大拉应力为左轮荷载下路面板底部纵向拉应力。

4 结论

本文介绍了基于半解析有限元方法的水泥混凝土路面力学响应计算理论，并根据相关理论和MATLAB开发了半解析有限元程序 CP-SAFEM，结合具体工况验证了 CP-SAFEM 具有良好的准确性与计算效率。

(1)针对单层板和双层板两种典型的水泥路面结构，按规范要求选择纵缝中部荷位，对比了规范、CP-SAFEM和ABAQUS的荷载应力计算结果，单层板和双层板下层的计算误差约为 1%，双层板上层计算误差约为6%。

(2)对于板中荷位下的水泥路面，CP-SAFEM和ABAQUS各项力学指标误差均在5%以内，进一步验证了CP-SAFEM的准确性。

(3)针对同一路面结构，采用相同的横截面网格划分形式，对比了CP-SAFEM和ABAQUS的模型复杂程度和计算时间，CP-SAFEM模型的单元数目仅为ABAQUS的1/50，运行时长仅为1/6，证明了CP-SAFEM良好的计算效率。

(4)通过CP-SAFEM对比了不同荷位下水泥混凝土路面的荷载应力，发现荷位为纵缝中部时出现最大荷载应力，证明了临界荷位应选择纵缝中部，这与现有研究结论相同。