基于并行FDM模型和流网计算的土体渗流

摘要

常规二维稳态流网计算得到的渗流安全系数接近于临界值，将导致设计过程中采取超量设计的方法对结构进行冗余防护。因此，提出了基于有限差分法（FDM）模型对流网进行二次修正的方法，将计算截面离散化，利用FDM步长单元小、迭代计算次数多的优势获得新的等势线和流线，并用修正后的水头梯度和渗流路径重新计算渗流量和安全系数。结果表明：基于FDM模型修正的每日渗流量为2.77 m³/m，与二维稳态流网的计算值3.29 m³/m基本一致，但经FDM模型修正后的渗流安全系数为2.7，远大于二维稳态流网的计算值1.6。由此说明，经FDM模型修正后的流网计算结果精度较高，结构安全性更可靠。

关键词

流网; 有限差分法; FDM模型; 渗流安全系数

0 引言

土壤是岩石经过风化、破碎和剥蚀后自然堆积成层的产物。由于土层内部的土粒是各项异性的非均匀体，因此，土层内部水平方向的渗透系数与竖直方向的渗透系数在数值上不相等，造成多孔介质内流体的渗流路径具有非唯一性，如单相渗流、多相渗流和双重介质渗流等。长期的渗流作用将改变土体内部孔隙水压力和浮力，从而使路堤边坡发生滑动失稳、管涌等工程灾害。早在半个世纪前，研究人员便对土的渗流安全系数和影响因素进行了一系列的研究。Barnes^［1］在1986年提出地下水流动的方式会改变流体内部可溶性硫酸盐浓度，从而腐蚀坝体和路堤边坡。Knappett等^［2］在1992年利用流网法中的等势线和流线对某水库土坝的二维渗流问题进行分析，发现通过流网法计算得到的安全系数能够满足规范要求。Knappett等^［3］在2019年描述了如何通过数学模型从拉普拉斯方程的有限差分形式中获得解，并根据有限差分法的解进一步计算了二维稳态渗流安全系数，发现传统的流网法计算值偏于保守，有限差分法计算得到的安全系数在满足规范要求的前提下，小于流网法计算值，使工程的经济性更加合理。

张超^［4］研究得出浑水会改变土的渗透特性，水头、浓度和不均匀系数越大，土的渗透系数越小。刘姚军^［5］提出土-岩界面的接触特性对层面渗流存在一定的影响，当土-岩界面剪胀效应越明显时，接触界面渗流通道相对越大，土-岩界面过流能力越强，非线性渗流特征也更加明显，导致的不稳定因素也愈多。李海鹏等^［6⁃7］提出上游土体的分期加载和流固耦合会增加坝体的渗流作用，同时随着上游土层顶角角度的增加，渗流对路堤或者挡墙墙体的作用越明显。徐令宇等^［8］提出一种数学方程和数值模型相结合的方法，对渗透固结理论的偏微分计算方程进行改进，并通过有限元软件验证了计算结果的精度，结果表明改进后的方程计算结果优于传统的太沙基固结微分方程结果。

本文通过试错法绘制流网，求得坝体边坡初步的渗流安全系数后，对边界条件进行简化，然后利用有限差分法将土体单元离散并建立数值模型；再通过拉普拉斯方程进一步求解数值模型的各边界值，推导出修正后新的流网形式，计算出修正后的渗流安全系数，从而判断结构的安全性和稳定性。

1 流网

流网是由互相正交的流线和等水头线（等势线）所组成的网，用于反映地下水在渗流场中的运动方向和渗流速度。流网法是讨论二维稳态平面渗流问题最直接、最全面的方法，从流网图中可获得用于评估渗流场特性和工程渗流控制设计所需的渗流速度、渗透压力、水头、水力坡降、过流断面上的渗流量等参数，也可根据流网的变化特征、流线以及等水头线的变化形态，来了解和判断渗透水的渗径、历程、场中各子域间的水量互补关系、子域的相对透水性和相对潜在的可能渗透变形区等。

