AAM  >> Vol. 6 No. 9 (December 2017)

    基于有限差分法的承压稳定井流计算方法的比较
    Comparison of Calculation Methods for Confined Steady Well Flow Based on Finite Difference Method

  • 全文下载: PDF(1214KB) HTML   XML   PP.1201-1206   DOI: 10.12677/AAM.2017.69145  
  • 下载量: 637  浏览量: 1,466  

作者:  

张思宇,周德亮,孙婧祎:辽宁师范大学数学学院,辽宁 大连

关键词:
有限差分法井水位计算比较奇点磨光Finite Difference Method Well Level Calculation Comparison Polishing Singular Point Method

摘要:

分别用有限差分法的早期方法、井水位优化方法和奇点磨光法计算实例中井水位。由实例计算并比较显示,奇点磨光法计算得到的井水位与井处真实水位最接近。

Based on finite difference method, we use early method, well level optimization method and pol-ishing singular point method to calculate the well level in example respectively. The example is calculated and compared. The well level calculated by polishing singular point method is closest to the real well level.

1. 引言

在水文地质实践中,很多情况下都涉及到井水位计算问题。它不仅关系到地下水资源的合理开发与利用,而且还关系到水源地建设的投资。对于井水位的计算,使用最广泛的方法为有限差分法 [1] [2] ,由地下水动力学理论和实际经验可知,井流是径向流动的,所以按有限差分法处理井流时,模拟得到的井点水位并不能代表井中的真实水位,这为动态分析地下水井水位运动的变化规律及水资源的开发与利用带来了负面影响。为解决这一问题,1967年Prickett提出了针对正方形网格井点处水位的校正公式。1978年Peacemen基于Theim公式提出了Well Index校正方法。1984年李贻祖,陈平等基于有限元的模拟,通过井周的特殊划分求导出井中水位修正公式。1990年陈崇希,唐仲华等根据不稳定流理论提出了正方形差分网格的井中水位校正公式。2003年曹洪,尹小玲等在李贻祖方法的基础上,将单元按等分圆周划分得出井中水位修正公式。本文针对有限差分法的早期方法 [3] 、井水位优化方法 [4] 和奇点磨光法 [5] 这三种井水位计算方法进行描述,并将它们在一具体实例中的计算结果与井中真实水位进行比较。

2. 井水位早期计算方法

考虑如下承压稳定井流问题

