轨道交通运行过程中，列车、轨道结构和场地土层三个子系统之间始终存在着复杂的动态变化的相互作用. 针对三者的耦合振动，一些学者(如 Sheng et al.，1999a，1999b，2003，2004; Lombaert et al.，2006，2009)建立了精妙的分析模型.其中场地对简谐荷载的动力响应函数，即频域格林函数，是解耦列车-轨道-场地系统的必要的关键因素.模型中的场地土被抽象为单相的黏弹性连续介质，没有考虑孔隙水与土骨架的流固两相耦合振动问题.天然土的孔隙水可能会分担土中的应力，也可能相对于土骨架运动而产生内摩擦.很多情况下，土中水对场地振动的影响相当大，基于单相介质模型的方法难以给出满意 的结果(Cai et al.，2009; Beskou and Theodorakopoulos，2011).

Biot(1956a，1956b，1962)建立了流固两相介质运动方程，被广泛应用于饱和土场地对地表移动荷载的动力响应问题.其中很多研究(Jin et al.，2004; Lu and Jeng，2007; Lefeuve-Mesgouez and Mesgouez，2008; Sun et al.，2010; Cao et al.，2011)关注均匀半空间两相介质的波动，忽略了场地土的分层现象.均匀半空间模型易于获得闭合形式的解析解，有助于理解饱和土场地振动的一些基本特性.但是由于没有分层界面，均匀半空间不能反映土层界面上反射、折射和透射等一些典型的波动现象.

比较而言，两相介质的层状半空间更接近场地的物理实际.针对移动荷载激励问题，一些学者(Xu et al.，2007，2008; Lu et al.，2009; Mesgouez and Lefeuve-Mesgouez，2009; Lefeuve-Mesgouez and Mesgouez，2012)采用透射反射矩阵方法(Lowe，1995; Luco and Apsel，1983; Rokhlin and Wang，2002)计算了场地动力响应.少数的研究采用了2.5维有限元(Gao et al.，2012)和薄层法(高广运等，2013)等数值模型.

本文可以看做是Sheng et al.(2003，2004)模型的进一步拓展，目的是构建饱和层状半空间的竖向位移格林函数，用以取代该文献中的单相介质格林函数，使其能够模拟列车引起的饱和场地振动.论文首先给出频率-波数域内Biot方程的三维解答，基于传递矩阵方法推导格林函数矩阵，然后探讨改善格林函数数值计算效率的措施，最后通过具体算例验证格林函数的可靠性并展示饱和层状半空间的振动特征.

在流体饱和的多孔弹性介质中建立笛卡尔坐标系，用U_i(i=x，y，z)表示固体骨架位移张量，用w_i表示流体相对于固体骨架的位移张量，依据Biot(1962)理论，在忽略孔隙形状影响的条件下，流固两相的运动方程为

在频率为ω的谐振条件下，对(1)式两边点乘笛卡尔空间梯度算子，得到体变方程：

对式(1)应用式(4a)定义的Fourier变换，将式(5)的体应变解代入，得到位移解，

如图 1所示，层状半空间由n个水平覆盖层和一个下卧均匀半空间构成.自上而下依次编号，第n+1“层”代表下卧半空间.层状半空间的材料均为饱和多孔介质，不同层物理力学参数有所不同，层的内部则处处相同.每层具有自己的局部坐标系，原点位于各自的顶面.第1层的局部坐标系也为全局坐标系，原点位于层状半空间顶面.设第j层的厚度为h_j，则在局部坐标系中，z_j=0表示顶面，z_j=h_j表示底面.分别采用${{S}_{j,0}}={{S}_{j}}\left( {{z}_{j}}=0 \right),{{\bar{\tilde{u}}}_{j,0}}={{\bar{\tilde{u}}}_{j}}\left( {{z}_{j}}=0 \right),{{\bar{\tilde{\sigma }}}_{j,0}}={{\bar{\tilde{\sigma }}}_{j}}\left( {{z}_{j}}=0 \right)$的记号方式表示顶面的状态向量、位移向量和应力向量，采用 ${{S}_{j\text{,}1}}={{S}_{j}}\left( {{z}_{j}}={{h}_{j}} \right)\text{,}{{\bar{\tilde{u}}}_{j\text{,}1}}={{\bar{\tilde{u}}}_{j}}\left( {{z}_{j}}={{h}_{j}} \right)\text{,}{{\bar{\tilde{\sigma }}}_{j\text{,}1}}={{\bar{\tilde{\sigma }}}_{j}}\left( {{z}_{j}}={{h}_{j}} \right)$，表示底面处的各向量.根据式(8)，

