張 敏，余震果

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

有限元法是解決地下水問(wèn)題的主流方法，利用此方法求解問(wèn)題需選定單元的形狀;對(duì)求解域作剖分;對(duì)節(jié)點(diǎn)編號(hào);構(gòu)造基函數(shù);形成有限元方程;求解方程。有些實(shí)際問(wèn)題求解區(qū)域比較復(fù)雜，在對(duì)幾何區(qū)域進(jìn)行單元網(wǎng)格剖分、節(jié)點(diǎn)編號(hào)時(shí)比較繁瑣，直接編寫(xiě)程序具有一定的難度。偏微分方程(PDE)工具箱提供了研究和求解空間二維偏微分方程問(wèn)題的一個(gè)強(qiáng)大而又簡(jiǎn)單便捷、靈活實(shí)用的環(huán)境。利用PDE工具箱求解偏微分方程問(wèn)題有兩種方法，一種是直接使用圖形用戶(hù)界面(GUⅠ)，可以從繁雜的編程中解脫出來(lái)，很容易實(shí)現(xiàn)復(fù)雜的幾何區(qū)域剖分、直接加密網(wǎng)格，也可以進(jìn)行復(fù)雜的定解條件下的偏微分方程求解，直接生產(chǎn)M代碼。另一種是采用命令函數(shù)編寫(xiě)程序計(jì)算。應(yīng)用GUⅠ解決地下水問(wèn)題步驟為:1.區(qū)域選擇;2.模型建立;3.定義邊界;4定義PDE類(lèi)型和系數(shù);5.三角形網(wǎng)格剖分;6.求解PDE;7.解的圖形表達(dá);8.數(shù)據(jù)輸出。本文將利用MATLAB PDE工具箱的基本原理及步驟解決承壓和非承壓的穩(wěn)定流問(wèn)題，只具體應(yīng)用了GUⅠ沒(méi)有應(yīng)用命令行編輯去解決問(wèn)題。

1 承壓含水層定降深完整井穩(wěn)定流計(jì)算

如圖1［1］，一承壓含水層完整井，含水層為均質(zhì)等向，厚度M，滲透系數(shù)K，含水層中的原始水位H0，抽水穩(wěn)定后，測(cè)得井中水位hw，濾水管半徑rw，抽水井的影響半徑R。取 K=0.8 m/d;M=100 m;R=300 m;H0=170 m;rw=80 m;hw=120 m。

圖1

此方程屬于橢圓型方程，它的基本形式為:-▽(c▽u)+au=f，inΩ。其中Ω是平面有界區(qū)域。對(duì)應(yīng)此題c=KM，a=0，f=0。Ω 為環(huán)形，邊界為 Dirichlet條件(hu=r)，在邊界x2+y2=R2上時(shí)，h=1，r=170，在邊界 x2+y2=上時(shí)，h=1，r=120。

應(yīng)用PDE工具箱來(lái)解決這個(gè)問(wèn)題的具體步驟:

1)在MATLAB命令窗口輸入 pdetool并運(yùn)行，出現(xiàn) PDE圖形用戶(hù)界面(GUⅠ)。

2)在GUⅠ中對(duì)區(qū)域進(jìn)行選擇，此題區(qū)域?yàn)榄h(huán)形，分步完成幾何區(qū)域C1、C2。用菜單或快捷工具畫(huà)出矩形 C1、C2，然后再Set formula欄直接鍵入C1-C2。

3)定義邊界條件，可單擊Boundary菜單中Specify Boundary Conditions選項(xiàng)或直接單擊Ω，打開(kāi)Boundary Conditions對(duì)話框，輸入邊界條件，此題全部設(shè)為Dirichlet條件，邊界顏色都為紅色。

4)定義PDE類(lèi)型和PDE系數(shù)，選擇PDE菜單中的PDE Mode命令，進(jìn)入偏微分方程模式，再單擊PDE Specification選項(xiàng)，設(shè)置方程為橢圓型，在type of PDE中選擇Elliptic，再分別輸入?yún)?shù) c，a，f，c=160，a=0，f=0。

5)對(duì)區(qū)域進(jìn)行三角形網(wǎng)格剖分，可用initmesh和refinemesh進(jìn)行剖分，還可以選擇Mesh菜單中的Ⅰnitialize Mesh命令進(jìn)行剖分、Refine Mesh命令對(duì)網(wǎng)格加密。

