潜水稳定井流的剖面二维数值模拟方法.pdf

第4 2卷 第4期2 0 2 3年 7月 地 质 科 技 通 报B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y V o l.4 2 N o.4J u l.2 0 2 3王旭升,谢永桦,陈崇希.潜水稳定井流的剖面二维数值模拟方法J.地质科技通报,2 0 2 3,4 2(4):2 7-3 6.W a n g X u s h e n g,X i e Y o n g h u a,C h e n C h o n g x i.S e c t i o n a l 2 D n u m e r i c a l m o d e l l i n g m e t h o d f o r s t e a d y s t a t e w e l l-f l o w i n a n u n c o n-f i n e d a q u i f e rJ.B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y,2 0 2 3,4 2(4):2 7-3 6.基金项目:国家自然科学基金项目(4 1 9 7 2 2 6 3;4 1 7 7 2 2 4 9)作者简介:王旭升(1 9 7 4),男,教授,博士生导师,主要从事地下水动力学和水文模型研究工作。E-m a i l:w x s h c u g b.e d u.c n E d i t o r i a l O f f i c e o f B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y.T h i s i s a n o p e n a c c e s s a r t i c l e u n d e r t h e C C B Y-N C-N D l i c e n s e.潜水稳定井流的剖面二维数值模拟方法王旭升1 a,谢永桦1 b,陈崇希2(1.中国地质大学(北京)a.水利部地下水保护重点实验室(筹);b.水资源与环境学院,北京 1 0 0 0 8 3;2.中国地质大学(武汉)环境学院,武汉 4 3 0 0 7 8)摘 要:经典D u p u i t井流模型与考虑入渗的改进D u p u i t井流模型,都受到D u p u i t假定的影响,可能存在系统误差。建立反映三维流或轴对称二维流特性的潜水井流数值模型,是检验D u p u i t模型可靠性的重要手段。提出一种模拟潜水稳定井流的剖面二维数值模拟方法,通过参数转换把柱坐标系的渗流方程变换为等效的直角坐标系方程,利用MO D F L OW方体网格有限差分模型实现剖面二维流场模拟。数值模型以抽水井的井中水位为已知条件,渗出面的排水量按照D a r c y定律差分公式计算,潜水面则通过MO D F L OW对干-湿单元的处理加以识别,抽水流量经水均衡计算得到。通过采用精细化网格建立典型案例模型,获得模拟精度很高的结果,使反算抽水井流量的相对误差不超过0.2%。以此检验D u p u i t井流模型,发现解析公式得到的水位线总体与数值模拟结果一致,仅在抽水井附近由于没有考虑水跃而偏低,且误差受到含水层渗透系数各向异性的影响。在有入渗的情况下,分水岭附近的渗流违反D u p u i t假定。然而,改进的D u p u i t井流公式计算的分水岭水位相对误差低于1%。这一数值模拟方法简单实用,但也受到MO D F L OW本身局限性的约束。关键词:D u p u i t假定;抽水井;入渗补给;渗出面;潜水面;各向异性;有限差分法2 0 2 3-0 1-1 3收稿;2 0 2 3-0 4-1 4修回;2 0 2 3-0 4-1 9接受中图分类号:P 6 4 1 文章编号:2 0 9 6-8 5 2 3(2 0 2 3)0 4-0 0 2 7-1 0d o i:1 0.1 9 5 0 9/j.c n k i.d z k q.t b 2 0 2 3 0 0 2 4 开放科学(资源服务)标识码(O S I D):S e c t i o n a l 2 D n u m e r i c a l m o d e l l i n g m e t h o d f o r s t e a d y s t a t e w e l l-f l o w i n a n u n c o n f i n e d a q u i f e rW a n g X u s h e n g1 a,X i e Y o n g h u a1 b,C h e n C h o n g x i2(1 a.K e y L a b o r a t o r y o f G r o u n d w a t e r C o n s e r v a t i o n o f M i n i s t r y o f W a t e r R e s o u r c e s(i n p r e p a r a t i o n);1 b.S c h o o l o f W a t e r R e s o u r c e s&E n v i r o n m e n t,C h i n a U n i v e r s i t y o f G e o s c i e n c e s(B e i j i n g),B e i j i n g 1 0 0 0 8 3,C h i n a;2.S c h o o l o f E n v i r o n m e n t S t u d i e s,C h i n a U n i v e r s i t y o f G e o s c i e n c e s(Wu h a n),Wu h a n 4 3 0 0 7 8,C h i n a)A b s t r a c t:O b j e c t i v eB o t h t h e c l a s s i c a l D u p u i t m o d e l a n d t h e m o d i f i e d D u p u i t m o d e l i n c l u d i n g i n f i l t r a t i o n a r e i n f l u e n c e d b y D u p u i ts a s s u m p t i o n a n d h a v e p o t e n t i a l s y s t e m a t i c e r r o r s.B u i l d i n g a n u m e r i c a l m o d e l o f w e l l f l o w i n a n u n c o n f i n e d a q u i f e r b y c h a r a c t e r i z i n g t h e t h r e e-d i m e n s i o n a l o r a x i s y mm e t r i c t w o-d i m e n s i o n a l(2 D)f l o w i s a n e s s e n t i a l a p p r o a c h t o v e r i f y t h e p e r f o r m a n c e o f D u p u i t-t y p e m o d e l s.M e t h o d sI n t h i s s t u d y,a 2 D n u m e r i c a l m o d e l i s p r o p o s e d f o r t h e s t e a d y s t a t e w e l l-f l o w i n a n u n c o n f i n e d a q u i f e r,i n w h i c h t h e c o n t r o l e q u a t i o n o f s e e p a g e i n t h e c y l i n d r i c a l c o o r d i n a t e s i s t r a n s f o r m e d e q u i v a l e n t l y t o t h e C a r t e s i a n c o o r d i n a t e s t h r o u g h p a r a m e t e r t r a n s f o r m a t i o n,a n d t h e s e c t i o n a l 2 D m o d e l l i n g i s i m p l e m e n t e d v i a t h e MO D F L OW f i n i t e-d i f f e r e n c e g r i d o f c u b i c b l o c k s.I n t h e n u m e r i c m o d e l,t h e w a t e r l e v e l i n t h e p u m p i n g h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 w e l l i s a g i v e n c o n d i t i o n,t h e f l u x a c r o s s t h e s e e p a g e f a c e i s e s t i m a t e d b y d i f f e r e n c e e q u a t i o n a c c o r d i n g t o D a r c ys l a w,t h e p h r e a t i c s u r f a c e i s i d e n t i f i e d b y t h e t r e a t m e n t o f d r y a n d w e t c e l l s i n MO D F L OW,a n d t h e p u m p i n g r a t e i s d e t e r m i n e d f r o m t h e w a t e r d e b u g c a l c u l a t i o n.R e s u l t sF i n e g r i d s a r e c o n s t r u c t e d i n-n u m e r i c a l m o d e l s o f t y p i c a l c a s e s t o o b t a i n h i g h-p r e c i s i o n r e s u l t s,i n w h i c h t h e r e l a t i v e e r r o r o f t h e b a c k-w a r d s e s t i m a t e d p u m p i n g r a t e i s n o m o r e t h a n 0.2%.T h i s n u m e r i c a l m o d e l i s u s e d t o c h e c k t h e D u p u i t-t y p e w e l l-f l o w m o d e l s.A s i n d i c a t e d,t h e g r o u n d w a t e r l e v e l e s t i m a t e d f r o m t h e a n a l y t i c a l f o r m u l a s g e n e r-a l l y a g r e e s w e l l w i t h t h e n u m e r i c a l m o d e l l i n g r e s u l t s,e x c e p t t h a t n e a r t h e w e l l,w h e r e t h e a n a l y t i c a l s o l u-t i o n u n d e r e s t i m a t e s t h e g r o u n d w a t e r l e v e l d u e t o i g n o r i n g t h e w a t e r j u m p a n d t h e e r r o r s d e p e n d o n t h e a n i-s o t r o p i c p e r m e a b i l i t y o f t h e a q u i f e r.Wh e n i n f i l t r a t i o n e x i s t s,t h e f l o w i n t h e v i c i n i t y o f w a t e r s h e d d o e s n o t f o l l o w D u p u i ts a s s u m p t i o n.H o w e v e r,t h e e s t i m a t e d g r o u n d w a t e r l e v e l o n t h e w a t e r s h e d b y m o d i f i e d D u-p u i t w e l l-f l o w e q u a t i o n h a s a l o w l e v e l o f r e l a t i v e e r r o r,w h i c h i s l e s s t h a n 1%.C o n c l u s i o nT h i s n u m e r i-c a l m e t h o d i s s i m p l e a n d p r a c t i c a l h o w e v e r i t i s a l s o i n f l u e n c e d b y l i m i t a t i o n s i n MO D F L OW.K e y w o r d s:D u p u i ts a s s u m p t i o n;p u m p i n g w e l l;i n f i l t r a t i o n r e c h a r g e;s e e p a g e f a c e;p h r e e t i c s u r f a c e;a n i-s o t r o p i c;f i n i t e d i f f e r e n c e m e t h o dR e c e i v e d:2 0 2 3-0 1-1 3;R e v i s e d:2 0 2 3-0 4-1 4;A c c e p t e d:2 0 2 3-0 4-1 9 地下水动力学中的井流模型具有轴对称特征,使用柱坐标系进行描述。经典D u p u i t井流模型1引入D u p u i t假定,把地下水流简化为单纯的径向一维流问题,便于获得控制方程的解析解。D u p u i t井流模型忽略了降水入渗的影响。为此,有研究者提出了改进模型2-4,将入渗补给考虑在径向一维流问题的控制方程中。这些模型能够提供抽水井周围潜水面形态的理论公式。轴对称状态的潜水井流实际属于柱坐标系下的二维渗流过程,采用D u p u i t假设是将其简化为径向一维流问题的手段。在含水层的局部地段(如地下分水岭处),地下水的垂向流速可能显著大于水平流速,不满足D u p u i t假定。那么,如何检验D u p u i t假定对潜水稳定井流问题的有效性?这是地下水动力学中一个受到关注的问题。2 0 2 0年陈崇希4在提出考虑入渗的改进D u p u i t模型时,也指出所建立的方程包含D u p u i t假定,可能在地下分水岭和抽水井附近带来误差,但只做了简单定性讨论,未定量分析D u p u i t假定对解析结果的影响。检验D u p u i t假定的影响,需要求解柱坐标系的二维潜水井流问题。不过,由于包含潜水面这样的自由面边界,获得解析解的难度很大。潜水井流问题还涉及另一个需要考察的实际因素,就是渗出面。地下水流向抽水井时,井壁水位高于井中水位(即水跃),介于两者之间的一段井壁就是渗出面。这是地表水等势体附近流线和等水头线弯曲状态(不满足D u p u i t假设)所产生的必然结果,其存在性与流场是否稳定无关2,5-7。前苏联学者已经证明5-6:即使存在水跃,D u p u i t模型的流量公式仍然是可靠的(无入渗补给条件),尽管在抽水井附近实际地下水位高于D u p u i t公式计算的结果。对于潜水井流的水跃现象,B o u l t o n7利用L a p l a c e方程的有限差分解法进行了系统的研究。张有龄8在阐述潜水井流问题和水跃现象时,详细引用了杨式德所做的有限差分模型,其方法与B o u l t o n所用方法类似。其他学者也用有限差分法9或有限单元法1 0对此做了数值模拟研究。由于非线性带来的困难,目前尚未研究得到水跃的严格解析解。前人所做的潜水井流数值模拟没有考虑降水入渗补给,而且需要采用柱坐标系进行网格离散,不容易在现有的通用数值模拟工具中实现。如果能够在直角坐标系建立与柱坐标系等价的剖面二维渗流数值模型,将为模型网格离散以及处理非均质含水层、入渗补给、渗出面等水文地质要素带来便利。目前国际上流行的三 维地下水流 有 限 差 分 模 拟 程 序MO D F L OW1 1采用的是直角坐标系,可在一些专业软件进行快速建模。笔者拟利用控制方程的等效原理,在直角坐标系下采用等效方法模拟潜水稳定井流剖面二维流场,与D u p u i t模型或改进D u p u i t模型解析解进行对比,对有关问题进行讨论。1 直角坐标系模拟轴对称渗流方法1.1两种坐标系的渗流方程在直角坐标系中(图1-a),二维地下水稳定渗方程1 2的一般形式为:xKxHx +zKzHz =0(1)式中:H为水头;x,z分别为水平与垂直坐标,而Kx和Kz分别为对应坐标方向的渗透系数。该方程成立的条件是渗透张量的主轴与坐标轴一致,且流场内部没有源汇项。潜水面作为二维流场的上边82 第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法界,其稳态约束条件可表示为1 2:KxHxzw tx-KzHz+w=0(2)式中:zw t为潜水面的高度(随x坐标变化);w为直角坐标系入渗补给强度。在柱坐标系中(图1-b),轴对称渗流方程1 3的一般形式为:1rrKhrHr +zKvHz =0(3)式中:r为径向距离;Kh和Kv分别为水平和垂直方向的渗透系数。这种柱坐标系下的潜水面稳态边界条件表示为:KhrHrzw tr-KvrHz+r=0(4)式中:为柱坐标系下的入渗补给强度。1.2等效变换方法已有一些研究者提出了在直角坐标系等效模拟轴对称地下水流的方法。S a m a n i等1 4将直角坐标系的水平距离变换为x=l n(r)建立MO D F L OW模型,然后按照这种对数变换,保持水平渗透系数不变,而对垂向渗透系数实施变换,产生等效结果。L a n g e v i n1 5提出另外一种在MO D F L OW中转换参数的方法,他发现MO D F L OW用对数加权平均法计算的块体单元界面等效渗透系数正好与D u-p u i t公式的效果一致。L o u w y c k等1 6建立x方向均匀划分的MO D F L OW模型,然后通过参数变换获得等效模拟结果,将水平渗透系数变换为:图1 直角坐标系(a)和柱坐标系(b)的潜水面示意图(z为相对高度;x为横向距离;r为径向距离;后同)F i g.1 S c h e m a t i c s o f t h e p h r e e t i c s u r f a c e i n t h e c a r t e s i a n(a)a n d c y l i n d r i c a l(b)c o o r d i n a t e s y s t e m s Kx=2 Khl n(rj+1/2/rj-1/2)(5)式中:j为 单 元 体 在x方 向 的 离 散 坐 标;rj+1/2,rj-1/2分别为单元外侧和内侧的径向距离。这种方法实际上也是取x与l nr成正比。本研究从控制方程等价的角度建立变换关系,即保持直角坐标系的水平距离x与柱坐标系的径向距离r相同,通过参数变换,使式(1)与式(3)等价,而式(2)与式(4)等价。满足这种等价关系的变换式为:Kx=2 x Khy,Kz=2 x Kvy,w=2 x y,x=r(6)式中:乘积因子2 x为径向距离r处完整过水断面的周长;肀y为模型在x,z正交方向的宽度。在式(6)中变量与肀y相除,可以保持参数的量纲不变。从物理意义上来讲,肀y是剖面模型所代表的过水断面单位宽度。MO D F L OW模型默认单元都是三维方体,因此必然会有y方向的宽度肀y,将其设定为一个任意单位宽 度(例如1 m)代入式(6)即可。L a n g e v i n1 5提出的方法本质也是如此,但他没有从控制方程等价的角度进行解释,也没有处理渗出面。1.3渗出面的处理渗出面是地下水流的特殊边界,在忽略流速水头的情况下,其水头与渗出点的位置高度相等。对于有入渗补给的潜水稳定井流模型,渗出面既可能出现在抽水井的井壁,也可能出现在外周含水层与大气的交界面上。假设井壁和含水层外周均沿铅直方向延伸,分别位于r=rw(井半径)和r=R(含水层外周半径)处,则渗出面的边界条件可表示为:H(x,z)=z,x=rw,hwz(hw+Lw)(7)H(x,z)=z,x=R,hRzhw)z92h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 或 H(R-x,zhR)z(9)Vx=0,H(rw+x,zhw)z或 H(R-x,zhR)z(1 0)式中:x0,表示偏移一小段距离。采用偏移小段距离的水头大小来判断是否处于渗出面附近,符合数值模型中边界单元的特性。以井壁附近的单元为例,如图2-a所示,第k层和k+1层的相邻单元(宽度为肀x)计算水头分别为H(k)和H(k+1),它们均与x=rw+0.5肀x处的潜水面高度zw近似相等,但是只有第k层的水头满足H(k)z(k)而处于渗出面附近。这里有x=0.5肀x。图2-a中,第k层单元与渗出面接触,并且排泄的侧向流量Qx按照D a r c y定律计算,即将式(9)的差分格式与单元界面的面积相乘,计算为:Qx=-Kx2(Hk-zk)xzy,Hkzk(1 1)Qx为负值表示流向与x轴相反。图2-a中第k+1层单元格点高于潜水面,符合式(1 0)所述条件,即Qx=0。在含水层外周,可以证明存在类似关系。因此,可以将渗出面的边界条件替换为第三类边界,即:Qd=Cd(Hk-zk),Hkzk;Qd=0,Hkzk(1 2)式中:Qd为模型单元向外部排泄的流量;Cd称为模型单元的排水系数。根据式(1 1),并把式(6)代入,可得:Cd=2Kxzyx=4 x Khzx,x=rw或x=R(1 3)a.井壁渗出面相邻单元;b.井孔饱水带相邻单元;rw为抽水井半径;x为单元长度;z为单元厚度;y为单元宽度;Qx为侧面流量;zw为井中水位;z(k)为单元中点高度;H(k)为单元格点水头图2 含水层网格单元与井孔的关系示意图F i g.2 S c h e m a t i c s o f t h e r e l a t i o n s h i p b e t w e e n g r i d b l o c k s o f t h e a q u i f e r a n d t h e w e l l 显然,排水系数的取值与模型网格剖分有关,并非自然界含水层的固定参数。实际上,井中水位以下的井壁(水头值为hw),也可以用类似的方法处理。如图2-b所示,设第k层有一个模型单元与井孔单元相邻,其宽度为肀x,2个单元交界面上的水头固定为hw,则侧向流量的差分公式为:Qx=-Kx2(Hk-hw)xzy(1 4)式(1 4)与处理渗出面的式(1 1)是类似的。因此也可使用类似式(1 2)的公式计算单元的排泄流量:Qd=Cd(Hk-hw)(1 5)其中排水系数Cd的计算方法与式(1 3)相同。2 模型设计与MO D F L OW用法2.1含水层模拟情景本研究模拟均质含水层中的潜水井流问题,以便与D u p u i t模型和改进的D u p u i t模型4解析解进x,r为距离;R为圆岛半径;LR为外侧边界水跃;H(r,z)为水头分布函数;hw为井中水位;Lw为井壁水跃图3 潜水井流模型示意图F i g.3 D i a g r a mm a t i c m a p o f t h e m o d e l f o r t h e w e l l-f l o w i n a n u n c o n f i n e d a q u i f e r行对比。模型如图3所示,有一半径为R=2 0 0 0 m03 第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法的圆岛状均质含水层,外周边界水位为hR=5 0 m,圆岛中心处有一个半径为rw=0.1 m的抽水井,以恒定流量Qw抽水,最终形成井中水位hw。含水层从隔水底板到地面的总厚度为7 0 m,地下水接受大气降水入渗补给,强度为,最大可达 1 mm/d。地下水流场达到稳定,形成非线性的潜水面形态,并可能在抽水井的井壁发育长度为Lw的渗出面(即水跃值为Lw),在含水层外周边界发育长度为LR的渗出面。按照含水层渗透系数与入渗补给强度的变化,设置不同情景进行模拟。情景编号和参数见表1。表1 模拟情景编号和参数T a b l e 1 C o d e a n d p a r a m e t e r s o f m o d e l l i n g s c e n a r i o s模拟情景水平方向渗透系数Kh/(md-1)各向异性比值Kv/Kh入渗补给强度/(mmd-1)井中水位hw/mA1 010.03 5B1 011.03 5C1 00.1,1 0.01.03 5D10.11.02 02.2离散网格模拟工作是在A q u a v e o公司发布软件G r o u n d-w a t e r M o d e l i n g S y s t e m(GM S)中完成的,构建三维方体有限差分网格(y方向只取一个单元宽度),调用MO D F L OW模块处理水文地质要素。模型网格的初始剖分方式为:在垂向上分为7 0层,每层厚度1 m;水平方向也进行密集划分,总体每个网格单元的宽度为2 m,抽水井本身占一排宽度为0.1 m的单元,其附近含水层单元加密,宽度缩小到0.0 1 m。模型中井孔所在的这一列单元设置为无效单元,即并不直接参与有限差分计算,而是以边界的方式作用到模型中(具体处理见2.3节)。模型最右侧一排单元与含水层外周边界相邻。在模拟过程中,根据渗出面和潜水面分水岭的结果进行模型网格的微调,提高精度。每个含水层单元的水平渗透系数Kx和垂向渗透系数Kv,根据式(6)计算,其中x为单元形心与抽水井中心线的距离。2.3抽水井和外周边界的处理本次模拟并不像解析模型那样把抽水井处理为已知流量的边界,因为抽水井本身也被离散为很多个单元格(类似混合井),事先并不能准确掌握每个单元分配的流量。为此,本研究将抽水井的水位处理为已知条件,即给定hw的数值,然后在模拟形成稳定流之后反求出抽水井的流量。在井壁处,z坐标高于hw的点可能位于渗出面,其水头值与z坐标一致。与井壁相邻的一排模型单元,采用MO D F-L OW中的D r a i n模块1 1处理地下水向抽水井的排泄。高于井中水位的边界单元,排泄量用式(1 2)计算,以zk为D r a i n模块的排水高程。低于井中水位的边界单元,排泄量用式(1 5)计算,以hw为D r a i n模块的排水高程。D r a i n模块的排水系数统一采用式(1 3)计算,取x=rw。调用MO D F L OW运行稳定流模拟之后,再根据模拟结果确定抽水井的流量Qw。它实际上由2个部分构成的:一部分是从渗出面排出的地下水量,一部分是从井水位下方井壁排出的地下水量。因此,抽水井流量的统计可以表示为:Qw=Nk=1Cd(Hk-hw)+Mk=N+1Cd(Hk-zk)(1 6)式中:N为井中水位以下的最高井壁单元层号;M为对应渗出面的最高井壁单元层号。从模拟结果计算Qw的方法,是将井壁相邻的饱和带单元组合设置为一个均衡区,进行编码,调用MO D F L OW的B u d g e t模块计算均衡区流进和流出的水量。没有计算误差的情况下,两者的绝对值应相等,并且与Qw一致。如果潜水面存在分水岭,也可以将分水岭半径范围以内的全部入渗补给水量计为Qw。这2种计算方法的结果不一定相等,只要控制相对误差不超过1%,就可以认为数值模拟结果是可靠的。在模型的外周边界处,渗出面的情况与井壁渗出面是类似的,也用D r a i n模块处理。对高于水位hR的边界 单 元,排 泄 量 用 式(1 2)计 算,以zk为D r a i n模块的排水高程,取x=R计算排水系数Cd。低于hR的边界单元则将水头固定为hR。2.4潜水面的处理潜水面属于自由水面,各点的水头值与其高度相等,这可以用来判断潜水面的位置。与此同时,潜水面还需要满足式(2),作为渗流场的边界条件。MO D F L OW并不直接求解式(2),而是以入渗补给量作为源、建立潜水面所在网格单元的水均衡方程,与饱和带的方程组一起求解。从数值计算的角度来讲,这与求解式(2)是等价的。如果计算水头高于单元的底部,则该单元为“湿单元”,是有效的。否则,就判断为“干单元”,相当于无效单元,不再参加下一次计算,而是将潜水面移入下部的相邻单元1 1,1 7。MO D F L OW的R e c h a r g e模块可以自动将入渗补给加入到最高处的湿单元上。模型入渗强度w根据式(6)计算,其中x为单元形心的r值。模拟过程中,潜水面的形成实际上表现为网格单元的逐渐疏干过程。先将所有不属于边界的网格单元水头设置为模型顶部高程(7 0 m)。然后启动渗流场计算和“干”、“湿”单元判断,反复迭代,逐渐将水头高于潜水面的单元疏干,而余下的网格单元13h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 构成饱和带的渗流空间。迭代过程中可能会出现疏干过快或振荡现象,导致计算不稳定,因此需要对迭代算法进行控制,使疏干过程平稳进行。2.5求解算法和收敛精度MO D F L OW形成有限差分方程组后,可采用多种方法进行求解。本次模拟采用预共轭梯度法(P C G 2模块1 8),这是一种迭代计算效果较高的算法。为保证数值模拟精度,P C G 2模块中水头计算的收敛标准设置为0.0 0 1 m,水均衡计算的收敛标准设置为0.0 1 m3/d。由于潜水面的形成是一个将包气带单元逐步疏干的过程,不宜太快,将P C G 2模块中迭代计算的松弛因子(r e l a x a t i o n p a r a m e t e r)从默认的1.0调整为0.5。3 典型模拟结果分析3.1对比D u p u i t模型解析解表1中的情景A是无入渗、均质各向同性含水层条件,剖面二维渗流模拟结果可以与D u p u i t公式进行对比。图4给出了抽水井附近的模型网格与模拟结果。显然,潜水面的模拟水位高于解析解,在井壁处模拟得到渗出面的顶部高为4 0.5 m,比井中水位3 5.0 m高出了5.5 m,这就是水跃值。圆岛外周的水跃值小于1 m。图4 情景A模拟水位与经典D u p u i t模型解析解的局部对比F i g.4 L o c a l c o m p a r i s o n b e t w e e n t h e m o d e l l e d g r o u n d-w a t e r l e v e l o f t h e s c e n a r i o-A a n d t h e a n a l y t i c a l s o l u t i o n o f t h e c l a s s i c a l D u p u i t m o d e l根据模拟结果,对抽水井流量进行了计算,得到Qw=4 0 5 1.8 m3/d,其中从渗出面排泄的地下水量为3 6 4.2 m3/d,占9.0%。模型的总体水均衡误差(流入、流出之差)绝对值小于0.0 1 m3/d,说明有很高的差分计算精度。解析方法计算流量的D u p u i t公式1为:Qw=Kh(h2R-h2w)l n(R/rw)(1 7)把含水层参数和边界水位代入式(1 7),得到Qw=4 0 4 4.6 m3/d。与之相比,数值模型得到的流量偏高0.2%。D u p u i t模型的流量公式对三维流也是适用的,而且不受水跃的影响,这一点已经得到了严格的证明5。因此,上述0.2%的相对误差代表了数值模型在流量计算上的误差水平。3.2对比改进的D u p u i t模型解析解表1中的情景B是有入渗、均质各向同性含水层条件,剖面二维渗流模拟结果可以与改进的D u-p u i t公式4进行对比(图5)。宏观上看,模拟的潜水面曲线与解析解很接近,在分水岭处几乎重合(图5-a)。在抽水井附近,解析解则明显低于模拟结果(图5-b),渗出面产生的水跃值Lw7.1 m。对于情景B,抽水井流量与井中水位的解析关系为4:Qw=Kh(h2R-h2w)l n(R/rw)+R22 l n(R/rw)(1 8)把参数和边界条件代入式(1 8),得到抽水井流量的解析解为Qw=4 6 7 9.0 m3/d,高于无入渗的情景。数值模型反算的结果为Qw=4 6 8 7.9 m3/d,其中渗出面流量5 1 0.7 m3/d,占1 0.9%。模型的水均衡误差仍然小于0.0 1 m3/d,数值模拟得到的抽水流量比解析解偏大0.2%,相对误差很小。这种一致性也反映在分水岭上,解析解确定的分水岭在x=1 2 2 0.4 m处,数值模拟的分水岭在x=1 2 2 1 m处,相对误差仅有0.0 5%,而且数值模拟的分水岭水位是5 0.5 2 m,与解析解一致。分水岭处地下水主要呈现垂向流动,明显违反D u p u i t假设。然而,从本例对比结果来看,引入D u p u i t假设并没有影响解析解计算分水岭水位的准确性。a.整体范围;b.抽水井附近图5 情景B模拟水位与改进D u p u i t模型解析解的对比F i g.5 C o m p a r i s o n b e t w e e n t h e m o d e l l e d g r o u n d w a t e r l e v e l o f s c e n a r i o-B a n d t h e a n a l y t i c a l s o l u t i o n o f t h e m o d i f i e d D u p u i t m o d e l3.3有入渗的各向异性含水层情景表1中的情景C和D都是各向异性含水层,可以考察含水层各向异性特征对潜水井流的影响。改23 第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法进D u p u i t模型4并没有包含对各向异性的处理。情景C与情景B的基本条件是相同的,唯一区别在于情景C的垂向渗透系数与水平渗透系数不同。根据模拟结果(图6),各向异性对水位的影响主要表现在抽水井附近(图6-a),而距离大于1 0 0 m时的水位线对各向异性比(Kv/Kh)不敏感(图6-b)。分水岭位置和水位的模拟结果与情景B几乎一致,说明各向异性对抽水井流