图 1

Fig. 1

图 1 层状半空间场地模型 Fig. 1 Layered half-space model for ground

消去 b _j得

在下卧半空间中，当z趋于无穷远时，位移必趋近于0.据此，可以建立下卧半空间顶面(也就是第n层底面)处的位移-应力关系:

将式(13)代入式(14)得

一般轨道周边地面运动的竖向分量远大于水平向分量，通常建模中仅考虑竖向激励而忽略两个水平向的剪切荷载(Sheng et al.，1999; Sheng et al.，2003; Sheng et al.，2004; Lombaert et al.，2006; Lombaert and Degrande，2009).在地面可渗水的条件下，式(17)的竖向位移分量

根据物理意义可知，竖向位移格林函数${{\bar{\tilde{\delta }}}_{44}}$在(β，γ)坐标平面上具有轴对称性.这个属性能够有效减少数值计算量.只须沿着某方向的极轴r计算${{\bar{\tilde{\delta }}}_{44}}$，则整个(β，γ)坐标平面上的值都可获得.不失一般性，可以将r的方向选为γ轴的正向.换言之，如果仅考虑竖向激励和竖向响应，本文附录矩阵的计算过程中可将β设为0.这将使矩阵中的一些元素成为零元素.令 K = A _u A_τ ^－1 T ₂₁－ T ₁₁，令 Q = T ₁₂－ A _u A_τ ^－1 T ₂₂，则式(15)可写为

Sheng等(2003，2004)建立了列车-轨道-场地耦合动力分析模型，全系统如图 2所示.列车子系统用弹簧阻尼器连接的多刚体模拟，质量矩阵为 M _T，复阻尼刚度矩阵 K _T.无限长轨道与场地接触面宽度为2b_t.钢轨用Euler梁模拟，单位长度质量为m_R，截面抗弯刚度EI.轨枕连续化成为沿轨道方向延伸的无抗弯刚度扁梁，单位长度质量m_s.轨垫板用竖向弹簧代表，刚度为k_p.道碴路堤用线密度m_B的阻尼弹性层代表，用一致质量近似方法简化，仅考虑竖向刚度k_B.EI、k_B和k_p均为复数，虚部包含各自的耗散系数.假定轨道与场地接触面上只存在竖向相互作用力，并且沿着垂直轨道的水平方向均匀分布.

图 2

Fig. 2

图 2 列车-轨道-场地耦合振动体系 Fig. 2 The train-track-ground coupling system

采用本文推导的两相介质格林函数${{\bar{\tilde{\delta }}}_{44}}$替代原模型的单相介质格林函数，即可模拟列车运行引起的饱和层状场地波动.限于篇幅，本文不再赘述模 型的计算流程.力学建构和公式推导详见文献Sheng等(2003，2004).

首先计算层状半空间上直接作用矩形分布荷载的简单算例，然后模拟包括列车、轨道和饱和场地全系统的工况.

Jones等(1998)和Lefeuve-Mesgouez and Mesgouez(2008)分别计算了矩形分布简谐荷载对单相介质层状半空间和两相介质均匀半空间的激励.前者的荷载作用位置不变，即非移动激励; 后者的荷载在地面匀速直线运动，为移动激励.现采用前面推导的格林函数分别计算这两个简单算例，与文献结果对比.

如图 3所示，竖向简谐荷载均匀分布于地表矩形区域，坐标原点设在矩形的中心.频率-波数域内荷载可表示为

图 3

Fig. 3

图 3 四个场地模型及其矩形分布简谐荷载 Fig. 3 Four ground models with rectangular harmonic load acting on surfaces

将式(23)代入式(18)得到(β，γ，ω)域中的地表竖向位移表达式：

Jones等(1998)采用单相介质模型，利用表 1中的两套参数区分A、B两种土.A、B土的不同组合形成4个场地模型，如图 3所示.场地1为A土构成的均匀半空间; 场地2为7 m厚的A土覆盖层加B土下卧半空间; 场地3为7 m厚A土覆盖层加下卧刚性基岩; 场地4为B土构成的均匀半空间.Biot(1962)指出，当参数ρ_f，M和b趋近于0时，固相的动力特性将趋近于相同参数值的单相介质.为验证本文推导的格林函数，固相参数值设为表 1 单相介质的相应值; 其余参数ρ_f=0.0001 kg·m^-3，M=0.0001 Pa，b=0，α=1.0和φ=0.6.

表 1 Table 1

