基于考虑变形效应的土-水特征曲线求解非饱和地基动力响应

胡静; 金林廉; 吕志豪; 张家康; 边学成

doi:10.11779/CJGE20230800

基于考虑变形效应的土-水特征曲线求解非饱和地基动力响应 English Version

1.
福州大学土木工程学院, 福建福州 350108
2.
浙江大学建筑工程学院, 浙江杭州 310058

基金项目:

国家自然科学基金项目 52108308

国家重点研发计划青年科学家项目 2024YFC2911000

详细信息

作者简介:
胡静(1991—)，博士，副研究员，从事交通荷载作用下软土地基动力响应方面的研究工作。E-mail：jingh@fzu.edu.cn

中图分类号: TU435
计量
- 文章访问数: 238
- HTML全文浏览量: 21
- PDF下载量: 44
出版历程
- 收稿日期: 2023-08-17
- 网络出版日期: 2024-06-04
- 刊出日期: 2025-01-31

Dynamic response solutions of unsaturated soil foundation using soil-water characteristic curve considering deformation effects

1.
College of Civil Engineering, Fuzhou University, Fuzhou 350108, China
2.
College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310508, China

摘要

摘要: 为研究移动荷载作用下的非饱和地基动力响应，在传统V-G土-水特征曲线模型的基础上，建立了考虑变形效应的土-水特征曲线模型，并基于该模型推导了新的非饱和土动力控制方程，从而完整描述非饱和土在受动力作用下的水-力耦合作用。进一步的，采用2.5维有限元法对控制方程进行求解，求解结果分别与单相弹性介质，双相饱和介质和三相非饱和介质的解析解进行对比，均验证了该求解方法的准确性；不同介质模型的计算耗时分析表明2.5维有限元法是目前求解多孔介质动力问题的一种优势算法。通过数值分析发现，采用传统的、未考虑变形效应的土-水特征曲线会低估非饱和地基的振动强度。
- 非饱和土 /
- 2.5维有限元法 /
- 动力响应 /
- 土-水特征曲线 /
- 水-力耦合
Abstract: To investigate the dynamic response of unsaturated soils under moving loads, a new soil-water characteristic curve (SWCC) model is established considering the deformation caused by applied load based on the traditional V-G SWCC model. Using this modified SWCC model, a dynamic governing equation for unsaturated soils, which fully describes the water-force coupling effects of unsaturated soils under dynamic loading, is derived. The governing equation is solved using the 2.5-dimensional finite element method (2.5D FEM). The obtained solutions are compared with the analytical ones for single-phase medium, double-phase saturated medium and three-phase unsaturated medium, respectively, which all confirm the accuracy of the proposed solution method. The computational time analysis for different medium models demonstrate that the 2.5D FEM is an advantageous algorithm for solving the dynamic problems of porous media. The numerical analysis reveals that using the traditional SWCC without considering the deformation will underestimate the vibration intensity of unsaturated foundations.
- unsaturated soil /
- 2.5-dimensional finite element method /
- dynamic response /
- soil-water characteristic curve /
- water-force coupling effect

HTML全文

0. 引言

近年来，全球区域性暴雨频发，降雨量远高于以往水平^[1-2]，对铁路线路的安全运营带来挑战。1965年，日本新干线就曾因连续降雨等原因导致线路沉降过大，不得不将列车运行速度由最初设计的210 km/h下降至110~180 km/h^[3]。1997年3月份，哥伦比亚Conrad地区的铁路线路，由于强降雨导致路基和地基土体含水率迅速增大，在列车仅以平均车速43.45 km/h通过时，线路突然坍塌，造成重大安全事故^[4]。对于土体，水的存在既减小了土颗粒之间的摩擦力，影响线路的运营。由此可见，实际铁路线路的稳定性受移动荷载与含水路基/地基形成的水-动力耦合作用的影响。由于自然界中土体的含水率处于变化中，其饱和度也随之改变；因此，将路基/地基视为土颗粒、水、气组成的三相非饱和土，研究不同饱和度下，移动荷载作用时的动力响应可以对运行线路的稳定性进行更准确的预判。

移动荷载作用下的地基动力响应问题最初是将地基视为单相弹性介质^[5-7]。近年来，国内外许多学者基于Biot饱和多孔介质波传播理论^[8]，将单相介质扩展到两相介质，建立了将地基视为饱和两相介质的动力求解方法和模型。Lu等^[9]建立了三维解析解模型，考虑了饱和多孔介质及移动点荷载，研究了动力响应与土体剪切波波速的关系。Cai等^[10]采用双重傅立叶变换和逆变换，考虑了流固耦合作用，研究了荷载加载方式和饱和地基土参数对动力响应的影响。Gao等^[11]、Bian等^[12]和胡静等^[13]基于2.5维有限元法研究了饱和土渗透系数和列车速度对地基动力响应的影响。考虑到实际土体的三相性，徐明江^[14]基于连续介质力学理论，结合土-水特征曲线（soil-water characteristic curve，简称SWCC）以及Mualem^[15]理论，推导了移动荷载作用下的非饱和土波动方程。Lu等^[16]、Fang等^[17]采用积分变换的方式求解了该方程。同样基于该方程，Gao等^[18]、李绍毅^[19]研究了列车荷载作用下的非饱和地基振动。

在非饱和土动力响应的研究中^{[14, 16-19]}，需要基于SWCC推导不同饱和度下的孔隙流体渗透性和土的强度，因而SWCC对土体的动力响应具有决定性作用。SWCC描述的是非饱和土土体含水率和基质吸力之间的关系，早期的SWCC研究中，通常假设基质吸力和含水率之间存在着唯一关系^[20-22]。实际上，除了含水率，变形也会使土体饱和度发生改变。非饱和土在荷载作用下，土体必然会发生相应的变形，因而，将变形效应考虑进SWCC模型中更为合理。遗憾的是，已有关于非饱和土动力响应的研究均采用早期的V-G模型^[21]来描述土-水特征曲线，未考虑变形对SWCC的影响，不能完整的描述荷载作用下的非饱和土持水特性。

综上所述，目前还未出现基于考虑变形效应的SWCC的非饱和土动力响应求解方法。为此，本文在现有研究的基础上，建立了考虑变形效应的SWCC模型，并据此推导了考虑SWCC变形效应的非饱和土波动方程，利用2.5有限元法对该方程进行求解；求解结果与现有文献中的解析解进行了验证；最后对该方法的计算效率进行了对比，并对引入变形效应的影响进行了初步分析。

1. 理论求解

1.1 考虑变形的土-水特征曲线

孙德安^[23]的研究表明，非饱和土的土-水特征曲线与应力历史和应力状态无直接关系，在土体受力变形过程中，可以用孔隙比来表示孔隙的变化。蔡国庆等^[24]、胡冉等^[25]及张雪东等^[26]同样也是采用孔隙比来表征变形。图 1整理了不同初始孔隙比的土在等吸力下的等向压缩试验结果；其中，孙德安^[23]对珍珠黏土进行了基质吸力为147 kPa的等向压缩试验；蔡国庆等^[24]给出了膨润土/高岭土混合土体在基质吸力为100 kPa的等向压缩试验结果。图 1中s为基质吸力，e为孔隙比，S_r为孔隙水饱和度，简称饱和度。

图 1 饱和度与孔隙比的关系图

Figure 1. Relationship between saturation and void ratio

