渗流数值计算的有限单元法

渗流数值计算的有限单元法

渗流问题常用的数值计算方法主要的是有限差分法和有限单元法，其中有限差分法出现较早，随着计算机和计算技术的发展，有限单元法在这一领域的应用日益广泛，并在计算复杂渗流工程问题中占有较大优势，下面简要介绍渗流问题有限单元法的基本概念。 （1）控制方程和边界条件

本章介绍的渗流仅限于饱和土中的渗流，且假定渗流过程中土的孔隙比不变，即土的渗透系数不随时间变化。前面已推导出二维渗流问题的控制方程为

02222=∂∂+∂∂y

h

k x h k y x （3-64） 渗流问题数值计算的边界条件有两类。第一类边界条件是给定水头边界，这种边界常出现在渗流区域与地表水的连接处。对于这种边界上的所有点，每一时刻水头h 是给定的，即

),,(),,(1t y x t y x h ϕ=Γ，1,Γ∈y x ，0>t （3-65）

式中：h －边界1Γ上某点),(y x 在t 时刻的给定水头；ϕ－已知函数。

第二类边界条件是给定水流通量（流入或流出）边界，在这种类型的边界上，单位面积流入（或流出）的通量是已知的，即

),,(),cos(),cos(2

t y x q y n y h k x n x h k w

y w x =⎥⎦⎤⎢⎣⎡∂∂+∂∂Γ

，2,Γ∈y x ，0>t （3-66）

式中：),cos(x n ,),cos(y n —边界外法线向量与坐标轴正向之间夹角的余弦；

),,(t y x q —t 时刻边界2Γ上某点),(y x 处水流量，为已知函数。

除了上述两类边界条件外，渗流问题的边界条件也可以是混合边界条件，即部分边界上的水头为已知、部分边界上的流量为已知。 （2）泛函和变分

式（3-64）所示微分方程在复杂的边界条件下无法得到解析解，采用数值方法计算时，首先建立h 的泛函，一定边值问题的解就是这个泛函的极小值，这个求解过程就是变分。

对二维渗流情况（图3.18），在x 方向，t d 时间内，外力在单位重量流体上所做的功的增量为

*

-=x

x x h q dA d d （3-67） 其中，x q d 为x 方向的流量增量；*x h d 为在x 方向上的近似水头差，上标*表示近似，*

x h d 可

以表示为

x x

h h x

d d ∂∂=*

* （3-68）

图3.18 单元流体做的功

则

x x

h q A x x d d d ∂∂-=*

（3-69）

由y x h k q x x d ∂∂-=*可得y

k q x h x x d =∂∂-*

，代入式（3-69），整理后得 x x x x q q y

k x

A d d d d =

（3-70） 到时间0t 外力所做的总功为

2

d d d d d 2

0x

x Q x x x x Q y k x q q y k x A x =

=⎰ （3-71） x Q 是在某时间0t 内，水头为h 时的总渗流量

y x

h

k Q x

x d ∂∂-= （3-72） 则

y x x

h k A x x d d )(22

∂∂=

（3-73） 单位体积外力所做的功

2

)(2x

h k a x x ∂∂=

（3-74） 由于外力做功等于土体内存储的能量，设渗流的能量密度为x ω、y ω，则

2

)(2x

h k a x x x ∂∂-

=-=ω （3-75a ）

2

)(2y

h k a y y y ∂∂-=-=ω （3-75b ）

同样，在某一渗流域Ω中，忽略流体的可压缩性，其渗流能的表达式为

y x y h k x h k h I y x d d )()(21)(22⎥⎦⎤

⎢⎣⎡∂∂+∂∂=⎰⎰Ω

（3-76）

对于非稳定渗流，存在自由水面的情况，边界上能量为⎰⎰

Γ∂∂=ΓΓd cos d 2

h t

h

qh θμ

，则上述渗流能为

⎰⎰⎰Ω

ΓΓ∂∂+⎥⎦⎤⎢⎣⎡∂∂+∂∂=2

d cos d d )()(21)(22h t h

y x y h k x h k h I z x θμ （3-77）

)(h I 是一个泛函，求其极小值，对应的),(y x h 就是式（3-64）的解。

（3）渗流场的离散与插值函数

将渗流场划分成有限个单元（常用三角形单元），其中任意一个三角形单元ijm 内任一点的水头),(y x h 可通过线性插值求得，即

y a x a a z x h 321),(++= （3-78）

式（3-78）中的三个常数1a 、2a 、3a 可由该单元三个节点i 、j 、m 的坐标及相应的水头i h 、j h 、m h 来表示，则此单元内任意点的水头可用矩阵形式表示为（右上角e 指单元而言）

{}e

m j i m j

i h N h h h N N N h ][][=⎪⎭

⎪⎬⎫

⎪⎩⎪⎨⎧= （3-79）

式中N 只是节点坐标x 、y 的函数，称为形函数，其表达式为

e i i i i A y c x b a N 2/)(++= （3-80a ）

e j j j j A y c x b a N 2/)(++= （3-80b ） e m m m m A y c x b a N 2/)(++= （3-80c ）

式中的系数a 、b 、c 和三角形的面积e

A 都是节点坐标x 、y 的函数，其表达式为

j m m j i y x y x a -= （3-81a ）