《南京工业大学学报》--承压含水层地下水稳定流的有限层分析-诸宏博

收稿日期：基金项目：国家自然科学基金项目（50278042）；江苏省高校自然科学研究计划项目（03KJB560044）作者简介：诸宏博（1981-），男，浙江杭州人，硕士生，主要研究方向为环境岩土工程及数值分析；王旭东（联系人），教授，E-mail：[e-mail]承压含水层地下水稳定流的有限层分析诸宏博,王旭东,宰金珉(南京工业大学土木工程学院，江苏南京210009)摘要：有限层法是一种对空间某一方向进行数值离散，而在其余两方向采用连续函数的半数值半解析方法，该方法能有效地将三维问题简化为一维问题求解，从根本上解决了常用数值计算方法在模拟三维地下水运动时存在的计算工作量大，占用内存多、耗时大等缺点。本文基于有限层法的优点，推导了以伽辽金法结合贝塞耳函数为基础的层状非均质各向异性承压含水层的稳定流有限层方程，并编制了相应的计算程序。通过对两个经典算例的数值解与解析解对比分析，验证了本文方法的正确性。关键词：有限层法；各向异性；承压含水层；稳定流；降深中图分类号：P641文献标识码：A文章编号：宏博（1981－），男，硕士研究生，主要研究方向为环境岩土工程及数值分析。AnalysisofSteadyGroundwaterFlowintheConfinedAquiferbyFiniteLayerMethodZHUHong-bo,WANGXu-dong,ZAIJin-min(CollegeofCivilEngineering,NanjingUniversityofTechnology,Nanjing210009,China)Abstract:Thefinitelayermethod(FLM)isaquasi-numericalandquasi-analyticmethod,whichistodiscretizeonedimensionofthespatialdomainusingfiniteelements,approximatingvariationsintheothertwodimensionsusingcontinuousfunction.Themethodthereforereducesthree-dimensionalproblemstoone-dimensionalproblemseffectively.TheFLMcanresolvethedisadvantagesradically,suchasthelargeamountofworks,themuchmemoryandmoretime,etal,byusingcommonnumericalmethodstosimulatethree-dimensionalgroundwaterflows.TakingadvantagesoftheFLM,finitelayerformulationsarederivedforcalculatingthegroundwaterdrawdownintheheterogeneousandanisotropicconfinedaquifer,whichbasedonGalerkinmethodandBesselfunction,andacomputerprogramusingFORTRANlanguageisdeveloped.Twonumericalexamplesarepresentedandcomparedwithanalyticalsolutionstodemonstratethevalidityoffinitelayermethod.Keywords:finitelayermethod;anisotropic;confinedaquifer;steadygroundwaterflow;drawdown地下水运动由于受含水层的非均质、各向异性和边界条件等因素的影响，其运动状态十分复杂，对于非均质、各向异性问题很难用解析方法进行求解，必须借助于数值模拟方法求其近似解［1-2］。国内外学者不断尝试寻找有效的数值方法用于地下水模拟计算，薛禹群等将多尺度有限元法应用于地下水模拟，节省了计算量，提高了计算精度［3］。朱学愚等采用边界元法对二维稳定流进行了分析计算［4］。但是，数值法在进行三维问题的计算分析时，常常会遇到计算分析占用内存多、计算工作量大等问题。Cheung于1968年首先提出了有限条法，并成功地将有限条法应用于结构工程问题的分析［5］。该方法能有效地将三维问题简化为一维问题求解，从而大大减少单元数、节省计算工作量和内存的需要量。有限层法正是利用了有限条法的特点，进一步发展和完善起来的。它是一种对空间某一方向进行数值离散，而在其余两方向采用连续函数的半数值半解析方法。Small和Booker等运用了有限层法对土的固结与流变耦合问题做了较为系统的研究［6］；Stanley、Myron等将有限层理论应用于地下水分析［7］；Southcott也采用该方法对竖向荷载下的桩土相互作用问题进行了探讨［8］；支喜兰、王秉纲利用有限层方法对弹性地基上受对称荷载作用的圆形厚板进行了力学分析［9］；宰金珉、王旭东等将有限层法应用于土与结构物的共同作用分析［10-11］。有限层法已在岩土工程数值计算中得到了广泛的应用。本文利用有限层法的优点，以伽辽金法结合贝塞耳函数为基础，推导了承压含水层稳定流有限层方程，编制了相应的计算程序，通过算例验证了该方法的正确性。1承压含水层稳定流模型的数学描述假设层状非均质、各向异性、等厚承压含水层，初始水位水平，定水头边界条件，如图1所示。根据水量均衡原理和达西定律，柱坐标下承压含水层稳定流的数学模型可表示为10,00rRrzssKrKWr,zrRrrrzzsrz()≤≤(,)(1)式中：s(r,z)为降深；Kr、Kz分别为r径向、z方向的渗透系数；W(r,z)为源汇项，补给源取正值，排泄汇取负值；R为影响半径。R承压含水层完整井隔水层含水层隔水层srzR定水头边界定水头边界Q定水头边界zR定水头边界隔水层含水层R隔水层srQ图1承压含水层完整井物理模型Fig.1Configurationoffullypenetratingwellphysicalmodelinconfinedaquifer2有限层格式推导有限层计算简图如图2所示。圆形计算区域径向r的计算尺寸取为R，竖向z的计算尺寸取为c。同时沿z轴方向将承压含水层离散成L个层元。层元1z含水层厚度cRr节面2节面1层元l层元L节面L+1O图2有限层计算简图Fig.2Finitelayermodels2.1试探函数设降深在径向用解析函数表示，而z方向采用有限层元离散，层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数111)()(,~LjMmjmmjzNrAΦzrs(2)式中：zrs,~为降深试探函数；Φmj为待求系数；rkrAmm0J)(为满足定水头边界的第一类零阶贝塞耳函数，其中Rkmm/0，0m为零阶贝塞耳函数的第m个正零点；M为级数项数；L为有限层层元数。Nj(z)为线性插值基函数，其形式如下：11111111//0jjjjjjjjjjjjjzzzzzzzNzzzzzzzzzzzz≤≤≤≤≤≥或(3)2.2伽辽金方程根据加权余量法原理，可得有限层法伽辽金方程1,1,2,,0dd~LizrrWsLDim其中(4)式中：D为含水层计算区域；W为源汇项；sL~为微分算子，zsKzrsrKrrsLzr~~1~；im为权函数，)(,rAzNzrmiim；mAr为第一类零阶贝塞耳函数；Ni(z)同样为线性插值基函数，其形式同式(3)。将含水层离散为L个层元，则式(4)可改写为1,1,2,,0dd~1LizrrWsLLeDimee其中(5)式中：De为层元计算区域；We为层元源汇项。将式(2)代入式(5)得1111dd01,2,,1eLLMeeemjmjimDejmLΦANWNArrziL其中(6)式中：ejNz，eiNz为层元线性插值基函数。2.3层元有限层方程取一典型有限层层元l，其节面分别为l，l+1，层厚cl，如图3所示。rz层元lcl节面l节面l+1RO图3典型有限层层元Fig.3Finitelayerelementmodels利用贝塞耳函数Am(r)和Am′(r)的正交性，可得])(J[20d)(J)(J021200RmmmmmRkRmmrrrkrk(7)式中：Rkm1J为一阶贝塞耳函数。故，当mm时，由式(6)得：111101200200ddddddd0,1lllllllMZReemmjrijmZjlmeRjecmjzimeeZRjimjzmZZReejmZAΦKNNzrArrrNΦKNArrzNNΦKzArrzzWNArrzil,l其中(8)结合基函数的性质，由式(8)积分可得到层元l的有限层方程：0111,1,11,,meemllmellellellellllqqgggg(9)式中：emijg为层元渗透矩阵系数；emjq为水量矩阵系数。jicKckKRjicKckKRgllzlmlrmllzlmlrmemij6)(])(J[23)(])(J[22201222012(10)对于完整井抽水：1000dJd2πJ(0),14πllzRemjjmzlrqNzQkrrrrcQjll其中(11)其中，Q为抽水率，lcQQ/，Q为井流量。对于半球状点源：000Jd2π1J(0),12πRemjmrqQkrrrrQjll其中(12)式(9)用矩阵可简化表示为：0eemmmGΦq(13)式中：emG为层元渗透矩阵；emq为水量矩阵；mΦ为待求系数矩阵。2.4整体有限层方程由层元有限层方程组装可得整体有限层方程：0mmmGΦq(14)式中：mG为整体渗透矩阵，1LemmeGG；T121jLLmmΦ，，，，，，；mq为整体水量矩阵，1Lemmeqq。整体矩阵mG由层元矩阵叠加而成，叠加方法如图4所示。L,L+1L+1,L+1x3223x21x12x+22xx2211x1111222LL+xLL,Lx332+ll+1,l+1lxl+1,llxl,l+1l,lxl-1+lxl,ll+1,l+1l+1xxxx+L+1,LLx图4[G]m矩阵的叠加Fig.4Matrixassemblyfor[G]m3有限层方程求解式(14)整体有限层方程为L+1阶线性代数方程组，可采用追赶法［12］求解方程，求得未知系数Φmj后，由式(2)计算任意位置的降深。上述计算过程通过编制相应的计算程序实现。4数值算例及分析4.1承压含水层中完整井假设承压含水层厚c=20m，为了便于对比分析，渗透系数取Kr=Kz=5m/d，计算区域取为R=2500m的圆形区域。在计算区域中心布置一完整井，抽水量Q=400m3/d，如图1所示。采用有限层法模拟计算完整井抽水引起的降深。为了验证承压含水层稳定流有限层分析方法的正确性，引入无限承压含水层稳定流完整井解析解［2］进行对比分析。0.00.10.20.30.40.50.60.70.80.91.08765432102sKc/Qr/R有限层解解析解图5