利用初流量法分析有自由面渗流问题之改进.pdf

第 34 卷 第 2 期 岩 土 工 程 学 报 Vol.34 No.2 2012 年 .2 月 Chinese Journal of Geotechnical Engineering Feb. 2012 利用初流量法分析有自由面渗流问题之改进 潘树来 1，王全凤1，俞 缙1,2 1. 华侨大学土木工程学院，福建 厦门 361021；2. 中国矿业大学深部岩土力学与地下工程国家重点实验室，江苏 徐州 221008 摘 要在利用初流量法分析带有自由面的渗流问题时，常因采用高斯积分点作为结点初流量计算判断依据所带来的 误差而令解出现振荡。为了减低这一因素的影响及提高计算的效率，针对常用的 4 结点平面单元和 8 结点 6 平面三维 单元的自由面方程进行分析，提出采用坐标变换和等参变换等技术来改造结点初流量的计算，并建议按自由面穿越单 元之状况引用分部积分，使其积分上、下限符合高斯积分法的格式化要求，从而可利用高斯积分法求出精确的结点初 流量，而精确的结点初流量则有助于改善解的振荡及提高计算结果的精度。通过一矩形均质坝实例分析，表明该方法 的稳定性和收敛性良好。 关键词初流量法；有自由面渗流；高斯积分法；固定网格 中图分类号TV139.11 文献标识码A 文章编号1000–4548201202–0202–08 作者简介 潘树来1967– ， 男， 澳门人， 博士研究生， 从事计算方法和岩土工程研究。 E-mail punsuloi.hk。 Improvement of analysis of free surface seepage problem by using initial flow PUN Su-loi1, WANG Quan-feng1, YU Jin1,2 1. College of Civil Engineering, Huaqiao University, Xiamen 361021, China; 2. State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221008, China Abstract When the initial flow is applied in solving problems of free surface seepage, we often find an oscillating solution determined due to the error caused by taking the Gaussian integration points as the assumed nodes based on calculation of nodal initial flow. In order to reduce the impact of this factor and to improve the efficiency of computation, an analysis is made with the free surface equation commonly used by 4-node plane elements and 8-node 6-plane three-dimensional elements, a technique of transation of coordinates as well as isoparameters is proposed to improve the calculation of the nodal initial flow, and an introduction of partial integration is also suggested according to the state of free surface penetrating elements, making the integration upper limit and lower limit in consistence with the atting of the Gaussian requirements. Accordingly, an exact nodal initial flow integration will be obtained by using the Gaussian , and this accurate nodal initial flow will then benefit the itself improvement of the oscillating solution and increase the precision of the calculation. Through the analysis of a homogeneous rectangular case, it is shown that a fair state of stability and convergence has been achieved by using the proposed improved . Key words initial flow ; unconfined seepage; Gaussian integration ; fixed mesh 0 引 言 利用有限元技术分析基坑、边坡、堤坝等岩土工 程时，常会遇到地下水渗流的问题。处理这类有自由 面渗流问题的难点在于其自由面（浸润面）位置为未 知，一般需要通过迭代求解。目前求解的方法可分为 变网格法和固定网格法两大类。 变网格法是确定自由面的一种传统方法，虽然已 取得了不少成功经验，但此方法也存有致命缺点每 次迭代都要先确定自由面的位置及调整总体传导矩 阵，计算量大；同时，在调整网格时容易造成单元网 格畸变或交错重叠，影响计算精度；另外，由于非渗 流区域没有被考虑到分析系统中去，故不适用于作应 力场与渗流场耦合分析，局限了它在工程中的应用。 因此，变网格法已逐步被淘汰。 为了解决变网格法的缺点，学者们开始尝试采用 ─────── 基金项目国家自然科学基金项目（51109084） ；高等学校博士学科点 专项科研基金项目（200803850001） ；福建省自然科学基金项目 （2011J01317） ； 中国矿业大学深部岩土力学与地下工程国家重点实验 室开放基金项目 （SKLGDUEK1012） ； 泉州市科技计划项目 （2010Z81） 收稿日期2010–11–09 第 2 期 潘树来，等. 利用初流量法分析有自由面渗流问题之改进 203 扩大分析系统范围的办法，将渗流区与非渗流区结合 起来形成总分析系统，并引入总系统的边界条件作分 析，以达到求解过程中无需变更单元网格的目的，从 而形成了固定网格法。自 Neumann 于 1973 年提出用 不变网格分析有自由面渗流的 Galerkin 法以来，先后 出现了多种固定网格法，如剩余流量法[1]、单元渗透 矩阵调整法[2]、初流量法[3]、节点虚流量法[4]、虚节点 法[5]、改进节点虚流量法[6]等。在众多固定网格法中， 由于初流量法具有理论明确、只需一次形成总体传导 矩阵、无需在每次迭代时确定自由面的近似位置和判 别自由面与单元相交状况等优点，因此，初流量法是 目前利用固定网格法求解有自由面渗流问题的一种较 为有效方法。 在应用初流量法时，关键是要算出自由面每次迭 代时由初流量所引起的结点初流量值 Q0（记作结点 “初” 流量， 用作区分由边界条件所产生的结点流量） ， 而 Q0的误差与选用的积分计算方法密切相关。 对于那 些被自由面穿过的单元，初流量法建议只取位于自由 面之上的积分点作计算，明显地，此方法未能充分考 虑由自由面所分割出非渗流区的实际形状及此位置上 被积函数的大小和分布， 因此计算得出的 Q0带有很大 的主观性和误差，对于那些体积较大或自由面刚好位 于积分点附近的单元来说，误差就更大。而过大的误 差将影响渗流场有限元方程的迭代计算，令解出现振 荡。 为了减小结点初流量之误差，人们提出了多种不 同的改进方法，例如朱伯芳[7]建议自由面附近单元 每排高斯积分点的数量应不小于 3，使每个高斯积分 点控制的面积尽可能小； 王媛[8]提出的改进初流量法， 针对影响初流量法解不稳定的主要因素提出了区域识 别函数的新概念，对不连续的区域识别函数进行了线 性微调，使之连续化，从而降低了初流量法解的不稳 定性；李新强[9]提出引入连续调整系数法的概念来改 进利用初流量法求解自由面的计算精度，实例表明， 经改进后的初流量法稳定性和收敛性良好。 总的来说，上述改进初流量法均建立在减小结点 初流量的误差的基础上，但这些方法均未能彻底解决 计算方法本身所造成的误差。而本文提出的改进，主 要通过对自由面方程进行研究，利用换元技术将自由 面由曲线变为直线，然后再透过等参变换及高斯积分 法直接求出结点初流量。基于高斯积分法理论，若选 取的高斯积分点数量足够，则计算所得的数值解为精 确解， 因此， 本文改进方法只会出现计算精度之误差， 并不会出现计算方法上之误差。 1 初流量法简介 初流量法的基本原理类似于非线性应力分析中的 初应力法，即在达西定律中增加一个初流量项，通过 对初流量的调整， 将非线性分析化为一系列线性分析， 相关方程的推导可参见文献[3，8]，其结果为    0 MhQQ ， 1 式中，[M]为总体传导矩阵， h为结点水头列矩阵，  Q为由边界条件引起的结点流量列矩阵， 0 Q为由 初流量引起的结点初流量列矩阵，       e T T 00e e e dQABKKB h    ，2 式中， eA为整体组装的单元选择矩阵，[B]为几何矩 阵，其关系如式（3）所示， e h为单元内结点水头，[K] 为渗透矩阵，[K0]为渗透修正矩阵， ij K， 0ij K为介质 渗透张量，饱和区 0ijij KK，非饱和区 0 0 ij K，   128 ,,,BB BB, T ,, iii i NNN B xyz      ， 3 其中，Ni为8结点单元的形函数。 为 便 于 表 述 ， 记    T 0 , ,F x y zBKK   e B h，则结点初流量列矩阵方程可表示为   e T 0 e e , ,dQAF x y z    。 4 从式（2）可以看出，饱和区单元的结点初流量为 0，非饱和区单元的则可通过数值积分计算得出。 2 初流量法之改进 2.1 二维4结点四边形渗流单元 对于图1（a）所示的平面4结点四边形单元，经 等参变换转为图1（b）时，单元的水头h和高度坐标 z、结点初流量列矩阵可分别表示为 1 1223344 1 1223344 hN hN hN hN h zN zN zN zN z     ， ， 5  T 0e e e ,d dQAF x zJ      。 6 式（6）中，由于F（x，z）要用到形函数对旧坐 标（x，z）的微分，但形函数是以新坐标（ξ，ζ）表 示，两者需通过雅可比矩阵作变换。另外，由自由面 定义（hz）可得出自由面曲线方程为  111222333444 0NhzNhzNhzNhz 。 7 式（7）是一个带有常数，ζ和ξζ项的二次多项 式，基于组成特点可改写为  1234 0aaaa  ， 8 式中， 204 岩 土 工 程 学 报 2012 年     111223344 211223344 311223344 411223344 ahzhzhzhz ahzhzhzhz ahzhzhzhz ahzhzhzhz               ， ， ， 。 9 式（8）经整理后，可得  23142 4434 a aa aa aaaa       。 10 令  2314 434 a aa a p aaa    进行坐标变换， 则可将自由面 方程由曲线变为式（11）表示的直线，它与单元的关 系如图1（c）所示，而位于自由面上结点3，4的坐 标值可通过式（11） 、 （12）求出 2 4 a p a   ， 11  2314 434 23143 2 44 a aa a p aaa a aa aa a pa              ， 。 12 图 1 带自由面的平面 4 结点单元 Fig. 1 A 4-node plane element with free surface 将式（11） 、 （12）代入式（6） ，便可得出经坐标 变换后的结点初流量方程为   s4 2 3 4 1T 0 e e , s xx pp a ppp a QAF x z       1423 22 4 d d a aa a Jp a p    。 13 式（13）中，由于积分区间是一个如图1（c）所 示的四边形，故可再进行一次等参变换，使之变成如 图1（d）的标准单元，然后便可通过高斯积分法算出 结点初流量的值。在等参变换时，新旧坐标关系如式 （14）所示，另外，由于旧单元有一组对边与坐标轴 平行，其雅可比矩阵必有元素vp 或u为0，因 此，新旧坐标面积元的关系为     12 34 12 34 1 1111 4 1111 1 1111 4 1111 puv puv p uv puv p uvuv uvuv                 ， ， 14   12341234 12341234 1 d d 16 d d pppppppppv uu v         。 15 将式（15）代入式（13） ，经整理后便可得出最终 的结点初流量方程为   11T 0p e 11 e ,d dQAF x zJu v        ， 16 式中， p 可被称为二维转换因子，  1423 p1234 22 4 16 a aa a pppp a p      1234 ppppv   12341234 u  。 17 式（16）的积分上、下限已符合高斯积分法的格 式化要求，可通过式（18）高斯积分法求出其数值解。   T 0p e e11 , nn ijij ij QAH HF x zJ      ，18 式中的变量x和z，应先透过雅可比矩阵变为ξ和ζ， 然后再通过式（12） 、 （14）转换为最终的变量u和v。 （1）单元结点初流量计算结果之验证 为了验证本文的改进方法，现以图1所示的均质 单元（KijK）为例，分别采用积分法、以积分点是否 位于非渗流区为判断准则的初流量法及本文改进方法 进行结点初流量计算，结点初流量的积分解为 01 Q 1.152704078K， 02 Q0.930373421K， 03 Q -0.650222894K， 04 Q-1.432854606K，其它方法的计 算结果列于表1，2。 从表1可看出，初流量法算出的结点初流量误差 很大，达到以倍数计，而随着积分点数量增加，误差 并没有改善。从表2可看出，本文改进方法能有效地 降低误差，随着积分点数量增加，误差改善明显，当 每行积分点取为5时，计算结果的相对误差只有百万 分之几，说明本文方法计算效果理想。 （2）自由面穿越平面单元的可能性及其模拟 从式 （8） 可看出， 自由面是一条光滑的二次曲线， 它穿越单元的种类可归纳为图2所示的4种情况。对 于图2（a） 、 （b）的情况，可按前述方法直接求出结 点1234范围的结点初流量。至于图2（c） 、 （d）之情 第 2 期 潘树来，等. 利用初流量法分析有自由面渗流问题之改进 205 表 1 由初流量法算出之结点初流量数值解及其与积分解之比较 Table 1 Numerical solutions calculated by initial flow and comparison between them and integral ones 01 Q 02 Q 03 Q 04 Q 每行积 分点数 数值解 ε 数值解 ε 数值解 ε 数值解 ε 4 0.483804197K -0.5803 1.742284226K0.8727 -0.046087282K-0.9291 -2.180001141K0.5214 5 0.658980423K -0.4283 1.942164814K1.0875 -0.114630806K-0.8237 -2.486514430K0.7354 6 0.535642414K -0.5353 1.747497602K0.8783 -0.080516974K-0.8762 -2.202623042K0.5372 注表中 ε数值解-积分解/积分解。 表 2 由本文改进方法算出之结点初流量数值解及其与积分解之比较 Table 2 Numerical solutions calculated by present improved and comparison between them and integral ones 01 Q 02 Q 03 Q 04 Q 每行积 分点数 数值解 ε 数值解 ε 数值解 ε 数值解 ε 4 1.152805902K 8.83310 -5 0.930271439K-1.09610 -4 -0.650285165K 9.57710 -5 -1.432792176K-4.35710 -5 5 1.152708965K 4.24010 -6 0.930367909K-5.92510 -6 -0.650225256K 3.63310 -6 -1.432851618K-2.08510 -6 6 1.152704234K 1.35310 -7 0.930373205K-2.32210 -7 -0.650222911K 2.61410 -8 -1.432854528K-5.44410 -8 注表中 ε数值解-积分解/积分解。 况，由于自由面曲线光滑及其上、下方必为非渗流区 及渗流区，故自由面与ξ1边界必有交点，基于分 部积分原理， 可分别求出结点1234和1 2 3 4  范围的结 点初流量，然后相减得出所需的解。 图 2 自由面穿越单元之状况 Fig. 2 Conditions of free surface through element 2.2 三维8结点6平面渗流单元 对于图3（a）所示的等参单元，由自由面定义可 得出自由面的曲面方程为  111222888 0NhzNhzNhz。19 经整理后，式（19）可改写为 1234 5678 bbbb bbbb        ， 20 式中， 111223344 55667788 211223344 55667788 311223344 55 bhzhzhzhz hzhzhzhz bhzhzhzhz hzhzhzhz bhzhzhzhz hzh          ， ， 667788 411223344 55667788 511223344 55667788 6112233 zhzhz bhzhzhzhz hzhzhzhz bhzhzhzhz hzhzhzhz bhzhzhz          ， ， ， 44 55667788 711223344 55667788 811223344 55667788 hz hzhzhzhz bhzhzhzhz hzhzhzhz bhzhzhzhz hzhzhzhz            ， ， 。                         21 在式（20）中，由于其分母带有ξ和η的相加和 相乘项， 无法将公式右侧部分分拆为ξ和η的独立项， 故不能直接采用2.1节中坐标变换方法将自由面方程 由曲面变为平面。然而，基于三维高斯积分原理，可 先对某一坐标进行处理，将原属三维积分问题变为二 维问题，从而可引用前述方法作计算。 由于自由面方程具有光滑及其上、下方必为非渗 流区和渗流区的特点，当自由面同时穿过单元的顶、 底面时，切割出顶面的范围必大于底面，因此，在选 定先进行高斯积分计算的坐标时，应以自由面与单元 顶面相交的状况作参考。另外，为了确保符合高斯积 分的格式化要求，还应遵从以下两项原则①选定的 206 岩 土 工 程 学 报 2012 年 坐标，其积分区域应为[-1，1]或可通过单一坐标变 换变为[-1，1]；②在选定坐标的积分区域内，其它 两个坐标的积分下限应只有自由面一组方程，若不能 满足时，可参考分部积分原理将单元划分为不同部分 去处理。 对于图3（a） 所示的等参单元， 其顶面如图3（b） 所示，应取坐标先行处理。在此实例中，的积分 区域刚好为[-1，1]，否则，可以以方向下方的坐 标 x 及上方坐标 s 作变换  1 xsxs2 q   ， 22 而相应的自由面方程变为   1xs32x4 5xs76xs8 22 22 s bbbb bbbb            xs3xs4 xs7xs8 b qbq b qbq        。 23 图 3 带自由面的 8 结点 6 平面等参单元 Fig. 3 8-node with 6-plane isoparametric element with free surface 以式（22）代入结点初流量方程，并先对坐标q 进行高斯积分处理后，结点初流量方程及自由面方程 可改写为   T xs 0 e e1 ,,d d 2 i i n ii i QAHF x y zJ           ， 24 12 34 cc cc        。 25 式中， 11xs3xs322 22 i cbbb qcb ，  xs4xs435 2 i bb qcb ， xs7 b  xs746 2 i b qcb， xs8 b xs8i b q。 在对坐标q先进行高斯积分处理后，原三维积分 问题会变成在高斯积分点qi上的二维积分问题，而由 于自由面在qi平面上变成一条如式（8）所示的二次 多项式，因此，完全可参照2.1节所介绍的平面问题 方法作处理，最终可得出三维8结点6平面渗流单元 的结点初流量方程为   0s e e111 ,, jik nnn T ijkjik ijk QAH H HF xy zJ        ， 26 式中， s 是三维转换因子，  1 42 3 s12341234 22 4 32 c cc c ppppppppv c p       12341234sx u   。 27 与平面单元相比，自由面穿越三维空间单元的状 况要复杂得多，单以被自由面切割出位于非渗流区的 结点就有1～7个这7种情况， 而且每种情况也有不同 形式的组合。因此，需要找出这些情况的规律，才能 用简洁的方法进行模拟。 在使用高斯积分法时，必须确保其积分边界由唯 一的一道方程组成，否则，应采用分部积分处理，基 于此原则性要求，可按积分边界的实际情况来划分分 部积分的范围。 在参照2.2节的方法选定水平面坐标η 先进行高斯积分计算时，可以以自由面与单元顶、底 面相交的情况作判断。当自由面与单元底面相交而使 积分下限由自由面和单元底面两部分组成时，可先计 算由单元顶面与自由面所组成的部分，然后再扣除由 单元底面与下方自由面所组成的多算部分。自由面与 单元顶、底面相交大致有图4所示的4种情况，对于 （a） 、 （c）两种情况，可直接求出由结点1256和自由 面所组成范围的结点初流量，至于（b） 、 （d）两种情 况，可分别计算由结点1256和自由面、及1 2 5 6   和 自由面所组成范围的结点初流量，然后相加便可得出 答案。总的来说，不管自由面穿越单元的状况有多复 杂，按以上方法最多只需分四部分积分便可得出单元 的结点初流量。 图 4 自由面穿越单元顶面或底面之状况 Fig. 4 Conditions of free surface across top or bottom of element 3 算 例 为考察本文的改进方法，现以图5所示的尺寸为 第 2 期 潘树来，等. 利用初流量法分析有自由面渗流问题之改进 207 表 3 本文改进方法与电模拟试验及其它数值方法算出的自由面位置之比较 Table 3 Comparison among free surface locations calculated by present improved , electrical simulation tests and other numerical s X/m 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 电模拟试验解[4] 10.00 9.78 9.31 9.05 8.31 7.91 7.64 6.93 5.86 5.47 4.73 节点虚流量法[4] 10.00 9.74 9.38 9.01 8.51 8.03 7.41 6.77 6.05 5.19 4.33 初流量法[3] 10.000 9.745 9.3949.0108.520 8.0307.4296.807 6.085 5.179 4.373 改进初流量法[8] 10.000 9.745 9.3939.0078.517 8.0257.4286.805 6.076 5.189 4.375 本文方法 10.000 9.745 9.3909.0068.517 8.0267.4276.807 6.073 5.190 4.372 表 4 单元内取不同积分点数对自由面位置计算结果之影响 Table 4 Effect on free surface locations as calculated by different integral points in an element X/m 1.0 2.0 4.0 5.0 6.0 8.0 9.0 10.0 n4 9.74456422 9.38660153 8.526680538.019625047.423057216.04434523 5.19150219 4.35736189 n6 9.74708443 9.39383819 8.508312008.028842657.428986086.08077863 5.17896002 4.37434586 n8 9.74526472 9.39423747 8.519765298.030425847.428742486.08482345 5.17892745 4.37291666 n12 9.74505825 9.38907987 8.515099118.027985207.426504196.06731194 5.19102046 4.37092718 初流 量法[3] n20 9.74543610 9.38906513 8.516506798.024234047.424850146.06825654 5.18957853 4.37307322 n4 9.74411992 9.39383364 8.514325098.032628017.426838496.06099322 5.18786827 4.36537094 n6 9.74217857 9.39353273 8.513154618.030446277.431668396.08984415 5.19330771 4.37448693 n8 9.74465622 9.39300503 8.517282758.024845757.427631786.07570802 5.18872026 4.37489685 n12 9.74561049 9.38877210 8.514990168.026711127.426880406.06859015 5.19193588 4.37178340 改进 初流 量法[8] n20 9.74524344 9.38964076 8.516754728.025265647.425600066.06989673 5.18988792 4.37208832 n4 9.74483860 9.39000990 8.517007198.026050937.427179946.07300624 5.19015159 4.37167357 n6 9.74483847 9.39000957 8.517006318.026049897.427178226.07300392 5.19017199 4.37165322 n8 9.74483847 9.39000957 8.517006308.026049897.427178226.07300391 5.19017214 4.37165306 n12 9.74483847 9.39000957 8.517006308.026049897.427178226.07300391 5.19017214 4.37165306 本文 方法 n20 9.74483847 9.39000957 8.517006308.026049897.427178226.07300391 5.19017214 4.37165306 注表中 n 为单元每行的高斯积分点个数。 10 m10 m1 m均质土坝作分析，土坝上、下游水位 分别为10 m和2 m，底部为不透水边界，下游渗出面 边界取为hz。建立计算模型时网格划分为100个1 m1 m1 m的8结点6平面单元，共有242个结点。 在按自由面位置算出结点初流量后，再使用迭代法分 析时，计算精度取1.010 -10。利用自编的三维有限 元程序，在取单元每行的积分点为8个及自由面经5 次迭代后，得出的自由面位置与其它几种方法的对比 见于表3， 自由面及等势线关系见图5， 从表3中的计 算结果与试验解及其它数值解之吻合状况可说明本文 改进方法正确可靠。 对于图5的实例，当采用不同的积分点时，自由 面经5次迭代计算后得出的自由面位置见于表4。从 表4可看出，随着积分点数量增加，由初流量法及改 进初流量法算得的自由面位置并不稳定，而本文方法 在取单元每行的积分点为8个及更多后，计算出的结 果已基本相同，此情况只有在结点初流量值相等或很 接近时才会出现，说明本文方法在取一定积分点后计 算得出的结点初流量精确可靠。 图 5 本文改进方法分析出的自由面及等势线位置 Fig. 5 Free surface and equipotential line introduced from analysis of present improved 为了解本文方法的收敛性和稳定性，现仍以图5 之实例及取单元每行的积分点为8个进行分析，得出 相连两次迭代后自由面的差值如表5所示，经不同次 数迭代后自由面位置平均振幅见于表6。从表5可看 出，随着迭代次数增加，改进初流量法及本文方法的 208 岩 土 工 程 学 报 2012 年 表 5 迭代过程中自由面位置之差值 ΔZx Table 5 Relative variation of positions of free surface during process of iteration 10 -3 m X/m 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 m1 0.48960 -3.76474 7.28274 -2.962927.69665-1.273083.7108815.37869 -13.878606.99741 m2 -1.79971 -3.97404 -6.56062 -7.14040-1.95096-7.02647-9.36627-18.62471 -7.618128.64576 m3 -0.16075 -0.20919 -0.18919 -0.