下载: 全尺寸图片幻灯片

从中可以发现，吸力一定时，土体受外力导致的孔隙比减小会使饱和度增加。不同初始孔隙比条件下，等吸力时的孔隙比与饱和度近似呈线性关系且直线斜率相近。可见，考虑变形效应的SWCC可近似认为是将原有SWCC进行偏移得到。假定中的直线斜率为 ${\lambda _{{\text{se}}}}$ ，结合蔡国庆等^[24]关于孔隙比与SWCC的研究，不同孔隙比下的SWCC可总结为图 2。

图 2 不同孔隙比下的土-水特征曲线示意图

Figure 2. Schematic diagram of soil-water characteristic curve under different void ratios

下载: 全尺寸图片幻灯片

根据图 2，考虑变形效应的SWCC模型可以表示为

${S_{\text{r}}} = {f_{{S_{\text{r}}}}}(e) + {f_{{S_{\text{r}}}}}(s) \text{，}$

(1)

${f_{{S_{\text{r}}}}}(e) = {\lambda _{{\text{se}}}}(e - {e_0}) 。$

(2)

式中： ${e_0}$ 为初始孔隙比； ${f_{{S_{\text{r}}}}}(e)$ 表示饱和度与孔隙比之间的关系； ${f_{{S_{\text{r}}}}}(s)$ 为 ${e_0}$ 时的SWCC。

在众多SWCC模型中，V-G模型适用范围广，并且用于求解非饱和动力问题时，能够形成关于时间的一阶常微分方程，实现求解，因此现有关于非饱和土的动力求解均是基于V-G模型^{[14, 16-19]}。当 ${f_{{S_{\text{r}}}}}(s)$ 为V-G模型，即 $\frac{{S_{\text{r}}^0 - {S_{{\text{w0}}}}}}{{(1 - {S_{{\text{w0}}}})}} = {\text{ }}{\left\{ {1 + {{(\beta s)}^k}} \right\}^{ - m}}$ ， $S_{\text{r}}^0$ 为 ${e_0}$ 时V-G模型中的饱和度， ${S_{{\text{w0}}}}$ 为残余饱和度，s为基质吸力， $\beta$ ，m，k为V-G模型的参数，其中m=1-1/k。则，考虑变形效应的SWCC模型可以写为

${S_{\text{r}}} = {S_{{\text{w0}}}} + {\lambda _{{\text{se}}}}(e - {e_0}) + (1 - {S_{{\text{w0}}}}){\left\{ {1 + {{(\beta s)}^k}} \right\}^{ - m}} 。$

(3)

此时，有效饱和度 ${S_{\text{e}}}$ 为

${S_{\text{e}}} = \frac{{[S_{\text{r}}^{\text{0}} + {\lambda _{{\text{se}}}}(e - {e_0})] - {S_{{\text{w0}}}}}}{{(1 - {S_{{\text{w0}}}})}} 。$

(4)

根据Mualem理论^[15]，非饱和土中流体的渗透系数是关于有效饱和度的函数为

${k_{\text{w}}} = \frac{{{\rho _{\text{w}}}g\kappa }}{{{\eta _{\text{w}}}}}\sqrt {{S_{\text{e}}}} {\left\{ {1 - {{\left[ {1 - {{({S_\text{e}})}^{\frac{1}{m}}}} \right]}^m}} \right\}^2} \text{，}$

(5)

${k_{\text{a}}} = \frac{{{\rho _{\text{a}}}g\kappa }}{{{\eta _{\text{a}}}}}\sqrt {1 - {S_{\text{e}}}} {\left[ {1 - {{({S_{\text{e}}})}^{\frac{1}{m}}}} \right]^{2m}} 。$

(6)

式中： ${k_{\text{w}}}$ ， ${k_{\text{a}}}$ 分别为孔隙水和孔隙气的达西渗透系数； ${\rho _{\text{w}}}$ ， ${\rho _{\text{a}}}$ 为孔隙水和孔隙气的密度； ${\eta _{\text{w}}}$ ， ${\eta _{\text{a}}}$ 分别为孔隙水和孔隙气的黏滞系数（Pa·s）；g为重力加速度； $\kappa$ 为固有渗透系数（m²）。

当考虑变形的影响时， ${S_{\text{e}}}$ 的表达式发生变化，流体渗透系数也发生相应的改变。

与此同时，徐明江^[14]根据试验结果和经验公式，建立了非饱和土的动剪切模量与有效饱和度之间的联系：

$G = {G_{\text{s}}} + 2050\tan \varphi '\int\limits_0^s {{S_{\text{e}}}} {\text{d}}s 。$

(7)

式中： $G$ 为非饱和土的动剪切模量； ${G_{\text{s}}}$ 为土体饱和状态下的动剪切模量， $\varphi '$ 为土体的有效内摩擦角。

1.2 非饱和多孔介质动力控制方程

固、液、气三相的质量守恒方程分别为

$\frac{{(1 - n)\partial {\rho _{\text{s}}}}}{{\partial t}} - {\rho _{\text{s}}}\frac{{\partial n}}{{\partial t}} + {\rho _{\text{s}}}(1 - n)\nabla \dot u = 0 \text{，}$

(8)

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

(9)

$\frac{{{S_{\text{a}}}{\rho _{\text{a}}}\partial n}}{{\partial t}} + \frac{{{\rho _{\text{a}}}n\partial {S_{\text{a}}}}}{{\partial t}} + \frac{{n{S_{\text{a}}}\partial {\rho _{\text{a}}}}}{{\partial t}} + {\rho _{\text{a}}}n{S_{\text{a}}}\nabla {\dot u^{\text{a}}} = 0 。$

(10)

式中： $n$ 为孔隙率；t为时间； ${\rho _{\text{s}}}$ 为土骨粒密度； $u$ ， ${u^{\text{w}}}$ ， ${u^{\text{a}}}$ 分别为土骨架位移，孔隙水位移，孔隙气位移，上标 $\cdot$ 表示对时间的导数； $\nabla$ 为汉密尔顿算子； ${S_{\text{a}}}$ 为孔隙气饱和度， ${S_{\text{a}}} = 1 - {S_{\text{r}}}$ 。

固、液、气三相应力引起土单元压缩变形的本构方程可以表示为

$\frac{{\partial {\rho _{\text{s}}}}}{{{\rho _{\text{s}}}\partial t}} = - \frac{{\partial \sigma _{ij}^{\text{s}}}}{{3{K_{\text{s}}}{\text{d}}t}} \text{，}$

(11)

$\frac{{\partial {\rho _{\text{w}}}}}{{{\rho _{\text{w}}}\partial t}} = \frac{{\partial {P^{\text{w}}}}}{{{K_{\text{w}}}{\text{d}}t}} \text{，}$

(12)

$\frac{{\partial {\rho _{\text{a}}}}}{{{\rho _{\text{a}}}\partial t}} = \frac{{\partial {P^{\text{a}}}}}{{{K_{\text{a}}}{\text{d}}t}} 。$

(13)

式中： $\sigma _{ij}^{\text{s}}$ 为作用在土颗粒上的应力； ${K_{\text{s}}}$ ， ${K_{\text{w}}}$ ， ${K_{\text{a}}}$ 分别为土颗粒、孔隙水、孔隙气的压缩模量； ${P^{\text{w}}}$ 为孔隙水压， ${P^{\text{a}}}$ 为孔隙气压。

