基于FLAC3D平台的边坡非饱和降雨入渗分析

蒋中明 ，熊小虎 ，曾 铃

(1. 长沙理工大学 水利工程学院，长沙 410004；2. 长沙理工大学 水沙科学与水灾害防治湖南省重点实验室，长沙 410004；3. 长沙理工大学 土木与建筑学院，长沙 410004)

1 引 言

降雨入渗是诱发边坡失稳的主要因素之一[1]。许多工程实践揭示了边坡失稳多发生在雨季[2]。降雨导致边坡原有非饱和区负孔隙水压力发生强烈的变化，表层岩土体不同部位出现暂态饱和区，增大该区域边坡土体自重，导致下滑力增加，对边坡稳定性不利。长期降雨情况下，坡体内地下水位还会出现大幅度的上升[3]，增大边坡岩土体中的孔隙水压力，降低有效应力，导致该区域抗剪强度降低。因此，对边坡进行降雨入渗过程分析是十分必要的。为了分析边坡在降雨条件下的饱和非饱和渗流过程，许多学者采用自编软件或各种商业软件对边坡降雨入渗饱和非饱和渗流过程进行了研究[4－7]。近年来，基于有限差分方法的FLAC3D软件因其专业性针对性强，在岩土工程分析中得到越来越广泛的应用，学习和使用的科技人员越来越多。然而FLAC3D软件在非饱和渗流分析方面存在很大的局限性，即FLAC3D在数值计算过程中始终将负压区的饱和度强制置为1，从而使得计算过程中，非饱和区的渗透系数仍然采用饱和区的渗透系数，这与非饱和渗流理论是相悖的。为弥补这一不足，部分学者采用不同的技术对 FLAC3D的非饱和渗流计算方法进行了改进，取得了较好的效果[8]。为更加充分利用FLAC3D软件平台，并进行流-固-热等多功能计算（部分科技人员采用其他软件进行非饱和渗流分析，然后将渗流计算结果导入 FLAC3D进行力学计算），同时降低软件购置费或避免软件功能的重复购置，本文在非饱和渗流分析基本理论的认识基础上，结合对FLAC3D软件渗流分析模块的解析，采用FISH语言编写降雨入渗边界处理函数和非饱和区渗透系数实时更新函数，从而实现了降雨条件下边坡饱和非饱和渗透过程的模拟，为大型复杂三维边坡降雨入渗饱和非饱和渗流场分析及降雨条件下复杂三维边坡的稳定性分析提供一条可行的途径。

2 饱和非饱和渗流基本原理[9]

多孔连续介质饱和非饱和渗流方程如下：

式中：qi为单位流量向量；Kij为渗透系数张量；kr与饱和度s相关的相对渗透系数，饱和区，kr=1，非饱和区0

由式（1）可知，饱和渗流可以看成是饱和度为1.0时的非饱和渗流。非饱和渗流分析的关键在于非饱和区渗透系数依赖于饱和度变化而变化，同时，非饱和区的孔隙水压力可以表述为负值，也就是基质吸力。

3 FLAC3D非饱和降雨入渗分析功能的二次开发

3.1 FLAC3D软件饱和非饱和渗流分析方法[10]

FLAC3D软件在进行饱和非饱和计算时，都采用达西定律，即

式中：Kij为饱和渗透系数张量； kr(s)= s2(3 - 2s)；ρf为流体密度；xj为笛卡尔坐标分量；gj为重力加速度分量； i = 1,3, j = 1,3,l = 1,3。

在渗流计算过程中，流体压力、饱和度与流体体积改变量的关系式如下：

式中：M为比奥模量；n为孔隙率；t为时间；ζ为流体扩散作用引起的单位体积多孔介质中流体体积改变量，当流体不可压缩时，其值等于单元流量改变量。

在进行饱和渗流计算情况下，FLAC3D软件直接将节点上的饱和度置为 1.0[10]，于是由式（3）得到下一计算时刻节点上孔隙水压力的更新计算公式为

在进行非饱和渗流计算时，FLAC3D将节点上孔隙水压力等于0，然后由式（3）得到计算节点饱和度的更新值为

由此可见，FLAC3D平台非饱和状态（饱和度）变化计算的主要依据是计算微元中流体体积改变量的变化值，而这与饱和度的定义相一致。

3.2 饱和非饱和渗流分析功能的FISH函数开发