6)選擇Solve菜單中Solve PDE命令，或直接單擊“=”求解偏微分方程且顯示圖形解。

7)解的圖形表達(dá)，單擊Plot菜單中的Parameters選項(xiàng)，打開(kāi)Plot Selection對(duì)話框，選中相應(yīng)的選項(xiàng)，然后單擊 Plot，顯示三維圖形解。

此類(lèi)地下水問(wèn)題，對(duì)網(wǎng)格進(jìn)行剖分，直接編寫(xiě)原始程序都較復(fù)雜。利用PDE工具箱對(duì)區(qū)域進(jìn)行網(wǎng)格剖分直觀簡(jiǎn)捷。本題對(duì)區(qū)域進(jìn)行6次初始化網(wǎng)格剖分，使解的精度得以提高。圖2為承壓含水層定降深完整井的降水漏斗曲面圖。

圖2

2 非承壓含水層穩(wěn)定流計(jì)算

如圖3［1］所示是一平面矩形河間地區(qū)，河口沖積扇平面中設(shè)置的灌，排河渠網(wǎng)的一典型地塊，由于主河渠水位均較高，所以造成河渠間的地下水水位過(guò)高，這一高位地下水長(zhǎng)期侵漬，使土地鹽化，堿化。為了排鹽排堿，進(jìn)行耕種，在垂直于河渠方向沒(méi)隔間距挖一條深排水溝進(jìn)行排油，以降低地下水位。

圖3

此方程為非線性橢圓型方程，此類(lèi)方程的基本形式為-▽(c(u)▽u)+a(u)u=f(u)，inΩ 。其中 c，a，f可以是解 u的函數(shù)。對(duì)應(yīng)此題c=Kh，a=0，f=0。Ω為矩形區(qū)域，邊界條件為Dirichlet條件(hu=r)，h和r是定義在αΩ上的函數(shù)，對(duì)于非線性情形，h和r可以依賴(lài)于u。滲透系數(shù)為K;已知水頭為h0;主渠內(nèi)水位為h1;主河內(nèi)水位為h2。分別取K=-0.8 m/d;h0=10 m;h1=30 m;h2=40 m;l1=250 m;l2=120 m。

求解步驟與上例有一點(diǎn)區(qū)別，即對(duì)上例所述第六步進(jìn)行設(shè)置，選擇 Solve菜單中 Parameters的 Use nonlinear solver，如圖4。其余大致相同。圖5是利用PDE工具箱對(duì)上面二維穩(wěn)定流問(wèn)題做出的水頭函數(shù)圖形。

3 結(jié)語(yǔ)

本文是直接利用GUⅠ解決地下水計(jì)算問(wèn)題，應(yīng)用此種方法簡(jiǎn)單便捷。除了用GUⅠ解決還可以用工具箱中的命令來(lái)創(chuàng)建描述幾何圖形的M文件來(lái)解決問(wèn)題。有時(shí)會(huì)遇到利用GUⅠ無(wú)法解決的問(wèn)題，例如幾何圖形不是由直線、圓弧、橢圓弧及其組合而成的圖形時(shí)，就只能用工具箱中的命令函數(shù)編寫(xiě)程序計(jì)算，此法比直接編寫(xiě)原始程序簡(jiǎn)單快速。

圖4

圖5

［1］王君連.工程地下水計(jì)算［M］.北京:中國(guó)水利水電出版社.2004.

［2］薛禹群，葉淑君，謝春紅，張?jiān)?多尺度有限元在地下水模擬中的應(yīng)用［J］.水力學(xué)報(bào).2004(7).

［3］陸君安，尚濤，謝進(jìn)，等.偏微分方程的 MATLAB解法［M］.武漢:武漢大學(xué)出版社.2001.

［4］陸金甫，關(guān)治.偏微分方程的數(shù)值解法［M］.北京:清華大學(xué)出版社.2004.

［5］蘇金明，張蓮花，劉波.MATLAB工具箱應(yīng)用［M］.北京:電子工業(yè)出版社.2004.

［6］陳群，謝家燁.PDETOOL工具箱在磁共振技術(shù)中的應(yīng)用［J］.江西科學(xué).2009，27(2).

［7］李明.偏微分方程的 MATLAB 解法［J］.湖南農(nóng)機(jī).2010，37(3):89—91.