大口径辐射井抽水效果的三维可视化渗流有限元分析

大口径辐射井抽水效果的三维可视化渗流有限元分析

0非稳定渗流场数值模拟技术大口径辐射井因其良好的排水效果，已广泛应用于许多节水工程中。但辐射井因辐射管长度、层数和数量的不同,其排水效果差异很大。如何进行辐射井的合理设计,包括辐射井间距、辐射管的几何尺寸(长度、管径和层数)等是应用辐射管进行排水设计的关键。有限元分析是解决该问题的重要手段,但辐射井的径向尺寸很小而横向尺寸很大、数量众多,且空间分布复杂,给有限元模拟其精细渗流场带来一定的难度。目前对堤坝渗控分析中密集排水孔的模拟方法较多。但目前还没有相关文献对辐射井的排水效果和渗流场特性进行数值模拟方面的研究。为了解辐射井的抽排水效果,在充分研究了辐射井工作原理和渗流行为的基础上,提出了辐射井子结构法和子结构渗流开关器等关键技术进行辐射井渗流场的有限元分析方法,并联合应用迭代增量法和求解大型稀疏矩阵的预处理共轭梯度算法进行非稳定渗流场的分析。该方法既解决了辐射井附近渗流场的精细模拟,又解决了非稳定渗流场自由面的迭代问题,可用于多种排水措施和辐射井不同设计方案的计算比较。在此基础上,采用IDL语言开发了三维可视化渗流有限元分析软件GWSS。将该软件应用于台兰河地下水库辐射井抽水过程和水位回复过程的的非稳定渗流场的计算分析中,并利用各观测孔的现场实测值对作者提出的辐射井子结构法的计算理论和编制软件的可靠性进行验证。1辐射井的工作原理和渗透1.1蓝波辐射井系统辐射井排渗系统一般由竖井、辐射管等组成,见图1。实际运行中辐射井常常和虹吸井一起使用,即虹吸辐射井系统。虹吸辐射井系统除竖井和辐射管外还包括虹吸管和水封槽等。本文重点研究辐射井工作原理和方法,有关虹吸辐射井的研究成果将另文发表。辐射管是关键的降水设施,其成井质量要求较高,另外辐射管的布置形式灵活,可以水平布置、也可以斜向布置,可根据地形地势进行选择。1.2辐射管安全边界的计算辐射井系统的工作模式有2种:一种是抽水,另一种是自流排水。对辐射井的渗流行为在算法上可分2种情况:一是自由面穿过辐射管;二是自由面不与辐射管相交。有自由面穿过的辐射管内边界面上的渗流行为如图2所示。1)如果自由面水位低于辐射管底部aa′,则该辐射管不起到排水作用,辐射管失效。因此可以在每个辐射管底部高程处虚构一个数学开关器,详细见文献。进行辐射管计算时,先假设开关器打开,管内边界全部作为可能逸出边界,每一步迭代后对孔内结点压力进行甄别,记孔内边界各结点压力为p。若辐射管底面压力p<0,认为辐射管全部位于非饱和区,此时辐射管完全失效。辐射管全部在非饱和区时,其边界条件数学表达式为kij∂h∂xjni|Γ1=0‚且h<x3。式中:kij为渗透材料的渗透张量矩阵;h为水头;xj为3个方向的坐标;ni为内边界面外法线方向的方向余弦,i为3个坐标主方向;x3为z方向坐标。若辐射管bb′面上的结点压力p>0,认为辐射管全部位于饱和区,此时边界为定水头边界,其数学表达式为kij∂h∂xjni|Γ1>0‚且h=h0。式中,h0为辐射井竖井内的水位。2)如果自由面水位位于aa′或bb′内,则辐射管内边界为可能逸出边界,其数学表达式为kij∂h∂xjni|Γ2≥0‚且h=z0。式中,z0为aa′或bb′逸出面上结点位置高程。2计算非稳定渗透流的理论2.1土体压缩和水的贮存水量的定解非均质各向异性有自由面稳定渗流的基本微分方程如下:(kijh,j),i+Q=SS∂h∂t。(1)式中:Q为源或汇项(1/s);SS为贮水率(1/m),即单位体积的饱和土体内,当下降1个单位水头时,由于土体压缩和水的膨胀所释放出来的贮存水量。式(1)的定解边界条件如下:水头条件为h|Γ1=h(x,y,z,t);流量边界条件为(kijh,jni)|Γ2=q(x,y,z,t);自由面边界条件为{h|Γ3=z(kijh,jni)|Γ3=q=-μ∂h∂tcosθ；初始条件为h|t=0=h(x,y,z,0)。式中:h(x,y,z,t)为已知水头函数;q为自由面上因自由变动而引起的流量补给;μ为饱和差(自由面上升时)或给水度(自由面下降时),它表示自由面改变单位高度下,含水层单位截面积上吸收(自由面下降)或排出(自由面上升时)的水量,是无量纲数;θ为自由面外法线方向与垂线的交角;t为时间。2.2整体有限元模型对2.1节提出的数学模型采用Galerkin法和等参8结点六面体单元在空间域上进行合理离散,将单元支配方程进行集成,可得整体有限元支配方程:Κp+(S+G)∂p∂t=F。(2)在此基础上,对时间域进行合理离散,可得到有限元时间隐式差分迭代格式:[Κ+1Δt(S+G)]pt+Δt=F+1Δt(S+G)pt。(3)式中:K为渗透系数矩阵;S为总质量矩阵;G为自由面贡献矩阵;p为结点压力列阵;F为等效结点流量列阵;Δt为时间步长。2.3筛选适用的迭代格式由于自由面位置是未知的,式(3)无法一次求解,这里采用文献介绍的增量迭代法,可推导得如下适于计算的迭代格式:AΔpk+1t+Δt=ΔBk+1t+Δt。(4)式中:A为非线性方程的系数矩阵,A=Κ+1Δt(S+G)；ΔBk+1t+Δt=ΔFkt+Δt-1Δt(S+G)(pkt+Δt-pt)。2.4预处理共轭梯度算法采用预处理共轭梯度算法(preconditionedconjugategradientsmethod,PCG)对线性方程组式(4)进行联立求解。PCG法假设矩阵A被分为2个矩阵,即:A=Μ+Ν。(5)式中:M为预优矩阵;N为剩余矩阵。M为A的预条件形式,稀疏性与A相同,总是对称且正定的。有关预处理共轭梯度算法的具体求解方法见文献。2.5流子结构总传导矩阵按照上述介绍的辐射井的渗流行为,结合有限元法求解原理,将子结构作为主网格整体结构的一部分一起来考虑,则辐射井的非稳定渗流子结构总传导矩阵及相应流量列阵写成分块形式,如式(6)所示:Κ0iΔΡ0i=ΔF0i,i=1,2,⋯,n。(6)式中:Κ0i=(Κ+S+GΔtk)0i;ΔF0i=(ΔF+(S+G)(Δpk)0Δtk)0i;i为子结构编号,n为子结构总数;k为迭代时步;式中各量的计算域和边界被限制在母单元的模型中。对于子结构内部边界的处理与外部边界相同即可,详见参考文献。3工程应用3.1台兰河流域地下水台兰河地下水库工程位于新疆维吾尔自治区阿克苏地区温宿县境内的台兰河灌区,距温宿县城33km,距阿克苏市46km,距乌鲁木齐市955km,研究区位于314国道穿过2km范围内。台兰河地下水库工程主要由渗库工程、取水工程和输水工程三大部分组成。渗库工程是通过在上游修建简易人工分洪工程(蓄洪库、挡水堤坝或坑、塘等),加大地表水地下水的回灌量,形成可灌的地下水库;另外在下游选择合适的地形地质条件修建地下暗坝或天然暗坝,使上游地下松散地层内一定空间的地下水位壅高,形成地下水库。取水工程是在地下水库取水点修建大口径取水井截取第四纪松散地层中的地下水,采用低投入、高产出的方法提取(或自流)地下水,以起到尽量减少蒸发、以丰补歉、下游防洪等调节水库的作用。台兰河地下水库工程是将台兰河洪水期的洪流和冬闲期河水转换成地下水,再使地下水转换为地表水资源,以达到增加台兰河流域灌溉水量的目的;利用地表水和地下水联合调度台兰河流域需水过程,创新和优化配置流域内的水利工程布局。目前取水井工程的一期已完成4眼大口径辐射井及输水管道的建设,可进行自流式供水,流量达到0.3～0.5m3/s。台兰河地下水库取水工程的二期完成“自流虹吸”输水系统的建设。在一期建设的基础上完成输水廊道、集水廊道以及在输水廊道距F1号井800m的下游布置F0号辐射井的施工建设。其中输水廊道和集水廊道的最大降深均为10m,F0号辐射井和F1号辐射井的虹吸水位降深达到12m,设计流量1.0～1.2m3/s;并且利用控制闸调控F0号井水位,使其在高水位状态下将地下水引入高位水池,进入自压式滴灌系统,示范区面积1万ha。整个取水工程的平面布置见图3。3.2冲积砾质倾斜沉积单次进入洪积平原图3台兰河由北向南穿过洼地进入古木别致低山丘陵地带,海拔高程1498～1532m;又经老龙口豁口下泄进入洪积平原,在下游依次形成冲积砾质倾斜平原,地形坡降10‰～16‰,东西向宽度20km;再向南为洪积细土平原,地面高程1185～1220m,地面坡降5‰～8‰。台兰河地下水库示范区选择在314国道以北山前冲洪积扇储水构造的小阿吉沟谷中,地下水位埋藏深度1.5～3.0m,含水层为砂卵砾石层。3.3辐射井的巷道布置取水工程系统为“横坎儿井”式的大口径竖井和水平辐射井的集水系统,垂直地下水流向布设3眼大口径辐射井。由于地形坡度较大,在大口径井下游约0.8km处布置1眼大口径辐射集水井,采用“自流虹吸”原理,用钢管将上述4眼井联成一体,采用闸门进行控制。垂直地下水流向的大口径井井距为250m左右,井径3.5m,井深30m,水平辐射井每井3层,每层8～10个辐射井。辐射井采用裸孔钢管,管径159mm×10mm,单根长10～15m,最长达24m。每个辐射井最大降深为10m,单井出水量大于1200m3/h。在集水井下游修建引水管道进行自流式引水。引水管道直径1.2m,长1.45km,自压式进入农田。3.4井位平面位置分布台兰河地下水库目前取水工程共布置4眼大口径辐射井(F0-F3)和12口观测井(G01-G12),其平面位置分布见图4。抽水试验选定F2辐射井为抽水井(其余辐射井作为观测井),抽水时间从2010年6月14日14:30开始,一直持续到2010年7月26日4:00。其中经历了10次停电,具体停电时间见表1。3.5研究区初始渗流场计算计算模型范围本次重点研究的是F1号辐射井抽水过程和水位恢复的渗流场及示范工程附近地下水渗流场变化。从抽水试验中各观测井的观测资料可看出,F1号辐射井在连续抽水25d后,地下水影响范围约为500m。因此,为减少边界水位对地下水渗流场的影响,将计算模型平面范围取为以F1号辐射井为中心,半径1.5km。研究区含水层深度为300～500m,垂直向为高程900m的平面至地表面。边界条件及初始条件研究区位于台兰河洪积扇与细土平原交界带上的洪水冲沟内,地下水流沿洪水冲沟由北向南流,地面坡降5‰～8‰,地下水的水力坡降约为4.4‰。根据观测孔和辐射井抽水前的水位拟合得到了研究区初始稳定渗流场(图5),根据初始渗流场确定计算范围内北部和南部边界为定水头边界,水头值分别为1202.5m和1187.0m。由于抽水试验历时较短(约26d),在此期间没有大的降雨,只下过2d小雨,且该地区蒸发量较大,降雨对地下水位几乎没有影响,因此计算模型的地表面边界为隔水边界,不考虑降雨入渗的影响。为了验证初始渗流场的合理性,将边界条件代入模型,利用稳定渗流计算得到了研究区的初始渗流场,各钻孔的计算水位和观测水位拟合情况见图6。由图6可看出,稳定渗流场计算得到的各观测井的计算水位和观测水位较为接近,且均匀分布在45°线上,由此可说明模型上下游边界选择较为合理。时间步长选取整个抽水试验过程共进行了约40d,由于抽水试验各观测井均采用了电子自动记录水位(每30min记录一次),数据较大,为了能够反映整个抽水过程的水位变化,计算模型的时间步长选为1d,每个时间步长末端时刻为12:00;为充分反映停电过程对地下水位的影响,每次停电时刻和停电后水泵开启的时刻单独作为一个时间步长。计算时间从2010年6月14日14:30到2010年6月30日14:30结束,共40个时间步长。计算时间步长内辐射井F2水位随时间变化曲线见图7。研究区网格剖分及辐射井子结构网格剖分图8给出了辐射井子结构单元的三维网格图。按照上述计算模型范围,结合辐射井和观测井的位置建立了研究区的有限元分析模型,对辐射井F2进行子结构模型处理,其他观测井都布置在单元结点上,计算区共剖分单元12000个,结点10487个,剖分网格见图8a。根据辐射井F2的结构设计(图1),确定了辐射井的子结构母单元网格的平切图见图8b,剖分形成的辐射井的三维子结构单元见图8c。渗流场计算及计算结果分析图9给出各观测井计算水位和观测水位对比图。图10给出了第1次停电时的第2层辐射管平面的地下水头等值线图和剖面水头等值线图。图11给出了第10次停电后水泵启动时刻的第2层辐射管平面的地下水头等值线图和剖面水头等值线图。利用观测井计算结果对模型参数进行了校正分析,表2给出了采用改进的潜水斜率解析法、Neuman算法和辐射井子结构法的计算反演等方法计算得到研究区平均渗透系数和给水度值。从表2可看出,反演得到的渗透系数比前2种方法稍大,给水度较为接近。主要原因是前2种方法不适合于辐射井抽水试验渗透系数反演。另外,由于数值分析模型对地层进行了概化,没有反映地层的各向异性,反演值存在一定的误差。从图9结合计算结果可知:所有观测井计算水位的平均绝对误差为0.22m,单井水位的平均误差最大值为0.40m,最小值为0.02m,各观测井地下水位的模拟值与实测值的变化趋势吻合较好,说明模型拟合的地下水位变化趋势与研究区的水