铁路路基降雨入渗过程中水-气两相渗流特性研究

2022-12-13贾羽王晅丁瑜张家生陈晓斌

铁道科学与工程学报 2022年11期

贾羽，王晅, ，丁瑜，张家生, ，陈晓斌,

(1.中南大学土木工程学院，湖南长沙 410075；2.中南大学高速铁路建造技术国家工程研究中心，湖南长沙 410075)

铁路路基作为列车和轨道结构的支承体，直接暴露在大气环境下，易受到降雨入渗的长期作用[1−2]。通常情况下，路基填料处于非饱和状态，降雨入渗到路基中，导致路基孔隙内的水和气体的流动，因此，降雨入渗的实质是路基内部的水−气两相渗流过程。以往的岩土工程非饱和渗流计算通常认为孔隙气压恒等于大气压力，忽略了孔隙气压力变化对水渗流特性的影响。KIM等[3]研究了降雨过程中土体的水−力耦合行为对非饱和土边坡稳定性的影响，计算过程中，将孔隙气压力假定为大气压力。YANG等[4]在Richards方程的基础上，建立了降雨入渗过程中边坡的水−力耦合分析模型。HAMDHAN等[5]分析了土体水力参数和降雨强度对边坡稳定性的影响，计算过程中假定孔隙气压力等于大气压力。大量实验证实，非饱和土体中孔隙气压力的变化对孔隙水渗流特性具有显著的阻滞作用[6−7]。因此，为了更好地模拟非饱和土体中水的渗流规律，国内外学者开始利用水-气两相渗流理论进行研究。LIU等[8]基于水−气两相流理论，建立了边坡降雨入渗过程中非饱和渗流模型。CHO[9]采用水气两相流分析方法，研究了强降雨引起的空气与水流相互作用对非饱和土边坡力学稳定性的影响，结果表明：由于雨水置换了非饱和带内的空气，孔隙气压力逐渐增加；孔隙气压力的增加对水的流动产生了显著的延迟效应。孙冬梅等[10−11]基于多孔介质中水−气两相渗流理论，建立了土体水−气两相渗流的数学模型，利用该模型对非饱和土质边坡的降雨入渗过程进行了模拟。上述研究多是针对降雨入渗对边坡稳定性的影响，边坡多是某种特定类型的土体，而路基属于层状结构物，每层填料的类型以及物理力学性质都有很大的差异。因此，有必要针对降雨入渗对路基各层结构水−气两相渗流特性进行研究。本文基于两相渗流理论，建立水−气两相渗流数学模型；利用该模型研究铁路路基在特大暴雨条件下水−气两相渗流规律，探讨降雨入渗对路基各层结构孔隙水压力以及孔隙气压力的影响。

1 路基降雨入渗过程中水-气两相渗流数学模型的建立

1.1 基本假设

1) 路基填料具有多孔介质性质，是由固相(s)、液相(水，w)以及气相(空气，a)组成的三相体，且孔隙间具有连通性；

2) 孔隙内水的流动以及气体的流动是由其各自的压力梯度引起的，并且均遵从Darcy定律[12]；

3) 不考虑水和气体之间的物质交换。

1.2 水−气两相渗流的控制方程

考虑水−气两相渗流时，多孔介质内水相的质量守恒方程以及气相的质量守恒方程[12]可以表示为：

式中：ϕ为多孔介质的孔隙率；Sw和Sa分别为水相和气相的饱和度；ρw和ρa分别为水相和气相的密度，kg/m3；qw和qa分别为水相和气相的源(汇)项，kg/m3；uw和ua分别为水相和气相的渗流速度矢量，m/s，并且分别满足Darcy定律[13]：

式中：k为多孔介质的绝对渗透率，m2；krw和kra分别为水相和气相的相对渗透率，0≤krw(kra)≤1；μw和μa分别为水相和气相的动力黏滞系数，(N∙s)/m2；pw和pa分别为水相和气相的孔隙压力，Pa；g为重力加速度矢量，N/kg。

考虑水相以及气相的压缩性，其密度可以分别表示为：

式中：Cw和Ca分别为水相和气相的体积压缩系数，Pa−1。

将式(3)和式(5)代入式(1)，式(4)和式(6)代入式(2)分别得到水-气两相渗流的控制方程：

1.3 辅助方程

上一节得到的水-气两相渗流的基本控制方程中有4个基本未知量(Sw，pw，Sa，pa)，求解结果还需要一些辅助方程。

1) 水−气两相流饱和度方程

2) 基质吸力−饱和度关系曲线

式中：pc为基质吸力，Pa，是饱和度的函数。

本文采用Van Genuchten模型[14]表征土体基质吸力−饱和度关系曲线：

式中：m和n为拟合参数，并且满足m=1-1/n；α为拟合参数；Se为有效水饱和度，定义为：

式中：Swr为残余水饱和度；Sws为饱和水饱和度。

3) 相对渗透率−饱和度关系曲线

表征水相和气相的相对渗透率−饱和度关系曲线分别采用Van Genuchten-Mualem模型[15]和Brook-Corey模型[16]：

