承壓穩(wěn)定井流計(jì)算的單位分解配點(diǎn)法

郭冬冬，周德亮

(遼寧師范大學(xué)數(shù)學(xué)學(xué)院，遼寧 大連 116029)

?

承壓穩(wěn)定井流計(jì)算的單位分解配點(diǎn)法

郭冬冬，周德亮

(遼寧師范大學(xué)數(shù)學(xué)學(xué)院，遼寧 大連 116029)

[摘要]用單位分解配點(diǎn)法求解承壓含水層中地下水向井的穩(wěn)定流動(dòng)問題，該方法擺脫了背景網(wǎng)格的束縛，是一種真正的無網(wǎng)格方法，根據(jù)具體模型計(jì)算后發(fā)現(xiàn)其不僅實(shí)施簡(jiǎn)單，而且計(jì)算精度高。

[關(guān)鍵詞]單位分解；配點(diǎn)法；穩(wěn)定井流；無網(wǎng)格

無網(wǎng)格(Meshfree)方法是近幾十年來迅速興起的一種數(shù)值計(jì)算方法，它很好的克服了有限差分和有限元法等方法需要網(wǎng)格的缺陷，并利用一組散布在求解域中及其域邊界上的節(jié)點(diǎn)來表示(而非離散)該求解域和其邊界，具有收斂快，后處理方便等優(yōu)點(diǎn)。其中單位分解配點(diǎn)法[1]繼承了無網(wǎng)格方法的性質(zhì)，一方面可以在子域內(nèi)調(diào)節(jié)節(jié)點(diǎn)的密度，也可以改變某個(gè)子域的形狀；另一方面局部逼近空間可以包含非多項(xiàng)式函數(shù),從而很好地局部逼近未知解。因此單位分解配點(diǎn)法具有靈活性且近似解有較高的精度。本文將單位分解配點(diǎn)法用于求解地下水的承壓穩(wěn)定井流動(dòng)的問題。

1配點(diǎn)法

考慮如下定解問題

L{u(x)}=f(x)，x∈D

(1)

B{u(x)}=g(x)，x∈?D

(2)

其中設(shè)D是有界區(qū)域，?D是D的邊界，L是偏微分算子，B是邊界微分算子，f(x)，g(x)是給定的函數(shù)。

(3)

一般情況下中心點(diǎn)取配置節(jié)點(diǎn)，αj是待定系數(shù)，N是配點(diǎn)總數(shù)，φ(‖x-ξj‖2)表示徑向基函數(shù)[2]，常見的徑向基函數(shù)有

高斯徑向基函數(shù)：φ(r)=e-(εr)2(ε>0)；

緊支撐徑向基函數(shù)：φ(r)=(1-εr)4+(4εr+1) (ε>0)；

將(3)代入(1)(2)中，使其分別在內(nèi)節(jié)點(diǎn)，一類邊界節(jié)點(diǎn)和二類邊界節(jié)點(diǎn)上成立，這樣線性方程組便可寫成

Aα=b

(4)

其中α=(α1，α2，…，αN)T是待定的未知向量，

b=(f(x1)，…，f(xN))，g(xNi+1)，g(xNi+N1+1)，g(xN)T

s(x)=φ(x)α=φ(x)A-1u=Ψ(x)u，其中Aα=u，u=[u(x1)…u(xN)]T，

φ(x)=[φ(‖x-x1‖)，…，φ(‖x-xN‖)]，Ψ(x)=φA-1

這里的Ψj(x)具有如下性質(zhì)

2單位分解配點(diǎn)法

?x∈DI(x)={j|x∈Dj}，I(x)≤K

集合I(x)表示的是包含x的子域的編號(hào)，包含x的子域個(gè)數(shù)不超過K個(gè)，這里常數(shù)K與子域的數(shù)量M無關(guān)。在單位分解配點(diǎn)法[1]中,定解問題的解函數(shù)u(x)的近似函數(shù)s(x)可以寫成下面形式：

(2.1)

其中sj(x)是u(x)在子域Dj內(nèi)的近似函數(shù)，wj是定義在Dj上的緊支正定權(quán)函數(shù),它可以通過Shepard的方法構(gòu)造,即表示為

(2.2)

其中

(2.3)

是子域Dj上的緊支徑向基函數(shù).在(2.2)中我們選擇Wendland緊支徑向基函數(shù)

(2.4)

