4.4.1 饱和渗流数值计算方法1在线视频

同学们好

下面我们学习第4节

“饱和渗流数值计算方法”

这一节主要包括

渗流分析的方法

饱和稳定渗流的有限元方法

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

以及饱和非稳定渗流有限元方法

共4个部分

首先来看渗流分析的方法

在工程实践中

对渗流问题通常可采用

如下四类方法进行求解

第一 数学解析法或近似解析法

就是采用数学的方法

求取相应边界条件下

拉普拉斯方程的解析解

或者近似解析解

通常适用于一些渗流域相对规则

和边界条件简单的渗流问题

二 流网法

流网法是一种图解法

具有简便迅速的优点

并能应用于渗流场边界

较为复杂的情况

因而该法在工程上

得到了比较广泛的应用

三 电比拟试验法

利用渗流场与电场

都满足拉普拉斯方程

采用电流场来模拟渗流场

该法简便直观

可以用于二维和三维的渗流问题

四 计算机数值解法

随着计算机

和数值计算技术的迅速发展

以有限单元法为代表的各种数值方法

在渗流问题的模拟计算中

得到了越来越广泛的应用

数值解法

不仅可用于各种二维或三维问题

也可很好地处理各种复杂的边界条件

近年来已成为了

求解渗流问题的主要方法

本节下面详细介绍有限元方法

有限单元方法是目前

最常用的数值求解方法

对某种特定的问题

有多种的途径建立相应的有限元方程

第1种 对于力学问题

可以直接利用虚功原理

建立有限元方程

第2种变分原理

也是建立有限元方程的常用方法

此时对于力学问题

可以利用最小位能或最小余能原理

而对于其它问题

则需要首先寻找

相应的微分方程的泛函

针对该泛函再应用变分原理

建立相应的有限元方程

第3种加权余量法

也是建立有限元方程的

一种普遍和常用的方法

该法可用于变分的泛函

找不到或者根本不存在的情况

本节介绍如何利用泛函和变分

来建立饱和稳定渗流

有限元方程的方法

为此下面首先给大家回顾一下

泛函和变分的基本概念

泛函通俗的讲就是函数的函数

而求泛函的极值问题称为变分

历史上有许多有关泛函和变分的

有名的极值问题

这里为了说明泛函与变分的概念

我们举其中的一个例子

这就是最速降线问题

如图中所示

A点和B点不处在同一条的铅垂线上

我们将坐标原点设在A点

假设在AB之间

连接着某种光滑的曲线轨道y=f(x)

有一个小球在重力作用下

由A点沿该曲线自由下滑

如果不考虑摩擦力的作用

显然小球从A下滑到B的时间

会随曲线形状的不同而各有不同

现在我们的问题是

从A下滑到B

所需时间最短的曲线为哪一条

这就是历史上有名的最速降线命题

他是伯努利在1696年

以公开信的形式提出来的

曾引起广泛的注意

而经许多著名科学家如莱布尼兹

牛顿等大约100年的努力

才获得较为完善的解答

今天我们可以通过高等数学中

有关泛函与变分的方法

很方便的求解这个问题

如图中所示根据能量守恒

可以计算出小球

达到某点时的速度v=√2gy

然后可以计算出

滑过微弧长ds所需的时间dt

再将dt沿整个弧长积分

可得到小球从A点下滑到B点

所需的总时间T为

T等于从0到xB积分

