移动荷载引发的道路结构力学响应一直是工程界所关心的课题。典型的道路结构具有横向有限，纵向无限的几何特点，同时地基的层状分布特性及各向异性属性也天然存在。对于柔性路面结构而言，在车辆移动荷载作用下，产生的振动波会通过近场道路传向远场地基。道路的结构参数会影响振动波的传播，而地基的结构特性也影响着道路结构的动力行为。因而，柔性路面与地基共同构成一种动力耦合体系，建立道路与地基的耦合结构模型是十分必要的。

柔性路面力学响应的研究常基于弹性(黏弹性)层状体系理论。马宪永等^[1]从横观各向同性黏弹性多层体系动力学方程出发，基于波传递方法研究了层间接触、材料及荷载对沥青柔性路面动力响应的影响。卢正等^[2]基于传递-反射矩阵法，推导了双圆轮胎振动荷载作用下弹性层状公路的动力响应解。LIU等^[3]建立了基于正交向量函数系的DVP方法，研究了层状柔性路面结构在表面荷载下的力学响应问题。整体层状体系的柔性道路力学分析模型能够很好地考虑振动波向外部地基的辐射效应，但是当需要考虑道路的实际几何特征时，解析类算法难以满足需要。

有限元法在模拟道路结构方面，有着绝对的优势，可以很方便地考虑道路结构的几何形态及材料的各向异性等。近来已有学者将有限元模型与层状地基求解理论相结合^[4-5]用来分析地震动入射波对于近场的影响，但其主要考虑二维平面结构问题，并不适用于分析移动荷载作用的道路结构问题。鉴于三维(3D)道路有限元模型需要耗费大量的计算成本，许多半解析有限元法得到广泛应用。沈凯仁等^[6-7]利用Fourier级数表征位移及荷载函数，建立了水泥路面的半解析有限元模型，并采用人工粘弹性边界来实现地基对振动波的吸收。LIU等^[8]同样借助Fourier级数展开建立了柔性路面的半解析有限元模型，并将无限元与其结合求解了柔性沥青路面在静力圆形荷载下的力学响应。BIAN等^[9]建立了车-轨道-地基系统的2.5D有限元-薄层单元模型，分析了列车荷载下高铁轨道的动力响应，但仅考虑底部为刚性基岩的情况。边界元法可以模拟波在结构边界向无穷远处传播的特性。SHENG 等^[10]和DEGRADEN等^[11]应用2.5D有限元法和2.5D边界元法求解了列车荷载诱发的大地振动问题。边界元法的应用需要基本解，当结构层数较多时，计算量会大幅增加。YASERI等^[12]应用2.5D有限元/比例边界有限元法研究了列车引起的地面振动，但只考虑远场为均质半空间的情况。YANG等^[13]应用2.5D有限元/无限元法^[14]研究了地下隧道在移动列车荷载下的动力响应问题。值得指出的是，无限元可以用来模拟波的辐射效应，但由于结构体内同时存在体波和面波，各类型波在体系内不同位置的贡献程度不同，每种波的位移振幅衰减因子也不相同，因此，需要人为地划分剪切波区、压缩波区和瑞利波区来近似模拟波混合的情况^[13]，含有一定的不确定性。此外，针对道路底部为层状地基的情况，2.5D有限元/无限元模型需要扩大近场道路的有限域模型尺寸，以适应无限元对层状地基的模拟。

综合来看，道路结构近场的半解析有限元模拟技术有着良好的应用前景，但有关近场边界条件的处理，即实现波向无穷远处辐射的模拟，仍存在着一些不足，或者不能考虑地基的半无限层状分布特性，或者存在地基层数增多所带来的求解不便。由此可见，道路-地基耦合模型的建立，关键在于实现对远场层状地基的高效模拟，以便于确定耦合界面处的边界条件。

有关层状地基的动力响应分析方法，主要包括传递矩阵法^[15-16]、刚度矩阵法^[17-18]、薄层法^[19-20]、解析层元法^[21]、传递-反射矩阵法^[2]、波函数法^{[1, 22]}、混合变量法^[23]和谱元法^[24]等。其中，传递矩阵法在处理高频或厚层问题时，会存在数值溢出的问题。刚度矩阵法和解析层元法的计算原理相近，基于解析类的方法可以获得层状体系动力响应的高精度解答。薄层法^[25]对层状介质进行垂直向的有限元离散，而在水平向采用解析的方法分析。薄层法在处理底部半无限域时，需要设置较厚的缓冲层或者吸收边界来近似模拟半空间体^[26]。同时薄层法的计算精度仍会受到薄层厚度的影响。传递-反射矩阵法及波传递法原理相近，通过构造层的上、下行波向量，可以消除传递矩阵中的正指数项，避免了数值溢出的问题^[1]。谱元法基于有限元思想，在变换域内将层状体系离散为一系列谱单元的组合，通过叠加有限个波数下的响应值即可获得时空域内的响应解答。谱元法所需的离散单元少，可以获得满意的求解精度^[27]。值得注意的是，解析类算法虽然可以获得问题的高精度解答，但在处理各向异性程度较高的层状介质时，需要重新推导求解，同时精确刚度矩阵的表达也会变得更加繁冗复杂^[28-30]。混合变量法是一种基于精细积分^[31]的层状地基动力响应求解方法，其求解精度高且具有无条件稳定性，基于矩阵代数的推导过程使得处理各向异性介质问题更加方便快捷。