表 1 单相介质参数值(Jones et al.(1998)) Table 1 Parameters of single phase media in Jones et al.(1998) 参数 单位 A土 B土 剪切模量μ0 MPa 107 428 拉梅常数λ0 MPa 113 453 密度ρ kg·m-3 1550 2000 阻尼比η 0.05 0.05

表 1 单相介质参数值(Jones et al.(1998)) Table 1 Parameters of single phase media in Jones et al.(1998)

图 4为计算得到的地表位移波数幅值谱${{\bar{\tilde{U}}}_{z}}\left( \beta \text{,}\gamma \right)$和空间幅值 ${{\tilde U}_z}\left({x，y} \right)$，可见前述采样点数和采样间隔足够，有效避免了混叠和泄漏.图 5为地表振动的绝对幅值|${{\tilde U}_z}$(x=0)| 沿y轴的衰减关系，其中实线为Jones等(1998)单相介质模型的计算结果，*号代表本文方法得到的固相位移.可见本文结果与文献结果良好符合，初步验证了格林函数的可靠性.

图 4

Fig. 4

图 4 当ρf，M和b取值极小时场地1的竖向位移 (a)波数域位移谱;(b)空间域位移幅值. Fig. 4 Vertical displacement of the Ground 1，when values of ρf，M and b are very small (a)Wave number domain;(b)Space domain.

图 5

Fig. 5

图 5 x=0处的竖向位移绝对值 Fig. 5 Amplitudes of vertical displacements along x=0

前述矩形分布简谐荷载沿x轴以速度c移动，式(24)变为

依据Lefeuve-Mesgouez and Mesgouez(2008)，场地土骨架无阻尼拉梅常数μ₀=462 MPa，λ₀=692 MPa，阻尼比η=0.05，土体密度ρ=1960 kg·m^-3，Biot第一参数M=5100 MPa，第二参数α=0.975，流体密度ρ_f=1000 kg·m³，孔隙率φ=0.4，黏滞耦合系数b=1×10⁹(Pa·s·m²).采用本文推导的格林函数依据式(25)计算地表位移，做关于(β，γ)的逆Fourier变换，得到(x－ct，y)移动空间域的幅值.图 6 展示了y=0处移动坐标系中地表竖向位移的瞬 间形态.对比图 4b可见，由于荷载移动，位移幅值不再关于荷载中心对称.本例中荷载移动速度122.5 m·s^-1，约为半空间瑞利波速的1/4.如果荷载速度进一步增加，非对称性将会更为明显(Lefeuve-Mesgouez and Mesgouez，2008).图 7 展示了移动坐标轴上的振动绝对幅值|${{\tilde U}_z}$(y=0)|，可见本文结果与文献基本符合.

图 6

Fig. 6

图 6 在随矩形分布荷载中心点一起移动的坐标系中，y=0位置处的地面竖向位移(荷载速度为122.5 m·s-1) Fig. 6 Vertical displacement of ground surface at y=0 in the moving frame of reference bound with the center of a harmonic square load traveling at the speed of 122.5 m·s-1

图 7

Fig. 7

图 7 移动坐标系中，地面y=0处的竖向位移绝对值 Fig. 7 Modulus of the vertical displacement of ground surface at y=0 in the moving frame of reference

笔者所在课题组曾对北京城轨13号线地面运行区间进行过现场实测(王福彤等，2011b)，地表测点的平面布局如图 8所示.列车4节编组，行驶速度 为60 km·h^-1(或17 m·s^-1).利用强震仪以200 Hz 的采样频率记录地表测点的竖向加速度; 根据加速度时程进行功率谱估计，并去除了本底振动成分(王福彤等，2011a; 郑鑫等，2013); 按照国家标准《城市区域环境振动量测方法》(GB10071-88)将功率谱转换为振动加速度级VAL谱，如图 9所示.

图 8

Fig. 8

图 8 地表观测点平面布置图 Fig. 8 Plane layout of the observation points on the ground surface

图 9

Fig. 9

图 9 地表测点上振动加速度级的模型计算值与观测值的对比 Fig. 9 Comparison between the simulated values and the observed values of the vibration acceleration levels (VALs) at the observation points on the ground surface

采用第5节模型模拟计算地表振动.表 2、表 3分别为车辆和轨道结构的模型参数表，数值基于设计资料综合确定.表 4数据为场地土介质参数，源于钻探取样和原位波速测结果.现场地下水位埋藏较深，因此Biot模型参数ρ_f，M和b近于0值.土骨架泊松比取为0.3.图 10为计算所需的轮轨不平顺功率谱密度，根据现场实测数据反演得来(王福彤等，2012).为便于比较，模型计算结果也展示在图 9中.从中可见，在远近不同的位置处，在环境振动主要频率范围内，计算结果与试验结果都很接近，说明本文方法能够很好模拟预测场地振动水平.