如果對(duì)每個(gè)節(jié)點(diǎn)Xi∈Dj，等式(2.4)中的函數(shù)sj(x)，j=1，…，M，都具有局部插值性質(zhì)sj(xi)=u(xi)，那么全局單位分解近似值繼承了局部插值的性質(zhì)，s(xi)也就可以寫成下面的形式：

利用單位分解配點(diǎn)法解決偏微分方程問題時(shí)，需要計(jì)算在空間微分算子在內(nèi)部節(jié)點(diǎn)上的值和邊界微分算子在邊界節(jié)點(diǎn)上的值，為此需要計(jì)算α階的導(dǎo)數(shù)，應(yīng)用Leibniz定理，可以得到近似解的α階導(dǎo)數(shù)為

其中α和β是多重階數(shù)。

3對(duì)井的處理

考慮平面承壓穩(wěn)定井流混合問題[4]

其中,T為導(dǎo)水系數(shù),h(x，y)為水頭函數(shù),ε(x，y)為越流補(bǔ)給強(qiáng)度，q為補(bǔ)給量，qk為第k口井的開采量，ρk為第k口井的半徑。

通過實(shí)際計(jì)算可發(fā)現(xiàn)，和傳統(tǒng)數(shù)值方法一樣，如果直接從問題中進(jìn)行離散，用徑向基函數(shù)配點(diǎn)法求解將會(huì)在井口附近產(chǎn)生很大的誤差，使結(jié)果不可用。原因是井口相對(duì)于求解域非常接近為一個(gè)點(diǎn)，從數(shù)學(xué)角度看就是一個(gè)奇點(diǎn)，目前解決這一問題的常規(guī)方法是將問題轉(zhuǎn)變?yōu)槿缦碌牡葍r(jià)問題

由此求得w,立即得到定解問題的解h=w+v。

圖1　計(jì)算域及節(jié)點(diǎn)配置圖

圖2　近似解

圖3　解析解

設(shè)ε(x，y)光滑無奇性。要保證數(shù)值求解的精度，函數(shù)vk除了要滿足文[6]給出的前述4個(gè)條件

4數(shù)值計(jì)算

在計(jì)算區(qū)域上按圖1所示布置節(jié)點(diǎn)。單位分解配點(diǎn)法得到的水頭函數(shù)近似解、解析解以及絕對(duì)誤差如圖2、圖3和圖4所示，本文中選用的是高斯徑向基函數(shù)，權(quán)函數(shù)選用緊支徑向基函數(shù)，參數(shù)取160，與解析解比較可知使用單位分解配點(diǎn)法配求解地下水向井流動(dòng)的問題時(shí)具有靈活性，且近似解有較高的精度，由于該方法是在局部區(qū)域上應(yīng)用配點(diǎn)法，因而形成的矩陣是稀疏的，從而減少了計(jì)算的費(fèi)用。

圖4　絕對(duì)誤差

[1]Safdari-Vaighani A,Heryudono A,Larsson E.A Radial Basis Function Partition of Unity Collocation Method for Convection-Diffusion Equations Arising in Financial Applications[J]. Sci.comput,2014.10.4-7.

[2]Sharan M, Kansa E J, Gupta S. Application of multiquadric method fornumerical solution of elliptic partial differential equations[J]. Applied Mathematics Computation, 1997, 84: 275-303.

[3]Franke C, Schaback R. Solving partial differential equations by collocation using radial basis functions[J]. Applied Mathematics Computation, 1998, 93: 73-82.

[4]孫訥正.地下水流的數(shù)學(xué)模型和數(shù)值方法[M]. 北京:地質(zhì)出版社.1981.7-48,110.

[5]楊天行.非穩(wěn)定流的有限元法中承壓水析奇與潛水逐步析奇的一個(gè)方法[J].長(zhǎng)春地質(zhì)學(xué)院學(xué)報(bào).1978(4): 37-54.

[6]周培德.計(jì)算幾何——算法分析與設(shè)計(jì)[M]. 北京;清華大學(xué)出版社.1995.88-92.

[收稿日期]2015-11-16

[作者簡(jiǎn)介]郭冬冬(1991-)，女，遼寧大連人，在讀碩士研究生，主攻方向：偏微分?jǐn)?shù)值解法。

[中圖分類號(hào)]P641.139

[文獻(xiàn)標(biāo)識(shí)碼]B

[文章編號(hào)]1004-1184(2016)01-0051-02