通过上述分析可知，非饱和渗流计算的关键是正确的获得饱和度与负压关系，并根据饱和度计算出非饱和区的渗透系数，然后在求解过程中适时调整单元渗透系数即可实现非饱和渗流过程的分析。

对于非饱和土，1980年，Van Genuchten提出了土体中体积含水率与负压的四参数关系方程[11]，即

式中：θ为体积含水率；θr为残余体积含水率；θ为饱和体积含水率；a、m′、n′拟合参数，对于一般土体，GEO-Slope软件采用默认值a=100，m′= 1，n′=2，对岩体，a值可取10。

根据体积含水率与饱和度关系θ=n s，由（6）式得到饱和度与负孔隙水压力的关系式为

式中：sr为残余饱和度。

由此可见，非饱和渗流分析的实质是：①计算过程中非饱和区的渗透系数小于饱和区渗透系数，饱和度越小，渗透系数越小；②负孔隙水压力（基质吸力）与饱和度存在对应函数关系，对不同岩土体，吸湿和脱湿过程的函数关系可能不同。为简化起见，也有在吸湿和脱湿过程中采用相同函数的做法。

FLAC3D通过设置流体抗拉强度来允许负孔隙水压力的产生与发展，这就为利用FLAC3D进行非饱和渗流分析成为可能。在 FLAC3D计算中，由式（5）可知，在非饱和状态下，FLAC3D软件直接将负孔隙水压力置 0，然后根据节点流体体积的改变量来计算饱和度的增量；非饱和负压的计算也根据节点代表流体体积的改变量用式（4）来完成。当流出节点的流体体积大于流入节点的流体体积时，流体体积改变量为负值（导致孔隙空间不能被流体全部充填，出现非饱和状态），从而计算出的孔隙水压力为负值。从上述过程看，FLAC3D的负压形成机制是合理的、正确的。正确进行非饱和渗流分析的关键是确保非饱和区渗透系数随饱和度的变化而变化。因此，如果能做到FLAC3D在非饱和渗流计算时段内，能自动根据上一计算增量时间步的结果来调整非饱和区的渗透系数，就能够完全实现利用FLAC3D进行非饱和渗流计算。FLAC3D的FISH语言正好提供了 FISHCALL命令来实现根据计算步的结果来自动调整各种FISH变量和FLAC3D内置变量（诸如各种力学参数和渗流计算参数等）的途径。这个过程与增量法有限元软件处理相同问题的思路完全一致。

根据上述思想，利用FISH语言在FLAC3D中实现非饱和渗流过程的方法如下：

（1）设置流体抗拉强度，允许渗流计算过程中因节点流量负流入而形成负压区。

（2）通过FISH内置变量z_pp(pz)获取单元的负孔隙水压力值（该值根据节点孔隙水压力，在FLAC3D中自动插值获得），然后用式（7）计算单元饱和度。

（3）根据式（7）计算得到单元饱和度，利用相对渗透系数计算公式计算单元相对渗透系数。

（4）对负压区单元的渗透系数进行修正：即对非饱和区单元的饱和渗透系数（计算输入值）乘以第3步计算得到的kr(s)值，然后再赋值给该单元，从而实现非饱和区渗透系数的修改。值得注意的是，通过式（7）计算得到的饱和度必须另外开辟存储单元。

经过上述4个步骤，即可在FLAC3D中实现一般意义上的非饱和渗流计算过程。值得一提的是，上述过程仅对渗流模块有效，而负压引起的土体有效应力的增加，需要在力学计算步中通过自编FISH函数另行处理。图1为FLAC3D在增量计算时步中的非饱和单元渗透系数修正计算FISH程序框图。

图1 非饱和单元渗透系数修正计算框图Fig. 1 Chart for calculation of unsaturated coefficient of permeability

3.3 降雨入渗及出渗边界功能的FISH函数开发

降雨入渗边坡条件实际上为一个动态的边界条件。地面产流前，在降雨入渗过程中，当降雨强度小于边坡表层土体最大入渗率时，入渗量取降雨强度，当降雨强度大于边坡表层土体最大入渗率时，入渗量取最大入渗率。本次研究最大入渗率按饱和渗透系数考虑。地表产流前的边坡表面施加流量边界，为第2类边界条件。地面产流后，土体表层单元处于饱和状态，此时在表层饱和单元处施加一个压力很小的压力边界或者取表面压力为 0，即第 1类边界条件。