考虑土颗粒由于孔隙流体压力作用产生的变形时，由弹性理论，有效应力分量可表示为

${\sigma '_{ij}} = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} + \frac{{{K_{\text{b}}}}}{{{K_{\text{s}}}}}P{\delta _{ij}} 。$

(14)

式中： $\lambda$ ， $\mu$ 为Lame常数； $\Theta$ 为土骨架的体积应变； ${\varepsilon _{ij}}$ 为土颗粒的应变； ${\delta _{ij}}$ 为克罗内克符号； ${K_{\text{b}}}$ 为土骨架的压缩模量且 ${K_{\text{b}}} = \lambda + \frac{{2\mu }}{3}$ ；P为孔隙流体压力，包括孔隙水压及孔隙气压， $P = {S_{\text{r}}}{P^{\text{w}}} + {S_{\text{a}}}{P^{\text{a}}}$ 。

根据应力空间平均化方法，土单元的总应力可以表示为

${\sigma _{ij}} = (1 - n)\sigma _{ij}^{\text{s}} - n{S_{\text{r}}}{P^{\text{w}}}{\delta _{ij}} - n{S_{\text{a}}}{P^{\text{a}}}{\delta _{ij}} 。$

(15)

根据Bishop单一变量非饱和土有效应力原理：

${\sigma '_{ij}} = {\sigma _{ij}} + {\delta _{ij}}P \text{，}$

(16)

将式（16）代入式（14），可以得到

${\sigma _{ij}} + {\delta _{ij}}P = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} + P{\delta _{ij}} 。$

(17)

将式（15）代入式（17），移项后可以得到

$(1 - n)\sigma _{ij}^{\text{s}} = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} - (\alpha - n)P{\delta _{ij}} 。$

(18)

式中， $\alpha = 1 - \frac{{{K_{\text{b}}}}}{{{K_{\text{s}}}}}$ 。

将式（11）代入式（18）并展开，可以得到

$\begin{aligned} &-\frac{\partial \sigma_{i j}^{\mathrm{s}}}{3 K_{\mathrm{s}} \mathrm{~d} t}=\frac{1}{3 K_{\mathrm{s}}} \frac{1}{n-1}\left[2 \mu \frac{\partial \varepsilon_{i j}}{\partial t}+\lambda \frac{\partial \Theta}{\partial t} \delta_{i j}-\right. \\ &\left.(\alpha-n) S_{\mathrm{r}} \frac{\partial P^{\mathrm{w}}}{\partial t}-(\alpha-n) S_{\mathrm{a}} \frac{\partial P^{\mathrm{a}}}{\partial t}\right] \end{aligned}$

(19)

将式（19）代入式（8）中，最终可以求解出孔隙率的偏导数：

$\frac{{\partial n}}{{\partial t}} = (\alpha - n)\nabla \dot u + \frac{{(\alpha - n){S_{\text{r}}}}}{{{K_{\text{s}}}}}\frac{{\partial {P^{\text{w}}}}}{{\partial t}} + \frac{{(\alpha - n){S_{\text{a}}}}}{{{K_{\text{s}}}}}\frac{{\partial {P^{\text{a}}}}}{{\partial t}}。$

(20)

结合式（9），（10），（20），可得

$\begin{aligned} & (\alpha-n) \nabla \dot{u}+n S_{\mathrm{r}} \nabla \dot{u}^{\mathrm{w}}+n S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{a}}+ \\ & {\left[\frac{(\alpha-n) S_{\mathrm{r}}}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{r}}}{K_{\mathrm{w}}}\right] \frac{\partial P^{\mathrm{w}}}{\partial t}+\left[\frac{(\alpha-n) S_{\mathrm{a}}}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{a}}}{K_{\mathrm{a}}}\right] \frac{\partial P^{\mathrm{a}}}{\partial t}=0 。} \end{aligned}$

(21)

从式（3）推导出饱和度对时间的导数：

$\begin{aligned} \frac{\partial S_{\mathrm{r}}}{\partial t}= & \lambda_{\mathrm{se}} \frac{\partial e}{\partial t}-\beta k m\left(1-S_{\mathrm{w} 0}\right) \\ & \left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{\frac{m+1}{m}}\left[\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{-\frac{1}{m}}-1\right]^{\frac{k-1}{k}} \frac{\partial s}{\partial t}。 \end{aligned}$

(22)

根据定义，孔隙比 $e = \frac{n}{{1 - n}}$ ，基质吸力 $s = {P^{\text{a}}} - {P^{\text{w}}}$ ，则式（22）可以改写为

$\begin{aligned} \frac{\partial S_{\mathrm{r}}}{\partial t}= & \lambda_{\mathrm{se}} \frac{1}{(1-n)^2} \frac{\partial n}{\partial t}-\beta k m\left(1-S_{\mathrm{w} 0}\right)\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{\frac{m+1}{m}}\cdot \\ & {\left[\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{-\frac{1}{m}}-1\right]^{\frac{k-1}{k}}\left(\frac{\partial P^{\mathrm{a}}}{\partial t}-\frac{\partial P^{\mathrm{w}}}{\partial t}\right) 。} \end{aligned}$

(23)

令 ${A}_{\text{e}}={\lambda }_{\text{se}}\frac{1}{{(1-n)}^{2}}，{A}_{\text{s}}=-\beta km(1-{S}_{\text{w0}}){\left(\frac{{S}_{\text{r}}^{\text{0}}-{S}_{\text{w0}}}{1-{S}_{\text{w0}}}\right)}^{\frac{m+1}{m}}\cdot$ ${\left[ {{{\left( {\frac{{S_{\text{r}}^{\text{0}} - {S_{{\text{w0}}}}}}{{1 - {S_{{\text{w0}}}}}}} \right)}^{ - \frac{1}{m}}} - 1} \right]^{\frac{{k - 1}}{k}}}$ ；结合式（9），（10），（20）和（23），可得

$\begin{aligned} & A_{\mathrm{e}}(\alpha-n) \nabla \dot{u}+S_{\mathrm{r}} S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{w}}-S_{\mathrm{r}} S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{a}}+ \\ & {\left[\frac{A_{\mathrm{e}}(\alpha-n) S_{\mathrm{r}}}{K_{\mathrm{s}}}-A_{\mathrm{s}}+\frac{S_{\mathrm{r}} S_{\mathrm{a}}}{K_{\mathrm{w}}}\right] \frac{\partial P^{\mathrm{w}}}{\partial t}+} \\ & \quad\left[\frac{A_{\mathrm{e}}(\alpha-n) S_{\mathrm{a}}}{K_{\mathrm{s}}}+A_{\mathrm{s}}-\frac{S_{\mathrm{r}} S_{\mathrm{a}}}{K_{\mathrm{a}}}\right] \frac{\partial P^{\mathrm{a}}}{\partial t}=0。 \end{aligned}$

(24)

式（21），（24）组成求解 $\frac{\partial {P}^{\text{w}}}{\partial t}和\frac{\partial {P}^{\text{a}}}{\partial t}$ 的方程组：