由公式(7)，(11)，(12)和(13)可得水相渗流的控制方程，由公式(8)，(9)，(11)，(12)和(14)可得气相的渗流控制方程，分别为：

式中：

1.4 边界条件和初始条件

上述数学模型是高度非线性的，在求解的过程中需要给定边界条件和初始条件。求解时，将孔隙水压力pw和孔隙气压力pa作为基本变量。

边界条件：边界条件通常有第1类边界条件，也称为Dirichlet boundary；以及第2类边界条件，也称为Neumann boundary。

第1类边界条件可以表示为：

第2类边界条件可以表示为：

qw=-vws∙nonΓqw；qa=-vas∙nonΓqa。

初始条件：t=0时，应满足以下条件：pw=p0w，pa=p0a。

1.5 数学模型的求解

利用多物理场仿真软件COMSOL Multiphys‐ics®对上述模型进行求解，采用Galerkin法对偏微分方程组进行离散化，使用内置非线性微分代数方程求解器作为时间积分器，控制方程组中的时间项离散采用欧拉向后差分法。将上述数学模型写入PDE模块中，系数分别为：

2 模型验证

2.1 Liakopoulos排水实验

Liakopoulos[17]排水实验曾得到众多学者的广泛验证和模拟[18−21]，本节采用上一节建立的水-气两相渗流数学模型对Liakopoulos排水实验进行模拟，以验证数学模型的正确性。

Liakopoulos排水实验采用的是Del Monte 砂，试样尺寸为：高度1 m，直径0.1 m，如图1所示。试样底部安装高渗透性透水石，侧面为不透水圆筒。实验开始前，试样上部有稳定水流供应，整个砂柱处于完全饱和状态。实验开始时，切断试样上部的供水。砂柱内的水在重力作用下逐渐从下部排出。在实验过程中，记录了砂柱内部的孔隙水压力以及砂柱底部排水速率等数据。

图1 Liakopoulos排水实验示意图Fig.1 Schematic of the Liakopoulos drainage test

2.2 Liakopoulos 排水实验的数值计算模型

数值计算模型如图1所示，模型计算采用物理参数如表1 所示[22−23]。

表1 数值计算物理参数Table 1 Physical parameters for numerical calculation

本文在计算中以标准大气压力值作为计算参考值，孔隙水压力和孔隙气压力均采用相对压力值。根据Liakopoulos排水实验，在数值模型计算中采用的初始条件为：pw=0，pa=0。边界条件定义为：1) 试样上部为不透水、透气边界，即qw=0,pa=0；2) 试样底部为排水、透气边界，即pw=0，pa=0；3) 试样侧面为不透水、不透气边界，即qw=0，qa=0。

2.3 计算结果的验证分析

单相渗流以及水-气两相渗流模型计算得出的砂柱不同高度处孔隙水压力演化曲线如图2所示。水−气两相渗流计算模型得出的砂柱底部孔隙水渗流速度演化曲线如图3所示。

图2 2个不同高度处孔隙水压力演化曲线Fig.2 Evolutions of pore water pressure at two different heights

图3 砂柱底部孔隙水渗流速度演化曲线Fig.3 Evolutions of water outflow rate through the bottom surface of the sand column

从图2中看出：1) 随着时间的增长，孔隙水压力非线性减小：前期减小的速率较大，中期减小的速率逐渐变缓，后期孔隙水压力逐渐趋于定值；并且高度越高，孔隙水压力的值越小。2) 单相渗流计算模型中，假设孔隙器压力恒定为0，为考虑其变化对孔隙水压力的影响，因此，同一时间，单相渗流计算得到的孔隙水压力变化更迅速；与单相渗流计算结果相比，水−气两相渗流计算结果与实测数据吻合更好，尤其是在t<40 min内。说明在进行非饱和渗流计算时，应该考虑孔隙气压力对孔隙水压力变化的影响。

从图3中看出：1) 随着时间的增加，砂柱底部孔隙水渗流速度先快速减小，后缓慢减小。2) 与实验数据相比，模型计算值在初期(即t<40 min)与试验数据吻合较好；随着时间的增大，计算值与实验数据存在一定的误差。

综合分析图2和图3可知：本文的水-气两相渗流数学模型可以很好地预测非饱和渗流过程。

3 路基降雨入渗过程中水-气两相渗流特性分析

3.1 有限元模型

根据我国《高速铁路设计规范》(TB10621—2014)[24]，参照京沪高速铁路路基横断面结构，建立如图4所示的二维路基降雨入渗计算模型。路基结构从上至下依次为：基床表层(厚度0.4 m)，基床底层(厚度2.3 m)，路基本体(厚度2.0 m)以及地基(厚度5.0 m)。基床表层和基床底层的边坡坡度为1:1.5。为分析路基不同位置处的渗流特性，设置了4个观测点，如图4所示。模型采用自由三角形网格，共划分了6 251个网格，平均单元质量0.92。

3.2 路基各层填料力学参数