根据达西定律^［9］，可分解和求出流网中势函数和流函数的基本方程。

设势函数为$\varphi =kh$，则有：

$$V_x=ki=\frac{\partial \varphi }{\partial x}$$

（1）

$$V_z=ki=\frac{\partial \varphi }{\partial z}$$

（2）

同时必存在一个流函数：

$$\frac{\partial \psi }{\partial z}=-\frac{\partial \varphi }{\partial x}$$

（3）

$$\frac{\partial \psi }{\partial x}=-\frac{\partial \varphi }{\partial z}$$

（4）

且势函数和流函数为共轭调和函数：

$$\frac{{\partial }^{\mathrm{2}}\psi }{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}\varphi }{\partial z^{\mathrm{2}}}=\mathrm{0}$$

（5）

流网的绘制步骤^［9］大致如下：（1）首先根据边界条件确定边界上的流线和等势线；（2）根据流网的正交性反复试画、修正；（3）当存在渗流的自由水面时，需反复试画；（4）当无法分割为整数时，流槽可为非整数。虽然流网对于估算简单的由渗流引起的孔隙压力和体积流速是有用的，但对于复杂、大型工程参数计算的效率仍有待提高。

2 FDM模型（Finite Difference Method）

有限差分法是通过近似数值解法，将微分问题直接转变为代数问题。它将求解域划分为差分网格模型，使既有连续的求解域被有限个网格节点代替，借助Taylor级数等展开方式，用网格节点函数值的差商去离散控制方程中的导数，从而建立以网格模型节点上的值为未知数的线性方程组。

2.1 模型建立

模型的设置通过对时间和空间不同差分方式的组合，从网格划分的精度、差分的空间形式和时间因子等因素出发，建立1阶、2阶和高阶，或中心格式和逆风格式等，并结合实际土层、地形因素和科朗稳定条件来确定离散元的步长。Taylor级数的展开仅需着眼于节点的微分平衡，不需要求出其分布函数，用节点上待求变量的数值作为控制容积上待求变量的代表，从而得到离散化的方程。

2.2 节点分配计算

FDM模型中的离散节点包括：上-左节点、上边界、上-右节点、右边界、下-右节点、下边界、下-左节点、左边界和连接节点等。主要离散节点的时空位置和计算控制方程如表1所示。

表1  主要节点分配和计算公式表

节点类型	节点图示	控制方程
上-左		h =$\frac{h_1+h_4}{2}$
上边界		h =$\frac{h_1+h_3+2h_4}{4}$
上-右		h =$\frac{h_3+h_4}{2}$
右边界		h =$\frac{h_2+2h_3+h_4}{4}$
下-右		h =$\frac{h_2+h_3}{4}$
下边界		h =$\frac{h_1+2h_2+h_3}{4}$
下-左		h =$\frac{h_1+h_2}{2}$
左边界		h =$\frac{2h_1+h_2+h_4}{4}$
连接节点		h =$\frac{h_1+h_2+h_3+h_4}{4}$

3 模型实例计算

以深圳市龙岗区龙口水库现状使用的堤坝为例建立FDM模型，得到土层孔隙水压力的分布，并用2种方法计算下游的渗流量。根据地质勘察报告，堤坝两侧上、下游均为均质的土层，上层粉砂土厚度D₁=4 m，渗透系数k₁=1.7×10^-5 m/s；下层粉土厚度D₂=5 m，渗透系数k₂=1.4×10^-5 m/s；坝体宽度L=8 m，上游水位H₁=5 m，下游水位H₂=0 m，坝底距上游河床高度h=1.5 m。水库坝体断面图见图1；土层地质剖面图见图2。

图2  土层地质剖面图（单位：m）

3.1 流网计算