$\left.\begin{array}{l} B_{11} \frac{\partial P^{\mathrm{w}}}{\partial t}+B_{12} \frac{\partial P^{\mathrm{a}}}{\partial t}+B_{13} \nabla \dot{u}+B_{14} \nabla \dot{u}^{\mathrm{w}}+B_{15} \nabla \dot{u}^{\mathrm{a}}=0, \\ B_{21} \frac{\partial P^{\mathrm{w}}}{\partial t}+B_{22} \frac{\partial P^{\mathrm{a}}}{\partial t}+B_{23} \nabla \dot{u}+B_{24} \nabla \dot{u}^{\mathrm{w}}+B_{25} \nabla \dot{u}^{\mathrm{a}}=0。 \end{array}\right\}$

(25)

式中： ${B_{11}} = \frac{{(\alpha - n){S_{\text{r}}}}}{{{K_{\text{s}}}}} + \frac{{n{S_{\text{r}}}}}{{{K_{\text{w}}}}}$ ， ${B_{12}} = \frac{{(\alpha - n){S_{\text{a}}}}}{{{K_{\text{s}}}}} + \frac{{n{S_{\text{a}}}}}{{{K_{\text{a}}}}}$ ， ${B_{13}} = \alpha - n$ ， ${B_{14}} = n{S_{\text{r}}}$ ， ${B_{15}} = n{S_{\text{a}}}$ ； ${B_{21}} =$ $\frac{{A}_{\text{e}}(\alpha -n){S}_{\text{r}}}{{K}_{\text{s}}}-{A}_{\text{s}}+\frac{{S}_{\text{r}}{S}_{\text{a}}}{{K}_{\text{w}}}，{B}_{22}=\frac{{A}_{\text{e}}(\alpha -n){S}_{\text{a}}}{{K}_{\text{s}}}+{A}_{\text{s}}-\frac{{S}_{\text{r}}{S}_{\text{a}}}{{K}_{\text{a}}}，$ ${B_{23}} = {A_{\text{e}}}(\alpha - n)$ ， ${B_{24}} = - {B_{25}} = {S_{\text{r}}}{S_{\text{a}}}$ 。当不考虑变形对SWCC的影响，即A_e为0时，方程（25）退化后的结果与Lu等^[16]和徐明江^[14]的结果一致。

令 $M=\frac{\alpha -n}{{K}_{\text{s}}}，{M}_{\text{w}}=\frac{n}{{K}_{\text{w}}}，{M}_{\text{a}}=\frac{n}{{K}_{\text{a}}}，{C}_{1}=\frac{{A}_{\text{e}}n+{S}_{\text{r}}}{n}，$ ${C_2} = \frac{{{A_{\text{e}}}n - {S_{\text{a}}}}}{n}$ ；采用相对位移， ${w}_{i}=n{S}_{\text{r}}({u}_{i}^{\text{w}}-{u}_{i}) ，$ ${v_i} = n{S_{\text{a}}}(u_i^{\text{a}} - {u_i})$ ，可求得

$- \frac{{\partial {P^{\text{w}}}}}{{\partial t}} = {D_{11}}\nabla \dot u + {D_{12}}\nabla \dot w + {D_{13}}\nabla \dot v \text{，}$

(26)

$- \frac{{\partial {P^a}}}{{\partial t}} = {D_{21}}\nabla \dot u + {D_{22}}\nabla \dot w + {D_{23}}\nabla \dot v 。$

(27)

式中： ${w_i}$ 为液相与土骨架的相对位移； ${v_i}$ 为气相与土骨架的相对位移，D₁₁，D₁₂，D₁₃，D₂₁，D₂₂，D₂₃为方程系数，具体表达式见附录1。

根据动量守恒原理，土体的运动方程可以表示为

${\sigma _{ij,j}} = {\rho _{\text{s}}}{\ddot u_i} + {\rho _{\text{w}}}{\ddot u_i}^{\text{w}} + {\rho _{\text{a}}}{\ddot u_i}^{\text{a}} 。$

(28)

式中：上标 $\cdot \cdot$ 表示对时间的二阶导数。

根据广义达西定律，孔隙水和孔隙气的渗流运动方程可以表示为

$n{S_{\text{r}}}(\dot u_i^{\text{w}} - {\dot u_i}) = \frac{{{k_{\text{w}}}}}{{{\rho _{\text{w}}}g}}( - P_{,i}^{\text{w}} - {\rho _{\text{w}}}\ddot u_i^{\text{w}}) \text{，}$

(29)

$n{S_{\text{a}}}(\dot u_i^{\text{a}} - {\dot u_i}) = \frac{{{k_a}}}{{{\rho _{\text{a}}}g}}( - P_{,i}^{\text{a}} - {\rho _{\text{a}}}\ddot u_i^{\text{a}}) 。$

(30)

将式（17）代入式（28）~（30），采用式（26）中的相对位移，可以得到三相介质运动方程：

$(\lambda + \mu )\nabla (\nabla u) + \mu {\nabla ^2}u - {\delta _{ij}}\alpha P = \rho \ddot u + {\rho _{\text{w}}}\ddot w + {\rho _{\text{a}}}\ddot v \text{，}$

(31)

$- P,_i^{\text{w}} = {\rho _{\text{w}}}{\ddot u_i} + \frac{{{\rho _{\text{w}}}}}{{n{S_{\text{r}}}}}{\ddot w_i} + \frac{{{\rho _{\text{w}}}g}}{{{k_{\text{w}}}}}{\dot w_i} \text{，}$

(32)

$- P,_i^{\text{a}} = {\rho _{\text{a}}}{\ddot u_i} + \frac{{{\rho _{\text{a}}}}}{{n{S_{\text{a}}}}}{\ddot v_i} + \frac{{{\rho _{\text{a}}}g}}{{{k_{\text{a}}}}}{\dot v_i} 。$

(33)

式中： $\rho$ 为非饱和土的总密度， $\rho = (1 - n){\rho _{\text{s}}} + n{S_r}{\rho _{\text{w}}} +$ $n{S_{\text{a}}}{\rho _{\text{a}}}$ 。

将式（26），（27）代入式（31）~（33），得到最终考虑变形效应的非饱和土动力控制方程为

$\begin{aligned} & \mu u_{i, j j}+\left(\lambda+\mu+\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) u_{j, j i}+ \\ & \left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_a D_{22}\right) w_{j, j i}+ \\ & \left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) v_{j, j i}=\rho \ddot{u}+\rho_{\mathrm{w}} \ddot{w}+\rho_{\mathrm{a}} \ddot{v}, \end{aligned}$

(34)

${D_{11}}{u_{j,ji}} + {D_{12}}{w_{j,ji}} + {D_{13}}{v_{j,ji}} = {\rho _{\text{w}}}{\ddot u_i} + \frac{{{\rho _{\text{w}}}}}{{n{S_{\text{r}}}}}{\ddot w_i} + \frac{{{\rho _{\text{w}}}g}}{{{k_{\text{w}}}}}{\dot w_i} \text{，}$

(35)

$\begin{array}{*{20}{c}} {{D_{21}}{u_{j,ji}} + {D_{22}}{w_{j,ji}} + {D_{23}}{v_{j,ji}} = {\rho _{\text{a}}}{{\ddot u}_i} + \frac{{{\rho _a}}}{{n{S_{\text{a}}}}}{{\ddot v}_i} + \frac{{{\rho _{\text{a}}}g}}{{{k_{\text{a}}}}}{{\dot v}_i}} \end{array}。$

(36)

1.3 2.5维有限元求解

