地下工程渗流排水孔数值模拟的隐式复合单元法

第29卷第6期岩土力学Vol.29No.62008年6月RockandSoilMechanicsJun.2008收稿日期：2006-09-06基金项目：国家自然科学基金重大研究计划支助项目（No.90715042No.50725931作者简介：倪绍虎，男，1982年生，博士研究生，主要从事地下工程渗流数值模拟方面的研究。E-mail:[e-mail]文章编号：1000－7598－(2008)06－1659－06地下工程渗流排水孔数值模拟的隐式复合单元法倪绍虎，肖明（武汉大学水资源与水电工程科学国家重点实验室，武汉430072）摘要：根据渗流基本原理，提出了排水孔在地下工程中渗流数值模拟的隐式复合单元法。将排水孔隐含于模型网格中，得到排水孔在模型中的数据信息。判断所有排水孔与单元的相对位置，通过修正排水孔所穿过单元的渗透传导矩阵来模拟排水孔的强排水效果，并通过实例对排水子结构法和隐式复合单元法进行了比较。运用编制的三维有限元计算程序对某地下水电站厂房渗流场及排水孔进行了模拟，结果表明，采用隐式复合单元法对排水孔进行模拟是可行的。关键词：渗流；数值模拟；排水孔；隐式复合单元法中图分类号：O357文献标识码：AStudyonnumericalsimulationofdrainageholeofseepagefieldinundergroundengineeringwithimplicitcompositeelementmethodNIShao-hu,XIAOMing(StateKeyLaboratoryofWaterResourcesandHydropowerEngineeringScience,WuhanUniversity,Wuhan430072,China)Abstract:Anumericalsimulationmethodofdrainageholeinthree-dimensionalseepagefieldinundergroundengineeringispresentedaccordingtothetheoryofseepageflow.Drainageholesareimplicitinmodelandtheirdataaregainedaccordingtopre-process.Therelativepositionsofdrainageholesarejudgedandtheirimpactsaresimulatedthroughmodifyingthepenetrabilityconductionmatrixesofelementsperforatedbythedrainageholes.Andthendrainagesubstructuremethodandimplicitcompositeelementmethodarecomparedviaanexample.Finally,athree-dimensionalfiniteelementmethodcalculationprogramcompiledbyourselvesisusedtosimulatetheseepagefieldofanundergroundhydropowerhouse.Thecalculationresultshowsthatthenewmethodintroducedisfeasibletosimulatedrainageholes.Keywords:seepage;numericalsimulation;drainagehole;implicitcompositeelementmethod1引言在地下工程、大坝和边坡等结构工程设计中，水的渗流问题是必然面临和必须解决的问题。排水孔是这些工程中最主要的渗流控制措施之一，应用非常广泛，正确、有效地模拟排水孔的排水效果是设计中的重要工作。国内外对排水孔的模拟方法很多，如排水子结构法、等效杆单元法、以管代孔法、以缝代井法、空气单元法和改进截止负压法等[1－3]。排水子结构法不增加整体自由度，只增加子结构的凝聚计算，可以较精确地模拟排水孔，可以运用于排水孔较少的工程。但一般工程中排水孔数目众多，采用子结构法对单个排水孔一一模拟会使工程量加大，使计算复杂，也不经济，因此，很难推广运用。等效杆单元法计算简单，适用性好，但因为它是在平均意义上模拟的结果，所以计算精度不够高。有的计算方法因为人为因素较大，会形成较大误差，计算精度不高。还有一些计算方法主要是对排水孔的节点进行水头约束，这样会加大计算难度，增加了人为的非计算误差，对准确模拟排水孔不利。本文主要针对以上问题，提出了排水孔的隐式复合单元法。运用打“锚杆”的处理方法直接在模型中引入排水孔，此类排水孔隐含在模型网格中而没有形成单元、节点，因此，没有改变模型，只需修正复合单元的渗透传导矩阵。此法只需增加少量的数据准备工作和计算量，可以较好地模拟排水孔岩土力学2008年以满足实际工程的要求，并运用自行编制的三维有限元计算程序对某地下水电站厂房渗流场及排水孔进行了模拟。2基本原理岩土工程中岩土材料的渗透特性一般表现为各向异性和不均匀性。对于岩土的各向异性，可以采用渗透系数矩阵进行模拟。而对于岩土材料的不均匀性则可以采用不同材料单元的单元渗透系数来模拟[4]。根据复合单元法的基本思想，在同一单元内部也会存在不同渗透特征的渗透子域，它们的渗透系数各不相同，此时单元的渗透系数需用复合单元法来确定。排水孔实际是中空的，为了计算方便，可以将其作为强渗透介质。排水孔会穿过单元或包含在单元内，可以将排水孔视为单元内的子域，单元内包含岩体、排水孔、断层等，此时的单元变成了复合单元。在渗流场离散化建立模型网格时，首先不考虑排水孔剖分网格。模型网格建成后，按照打“锚杆”的方法加入排水孔，然后输出排水孔的信息数据，包括排水孔的起点坐标、终点坐标、起点单元号、终点单元号、直径等信息。在进行有限元计算时，将排水孔和它穿过的岩体单元合并考虑为复合单元，对该单元的渗透传导矩阵进行等效处理加以修正，这样就可以对排水孔进行模拟。该法只修正了复合单元的渗透传导矩阵，没有对排水孔进行实际单元模拟。排水孔只是隐含于岩体中，并没有改变网格，建模简单，加入排水孔和输出排水孔信息数据容易，方法简单可行，便于运用。3基本方法和步骤3.1渗流计算方法本文采用改进的丢单元法[5]进行渗流计算，根据前一次有限元计算所求得的节点势，用插值法求出自由面与单元的交点和相对位置。自由面将单元分成上下两个区域：自由面以上为渗流虚区，将此区单元（空单元）和单元的节点（虚节点）固定，不参与下一次迭代计算；自由面以下为渗流实区，该区域为渗流计算区域。为了避免含自由面单元处复合缝面材料发生突变，采用分段线性函数过渡对含自由面单元的渗透矩阵进行调整[6]。通过迭代逐步逼近渗流场真实解。采用丢单元法不需要改变模型网格，迭代速度快，计算精度相对较高，渗流计算中累积误差小。3.2复合单元的判断在修正复合单元的渗透传导矩阵前，首先要判断排水孔穿过的单元，即复合单元。根据读入的单元信息和排水孔信息在每个单元内对所有排水孔与它的相对位置情况分别进行判断，求解排水孔的中轴线与单元6个面的相交情况，若交点在排水孔的始点和终点之间，则排水孔与单元存在交点（视排水孔为线段）。排水孔与单元的相对位置有3种：第1种是排水孔与单元相离，没有交点；第2种是排水孔贯穿单元，有两个交点；第3种是排水孔一端位于单元内，有1个交点。记录所有单元的交点情况，有1个或2个交点的单元即为复合单元。3.3复合单元内排水孔主渗透系数方向的确定在判断复合单元的同时，分别判断所有排水孔与每个单元的相对位置。如果有交点，则求解出交点坐标，并记录每个单元内的排水孔个数，特别是当只有一个交点时，要准确判断排水孔的哪一端位于单元内。而排水孔的走向、倾角不同，会使单元内的渗透情况有差异，排水孔沿轴向的渗透性最强，沿径向的渗透性也有所加强。将排水孔视为一个圆柱体，确定每个排水孔轴向方向向量，取向下为正。将所有排水孔的方向向量通过矢量叠加后得出复合单元内所有排水孔等效后的轴向主渗透系数（简称主渗）方向向量。如图1为一个8节点单元内有两个排水孔A和B，排水孔A穿过复合单元的面2和3，排水孔B穿过面1和3。图2为图1中2个排水孔的轴向方向向量平移后在2个向量组成的平面上的投影，根据向量相加得到此复合单元内等效后排水孔轴向的方向。假定复合单元内等效后的排水孔轴线方向向下为它的x主渗方向。渗透系数为xK′。垂直于轴线方向的任意一个方向为它的y主渗方向。为方便起见，可将y主渗方向取在x主渗方向与y坐标轴组成的平面内。x主渗方向向量与y主渗方向向量的叉积即可得到z主渗方向向量，如图3所示。同理，可以求出复合单元内排水孔等效后的y和z主渗yK′、zK′的方向。图1复合单元内排水孔位置示意图Fig.1Thedrainageholesincompositeelement1660第6期倪绍虎等：地下工程渗流排水孔数值模拟的隐式复合单元法图2排水孔向量相加示意图Fig.2Thevectorsofdrainageholes图3等效排水孔主渗方向示意图Fig.3Themainseepageorientationsofdrainagehole3.4排水孔等效渗透系数矩阵p[]K的取值排水孔的等效主渗透系数根据排水孔的数目、间距、长度等按照文献[7]的等效法取值。事实上，排水孔沿轴向和垂直于轴向的渗透性是有差异的，表现出各向异性，沿轴向方向是中空的渗透系数大于垂直于轴线方向的渗透系数。排水孔的主渗透系数按文献[7]的建议取值，这样就可以得到排水孔的主渗透系数矩阵[K′]。根据所有排水孔等效后在复合单元内的位置可以计算出3个主渗透系数与坐标轴的方向余弦xl，yl，zl。由方向余弦组成坐标转换矩阵[Rp]：111p222333[]xyzxyzxyzlllRllllll⎛⎞⎜⎟=⎜⎟⎜⎟⎝⎠（1）根据下式计算出整体坐标下的排水孔等效渗透系数矩阵p[]K：Tppp[][][][]KRKR′=（2）3.5复合单元的渗透系数矩阵本文采用6面体8节点等参单元。复合单元如图4所示。根据复合单元内所有排水孔的信息可以求得排水孔在单元内的体积pV，单元体积为V，单元内岩体体积为rV。排水孔的渗透系数矩阵为p[]K，岩体的渗透系数矩阵为r[]K。则有：图4复合单元示意图Fig.4ThecompositeelementrpVVV=−（3）复合单元的渗透系数矩阵按有限元法单元渗透传导矩阵计算原理来计算。复合单元渗透传导矩阵可以表示为[][][][][][]TTerrppddrpABKBvBKBvΩΩ⎡⎤=+⎣⎦∫∫∫∫∫∫（4）式中：[B]为单元几何矩阵；r[]K，p[]K分别为岩体和排水孔的渗透系数矩阵。复合单元具有整体性，那么就需要将不同材料渗透系数等效成单元的等效渗透系数。假设复合单元的等效渗透系数矩阵为c[]K，根据有限元方法计算出复合单元的渗透传导矩阵，则有：Tec[][][][]dABKBvΩ=∫∫∫（5）联立式（4）、式（5）两式，即可求出复合单元等效渗透系数矩阵c[]K。4算例为了检验用隐式复合单元法模拟排水孔的正确性，笔者同时采用排水孔子结构法和复合单元法模拟了排水孔的排水效果，并对两种方法的计算结果进行了比较。有限元计算采用改进的丢单元法，计算区域为（0～30，0～24，0～30m），上游水位为27.00m，下游水位为15.00m。模型中布置有8个直径为100mm的排水孔，长度为13.5m，间距为3.0m。如图5、图6所示为排水子结构法的计算模型，图7、图8所示为隐式复合单元法的计算模型。其中为了控制计算误差，采用隐式复合单元法模拟时排水孔附近的单元应适当加密，使单元体积不会太大，并使网格从排水孔向四周渐变。经笔者反复模拟，若排水孔附近单元体积控制在10～20m3甚至更小，计算误差很小。排水孔的渗透系数按文献[7]的建议取值，x主渗透系数取为围岩的5000倍，y和z主渗透