本文基于谱元法思想和混合变量法原理，从各向异性介质的波动方程出发，借助位移谱向量和积分变换，推导出各向异性层状地基在任意荷载作用下的动力响应解答。然后依据子结构原理，建立起近场道路半解析有限元模型的边界条件，使得近、远场结构相耦合。最后通过快速Fourier变换(FFT)来获得柔性路面在时空域内的动力响应解答。该计算模型可以实现对道路近场几何形态的考虑，可以计入远场各向异性层状地基的影响，并且近场模型尺寸不会受到地基层数及层厚的影响。基于本文方法，分析了耦合柔性路面-层状地基模型与整体层状道路模型所得计算结果的差异，并分析了道路材料水平向及竖直向弹性模量对路面动力响应的影响。

1   计算模型与求解方法

图1为3D柔性路面-层状地基的耦合结构模型。将坐标原点O置于柔性道路表面，z轴竖直向下指向无穷远处。假定耦合体系的横截面几何形态及材料参数沿x轴保持不变，车辆移动荷载作用于路表，并沿x轴正向以速度c匀速移动。利用结构的特殊性，基于Fourier变换将问题转换到x向波数${\kappa _x}$域内，由此可只对波数${\kappa _x}$下形如体系横截面的2D耦合结构进行分析，于是3D道路问题的求解转化为多个独立的2D问题与解析问题的组合，而通过叠加不同波数${\kappa _x}$下的2D问题求解结果，即可获得原3D耦合结构的动力响应解答。

值得提及的是，道路的近场可以具有任意结构形式，路面可以高出原地面，也可以低于原地面，只要规划出待求解的近场部分即可。本文仅以图1所示简化的道路近场为例，阐述具体求解方法，其他结构形式的近场同理。

图  1  道路-层状地基耦合结构模型

Figure  1.  Road-layered foundation coupling model

下载: 全尺寸图片 幻灯片

1.1   近场道路的波数域有限元模型

建立近场道路的波数域有限元模型。如图2所示，原x向对应于波数域${\kappa _x}$。将波数${\kappa _x}$下的路面结构体离散成y-z平面内的若干个8结点曲边等参单元，单元内每个结点包含x、y及z三个方向的自由度，则波数$\kappa _{x}$下单元的位移表示为：

单元的坐标变换式表示为：

式中，${N_i}\left( {\xi ,\eta } \right)$为8结点曲边等参单元的形函数。

图  2  近场道路的波数域有限元模型

Figure  2.  Wavenumber finite element model of near-field road structure

下载: 全尺寸图片 幻灯片

假设移动车辆仅产生竖直向的简谐激振力，其振动幅值为${F_{\textit{z}}}$，激振频率为${\omega _0}$，将荷载的空间分布形式记为$\varPhi (x)\varPsi (y)$，则移动车辆荷载表示为：

将式(3)移动荷载写成波谱分量的形式：

式中：$\omega = {\omega _0} - {\kappa _x}c$；$\tilde \varPhi ({\kappa _x})$为$\varPhi (x)$的Fourier变换。在波数${\kappa _x}$下，移动荷载${\underline{\tilde f}{} _{\textit{z}}}$作用于有限元模型所产生等效结点载荷向量$\ {{{\underline{\tilde {{\boldsymbol{F}}}} }}} \ $可表示为：

式中，$\ {{{\tilde { {\boldsymbol{F}}}}}} \ $为结点力指示向量，用来指示荷载${F_{\textit{z}}}\varPsi (y)$产生的等效结点力及其所在结点自由度的位置。考虑简谐荷载$\ {{{\tilde { {\boldsymbol{F}}}}}}{{\text{e}}^{{\text{i}}{\kappa _x}x}}{{\text{e}}^{{\text{i}}\omega t}}$作用下的稳态响应问题，将结构位移写成波谱分量的形式：

式中：$\omega = {\omega _0} - {\kappa _x}c$；${\boldsymbol{\tilde u}} = {( {{{{\tilde{u}}}_x}}\;\;\;{{{\tilde u}_y}}\;\;\;{{{\tilde u}_{\textit{z}}}} )^{\text{T}}}$。根据有限元理论，将位移式(6)代入动力有限元方程可得：

式中：${\boldsymbol{\tilde S}}_{}^{{\text{fin}}}$为近场结构总动力刚度矩阵，${\boldsymbol{\tilde S}}_{}^{{\text{fin}}} = {\boldsymbol{\tilde K}}_{}^{{\text{fin}}} - {\omega ^2}{\boldsymbol{\tilde M}}_{}^{{\text{fin}}}$；${{\boldsymbol{\tilde U}}}$为近场结构整体结点位移波谱向量；${\boldsymbol{\tilde K}}_{}^{{\text{fin}}}$和${\boldsymbol{\tilde M}}_{}^{{\text{fin}}}$分别为近场结构总静力刚度矩阵和总质量矩阵，具体表示为：

