渗流有限元分析理论2

1、渗流有限元分析渗流有限元分析理论理论 （1）设想将连续系统分割为数量有限的单元，单元之间通过指定）设想将连续系统分割为数量有限的单元，单元之间通过指定点即节点进行连接，用这样形成的单元集合体来替代原来的连续系点即节点进行连接，用这样形成的单元集合体来替代原来的连续系统。作用于系统上的外荷载用节点上的等效荷载来代替。统。作用于系统上的外荷载用节点上的等效荷载来代替。 （2）按一定的规则，对每个单元建立起求解未知量与节点相互作）按一定的规则，对每个单元建立起求解未知量与节点相互作用力之间的关系。用力之间的关系。 （3）按一定的条件集合全部单元，边界条件引进后就形成了一组）按一定的条件集合全部单元，

2、边界条件引进后就形成了一组以节点变量为未知量的代数方程组，进行求解就可得到某个节点处以节点变量为未知量的代数方程组，进行求解就可得到某个节点处的待求变量。因此，有限元法的实质是将拥有无限个自由度的连续的待求变量。因此，有限元法的实质是将拥有无限个自由度的连续系统，抽象为仅有有限个自由度的单元集合体，使问题更适于数值系统，抽象为仅有有限个自由度的单元集合体，使问题更适于数值求解。求解。2.1.2 2.1.2 有限元法的有限元法的基础理论基础理论 由上述可知，有限单元法的重要思想是将一个连续域离散化为有限个单由上述可知，有限单元法的重要思想是将一个连续域离散化为有限个单元，这些单元通过有限的节点相

3、互连接形成一个等效的集合体。这些单元之元，这些单元通过有限的节点相互连接形成一个等效的集合体。这些单元之间的连接方式、组合方法均可以不相同，单元自身的形状又是多种多样，因间的连接方式、组合方法均可以不相同，单元自身的形状又是多种多样，因此可以对几何形状比较复杂的模型进行数值模拟求解。该法是通过将全部求此可以对几何形状比较复杂的模型进行数值模拟求解。该法是通过将全部求解域内将要求解的未知区域函数分离到所有单元内假定的逼近相似函数分区解域内将要求解的未知区域函数分离到所有单元内假定的逼近相似函数分区域表示来实现的。与此同时，全部单元里的逼近相似函数由未知区域里的函域表示来实现的。与此同时，全部单元

4、里的逼近相似函数由未知区域里的函数在单元的各个结点的数值以及插值函数表示。随着单元个数的增加单元数在单元的各个结点的数值以及插值函数表示。随着单元个数的增加单元尺寸不断缩小、其自由度不断增加，插值函数更加精确。所得到的解与真实尺寸不断缩小、其自由度不断增加，插值函数更加精确。所得到的解与真实值的差距将会不断缩小，假如该单元的收敛条件得以满足，那么所求解最终值的差距将会不断缩小，假如该单元的收敛条件得以满足，那么所求解最终将趋于我们所希望得到的精确解。将趋于我们所希望得到的精确解。 有限元法有限元法对于单元应力应变的求解只要单元位移确定，就可以利用几何方程和对于单元应力应变的求解只要单元位移确定

5、，就可以利用几何方程和物理方程就可以求单元的应力和应变。下面仍以平面四节点矩形单元为例推导单元刚物理方程就可以求单元的应力和应变。下面仍以平面四节点矩形单元为例推导单元刚度矩阵。根据弹塑性力学中平面问题几何方程能得到单元里任何一点的应变表示式如度矩阵。根据弹塑性力学中平面问题几何方程能得到单元里任何一点的应变表示式如下：下：式式中中B称为应变矩阵，它的分块子矩阵可表示成：称为应变矩阵，它的分块子矩阵可表示成： xeeyijlmzBBBBB （2.1）000000iiiiiiiNxxNNByyNNNxyxy（2.2）代入插值函数后，上式可写为代入插值函数后，上式可写为根据弹性力学平面问题物理方程