2.5维有限元法是由Yang等^[27]用于求解移动荷载作用下完全弹性半空间的响应。该方法的特点在于需要假定研究对象在沿荷载移动方向的几何形状以及材料属性都是连续一致的，通过傅里叶变换将荷载移动方向的空间坐标（假定为x方向）变换成波数，只在频域-波数域内求解出y-z二维平面的响应，再通过傅里叶逆变换就可以得到x方向的响应，从而获得三维空间响应的解答。由于轨道结构大都具备沿列车行进方向连续一致的特性，因此，2.5维有限元法近年来在交通岩土领域得到了广泛的应用^[11-13]。

在应用时，首先定义关于荷载移动方向x以及时间t的傅里叶变换：

$\tilde{\bar{u}}\left(\xi_x, y, z, \omega\right)=\int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} u(x, y, z, t) \mathrm{e}^{\mathrm{i} \xi_x x} \mathrm{e}^{-\mathrm{i} \omega t} \mathrm{~d} x \mathrm{~d} t。$

(37)

相应的傅里叶逆变换为

$u(x, y, z, t)=\frac{1}{4 \pi^2} \int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} \tilde{\bar{u}}\left(\xi_x, y, z, \omega\right) \mathrm{e}^{\mathrm{i} \xi_{x^x}} \mathrm{e}^{\mathrm{i} \omega t} \mathrm{~d} \xi_x \mathrm{~d} \omega。$

(38)

式中： ${\xi _x}$ 为沿x方向的波数（m^-1）； $\omega$ 为圆角频率（rad/s）；上标‘–’和‘~’分别为波数域和频域内的量。

对式（34）~（36）进行傅里叶变换后，非饱和土的动力控制方程在频域内的表达式为

$\begin{aligned} & \mu u_{i, j j}+\left(\lambda+\mu+\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) u_{j, j i}+\rho \omega^2 u_i+ \\ & \left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_{\mathrm{a}} D_{22}\right) w_{j, j i}+\rho_{\mathrm{w}} \omega^2 w_i+ \\ & \quad\left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) v_{j, j i}+\rho_{\mathrm{a}} \omega^2 v_i=0, \end{aligned}$

(39)

$\begin{gathered} D_{11} u_{j, j i}+\rho_{\mathrm{w}} \omega^2 u_i+D_{12} w_{j, j i}+\frac{\rho_{\mathrm{w}}}{n S_{\mathrm{r}}} \omega^2 w_i-\frac{\rho_{\mathrm{w}} g}{k_{\mathrm{w}}} \mathrm{i} \omega w_i+ \\ D_{13} v_{j, j i}=0 \quad, \end{gathered}$

(40)

$\begin{gathered} D_{21} u_{j, j i}+\rho_{\mathrm{a}} \omega^2 u_i+D_{22} w_{j, j i}+\frac{\rho_{\mathrm{a}}}{n S_{\mathrm{a}}} \omega^2 v_i- \\ \frac{\rho_{\mathrm{a}} g}{k_{\mathrm{a}}} \mathrm{i} \omega v_i+D_{23} v_{j, j i}=0 \end{gathered}$

(41)

经过波数变换后整理得到的矩阵形式：

$\left[\begin{array}{ccc} \boldsymbol{K}_\boldsymbol{1}+\boldsymbol{K}_\boldsymbol{2}-\boldsymbol{M}_\boldsymbol{1} & \boldsymbol{L}_\boldsymbol{1}-\boldsymbol{M}_\boldsymbol{2} & \boldsymbol{G}_\boldsymbol{1}-\boldsymbol{M}_\boldsymbol{3} \\ \boldsymbol{K}_\boldsymbol{3}-\boldsymbol{M}_\boldsymbol{2} & \boldsymbol{L}_\boldsymbol{2}-\boldsymbol{M}_\boldsymbol{4} & \boldsymbol{G}_\boldsymbol{2} \\ \boldsymbol{K}_\boldsymbol{4}-\boldsymbol{M}_\boldsymbol{3} & \boldsymbol{L}_\boldsymbol{3} & \boldsymbol{G}_\boldsymbol{3}-\boldsymbol{M}_\boldsymbol{5} \end{array}\right]\left[\begin{array}{l} \boldsymbol{U} \\ \boldsymbol{W} \\ \boldsymbol{V} \end{array}\right]=\left[\begin{array}{l} \boldsymbol{F}^{\mathrm{\mathbf{s}}} \\ \boldsymbol{F}^{\mathrm{\mathbf{f}}} \\ \boldsymbol{F}^{\mathrm{\mathbf{g}}} \end{array}\right] 。$

(42)

式中：K，L，G分别为固相、液相和气相的刚度矩阵；M为各相的质量矩阵；U，W，V分别为固相、液相和气相的位移矩阵； $\boldsymbol{F}^\mathbf{s}，\boldsymbol{F}^\mathbf{f}，\boldsymbol{F}^\mathbf{g}$ 分别为固相、液相和气相的外力矢量矩阵。各矩阵的具体表达式见附录1。

求解时，首先建立y-z剖面的二维有限元模型，采用4节点等参单元划分模型网格。参考王勖成《有限单元法》^[28]，采用拉格朗日插值函数和高斯积分，根据式（42）得到单元刚度矩阵，最后集成总体刚度矩阵后进行求解。整个求解过程是根据式（42）在Matlab中通过编程实现。

2. 求解精度及求解效率

2.1 模型介绍