式中：${[ {{\boldsymbol{\tilde K}}_{}^{{\text{fin}}}} ]^{{\text{ele}}}}$、${[ {{\boldsymbol{\tilde M}}_{}^{{\text{fin}}}} ]^{{\text{ele}}}}$分别为单元刚度矩阵和单元质量矩阵；$\displaystyle\sum\limits_{{\text{ele}}} {}$表示对各单元矩阵进行组装；$\rho$为单元材料密度；${\boldsymbol{D}}$为单元弹性矩阵，当需要考虑材料阻尼效应时，可将弹性矩阵${\boldsymbol{D}}$乘以系数$( {1 + 2{\text{i}}\zeta } )$，$\zeta$为材料阻尼比；$|{\boldsymbol{J}}|$为雅可比矩阵${\boldsymbol{J}}$的行列式，${\boldsymbol{J}}$可写为：

${{\tilde {\boldsymbol{B}}}} = [ {{{{{\tilde {\boldsymbol{B}}}}}_1}}\;\;{{{{{\tilde {\boldsymbol{B}}}}}_2}}\;\; \cdots \;\;{{{{{\tilde {\boldsymbol{B}}}}}_8}} ]$为波数$\kappa_{x}$下的单元应变矩阵，${{{\tilde {\boldsymbol{B}}}}_i}{{ (i = 1, 2, }} \cdots {\text{, 8)}}$可表示为：

式(7)的求解还需确定波数域有限元模型的边界条件。

1.2   远场各向异性层状地基的模拟

1.2.1   各向异性介质的波动方程

本节基于弹性波动理论，依据谱元法求解思想及精细积分技术，建立远场层状地基的动力响应计算模型。对于各向异性层状地基(图1)中的一般典型层来说，其满足波动方程(不计体力)：

式中：应力向量${{\underline{{\boldsymbol{\sigma}} } }} = ( {{{\underline{\sigma } }{}_{xx}}}\;\;{{{\underline{\sigma } }{}_{yy}}}\;\;{{{\underline{\sigma } }{}_{{\textit{z}}{\textit{z}}}}}\;\;{{{\underline{\tau } }{}_{xy}}}\;\;{{{\underline{\tau } }{}_{x{\textit{z}}}}}\;\;{{{\underline{\tau } }{}_{y{\textit{z}}}}} {)^{\text{T}}}$；位移向量${{ {\underline {{\boldsymbol{u}}}} }} = {( {{{ {\underline {u}} }{}_x}}\;\;\;{{{ {\underline {u}} }{}_y}}\;\;\;{{{ {\underline {u}} }{}_{\textit{z}}}} )^{\text{T}}}$；${{\underline{\ddot {\boldsymbol{u}}} }}$为位移对时间的二阶导数；${\boldsymbol{L}}$为微分算子。

弹性体的应变-位移关系表示为${\boldsymbol{\varepsilon }}{\text{ = }}{\boldsymbol{Lu}}$，应力-应变关系表示为${\boldsymbol{\sigma }}{\text{ = }}{\boldsymbol{D\varepsilon }}$，${\boldsymbol{D}}$为各向异性介质的弹性矩阵，为对称矩阵，当需要考虑地基的阻尼效应时，可用复弹性矩阵${{\boldsymbol{D}}^{\text{*}}}{\text{ = }}{\boldsymbol{D}}\left( {1 + 2{\text{i}}\zeta } \right)$来替代，其中$\zeta$为介质阻尼比。应用Fourier变换对可实现波动方程式(11)时间域与频域的相互转换：

式中，${\boldsymbol{\varGamma }}$为时空域内的任意变量。利用应力-应变及应变-位移关系，式(11)可转化成频域下位移表示的波动方程：

其中：

式中，${d_{ij}}{\text{ }}(i,j = 1,2, \cdots ,6)$为弹性矩阵${\boldsymbol{D}}$中第i行第j列的弹性常数。式(14)为包含各向异性介质弹性常数的二阶偏微分波动方程，其解为以e为底数的指函数。

1.2.2   谱单元刚度矩阵

1)双节点谱单元

使用双结点谱单元来模拟层状地基中非底部半空间的结构层。谱单元的上、下结点分别代表结构层的上、下两个界面。以第j号结构层为例，其对应于第j号双结点谱单元(如图3所示)，设该结构层的密度为${\rho _j}$，阻尼比为${\zeta _j}$，厚度${h_j}{\text{ = }}{{\textit{z}}_j} - {{\textit{z}}_{j - 1}} (j = 1,2, \cdots, l)$，弹性矩阵为${{\boldsymbol{D}}_j}$。将谱单元结点位移写成波谱分量的形式：

式中：${{\boldsymbol{\tilde u}}^{^j}} = {( {\tilde u_x^j}\;\;{\tilde u_y^j}\;\;{\tilde u_{\textit{z}}^j} )^{\text{T}}}$；${\kappa _x}$和${\kappa _y}$分别为沿x向和y向的波数。将式(16)代入式(14)可得：

其中：

式(17)为位移关于坐标z的二阶常微分方程。引入结点应力波谱向量${{\boldsymbol{\tilde p}}^j} = - {( {\tilde \tau _{x{\textit{z}}}^j}\;\;{\tilde \tau _{y{\textit{z}}}^j}\;\;{\tilde \sigma _{\textit{z}}^j} )^{\text{T}}}$，其与结点位移波谱向量${{\boldsymbol{\tilde u}}^j}$共同构成谱单元结点状态向量${{\boldsymbol{\tilde S}}^j} = {( {{{{\boldsymbol{\tilde u}}}^j}}\;\;{{{{\boldsymbol{\tilde p}}}^j}} )^{\text{T}}}$。根据应力-位移关系，应力向量满足关系式：

