張 敏，余震果

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

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

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

如圖1［1］，一承壓含水層完整井，含水層為均質(zhì)等向，厚度M，滲透系數(shù)K，含水層中的原始水位H0，抽水穩(wěn)定后，測得井中水位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ū)域。對應(yīng)此題c=KM，a=0，f=0。Ω 為環(huán)形，邊界為 Dirichlet條件(hu=r)，在邊界x2+y2=R2上時，h=1，r=170，在邊界 x2+y2=上時，h=1，r=120。

應(yīng)用PDE工具箱來解決這個問題的具體步驟:

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

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

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

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

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

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

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

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

圖2

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

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

圖3

此方程為非線性橢圓型方程，此類方程的基本形式為-▽(c(u)▽u)+a(u)u=f(u)，inΩ 。其中 c，a，f可以是解 u的函數(shù)。對應(yīng)此題c=Kh，a=0，f=0。Ω為矩形區(qū)域，邊界條件為Dirichlet條件(hu=r)，h和r是定義在αΩ上的函數(shù)，對于非線性情形，h和r可以依賴于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。

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

3 結(jié)語

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

圖4

圖5

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

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

［3］陸君安，尚濤，謝進，等.偏微分方程的 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)機.2010，37(3):89—91.