Yang^[27]讨论了网格尺寸与求解精度以及计算经济性的关系，认为计算区域2.5维有限元网格尺寸的最大值L(m)必须要小于等于 ${\lambda '_{\text{s}}}/6$ ，才能满足计算精度要求，网格尺寸的最小值大于等于0.5 ${\lambda '_{\text{s}}}$ 能兼顾计算的经济性；其中， ${{\lambda }^{\prime }}_{\text{s}}=2\text{π}/{{k}^{\prime }}_{\text{s}}，$ ${k'_{\text{s}}} = \sqrt {{k_{\text{s}}}^{\text{2}} - {{\left( {\frac{{\omega - {\omega _0}}}{c}} \right)}^2}}$ ， ${k_{\text{s}}}$ 为弹性半空间土体沿x方向上的剪切波波数， ${k'_{\text{s}}}$ 为移动荷载作用下土体沿x方向的剪切波波数， ${\lambda '_{\text{s}}}$ 为移动荷载作用下剪切波的波长，单位为m， ${\omega _0}$ 为移动荷载初始自振频率，c为外部荷载的移动速度。

本节建立了一个宽100 m，深50 m的2.5维有限元模型用以模拟非饱和半无限空间，用四边形等参单元模拟地基中的土单元，其中每个节点包含空间3个方向的位移自由度，如图 3所示，其中，由于x方向也被波数代替，因此只需在y-z平面进行网格划分。文献[16]中非饱和土的最小剪切波速为38 m/s，如外加荷载最大运行速度为200 m/s，则2.5维有限元网格的最大值不能超过0.2 m。因此，在建模时，本模型计算核心区域的网格尺寸为0.08 m×0.08 m，0.2 m×0.2 m；随后逐渐过渡至0.5 m×0.5 m，边缘网格尺寸为1.0 m×1.0 m。Bian等^[12]和胡静等^[13]在求解饱和地基动力响应时采用多层阻尼边界，很好地模拟了饱和半无限空间的动力响应。因此，本文同样采用多层阻尼边界对模型边界的波动进行吸收。整个2.5维有限元模型共计9052个网格，9250个节点。

图 3 模型网格分布示意图

Figure 3. Diagram of mesh distribution

下载: 全尺寸图片幻灯片

2.2 退化验证

（1）弹性地基验证

首先，将参数 $\alpha$ ， ${\rho _{\text{w}}}$ ， ${\rho _{\text{a}}}$ ，S_r和S_a的值都取无限小，则非饱和动力控制方程退化为单相弹性介质的Naiver波动方程。采用文献[12]中的单相弹性地基的参数（表 1所示）赋至上述2.5维有限元模型中进行计算。

表 1 单相弹性地基计算参数^[12]

Table 1. Parameters of single-phase elastic medium^[12]

K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r
35	2.25	145	25	0.125	0.001	2500	0.001	0.001	0.001

下载: 导出CSV

| 显示表格

图 4为单位荷载（1 N）以时速75 m/s移动时，观察点(x, 1, 1)处的位移响应时程曲线，x为荷载移动方向。由图 4可知，基于退化至弹性地基参数的2.5维非饱土模型计算得到的结果与半解析解结果吻合。

图 4 单相弹性地基验证

Figure 4. Verification of single-phase elastic medium

下载: 全尺寸图片幻灯片

（2）饱和地基验证

令S_r等于1，参数A_e，A_s，ρ_a和S_a的值都取无限小，非饱和土的动力控制方程退化为Biot波动方程。将文献[9]中双相饱和地基的参数（见表 2）赋至2.5维非饱和半空间模型中进行计算。

表 2 双相饱和地基计算参数^[9]

Table 2. Parameters of two-phase saturated medium^[9]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r	k_D/(m·s^-1)
8.77	11	2.25	145	2997.6	0.125	0.3	2500	1000	0.001	1	1×10^-6

下载: 导出CSV

| 显示表格

图 5为单位荷载（1 N）以121 m/s移动时观察点(x, 1, 1)处的位移及超静孔压响应。中，位移乘以 $G{a_{\text{R}}}{\text{/}}F$ ，超静孔压乘以 $a_{\text{R}}^2/F$ 进行了归一化，其中F=1 N，a_R=1.0 m，G为表 2中的剪切模量。从图 5中可以发现，基于退化至饱和地基参数的2.5维非饱土模型计算得到的位移响应和超静孔压响应都与半解析解吻合。

图 5 饱和地基验证

Figure 5. Verification of two-phase saturated medium

下载: 全尺寸图片幻灯片

2.3 与非饱和解析解验证

Lu等^[16]采用解析法求解了矩形荷载作用下非饱和半空间的动力响应，但该研究并未考虑变形对土-水特征曲线的影响。此处选取该研究中的参数进行计算，并与其结果进行对比以验证本模型对非饱和动力问题求解的准确性。模型中，荷载为0.32 m×0.32 m的矩形荷载，大小为30 kN，频率为2 Hz；地基土参数见表 3，地基上部为完全透水边界。由于非饱和土的剪切模量与饱和度有关，因而土体的剪切波波速也随着饱和度的变化而变化。为此，定义参数v_s^sat为不同饱和度下的土体剪切波波速，其上标sat代表饱和度。

表 3 非饱和地基计算参数^[16]

Table 3. Parameters of unsaturated medium^[16]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G_s/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)
8.33	35	2.25	145	3.85	0.35	0.45	2650	1000	1.28

S_r	S_a	η_w/(Ns·m^-2)	η_a/(Ns·m^-2)	S_w0	κ/m²	β/(Pa^-1)	m	k	φ'
0.5，0.8	0.5，0.2	0.001	18×10^-6	0.05	5.3×10^-13	1.0×10^-4	0.5	2	10^°

下载: 导出CSV

| 显示表格

图 6为观察点（0, 0, 0.5）在不同饱和度和不同荷载移动速度下的响应。从图 6中可以看出，在低速和高速条件下，2.5维有限元在不同饱和度下的计算结果均与Lu等^[16]的解析解吻合良好。

图 6 非饱和地基验证

Figure 6. Verification of three-phase unsaturated medium

下载: 全尺寸图片幻灯片

2.2和2.3节共进行了3种介质模型的验证，数值结果都与已有解析/半解析结果吻合一致，证明了2.5维有限元方法求解多孔介质动力问题的准确性。

2.4 计算效率

采用Intel Core i9-10900k 10核20线程，基本频率为3.7 GHz的电脑对上述2.5维有限元模型进行计算。移动荷载为8节编组的列车荷载，移动速度为75 m/s，移动距离为384 m，结果表明，仅考虑每节点3个自由度的单相弹性介质模型所需的计算时间为60.43 s，考虑每节点6个自由度的双相饱和介质模型所需的计算时间为447.54 s，考虑每节点9个自由度的三相非饱和介质模型所需的计算时间为1807.24 s。可见，随着每节点自由度数的增加，模型所需的计算时间显著增大。Hall^[29]曾用ABAQUS建立了一个65 m×30 m的三维完全弹性地基模型，计算8节编组列车荷载下的动力响应，耗时约23 h。由此可见，采用有限元法求解多孔介质的动力响应问题不仅面临着二次开发困难，还必将存在计算用时过长的问题。与此同时，非饱和土波动方程求解较为复杂，将变形对土-水特征曲线的影响引入波动方程后，需要对关键控制项进行重新推导，解析和半解析求解难度增大。如果再考虑到路基/地基中饱和度横向及纵向的不均匀分布，解析和半解析方式则较难实现求解。对比而言，2.5维有限元方法体现出了高效的计算效率，与此同时还具备良好的准确性，因此，2.5维有限元方法是现阶段求解多孔介质动力问题的一种优势方法。

3. 数值分析

为研究考虑变形效应的SWCC对非饱和土动力响应的影响，基于孙德安^[23]研究，将孔隙比-饱和度直线图的斜率 ${\lambda _{{\text{se}}}}$ 取为-0.35，代入2.3节中的模型进行计算。对比了观察点（0, 0, 0.5）处，基于考虑变形与未考虑变形两种SWCC模型计算出的位移响应结果，如图 7所示。

图 7 SWCC考虑变形对位移响应的影响

Figure 7. Displacement responses under SWCC considering effect of deformation

下载: 全尺寸图片幻灯片

从图 7中可以看出，采用考虑变形效应的土-水特征曲线计算得到的位移响应大于未考虑变形的结果。根据图 2和式（3），（7）可知，考虑变形的SWCC相较于未考虑变形的曲线向右上方移动，导致同一基质吸力所对应的饱和度增大，土体剪切模量减小；因而在相同的荷载作用下其响应值更大。

图 8给出了多个饱和度和速度下，考虑变形效应所产生的位移增幅。

图 8 不同速度下考虑变形效应所产生的位移增幅随饱和度的变化

Figure 8. Amplification of deformation and saturation considering different velocities

下载: 全尺寸图片幻灯片

由图 8可知，当荷载移动速度没有超过土体的剪切波波速时，考虑变形带来的位移响应增幅与饱和度和列车速度均成正比。采用传统的土-水特征曲线模型去求解非饱和土的动力问题将低估土体的振动强度。

4. 结论

（1）本文在传统V-G模型的基础上，根据饱和度与孔隙比的关系，建立了考虑变形效应的土-水特征曲线模型，并据此推导了新的非饱和土动力控制方程，该方程能够完整描述动力作用下非饱土的水-力耦合作用。