图  3  双结点谱单元

Figure  3.  2-noded spectral element

下载: 全尺寸图片 幻灯片

根据式(17)和式(19)，可推导出关于结点状态向量的一阶常微分方程：

其中：

常微分方程式(20)的解可写成：

由式(22)可建立谱单元的上、下结点的状态向量关系式：

式中，${{\boldsymbol{T}}^j}$为第j号谱单元的结点状态传递矩阵，其表达式为：

结点状态转移矩阵式(24)为矩阵指数函数，针对此类函数的计算，ZHONG等^[31]提出了一种基于2^N类算法的精细积分求解法，该方法可获得计算机意义下的高精度解，在合理的积分步长内，不会发生稳定性问题。本文应用精细积分法来求解谱单元内的结点状态转移矩阵。先将典型层在竖直向均分成厚度为$h_{{\text{mini}}}^j$的${2^N}$个细微薄层，则式(24)表示成：

一般N取20，于是$h_{{\text{mini}}}^j \approx {10^{ - 6}}h_{}^j$，$h_{{\text{mini}}}^j$一般远小于结构的振动波长，因此可利用Taylor展开将式(25)表示为：

令$({\boldsymbol{I}} + {\boldsymbol{T}}_n^j) = {({\boldsymbol{I}} +{\boldsymbol{T}}_{n - 1}^j)^2},{\text{ }}n = 1,2, \cdots ,N$，则可得${{\boldsymbol{T}}^j} = {\boldsymbol{I}} + {\boldsymbol{T}}_N^j$，其中${\boldsymbol{T}}_N^j$可通过递推公式${\boldsymbol{T}}_n^j = 2{\boldsymbol{T}}_{n - 1}^j + {({\boldsymbol{T}}_{n - 1}^j)^2},\; ({n = 1,2, \cdots ,N} )$来计算。式(26)的Taylor展开若略去${( {h_{{\text{mini}}}^j} )^5}$及以上的高阶项，则截断误差的量级为$O( {{{(h_{\min {\text{i}}}^j)}^4}} ) \approx {10^{ - 24}}O( {{{({h^j})}^4}} )$，已达到计算机的舍入误差范围，可认为是计算机意义下的精确解。

根据计算得到的${{\boldsymbol{T}}^j}$，将式(23)展开写成谱单元结点应力-位移关系式：

式中，${\boldsymbol{\tilde{u}}}、{\boldsymbol{\tilde{p}}}$的下标代表谱单元的结点编号。根据柯西应力法则和有限元理论，定义谱单元结点力向量和位移向量分别为：

则谱单元的结点力与结点位移关系可表示成：

其中：

将式(29)简写成：

式中，${{\boldsymbol{\tilde K}}^j}$即为双结点谱单元的刚度矩阵。

2) 单结点谱单元

层状体系底部的半空间只包含一个上界面，波在半空间内向下传播，不发生反射。采用单节点谱单元来模拟无限半空间，如图4所示。与双结点谱单元中的推导类似，单结点谱单元的结点状态向量满足关系式：

图  4  单结点半空间谱单元

Figure  4.  1-noded half-space spectral element

下载: 全尺寸图片 幻灯片

于是对${{\boldsymbol{H}}^{l + 1}}$求解特征值问题，定义${\boldsymbol{ \varPhi}} 、{\boldsymbol{\varLambda}}$分别为${{\boldsymbol{H}}^{l + 1}}$的特征向量和特征值矩阵，那么：

式(33)中的特征值成对出现，将${\boldsymbol{\varLambda }}$中的特征值按实部进行降序排列，并调整${\boldsymbol{\varPhi }}$中对应的特征向量位置，可得到：

为求解式(32)，引入向量${\boldsymbol{Q}} = {{\boldsymbol{\varPhi }}^{ - 1}}{{\boldsymbol{\tilde S}}^{l + 1}}$，其满足关系式：

式(35)为一阶常微分方程，其解可表示为：

式中，${{\boldsymbol{c}}_{g1}}$和${{\boldsymbol{c}}_{g2}}$为常数，则式(32)的解可表示为：

式中：${{\boldsymbol{\lambda }}_i}$为波从远处传播过来；$- {{\boldsymbol{\lambda }}_i}$为波向无穷远处传播出去。因此，式(37)中状态向量不应包含$\exp \left( {{\textit{z}}{{\boldsymbol{\lambda }}_i}} \right)$项，则${{\boldsymbol{c}}_{g1}} = 0$。于是可得到单结点谱单元结点应力-位移关系：

将式(38)简写为：

式中，${{\boldsymbol{R}}^\infty }$即为单结点谱单元的刚度矩阵。

1.2.3   层状体系整体谱刚度矩阵

根据传统有限元理论，图1中的层状地基整体谱刚度矩阵应由l个双结点谱单元和一个单结点谱单元组装而成。假定层间接触条件为完全连续，那么相邻两结构层的谱单元在公共结点处的位移相同，则有：

式中，${\boldsymbol{\tilde u}}_j^{}$为第j号结点的位移向量。将每个谱单元的刚度矩阵按照整体结点编号顺序“对号入座”，进而建立层状地基的整体平衡方程：

