裂隙岩体渗流场与应力场耦合分析方法

第12卷　第2期

1997年　6月云南农业大学学报JournalofYunnanAgriculturalUniversityVol.12No.2Jun.1997

裂隙岩体渗流场与应力场耦合分析方法*

高海鹰

(河海大学结构中心,南京210024)

摘要　从工程观点出发,忽略不计完整岩块的渗透系数,借助于层面缝隙流运动规律(立方

定律)和变分基本原理,建立了三维裂隙岩体渗流场和应力场耦合模型.并推出了耦合的有

限元支配方程,提出了运用同一套单元网络对两场进行区域离散的观点.以重力坝坝基渗流

场与应力场耦合计算为例,验证了模型的可行性.

关键词　耦合;立方定律;裂隙岩体;变分原理

中图分类号　TV139.14

近年来裂隙岩体渗流与应力耦合分析方法主要有两种:一是以应力为直接耦合对象,

1〕耦合关系采用Louis〔方法,此法只表明垂直应力对与之正交方向上的渗透性的影响,未

能反映实际的复杂应力状态对渗透特性的影响.

2〕另一种方法是以应变为直接耦合对象,如Erichsen〔给出了平面问题且两组正交的

裂隙的算例.较第一种方法的经验关系更符合实际.但Erichsen中未能按实际的由多组裂隙构成复杂的裂隙网路进行三维分析.本文忽略完整岩块的渗透性,认为水流只在岩体的裂隙网络中流动,可由立方定律给出裂隙的变形和渗透系数之间的关系,耦合的机理非常清楚,在此基础上建立三维裂隙岩体渗流与应力耦合分析方法.

1　渗流场和位移场(应力场)的耦合变分原理

当仅进行渗流场分析时,可以找到对应的泛函

I(h)=8{〔kx()2+ky()2+kz()2〕+hU0}d822+p2h-Udp+p3(Kh+ch)dp2

令DI(h)=0,可得到渗流场支配方程.

当仅要求位移场时,可以取势能泛函∫∫

8∫(1)Pp(u)={A(e)-∫ijfiui}dS-u∫sRii(2)

令DPp(u)=0,可得到位移场的支配方程,进而可求应力场.

在正常荷载作用下,岩体绝大部分变形发生在裂隙.裂隙的过水能力又和裂隙隙宽的3次方成比例.因此,渗流场受应力环境的影响很大,而渗流场的改变将改变渗透压力

*1996—06—15收稿

　138云南农业大学学报第12卷的分布,后者又必将对应力场产生影响.因此,泛函Ⅰ的宗量是水头h,参变量是位移u.泛函Pp的宗量是位移u,参变量是h.待求的水头函数h和位移函数u必须使得方程

(1)(2)泛函同时取得驻值,可使下式成立.

DI=0(3)DPp=0

　　方程(3)式中的两个变分方程是两个含参变量的方程,按位移求解的有限元支配方,可得到最终的耦合支配方程组:

〔A(u)〕{h}={U}(4)　　(5)〔k〕{u}={R(h)}

其中{R(h)}={R0}+{RX},{R0}是由除去渗透压力以外的外荷载形成的等效结点外荷载列阵,{RX}是由渗透压力产生的等效结点渗透荷载列阵.〔A(u)〕为总体水力传导矩阵,{h}是结点水头列阵,{U}为等效列点流量列阵,〔K〕是结构的劲度矩阵,{u}是结点位移列阵.

2　耦合支配方程的求解

耦合方程的求解实质上是以耦合本构关系为基础,迭代求解非线性方程组(4)和(5).

2.1　耦合本构关系　裂隙水压以面力形式作用在岩块上,通过裂隙的变形,即裂隙开度的变化来反映两场的耦合关系,裂隙的变形和渗透系数的关系,由立方定律给出:

　　K=12C2(6)其中e为裂隙宽度,$e为裂隙张开度的改变量,可由裂隙单元各结点位移求得.k为渗透系数,C为动力粘滞系数,g为重力加速度.

2.2　求解三维裂隙渗流场　首先将耦合区域离散化,为了耦合的需要,分析渗流场和应力场必须采用同一套有限元网络,而这一套网络是基于应力场分析的需要建立的.裂隙单元在结构分析中采用空间八结点夹层单元,在渗流分析中采用三维空间中的平面四结点等参数单元模拟,这样通过相应的结点对压缩就可得到渗流分析有限元网络信息.

