高速铁路以其速度快、低排放、方便快捷的优点近年得到迅速发展.截至2016年底，我国处于运营状态的高速铁路里程超过2.5万km，占世界67%以上.高速铁路带来便捷的同时，其运行引起的地面振动日益影响到人们的生活、工作以及健康^[1].建立合理的列车荷载下轨道-路基动力分析模型、有效评价高铁运行诱发地面振动对高铁线路的设计和建设具有重要意义，而其中的路基模型很大程度上决定着数值模型的计算量和精度.

最初路基常被视为单相弹性体，学者们对其在移动荷载下的动力响应进行研究.Cole和Huth^[2]最早研究了二维均质弹性地基上作用超高速移动线荷载的动力响应.Eason^[3]研究了移动点荷载、圆形和矩形荷载下三维均质半空间稳态问题.Xia等^[4]采用列车-轨道-路基动力模型分析了列车运行对轨道及周围地面振动的影响.王永刚等^[5]基于移动荷载作用下弹塑性半无限空间的动应力分布和残余应力场分析了交通动应力下路基土体安定性下限值及车辆移动速度的影响.周凤玺等^[6]基于线弹性动力学理论，结合坐标变换用半解析法研究了移动荷载作用下二维非均匀路基动力响应.为了考虑饱和土体中孔隙水对路基振动的影响，Burke等^[7]基于Biot理论^[8-9]首先给出了忽略惯性项的二维多孔饱和半平面在表面移动荷载作用下的解析解.Theodorakopoulos等^[10-11]考虑了水土耦合及平面应变条件，采用半解析法求解了移动荷载下饱和路基动力响应，Jin等^[12]使用半解析方法研究了高速移动荷载下饱和多孔介质的应力和超静孔压.Lefeuve-Mesgouez等^[13]采用Fourier变换的方法分析了竖向矩形简谐移动荷载引起的地面振动.Cai等^[14-16]考虑了路轨系统的影响，获得了列车移动荷载作用下饱和半空间的动力响应，并讨论了荷载移动速度、土体物理参数和渗透系数对路面振动的影响.高广运等^[17-22]利用2.5D有限元系统研究了饱和路基(均匀、分层均匀、横观各向同性和分层横观各向同性)在高铁移动荷载下的动力响应问题.

然而由于大量的蒸发和蒸腾作用路基常处于非饱和状态^[23]，Wu^[24]、Yang等^[25]指出非饱和土饱和度对其动力特性有着明显的影响，因此有必要研究饱和度对移动荷载下地面振动及超静孔隙水压力的影响.相比于三维的有限元分析，2.5D有限元方法在求解移动荷载作用下动力响应问题速度快、精度高，具有很大优势.Yang和Huang^[26]、Takemiya^[27]、Gao等^{[17-22, 28-30]}分别利用2.5D有限元研究了弹性、成层以及饱和路基在移动荷载下的动力响应，指出在求解移动荷载下土体动力响应过程中2.5D有限元法在保证精度的情况下可极大地减小计算量.本文视路基为更符合工程实际的非饱和介质，开发一种非饱和路基的2.5维有限单元方法，研究高速列车移动荷载引起的地面振动.对控制方程时间进行Fourier变换，沿轨道方向波数变换将三维问题降为平面问题，结合边界条件和Galerkin法推导出频域内的2.5维有限元方程.轨道结构视为非饱和路基上的Euler梁，所得频域-波数域内解答通过快速Fourier逆变换得到三维时域-空间域内结果.数值实验研究了车速和路基液体饱和度对地面振动和超静孔隙水压力的影响.

1 计算理论与模型 1.1 路基控制方程的2.5D有限元格式

对非饱和土建模采用如下假定：①表征单元大小满足连续介质尺寸要求；②孔隙水和空气的运动符合广义Darcy定律; ③土颗粒密度与流体压力和体积变形有关，孔隙中水和气的密度分别是孔隙水压力和孔隙气压力的函数.首先定义对t和x的傅里叶变换为：$\tilde{\bar{u}}\left( {{\varepsilon }_{x}},y,z,\omega \right)=\int_{-\infty }^{+\infty }{\int_{-\infty }^{+\infty }{u}}(x,y,z,t){{\text{e}}^{\text{i}{{\varepsilon }_{x}}x}}{{\text{e}}^{-\text{i}\omega t}}\text{d}x\text{d}t $.土骨架质量守恒方程可表示为^[31]

$ \left( {1 - n} \right)\frac{{\partial {\rho _{\rm{s}}}}}{{\partial t}} - {\rho _{\rm{s}}}\frac{{\partial n}}{{\partial t}} + {\rho _{\rm{s}}}\left( {1 - n} \right)\nabla \cdot \dot u = 0 $

(1)

式中：ρ_s为土颗粒密度；n为孔隙率；u为土颗粒位移；“·”表示对时间一阶导数；$\nabla $为拉普拉斯算子.流体压力和体积变形对固体介质密度的影响可表示为^[31]

$ \frac{1}{\rho_{\mathrm{s}}} \frac{\partial p_{\mathrm{c}}}{\partial t}=\frac{\alpha-n}{K_{\mathrm{g}}} \frac{\partial p_{\mathrm{c}}}{\partial t}-(1-\alpha) \frac{\partial \varepsilon_{\mathrm{s}}}{\partial t} $

(2)

式中：K_g为土体颗粒压缩模量，α=1－K_sk/K_g为Biot系数，其中K_sk为土骨架压缩模量；ε_s为土体体积变形; p_c=S_rp^w+(1－S_r)p^a为非饱和流体压力，p^w、p^a分别为孔隙水压力和孔隙气压力，S_r为流体饱和度.将式(2)代入式(1)有