√[(1+y'^2)/(2gy)]·dx

这里时间T是一个标量

从这个式子可以看出

标量T是函数y的一种广义函数

我们称之为泛函

那么这里的最速降线问题

实质上也就变成了

求取这样一个泛函的极值问题

我们称之为变分

下面我们简单讨论求泛函的极值

也就是变分问题

与求解一个微分方程的关系

我们知道求解一个泛函的极值

一般总能通过欧拉方程

将其转换为求解一个

微分方程的边值问题

但是反过来

如果我们要求解一个微分方程

是否也可以将其转换为

求解一个泛函的极值问题呢

答案是这对于许多问题

但并不是所有的情况都是可以的

这时需要寻找该微分方程

所对应的泛函是什么

如果找到了相应的泛函

则求解该微分方程的边值问题

也就可以用求解

其相应泛函的极值问题来等价

但可惜并不是所有的微分方程

都能找到相对应的泛函

求解泛函的极值问题

对于有限元方法是非常合适的

是建立有限元方程的有力的工具

非常方便

图中是饱和稳定渗流的运动方程

和常用的4类边界条件

所以我们现在求解的

是这个偏微分方程的边值问题

如果我们想把它转换为

求解一个泛函的极值问题

就需要首先找到

渗流微分方程相对应的泛函

如果我们想把它转换为

求解一个泛函的极值问题

就需要首先找到

渗流微分方程相对应的泛函

很幸运有学者找到了相应的泛函

这个泛函的具体表达式如图中所示

这个泛函包含了渗流的运动方程

和我们前面讨论过的

第2类流速边界条件

也就是说由渗流运动方程

和第2类边界条件所决定的渗流场

等价于让这个泛函取极值

所对应的渗流场

其它的边界条件未包括在这个泛函中

需要在进行有限元求解时

强制得到满足

这样我们就把

求解一个偏微分方程的边值问题

转换成了一个求解泛函的极值问题

而求解泛函的极值问题

可用有限元法方便地来进行求解

为了简化讨论

下面具体以图中的

二维情况为例进行讨论

由变分原理建立一个有限元方程

一般可按照下面的标准的步骤来进行

一 离散计算域

将计算域划分为有限个的单元

二 建立单元的描述

对渗流问题

可取结点水头为基本未知量

通过单元插值函数建立单元的描述

包括求取各种的导数等

三 求取单元上的泛函值

根据单元上的描述

可以计算各个单元上的泛函值

实际上是将单元上的泛函

表示成为了一个结点水头的多元函数

四

泛函对单元求和得到域上总的泛函值

每个单元的泛函值就是

单元结点水头的多元函数

求和之后总的泛函值

也就是域内所有结点水头的多元函数

五 求泛函的极值

这时求泛函的极值

就变成了一个求多元函数的极值问题

这个多元函数的未知量

就是域内所有的结点水头

而多元函数求极值的条件

就是将其对每一个未知量求导

也就是对每一个结点水头求导

并让导数等于0

这样就可以得到一个

以结点水头为未知量的方程组

六 边界条件的处理

对没有包括在这个泛函中的

边界条件进行处理

让它们强制得到满足

七 求解方程组

得到域内所有结点上的水头

同学们好

下面我们学习建立渗流的

有限元方程的具体过程

建立有限元方程的第1步是离散计算域

将计算域划分为有限个的单元

然后第2步是建立单元的描述

以图中的单元e为例

设单元e的结点分别为i、j、k等

结点水头分别hi、hj、hk等

单元形函数分别为Ni、Nj、Nk等

则单元内任一点的水头

可用形函数表示为图中的(1)式

式中[N]矩阵为单元的形函数矩阵

对(1)式求导数并取负号

可得单元的水力坡降得到(2)式

其中[B]矩阵为单元形函数的导数矩阵

将水力坡降再乘以渗透系数矩阵[k]

可得单元的渗透流速得到(3)式

建立了单元的描述之后

第3步是求取单元上的泛函值

计算的公式为(4)式

显然在建立了单元的描述之后

公式(4)中的各项都是可以计算的

且由积分计算所得到的

单元e子域上的泛函Ie(h)

是单元结点水头的多元函数

下面是第4步

对整个域上所有单元求和

可得域上的总泛函I(h)

由公式(5)表示

公式(5)是对所有的单元求和

其中N为单元总数

第5步求总泛函I(h)的极值

显然前面积分得到的总泛函I(h)

应是域内所有结点水头的多元函数

多元函数取极值的条件

就是将其对每一个结点水头求导

并让导数等于零

也就是∂I/∂hl=0 l=1,…,M

式中M为结点总数

将前面的公式(4)和公式(5)

代入到 上式

可得图中的公式(6)

公式(6)中的各项都可以由前面

所建立的单元描述得到

图中给出了一些项的具体示例如

单元内一点水头he(x,y)的表达式

水头的一阶导数∂he/∂x的表达式

水头的一阶导数

再对水头求偏导的公式等

把这些式子代入公式(6)

并对单元e进行单独整理

可得单元e的泛函Ie

对某个结点水头

hi导数的计算公式 公式(7)

这个式子虽然看起来很长

但各个项都很有规律

同样可得Ie

对其它结点水头导数的计算公式

将单元e上的这些计算式

按单元e的上的结点水头

hi、hj等排列成列

组成一个行列式

可得到图中的式(8)

公式(8)实际上就是

单元e的有限元方程

其中[C]e为单元渗透矩阵

相当于结构计算中的单元刚度矩阵

它的计算公式如图中所示

为[B]矩阵的转置乘渗透系数矩阵[k]

再乘上[B]矩阵

然后在单元e上进行积分

Qe为单元结点流量列阵

相当于结构计算中的单元结点力列阵

它的计算公式如图中所示

包含有两项

第1项为由单元内源项q

引起的结点流量

在结构计算中

相当于由体积力引起的结点力

第2项为由单元边界上

法向流速渗入引起的结点流量

在结构计算中相当于

边界上表面力所引起的结点荷载

将得到的单元e的有限元方程

对渗流域内所有的单元

进行集成求和

对于渗流域中的全部结点

可得到如图中所示的

总的有限元方程组

这是一个方程数目

和域内结点总数M相同的线性方程组

其中[C]称为渗透矩阵

相当于结构计算中的刚度矩阵

如果域内总共有m个结点

则[C]为m×m的方阵

h为结点水头列阵

Q为右边项列阵

相当于结构计算中的

右边项结点力列阵

从本质上看

我们所得到的这个渗流的

有限元方程组是基于渗流的

连续性条件推导所得到的

因此每个方程对应的物理意义

是相应结点上的渗透流量的平衡

每个方程的右边项对应的是

相应结点上的结点流量

下面我们将渗流有限元方程

与应力变形有限元方程进行对比

尽管所涉及的物理问题和微分方程

具有很大的不同

两个有限元方程在许多方面

却具有相似性

我们对比下面的4个方面

一 基本未知量和单元插值

在应力变形分析中

基本未知量为结点位移

对三维问题每个结点有3个未知量

单元插值采用形函数[N]

在渗流分析中

基本未知量为结点水头

每个结点1个未知量

单元插值也采用形函数[N]

从未知量数目看

对相同有限元计算网格的计算分析

渗流分析的计算量相对要小很多

二 1阶导数

在应力变形分析中

基本未知量的一阶导数为单元应变

可写成乘以[B]矩阵的形式

如果再乘以弹性矩阵[D]

可得到单元的应力

在渗流分析中

基本未知量的一阶导数

为单元的水力坡降

也可写成乘以[B]矩阵的形式

如果再乘以渗透系数矩阵[k]

可得到单元的渗透流速

三 方程右边项

在应力变形分析中

我们知道每一个方程的右边项

都是对应有具体的物理含义的

那就是在相应结点某个位移方向上

所作用的结点力的大小

在渗流分析中

每个方程的右边项

也同样是对应有具体的物理含义的

那就是相应结点上的渗透流量

四 有限元方程的物理意义

在应力变形分析中

每一个方程也都对应有

具体的物理含义

那就是在相应结点位移方向上的

结点力的平衡

在渗流分析中

每一个方程也同样都是对应了

具体的物理含义

那就是相应结点上的渗透流量的平衡