降雨停止后，由于降雨入渗进入斜坡体内的水体在重力作用下，继续保持向边坡体内和下部渗流的趋势，在边坡斜坡表面，坡体内的水分既可向边界外溢出，也可能在重力作用下继续下渗。当斜坡上点因坡体内水分溢出而达到正孔隙水压力时，边坡表层节点的边界调整为压力为0的压力边界。降雨停止后边坡表面上的孔隙水压力在重力作用下变为负值时，渗流边界改为零流量边界，即坡体内的水分不再溢出。

FLAC3D可以直接利用的边界有 3种：①固定流量边界；②固定压力边界；③渗漏边界。为了实现边坡降雨条件的模拟，需要对边坡表面的边界状况进行实时调整，以实现降雨入渗和边坡坡面在降雨停止后的出渗模拟。FISH函数的实施过程如下：

（1）在地表面生成interface单元，利用 interface单元信息自动寻找表层实体单元(zone)，获得地表不同部位的饱和渗透系数ksat（控制最大入渗率）。

（2）通过interface单元信息，寻找地表表层每一个实体单元(zone)的外表面，并在该表面上施加流量边界，边界上的流量数值可以随计算时间实施调整，以模拟不同时期降雨强度的变化；在每一个计算时间步，对比降雨强度值与边坡岩土体的入渗率（本研究取边坡表层实体单元饱和渗透系数）之间的大小关系，取两者之间的小值作为边界上的流量输入值qs。

（3）通过interface节点信息，定位地表单元的节点位置，并在每一个计算步中，判断地表单元节点的孔隙水压力ps是否大于0。如地表节点的孔隙水压力大于 0，表明该节点达到饱和状态，地表将出现积水，因此，在后一计算时间步中，将该节点的流量边界修改为压力边界（可取压力为 0；考虑边坡表层积水驻留时间短，故可以按照一薄层水厚度换算成水压力，该值一般较小）。

（4）降雨停止后，移除地表单元表面的流量边界。对地表节点的孔隙水压力进行实时监控，当地表节点孔隙水压力大于0时，将该节点孔隙水压力强制置零。节点孔隙水压力为0的部位即为坡体内的地下水溢出点。

将上述过程的编写为 FISH函数。然后利用FISHCALL命令，即可以实现FALC3D在每一流体计算时步内对降雨边界的模拟。图2为边坡降雨入渗边界及出渗边界模拟计算的FISH程序框图。

图2 降雨入渗边界及出渗边界程序框图Fig. 2 Chart for determination of rainfall boundary

4 算例分析

4.1 计算模型

选取某实际路堑边坡进行渗流场分析。图3为网格模型，X坐标在坡脚处沿水平方向指向坡内，Z方向在坡脚处沿竖直方向向上为正。计算网格单元数量为1 591，节点数量为3 356。根据现场压水试验成果，边坡岩体饱和渗透系数为2.5×10-4cm/s，非饱和区渗透系数按照式（8）所示的Van-Genuchten模型确定，模型计算中的拟合参数取a=10，m=0.5，n= 2，残余饱和度为0.05。初始孔隙水压力场根据边坡体上稳定的钻孔水位和地表处基质吸力（取-200 kPa，地表各点相同），采用改进后饱和非饱和分析模块进行拟合生成。图3同时给出了饱和非饱和初始渗压场分布，图中虚线表示初始状态下的非饱和区负孔隙水压力等值线，孔隙水压力为0的等值线为初始地下水位线。

计算条件：连续降雨10 d，降雨强度为3×10-6m/s；降雨停止后，继续计算5 d，以查看雨后边坡孔隙水压力的变化。

由于 FLAC3D只能给出等值云图，故本文将FLAC3D计算结果导出后，用Tecplot软件进行后处理，并绘制等值线。

图3 边坡计算网格及初始孔隙水压力场分布图(单位: kPa)Fig.3 Calculation grid of slope and initial pore water pressure distribution(unit: kPa)

4.2 饱和非饱和孔隙水压力场变化过程

图4为边坡在降雨过程中和雨后5 d时间内的孔隙水压力分布变化图。连续降雨5 d后，边坡坡体、坡顶和坡脚土体已经由降雨前的负孔隙水压力状态转化为正孔隙水压力状态，且与坡体下部中的正孔隙水压力区连通；坡体以上正孔隙水压力区范围相对较小。降雨10 d后，坡体中的负孔隙水压力区继续减少，但相对于前5 d的负压区减少速度要缓慢；边坡坡体上部中的正压区范围继续扩大。降雨结束时，边坡中的孔隙水压力场分布规律与文献[12]相同。