（2）采用2.5维有限元方法对考虑SWCC变形效应的非饱和土动力控制方程进行了求解，求解结果分别与单相弹性介质，双相饱和介质和三相非饱和介质的解析/半解析结果进行对比，均验证了该求解方法的准确性。

（3）随着多孔介质模型的复杂程度即考虑相数的增加，动力求解模型的计算耗时会显著增长；2.5维有限元方法具备良好的准确性和高效的计算效率，并且能够根据实际地层特征进行建模，是目前求解多孔介质动力问题的一种优势算法。

（4）采用考虑变形效应的土-水特征曲线计算得到的位移响应大于未考虑变形的结果，说明传统的计算方法会低估非饱和土动力响应的强度。

（5）由于非饱和三相介质动力计算需要涉及众多参数，全套参数的测定也较耗时，因此文中进行的分析是基于前人已发表的参数。后续可进行非饱和土的试验研究，测定SWCC、渗透系数及强度，获得非饱和三相介质动力计算所需的参数，更加贴近实际应用。

附录

$\begin{aligned} & D_{11}=\left[A_{\mathrm{e}} n S_{\mathrm{a}}\left(M+M_{\mathrm{a}}\right)-\alpha C_1 S_{\mathrm{a}} M_{\mathrm{a}}+\alpha A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{12}=\left[C_2 S_{\mathrm{a}} M+\left(C_2-C_{\mathrm{r}}\right) S_{\mathrm{a}} M_{\mathrm{a}}+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{13}=\left[C_1 S_{\mathrm{a}} M+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+\right. \\ & \left.M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{21}=\left[-A_{\mathrm{e}} n S_{\mathrm{r}}\left(M+M_{\mathrm{w}}\right)+\alpha C_2 S_{\mathrm{r}} M_{\mathrm{w}}+\alpha A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+\right. \\ & \left.A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{22}=\left[-C_2 S_{\mathrm{r}} M+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+\right. \\ & \left.M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{23}=\left[-C_1 S_{\mathrm{r}} M-\left(C_1-C_2\right) S_{\mathrm{r}} M_{\mathrm{w}}+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{K}_1=\sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_2=\left(\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_3=D_{11} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_4=D_{21} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \mathbf{B} \mathbf{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma。 \end{aligned}$

$\begin{aligned} & \boldsymbol{L}_1=\left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_{\mathrm{a}} D_{22}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{L}_2=D_{12} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{L}_3=D_{22} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{G}_1=\left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{G}_2=D_{13} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{G}_3=D_{23} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{M}_1=\omega^2 \rho \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_2=\omega^2 \rho_{\mathrm{w}} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_3=\omega^2 \rho_{\mathrm{a}} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_4=\left(\omega^2 \frac{\rho_{\mathrm{w}}}{n S_{\mathrm{r}}}-\mathrm{i} \omega \frac{\rho_{\mathrm{w}} g}{k_{\mathrm{w}}}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \zeta 。 \end{aligned}$

$\boldsymbol{D}=\left(\begin{array}{cccccc} \lambda+2 \mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda+2 \mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda+2 \mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{array}\right),$

$\boldsymbol{B}^{\mathbf{T}}=\left[\begin{array}{cccccc} -\mathrm{i} \xi_x & 0 & 0 & \frac{\partial}{\partial_y} & 0 & \frac{\partial}{\partial_z} \\ 0 & \frac{\partial}{\partial_y} & 0 & -\mathrm{i} \xi_x & \frac{\partial}{\partial_z} & 0 \\ 0 & 0 & \frac{\partial}{\partial_z} & 0 & \frac{\partial}{\partial_y} & -\mathrm{i} \xi_x \end{array}\right],$

$\boldsymbol{m}=\left(\begin{array}{l} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{array}\right)|\mathbf{J}|=\left|\begin{array}{ll} \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \eta} y_i^{\text {element }} & \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \xi} y_i^{\text {element }} \\ \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \eta} z_i^{\text {element }} & \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \xi} z_i^{\text {element }} \end{array}\right| 。$

式中：B为应变矩阵；D为弹性矩阵； $N$ 为插值函数矩阵； ${\mathbf{J}}$ 为Jocobi矩阵； $\eta ,\zeta$ 表示局部坐标系中的变量；y_i，z_i为总体坐标系中的坐标值。

图 1 饱和度与孔隙比的关系图

Figure 1. Relationship between saturation and void ratio

下载: 全尺寸图片幻灯片

图 2 不同孔隙比下的土-水特征曲线示意图

Figure 2. Schematic diagram of soil-water characteristic curve under different void ratios

下载: 全尺寸图片幻灯片

图 3 模型网格分布示意图

Figure 3. Diagram of mesh distribution

下载: 全尺寸图片幻灯片

图 4 单相弹性地基验证

Figure 4. Verification of single-phase elastic medium

下载: 全尺寸图片幻灯片

图 5 饱和地基验证

Figure 5. Verification of two-phase saturated medium

下载: 全尺寸图片幻灯片

图 6 非饱和地基验证

Figure 6. Verification of three-phase unsaturated medium

下载: 全尺寸图片幻灯片

图 7 SWCC考虑变形对位移响应的影响

Figure 7. Displacement responses under SWCC considering effect of deformation

下载: 全尺寸图片幻灯片

图 8 不同速度下考虑变形效应所产生的位移增幅随饱和度的变化

Figure 8. Amplification of deformation and saturation considering different velocities

下载: 全尺寸图片幻灯片

表 1 单相弹性地基计算参数^[12]

Table 1 Parameters of single-phase elastic medium^[12]

K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r
35	2.25	145	25	0.125	0.001	2500	0.001	0.001	0.001

下载: 导出CSV

表 2 双相饱和地基计算参数^[9]

Table 2 Parameters of two-phase saturated medium^[9]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r	k_D/(m·s^-1)
8.77	11	2.25	145	2997.6	0.125	0.3	2500	1000	0.001	1	1×10^-6

下载: 导出CSV

表 3 非饱和地基计算参数^[16]

Table 3 Parameters of unsaturated medium^[16]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G_s/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)
8.33	35	2.25	145	3.85	0.35	0.45	2650	1000	1.28

S_r	S_a	η_w/(Ns·m^-2)	η_a/(Ns·m^-2)	S_w0	κ/m²	β/(Pa^-1)	m	k	φ'
0.5，0.8	0.5，0.2	0.001	18×10^-6	0.05	5.3×10^-13	1.0×10^-4	0.5	2	10^°

下载: 导出CSV

参考文献(29)

[1]	黄长生, 周耘, 张胜男, 等. 长江流域地下水资源特征与开发利用现状[J]. 中国地质, 2021, 48(4): 979-1000. HUANG Changsheng, ZHOU Yun, ZHANG Shengnan, et al. Groundwater resources in the Yangtze River Basin and its current development and utilization[J]. Geology in China, 2021, 48(4): 979-1000. (in Chinese)
[2]	吴庆华, 汪啸, 范越. 长江中下游地下水资源战略储备选址适宜性评价指标体系[J]. 长江科学院院报, 2022, 39(8): 145-151, 158. WU Qinghua, WANG Xiao, FAN Yue. Evaluation index system for the suitability of groundwater strategic reserve site in middle and Lower Yangtze River[J]. Journal of Yangtze River Scientific Research Institute, 2022, 39(8): 145-151, 158. (in Chinese)
[3]	池田俊雄. 地盤地質と鉄道土木の 50年[J]. 応用地質, 1998, 39(1): 4-10. IKEDA T. 50 years of geology and railway engineering [J]. Applied Geology, 1998, 39(1): 4-10. (in Japanese)
[4]	Conrad, British Columbia. Canadian National Train No. Q-102-51-26 Derailment[R]. Canada: Transportation Safety Board of Canada, 1997: 1-24.
[5]	EASON G. The stresses produced in a semi-infinite solid by a moving surface force[J]. International Journal of Engineering Science, 1965, 2(6): 581-609. doi: 10.1016/0020-7225(65)90038-8
[6]	SHENG X, JONES C J C, PETYT M. Ground vibration generated by a load moving along a railway track[J]. Journal of Sound and Vibration, 1999, 228(1): 129-156. doi: 10.1006/jsvi.1999.2406
[7]	雷晓燕, 徐斌, 徐满清. 半无限弹性空间中移动荷载动力响应的频域-波数域比例边界有限元法分析[J]. 振动工程学报, 2017, 30(5): 798-805. LEI Xiaoyan, XU Bin, XU Manqing. Using the frequency-wave domain scaled boundary finite element method for the dynamic response of the elastic half space due to moving loads[J]. Journal of Vibration Engineering, 2017, 30(5): 798-805. (in Chinese)
[8]	BIOT M A. Mechanics of deformation and acoustic propagation in porous media[J]. Journal of Applied Physics, 1962, 33(4): 1482-1498. doi: 10.1063/1.1728759
[9]	LU J F, JENG D S. A half-space saturated poro-elastic medium subjected to a moving point load[J]. International Journal of Solids and Structures, 2007, 44(2): 573-586. doi: 10.1016/j.ijsolstr.2006.05.020
[10]	CAI Y Q, SUN H L, XU C J. Response of railway track system on poroelastic half-space soil medium subjected to a moving train load[J]. International Journal of Solids and Structures, 2008, 45(18/19): 5015-5034. http://core.ac.uk/download/pdf/82345375.pdf
[11]	GAO G Y, CHEN Q S, HE J F, et al. Investigation of ground vibration due to trains moving on saturated multi-layered ground by 2.5D finite element method[J]. Soil Dynamics and Earthquake Engineering, 2012, 40: 87-98. doi: 10.1016/j.soildyn.2011.12.003
[12]	BIAN X C, HU J, THOMPSON D, et al. Pore pressure generation in a poro-elastic soil under moving train loads[J]. Soil Dynamics and Earthquake Engineering, 2019, 125: 105711. doi: 10.1016/j.soildyn.2019.105711
[13]	胡静, 唐跃, 张家康, 等. 高速列车荷载作用下饱和软土地基动力响应研究[J]. 岩土力学, 2021, 42(11): 3169-3181. HU Jing, TANG Yue, ZHANG Jiakang, et al. Dynamic responses of saturated soft soil foundation under high speed train[J]. Rock and Soil Mechanics, 2021, 42(11): 3169-3181. (in Chinese)
[14]	徐明江. 非饱和土地基与基础的动力响应研究[D]. 广州: 华南理工大学, 2010. XU Mingjiang. Study on Dynamic Response of Unsaturated Soil Foundation and Foundation[D]. Guangzhou: South China University of Technology, 2010. (in Chinese)
[15]	MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research, 1976, 12(3): 513-522. doi: 10.1029/WR012i003p00513
[16]	LU Z, FANG R, YAO H L, et al. Dynamic responses of unsaturated half-space soil to a moving harmonic rectangular load[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42(9): 1057-1077. doi: 10.1002/nag.2780
[17]	FANG R, LU Z, YAO H L, et al. Study on dynamic responses of unsaturated railway subgrade subjected to moving train load[J]. Soil Dynamics and Earthquake Engineering, 2018, 115: 319-323. doi: 10.1016/j.soildyn.2018.08.037
[18]	GAO G Y, YAO S F, YANG J, et al. Investigating ground vibration induced by moving train loads on unsaturated ground using 2.5D FEM[J]. Soil Dynamics and Earthquake Engineering, 2019, 124: 72-85. doi: 10.1016/j.soildyn.2019.05.034
[19]	李绍毅. 土体饱和度对移动荷载引起多层非饱和铁路地基振动的影响[J]. 岩土力学, 2021, 42(1): 151-159, 167. LI Shaoyi. Influences of soil saturation on the vibration of multi-layer unsaturated railway ground induced by moving load[J]. Rock and Soil Mechanics, 2021, 42(1): 151-159, 167. (in Chinese)
[20]	Corey R H Brooks and A T. Hydraulic properties of porous media and their relation to drainage design[J]. Transactions of the ASAE, 1964, 7(1): 26-28. doi: 10.13031/2013.40684
[21]	VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated Soils1[J]. Soil Science Society of America Journal, 1980, 44(5): 892. doi: 10.2136/sssaj1980.03615995004400050002x
[22]	FREDLUND D G, XING A Q. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(4): 521-532. doi: 10.1139/t94-061
[23]	孙德安. 非饱和土力学特性及本构模型[J]. 岩土工程学报, 2023, 45(1): 1-23. doi: 10.11779/CJGE20221450 SUN De'an. Mechanical behaviors and constitutive model for unsaturated soils[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(1): 1-23. (in Chinese) doi: 10.11779/CJGE20221450
[24]	蔡国庆, 田京京, 李舰, 等. 考虑变形及滞回效应影响的三维土-水特征曲面模型[J]. 土木工程学报, 2019, 52(11): 97-107. CAI Guoqing, TIAN Jingjing, LI Jian, et al. A three-dimensional soil water characteristic surface model considering deformation and hysteresis effect[J]. China Civil Engineering Journal, 2019, 52(11): 97-107. (in Chinese)
[25]	胡冉, 陈益峰, 周创兵. 考虑变形效应的非饱和土相对渗透系数模型[J]. 岩石力学与工程学报, 2013, 32(6): 1279-1287. doi: 10.3969/j.issn.1000-6915.2013.06.023 HU Ran, CHEN Yifeng, ZHOU Chuangbing. A relative hydraulic conductivity model for unsaturated deformable soils[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(6): 1279-1287. (in Chinese) doi: 10.3969/j.issn.1000-6915.2013.06.023
[26]	张雪东, 赵成刚, 刘艳, 等. 变形对土水特征曲线影响规律模拟研究[J]. 土木工程学报, 2011, 44(7): 119-126. ZHANG Xuedong, ZHAO Chenggang, LIU Yan, et al. Modeling study of the relationship between deformation and water retention curve[J]. China Civil Engineering Journal, 2011, 44(7): 119-126. (in Chinese)
[27]	YANG Y B, HUNG H H. A 2.5D finite/infinite element approach for modelling visco-elastic bodies subjected to moving loads[J]. International Journal for Numerical Methods in Engineering, 2001, 51(11): 1317-1336. doi: 10.1002/nme.208
[28]	王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003. WANG Xucheng. Finite Element Method[M]. Beijing: Tsinghua University Press, 2003. (in Chinese)
[29]	HALL L. Simulations and analyses of train-induced ground vibrations in finite element models[J]. Soil Dynamics and Earthquake Engineering, 2003, 23(5): 403-413. doi: 10.1016/S0267-7261(02)00209-9