$ \frac{\partial n}{\partial t}=\frac{\alpha-n}{K_{g}} \frac{\partial\left(S_{r} p^{w}+\left(1-S_{r}\right) p^{a}\right)}{\partial t}+(\alpha-n) \nabla \cdot \dot{u} $

(3)

孔隙液体的质量守恒方程为^[31]

$ {S_{\rm{r}}}{\rho _{\rm{w}}}\frac{{\partial n}}{{\partial t}} + n{S_{\rm{r}}}\frac{{\partial {\rho _{\rm{w}}}}}{{\partial t}} + \eta {\rho _{\rm{w}}}\frac{{\partial {S_{\rm{r}}}}}{{\partial t}} + {\rho _{\rm{w}}}n{S_{\rm{r}}}\nabla \cdot {\dot u^{\rm{w}}} = 0 $

(4)

式中：ρ_w为液体密度；u^w为液体的位移.液相密度和液体压力关系为^[23]

$ \frac{{{\rm{d}}{\rho _{\rm{w}}}}}{{{\rho _{\rm{w}}}}} = \frac{{{\rm{d}}{p^{\rm{w}}}}}{{{K_{\rm{w}}}}} $

(5)

式中：K_w为液体压缩模量.基质吸力s=p^a－p^w与液体饱和度S_r的关系可以采用水土特征曲线(SWCC)表示^[23]

$ {S_{\rm{r}}} = {S_{\rm{r}}}\left( s \right) = {S_{\rm{r}}}\left( {{p^{\rm{a}}} - {p^{\rm{w}}}} \right) $

(6)

在小变形假定条件下，将式(3)和式(5)代入式(4)，并结合式(6)可得

$ {A_{11}}{\dot p^{\rm{w}}} + {A_{12}}{\dot p^{\rm{a}}} + {A_{13}}\nabla \dot u + {A_{14}}\nabla {\dot u^{\rm{w}}} = 0 $

(7)

参数A₁₁~A₁₄为$ {A_{11}} = \frac{{(\alpha - n)S_{\rm{r}}^{\rm{2}}}}{{{K_{\rm{g}}}}} + \frac{{n{S_{\rm{r}}}}}{{{K_{\rm{w}}}}} - {A_{ss}}\left( {n - \frac{{(\alpha - n){S_{\rm{r}}}s}}{{{K_{\rm{g}}}}}} \right), $ $ {A_{12}} = \frac{{(\alpha - n){S_{\rm{r}}}\left( {1 - {S_{\rm{r}}}} \right)}}{{{K_{\rm{g}}}}} + {A_{{\rm{ss}}}}\left( {n - \frac{{(\alpha - n){S_{\rm{r}}}s}}{{{K_{\rm{g}}}}}} \right), $ $ {A_{13}} = {S_{\rm{r}}}(\alpha - n), {A_{14}} = n{S_{\rm{r}}}, {A_{{\rm{ss}}}} = \frac{{\partial {S_{\rm{r}}}}}{{\partial s}}.$

孔隙气体的质量守恒方程为^[31]

$ \begin{array}{*{20}{c}} {\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}\frac{{\partial n}}{{\partial t}} + n\left( {1 - {S_{\rm{r}}}} \right)\frac{{\partial {\rho _{\rm{a}}}}}{{\partial t}} - n{\rho _{\rm{a}}}\frac{{\partial {S_{\rm{r}}}}}{{\partial t}} + }\\ {n\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}\nabla \cdot {{\dot u}^{\rm{a}}} = 0} \end{array} $

(8)

式中:ρ_a为气体密度；u^a为气体位移.气相压缩变形和气相压力关系为^[23]

$ \frac{{{\rm{d}}{\rho _{\rm{a}}}}}{{{\rho _{\rm{a}}}}} = \frac{{{\rm{d}}{p^{\rm{a}}}}}{{{p^{\rm{a}}}}} $

(9)

在小变形假定条件下，将式(3)和式(9)代入式(8)，并结合式(6)有

$ {A_{21}}{\dot p^{\rm{w}}} + {A_{22}}{\dot p^{\rm{a}}} + {A_{23}}\nabla \dot u + {A_{24}}\nabla {\dot u^{\rm{a}}} = 0 $

(10)

$ {A_{21}} = \frac{{\left( {\alpha - n} \right){S_{\rm{r}}}\left( {1 - {S_{\rm{r}}}} \right)}}{{{K_{\rm{g}}}}} + {A_{{\rm{ss}}}}\left( {n + \frac{{\left( {\alpha - n} \right)\left( {1 - {S_{\rm{r}}}} \right)s}}{{{K_{\rm{g}}}}}} \right), $

$ \begin{array}{*{20}{c}} {{A_{22}} = \frac{{\left( {\alpha - n} \right){{\left( {1 - {S_{\rm{r}}}} \right)}^2}}}{{{K_{\rm{g}}}}} + \frac{{n\left( {1 - {S_{\rm{r}}}} \right)}}{{{p^{\rm{a}}}}} - }\\ {{A_{{\rm{ss}}}}\left( {n + \frac{{\left( {\alpha - n} \right)\left( {1 - {S_{\rm{r}}}} \right)s}}{{{K_{\rm{g}}}}}} \right),} \end{array} $

$ {A_{23}} = \left( {1 - {S_{\rm{r}}}} \right)\left( {\alpha - n} \right),{A_{24}} = n\left( {1 - {S_{\rm{r}}}} \right). $

根据广义Darcy定律，孔隙流体(液相和气相)的渗流运动方程分别为^[23]

