4.4.3 饱和渗流数值计算方法3在线视频

同学们好

下面我们继续学习

“饱和稳定渗流边界条件的处理”

我们接着讨论第(3)种边界条件

自由水面的处理

如前所述自由水面也称为浸润面

在浸润面上应该同时满足两个条件

1）孔隙水压力为零

测管水头等于位置水头也即h=z

2）法向流速为零亦即vn=0

自由水面的位置通常无法事先确定

需在计算中通过迭代的方法计算确定

自由水面边界的迭代计算

是渗流计算中的难点问题

在饱和稳定渗流的有限元计算中

目前常采用的

迭代计算自由水面的方法主要有4种

包括，一 网格修正法

二 单元传导矩阵修正法

三 剩余流量法

和四 初始流量法

下面逐一进行介绍

第一种方法为网格修正法

网格修正法是一种最原始的方法

如图中所示采用该种方法时

第1步是要对全域

都划分有限元计算网格进行计算

第2步根据计算结果

用自由水面的第1个条件

测管水头等于位置水头

也就是h=z

确定出一个近似的自由水面

第3步

对h=z这个近似自由水面之下的区域

重新划分有限元计算网格并进行计算

第4步

再次用h=z确定一个近似的自由水面

如此反复直到前后两次

所确定的近似自由水面差别不大为止

显然网格修正法的计算效率很低

每次迭代均需人工重新划分

有限元计算网格

所需的人工工作量也非常巨大

第二种方法为单元传导矩阵修正法

这种方法是由著名学者巴特（Bathe)

在1979年提出的

如图中所示

该法和网格修正法的前两步是相同的

也即第1步对全域划分

有限元计算网格进行计算

第2步

用h=z确定出一个近似的自由水面

该自由水面将计算域

分成上下两个区域

其中自由水面上面的区域

是不应该包括在计算中的

在下面的迭代计算中

是应该拿掉的部分

但是巴特认为

没有必要重新划分计算网格

在下面的计算中

仍然可以采用整体的

计算网格进行计算

但是对自由水面上下的两个区域

可分别采取不同的渗透系数进行计算

所以对下面的第3步

巴特提出应改变渗透系数重新计算

具体方法为

对近似自由水面下面的区域

仍采用原土体

真正的渗透系数k进行计算

而对近似自由水面上面的区域

要将原渗透系数除以1000后进行计算

将渗透系数缩小1000倍进行计算

实际上是降低了近似自由水面

以上区域对总体渗透矩阵的贡献

相当于“在数值”上去掉了这部分的网格

这和在划分计算网格时

直接拿掉这些单元所起的作用是相同

但是由于不再需要重新划分

有限元网格

所以可以大大地减少人工的工作量

第4步

再次用h=z确定一个近似的自由水面

如此反复直到前后两次

确定的近似自由水面差别不大为止

巴特提出的单元传导矩阵修正法

只使用了1个有限元计算网格

所以可称为固定网格法

但是该法在每次迭代计算时

仍然需要进行

整体渗透矩阵的集成和分解

计算工作量仍然较大

第三种方法为剩余流量法

这种方法是由著名学者德赛（Desai)

于1976年提出的

剩余流量法也是固定网格法

按全域剖分计算网格

并按照整体网格进行计算

如图中所示

在求解出整体域内水头h的分布后

可根据h=z确定出

相应近似自由水面的位置AD

对近似自由水面AD

可计算其上的法向流速vn

如果在AD上均有vn=0

则此时在AD上将会满足

自由水面的全部两项条件

说明此时将会是问题的真解

在进行迭代计算的初期

上述情况是不可能发生的

也即在近似自由水面AD上

一般总有法向流速vn≠0

实际上该法向流速的大小

描述了上部ABCD区域

对下部计算域影响的大小

Desai建议

可通过在AD上

叠加一反向的法向流速-vn的方法

消除上部ABCD区域

对整体计算所造成的影响

具体方法为

首先计算该反向法向流速

所产生的结点流量

将计算所得到的结点流量

叠加到整体渗流方程的右边项

重新进行求解

可得到相对更加精确的近似解

如此反复进行迭代计算

直到得到满意的计算结果为止

从上述原理和计算过程可见

剩余流量法是固定网格法

在迭代过程中不再需要进行

整体渗透矩阵的集成和分解

计算效率得到了大大的提高

但是剩余流量法需要计算法向流速

在2维和3维的情况下

分别需要进行曲线和曲面的积分

计算过程相对较为繁琐

第四种方法为初始流量法

这是德国亚琛工业大学的

韦德卡(Wittke)教授等在1984年提出的方法

韦德卡教授是我在德国读博时的博士生导师

初始流量法同样是固定网格法

按全域剖分计算网格

并按照整体网格进行计算

如图中所示

在求解出整体域内水头h的分布后

同样可根据h=z确定出

第一次迭代对应的自由水面的位置

此时位于该自由面以上区域中的结点

都会有一个负的压力水头