{ T Δ u = f ( x , y ) , x , y G ( 2.1 ) u | Γ 1 = a ( x , y ) , Γ 1 ( 2.2 ) T u n | Γ 2 = b ( x , y ) , Γ 2 ( 2.3 ) T u n | C k = q k 2 π ρ k , C k , k = 1 , 2 , , n ( 2.4 )

其中, Δ u = 2 u 2 x + 2 u 2 y , T 为导水系数, u 为水头函数, f α β 是已知函数, q k 为第 k 口井的开采量, ρ k 为第 k 口井的半径。

对于有限差分法,先进行网格剖分,取定沿 x 轴, y 轴的步长 h x h y ,做两族与坐标轴平行的直线: x = i h x i = 0 , ± 1 , y = j h y j = 0 , ± 1 , ,两族直线的交点 ( i h x , j h y ) 称为节点,记为 ( i , j ) 。对于区域内的节点,分别用沿 x y 方向的二阶中心差商代替 u x x u y y ,即

Δ u = ( u i + 1 , j 2 u i , j + u i 1 , j h x 2 + u i , j + 1 2 u i , j + u i , j 1 h y 2 ) = 1 h x 2 u i + 1 , j 1 h x 2 u i 1 , j 1 h y 2 u i , j + 1 1 h y 2 u i , j 1 + 2 ( 1 h x 2 + 1 h y 2 ) u i , j

同理,对于在二类边界 Γ 2 C k 上的节点,也可以用二阶中心差商代替节点在其外法线方向导数,即将(2.1)~(2.4)转化为求解如下问题

{ T [ 1 h x 2 u i + 1 , j 1 h x 2 u i 1 , j 1 h y 2 u i , j + 1 1 h y 2 u i , j 1 + 2 ( 1 h x 2 + 1 h y 2 ) u i , j ] = f , G ( 2.5 ) u | Γ 1 = a ( x , y ) , Γ 1 ( 2.6 ) T ( u i , j + 1 2 u i , j + u i , j 1 h y 2 ) | Γ 2 = B ( x , y ) , Γ 2 C k ( 2.7 )

联立(2.5)、(2.6)和(2.7)解出 u i , j ( i = 0 , ± 1 , , j = 0 , ± 1 , ) ,即得到承压稳定井流问题的解。

井水位早期计算方法是对(2.1)~(2.4)进行离散,将原问题转化为如下问题

{ x ( T u x ) + y ( T u y ) = k = 1 n q k δ ( x x ¯ k , y y ¯ k ) + ε ( x , y ) , G ( 2.8 ) u | Γ 1 = a ( x , y ) , Γ 1 ( 2.9 ) T u n | Γ 2 = b ( x , y ) , Γ 2 ( 2.10 )

再用 q k S k 近似井点处平均流量,即(2.8)式中右端项量,其中 S k 是井点附近网格对偶剖分形成的面积。

3. 井水位优化方法

图1,依据有限差分法在井点四周布置间距为 Δ x Δ y 的矩形网格,假定抽水井所在节点 i 与周围节点间已经形成了拟稳定流动(即水位随时间变化,而水力坡度不随时间变化的不稳定流动)。由Theis公式的近似式可知, j 点水位 h j m 点水位 u m 与井中的水位 u w 之间满足下述关系:

u j = u w + Q 2 π T ln Δ x r w (3.1)

u m = u w + Q 2 π T ln Δ y r w (3.2)

其中, Q 为抽水井的流量; T 为导水系数; r w 为井的半径。

同理,得到节点 k n 处水位 u k u n

u k = u w + Q 2 π T ln Δ x r w (3.3)

u n = u w + Q 2 π T ln Δ y r w (3.4)

将有限差分法模拟得到的节点 i 处的水位记为 u i ,利用差分公式计算的从四周进入 i 点处的流量与水井的流量是相同,所以有

Q = T Δ x [ u m u i Δ y + u n u i Δ y ] + T Δ y [ u j u i Δ x + u k u i Δ x ] (3.5)

把上述 u j u m u k u n 的表达式都代入到(3.5)式,可得到如下优化公式

Figure 1. Plane difference grid near well point

图1. 井点附近平面差分网格

u w = u i Q 2 π T Δ x Δ y Δ x 2 + Δ y 2 [ Δ y Δ x ln Δ x r w + Δ x Δ y ln Δ y r w π ] (3.6)

当取 Δ x = Δ y 时,(3.6)式可简化为

u w = u i Q 2 π T ( ln Δ x r w π 2 ) (3.7)

一般情况下, Δ x 会远大于 r w ,所以 u w < u i ,这就是有限差分法得到的井中水位偏高的原因。

4. 奇点磨光法

奇点磨光法的基本思想是将问题(2.1)~(2.4)化成求解一个无井存在且方程右端充分光滑的问题。为此,在区域 G 内定义辅助函数 v = k = 1 n v k ,即相对于每一个井造一个定义在区域 G 上的辅助函数 v k ( k = 1 , 2 , , n ) ,为使辅助函数在抽水井处的奇性与水头函数 u 在抽水井处的奇性相同,函数 v k 应满足的条件为:(1) T v k n | c k = Q 2 π ρ k (2); v k | r k = R k = 0 ;(3) v k r k | r k = R k = 0 ;(4) v k | r k > R k = 0 r k = ( x x ¯ k ) 2 + ( y y ¯ k ) 2 R k 是以井点 ( x ¯ k , y ¯ k ) 为心的圆周半径,要求该圆周不能越过 Ω 的边界 Γ = Γ 1 + Γ 2

w = u v = u k = 1 n v k ,对于问题(2.1)~(2.4)就转变为求解下述问题

{ x ( T w x ) + y ( T w y ) = f ( x ( T v x ) + y ( T v y ) ) , Ω ( 4.1 ) w ( x , y ) | Γ 1 = a ( x , y ) , Γ 1 ( 4.2 ) T w n | Γ 2 = b ( x , y ) , Γ 2 ( 4.3 )

基于有限差分法求得 w ,即得到问题(2.1)~(2.4)的解 u = w + v

f 光滑无奇性。为保证求解精度,函数 v k 的选择除了要满足前述4个条件外,还要能使方程(4.1)的右端项 x ( T v x ) + y ( T v y ) 无奇性且连续。

5. 数值算例

为比较上述三种方法的求解效果,具体考虑如下混合边界井流计算问题。存在一个长为 a = 2000 m ,宽为 b = 1000 m 的矩形区域,东西边界 A D ¯ B C ¯ 为定水头边界,南北边界 A B ¯ D C ¯ 为隔水边界。承压含水层的导水系数 T = 100 m 2 / d h ( x , y ) 为水头函数。在区域内有一抽水井以 Q = 2000 m 2 / d 的流量抽水,抽水井的中心坐标 ( x 0 ¯ , y 0 ¯ ) = ( 1100 , 500 ) ,井的半径为 ρ = 0.9 m ,井的边界设为 L

针对奇点磨光法,根据辅助函数 v 满足的条件,取 v 的紧支域半径为 R = 400 m ,利用待定系数法可得到

v = ln r ( r 1 ) + 1 2 ( r 1 ) 2 + 2 3 ( r 1 ) 3 ( r < R ) v = 0 ( r R )

其中 r = ( x x ¯ 0 ) 2 + ( y y ¯ 0 ) 2 。三种计算方法得到的井水位的近似解与井中真实水位的绝对误差分别如图2图3图4所示。

Figure 2. Early method approximate water level and true water level error

图2. 早期方法近似水位与真实水位误差

Figure 3. Water level optimization method approximate water level and true water level error

图3. 井水位优化方法近似水位与真实水位误差

Figure 4. Approximate water level by singular point smoothing method and true water level error

图4. 奇点磨光法近似水位与真实水位误差

6. 结语

由实例计算得出,早期方法在井点处的计算误差较大;井水位优化方法可以达到很好的模拟效果但只适合用于矩形网格,对实际应用具有一定的局限性;奇点磨光法的求解精度高,同时很好的处理了井点处奇异性问题,对于实际应用没有网格划分上的限制。因此,对于井水位计算奇点磨光法是一种值得推广的方法。

文章引用:
张思宇, 周德亮, 孙婧祎. 基于有限差分法的承压稳定井流计算方法的比较[J]. 应用数学进展, 2017, 6(9): 1201-1206. https://doi.org/10.12677/AAM.2017.69145

参考文献

[1] 张宏仁. 解地下水渗流问题的有限差分方法(二) [J]. 水文地质工程地质, 1980(2): 58-61.
[2] 陈崇希, 等. 地下水流数值模拟理论方法及模型设计[M]. 北京: 地质出版社, 2014.
[3] 孙讷正. 地下水流的数学模型和数值方法[M]. 北京: 地质出版社, 1981: 110.
[4] 陈崇希, 王旭升, 胡立堂. 地下水数值模拟中抽水井水位校正[J]. 水利学报, 2007, 38(4): 481-485.
[5] 吉林大学数学系计算服务站地下水小组. 地下水非稳定流计算的有限元方法(I) [J]. 吉林大学学报, 1977(1): 15-27.