6、，单元里每一点的应力能表示为：根据弹性力学平面问题物理方程，单元里每一点的应力能表示为：0000(1)010(1)4(1)(1)iiiiibBaabab（2.3）00,ii式中式中 xeeyxyDDBS （2.4）式中式中s表示应力矩阵，表示应力矩阵，D为弹性矩阼其表达式为为弹性矩阼其表达式为E是杨氏模量，是杨氏模量，是泊松比。是泊松比。应力矩阵应力矩阵S的分块子矩阵可用下式表示：的分块子矩阵可用下式表示：2101011012ED （2.5）211122iiiiiiiNNNNESNNxyxyyx（2.6）将其代入无量纲形函数后，上式变为将其代入无量纲形函数后，上式变为对于平面应变问题，只需将上

7、式中的对于平面应变问题，只需将上式中的E和和分别用分别用代替即可代替即可 0000200111141111122iiiiiiibaESbaabaa（2.7）2,11E22 渗流分析理论渗流分析理论 达西达西（Darcy） 于于1852年通过实验研究，总结出了渗流年通过实验研究，总结出了渗流水头水头损失损失与渗流速度之间的关系，后人称之为达西定律。与渗流速度之间的关系，后人称之为达西定律。 设设一均匀渗流装置，水经过一均匀渗流装置，水经过达西通过实验得出，圆筒内达西通过实验得出，圆筒内的渗流量的渗流量 Q Q 与渗流模型过水断面与渗流模型过水断面面积面积 A A 及水力坡度及水力坡度 J J 成

8、正比，与成正比，与土壤的透水性能有关，即土壤的透水性能有关，即lh1h2Q，A长为长为 L L的的砂滤层后，测压管水头砂滤层后，测压管水头由由 h h1 1 减小为减小为 h h2 2 ，于是其间的水于是其间的水力坡度力坡度 J J 可表示为可表示为12lhhhJLL式中：式中：V断面断面A上的平均流速，或称达西流速；上的平均流速，或称达西流速； J渗透坡降，即沿流程渗透坡降，即沿流程S的水头损失率；的水头损失率； k渗透系数；渗透系数； h测压管水头，它是压力水头与位置高度之和。测压管水头，它是压力水头与位置高度之和。12hhQAkLQdhvkkJAdS (2.8)(2.9) 它它指出了渗透

9、速度指出了渗透速度v与水力梯度与水力梯度J，或渗透坡降的线性关系。故又称为线，或渗透坡降的线性关系。故又称为线性渗透定律性渗透定律。从。从公式中可以看到达西渗透定律是把流速公式中可以看到达西渗透定律是把流速v与渗透坡降与渗透坡降J的关系的关系作为正比关系来考虑的，通过许多学者的研究证明这一正比关系在一定的条作为正比关系来考虑的，通过许多学者的研究证明这一正比关系在一定的条件下才能成立，太沙基通过大量试验证明从砂土到黏土达西渗透定律在很大件下才能成立，太沙基通过大量试验证明从砂土到黏土达西渗透定律在很大的范围内都能适用，其适用范围是由雷诺数（的范围内都能适用，其适用范围是由雷诺数（Re）来决定的

10、，也就是说达西）来决定的，也就是说达西定律仅适用于线性阻力关系的层流动力。定律仅适用于线性阻力关系的层流动力。 在在水利水利水电水电建设建设工程工程中中，除了一些水工建筑物诸如堆石坝、堆石排水体，除了一些水工建筑物诸如堆石坝、堆石排水体这样的孔隙相当大的介质中的渗流被认为是紊流之外。相当一部分渗流都包这样的孔隙相当大的介质中的渗流被认为是紊流之外。相当一部分渗流都包括于层流这个范围之内，均可运用达西定律进行渗流分析。根据达西定律所括于层流这个范围之内，均可运用达西定律进行渗流分析。根据达西定律所建立起来的渗流分析的数值分析模型，与很多工程实际问题还存在一些差别。建立起来的渗流分析的数值分析模型

11、，与很多工程实际问题还存在一些差别。故而应用达西渗流理论分析解决实际渗流问题时，应多考虑工程实际情况。故而应用达西渗流理论分析解决实际渗流问题时，应多考虑工程实际情况。加以灵活运用，加以灵活运用，2.2.2 渗流连续性方程渗流连续性方程 渗流渗流理论基础之一是水力学中的连续性方程，该方程是自然界中质理论基础之一是水力学中的连续性方程，该方程是自然界中质量守恒定律应用于渗流问题的一个具体表现，它表明，运行于渗透介质量守恒定律应用于渗流问题的一个具体表现，它表明，运行于渗透介质中流体，在其运行变化的整个过程中，其质量始终保持恒定不变，即不中流体，在其运行变化的整个过程中，其质量始终保持恒定不变，即