也即这些结点的总水头小于

其相应的高程

然而这样的计算结果还存在着误差

因为该区域内的单元

对整体的渗流方程组还存在有影响

对于该区域中的某个单元e'而言

其对渗流方程的影响

就是该单元的有限元方程 Q8e'=[Ce']·he'

可见其影响相当于产生了结点流量Q8e'

因此可通过在方程的右边项

施加一个和它相反的结点流量

以校正它的影响

对所有这种的单元的结点流量

进行求和

可得到总体的结点流量的修正列阵

将其从渗流方程组的右边项中减去

并重新进行求解

可得到一个新的水头分布

此时所得到的计算结果

较前一次结果会有所改进

所确定的自由水面的位置

相对也会更加的准确

如此重复进行迭代计算

直到得到满意的计算结果为止

我们知道在有限元的数值计算中

在进行单元积分时

一般是通过对高斯点积分来进行的

因此对于和自由水面相交叉的单元

也即对仅部分位于渗流区域中的单元

可分别按照具体的高斯点

进行h

在数值积分流量修正列阵时

仅仅需要积分和累计具有负压力水头

也就是位于自由水面以上的

那些高斯点的数值

相反对于在自由水面下部的高斯点

则要忽略不计

这张图图示了

对一个具有8个高斯点的立方体单元

在计算其结点流量修正值时

如何用这种方法考虑自由水面的位置

对于图中的这个单元

只有1、2、3、4和7号高斯点

位于自由水面以上

在数值积分结点流量修正列阵时

也仅需要考虑这些高斯点的数值

图中还给出了具体的计算流程图

在某一的迭代步

对所有的单元进行循环

其中对每个单元

再对单元所有的高斯点循环

计算高斯点上的水头hj

并判断是否hj>zj

是则说明该高斯点位于自由水面以下

不需要计算

否则说明该高斯点位于自由水面以上

需计算该高斯点的结点流量

并累计到总修正列阵中

对所有的单元循环完成后

在渗流方程组的右边项中

减去所得到的总修正列阵

并重新进行求解

可得到一个新的水头分布

进行误差比较

满足计算精度要求后

自由水面迭代完成

最后我们讨论第(4)种边界条件

渗出面边界的处理

如前所述在渗出面边界上

也应该同时满足两个条件

1）孔隙水压力为零

测管水头等于位置水头

也就是h=z

2）vn≤0

也即渗透流速指向渗流域的外部

也就是说要有流量从域内渗出到域外

对于一个实际的渗流问题

渗出面的大致位置常常是已知的

或者可以根据日常经验

大体进行判断

但其精确的范围

却是需要在计算中通过迭代确定

例如对图中所示的隧道渗流

隧道周围的边壁都是可能的渗出面

再如对图中的

土石坝及岸坡渗流问题

土石坝下游边坡

两岸岸坡和河谷底部

都是可能的渗透面

渗出面区域的精确确定

需要在求解过程中

采用迭代计算来实现

图中以边坡渗流计算为例

给出了确定渗出面位置的迭代方法

在进行第一次迭代计算前

首先要给出一个假设的渗出面的区域

要足够大使得实际的渗出面区域

一定要包含在该假设区域之内

在第一次迭代计算时

首先假设所给该区域中的所有结点

均是渗出点

取h=z并按照给定水头边界进行计算

在求解方程组后

再对区域中所有的结点进行逐点检验

判断是否真实的渗出点

判断的依据

是渗出面边界的第2项条件

法向流速vn是否≤0

也即其结点流量

是否是由域内渗出域外

如果对某个可能的渗出点

计算所得到的结点流量小于零

例如图中的结点②

则说明该点可满足渗出面的

全部两项边界条件

该点就是真的渗出点

反之如果计算所得到的结点流量

大于零

如图中的结点①

则说明对于该点

如果要使其水头达到h=z

则必须向该结点补充水量才行

这表明该点是位于

高于实际渗出面位置的地方

它不是一个真的渗出点

而是一个假的渗出点

在下一次的迭代时

需要对这样的假渗出点进行修正

为了提高计算效率

应尽量通过修正有限元方程右边项

并保持渗透矩阵不变的方法

来进行迭代计算

具体方法可采用类似确定

自由水面位置的结点流量法进行修正

例如对于图中的结点①

经判断它是一个假渗出点

这时h=z这个假设就不能成立了

可将其边界条件改写为

hji+1=χ·hji

其中χ<1 如可取χ=0.95

修正了结点1的水头之后

需要相应对有限元方程的右边项

进行修正

修正由于将结点水头

由hji改为hji+1

所对应的方程右边项的改变

重新进行计算

并按上述的步骤进行迭代

直到所有渗出点都满足相应的条件

由此可见这种方法

与确定自由面位置的迭代方法很类似

实际上渗出面的位置

对自由水面的位置也是有影响的

因此在实际的有限元计算中

渗出面迭代和自由水面迭代

是需要相互嵌套交替进行的