《工程勘察》承压含水层非稳定流有限层分析-修改稿

承压含水层非稳定流有限层分析诸宏博1,2王旭东1宰金珉1(1.南京工业大学岩土工程研究所，江苏南京210009；2.浙江省建筑科学设计研究院，浙江杭州310012)摘要：推导了以伽辽金法结合傅立叶函数为基础的各向异性承压含水层非稳定流有限层方程，编制了相应的计算程序。通过算例的有限层解与解析解的对比分析，在验证有限层法正确性的基础上，深入探讨了计算区域、计算项数等因素对计算结果的影响，得到了计算区域尺寸影响有限层解的收敛速度等结论，最后给出了一定精度条件下计算项数选用的近似估算公式。关键词：有限层法；承压含水层；非稳定流；降深0引言有限层法是一种对某一方向进行数值离散，而在其余两方向采用连续函数的半数值半解析方法，这种方法能有效地将三维问题简化为一维问题求解，从而大大减少单元数，节省计算工作量和内存的需要量。基于这些特点，Cheung在1968年首先提出了有限条法，并成功地将有限条法应用于结构工程问题的分析［1］。Small和Booker等运用了有限层法对层状弹性地基和土的固结问题做了较为系统的研究［2-4］。A.M.oskoorouchi和B.Novrouzian将有限条方法推广应用于岩土工程分析中，对基坑开挖、隧道和土坝等问题进行了研究分析［5］。Stanley,Myron等将有限层理论应用于地下水分析［6］。张问清、赵锡宏、宰金珉等采用有限层方法探讨了任意力系作用下的层状弹性半空间问题［7］。宰金珉、王旭东等将有限层法应用于土与结构物的共同作用分析［8-9］。地下水计算方法主要有解析法和数值法。解析法可以将描述地下水运动的各种物理量反映在一个表达式中，使用方便，但自然界地下水运动状态复杂，往往与解析解的假设条件相差甚远，其应用范围也受到一定的限制。对于一些复杂的地下水问题，就必须借助于数值模拟方法求其近似解。目前，地下水数值分析中主要有有限差分法、有限元法、边界元法等数值方法［10-12］。数值法能较为方便地处理解析解难以解决的复杂边界等问题，因此，得到广泛的应用。但是，数值法计算工作量大，尤其是在考虑三维模型的计算分析时，占用内存多。国内外学者不断尝试寻找有效的数值方法用于地下水模拟计算，薛禹群，叶淑君，谢春红等将多尺度有限元法应用于地下水模拟，节省了计算量提高了计算精度［13］。本文吸取有限层法的优点，以伽辽金法结合傅立叶函数为基础，推导了承压含水层非稳定流有限层方程，编制了相应的计算程序，在验证方法的正确性基础上，探讨了计算区域、计算项数等因素对计算结果的影响。1承压含水层非稳定流模型的数学描述假设均质、各向异性、等厚的有限区域承压含水层，在区域中心(x0,y0)处有一个以恒定流量Q抽水的完整井，初始水位水平，定水头边界条件，如图1所示。根据质量守恒定律和达西定律，承压含水层非稳定流的数学模型用降深可表示为：--(,)0(,,,)00(,,,)00(,,,)00xaxaybybxyzsssssKKKRx,y,ztSxxyyzztt,axa,bybsxyztt,axa,bybsxyztt,bybsxyztt,axa＞≤≤≤≤≤≤≤≤＞≤≤＞≤≤(1)式中，s(x,y,z,t)为降深，Kx、Ky、Kz分别为x，y，z方向的渗透系数，Ss为单位贮水系数，R(x,y,z,t)为源汇项，补给源取正值，排泄汇取负值。(1)式可用微分算子表示为：0RsL(2)式中，tsSzsKzysKyxsKxsLszyx，L为微分算子。───────收稿日期:基金项目：国家自然科学基金项目（50278042）；江苏省高校自然科学研究计划项目（03KJB560044）作者简介：诸宏博（1981－），男（汉族），浙江杭州人，硕士研究生，主要研究方向为环境岩土工程及数值分析。E-mail:[e-mail]层元iQaboxyz(x0,y0)隔水层含水层隔水层cb1a图1承压含水层物理模型与有限层层元Fig.1Configurationofconfinedaquiferphysicalmodelandfinitelayermodels2有限层格式推导有限层计算简图如图1所示。x和y方向的计算尺寸分别取为2a和2b，沿z轴方向将承压含水层离散成L个层元。2.1试探函数设降深在x和y方向分别用解析函数表示，而z方向采用有限层元离散，层元内降深线性分布。则试探函数可表示为下列分离变量的函数项级数：1111,,,(,)()LMNmnjmnjjmnsxyztΦtAxyNz(3)式中，tzyxs,,,~为降深试探函数；Φmnj(t)为待求系数；(,)sinsinmnmnAxykxakyb为满足定水头边界的已知函数，bnkamknm2/,2/；M，N为级数项数；L为有限层层元数。Ni(z)为线性插值基函数，其形式如下：11111111//0iiiiiiiiiiiiizzzzzzzNzzzzzzzzzzzz≤≤≤≤≤或≥(4)2.2伽辽金方程根据加权余量法原理，可得有限层法伽辽金方程：1,1,2,,0~LidzdydxRsLDinm其中(5)式中,inm为权函数,),(,,yxAzNzyxnmiinm。将含水层离散为L个层元，则(5)式可改写为：1,1,2,,0~1LidzdydxRsLLeDinme其中(6)将(3)式代入(6)式得：1111101,2,,1eLLMNeeemnjmnjimnDejmnLΦtANRNAdxdydziL其中(7)2.3层元有限层方程取一典型有限层层元l，其节面分别为l，l+1，层厚cl，如图2所示。a节面l+1lc层元l节面lxzyaobb图2典型有限层层元Fig.2Configurationoffinitelayerelementmodels由Amn(x,y)和Am′n′(x,y)的正交性，可得：0sinsinaammmmammxdxkxk(8)0sinsinbbnnnnbnnydykyk(9)故，当mm和nn时，由(7)式得：1112121122202lllllllMNZabeemnmnjxijmnZabjlmnZabeemnmnjyijmnZabeabjecmnjzimnabeeZjimnjzmnZAΦKNNdzAdxdyxAΦKNNdzAdxdyyNΦKNAdxdyzNNΦKdzAdxdyzz1120,1llllababZabmnjeesijmnZabZablejmnZabΦSNNdzAdxdytRNAdxdydzil,l其中(10)由基函数的性质可得：631llZZllejeijicjicdzNN(11)由(11)式可分别得到m、n阶的层元渗透矩阵系数emnijg、贮水矩阵系数emnijb和水量矩阵系数emnijq：aabbmnZZejeizaabbmnZjeizaabbmnmnZZejeiyaabbmnmnZZejeixemnijdydxAdzzNzNKdydxANzNKdydxAyAzdNNKdydxAxAzdNNKgllllll2202222111(12a)121lliaabbmnZZejeisemnijdydxAdzNNSbll(12b)aabbmnejlZZemnjdzdydxANRqll1(12c)式中，i=l,l+1，j=l,l+1。由(12a)式、(12b)式积分可得：jicKckKkKabjicKckKkKabgllzlnlymlxllzlnlymlxemnij632222(13a)jiScabjiScabbslslemnij63(13b)对于位于区域中心处的完整抽水井，由(12c)式积分可得：1,,sinsin,2,,00001lljbkakyxQcdydxAyxyxQdzNqnmlaabbmnjzzemnjll其中(13c)式中，00,yxQ为抽水率lcyxQyxQ/,,0000。层元l的有限层方程为：0111,1,11,,11,1,11,,mneemnllmnellellellellmnllmnellellellellllqqdtdbbbbgggg(14)用矩阵可简化表示为：0eeemnmnmnmnmndGΦBΦQdt(15)2.4整体有限层方程由层元有限层方程组装可得整体有限层方程：0mnmnmnmnmndGΦBΦQdt(16)式中,LeemnmnGG1,LeemnmnBB1,LeemnmnQQ1，121TjLLmnmnΦ。整体矩阵mnG、mnB由层元矩阵叠加而成，叠加方法如图3所示。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图3[G]mn和[B]mn矩阵叠加Fig.3Matrixassemblyfor[G]mnand[B]mn3有限层方程求解对于非稳定流来说，降深是时间t的函数，则()mnjΦt也是时间t的函数，采用差分法对时间进行离散得：1kkmnmnmndΦΦΦdtt(17)假设在时间k和k+1之间，水头呈线性变化，则有：11kkkmnmnmnΦΦΦ(18)式中，kmnΦ和1kmnΦ分别为k时刻和k＋1时刻的值，t为时间步长，为加权因子，取值范围在0～1之间。将(18)式代入(16)式可得有限层方程的差分格式：11()1(1)kmnmnmnkkmnmnmnmnGBΦtBGΦQt(19)在(19)式中，当0时，即为显式差分格式；1时，为隐式差分格式；5.0时，为Crank－Nicolson差分格式，其中，隐式差分格式和Crank－Nicolson差分格式无条件稳定。(19)式有限层方程的差分格式为L+1阶线性代数方程组，求解方程组可得到未知系数Φmnj(t)，并由(3)式计算k+1时刻的降深。数值计算时初始降深为0，故00mnΦ，计算程序流程图如图4所示，并编制了相应的计算程序，采用追赶法［14］求解方程，以节省计算时间。开始输入计算参数计算单元矩阵[G]、[B]、[Q]求解方程得{}值mn输出计算值计算{}的初值结束mn0k+1m、nM、NNYk+1时间步数NY时间步数=T/Tm=1,2,3...M;n=1,2,3...Neee形成整体矩阵[G]、[B]、[Q]计算降深值kmn{}={}mnk+1图4有限层计算程序流程图Fig.4Programflowgraphoffinitelayermethod4数值算例及分析4.1有限承压含水层中抽水井算例：承压含水层厚c=100