12、不会自动增加也不会自动减少。会自动增加也不会自动减少。 我们我们假设渗透进介质中流动的水是不能被压缩的均质液体，而且仅假设渗透进介质中流动的水是不能被压缩的均质液体，而且仅是涉及竖直方向这一单个方向的压缩，那么由是涉及竖直方向这一单个方向的压缩，那么由质质量量守恒定律守恒定律，通过对公，通过对公式加以变换就能得到如下可压缩介质的渗流方程（即连续性方程）：式加以变换就能得到如下可压缩介质的渗流方程（即连续性方程）：yxzvvvhgnxyzt（2.10） 式式中，中，为多孔介质压缩系数；为多孔介质压缩系数； 为水的压缩系数；为水的压缩系数； 为渗透水的密度；为渗透水的密度； g(+n) 为单位贮水

13、量或贮存率为单位贮水量或贮存率; vx，vy，vz 分别为渗流沿坐标轴方向的分速度。分别为渗流沿坐标轴方向的分速度。 假设假设水体和土体均为不可压缩的，则上述公式可转化为；水体和土体均为不可压缩的，则上述公式可转化为； 该该公式是不可压缩液体在刚体介质中流动的渗流方程，设方程式是连公式是不可压缩液体在刚体介质中流动的渗流方程，设方程式是连续性的，说明在流体中任何一点的单位流量（即流速）的纯变化率是零。续性的，说明在流体中任何一点的单位流量（即流速）的纯变化率是零。0yxzvvvxyz(2.11)2.2. 3渗流微分方程渗流微分方程将达西渗透定律和连续条件结合起来将达西渗透定律和连续条件结合起来

14、 式式中：中：vx，vy，vz分别为分别为x，y，z 三个渗透主轴方向上的渗流流速；三个渗透主轴方向上的渗流流速； kx，ky，kz分别为分别为x，y，z方向上的渗透系数方向上的渗透系数。代入。代入式（式（2.10）中，得）中，得式中，式中，h为总水头：为总水头：kx，ky，kz分别为分别为x，y，z方向上的渗透系数；方向上的渗透系数；x=g(+n) 为单位贮水量或贮存率；为单位贮水量或贮存率；t为时间。为时间。,xxyyxzHHHvkvkvkxyz (2.12)xyzxhhhhhkkkgnxxyyzztt(2.13)当土和水为不可压缩时，当土和水为不可压缩时，x=0，上式，上式(213)变为

15、变为即稳定渗流情况下的基本微分方程式：即稳定渗流情况下的基本微分方程式：在各个方向上渗透系数为常数时，上式为在各个方向上渗透系数为常数时，上式为0 xyzhhhkkkxxyyzz(2.14)2222220 xyzhhhkkkxyz(2.15)如果每一方向有相同性质，即如果每一方向有相同性质，即kx=ky=kz时，可转换成为常见的拉普拉斯方程时，可转换成为常见的拉普拉斯方程：由以上推理可直接得出二维稳定渗流的微分方程式由以上推理可直接得出二维稳定渗流的微分方程式式中：式中：h为总水头；为总水头； kx，ky分别为分别为x和和y方向的渗透系数。方向的渗透系数。2222220hhhxyz(2.16)

16、0 xyhhkkxxyy(2.17)2.2.4定解条件定解条件 流体流体的运动总是发生在一定的流场内，要确定流场的分布的运动总是发生在一定的流场内，要确定流场的分布仅渗流基本微分方程还是不够的，同时还必须依靠边界条件与仅渗流基本微分方程还是不够的，同时还必须依靠边界条件与初始条件来确定。边界条件是指顺着流场边界发挥主导作用的初始条件来确定。边界条件是指顺着流场边界发挥主导作用的条件；而初始条件指的是研究分析开始时流场内部整体流场状条件；而初始条件指的是研究分析开始时流场内部整体流场状态或对流动起决定作用的条件（譬如水头所处区域等）。通常态或对流动起决定作用的条件（譬如水头所处区域等）。通常所说