已知渗流总是由高水头流向低水头，因此根据地质勘察报告提供的参数，流线应从左往右。根据李宏伟等^［10］的研究发现，不同的流道结构对流网的结果有着不同的影响，等水头线的间隔数量也会影响渗流的最终计算。因此，结合坝体的宽度和地址条件，初步绘制和分析时，采用10间隔等势线的流网，如图3所示。

图3  试错法初步绘制的流网（单位：m）

由于2层渗透性不同，且渗流在土体中存在水平方向（k_h）和垂直方向（k_v），因此在计算中应采用等效渗透系数k'。k_h、k_v、k'的计算式为：

$$k_{\mathrm{h}}=\frac{t_{\mathrm{1}}{k}_{\mathrm{1}}+t_{\mathrm{2}}{k}_{\mathrm{2}}}{t_{\mathrm{1}}+t_{\mathrm{2}}}=\frac{\mathrm{4}\times \mathrm{17}\times {\mathrm{10}}^{-\mathrm{4}}+\mathrm{5}\times \mathrm{14}\times {\mathrm{10}}^{-\mathrm{4}}}{\mathrm{4}+\mathrm{5}}=\mathrm{1.53}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}$$

$$k_{\mathrm{v}}=\frac{t_{\mathrm{1}}+t_{\mathrm{2}}}{\frac{t_{\mathrm{1}}}{k_{\mathrm{1}}}+\frac{t_{\mathrm{2}}}{k_{\mathrm{2}}}}=\frac{\mathrm{4}+\mathrm{5}}{\frac{\mathrm{4}}{\mathrm{17}\times {\mathrm{10}}^{-\mathrm{4}}}+\frac{\mathrm{5}}{\mathrm{14}\times {\mathrm{10}}^{-\mathrm{4}}}}=\mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}$$

$$k^{\text{'}}=\sqrt[]{k_{\mathrm{h}}\times k_{\mathrm{v}}}=\sqrt[]{\mathrm{1.53}\times {\mathrm{10}}^{-\mathrm{3}}\times \mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{3}}}=\mathrm{1.525}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}$$

根据流网计算公式，每日渗流量$q=k^{\text{'}}h\frac{N_{\mathrm{f}}}{N_{\mathrm{d}}}$。其中：流线间隔数N_f=5，等势线数N_d=10。因此，试错法下流网计算得到的每日渗流量为：

q = k'h$\frac{N_{\mathrm{f}}}{N_{\mathrm{d}}}$=1.525×10^-3×（5-0）×$\frac{\mathrm{5}}{\mathrm{10}}$=3.29 m³/m

3.2 FDM模型

根据坝体断面图，以上游高水位处为坐标原点，对结构建立FDM模型。沿深度方向和水平方向步长均为0.5，得到上游河床处水头为5 m；将模型离散成860个节点进行差分法计算，离散后的模型如图4所示。

图4  离散化的FDM模型

水头的值可通过FDM数值表的迭代计算得到，即∑∆h =1.84 m，模型计算结果如图7所示。其中，等效渗透系数与之前相同。利用达西定律，可以计算出每日渗流量：

q = ki = 1.525×10^-3×1.84 = 2.48 m³/m

图5  FDM模型计算结果（单位：m）

3.3 基于FDM模型修正的流网

根据FDM表中孔隙水压力的分布得到10条等势线，分别为5.00、4.50、4.00、3.50、3.00、2.65、2.25、1.85、1.50、1.00、0.50 m。根据FDM模型计算得到的等势线对上一步绘制的流网进行修正，并绘制新的流网图，如图6所示。

图6  基于FDM模型结果绘制的流网（单位：m）

根据修正后的流网，将流线间隔数调整为N_f=4.2；等势线数保持不变，N_d=10。因此修正后的每日渗流量为：

q = k'h$\frac{N_{\mathrm{f}}}{N_{\mathrm{d}}}$=1.525×10^-3×（5-0）×$\frac{\mathrm{4.2}}{\mathrm{10}}$=2.77 m³/m

3.4 水利梯度和渗流安全系数计算

3.4.1 土层初始状态下的水利梯度

水利梯度计算式为：

$$i_{\mathrm{e}}=\frac{∆h}{∆L}=\frac{\mathrm{5}}{\mathrm{8}}=\mathrm{0.625}$$

式中：△h为坝体水头差；△L为坝体界面的宽度，即渗流路径。假设临界水利梯度i_c为1，则在临界状态下抵抗渗流的安全系数F${}_{\mathrm{b}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}}$为：

$$F_{\mathrm{b}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}}=\frac{i_{\mathrm{c}}}{i_{\mathrm{e}}}=\frac{\mathrm{1}}{\mathrm{0.625}}=\mathrm{1.6}$$

F值接近于1，这是一个逼近于临界状态的破坏值，结构可能会不足以抵抗垂直向上的压力，大坝可能会不稳定。应采取稳定处理措施，如在坝体上游加设桩，以延长渗流通道。此外，还可在下游水位处采用井点降水法，通过抽水来降低地下水位，以保持坝体的稳定。

3.4.2 基于FDM模型修正后的水利梯度

经FDM模型修正后的水利梯度计算式为：

$$i_{\mathrm{e}}=\frac{∆h}{∆S}=\frac{\mathrm{1.84}/\mathrm{10}}{\mathrm{0.5}}=\mathrm{0.37}$$

式中：△S为步长，即离散单元体内的渗流路径；△h为离散单元体内的水头差。假设临界水力梯度为1，则在临界状态下抵抗渗流的安全系数F${}_{\mathrm{b}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}}$为：

$$F_{\mathrm{b}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}}=\frac{i_c}{i_{\mathrm{e}}}=\frac{\mathrm{1}}{\mathrm{0.37}}=\mathrm{2.7}$$