式(41)可简写成：

式中：$[{\boldsymbol{\tilde K}}(\omega )]$为层状地基的整体谱刚度矩阵；上标mn对应于当前波数${\kappa _{xm}}$和${\kappa _{yn}}$；${{\tilde {\boldsymbol{F}}}}_j^{}{\text{ }}(j = 1,2, \cdots ,l)$为施加于谱单元结点j处的结点载荷向量。根据结点平衡关系，结点外载荷${{\tilde {\boldsymbol{F}}}}_j^{}$应与环绕结点j的各谱单元在结点j处的结点力之和相等，即：

通常结点外载荷向量为已知量，则结点位移可通过下式求解：

式中，$[{\boldsymbol{\tilde G}}(\omega )] = {[{\boldsymbol{\tilde K}}(\omega )]^{ - 1}}$为层状地基的整体谱柔度矩阵。为与波数域有限元法相适应，根据谱元法理论，波数${\kappa _{xm}}$下层状地基位移响应${\boldsymbol{\tilde U}}({\kappa _{xm}},y,t)$可通过叠加有限个波数${\kappa _{yn}}$下的波谱解来获得，于是：

式中：$\omega = {\omega _0} - {\kappa _{xm}}c$；m和n的取值范围分别为$\left[ { - M,M} \right]$和$\left[ { - N,N} \right]$，N的取值需要足够大以涵盖激振产生的大多数波形。在计算式(45)时，本文采用快速Fourier变换(FFT)来加快求解速度。

1.3   近-远场结构的耦合

1.3.1   边界条件的确定

道路近场波数域有限元模型通过边界结点与周围地基连接，假定连接点处的位移满足连续条件，则近场模型的边界条件，可通过建立周围地基在边界结点处的动力平衡关系而获得，即边界结点动刚度的求解。本文利用容积法^[25]，并结合子结构理论来确定周围地基的边界结点动力刚度。如图5所示，不规则地基的边界结点动刚度，可等效为规则原场地基的结点动刚度与开挖土体的动刚度之差。具体步骤为：

图  5  道路-地基结构分解示意图

Figure  5.  Schematic diagram of road-foundation structure decomposition

下载: 全尺寸图片 幻灯片

1)确定开挖土体的动刚度${{\boldsymbol{\tilde S}}^{\text{e}}}({\kappa _x})$。建立开挖土体的波数域有限元模型，方法同1.1节。在对开挖土体进行有限域离散时，其边界结点的位置应与近场道路的波数域有限元模型边界结点位置相同，这样就可以保证开挖地基边界条件对近场道路模型的适用性。在波数${\kappa _x}$下，设开挖土体的刚度矩阵和质量矩阵分别为${{\boldsymbol{\tilde K}}^{\text{e}}}({\kappa _x})$和${{\boldsymbol{\tilde M}}^{\text{e}}}({\kappa _x})$，则开挖土体的动刚度矩阵可表示为${{\boldsymbol{\tilde S}}^{\text{e}}}({\kappa _x}) = {{\boldsymbol{\tilde K}}^{\text{e}}}({\kappa _x}) - {\omega ^2}{{\boldsymbol{\tilde M}}^{\text{e}}}({\kappa _x})$。

2)确定规则原场地基的动刚度${\boldsymbol{\tilde S}}({\kappa _x})$。在原场地基中设置与开挖土体相同的离散结点群，在第i个结点处施加合外力为${{\boldsymbol{p}}_i}(x,t) = {( {{p_{ix}}}\;\;{{p_{iy}}}\;\;{{p_{iz}}} )^{\text{T}}}$的结点外荷载，${{\boldsymbol{p}}_i}(x,t)$均匀地作用在结点i处宽度为$2{b_i} = \left| {y_i^1 - y_i^2} \right|$的区域内，其中${y_i}$表示i结点的横向坐标，$y_i^1$和$y_i^2$为作用区域的两侧界限值，一般取为相邻两结点的中点，于是结点i处作用于回填地基的荷载可表示为：

对式(46)作Fourier变换，得到变换域内的结点荷载：

将式(47)代入式(45)可得到i结点荷载引起的j结点的位移响应：

那么在波数${\kappa _x}$下，规则原场地基结点群荷载引起的结点群位移响应可表示为：

将式(49)简写成${\boldsymbol{\tilde U}}({\kappa _x}) = {\boldsymbol{\tilde Q}}({\kappa _x}){\boldsymbol{\tilde P}}({\kappa _x})$，${\boldsymbol{\tilde Q}}({\kappa _x})$即为波数${\kappa _x}$下规则原场地基离散结点群的动力柔度矩阵，则规则原场地基结点群的动刚度矩阵可表示为${\boldsymbol{\tilde S}}({\kappa _x}) = {[{\boldsymbol{\tilde Q}}({\kappa _x})]^{ - 1}}$。