17、的定解条件就是指边界条件及初始条件这两个的合称。当所说的定解条件就是指边界条件及初始条件这两个的合称。当研究稳定渗流场问题时，无需考虑初值问题；而当研究非稳定研究稳定渗流场问题时，无需考虑初值问题；而当研究非稳定渗流场问题时，需要考虑边界条件，此条件有可能是变化的。渗流场问题时，需要考虑边界条件，此条件有可能是变化的。 因为因为稳定渗流场问题不需考虑初值条件稳定渗流场问题不需考虑初值条件，要求解要求解，就必须采用，就必须采用基本微分方程的定解条件之一即边界条件才可解决，我们称这样的基本微分方程的定解条件之一即边界条件才可解决，我们称这样的定解问题为边界问题。假设知道了所研究区域边界上的水头值，

18、那定解问题为边界问题。假设知道了所研究区域边界上的水头值，那么此边界条件即可以下式表示：么此边界条件即可以下式表示：式中，式中，1为渗流区域边界：为渗流区域边界： f(x，y，z)为已知函数；为已知函数；x、y、z处于边界处于边界1上。称此种边界条件为第一类边界条件，上。称此种边界条件为第一类边界条件，即水头边界条件。即水头边界条件。1, , ,h x y zf x y z(2.18) 假设假设不能确定研究的渗流边界上的水头，然而却知道边界上单位面积不能确定研究的渗流边界上的水头，然而却知道边界上单位面积流出或者流入的流量，那么这样的边界条件问题可用下式表示为：流出或者流入的流量，那么这样的边

19、界条件问题可用下式表示为： 式式中中， 2为已知流出或流入流量的边界段；为已知流出或流入流量的边界段； n为为2的外法线方向。的外法线方向。 当当边界不透水时，边界不透水时，q=0，上式可转变为，上式可转变为 称称此种边界条件为第二类边界条件，即流量边界条件。此种边界条件为第二类边界条件，即流量边界条件。2, ,hkq x y zn(2.19)0hn 当当研究非稳定的渗流问题时，则要将边界条件及初始条件全研究非稳定的渗流问题时，则要将边界条件及初始条件全部加以考虑。部加以考虑。水头边界：水头边界：流量边界：流量边界： 不透水边界我们认为是流量边界条件的特例，即不透水边界我们认为是流量边界条件的

20、特例，即1, , , , ,h x y z tf x y z t0t (2.20)2, , ,hkq x y z tn0t (2.21)0hn此外，还有以下定解条件此外，还有以下定解条件：混合边界：混合边界：初始条件：初始条件：其中其中h0是初始水头。是初始水头。hht00, , , ,th x y z thx y z(2.24)2.2.5 渗流分析的有限元法原理渗流分析的有限元法原理 渗流渗流分析的有限元法是通过利用变分原理将渗流基本微分析的有限元法是通过利用变分原理将渗流基本微分方程和它的边界条件变换成为泛函数的极值问题加以解决分方程和它的边界条件变换成为泛函数的极值问题加以解决的，它是一

21、种分块近似里兹的，它是一种分块近似里兹(Ritz)法的应用。首先将所研究区法的应用。首先将所研究区域以适当的方法分割成有限个适合求解的小区域，即单元，域以适当的方法分割成有限个适合求解的小区域，即单元，单元的角点即为结点，单元里建立单元局部方程的方法是借单元的角点即为结点，单元里建立单元局部方程的方法是借助于连续分片的插值函数，之后通过结点之问的关系装配全助于连续分片的插值函数，之后通过结点之问的关系装配全部单元成为一个整体，利用计算机，再通过线性代数方程组部单元成为一个整体，利用计算机，再通过线性代数方程组的形式进行求解。有限元法求解渗流场的水头函数的形式进行求解。有限元法求解渗流场的水头函