表 2 Table 2

表 2 机车车辆参数值 Table 2 Parameter values for the vehicles 参数 单位 拖车 动车 车辆质量 kg 4.5×104 5.0×104 转向架簧上质量 kg 2.8×103 3.0×103 转向架簧下质量 kg 1.8×103 1.8×103 车体点头惯量 kg·m2 1.4×106 1.5×106 转向架点头惯量 kg·m2 2.2×103 2.1×103 车辆一系悬挂刚度 N/m 2.5×106 2.5×106 车辆二系悬挂刚度 N/m 1.0×106 1.0×106 车辆一系悬挂阻尼 N·s/m 1.2×105 8.0×104 车辆二系悬挂阻尼 N·s/m 1.1×105 6.0×104 车辆定距 m 12.6 12.6 转向架固定轴距 m 2.0 2.0 车轮半径 m 0.42 0.42 车辆长度 m 19.0 19.0

表 2 机车车辆参数值 Table 2 Parameter values for the vehicles

表 3 Table 3

表 3 轨道系统基本参数 Table 3 Parameter values for the track structure 轨道部件 符号 含义 数值 钢轨 mR 线密度 121 kg·m-1 E0I 弯曲刚度 1.3×107 N·m2 ηR 滞回阻尼比 0.005 轨垫板层 kp0 无阻尼线刚度 3.5×108 N·m-2 ηp 滞回阻尼比 0.075 轨枕层 ms 线密度 416 kg·m-1 道砟路基层 mB 线密度 2400 kg·m-1 kB0 无阻尼线刚度 1.45×107 N·m-2 ηB 滞回阻尼比 0.05 2bt 与场地接触宽度 4.50 m

表 3 轨道系统基本参数 Table 3 Parameter values for the track structure

表 4 Table 4

表 4 地层模型基本参数 Table 4 Parameter values for the ground 层厚/m 剪切波速/(m·s-1) 杨氏模量/MPa 质量密度/(kg·m-3) 3.0 158 117 1800 4.0 169 145 1950 1.0 177 163 2000 1.0 182 168 1950 1.0 187 182 2000 2.0 197 192 1900 1.0 202 218 2050 1.0 209 222 1950 ∞ 225 257 1950

表 4 地层模型基本参数 Table 4 Parameter values for the ground

图 10

Fig. 10

图 10 北京城轨13号线轮轨不平顺谱密度(王福彤等，2012) Fig. 10 Power spectral density for the wheel-rail unevenness of the Line 13，Beijing Urban Railway(Wang et al.，2012)

前述观测试验的局限是没有得到饱和土场地的振动数据，现采用本文模型模拟计算饱和土场地的环境振动.众所周知，饱和软土与饱和密砂的动力学特性存在很大差异.陈龙珠等(1998)提供了饱和软土和饱和砂土的典型参数值，如表 5所示.本算例采 用这两套数值表示两种覆盖层，分别与一个模量较高、孔隙率极小的近似不透水坚硬下卧层搭配，形成两个典型的饱和层状场地.列车和轨道结构参数仍旧保持表 2和表 3中数值不变.

表 5 Table 5

表 5 饱和软土场地和饱和砂土场地的模型参数 Table 5 Parameter values in the models for the grounds of saturated soft clay and saturated sand 参数 单位 软土 砂土 下卧层 排水剪切模量μ0 MPa 19.5 183 462 排水拉梅常数λ0 MPa 77.6 274.5 692 密度ρ kg·m-3 1680 2031 1960 骨架的材料阻尼比η 0.015 0.015 0.015 Biot第一参数M MPa 3500 5600 5600 第二参数α 1.0 1.0 1.0 流体密度ρf kg·m-3 1000 1000 2000 孔隙率φ 0.600 0.375 0.001 黏滞耦合系数b Pa·s·m-2 1×108 1×106 1×109 层厚 m 7 7 ∞

表 5 饱和软土场地和饱和砂土场地的模型参数 Table 5 Parameter values in the models for the grounds of saturated soft clay and saturated sand

图 11和图 12展示了距离轨道中线10、20 m和30 m处各地表点的振动级模拟计算值.图 11的车速为17 m·s^-1，模拟慢速列车; 图 12的车速为101 m·s^-1(约等于表 5软土的120 Hz以内的瑞利波速)，模拟高速列车.对比饱和砂土场地和软土场地可见，无论距轨道远近，在几乎所有频段上，软土场地的振动水平均高于饱和砂土场地.

图 11

Fig. 11