$ n{S_{\rm{r}}}\left( {\dot u_i^{\rm{w}} - {{\dot u}_i}} \right) = - \frac{{{k_{\rm{w}}}}}{{{\rho _{\rm{w}}}g}}\left( {p_{,i}^{\rm{w}} + {\rho _{\rm{w}}}\ddot u_i^{\rm{w}}} \right) $

(11)

$ n\left( {1 - {S_{\rm{r}}}} \right)\left( {\dot u_i^{\rm{a}} - {{\dot u}_i}} \right) = - \frac{{{k_{\rm{a}}}}}{{{\rho _{\rm{a}}}g}}\left( {p_{,i}^{\rm{a}} + {\rho _{\rm{a}}}\ddot u_i^{\rm{a}}} \right) $

(12)

式中“··”表示对时间二阶导数.采用傅里叶变换处理式(11)和式(12)，得孔隙水和孔隙气平均位移在频率域内表达式分别为

$ \tilde u_i^{\rm{w}} = \left( {{F_{\rm{w}}}{{\tilde u}_i} - \tilde p_{,i}^{\rm{w}}} \right)/\left( {{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}} \right) $

(13)

$ \tilde u_i^{\rm{a}} = \left( {{F_{\rm{a}}}{{\tilde u}_i} - \tilde p_{,i}^{\rm{a}}} \right)/\left( {{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}} \right) $

(14)

式中“~”代表参数频域内的解；$ F_{\mathrm{w}}=n S_{\mathrm{r}} \rho_{\mathrm{w}} g \omega i / k_{\mathrm{w}}$；$F_{\mathrm{a}}=n\left(1-S_{\mathrm{r}}\right) \rho_{\mathrm{a}} g \omega i / k_{\mathrm{a}}, \omega $为振动圆频率.采用傅里叶变换处理式(7)和式(10)，结合(13)和式(14)得频域内液体和气体的质量守恒方程为