图4(c)、4(d)、4(e)表明：降雨停止后，边坡坡顶表层以下暂态饱和区中的水分在重力作用下，继续向边坡下部渗透；而在斜坡面表面以下暂态饱和区中的水分在重力作用下一方面继续向边坡体更低位置渗透，另一方面暂态饱和区中的水分通过坡面出渗，向坡体外部排出。降雨停止1 d，斜坡表层和坡顶以下1~2 m范围内重新出现负压区；而雨后3~5 d后，负压区大幅度增加，边坡原来连通的暂态饱和区被明显分割成两部分。

图4(b)、4(c)清晰地表明：在降雨过程中及降雨停止后，边坡体中出现了较大范围的暂态饱和区，饱和区内出现正孔隙水压力，即暂态水压力。暂态水压力分布规律：坡面孔隙水压力很小（近似为0）；坡面附近数值较小；暂态饱和区下缘与非饱和区接触面上的孔隙水压力为 0；最大值出现在暂态饱和区中部偏上部位。由于边坡坡脚处受坡脚平台降雨入渗补给和斜坡体下岩土体中孔隙水下渗补给的双重影响，因此，地下水位上升幅度较大。坡顶部位铅直下方的土体由于仅仅受坡顶面上的降雨入渗补给，降雨过程中始终存在非饱和区，使得暂态饱和区很难延伸到潜水面的位置，潜水面的升高也将会很小（通过非饱和区补给到初始潜水面上引起的潜水面升高），因此，该部位饱和区的渗流场变化不大。

图5为不同时期边坡的饱和度。饱和度小于1.0的区域表示边坡中出现的非饱和区，其对应的孔隙水压力小于0。饱和度为1.0的位置与孔隙水压力为0的等值线相对应。由此可见，边坡中的地下水位线的动态变化过程是十分复杂的。

图4 孔隙水压力分布变化图(单位: kPa)Fig.4 Variation of pore water pressure contours(unit: kPa)

图5 边坡饱和度等值图Fig.5 Variation of saturation degree contours

4.3 地表测点孔隙水压力变化过程

图6为边坡测点孔隙水压力在降雨过程中及降雨停止后的变化过程图。降雨过程中，边坡表面坡脚（测点1）、坡腰（测点2）及坡顶（测点3）部位测点的孔隙水压力在降雨入渗作用下逐渐升高，且坡顶部位测点孔隙水压力的上升速度明显小于坡脚和坡腰，其原因是坡脚和坡腰测点的孔隙水压力除了受到该点降雨入渗作用的影响，同时，还受到通过斜坡坡体内部扩散来的水分的影响，因此，其孔隙水压力上升相对较快。但地表附近因饱和出现正孔隙水压力后，如果边坡表面积水能及时排出，那么地表孔隙水压力将维持在大气压力值，即零孔隙水压力值。降雨停止后，边坡顶部测点3的孔隙水压力在降雨停止2 h后开始由正压转变为负压，测点2在降雨停止10 h后，开始出现负压。测点出现负压后，孔隙水压力逐渐降低，大约2 d后，孔隙水压力降低速率越来越缓慢，并趋向与平稳。地表测点1压力在降雨停止5 d内，因该点一致处于饱和状态，所以其压力一直维持为0。

图6 边坡测点孔隙水压力时间过程线Fig.6 Pore pressure of monitor points vs. time

4.4 对比分析

FLAC3D默认的饱和非饱和渗流分析方法中，负压区中单元渗透系数是固定不变的[10]，这种处理方式对降雨非饱和渗流的影响是很大的。图7为利用FLAC3D程序默认处理方式计算得到的边坡在降雨结束时和雨后5 d情况下的孔隙水压力分布等值图。对比同一时刻条件下的孔隙水压力等值图（图7(a)和图4(b)、图7(b)和图4(e)）可知，在其他条件相同的情况下，边坡中暂态饱和区及负压区的孔隙水压力分布存在较大的差异。特别是降雨停止后，在不对非饱和区渗透系数进行折减的情况下，降雨入渗进入边坡的雨水快速下渗到了坡体内部，因此，边坡中的地下水位快速大幅度上升，且地下水位线以上暂态饱和区完全消失。由此可见，FLAC3D默认的渗流分析模式不能正确反映边坡降雨入渗非饱和渗流场的分布与渗流场的发展与变化。

图7 负压区渗透系数固定时孔隙水压力场(单位: kPa)Fig.7 Pore pressure contours without modification of the coefficient of permeability in negative pressure zone(unit: kPa)