3)确定开挖地基边界结点群动刚度${{\boldsymbol{\tilde S}}}_{\text{b}}^{\text{g}}({\kappa _x})$。基于子结构原理，首先对规则原场地基及开挖土体的结点群自由度进行划分，其中开挖边界(即开挖土体与开挖地基的交界面)对应位置的结点自由度用下标b来表示，其余自由度用下标i来表示；然后对原场地基及开挖土的结点群动刚度${{\tilde {\boldsymbol{S}}}}$和${\tilde {\boldsymbol{S}}^{\text{e}}}$按照自由度类别b和i分别进行自由度重排，并进行矩阵分块处理；最后利用子结构的凝聚原理^[32]，将原场地基及开挖土体的结点群动刚度向交界面结点群凝聚。于是，可得凝聚后的原场地基交界面结点群动刚度矩阵为${\boldsymbol{\tilde S}}_{\text{b}}^{} = {\boldsymbol{\tilde S}}_{{\text{bb}}}^{} - {\boldsymbol{\tilde S}}_{{\text{bi}}}^{}{( {{\boldsymbol{\tilde S}}_{{\text{ii}}}^{}} )^{ - 1}}{\boldsymbol{\tilde S}}_{{\text{ib}}}^{}$和凝聚后的开挖土体交界面结点群动刚度矩阵为${\boldsymbol{\tilde S}}_{\text{b}}^{\text{e}} = {\boldsymbol{\tilde S}}_{{\text{bb}}}^{\text{e}} - {\boldsymbol{\tilde S}}_{{\text{bi}}}^{\text{e}} {( {{\boldsymbol{\tilde S}}_{{\text{ii}}}^{\text{e}}} )^{ - 1}}{\boldsymbol{\tilde S}}_{{\text{ib}}}^{\text{e}}$，那么开挖地基边界结点群动刚度可表示为：

式(50)即为近场道路波数域有限元模型的边界条件。

1.3.2   近-远场模型耦合及动力响应求解

对近场道路波数域有限元动刚度矩阵${\boldsymbol{\tilde S}}_{}^{{\text{fin}}}$(式(7))按照自由度类别b和i进行自由度重排，之后按照结点自由度对应关系，将开挖地基边界结点的动刚度矩阵${\boldsymbol{\tilde S}}_{\text{b}}^{\text{g}}$叠加到近场道路模型的动刚度矩阵${\boldsymbol{\tilde S}}_{}^{{\text{fin}}}$中，可得道路-层状地基耦合结构整体动力平衡方程：

将式(51)简写为$[ {{\boldsymbol{\tilde S}}_{\text{a}}^{{\text{fin}}}({\kappa _x})} ]\{ {{\boldsymbol{\tilde U}}({\kappa _x})} \} = \{ {{{\tilde { \boldsymbol F}}}} \}$，于是道路模型结点位移可表示为$\{ {{\boldsymbol{\tilde U}}( {{\kappa _x}} )} \} = {[ {{\boldsymbol{\tilde S}}_{\text{a}}^{{\text{fin}}}({\kappa _x})} ]^{ - 1}}\{ {{{\tilde { \boldsymbol F}}}} \}$。在式(3)所表示的移动荷载下，道路耦合结构在时空域内的结点位移响应可通过叠加不同波数${\kappa _x}$下的波谱响应值得出，即：

式中，$\omega = {\omega _0} - {\kappa _{xm}}c$。本文采用快速Fourier变换(FFT)来计算式(52)，其中m的取值范围为$\left[ { - M,M} \right]$，M取值需足够大以涵盖激振产生的大多数波形。

2   方法验证及算例分析

2.1   计算方法验证

为验证本文方法的正确性，将本文方法与精确类解法进行计算结果的对比验证。

2.1.1   均质各向同性体的对比验证

将本文方法与HUNG和YANG^[33]的解析方法进行计算结果的对比验证。考虑各向同性均质半空间体^[33]，其材料参数为：弹性模量$E = 50{\text{ MPa}}$，泊松比$\nu = 0.25$，阻尼比$\zeta = 0.02$，密度$\rho = 2000{\text{ kg/}}{{\text{m}}^3}$。外荷载沿x轴正向移动，其波数${\kappa _x}$域内表达式为${\tilde f_{1{\textit{z}}}} = 4T/(4 + \kappa _x^4{\alpha ^4})$，其中轮载$T = {\text{10 t}}$，特征长度$\alpha = 0.8{\text{ m}}$，荷载自振频率为10 Hz。为与解析解相比较，采用本文方法计算时，建立了宽度为20 m，厚度为1 m的近场道路波数域有限元模型，模型采用二维8结点有限单元来离散，形成了上、下两层共80个单元和325个结点，网格结点间距为0.25 m。坐标原点O位于模型上表面中心位置处(参考图2)。远场地基与近场道路的材料参数保持一致。外荷载${\tilde f_{\textit{z}}}$仅是关于纵轴x向波数${\kappa _x}$的函数，在横轴y向可视为集中荷载，因此将其直接施加在原点O上。FFT求解参数为：采样频率$- 32 < {\kappa _x} < 32$，采样点数范围值$M，N=2048$。选取$t = 0{\text{ s}}$时刻，移动载荷通过坐标原点时，竖直向位移${u_{\textit{z}}}$沿轴(y=0 m, z=1 m)的分布作为观察对象，对比两种荷载速率c=20 m/s和c=50 m/s下，本文结果与解析结果的差异。由图6可知，文本解与解析解在两种移动速度下都具有较高的吻合度，因此，本文方法是行之有效的。

图  6  本文方法与HUNG方法^[33]的计算结果对比

Figure  6.  Comparison between the results for the proposed method and those in HUNG^[33]