基于FDM模型修正后的流网计算得到的F值远大于1，可判定大坝是稳定的，无需采取特殊措施来解决流量问题。

很明显2种方法的计算结果是不同。采用FDM计算的安全系数大于采用初步流网法计算的安全系数，即大坝可能是稳定的，而流网法计算的结果可能是不安全的。Tanaka等^［11］采用不同的渗流计算方法进行了大量试验，其计算结果对流量的误差在13%以下，对压力的误差在3%以下；其平均值在流量和压力上的误差分别为6%和1%。

2种方法的主要偏差在于水头计算方法和渗流路径计算方法的差异。有限差分法节点多，网格间距小，对模型进行迭代计算的次数多，有利于获得可靠、准确的数据^［12］。然而，流网法因其具有较高的便捷性，在实际工程中得到了广泛的应用。在关键参数缺乏的通常情况下，采用流网法可以较好地获得工程的渗流情况。

4 结语

（1）有限差分法（FDM）采用成熟的理论，结果直观、精度可靠，易于编程和并行计算，但处理不规则区域比较复杂。尽管网格生成可以使FDM适用于不规则区域，但它对离散单元体的连续性有着严格的要求。

（2）FDM模型中离散单元的渗流水头和渗流路径长度的取值在计算中占据决定性因素，单元离散精度高，取值合理、精确，可以大大降低计算误差。

（3）由达西定律中水力梯度的计算值与临界水力梯度的对比结果可知，结构可能是稳定的，但由于安全系数较低，会造成冗余设计，使得项目的建安费超出既有概算。

（4）2种方法计算的每日渗流量不同。使用FDM计算的渗流安全系数为2.7，而使用流网法计算的渗流安全系数为1.6。使用FDM计算的每日渗流量为2.77 m³/m，而使用试错法绘制的初步流网计算的每日渗流量为3.29 m³/m。经FDM模型修正后计算的每日渗流量与初布流网计算的每日渗流量高度一致，但FDM计算得到的渗流安全系数远大于流网法的计算值，更接近真实情况，更可靠。