22、数h有其自有其自身的的方程身的的方程。其一般形式可表示如下其一般形式可表示如下： 式式中，中，K 为渗透矩阵；为渗透矩阵； h一为列向量；一为列向量； f为自由项列向量。为自由项列向量。 故此故此求解原有的渗流偏微分方程就等价于求解现在的代数求解原有的渗流偏微分方程就等价于求解现在的代数方程组。影响渗流有限元法求解渗流场精度的因素包括划分好方程组。影响渗流有限元法求解渗流场精度的因素包括划分好的计算区域对原来区域模拟的精度，分片插值对实际工程渗流的计算区域对原来区域模拟的精度，分片插值对实际工程渗流场模拟的精度和代数方程组自身求解的精度。场模拟的精度和代数方程组自身求解的精度。 Khf(2.2

23、5)2.2.6 有限元法计算渗流的基本公式有限元法计算渗流的基本公式 由由以上对渗流微分方程的分析，依据变分原理，将其和下面取以上对渗流微分方程的分析，依据变分原理，将其和下面取极小值的泛函进行等价。极小值的泛函进行等价。 式式中右端的最后一项为第二类边界积分，经过取泛函极小值后，中右端的最后一项为第二类边界积分，经过取泛函极小值后，在计算中便自动转化为第二类边界条件；而第一类边界条件，则是在计算中便自动转化为第二类边界条件；而第一类边界条件，则是在计算时直接赋给已知的边界水头值在计算时直接赋给已知的边界水头值。将渗流场划分成有限个单元。将渗流场划分成有限个单元后，便分解为各个单元的和，则泛函

24、式（后，便分解为各个单元的和，则泛函式（2.26）可以相对的分解为）可以相对的分解为某些单元泛函的和，如下式：某些单元泛函的和，如下式： 222212xyzshhhhI hkkkS hdxdydzqhdxyzt（2.26） 22221112mkxyzsejehhhhI hkkkS hdxdydzqhdxyzt（2.27） 以下以下顺次对式（顺次对式（2.28）中的各项取其导数，再求出这些项的最小值。）中的各项取其导数，再求出这些项的最小值。第一项第一项为为上上式对每个单元中各个节点的水头式对每个单元中各个节点的水头h1,h2,hm求导数，得：求导数，得： 222212312eeeexyzseh

25、hhhIhkkkS hdxdydzqhdIIIxyzt（2.28）222112exyzehhhIkkkdxdydzxyz（2.29）将单元将单元e内任一点表达式内任一点表达式代入代入式（式（2.30）得：）得：22212221212exyziiexyziiieIhhhkkkdxdydzhhxyzhhhkkkdxdydzhxhyhz（2.30）1MiiihN h11111122221,2,.,eMMMkikikixkykzkkkkieMkikikikxyzkeNNNNNNIkhkhkhdxdydzhxxyyzzNNNNNNhkkkdxdydz iMxxyyzz （2.31）令令则则：,jjjii

26、iijxyzeNNNNNNKkkkdxdydzxxyyzz 111112111212 2222121eMeeeMMMM MMeMIhKKKhIKKKhKhhKKKhIh (2.32)2esehIS hdxdydzt(2.33)求出单元求出单元e的的m个节点水头的导数，得：个节点水头的导数，得：令令，则，则:21111eMMkskkkkkieMkskikeMkskikehISN hNdxdydzhthSNN dxdydzthSN N dxdydzt（2.34）ijsijeSSN N dxdydz，其中，其中S是给水度，则：是给水度，则： 211112112212 2222122eMeeeMMMM

27、 MMeMIhSSShISSShShhSSShIh (2.35)hqSt23323111MMMek kkk kkkkhhIqhdShdSN hNdqN hdtt（2.36）232333111,eMkiikkiiiiMMIhqN dSNNdhthtqN dSN N dSN N dht（2.37）得：得： 313eeeeeMIhhPFtIh 2321;eeijijijMqN dPPPSN N dFqNd（2.38）（2.39）式中式中这样这样，对任意单元，对任意单元e，有：，有： 对所有单元的泛函求微分后，使其等于零，就得到泛对所有单元的泛函求微分后，使其等于零，就得到泛函对节点水头进行微分的函对节点水头进行微分的方程组即方程组即：式式中，中，n为节点总数。为节点总数。对对上面汇总的方程组写成矩阵形式，为：上面汇总的方程组写成矩阵形式，为： eeeeeeeeIhhKhSPFhtt（2.40） 0,1,2,.,neeiiIhI