基床表层填料为级配碎石，基床底层填料和路基本体填料为改良粗颗粒土，地基土为低液限粉土。模型计算所采用的物理参数如表2所示[2,22−23,25]。

表2 数值计算物理参数Table 2 Physical parameters for numerical calculation

3.3 初始条件和边界条件

假定在计算过程中不考虑地下水的影响，路基模型计算的初始条件为：pw根据每种填料的初始水饱和度计算，pa=0。

路基模型计算的边界条件如图4所示。由于基床表层上混凝土支撑层的覆盖作用，因此DE设置为不透水、不透气边界，地基下表面IJ也假定为不透水、不透气边界；上表面ABCD以及EFGH为降雨入渗边界以及透气边界；其余为透水、透气边界。

图4 路基降雨入渗计算模型Fig.4 Illustration of subgrade model under rainfall infiltration

路基降雨入渗的过程与降雨强度和土体的渗透性能相关[26]。假设降雨强度为ε(t)，边坡坡面外法线方向为n=(nx,ny,nz)，当降雨强度小于土体渗透系数时，认为雨水全部入渗，此时的降雨入渗条件简化为流量边界条件，那么，边坡坡面上雨水实际入渗量qs(t)与降雨强度的关系满足[27−28]：

当降雨强度大于土体渗透系数时，雨水在土体表面形成积水(积水水头不高)或形成表面径流，此时可将降雨入渗条件简化为压力条件。

本文计算过程中模拟特大暴雨情况，假设降雨强度为 3.47×10−6m/s(300 mm/24 h)，降雨持续时间为9 h[29]。

3.4 计算结果与分析

3.4.1 路基降雨入渗过程中孔隙水压力的变化

图5显示了水−气两相渗流模型计算得到的不同降雨时间路基内孔隙水压力分布云图。从图中可见，在降雨初期(如图5(a))，基床表层两侧上表面以及基床表层边坡处因为降雨的入渗，孔隙水压力最先增大为正值，说明此处达到饱和状态。随着降雨的进行(如图5(b))，基床表层两侧入渗的水分逐渐向基床表层下部以及中间部位渗流，使基床表层下部以及中间部位的孔隙水压力逐渐增大达到饱和状态；当降雨进行至9 h左右时(如图5(c))，基床表层两侧下部的饱和区域基本贯通，基床表层已基本处于饱和状态，此时基床表层和基床底层界面处滞留水分。若此时遇到在列车荷载的反复作用，则界面处填料中细颗粒会被软化进而形成泥浆，易产生翻浆冒泥等路基病害。随着降雨的继续进行，路基内部的饱和区域逐渐向下延展，由于基床底层、路基本体以及地基等处填料的渗透率只有基床表层的1/100甚至更小，因此，雨水入渗速率变缓慢，直至计算结束(如图5(d))，基床底层、路基本体以及地基两侧降雨入渗边界只有很小厚度范围内处于饱和状态。这也说明，降雨过程中水分很难直接迅速地渗入基床底层、路基本体以及地基内部。

图5 不同降雨时间路基内孔隙水压力分布云图Fig.5 Cloud maps of pore water pressure at different times

图6显示了路基内不同观测点处(见图4)孔隙水压力演化图。由图6可知：观测点1处随着降雨时间的增加由负孔隙水压力逐渐增加至正孔隙水压力，t=2 h左右达到饱和状态，随着时间的继续增长，孔隙水压力先逐渐增大后缓慢减小；观测点2处的孔隙水压力变化趋势与1处相似，t=20 h左右达到饱和状态。在雨水还没有入渗到观测点3和观测点4时，此处的孔隙水压力由于孔隙气压力的影响(如图8所示)呈先增大后减小的变化趋势[9]。

图6 路基内不同观测点处孔隙水压力演化Fig.6 Evolutions of pore water pressure at different observation points

图8 路基内不同观测点处孔隙气压力演化Fig.8 Evolutions of pore air pressure at different observation points

3.4.2 路基降雨入渗过程中孔隙气压力的变化

图7显示了不同降雨时间路基内孔隙气压力分布云图。在降雨初期，随着基床表层两侧以及基床表层边坡降雨的入渗，挤压土体孔隙内的气体，基床表层的孔隙气压力增大，随着时间的增长，受挤压的孔隙气体逐渐向下消散，基床表层的孔隙气压力逐渐减小。由于地基处填料的渗透率很小，雨水入渗速率缓慢，其孔隙气压力变化缓慢，直至计算结束，地基两侧边界处的孔隙气压力尚未消散完全。图8显示了路基内不同观测点处(见图4)孔隙气压力演化。由图8可知：4个观测点处的孔隙气压力的变化规律基本一致：随着降雨时间的增加，孔隙气压力由于入渗雨水的挤压作用首先呈现增长的趋势；随着受挤压的孔隙气体逐渐向下扩散，孔隙气压力逐渐减小直至趋近于大气压力。结合图6中观测点2处3～15 h内孔隙水压力产生下降的趋势，说明孔隙气压力的增大对水流的运动存在一定阻滞作用。