5 结 论

（1）利用 FLAC3D计算得到的负压来计算相对渗透系数，并在FLAC3D求解过程中的下一时间步增量中对负压区的渗透系数进行修改，从而完善FLAC3D的非饱和渗流分析功能。

（2）通过在当前计算时间步中引入降雨入渗边界，可实现了复杂三维地形条件下FLAC3D降雨入渗边界的自动模拟。

（3）利用编写的非饱和降雨入渗分析修正模块，分析了算例边坡在降雨条件下的非饱和渗流过程，并与FLAC3D内置渗流分析模块成果进行对比。计算成果表明，在对 FLAC3D软件现有功能修改的基础上，完全可以正确实现边坡降雨非饱和渗流过程的数值模拟过程，为利用FLAC3D软件实现复杂三维边坡降雨入渗分析提供了基础。

[1]NG C W, SHI Q. A numerical investigation of the stability of unsaturated soil slopes subjected to transient seepage[J].Computers and Geotechnics, 1998, 22(1): 1－28.

[2]孙广忠. 中国典型滑坡[M]. 北京: 科学出版社, 1998.

[3]张卓, 练继建, 杨晓慧. 不同降雨强度下岩体边坡的渗流场分析[J]. 水力发电学报, 2006, 25(5): 27－30.ZHANG Zhuo, LIAN Ji-jian, YANG Xiao-hui. The analysis of seepage field in rock slope under different rainfall degrees[J]. Journal of Hydroelectric Engineering,2006, 25(5): 27－30.

[4]戚国庆, 黄润秋, 速宝玉, 等. 岩质边坡降雨入渗过程的数值模拟[J]. 岩石力学与工程学报, 2003, 22(4): 625－629.QI Guo-qing, HUANG Run-qiu, SU Bao-yu, et al.Numerical simulation on rainfall infiltration on rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(4): 625－629.

[5]周桂云, 李同春. 饱和-非饱和非稳定渗流作用下岩质边坡稳定性分析[J]. 水电能源科学, 2006, 24(5): 79－82.ZHOU Gui-yun, LI Tong-chun. The rock slope stability analysis under saturated-unsaturated seepage[J]. Water Resources and Power, 2006, 24(5): 79－82.

[6]石露, 李小春, 白冰, 等. 降雨条件下露天平台边坡的稳定性研究[J]. 岩土力学, 2012, 33(5): 1519－1526.SHI Lu, LI Xiao-chun, BAI Bing, et al. Research on open-pit platform slope stability under rainfall condition[J]. Rock and Soil Mechanics, 2012, 33(5):1519－1526.

[7]付宏渊, 曾玲, 王桂尧, 等. 降雨入渗条件下软岩边坡稳定性分析[J]. 岩土力学, 2012, 33(8): 2359－2365.FU Hong-yuan, ZENG Ling, WANG Gui-yao, et al.Stability analysis of soft rock slope under rainfall infiltration[J]. Rock and Soil Mechanics, 2012, 33(8):2359－2365.

[8]李毅, 伍嘉, 李坤. 基于 FLAC3D的饱和-非饱和渗流分析[J]. 岩土力学, 2012, 33(2): 617－622.LI Yi, WU Jia, LI Kun. Saturated-unsaturated seepage analysis based on FLAC3D[J]. Rock and Soil Mechanics,2012, 33(2): 617－622.

[9]吴梦喜. 饱和-非饱和土中渗流Richards方程有限元算法[J]. 水利学报, 2009, 40(10): 1274－1279.WU Meng-xi. Finite element algorithm for Richards equation for saturated-unsaturated seepage flow[J].Journal of Hydraulic Engineering, 2009, 40(10): 1274－1279.

[10]Itasca Consulting Group, Inc. Online Manual of FLAC3DFast Lagrangian Analysis of Continua in 3 Dimensions[M].Version 3.0. Minneapolis: Itasca Inc., 2003.

[11]VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(5): 892－898.

[12]胡冉, 陈益峰, 周创兵. 降雨入渗过程中土质边坡的固-液-气三相耦合分析[J]. 中国科学: 技术科学, 2011,41(11): 1469－1482.HU Ran, CHEN Yi-feng, ZHOU Chuang-bing.Solid-liquid-gas three-phase coupling analysis of soil slope during rain fall infiltration[J]. Scientia Sinica Technologica, 2011, 41(11): 1469－1482.