下载: 全尺寸图片 幻灯片

2.1.2   层状横观各向同性体的验证

对于层状横观各向同性介质的情况，将本文方法与混合变量法(MVF)^[23]进行计算结果的对比验证。考虑文献[28]中的5层状横观各向同性体，各层材料参数详见文献[28]中表6内的工况2。外荷载为沿x轴正向以速度$c = 70{\text{ km/h}}$移动的点荷载。其波数${\kappa _x}$域内的表达式为${\tilde f_{2{\textit{z}}}} = 1$。近场道路波数域有限域模型与2.1.1节保持一致，远场地基材料参数与近场保持一致。外荷载${\tilde f_{2{\textit{z}}}}$仅是常函数点荷载，因此将其直接施加在近场模型的原点o上。FFT求解参数为：采样频率$- 32 < {\kappa _x} < 32$，采样点数范围值$M，N=1024$。选取$t = 0{\text{ s}}$时刻，移动载荷通过坐标原点时，竖直向位移${u_\textit{z}}$沿轴(y=0 m, z=0 m)和轴(y=0 m, z=1 m)的分布作为观察对象，对比本文结果与MVF计算结果的差异。由图7可知，本文方法与MVF法所求解的路表面及路基底面下的动位移，具有较高的吻合度。综上验证算例可以表明，无论对于均质体还是层状体，各向同性体还是各向异性体，本文方法均是行之有效的。

图  7  本文方法与混合变量法^[23]的计算结果对比

Figure  7.  Comparison between the results for the proposed method and the mixed variable formulation^[23] method

下载: 全尺寸图片 幻灯片

2.2   数值算例分析

2.2.1   耦合模型与整体层状模型的计算结果对比

道路-层状地基耦合模型与以往整体层状结构模型在结构形式上存在一定的差异，图8为两种简化模型的示意图。本节将对比分析这两种模型所得计算结果的差异。道路结构属于多层结构，面层和基层材料力学性能差异显著，为便于分析两种简化模型所得计算结果的差异这一问题，以及进一步减小近场模型的计算规模，根据《公路路线设计规范》^[34]和文献[2]中的道路试验，对两种模型的计算参数进行适当的配置。对于道路-层状地基耦合模型，建立宽度为20 m，厚度为1 m的近场道路模型。不失一般性，将多层道路结构^[2]简化为两个主要结构层：上层为包含面层、基层及底基层的路面结构层部分，处于标高0 m~0.5 m范围内；而下层为包含各路基层的路基结构层部分，处于标高0.5 m~1 m范围内。路面结构层的材料参数设为${E_1} = 1000{\text{ MPa}}$，${\nu _1} = 0.3$，${\zeta _1} = 0.02$，${\rho _1} = 2200{\text{ kg/}}{{\text{m}}^3}$；路基结构层的材料参数设为${E_2} = 60{\text{ MPa}}$，${\nu _2} = 0.35$，${\zeta _2} = 0.02$，${\rho _2} = 2000{\text{ kg/}}{{\text{m}}^3}$，地基参数和其他计算参数与2.1.1节中相同。对于道路整体层状结构模型，层状体系顶部两层分别设置为路面层和路基层，厚度均为0.5 m，相关计算参数与耦合模型的保持一致。

图  8  道路耦合结构模型与整体层状结构模型示意图

Figure  8.  Schematic diagram of coupling model and overall layered model for road structure

下载: 全尺寸图片 幻灯片

计算结果如图9所示，可以看到，无论在路表(z=0 m)还是在路基底部(z=1 m)，耦合模型所得的振动幅值要比整体层状模型的大，并且前者的波峰间距也比后者的小。由此可见，道路耦合模型下的波动情况更加剧烈，若单纯采用整体层状模型来求解道路动力响应的话，会低估道路的实际振动水平。

图  9  道路耦合模型与整体层状模型的竖向位移

Figure  9.  Vertical displacement of coupling model and overall layered model for road structure

下载: 全尺寸图片 幻灯片

2.2.2   道路结构层弹性模量的影响

本节基于道路-层状地基耦合模型，研究横观各向同性路面材料的弹性模量对路面动力响应的影响。横观各向同性材料包含5种独立的弹性常数，分为各向同性平面内的水平向弹性模量${E_{\rm{h}}}$，各向同性平面内的泊松比${\nu _{\rm{h}}}$，垂直于各向同性平面的竖直向弹性模量${E_{\rm{v}}}$，垂直于各向同性平面的剪切模量${G_{\rm{v}}}$，垂直于各向同性平面的泊松比${\nu _{\rm{v}}}$，而各向同性平面内的剪切模量${G_{\rm{h}}} = {E_{\rm{h}}}/2(1 + {\nu _{\rm{h}}})$^[1]。本节在2.2.1节道路耦合结构算例的基础上，对模型尺寸、网格、阻尼比、密度及泊松比等参数不变，而对路面层及路基层的材料弹性模量进行设定(如表1)，以研究道路材料模量${E_{\rm{h}}}$和${E_{\rm{v}}}$的影响规律。地基材料参数见表2。表1中：a为弹性模量${E_{\rm{h}}}$与${E_{\rm{v}}}$之间的比值，本文a分别取0.5、1.0、1.5、2.0四组对比值^[35]。工况1-1与工况1-2用来分析路面层内${E_{\rm{h}}}$和${E_{\rm{v}}}$对路面动力响应的影响，工况2-1与工况2-2用来分析路基层内${E_{\rm{h}}}$和${E_{\rm{v}}}$对路面动力响应的影响，同时对比路基层与路面层的影响效果。

