古树包滑坡坡面径流和饱和-非饱

古树包滑坡坡面径流和饱和-非饱

和渗流有限元分析

汪海洋1 谢新生1 岳桃丽2

(1.四川大学水利水电学院，四川 成都6100651；2广东省罗定市水务局，广东 罗定 527200）

摘要：本文阐述了坡面径流和降雨入渗关系，对滑坡坡面径流和非饱和渗流进

行了耦合分析。利用有限元的方法模拟古树包滑坡降雨入渗过程中、分析了降雨强度、初始含水率等因素对坡面径流产流的影响规律。 关键词:径流 饱和非饱和渗流 有限元 降雨入渗

1工程概况

古树包滑坡体位于清江水布垭工程区域坝区外台子上滑坡体下游侧的清江左岸岸坡上，距水布垭大坝1.6公里，在水布垭水利枢纽工程可行性研究阶段未进行1/1000的地勘工作。2000年4月坡体开始发生拉裂变形，随着雨季的来临，变形逐渐加剧，局部解体严重。随后的人工活动包括为贯通1#公路而开挖的较大型的人工高边坡、坡体之上的移民活动加剧、坡体之上地表水活动强烈，坡体进一步产生强烈变形，裂缝逐渐贯通至450m 高程，2001年1月7日，古树包巨厚的滑坡体终于滑坡化，开始全面解体。至2001年6月，坡体变形已发展到5000m 高程。

2渗流问题有限单元法求解原理

根据建立线性代数方程组的依据不同，求解渗流问题的有限单元法主要有变分法和伽辽金法，边界元法等。变分法将求解微分方程的问题用变分问题(即寻找泛函的极小值问题)来代替；Galinjin 伽辽金法[1]从问题的基本方程出发，结合线性剖分插值与加权剩余原理进行渗流分析，在计算中不依赖于泛函或变分原理，可用来求解方程已知、但泛函未知或没有泛函可采用的问题。伽辽金法的计算精度较高。当存在相应的泛函时，伽辽金法与变分法将得出相同的结果。因此，伽辽金法应用较广，下面介绍伽辽金法求解渗流问题。

引入权函数i W ，由伽辽金法建立饱和/非饱和、非恒定渗流有限单元法计 算方程(不考虑内源，不考虑水体和岩土体的压缩性)，即

()()()0w h

h

h h x x y y

z

z

p t u u u i

S dxdydz k k k W γ

∂

∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂Ω⎡⎤+++=⎢⎥⎣

⎦

⎰⎰⎰ （1）

对式（1）的前3项进行分部积分，经整理后得到

0=--⎪⎭

⎫ ⎝⎛++∂∂Ω∂∂∂∂∂∂∂∂∂∂∂∂∂∂Ω⎰⎰⎰⎰⎰⎰q dxdydz S dxdydz W W W W k i t h p z z h y y h x x h u w i i i γ （2） Γ=

∂∂Γ

⎰d q n

h

i u

W

k （3）

式中：q 为边界上给定的流量。 将计算域剖分为单元，式（3）化为

(

){}0i

i

i

w

e e

h h h h x

x

y

y

z

z

p

t

u

i e

e

W W W

S

dxdydz dxdydz q k

W γ

∂

∂

∂

∂

∂∂∂∂∂∂∂∂∂∂∂∂Ω

Ω

+

+--=⎡⎤⎣⎦

∑∑⎰⎰⎰⎰⎰⎰ （4）

设N

j

为形函数，则 ⎪⎪

⎪⎪

⎭

⎪⎪⎪

⎪⎬

⎫====

∑∑∑∑j j j j j j j

j j j

j

j

h N z N y N x

N h z y x （5） 令N W i i =，则式（4）变为

{}e

e

i j j

j

z j

z

j

y

j y

j

x

j x

e

u e

q d d Jd h S d d Jd N N N N N N N N h

h

h

k j w j j

j j

j j

∑∑⎰

⎰⎰

∑∑∑⎰⎰⎰∑=-⎪⎪⎭

⎫ ⎝

⎛+

+

∂---∂∂∂∂∂

∂∂

∂∂

∂∂

∂∂

---][11111

1

11111

1

ζηξγζ

ηξ （6） 可得如下代数方程组

{}{}{}{}e

e

e

e t

h e e

e

q B h A ∑∑=-∂∂][][ （7）

其中：

ζηξd d DJd N N k A i

j

e

u ij ⎰

⎰⎰---=1

1111

1

（8）

因此，用有限单元法求解三维饱和/非饱和、非恒定渗流的计算式为

[]{}[]{}{}q B h A t h =-∂∂ （9）

用向后差分求解式（9）设第m 时步的解已求得，则第1+m 时步的解为

[]

[]{}

[]{}{}

1

1

1

1

11

1)(++∆-++∆+-=

-

m m

m t m m t m q h B h B A m

m

（10）

首先求出21+m 时步的系数矩阵[]

2

1+

m A 、[]

2

1

+m B ，它们决定于

{}

2

1

+m h ，而

{}1+

m h 是未知的，故需要迭代求解。由 m-1时步m 时步的h 值外插可求得{}

2

1

+m h ，即