$ \begin{array}{*{20}{c}} {\left( {{A_{13}} + \frac{{{A_{14}}{F_{\rm{w}}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}}} \right){{\tilde u}_{i,i}} - \frac{{{A_{14}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}}\tilde p_{,ii}^{\rm{w}} + }\\ {{A_{11}}{{\tilde p}^{\rm{w}}} + {A_{12}}{{\tilde p}^{\rm{a}}} = 0} \end{array} $

(15)

$ \begin{array}{*{20}{c}} {\left( {{A_{23}} + \frac{{{A_{24}}{F_{\rm{a}}}}}{{{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}}}} \right){{\tilde u}_{i,i}} - \frac{{{A_{24}}}}{{{F_a} - {\rho _a}{\omega ^2}}}\tilde p_{,ii}^{\rm{a}} + }\\ {{A_{21}}{{\tilde p}^{\rm{w}}} + {A_{22}}{{\tilde p}^{\rm{a}}} = 0} \end{array} $

(16)

由非饱和多孔介质平衡方程$ {\sigma _{ij, j}} + {F_i} = {{\bar \rho }_s}{{\ddot u}_i} + {{\bar \rho }_{\rm{w}}}\ddot u_i^{\rm{w}} + {{\bar \rho }_{\rm{a}}}\ddot u_i^{\rm{a}}$、小变形条件下应力-应变关系式$ {\sigma _{ij}} = 2\mu {\varepsilon _{ij}} + \lambda {\delta _{ij}}{\varepsilon _{kk}} - {\delta _{ij}}\alpha {p_{\rm{c}}}$和弹性介质几何方程${\varepsilon _{ij}} = \frac{1}{2}\left( {{u_{i, j}} + {u_{j, i}}} \right) $可得非饱和多孔介质动量守恒方程为^[23]

$ \begin{array}{*{20}{c}} {\mu {u_{i,jj}} + \left( {\lambda + \mu } \right){u_{i,jj}} - {S_r}p_{,i}^{\rm{w}} - \left( {1 - {S_r}} \right)p_{,i}^{\rm{a}} = }\\ {{{\bar \rho }_{\rm{s}}}{{\ddot u}_i} + {{\bar \rho }_{\rm{w}}}\ddot u_i^{\rm{w}} + {{\bar \rho }_{\rm{a}}}\ddot u_i^{\rm{a}}} \end{array} $

(17)

式中:λ和μ为土体Lame常数；采用复数形式的Lame常数$ \bar \lambda = \left( {1 + 2i{\eta _{\rm{s}}}} \right)\lambda , \bar \mu = \left( {1 + 2i{\eta _{\rm{s}}}} \right)\mu $考虑土体阻尼系数的影响；η_s为阻尼系数；ε_ij为土体的应变分量；δ_ij为Kronecker符号；σ_ij和F_i为土体的应力分量和体力分量；$\overline{\rho}_{\mathrm{s}}=(1-n) \rho_{\mathrm{s}}、\overline{\rho}_{\mathrm{w}}=n \mathrm{S}_{\mathrm{r}} \rho_{\mathrm{w}}、\overline{\rho}_{\mathrm{a}}=n\left(1-S_{\mathrm{r}}\right) \rho_{\mathrm{a}} $分别为固相、液相和气相的相对密度.

采用傅里叶变换处理式(17)，结合式(13)和式(14)，用应力-应变关系式消去应力项，得到非饱和土体频域内动力方程如下:

$ \begin{array}{l} \bar \mu {{\tilde u}_{i,jj}} + \left( {\bar \lambda + \bar \mu } \right){{\tilde u}_{i,ji}} - {S_{\rm{r}}}\tilde p_{,i}^{\rm{w}} - \left( {1 - {S_{\rm{r}}}} \right)\tilde p_{,i}^{\rm{a}} + \\ {\omega ^2}\left[ {\left( {1 - n} \right){\rho _{\rm{s}}} + \frac{{n{S_{\rm{r}}}{\rho _{\rm{w}}}{F_{\rm{w}}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}} + \frac{{n\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}{F_{\rm{a}}}}}{{{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}}}} \right]{{\tilde u}_i} - \\ \frac{{{\omega ^2}n{S_{\rm{r}}}{\rho _{\rm{w}}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}}\tilde p_{,i}^{\rm{w}} - \frac{{{\omega ^2}n\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}}}{{{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}}}\tilde p_{,i}^{\rm{a}} = 0 \end{array} $

(18)

非饱和多孔介质控制方式(15)、式(16)和式(18)可以采用有限元法求解，对排水排气边界，应力、流体孔压及流量边界条件可分别表示为

$ \sigma_{i j}^{\mathrm{s}} n_{j}=f_{i}, p=p^{\prime},-k_{\mathrm{d}} p_{, j}^{\prime} n_{j}=\rho_{\mathrm{f}} g v_{\mathrm{n}}=q $

式中:f_i为边界处应力；p′为边界处孔隙流体压力；v_n, q分别表示边界处流速和流量；ρ_f为流量密度；k_d为渗透系数；g为重力加速度；n_j为方向向量.由式(18)和应力边界条件，引入虚位移δu_i^*和频域内有效应力${{\tilde \sigma '}_{ij, j}} $，利用虚功原理及高斯公式得多孔介质控制方程如下:

$ \begin{array}{l} \int {\left[ {\delta \varepsilon _i^ * {{\tilde \sigma '}_{ij}} - \delta u_i^ * {\omega ^2}\left[ {\left( {1 - n} \right){\rho _{\rm{s}}} + \frac{{n{S_{\rm{r}}}{\rho _{\rm{w}}}{F_{\rm{w}}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}} + } \right.} \right.} \\ \left. {\left. {\frac{{n\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}{F_{\rm{a}}}}}{{{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}}}} \right]{{\tilde u}_i}} \right]{\rm{d}}V + \int {\left( {\delta u_i^ * \frac{{{\omega ^2}n{S_{\rm{r}}}{\rho _{\rm{w}}}}}{{{F_{\rm{w}}} - {\rho _{\rm{w}}}{\omega ^2}}}\tilde p_{,i}^{\rm{w}} - } \right.} \\ \left. {\delta \varepsilon _i^ * {\delta _{ij}}{S_{\rm{r}}}\tilde p_{,i}^{\rm{w}}} \right){\rm{d}}V + \int {\left( {\delta u_i^ * \frac{{{\omega ^2}n\left( {1 - {S_{\rm{r}}}} \right){\rho _{\rm{a}}}}}{{{F_{\rm{a}}} - {\rho _{\rm{a}}}{\omega ^2}}}\tilde p_{,i}^{\rm{a}} - } \right.} \\ \left. {\delta \varepsilon _i^ * {\delta _{ij}}\left( {1 - {S_{\rm{r}}}} \right)\tilde p_{,i}^{\rm{a}}} \right){\rm{d}}V = \int {\delta u_i^ * {f_i}{\rm{d}}S} \end{array} $

(19)

对式(19)在列车运行方向采用波数展开，并采用4节点等参单元离散.单元形函数取$ {N_i}(\eta , \xi ) = \frac{1}{4}\left( {1 + {\eta _i}\eta } \right)\left( {1 + {\xi _i}\xi } \right)$，单元位移和孔隙压力与节点位移、节点孔隙压力关系为

$ {u_i} = \sum\limits_{i = 1}^4 {{N_j}\left( {\eta ,\xi } \right)u_{ij}^{\rm{e}}} ,p = \sum\limits_{i = 1}^4 {{N_j}\left( {\eta ,\xi } \right)p_i^{\rm{e}}} $

式中：η、ξ分别为单元局部坐标；η_i、ξ_i为节点局部坐标，则可得到矩阵形式2.5维有限元控制方程为

$ \begin{array}{*{20}{c}} {\left( {{{\mathit{\boldsymbol{K'}}}_{{\rm{up}}}} - {\mathit{\boldsymbol{M}}_{{\rm{up}}}}} \right)\mathit{\boldsymbol{\tilde {\bar u}}} + \left( {{{\mathit{\boldsymbol{Q'}}}_{{\rm{up}}}} - {\mathit{\boldsymbol{Q}}_{{\rm{up}}}}} \right){{\mathit{\boldsymbol{\tilde {\bar p}}}}^{\rm{w}}} + }\\ {\left( {{{\mathit{\boldsymbol{G'}}}_{{\rm{up}}}} - {\mathit{\boldsymbol{G}}_{{\rm{up}}}}} \right){{\mathit{\boldsymbol{\tilde {\bar p}}}}^{\rm{g}}} = \mathit{\boldsymbol{\tilde {\bar f}}}_{{\rm{up}}}^{\rm{s}}} \end{array} $

(20)

式中：“-”代表参数在波数域中的量；K′_up，M_up，Q_up，G_up，f_up^s分别为刚度矩阵、质量矩阵、液体贡献矩阵、气体贡献矩阵和等效节点力矩阵.各参量表达式如下:

$ \begin{gathered} {{\mathit{\boldsymbol{K'}}}_{{\text{up}}}} = \sum\limits_{\text{e}} {\iint {{{\left( {{\mathit{\boldsymbol{B}}^ * }\mathit{\boldsymbol{N}}} \right)}^{\text{T}}}\mathit{\boldsymbol{D}}\left( {\mathit{\boldsymbol{BN}}} \right)\left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} ;{\mathit{\boldsymbol{M}}_{{\text{up}}}} = \hfill \\ \;\;\;\;\;\;{\omega ^2}\left[ {\left( {1 - n} \right){\rho _{\text{s}}} + n{S_{\text{r}}}{\rho _{\text{w}}}\frac{{{F_{\text{w}}}}}{{{F_{\text{w}}} - {\rho _{\text{w}}}{\omega ^2}}} + } \right. \hfill \\ \;\;\;\;\;\;\left. {n\left( {1 - {S_{\text{r}}}} \right){\rho _{\text{w}}}\frac{{{F_{\text{a}}}}}{{{F_{\text{a}}} - {\rho _{\text{a}}}{\omega ^2}}}} \right]\sum\limits_{\text{e}} {\iint {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^{\text{T}}}\left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} ; \hfill \\ \end{gathered} $

$ {\mathit{\boldsymbol{Q'}}_{up}} = \frac{{{\omega ^2}n{S_{\text{r}}}{\rho _{\text{w}}}}}{{{F_{\text{w}}} - {\rho _{\text{w}}}{\omega ^2}}}\sum\limits_{\text{e}} {\iint {{\mathit{\boldsymbol{N}}^{\text{T}}}\overline {\mathit{\boldsymbol{BN}}} \left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} ; $

$ {\mathit{\boldsymbol{Q}}_{up}} = \alpha {S_{\text{r}}}\sum\limits_{\text{e}} {\iint {{{\left( {{\mathit{\boldsymbol{B}}^*}\mathit{\boldsymbol{N}}} \right)}^{\text{T}}}}} \mathit{\boldsymbol{m}}\mathit{\boldsymbol{\overline N}} |J|{\text{d}}\eta {\text{d}}\xi ; $

$ {\mathit{\boldsymbol{G'}}_{{\text{up}}}} = \frac{{{\omega ^2}n\left( {1 - {S_r}} \right){\rho _{\text{a}}}}}{{{F_{\text{a}}} - {\rho _{\text{a}}}{\omega ^2}}}\sum\limits_{\text{e}} {\iint {{N^{\text{T}}}\overline {\mathit{\boldsymbol{BN}}} \left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} ; $

$ {\mathit{\boldsymbol{G}}_{{\text{up}}}} = \alpha \left( {1 - {S_{\text{r}}}} \right)\sum\limits_{\text{e}} {\iint {{{\left( {{\mathit{\boldsymbol{B}}^ * }\mathit{\boldsymbol{N}}} \right)}^{\text{T}}}\mathit{\boldsymbol{m\bar N}}\left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} ; $

$ \mathit{\boldsymbol{\tilde {\bar f}}}_{up}^{\text{s}} = \sum\limits_{\text{e}} {\iint {{\mathit{\boldsymbol{N}}^T}\mathit{\boldsymbol{\tilde {\bar f}}}\left| \mathit{\boldsymbol{J}} \right|{\text{d}}\eta {\text{d}}\xi }} $

式中：“*”代表矩阵共轭；N, $ \overline{\boldsymbol{N}}$为空间域和波数域形函数；J为雅可比矩阵；D为应力矩阵；B为偏导矩阵.对式(15)和式(16)采用Galerkin方法并结合流量和孔压边界条件可得相应有限元控制方程如下:

$ \left( {{\mathit{\boldsymbol{H}}_{{\text{md}}}}} \right)\mathit{\boldsymbol{\tilde {\bar u}}} + \left( {{\mathit{\boldsymbol{Q}}_{{\text{md}}}} + {{\mathit{\boldsymbol{Q'}}}_{{\text{md}}}}} \right){\mathit{\boldsymbol{\tilde {\bar p}}}^w} + \left( {{{\mathit{\boldsymbol{G'}}}_{{\text{md}}}}} \right){\mathit{\boldsymbol{\tilde {\bar p}}}^{\text{g}}} = {\mathit{\boldsymbol{\tilde {\bar f}}}^{\text{w}}} $

(21)

$ \left( {{\mathit{\boldsymbol{H}}_{{\text{dw}}}}} \right)\mathit{\boldsymbol{\tilde {\bar u}}} + \left( {{{\mathit{\boldsymbol{Q'}}}_{{\text{dw}}}}} \right){\mathit{\boldsymbol{\tilde {\bar p}}}^w} + \left( {{\mathit{\boldsymbol{G}}_{{\text{dw}}}} + {{\mathit{\boldsymbol{G'}}}_{{\text{dw}}}}} \right){\mathit{\boldsymbol{\tilde {\bar p}}}^{\text{g}}} = {\mathit{\boldsymbol{\tilde {\bar f}}}^{\text{g}}} $

(22)

$ {\mathit{\boldsymbol{H}}_{{\text{md}}}} = \left( {{A_{13}} + \frac{{{A_{14}}{F_{\text{w}}}}}{{{F_{\text{w}}} - {\rho _{\text{w}}}{\omega ^2}}}} \right)\sum\limits_{\text{e}} {\iint {{{\mathit{\boldsymbol{\overline N}} }^{\text{T}}}}} {\mathit{\boldsymbol{m}}^{\text{T}}}\mathit{\boldsymbol{BN}}|\mathit{\boldsymbol{J}}|{\text{d}}\eta {\text{d}}\xi ; $

$ {\mathit{\boldsymbol{H}}_{{\text{dw}}}} = \left( {{A_{23}} + \frac{{{A_{24}}{F_{\text{a}}}}}{{{F_{\text{a}}} - {\rho _{\text{a}}}{\omega ^2}}}} \right)\sum\limits_{\text{e}} {\iint {{{\mathit{\boldsymbol{\overline N}} }^{\text{T}}}}} {\mathit{\boldsymbol{m}}^{\text{T}}}\mathit{\boldsymbol{BN}}|\mathit{\boldsymbol{J}}|{\text{d}}\eta {\text{d}}\xi ; $

$ \boldsymbol{Q}_{\mathrm{md}}=\frac{A_{14}}{F_{\mathrm{w}}-\rho_{\mathrm{w}} \omega^{2}} \sum\limits_{\mathrm{e}} \iint\left(\boldsymbol{B}_{\mathrm{s}}^{*} \overline{\boldsymbol{N}}\right)^{\mathrm{T}}\left(\boldsymbol{B}_{\mathrm{s}} \overline{\boldsymbol{N}}\right)|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi; $

$ Q_{\mathrm{md}}^{\prime}=A_{11} \sum\limits_{\mathrm{e}} \iint \overline{\boldsymbol{N}}^{\mathrm{T}} \overline{\boldsymbol{N}}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi; $

$ \boldsymbol{G}_{\mathrm{dw}}=\frac{A_{24}}{F_{\mathrm{a}}-\rho_{\mathrm{a}} \omega^{2}} \sum\limits_{\mathrm{e}} \iint\left(\boldsymbol{B}_{\mathrm{s}}^{*} \overline{\boldsymbol{N}}\right)^{\mathrm{T}}\left(\boldsymbol{B}_{\mathrm{s}} \overline{\boldsymbol{N}}\right)|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi; $

$ \boldsymbol{Q}_{\mathrm{dw}}^{\prime}=A_{21} \sum\limits_{\mathrm{e}} \iint \overline{\boldsymbol{N}}^{\mathrm{T}} \overline{\boldsymbol{N}}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi; $

$ \boldsymbol{G}_{\mathrm{md}}^{\prime}=A_{12} \sum\limits_{\mathrm{e}} \iint \overline{\boldsymbol{N}}^{\mathrm{T}} \overline{\boldsymbol{N}}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi $

$ \boldsymbol{G}_{\mathrm{dw}}^{\prime}=A_{22} \sum\limits_{\bf{e}} \iint \overline{\boldsymbol{N}}^{\mathrm{T}} \overline{\boldsymbol{N}}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{d} \xi; $

$ \mathit{\boldsymbol{\tilde {\bar f}}}_{{\text{up}}}^{\text{w}} = - \sum\limits_{\text{e}} {\iint {\frac{{{F_{\text{w}}}{\rho _{\text{w}}}g{{\mathit{\boldsymbol{\tilde v}}}_{\text{n}}}}}{{{k_{\text{w}}}}}}} |\mathit{\boldsymbol{J}}|{\text{d}}\eta {\text{d}}\xi ; $

$ \mathit{\boldsymbol{\tilde {\bar f}}}_{{\text{up}}}^{\text{g}} = - \sum\limits_{\text{e}} {\iint {\frac{{{F_{\text{g}}}{\rho _{\text{g}}}g{{\mathit{\boldsymbol{\tilde v}}}_{\text{n}}}}}{{{k_{\text{g}}}}}}} |\mathit{\boldsymbol{J}}|{\text{d}}\eta {\text{d}}\xi . $

综合式(20)、式(21)和式(22)可以得到非饱和地基2.5维有限元控制方程为

$ \boldsymbol{K U}=\boldsymbol{R} $

(23)

$ \begin{array}{l} \mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{K'}}}_{{\rm{up}}}} - {\mathit{\boldsymbol{M}}_{{\rm{up}}}}}&{{{\mathit{\boldsymbol{Q'}}}_{{\rm{up}}}} - {\mathit{\boldsymbol{Q}}_{{\rm{up}}}}}&{{{\mathit{\boldsymbol{G'}}}_{{\rm{up}}}} + {\mathit{\boldsymbol{G}}_{{\rm{up}}}}}\\ {{\mathit{\boldsymbol{H}}_{{\rm{md}}}}}&{{\mathit{\boldsymbol{Q}}_{{\rm{md}}}} + {{\mathit{\boldsymbol{Q'}}}_{{\rm{md}}}}}&{{{\mathit{\boldsymbol{G'}}}_{{\rm{md}}}}}\\ {{\mathit{\boldsymbol{H}}_{{\rm{dw}}}}}&{{{\mathit{\boldsymbol{Q'}}}_{{\rm{dw}}}}}&{{{\mathit{\boldsymbol{G'}}}_{{\rm{dw}}}} + {\mathit{\boldsymbol{G}}_{{\rm{dw}}}}} \end{array}} \right],\mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde {\bar u}}}}\\ {{{\mathit{\boldsymbol{\tilde {\bar p}}}}^{\rm{w}}}}\\ {{{\mathit{\boldsymbol{\tilde {\bar p}}}}^{\rm{g}}}} \end{array}} \right],\\ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde {\bar f}}}_{{\rm{up}}}^{\rm{s}}}\\ {{{\mathit{\boldsymbol{\tilde {\bar f}}}}^{\rm{w}}}}\\ {{{\mathit{\boldsymbol{\tilde {\bar f}}}}^{\rm{g}}}} \end{array}} \right]. \end{array} $

1.2 轨道模型

采用Euler梁模型模拟轨道系统，根据结构力学理论，波数-频率域轨道动力方程如下^[28]：

$ \left( {EI\xi _x^4 - m{\omega ^4}} \right)u_r^{xt} = f_{IT}^{xt}\left( {{\xi _x},\omega } \right) + p_0^{xt}\left( {{\xi _x},\omega } \right) $

(24)

式中：u_r^xt为波数-频率域位移；EI为轨道的弯曲刚度；m为轨道和枕木的综合质量；f_IT^xt为波数域-频率域路基接触反力；p₀^xt表示波数-频率域轨道上动荷载，其中p₀为列车车轮轴重，δ(x)为迪克拉函数；ξ_x为波数变量.列车作用在轨道上的连续轴重荷载同文献[27].

本文计算模型示意图如图 1所示，采用黏弹性吸波边界防止外行波在边界的反射.由于边界处距离振源较远，非饱和路基边界处孔隙气压力可忽略不计，黏弹性吸波边界推导详见文献[29].编制程序求解轨道路基耦合控制方程, 式(23)和式(24)可得到方程在频域-波数域内解答，将结果通过快速傅里叶逆变换(FFT)可以得到时间空间域内的解答.

2 模型验证与计算参数

本文地基地层分布和土体参数采用瑞典国家铁路局X2000列车现场实测数据，实测地层、列车轴重参数以及轨道参数见表 1和表 2^[27].值得注意的是，本文主要参考的非饱和土理论(文献[31])，在S_r=1的情况下，其退化为Biot波动方程.首先对模型进行退化验证，取饱和度S_r=0土体退化为弹性土体进行计算.图 2给出列车速度为70 km·h^-1和200 km·h^-1时轨道中心处路基地面振动结果，同时给出了与文献[27]的对比.由图可以看出两者相当一致，说明本文模型可靠.

非饱和路基动力响应计算中，土水土特征曲线采用V-G模型^[31]：$ {S_{\rm{e}}} = {\left[ {1 + {{\left( {{\alpha _1}s} \right)}^k}} \right]^{ - m}}$；孔隙气和孔隙水的渗透系数采用Mualem模型^[31]：$ {k_{\rm{w}}} = \frac{{{\rho _{\rm{w}}}g\kappa }}{{{\eta _{\rm{w}}}}}\sqrt {{{\rm{S}}_{\rm{r}}}} {\left\{ {1 - {{\left[ {1 - {{\left( {{{\rm{S}}_{\rm{e}}}} \right)}^{\frac{1}{m}}}} \right]}^m}} \right\}^2}, {k_{\rm{g}}} = \frac{{{\rho _{\rm{g}}}g\kappa }}{{{\eta _{\rm{g}}}}}\sqrt {1 - {{\rm{S}}_{\rm{r}}}} {[1 - {\left( {{S_{\rm{e}}}} \right)^{\frac{1}{m}}}]^{2m}}$，其中参数α₁=1×10^-5，k=4，m= 1－1/k=0.75，S_e=(S_r－S_wo)/(1－S_wo)为有效饱和度，S_wo=0.05为水的束缚饱和度，s为基质吸力.η_w=1.0 mPa·s为孔隙水黏性系数，η_g=0.015 mPa·s为孔隙气体的黏性系数；κ=1.0 μm²为土体的固有渗透率；孔隙水压缩模量K_w=2.1×109 N·m^-2.流体的物理参数在动力作用过程中假定为常数；地面边界设为排水边界和排气边界，地层底边和无穷远边设为不排水和不排气边界，以大气压力为参考气压力.高铁多分布于东南沿海地区，本文选取较高饱和度(100%、99%和90%)进行分析.

3 不同车速下非饱和路基地面振动位移时程分析

分析轨道中心处和距离轨道中心8 m远处(图 1中点A)路基地面振动位移随车速和饱和度变化.车速分别取200、250、300和350 km·h^-1，每个车速下路基饱和度取100%和90%.

3.1 轨道中心处分析

图 3为不同速度和饱和度下轨道中心处地面竖向振动位移时程曲线.在车速不超过250 km·h^-1时，非饱和路基地面位移振动幅值明显小于饱和路基；随着速度的增大，饱和与非饱和路基地面振动位移幅值的差距呈减小趋势.值得注意的是饱和度从100%降低为90%时，饱和土中进入少量气体，土体从完全饱和变为非饱和状态，带来土体相状态的彻底改变(三相变为两相)，进而带来土体动力特性的实质性改变，非饱和路基中的孔隙气体降低了轨道中心处的地面位移振幅，而完全饱和是振动位移最不利状态.

与速度低于250 km·h^-1时情况相反，300 km·h^-1和350 km·h^-1时饱和度的减小对位移幅值的增大作用有限，随饱和度的降低路基地面振动位移幅值稍有增大，但350 km·h^-1时饱和度的改变对振动位移幅值的影响已不明显：车速为300 km·h^-1时, 饱和度为100%和90%时地面振动峰值为3.5 mm和4.7 mm；车速为350 km·h^-1时, 饱和度为100%和90%时地面振动峰值为4.0 mm和4.7 mm，此时随速度增大，饱和和非饱和路基地面振动位移幅值接近相等.从时程图上也可以看到，相比于S_r=100%的饱和土，S_r=90%的非饱和土地面振动位移随时间衰减地更快；并且随车速增大，非饱和路基地面振动位移持时变短.

3.2 轨道中心8 m远处时程分析

图 4为不同车速和饱和度下距轨道中心8 m远处地面振动位移时程曲线.各车速下，饱和路基地面位移振幅峰值在0.7~1.2 mm之间，饱和度为90%的路基地面振动位移幅值相近，约为1~1.5 mm；与3.1节中轨道中心处的情况相反，在车速未超过250 km·h^-1时，距轨道中心8 m远处非饱和路基地面振动位移幅值大于饱和路基，这可能是由于非饱和路基中孔隙气的存在增加了振动波在传播过程中的散射和折射，孔隙气增大了非饱和路基远处的地面振动加速度；对于饱和路基，车速从200 km·h^-1增大到350 km·h^-1其地面振动位移幅值逐渐增大，但增大的幅度越来越小(车速为200、250、300和350 km·h^-1时，饱和路基地面位移振幅峰值为0.70、0.78、1.0和1.2 mm).车速高于300 km·h^-1后，随车速提高饱和及非饱和路基地面振动位移幅值均增大且两者逐渐趋于相等，且350 km·h^-1时饱和路基振动位移幅值有超过非饱和土路基的趋势.在同一速度下，饱和路基地面位移振动持时大于非饱和路基；不同速度下，随着车速提高，非饱和土振动持时变短，而饱和土地面振动持续时间变长.

4 不同车速下非饱和土路基地面振动衰减特性

图 5为饱和度为90%路基地面振动位移幅值(图 5a)和加速度幅值(图 5b)随轨道中心距离衰减曲线，车速取200、250、300和350 km·h^-1.200 km·h^-1车速下的地面振动位移峰值在近轨道处(距离小于5 m时)大于其余速度.近轨道处各速度下地面振动位移峰值均快速衰减且衰减速度相当；远轨道处(距离超过5 m)200 km·h^-1车速下衰减最快，其余车速下地面振动位移振幅大小相当且几乎不衰减.由图 5b可以看出，同一位置处低速下的加速度级大于高速下的加速度级，且距轨道中心2 m以内均快速衰减.距轨道中心2 m以外车速为300 km·h^-1时地面振动加速度衰减曲线平缓，衰减最慢，而200 km·h^-1车速下衰减最快.地面振动加速度级在某些列车速度下的衰减会出现反弹增大现象，250 km·h^-1时地面振动加速度反弹增大现象出现在距轨道中心约3 m远处，350 km·h^-1时出现多次地面振动加速度衰减的反弹增大现象，距轨道中心约4、6和8 m.可见加速度级的反弹增大现象出现与否及其位置与车速密切相关.

5 不同车速下非饱和土路基超静孔隙水压力

图 6为不同速度下轨道中心处超静孔隙水压力随深度变化曲线，车速为200、250、300和350 km·h^-1，路基土饱和度取100%和90%.可以发现不同车速下饱和土与非饱和土超静孔压变化规律类似，均在轨道中心下1.5~2.0 m之间达最大值，之后再3.5~4 m之间迅速减小，4.5 m之下超静孔隙水压力已经很小.即超静孔隙水压力主要在浅层分布，分布深度约为0~4.5 m.同一车速下饱和度对超静孔隙水压力具有较大影响，饱和土超静孔隙水压力明显大于饱和度为90%情况，随路基饱和度下降超静孔隙水压力峰值显著降低.饱和度90%的路基轨道中心下超静孔隙水压力峰值约为1 kPa，因此对于饱和度较低的路基，高铁荷载产生的超静孔隙水压力可以略去不计.饱和度为100%时随车速增大超静孔隙水压力有稍减小的趋势.4.5 m以下车速对超静孔隙水压力影响较小，各车速下饱和度100%路基轨道中心4.5 m深处孔压约为0.6 kPa，10 m深处约为0.2 kPa；饱和度90%路基4.5 m深处孔压约为0.1 kPa，10 m处仅约为0.02 kPa(分别为同一深度处饱和土路基的17%和10%).

6 结论

本文利用非饱和路基2.5维有限单元方法研究了高速列车移动荷载引起的非饱和路基地面振动，主要结论如下：

(1) 在轨道中心处：列车速度不超过250 km·h^-1时，非饱和土路基地面位移振动幅值明显小于饱和路基，路基从完全饱和到非饱和状态轨道中心处地面振动位移显著减小，完全饱和是振动最不利状态；在车速小于等于250 km·h^-1时，随着车速增大，饱和与非饱和路基地面位移幅值差值呈减小趋势；车速为350 km·h^-1时位移幅值几近相等.相比于饱和路基，非饱和路基地面振动随时间衰减更快，持时更短.

(2) 距离轨道中心8 m远处：车速不超过250 km·h^-1时非饱和路基振动位移幅值明显大于饱和路基，车速高于300 km·h^-1后随车速提高饱和与非饱和路基地面振动位移幅值增大且两者逐渐趋于接近，在350 km·h^-1时饱和路基地面振动位移幅值有超过非饱和土路基的趋势.在同一速度下，距离轨道中心8 m远处饱和路基路面振动持时大于非饱和路基；车速高于250 km·h^-1后随车速增大非饱和土振动持续时间变短，而饱和土地面振动持时变长.

(3) 近轨道处(距离小于5 m)位移幅值低速大于高速且均以大小相当的速率快速衰减；远轨道处(距离超过5 m)200 km·h^-1车速下衰减最快，其余车速下地面振动位移振幅大小相当且几乎不衰减.同一位置处低速下地面振动加速度级大于高速下的加速度级，且距轨道中心2 m以内均快速衰减.地面振动加速度级在某些列车速度下的衰减会出现反弹增大现象，反弹增大现象的出现与否及其位置与车速密切相关.

(4) 轨道中心下饱和与非饱和路基超静孔隙水压力分布深度约为0~4.5 m，在轨道中心下1.5~2.0 m之间达最大值，在3.5~4 m深度内急剧减小；同一车速下土体由饱和变为非饱和时超静孔隙水压力峰值显著减小，对于饱和度较低的路基，高铁荷载产生的超静孔隙水压力可以略去不计.