表  1  不同工况下道路材料的弹性模量

Table  1.  Elastic modulus of road materials for different cases

道路层	工况1-1 模量/MPa	工况1-2 模量/MPa	工况2-1 模量/MPa	工况2-2 模量/MPa
路面层	1000–$a{E_{\rm{h}}}$–600	$a{E_{\rm{v}}}$–1000–600	2000–1000–600	2000–1000–600
路基层	50–60–55	50–60–55	50–$a{E_{\rm{h}}}$–55	$a{E_{\rm{v}}}$-50-55
注：模量的格式为${E_{\rm{h}}} - {E_{\rm{v}}} - {G_{\rm{v}}}$，$a$为弹性模量${E_{\rm{h}}}$与${E_{\rm{v}}}$之间的比值，分别取0.5、1.0、1.5、2.0四组对比值。

下载: 导出CSV 

| 显示表格

表  2  地基材料参数

Table  2.  Material parameters of ground

地基层	模量/MPa	泊松比${\nu _{\rm{h}}}{\text{ = } }{\nu _{\rm{v}}}$	密度$\rho /({\text{kg}}\cdot{{\text{m}}^{-3}})$	层厚/$\text{m}$	阻尼比$\zeta$
地基1层	30−37.5−15	0.30	2000	2	0.02
地基2层	32−39−16	0.30	2200	4	0.02
地基3层	36−40−18	0.25	2200	−	0.02
注：模量的格式为${{E}_{{\rm{h}}}}-{{E}_{{\rm{v}}}}-{{G}_{{\rm{v}}}}$。

下载: 导出CSV 

| 显示表格

由图10知，路面层${E_{\rm{v}}}$的变化对路面动位移的影响总体不太明显，但仍存在一定的规律：随着路面层${E_{\rm{v}}}$的增大，最大位移幅值有所降低，而其他波峰幅值有所升高，且波峰间距减小。路面层${E_{\rm{h}}}$对动位移的影响较为明显，随着${E_{\rm{h}}}$的增大，路面动位移幅值明显降低，移动荷载前后的波峰间距明显增大，因而，路面层的${E_{\rm{h}}}$相较于${E_{\rm{v}}}$对路面动力响应的减弱作用更加明显。由图11可知，路基层的${E_{\rm{v}}}$越大，路表最大位移幅值越低，但降低率越来越小，而其他波峰处的幅值随着路基层${E_{\rm{v}}}$的增大而越来越大，同时波峰间距减小。路基层${E_{\rm{h}}}$对路面的动力响应的影响较为明显，较大的${E_{\rm{h}}}$值会引起路面动位移幅值显著增大。对比图10和图11可以看到，路基层弹性模量对路面动位移的影响相较于路面层要更加明显。综上分析，路基层模量${E_{\rm{h}}}$和${E_{\rm{v}}}$的升高都会增强路面的波动水平，而${E_{\rm{h}}}$的增强作用更加明显。路面层${E_{\rm{h}}}$的增大可以较大程度地降低路面的动力响应。

下载: 全尺寸图片 幻灯片

  10  工况1-1和工况1-2不同模量比率a下位移沿x轴和y轴的分布

  10.  Distribution of displacement along x-axis and y-axis under different modulus ratio a for case 1-1 and case 1-2

下载: 全尺寸图片 幻灯片

图  11  工况2-1和工况2-2不同模量比率a下位移沿x轴和y轴的分布

Figure  11.  Distribution of displacement along x-axis and y-axis under different modulus ratio a for case 2-1 and case 2-2

下载: 全尺寸图片 幻灯片

3   结论

本文介绍了基于半解析法的移动荷载作用下柔性路面-各向异性层状地基耦合结构动力响应的计算原理。根据这一原理进行了数值算例验证，同时对比分析了道路-地基耦合模型与整体层状道路模型所得动力响应的差异，并结合具体工况分析了道路结构层弹性模量的影响规律。

(1)数值算例表明：本文柔性路面-层状地基耦合结构的半解析解与解析结果具有较高的吻合度，可为工程路面结构设计提供可靠的计算依据。

(2)本文提出的基于有限元及谱元法的半解析方法具有广泛的适用性，可以考虑移动荷载作用的情况，可以考虑道路近场的复杂几何形态，还可以考虑地基的层状分布及各向异性特性等，并且近场有限域尺寸不会受远场地基层数及层厚度的影响，数值计算稳定，求解精度较高。

(3)道路-层状地基耦合模型与以往道路整体层状结构模型相比，前者求解的路面动力响应更为剧烈。对于横观各向同性道路而言，路基层弹性模量${E_{\rm{h}}}$或${E_{\rm{v}}}$的升高会增强路面的波动水平，其中${E_{\rm{h}}}$的增强作用更加明显。路面层${E_{\rm{h}}}$的升高可以显著地降低路面的动力响应，而${E_{\rm{v}}}$的升高则会提高路面整体振动水平。路基层弹性模量对路面振动的影响相较于路面层更加明显。