本文忽略岩块透水性,认为裂隙网络中各局部裂隙上的二维流动构成整体的三维流动.由于各裂隙单元不在同一平面上,需对每个裂隙单元建立局部坐标系o'x'y',在o'x'y'坐标系中将裂隙单元上的水流视为各向同性流,类似于二维连续介质渗流分析方法,可得到裂隙单元的传导矩阵中的元素.

　　Kij=Kfõee∫∫-111(-1iijjT　)(　)ûJûd§dG(7)

其中e为裂隙宽度,Kf为渗透系数,可由立方定律求得.

将各裂隙单元的传导矩阵总装成整体传导矩阵[A],可得三维裂隙网络渗流有限代数方程组,如公式(4)所示.根据渗流边界条件求解方程(4)即可得网络各结点水头

第2期

值.高海鹰:裂隙岩体渗流场与应力场耦合分析方法139　

裂隙面上某一点受到的渗透压力为P

4

P=CX(h-z)=CX∑Ni(hi-Zi)

i=1(8)

将P进行面积分可求得等效结点荷载.

2.3　位移(应力)场的非线性有限元分析

本文将裂隙单元采用空间八结点夹层单元模拟,岩块用空间八结点块体单元模拟,夹层取为不抗拉的理想弹塑性模型,岩块和坝体部分取为弹性材料.在分析弹塑性材料本构关系时,采用塑性流动理论,屈服准则取用D.C.Drucker条件.用“变Kp迭代法”求解有限元支配方程.详见文献〔3〕和〔4〕.

2.4　耦合方程组迭代步骤

(1)计算岩体初始位移(应力)场Li,Rij各裂隙的隙宽bi.

(2)由立方定律求出裂隙单元渗透系数,并求得单元传导矩阵,总装形成总体传导矩阵.

(3)由渗流边界条件求解渗流场,可得到各结点的水头和渗透压力分布,进而可求得等效结点渗透荷载列阵{RX}.

(4)将{RX}加上其他荷载形成的等效结点外荷载列阵{Ro}作非线性结构位移、应力分析.

(5)根据位移,应力状态确定各裂隙单元新的隙宽和渗透系数.

(6)判断迭代过程是否收敛,收敛了则输出结果,若没有收敛,则返回第(1)步,继续迭代.　收敛判据:前后两次迭代得到的各结点水头(或位移)保持不变或非常接近.　　　

附表　材料力学参数

Tab.　MechanicsofMaterialsparameters

裂隙

种类

水平裂隙

竖直裂隙

岩块

坝体

帷幕弹横E(mpa)50522.0×1042.4×1042.3×104泊松比L0.350.350.20.1670.167容重C(T/m3)1.91.92.62.42.4内摩擦角U15°15°凝聚力c(mpa)0.10.1(0)(0)(0)

3　算例

图1为某重力坝坝基渗流与应力耦合分析的计算域,坝高210m,底宽150m,坝段长20m,上游水位为200m,下游水位为0m,帷幕厚度为15m,深90m,渗透系数为5×10-6m/s.帷幕下游设排水孔,直径为10cm,间距为5m,孔深75m,排水孔出

　140云南农业大学学报第12卷的两组裂隙面,考虑到计算机容量的限制,隙宽分别取为3.0mm和2.5mm,间距取为15m.材料力学参数见附表

坝体渗透性与地基渗透性相比甚小,可不参加渗流分析.水的运动粘滞系数为1.5×10-6m2/s.渗流边界条件为:上游水位200m,下游水位0m,其余边界为不透水边界.整个结构受到的荷载有自重、水压力和渗透压力.本文共作了以下二组耦合分析;

(1)无帷幕和无排水设施

(2)既有帷幕又有排水设施

图2和图3为渗流场的等值线图,其中(a)为耦合前等势线,(b)为耦合后等势线.耦合后等势线均偏向下游,上游区等势线较疏,坝底及下游区等势线较密.这是由于建坝蓄水后坝基应力分布不均匀,不均匀的应力场使裂隙网络渗透性发生变化,上游处竖直裂隙有张开趋势,坝下水平裂隙有闭合趋势,因而上游处裂隙的渗透系数相对地比坝下裂隙的渗透系数大得多,上游水头损失小,下游水头损失大,等势线编向下游.图4为坝底扬压力分布,实线是渗流与应力耦合前结果,虚线为耦合后结果.从中可见考虑应力场影响下坝基扬压力值将增大

.

…… 此处隐藏：3508字，全部文档内容请下载后查看。喜欢就下载吧 ……