地下水流数值法.docVIP

地下水流数值法 一 地下水流的有限单元法 图1是某平面二维有限单元三角形网格的一个局部，其中3, 7, 11, 21, 32为节点编号，而(2), (5), (7), (9)为三角形单元编号，它们的关系如表1。假设含水层为均质等厚承压含水层，请设计这些节点的合适坐标和含水层参数，计算单元(5)的单元导水矩阵 [De]-(5)，给出对应节点11和21的总体导水矩阵元素aLp=a11,21。说明计算方法。 表1 单元 节 点 i j k (2) 3 21 7 (5) 7 21 11 (7) 32 11 21 (9) 3 32 21 各节点坐标如表2所设，含水层渗透系数K=5 m/d，厚度M=20 m。 表2 节点 坐标 3 (10,10) 7 (15,27) 11 (35,30) 21 (23,20) 32 (33,12) 计算(5)和(7)几何系数如表3： 表3 单元 bi bj bk ci cj ck △ (5) -5 3 7 12 -10 8 47 (7) 10 8 -18 -12 10 2 98 则(5)号单元的流量为： 其单元导水矩阵 [De]-(5)为 (7)号单元的流量为： 节点11和21的总体导水矩阵元素 二 地下水流的有限差分法 图2为某平面二维有限差分矩形网格的一个局部，网格单元的尺寸用a, b表示，采用格点法建立有限差分方程，格点编号分别为1, 2, 3, ……, 7, 8, 9等，在单元格5中含有一个抽水强度为Q=200 m3/d的抽水井，该承压含水层的厚度为M=15 m，渗透系数为 K=6 m/d，贮水率为Ss=0.00001 m(1。取时间步长为(t，已知n时步各个格点的水头hn求n+1时步的水头hn+1, 则单元格5的差分方程可写为 请设计合适的a, b值并求出上述公式中的C4, C2, C6, C8, C5, D5, E5。说明计算方法。 已知Q=200 m3/d，M=15 m，K=6 m/d，Ss=0.00001 m(1，设a=10 m，b=20 m。 对5号格点建立流量均衡方程，有： 根据达西定律： 带入流量均衡方程化简得： 与题目公式对应，则 三 地下水流数值模拟算例 请设计一个非稳定流的地下水流动模型，运用合适的软件对它进行模拟，并写出一份模拟报告。要求如下： 模型至少有一条定水头边界、一条隔水边界、一个抽水井和一个观测孔，模拟层不少于3层。其它条件自定。 写出模型的数学描述，包括水流偏微分方程、边界条件和初始条件。 说明使用的软件，模型的网格划分方法（包括总节点数、格点数或单元数），数值法的类型（如：有限元法、有限差分格点中心法、有限差分节点法），时间步长的控制，参数的设置等。 根据模拟结果说明抽水井附近水头分布的变化和观测孔水位的变化趋势，分析模拟结果与含水层渗透性参数和贮水性参数之间的关系。 阐述其它感兴趣的问题。 现对一个400×400 m的均质正方形区域利用Processing Modflow软件进行非稳定流的地下水流动模拟。 1、模型条件 区域面积400×400 m，厚10 m；东西两侧为定水头边界，东侧水位8 m，西侧水位10 m，南北两侧为隔水边界。 区域中心有一口抽水井，抽水量为10 m3/s，附近有一个观测孔。 2、数学描述 其中，, D为研究区域，ε为源汇项，Ss为储水率，μ为给水度。 3、模型建立 使用Processing Modflow软件进行有限差分法模拟，将区域剖分成20×20的网格，每个网格单元表示20×20 m。模拟层为3层，由上到下顶板标高分别为10 m、6 m、3 m。应力期为365 d，步长为12；水平和垂直渗透系数分别为10 m/s和1 m/s；储水率为0.0001 m-1，给水度为0.2。 抽水井位于（10，11）单元中心，观测孔位于（200，200）坐标处。 4、模拟结果 先按照模型条件进行稳定流模拟，将得到的水头值再作为初始水头条件继续进行非稳定流的模拟，抽水量增加为第1层100 m3/s，第2、3层50 m3/s（图2）。 a b c 图2 模拟水头值 a-第1层 b-第2层 c-第3层 抽水井附近水头减小明显，且其水力坡度西侧比东侧大。观测孔水位下降可以反应抽水井的水头变化趋势（图3），抽水开始，水头下降迅速，随时间步长增加，水位下降逐渐变缓，趋于稳定值。 图3 观测孔水头随时间变化曲线 当含水层渗透系数较大时，水平渗透系数50 m/s，垂直渗透系数5 m/s，模拟结果如图4和图5。抽水开始时水位下降迅速，但水位下降小且很快就趋于稳定。含水层渗透性参数对此模型有明显影响。 图4 渗透系数增大时第1层地下水水位 图5 渗透系数增大时观测孔水位 当储