（水工结构工程专业论文）土坝渗流三维有限元分析及其支持向量机监测模型研究.pdf

摘要 土石坝在渗流作用下往往会产生渗漏、渗透破坏、坝坡失稳等不利的事故， 直接关系到工程的经济效益与安全问题。因此，对土石坝的渗流计算和渗流监测 问题进行研究具有十分重要的意义。 土石坝的渗流计算是一个边界条件复杂的三维问题，通常采用的水力学方法 及平面闷题假设，因与实际状况有较大的差异而使计算结果与实际相差较大。本 文在介绍渗流计算的基础理论及其有限元理论的基础上，编制了三维有限元渗流 计算程序，可以更真实地模拟渗流边界条件。算例验证表明本程序具有良好的适 用性，而且实例计算的结果与实测值吻合良好。 土石坝的渗流监测在保障工程安全运行中起着重要作用，其中一个重要方面 是对测压管水位( 即浸润线位置) 的监测。本文在介绍测压管水位统计模型相关 理论和方法的基础上，首先建立了测压管水位的统计监测模型。此后，针对统计 模型的不足之处，引入了支持向量机方法的相关理论，建立了基于支持向量机方 法的测压管水位监测预报模型，与统计模型相比，结果拟合精度更高，预测趋势 更加符合实际情况。 关键词：渗流计算有限元分析支持向量机渗流监测统计模型 a b s t r a c t u n d e rt h ea c t i o no fs e e p a g e ，e a r t hd a ma l w a y so c c u rs o m ea d v e r s ea c c i d e n ts u c h a sl e a k a g e ，s e e p a g ef a i l u r e , s l o p ed e s t a b i l i z a t i o na n ds oo n , a n di td i r e c t l yi n f l u e n c e s e c o n o m i cb e n e f i ta n ds a f e t yo fp r o j e c t s oi ti sv e r yi m p o r t a n t1 0r e s e a r c hi ns e e p a g e c a l c u l a t i o na n dm o n i t o r i n go f e a r t hd a m t h es e e p a g ec a l c u l a t i o no f e a r t hd a mi sa3 dp r o b l e mw i t hc o m p l i c a 【e ib o u n d a r y c o n d i t i o n , a n di tu s u a l l yc a u s e ss o m ed i f f e r e n c et or e a l i t yb ya d o p t i n gt h eh y d r a u l i c m e t h o da n dh y p o t h e s i so fp l a n ep r o b l e m i nt h i sp a p e r , t h em o r ea c c u r a t eb o u n d a r y c o n d i t i o ni ss i m u l a w ab y c o m p i l i n g 3 df i n i t ed e m e n t s e e p a g e c a l c u l a t i o n p r o g r a m m i n g t h er e s u l to ft h es p e c i f i ce x a m p l ep r o v e st h a ti ti saf e a s i b l ew a yt o s o l v es e e p a g ep r o b l e m s e e p a g em o n i t o r i n go f e a r t hd a mp l a y sv e r yi m p o r t a n tp a r ti ns a f e t yo p e r a t i o no f e n g i n e e r i n g o n eo ft h ek e ya s p e c t so fs e e p a g em o n i t o r i n gi st h ep h r c a t i cl i n e ( p i e z o m e t r i cl e v e om o n i t o r i n g a c c o r d i n gt ot h es t a t i s t i cm o d e lo fp i e z o m e t r i cl e v e l ， m o d e lo fp i e z o m e t r i cl e v e lm o n i t o r i n go nt h eb a s i so fs u p p o r tv e c t o rm a c h i n e si s f o u n d e d c o m p a r i n gw i t ht r a d i t i o n a ls t a t i s t i cm o d e l t h en e wm o d e li sm o r ea c c u r a t e a n db u sb e t t e re x t e n s i o n f o r e c a s t i n gc a p a b i l i t y k e yw o r d ：s e e p a g ec a l c u l a t i o n ，f i n i t ee l e m e n ta n a l y s i s ，s u p p o r tv e c t o r m a c h i n e s ，s e e p a g em o n i t o r i n g ，s t a t i s t i cm o d e l 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果。除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得盘洼盘堂或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：砖晚采 签字日期：炒谬年，二月矽日 学位论文版权使用授权书 本学位论文作者完全了解鑫鲞盘堂有关保留、使用学位论文的规定。 特授权鑫鲞盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：砖吮禾 签字日期：，年阻月口日 导师签名： 签字日期：) 酗孵，上月却日 第一章绪论 1 1引言 第一章绪论 土石坝是最古老的坝型之一。几千年来，在全球所建的众多的水坝中，大多 数都是土石坝。特别是我国，在建国后修建的8 6 万多座水坝中绝大多数也都是 土石坝。坝工一般修建在河道或者渠道上。由于坝工的修建，建筑物的上游常常 形成壅水，并使建筑物的上下游产生水位差，形成水头，在这一水头作用下，建 筑物及其地基中将产生从水位高的一方面向水位低的一方面的渗流。渗流对于建 筑物，特别是土质建筑物及其地基会产生许多不利的影响，因此渗透稳定作为土 石坝工程中的主要问题，直接关系到大坝的经济效益与安全问题。世界大坝破坏 实例调查结果表明，由于渗透稳定造成的破坏与失事，占整个事故的4 0 以上， 例如美国提堂( t e t o a ) 土坝的失事就是由于渗流造成的。可见。对土石坝的渗 流计算和渗流监测问题进行研究具有十分重要的意义。 渗流计算的任务可归纳如下： 1 计算通过建筑物( 如土质堤坝) 及其地基的渗流量，以确定建筑物上游 水库的渗漏损失水量，分析其对水库蓄水的影响，以及是否需要采取其他防渗措 施，以减少渗漏损失。 2 计算土石堤坝坝体内的浸润线，以便分析土石堤坝坝坡的稳定性。 3 计算建筑物( 如土质堤坝) 及其地基内各点处的水头和压力，以确定在 建筑物及其地基内的分布和变化，用以预测产生渗透变形的可能性。 4 计算作用在建筑物底面上的扬压力，以便分析建筑物的稳定性及应采取 的防渗、排渗措施。 5 。计算建筑物及其地基内各点处的渗透流速和水力坡降，特别是渗流出逸 处的渗透流速和水力坡降，以便分析建筑物及其地基的渗透稳定性，所应采取的 防渗、排渗措施。 土石坝的渗流监测是运行监控中不可缺少的重要组成部分。由于土石坝的材 料特性以及库水位的频繁变化，坝内部渗流状况比较复杂。直接关系到前后坝坡 的稳定性及坝体的渗透稳定性，一直为大坝管理人员所关心，在坝体和坝基的适 当部位设置一定数量的测压管，并进行观测，可以了解坝体的浸润线和各点渗透 压力的大小，并对测压管水头进行预报，这对坝体的稳定分析和渗流分析，以及 为运行决策的提供依据等都有重要的意义。 第一章绪论 1 2 渗流计算方法概述 渗流力学是- f l 边缘学科。1 8 5 6 年，法国工程师达西( h d a r c y ) 通过实验 提出了线形渗透理论，为渗流理论的发展奠定了基础。1 8 8 9 年，h e 茹可夫斯 基首先推导了渗流的微分方程。此后。许多科学家和地下水动力学研究人员对渗 流力学及其解析解法迸行了广泛和深入的研究，并取得了一系列的研究成果。 1 9 2 2 年，h h 巴甫洛夫斯基正式提出了求解渗流场的电模拟法，为解决比 较复杂的渗流问题提供了一套有效的工具。为解决更复杂的渗流问题，逐步发展 和研究了电网络法。 随着电子计算机的迅速发展，数值方法即有限差分法、有限单元法和边界元 法在渗流分析中得到了越来越广泛的应用。有限差分法是在1 9 1 0 年由l f r i c h a r d s o n 首先提出来的。经过长期的研究和广泛的应用，目前己具有较完善 的理论基础和实践经验。有限单元法的基本思想是在1 9 4 3 年由r c o u r a n t 提出 的，1 9 6 5 年0 c z i e n k i e w i c z 和y k c h e u n g 提出有限单元法适用于所有按 变分形式计算的场问题，为该方法在渗流分析中的应用提供了理论基础。边界元 法是建立在经典力学理论b e t t i 互换定理和f r e d h o l m 积分方程基础上的，初见 于2 0 世纪6 0 年代后期，当时被称为边界积分方程法( b o u n d a r yi n t e g r a l e q u a t i o nm e t h o d ) 。 我国有关科学工作者和工程技术人员结合大规模的坝工建设对渗流理论、 数值方法和实验技术等都进行了广泛深入的研究，并取得了丰富的研究成果。 综上所述，渗流计算的方法归纳起来有解析法、数值分析法和实验分析法。 解析法是指利用有关数学手段直接求解基本微分方程的方法，包括直接解 法、复变函数法、组合法和水力学法。直接解法是采用直接解渗流基本方程来计 算渗流要素的一种方法；复变函数法是利用复变函数保角变换理论，将实际的渗 流问题变换为一个已有解答的区域，从而使问题获得解答；组合法是将复杂的渗 流区域划分为几个简单的渗流区段来进行计算的一种方法；水力学法是建立在对 渗流条件作某些简化假定基础上的一种方法，计算比较简单，能用于计算各种实 际渗流问题，但这种方法只能得出渗流场中某一截面上的平均渗流要素，而不能 计算出渗流场任意一点处的渗流要素，同时由于其基本假定与实际情况有一定出 入，所以计算存在一定误差。 数值方法目前主要有有限差分法、有限元法及边界元法。 有限差分法是把导数近似的用差商代替，边界条件、初始条件也相应的作类 似的代替，最后便把定解问题的求解转化为一组代数方程组的求解问题。有限差 分法的优点是其原理易懂，算式简单，有较成熟的理论基础。其缺点是往往局限 于差分网格，对曲线边界和渗透介质的各向异性模拟比较困难。 2 第一章绪论 有限单元法的基础是用有限个单元的集合体来代替渗流区。这些单元的结合 点称为结点。然后选择简单的函数来近似的表示每个单元上的水头分布。最后的 解将得出渗流区中离散点处的近似水头。即首先分割所研究的区域，然后用比较 简单的函数来构造每个子区域中的水头表达式，最后集合起来形成线性代数方程 组，从而得到原来的渗流场的解，而不是一下解出有关整个渗流区的问题。由于 建立线性代数方程组的依据不同，有限单元法本身又分为若干种方法，如里兹法、 迦辽金法、均衡法等等。有限单儿法在模拟曲线边界和各向异性的渗透介质方面 比有限差分法具有较大的灵活性。 边界元法是通过将求解域边界剖分为若干个单元，化边界积分为线性代数方 程组来求其数值解。建立所研究问题的边界积分方程是边界元法的基础。边界元 法与有限单元法相比，它便于处理无限和半无限渗透介质、渗流奇异( 如排水井 点) 和自由面等问题，而且由于它只对研究域的边界进行剖分，其数据信息显著 减少，一般说来其精度也高于有限单元法。其缺点是它所建立的系数矩阵是满阵， 而且不对称，即使结点较少，也将占据可观的内存。同时，它对三维非均质问题 的应用尚存在相当大的困难。 目前常用的实验方法有：砂槽( 土槽) 模拟法、粘滞流模型法( 缝隙水槽法) 、 水力网模型法、电模拟法和电阻网模型法等。 1 3 渗流分析的回归模型及支持向量机方法的引入 由于大坝渗流是一个三维问题，影响因素比较复杂，边界条件难以确定，致 使根据计算或实验而采用的防渗、导渗措施不可能是很完善的。实践证明，原型 观测是监控渗流状态最真实的行之有效的方法，因而被大多数工程所采用。国际 大坝委员会第2 3 号会刊的报告中指出“坝基扬压力和渗流量的观测是最直接的 也是最有意义的安全措旌。如果是重力坝，这些观测是头等重要的”。 土石坝的测压管水位是反映土石坝渗流状况的重要因素，配以巡视观测，可 以找出测压管水位变化的趋势规律，并根据新测值，评价判断大坝的渗透状况。 目前，测压管水位统计模型主要采用三种模型：统计模型、确定性模型及混合模 型。统计模型由于其计算方便，实用性较强，一直被工程界广泛运用。库水位、 降雨和坝后水位是影响测压管水位较明显的因素。随着工程运行时间的增长，土 颗粒材料逐步固结，使得土的渗透性能发生变化，这说明还存在着时间效应，测 压管水位变化还具有滞后性，它滞后于库水位的变化，前期库水位会影响到当前 的测压管水位，需根据回归及实际情况选择前期影响天数。 本文拟将机器学习领域中的支持向量机方法和理论应用到渗流的监测和预 测，尝试建立基予支持向量机方法的渗流监控方法和理论。支持向量机方法是建 第一章绪论 立在统计学习理论的v c 维理论和结构风险最小原理基础上的，根据有限的样本 信息，在模型的复杂性( 即对特定训练样本的学习精度) 和学习能力( 即无错造 识别任意样本的能力) 之甸寻求最佳折衷，以期获得最好的推广能力。支持向量 机方法的几个主要优点是： 1 它是专门针对有限样本情况的，其目标是得到现有信息下的最优解而不仅仅 是样本数趋于无穷时的最优值； 2 算法最终转化成为一个最优化问题； 3 算法将实际问题通过非线性变换转换到高维的特征空间，在高维空间中构造 线性判别函数，其特殊性质能保证机器有较好的推广能力，同时巧妙地解决 了维数问题，其算法复杂度与样本维数无关。 1 。4 本文的主要工作 渗流计算及渗流监测预报在土石坝设计、施工和管理中发挥着重要作用，也 是评价工程的社会效益、经济效益和环境效益的重要内容。解决工程中的诸多渗 流问题必须依据渗流的基本理论和方法，通过对渗流场中各个物理量的变化规律 的分析和研究，才能为工程实践提供足够的理论和监测依据。 本文根据土石坝的特点，首先将三维有限元分析引进到渗流计算当中，对引 滦入津工程中的尔王庄水库坝体的渗流要素进行了计算和分析，然后结合该水库 的渗流监测数据，引进了支持向量枫方法，对数据进行了回归，并得到其监测预 报模型，经检验，具有较高的精度，可以将之应用于水库安全信息系统的渗流监 测预报当中。 本文主要完成了以下几方面的工作： 1 介绍了渗流计算的基本理论及其有限单元法原理； 2 开发了一套适用对象较为广泛的三维有限元渗流计算程序，并对之进行 了算例验证。利用该程序，可确定在特定库水位下坝体的浸润线、等势 线、流线等渗流物理量； 3 介绍了渗流监测的基本理论及其多项式统计模型，并求得该模型； 4 介绍了支持向量机的基本理论及其回归方法，编制了基于m a t l a b 的支 持向量机回归程序，并运用之 寻到了渗流的支持向量机监测预报模型。 4 第二章渗流分析的基础理论和有限单元法原理 第二章渗流计算的基础理论和有限单元法原理 渗流计算的主要内容在本文的绪论里已做了陈述，在此不再赘述。在实际工 程设计或管理中，渗流计算的主要问题是稳定渗流，即在渗流场中任点的流速 大小和方向都不随时间变化的水流。本章主要介绍稳定渗流计算的基本理论及其 相应的有限单元法原理。 2 1 渗流计算的基础理论 2 1 1 渗流基本概念和基本定律 水或其他流体在多孔介质中运动，由于多孔介质中孔隙、孔隙大小，形状很 复杂，因此水质点在其中的运动毫无规律，有些地方甚至不连续，所以就不能像 研究地表水一样直接研究水质点的运动，而常用平均概念和综合性的参数来研究 其运动规律。这种方法的实质就是用和真实水流属于同一流体的、充满整个含水 层( 包括全部的孔隙空间和土体岩石颗粒所占据的空间) 的假想水流来代替仅仅 在孔隙中运动的真实水流，通过对这一假想水流的研究来达到了解真实水流平均 渗透规律的目的。这种假想水流同时还应具有下列性质： ( 1 它通过过水断面的流量与真实水流通过该断面的流量相同； ( 2 ) 它在断面上的水头及压力与真实水流的水头及压力相等： ( 3 ) 它在多孔介质中运动时所受到的阻力等于真实水流所受到的阻力。 满足这些条件的假想水流称为渗流。本文把渗流所占据的空间称为渗流场， 描述渗流的参数称为渗流要素，如压力p 、速度v 及水头h 等等；在计算时，把 它们看作空间坐标x , y ，z 和时间r 的连续函数。 2 1 1 1 渗流的几个基本概念： 1 、渗流作用力 ( 1 ) 静水压力与浮力 在饱和的多孔介质中，渗流对某一接触面上的静水压力符合流体的静水压力 分布，即任一点上的静水压力为p = 日( 日为水头) 。 固体颗粒淹没于水中，由于静水压力作用而产生浮力，使颗粒的重量减轻。 同样对于一定体积的多孔介质，只要孔隙彼此连通并全部充满水时，由于各点静 水压力的存在，使多孔介质整体也将受到浮力，且等于各颗粒所受浮力的累加总 和。根据土的三相比例关系得浮重度，。为： 第二章渗流分析的基础理论和有限单元法原理 ，= ，。一h ( 2 1 ) 其中，。、分别为土的饱和重度、水的重度。 ( 2 ) 渗透力 在饱和的多孔介质中，地下水在孔隙中的运动，对颗粒骨架的稳定性将发生 破坏作用。由动水压力作用在土体上产生的力称为渗透力。单位体积多孔介质沿 流线方向所受的渗透力为： z 。一h 筹：j ( 2 2 ) 以一h 面2 l 2 一圳 式中歹水力梯度； 以渗流路径长度，m ； 羽钇上的水头损失，m 。 浮力和渗透力的作用直接关系着土体的渗透稳定性，对土体的渗透变形研究 有着重要的意义。虽然静水压力产生的浮力不直接破坏土体，但能使土体的有效 重量减轻，降低其抵抗破坏的能力，因而也是一种消极破坏力。动水压力产生的 渗透力和渗流冲刷力，则是一种主动的破坏力，它与渗透破坏的程度成直接的比 例关系。 2 、渗流水头和水力梯度 地下水流中任意一点的总水头为 日：+ 三+ ! 二( 2 3 ) 凡2 9 式中p 研究点上的动水压力，k p a ； 凡水的重度，k n 一； ：研究点的位置高度，即研究点到任意选定基准面的垂直距离，m 流速水头，m 。 z g 由于自然界中地下水的运动速度很缓慢，流速水头 很小，可以忽略不计。 z g 因此在渗流计算中可以认为总水头日等于测压管水头王，即 以= ：+ 万p ( 2 - 4 ) 渗漉不是理想流体，具有粘滞性。因丽在运动过程中能量要不断消耗，反映 在水头上为沿程不断减小。因而在渗流场中各点的水头并不都是相同的。可标识 6 第二章渗流分析的基础理论和有限单元法原理 为h = h ( x ，y ，z ，f ) ，它构成一个标量场。由场论可知，标量场可构成个梯度场。 在渗流力学中，把太小等于梯度值，方向沿着等水头面的法线指向水头降低的方 向的矢量称为水力梯度，用表示，即 j = - j 塑万 ( 2 5 ) 一月 、一d ， 厕 式中万法线方向的单位矢量。 矢量，在直角坐标系中的三个分量为： 以= 一警；以= 一百o h ；以= 一警 ( 2 - 6 ) 3 、渗流量 通过渗流计算求得所研究渗流区域中各点的渗透速度后，就可以选择适当的 渗流控制断面计算通过相应区域的渗流量q ： q _ - v 。钌 ( 2 7 ) 式中s 过流断面，m 2 ； 叱过流断面的法向渗透速度，m s 。 上述流速、流量、水头、水力梯度以及渗透力等物理量是渗流计算中必不可 少的要素。 2 1 1 2 渗流的基本定律d a r c y 定律： 法国工程师d a r c y 通过渗透实验得到了渗流量q 与横截面4 、水头差 一h 2 ) 成正比，- q 渗透路径三成反比的规律，即： q = k - a - h 1 - h 2 ( 2 8 ) 或v ：里：七星! = 墨= k - ， 4上 ( 2 9 ) 上式称为d a r c y 定律。它指出了渗透速度v 与水力梯度，或渗透阻力的线性关 系，故又称为线性渗透定律。 d a r c y 定律描述的是层流状态下渗透速度与水头损失关系的规律，即渗透速度 v 与水力梯度，成线性关系只适用于层流范围。所鞋它的应用受到一定水力条件的 限制。 许多研究者都曾指出，随着渗透速度的增大，d a r c y 定律即渗透速度与水力梯 度j 之间的线性关系便不再成立。由于水沿着弯弯曲曲的途径运动，并且在不断 地改变它的运动速度、加速度和流动方向。这种变动有时是很剧烈的，受惯性力 的影响，使水流运动不再符合达西定律。当水流运动速度较小时，这种惯性力的 7 第二章渗流分析的基础理论和有限单元法原理 影响是不大的，粘滞力占优势，水流运动符合达西定律。随着速度的增加，惯性 力也相应的增大了。当惯性力占优势时，由于惯性力与速度的平方成芷比，达西 定律就不再适用了。 下列两种情况可认为超出d a r c y 定律的适用范围。 一种情况是在粒径为纯砾以上的介质中的渗流，如堆石体中的渗流，当水力 梯度较大时，流态已不再是层流而是紊流。这时d a r c y 定律不再适用。 另一种情况是在粘性很强的介质中。由于初始水力梯度的存在，使得当梯度， 小于初始水力梯度山时，没有渗流发生。只有当j 山时才会发生渗流。初始水 力梯度以的存在将会对渗流场以及介质的周结产生影响。 由于目前渗流计算的数学模型大都是以d a r c y 定律为基础建立起来的，这显 然对于许多实际情况带有近似的性质。因此，在设计中实际运用渗流计算资料时， 应对此予以适当考虑。 渗透系数k ，也称为水力传导系数，是一个重要的水力地质参数。它是表征 多孔介质输运流体能力的标量。在数值上为当水力梯度j = l 时的渗透速度。渗透 系数不仅仅取决于介质的性质( 如粒度成分、颗粒排列、裂隙性质和发育程度等) 。 而且和渗透流体的物理性质( 重度、粘滞性等) 有关。渗透系数一般可在实验室 或野外现场通过实验测定。 2 1 2 渗流基本方程 2 1 2 1 渗流连续性方程 渗流连续性方程是从质量守恒原理出发来建立的。在渗流场中取出一个单元 体，见图2 - i ，设各边长分别为缸，缈，缸，并和坐标轴平行。设渗流沿坐标轴方 向的分速度分别为匕，v ，叱，流体的密度为p ，则单位时间内通过垂直于坐标轴方 向的单位面积的水体质量分别为，：。单位时间内沿：f 方向流入单元体的 质量为以缈z ，流出单元体的质量为( p + 昙，血) 4 y a = 出，因此，沿j 方向流 入和流出单元体的质量差为： 8 第二章渗流分析的基础理论和有限单元法原理 z x 图2 一l 缪沉场中的单兀俸 哪y a z a t o 叩l 冬咿戡帑k 缸= 一! 譬缸弩缸心一1 0 ) 同理，沿，方向和z 方向流入和流出单元体的质量差分别为一! 墼a x a y a z a t 和一旦璺妙厶2 & ，因此，f 时间内流入和流出单元体总的质量差为： 一警a x a y a z a t 一孥舭f 一孥6 x a y a z a t 盘鲫出 根据质量守恒原理，它应等于单元体内流体质量随时间的变化量 詈l 册缸妙心k 。所以有 凹 一l 警a r a y a z a t + 警缸舭+ 孥o z 蚴& 斗昙b 舭k l 办卯l 甜 上式称为渗流的连续性方程。 若流体密度p 为常数且多孔介质不可压缩，则上式变为： 冬+ 冬+ 誓：0 。( z - 1 2 ) 该式表明在同时间内流入单元体的水的体积等于流出的水体积，即体积守 恒，此时把渗流当成刚性流体，即为不可压缩流体在刚性介质中流动的连续性方 程。当地下水为稳定渗流时就可得到上述结果。 9 第二章渗流分析的基础理论和有限单元法原理 2 1 2 2 渗流基本微分方程 根据d a r c y 定律在各向异性介质中有： = 一罢；哆= 也警；屹= 一也警( 2 - 1 3 ) 通过对( 2 1 1 ) 式进行一系列交换及根据d a r c y 定律有： 砜a ( k 。a 苏m j + ) 号一等) + 鲁i 警) = 以警 c z 州， 其中 1 t 。储水率或单位贮存量，其值表示对于单位体积多孔介质，当水 头降低一个单位时，有多孔介质压缩及水的膨胀所释放出来的 水量。 对于稳定渗流，由于掣：0 ，所以方程为： 丢卜警) + 参卜等 + 丢i 警) = 。 c z ， 上述方程就是稳定渗流基本微分方程。 2 1 3 渗流基本微分方程的定解条件 前面描述渗流的偏微分方程具有多解性，为了能从它们全部的解中选出一个 满足某个具体问题的确定解，就必须加上一些附加条件，这些附加条件就是通常 所说的定解条件。定解条件包括边界条件和初始条件。在稳定渗流中由于渗流状 态不随时间变化，可仅考虑边界条件。 2 1 3 1 边界条件 边界条件指渗流区域几何边界上的水力性质。又可分为第一类边界条件和第 二类边界条件。 ( 1 ) 第一类边界条件，又称为给定水头边界。 当渗流区域的某一部分边界( 如r i ) 上的水头为已知时，边界条件为： t t ( x ，y ，z ，f ) i r = 妒 ，y , z ，) ( 2 - 1 6 ) 应当注意，给定水头边界不是定水头边界，两者要分开。定水头边界是指在 边界上的水头函数日或势函数伊是不随时间变化的，是个常数。 ( 2 ) 第二类边界条件，又称为给定流量边界。 当渗流区域的某一部分边界( 如r ) 上的法向流速己知时，边界条件为： 1 0 第二章渗流分析的基础理论和有限单元法原理 七剿：g o 小毛f ) ( 2 - 1 7 ) i l 式中万厂。的外法线方向5 q ，的侧向补给量 最常见的这类边界为隔水边界，此时q = 0 。在介质各向同性的条件下，上式 简化为： 剖n - o 泌 ( 3 ) 自由面边界条件 在自由面( n ) 上有 h = 剿：o 锄lr | h ( x ，y ，z ，r ) lr | = z ( x ，y ，r ) 嚣1 日 2 1 2 2 塑：o 功 ( 2 - 1 9 ) ( 2 - 2 0 ) 在a b 上 在d c 上 电a e 上 ( 2 2 1 ) 在e d 上 桃、b c 上 确定士坝渗流的边界条件可参见图2 2 。- 图2 - 2 渗流区域和边界示意 2 1 3 2 初始条件 通常第一类边界条件( 即流场中的水头分布) ，它在开始时刻t = 0 时对整个流 场起支配作用。所以进行非稳定渗流计算时，必须先求得开始时刻流场的水头分 布作为初始条件。即 h ( x ，y ，z ，o l ，卸= 日y ，z ) ( 2 2 2 ) 第二章渗流分析的基础理论和有限单元法原理 初始条件可以根据需要，任意选择某一时刻作为初始条件。 2 2 渗流计算的有限单元法原理 2 2 1有限单元法渗流计算的基本思想 有限单元法是与数字电子计算机同时代的产物。因此，尽管这种方法和以前 的数值方法有很多共同的地方，但还是有它自己的一些特点，即利用了高速电子 计算机所提供的特殊方便。特别是它能系统的编成计算机程序，比较方便的处理 复杂的边界条件、非均质和厚度变化的岩土层、多个含水层组成的含水组等麻烦 和困难的问题。有限单元法对第二、第三类边界不必作专门处理，能自动满足， 因而便于处理复杂的边界条件。单元大小随便，形状也有多种，三角形、四边形、 四面体、六面体等，可以视需要选择，有较大的灵活性。 有限单元法是采用“分块逼近”手段来求解偏微分方程的一种数值方法。它 是把渗流区域划分成许多小的相互联系的亚区域，成为“单元”。这些单元的形状 不同，一维的单元为线段，二维的单元为三角形、四边形，还有三维的为多面体 单元等等。分割所研究的区域后，用比较简单的函数来构造每个子区域中的水头 表示式，最后集合起来形成线性代数方程组，从而得到原来渗流区域的解。 渗流问题的有限元分析大致步骤如下： ( 1 ) 把待求解区域划分为一系列数目的有限个单元，单元的顶点称为结点。 单元和单元之间通过结点相联系。这个过程称为“离散化”或“区域剖分”。如图 2 - 3 所示： 图2 - 3 渗流区域的离散化示意 ( 2 ) 有限单元内的待定近似函数( 单元水头分布函数) 由已知的若干插值函 数迭加组成。在多个单元所共有的结点上，其水头函数应相等，而且要保证在单 元交界面上水头函数也相等： ( 3 ) 建立每个结点的单元系数矩阵，也称单元渗透矩阵； 第二章渗流分析的基础理论和有限单元法原理 ( 4 ) 把单元渗透矩阵组合起来，形成一组描述整个渗流区域的代数方程组， 建立总的渗透矩阵： ( 5 ) 把给定的边界条件也归并到总矩阵中； ( 6 ) 求解线性代数方程组的解，最终得到问题的解答。 根据建立线性代数方程组的依据的不同，有限单元法又可分为四种方法： ( 1 ) 直接法：来源于结构分析的刚度矩阵法，一般只能用来处理简单的问题。 ( 2 ) 变分法：它的根据是把变分泛函区域离散化，然后对待定的结点场函数 值变分，使变分泛函达到极值，它可以用来处理比较复杂的问题。 ( 3 ) 加权余量法( 迦辽金法) ：它从问题的基本方程出发，在分析中并不依赖 于泛函和变分原理。其优点是可以用来处理方程已知，但泛函未知或者 没有泛函可采用的问题。因此，它是用途更为广泛的方法。 ( 4 ) 能量平衡法：它通常指某一系统的机械能或热能的平衡，它与迦辽金法 一样，不需要有现成的泛函。 对于处理渗流问题主要是后三种方法，而变分法更为常用。本文采用变分法。 2 2 2 稳定渗流问题的变分原理 所谓渗流问题的变分原理，就是把描述渗流场的偏微分方程的边值问题的求 解转化为求某个泛函的极值的问题。 描述符合d a r c y 定律，非均质各向异性土体的稳定渗流问题的偏微分方程的 边值问题为： 昙( 屯罢 + 昙( b 詈 + 丢i 警) + 形x , y , z e c 2 h ( x , y ，z ) | r l = 烈x ，y ，z ) 在j r 。上 ( 2 2 3 ) k o 饿h e o 嘲加屯掣o yc 0 如川+ 屯警毗力= g 札上 珊 。 也 一一 式中 日0 ，儿力水头函数，m 。 t ，k y ，t 当坐标轴与渗透主轴方向一致时，墨y ，z 方向的渗透系数， m s ； 一一入渗或蒸发水量，即单位时间内，单位体积上进出含水层的水 量，m v s ； q 渗流区域，即为r 1 和r 2 所围成的研究区域，其中r 1 为给定水头 边界，r ：为给定流量边界； 盯边界上单位面积上的流入流出的流量，m 3 t 月- 边界的外法线方向。 第二章渗流分析的基础理论和有限单元法原理 对于所研究的渗流场，根据变分原理，则渗流场上述边值问题等价于泛函 j ( 王，) 的极值问题，可表达如下： ，= j 1 嚣。 k 。o 缸h ，2 + k 。o 砂h ，2 + 吒( 罢) 2 一：船 一l l q 溜睁2 4 ) r 根据研究区域的结构特征，进行计算区域的离散化，即： q = q f ( 2 2 5 ) 某单元内的水头插值函数为； 日儿力= m ( f ，r ，f ) 岛 ( 2 - 2 6 ) 式中，m g ，r ， ) 为单元的形态函数；h ，为结点水头值；g ，玑a 为基本单元 的局部坐标。 由于渗流场整个区域已离散成单元，故( 2 - 2 4 ) 式可表达为： i ( h ) * 1 8 ( 奶= r a i n ( 2 2 7 ) 由上述可见，在一定的单元划分情况下，泛函，( 上d 完全确定结点水头，因此， 式( 2 - 2 7 ) 的极值条件即为： 罢孕：塑盟：0 ，川，2 叫 ( 2 - 2 8 ) 为h 1o h t 。 对各子区域迭加，可得由有限元求解渗流场的方程组： 【k 胆 = q ) ( 2 2 9 ) 式中 医】总体渗透矩阵，由单元渗透矩阵集合而成； 幻) 自由项向量，为域内源和已知边界流量对结点的贡献； 沮 待求结点水头向量。 在解微分方程的边值问题时，第二、第三类边界条件必须作为定解条件列出， 而在求解相应的变分问题时，它们会自动满足，无须作为定解条件列出。由于变 分原理的解能自动满足，所以基于变分原理的有限单元法就能使这个有利因素得 到发挥。对于第一类边界条件，必须作为定解条件作专门处理。 由于渗流微分方程定解问题的求解可以转化为求解某一泛函的极值的问题。 因此，对于复杂的一时难于求解的渗流微分方程定解问题可以不去直接求它的解， 而是通过求相应的泛函的极值来求它的近似解。 1 4 第二章渗流分析的基础理论和有限单元法原理 2 2 3 等参单元 在有限单元方法中，随着单元结点数目的增加，插值函数的方次随着增加， 因而用于实际问题的分析时可能达到的精度也随之提高了。对于一个给定问题的 求解域，预期用较少的单元即可获得需要精度的解答。但是用较少的形状规则的 单元离散几何形状比较复杂的求解域常会遇到困难，因此需要寻找适当的方法将 规则形状的单元转化为相应的其边界为曲线或曲面的单元。在有限单元方法中最 普遍采用的变换方法是等参变换，即单元几何形状的变换和单元内的场函数采用 相同数目的结点参数及相同的插值函数进行变换。采用等参变换的单元称之为等 参单元。借助于等参单元可以对一般的任意几何形状的工程问题和物理问题方便 地进行有限单元离散。 由于等参变换的采用使等参单元的刚度、质量及荷载等特性矩阵的计算仍在 规则的域内进行。因此不管各个积分形式的矩阵的被积函数如何复杂，都可以方 便的采用标准化的数值积分方法进行计算。用来离散三维空间的单元形状有许多 种，最简单的是四面体和六面体。采用四面体单元由于使分割后的空间各单元之 间相互关系复杂，一般很少单独使用。最常用的还是六面体单元，这种单元由于 分割简单，各单元位置之间相互关系清楚，公式推导简单，因而得到广泛的应用。 本文采用的即为8 结点等参单元。如图2 - 4 所示。 ：锣 整体坐标e 局部坐标e 图2 q 空间等参数单元 取形函数为： 奶= 1 ( 1 + 茧) ( 1 + ，7 玎，) ( 1 + 为) f = 1 ，2 ，8 ( 2 - 3 0 ) 其中 ，矾，白) 为相应结点f 在g 川，力上的坐标。此时各变量的变换公式为： 888 工= z x , 。；，) ，= e y , 。以，z = z l 。f ( 2 3 1 ) 则在整体坐标中日的插值函数为 霞：皂h ? n i i = l ( 2 - 3 2 ) 第二章渗流分析的基础理论和有限单元法原理 芘幽表达式燹为 ，c 册姐面= 三删t 。瞧县警卜b 陲e 。警 2w 瞎日警) 2 2 w c 毒e 。f 蝴一m 鼢 ( 2 3 3 ) 一e 。f 幽啮一m 鼢 ( 2 一 对于每一个单元有 等。喜研珥警警+ k ， a 钞n ia 砂n j 一+ 鬈警警k j 驴。m 叠吻仡 ，= 1 ，2 ，8 ( 2 3 4 ) 令 办= 名j f 。( l k ；, a 苏n , 动缸v j 。+ 联y w - 哪面- + k ；, i a n , i a n a r l 撒 e j = 1 1 1 w e n j a x a y a z 则得： 研( 费) 8 h i a i ；l m 鲫； o z f ( r ) a h ； 且 岛 岛 ( 2 - 3 5 ) 利用参数坐标系和原坐标系体积元素之间的关系： 西叻沈= i ，i 倒叩右 ( 2 3 6 ) 其中i 卅表示雅各比行列式的绝对值。再通过高斯积分即可求出以和e ，。 对于式( 2 - 3 3 ) 中的最后一项l g 胁，如果曲面是包含单元e 的一个面，则 参数方程为： 4 x - 善；n g ，7 ) 扭l 4 y = e y ；n ( 善，7 ) i = l 4 z = ：；q ，7 ) f = l 1 6 ( 2 3 7 ) 研雕磁 儿 靠如 九屯妃 第二章渗流分析的基础理论和有限单元法原理 1 - i 在曲面上的插值公式为： 膏：4 日；0 n g ，7 ) 于是 鼍= 蠢廷q 叠a s = 面紫翠4k 玉= 妒一 其中凼表示q 的面积微元，它可以利用下列公式计算 l l s ：压瓦予d 弘q 式中 e = ( 势盼( 毒 2 g = ( 黔附+ 蝌 f ：鱼鱼+ 堂鱼+ 鱼鱼 o 0 1 7d ；d q d o q 冬上述各式代入计算券的表达式，并以妨表示之，则有 鹞；l l q n a i a s = u - l q n n | 厨d 承珥3 = 1 ，2 3 4 于是可得任意单元ei - 的单元渗透矩阵为 甜；( 日) 8 h ； 研( 膏) 8 h ： 研( 青) 勰； = 巨睦 蜀 岛 最 翻 跌 珐 ， 级 - 0 ( 2 4 5 ) 由极小值条件把各单元的微商叠加，并使其和为零，即可得有关 玩，日2 ，日。( m 为未知结点总数) 的代数方程组 - 肛 = 护 ( 2 4 6 ) 阻】为总渗透矩阵。求解此方程组即可得本定解问题的近似解。 8 9 o l 2 3 4 3 3 4 4 4 4 4 一 一 一 一 一 一 一 2 2 2 2 2 2 2 第二章渗流分析的基础理论和有限单元法原理 2 2 4 无压渗流场自由面搜索技术 用有限单元法求解渗流场时，如何简捷高效地处理并确定自由面位置，同时 保证迭代过程收敛的稳定性是该问题的关键所在。多年来，从变网格法有限元数 值分析到固定网格的各种算法中，反映出该问题至今仍未得到完全解决。 无压渗流场自由面需要迭代确定，故属于边界非线性问题。采用有限元法求 解无压渗流场时，通常有变网格法和固定网格法两种。变网格法直接引用6 0 年代 中期0 c z i e n k i e w i c z 等人最初求解简单有压流问题时的方法。此法是将自由面 作为可动边界，在每一步迭代计算中需重新形成渗透矩阵。它对初始渗流自由面 位置的要求较高。如果初始位置与最终自由面位置相差很大，则极易造成单元严 重畸变，影响计算精度。当自由面附近渗透介质为非均质，特别是有水平分层的 不同渗流介质时，程序处理十分困难。另外在渗流场与应力场耦合分析中，应用 受到一定限制。因此在无压渗流场分析中，越来越倾向于采用固定网格法。 固定网格法即保持有限元网格不变的前提下求解渗流场。这些算法的核心是 在计算过程中网格不变。n e u m a n 于1 9 7 3 年提出求解自由面的g a l e r k i n 方法。d e s a i 于1 9 7 6 年提出剩余流量法，并于1 9 8 3 年发展到用于计算非稳定渗流。1 9 7 3 年 b a i o c c h i 基于变分不等式概念提出了新算法，后来o d e n 及b r u c h 又作了进一步的 发展。b a t h e 、s t e v e n 、a r a l 、l a i n 以及g e l l 等人在这方面都有所贡献。李春华提 出用固定网格法分析有自由面的二维渗流，并补充推导了公式，给出了算例。在 诸多的算法中，比较有影响的有单元渗透矩阵调整法、虚单元法、剩余流量法和 初流量法等。 2 2 4 1 单元传导矩阵调整法口” 利用对渗流场有限元计算的结果，根据单元永头与结点位置势的比较。自由 面把渗流场分成两个域，在自由面以上的子域中，各点流速为零；在自由面以下 的子域中，各点有大于零的流速。为了反映这一现象，可以在自由面以下子域中 采用实际的渗透系数，而在自由面以上的子域中令渗透系数为零。当然，为了计 算的稳定性，应取大于零的小数。b a t h e 建议； 当矿z 时：七：七 当p z 时；_ j ：k 1 0 0 0 这样，在迭代过程中，只要修改渗透系数，不必修改计算网格，因而可用固 定网格。实际计算步骤是；按全域剖分计算网格，首先假定全域均为实际渗透系 数，建立整体渗透矩阵解出水头；然后修改各单元的渗透系数，建立新的渗透矩 阵，解出新的水头。如此迭代计算，直至前后两次求出的水头充分接近。该算法 1 8 第二章渗流分析的基础理论和有限单元法原理 实际上是把边界不确定的非线形问题转化为材料非线形问题来考虑。 本法的优点是可用固定网格计算，但在每次迭代仍然要建立渗透矩阵并求逆， 计算量较大。而且该算法对含自由面单元组成区域的处理往往使迭代过程的收敛 性和解的稳定性不太好，结果会出现震荡现象。 2 2 4 2 虚单元法 虚单元法以上一次有限元计算求得的结点水头为基