图 11 慢速轨道交通地表点的振动加速度级 Fig. 11 VALs at some points on the ground surface caused by a low speed train

图 12

Fig. 12

图 12 高速铁路地表点的振动加速度级 Fig. 12 VALs at some points on the ground surface caused by a high speed train

多篇文献(Adolfsson et al.，1999; Madshus and Kania 2000; Takemiya，2003; Sheng et al.，2003)报道，瑞典X2000高速列车行经西部海岸线Ledsgard路段的软土场地时，车速接近瑞利波速，场地振动急剧增加.通过对比图 11和图 12的各条曲线可以发现，本文模拟结果反映了这个现象.对于软土场地来说，图 12的车速已达瑞利波速，几乎所有频段上的振动级都明显高于图 11中相应数值.

第5节模型可以计算单一波长的谐波不平顺激励下的地表位移场.设不平顺幅值为1 mm，波长2 m，则对于速度c=17 m·s^-1的慢速列车，轮轨动态相互作用力的频率为17/2=8.5 Hz.图 13a为慢速列车运行引起的砂土场地表面的竖向位移幅值，列车位置在移动坐标x-ct轴的0到-76 m处.整体看饱和砂土场地振动幅值极小，各组车轮正下方的位移略为明显.图 13b为慢速列车下的饱和软土场地，振动强度比砂土场地剧烈得多; 各个轮对激发环状波动，相互干涉，传播很远; 由于车速远小于瑞利波速，列车前方地表可见向前传播的振动.

图 13

Fig. 13

图 13 c=17 m·s-1的列车在轨道上运行引起的地表竖向位移幅值 (a)饱和砂土场地;(a)饱和软土场地. Fig. 13 Vertical displacement of ground surface when a train moving on a track with speed of 17 m·s-1 (a)S and y ground;(b)Soft clay ground.

图 14为高速列车计算结果，此时轮轨力频率为101/2=50.5 Hz.饱和砂土场地与软土场地的振动幅值仍然迥异.50 Hz的移动激扰，在砂土地面激发了向远处传播的环状的波动，如图 14a所示.尽管车速高达101 m·s，仍然远小于饱和砂土的瑞利波速(0~120 Hz频段为278~300 m·s^-1)，所以前方可见领先于列车的表面波.

图 14

Fig. 14

图 14 c=101 m·s-1的列车在轨道上运行引起的地表竖向位移幅值 (a)饱和砂土场地;(a)饱和软土场地. Fig. 14 Vertical displacement of ground surface when a train moving on a track with speed of 101 m·s-1 (a)S and y ground;(b)Soft clay ground.

无论从位移幅值，还是从影响范围看，图 14b的软土场地位移反应都比图 14a的砂土场地剧烈.最大峰值出现在列车前部第一转向架下面，振动向列车侧后方向传播，形成典型的马赫锥.振动局限于马赫锥之内，前方场地几乎没有地表波动.

众所周知，如果移动荷载直接作用于半空间表面，与瑞利波同行达到临界速度，马赫角应该等于90°.图 14b的马赫角明显小于90°，显示车速已经超越了临界速度.这说明，由于轨道的存在，临界速度并不等于场地瑞利波速，而是有所降低.保持其他条件不变，仅将车速降至90 m·s^-1，地表位移波场中仍然出现了马赫锥，如图 15所示.这个结果提醒软土场地工程实践注意，即使设计车速低于场地的瑞利波速，剧烈振动的马赫锥效应也可能出现.

图 15

Fig. 15

图 15 c=90 m·s-1时的地面竖向位移幅值 Fig. 15 Vertical displacement of ground surface when a train speed is 90 m·s-1

动力格林函数是表达场地的数学模型，在列车-轨道-场地系统的耦合振动分析中起到关键作用.本文构建的格林函数，能够考虑饱和岩土体的流固两相动力相互作用，场地分层沉积的几何物理特性以及弹性波在场地中的三维传播.数值算例的结果与文献结果十分接近，与现场观测试验基本符合.模型反映了高速列车的场地振动强度大于低速列车、软土场地振动强度大于饱和砂土场地的规律.当车速接近场地瑞利波速时，基于格林函数计算的地表振动幅值能够反映出马赫锥效应; 即使车速略低于瑞利波速，剧烈振动的马赫锥亦可能出现，这为软土场地高铁工程实践提供了有益的参考.

应该说明的是，还有很多影响因素值得进行数值分析，例如土层厚度，土体的渗透性、孔隙比、密度等对场地振动的影响.但本文意在提出格林函数的计算方法并数值验证，限于篇幅，影响因素将另行撰文讨论.