周旭 江山 羅陽 居加穎

【摘要】研究一維固定長度的水槽中化學(xué)污染物的穩(wěn)態(tài)分布情況，利用有限元法作數(shù)值模擬，將其有限元解與真解進(jìn)行比較，發(fā)現(xiàn)隨著網(wǎng)格剖分?jǐn)?shù)的不斷增大，絕對(duì)誤差的數(shù)量級(jí)可達(dá)10-6，L2范數(shù)的數(shù)量級(jí)可達(dá)10-5，且有高于1階收斂結(jié)果.

【關(guān)鍵詞】有限元解；誤差分析；水槽化學(xué)污染

【基金項(xiàng)目】江蘇省高等學(xué)校大學(xué)生實(shí)踐創(chuàng)新訓(xùn)練計(jì)劃項(xiàng)目

有限元方法［1-3］是求解各種復(fù)雜數(shù)學(xué)物理問題的重要方法，可以利用變分原理將實(shí)際工程計(jì)算中的微分方程問題轉(zhuǎn)化成與之等價(jià)的變分問題.因此，有些工程問題可以轉(zhuǎn)化為控制方程，由基本方程和平衡方程給出.從變分原理出發(fā)，給出相應(yīng)問題的有限元模型，采用有限元法給出相應(yīng)的微分方程邊值問題的近似解，進(jìn)而考慮滿足高精度要求的數(shù)值結(jié)果.

一、一維定長水槽問題

1.位勢流

在地下水問題中，設(shè)流體為定長不可壓流動(dòng)，由連續(xù)方程或質(zhì)量守恒方程描述此問題.設(shè)面積為常數(shù)，則一元?jiǎng)莺瘮?shù)可為φ(x)=-K(x)h(x)與u(x)=玠φ[]玠玿，其中u為速度，h為壓力頭，L為滲透率.基本方程由獶arcy定律［4］猽(x)=-K·玠玥(x)[]玠玿給出.定長不可壓縮流動(dòng)滿足玠玼[]玠玿=0關(guān)系，聯(lián)立上述三個(gè)方程有玠2φ[]玠玿2=0，其解為線性函數(shù)，因而一維流動(dòng)問題的速度為常數(shù).

2.質(zhì)量輸運(yùn)方程

一維水槽中介質(zhì)流動(dòng)是定常的，其控制方程將發(fā)生擴(kuò)散現(xiàn)象，質(zhì)量方程的平衡以稀釋混合項(xiàng)的形式給出.將位勢流理論與質(zhì)量擴(kuò)散理論相結(jié)合，得到問題的一個(gè)較全面的描述.設(shè)截面積為常數(shù)，則稀釋混合項(xiàng)的質(zhì)量平衡描述為

u(x)玠獵(x)[]玠玿+玠玧(x)[]玠玿+K璻C(x)=m.(1)

這里u(x)表示流動(dòng)速度，C(x),j(x)分別表示濃縮系數(shù)和流通量，K璻表示稀釋項(xiàng)與周圍物質(zhì)的反應(yīng)速度，m表示質(zhì)量源函數(shù)，其基本方程稱為獸ick定律［5-6］，表示為j(x)=-D(x)·玠獵(x)[]玠玿，這里D(x)是擴(kuò)散系數(shù)，聯(lián)立可得控制方程：

u(x)玠獵(x)[]玠玿-玠玔]玠玿D(x)玠獵(x)[]玠玿+K璻C(x)=m.(2)

設(shè)已知速度u(x)，C(x)邊界條件為本質(zhì)邊界條件，﹋(x)邊界條件為自然邊界條件.對(duì)方程(2)用擬變分法［2］，其擬變分函數(shù)可寫為

J(C)=А要璙1[]2D玠獵[]玠玿2+Cu玠獵[]玠玿+K璻C2-2mC玠玍.(3)

將體積積分改為關(guān)于定常截面積A和長度L的重積分，并將(3)改寫為基于形函數(shù)的矩陣形式，可得

J(C)=A[]2А要琇0{C}琓玠玁[]玠玿

琓［D］玠玁[]玠玿

{C}玠玿+A[]2А要琇0{C}琓·［Ｎ］Ｔ［u］玠玁[]玠玿{C}玠玿+A[]2

А要琇0{C}琓［N］琓［K璻］［N］{C}玠玿-И〢∫琇0{C}{N}琓m玠玿.

二、數(shù)值應(yīng)用及分析

設(shè)一維水槽長L=1 玬，各處的橫截面相同，內(nèi)有流體流動(dòng)，在x=0處，高度0.33 玬，在x=1處，高度0.11 玬，在x=0處測得某化學(xué)污染物的濃度為0，而在x=1處的濃度為10 玬g/m3.在流動(dòng)過程中，污染物與周圍物質(zhì)以K璻=6.6×10-4的常數(shù)比發(fā)生反應(yīng)，化學(xué)物質(zhì)以擴(kuò)散系數(shù)D=1.0×10-4擴(kuò)散，物質(zhì)的水力傳導(dǎo)率K=1.0×10-4，我們分析水槽中化學(xué)污染物的穩(wěn)態(tài)分布情況.

構(gòu)造有限元解，將剛度矩陣代入非零本質(zhì)邊界條件后，可得整體剛度矩陣：

1.0[]

[]22.0[]-11.0

[]-11.0[]22.0[]-11.0

[3]鱗]鱗]鳘

[4]鱗]鱗]鳘

[5]-11.0[]22.0[]-11.0

[6]-11.0[]22.0

[8]1.0

h1

h2h3う螵う螵h10h11h12

=0.33

3.630.0う螵う螵0.01.210.11

用列主元獹auss消去法求解得h1=0.33,h2=0.31,﹉3=0.29,h4=0.27,h5=0.25,h6=0.23,h7=0.21,h8=0.19,h9=0.17,h10=0.15,h11=0.13,h12=0.11,有﹗=-K(h2-h1)[]L[]n，

（10-4）11.09 -10.79

-11.01 11.31

C1C2

(4)

由式(4)邊界條件C﹛=0=0和C﹛=1=10，組成如下整體矩陣:

（10-4)·104

[]22.4[]-10.79

[]-11.01[]22.4[]-10.79

[3]鱗]鱗]鳘

[4]-11.01[]22.4[]-10.79

[5]-11.01[]22.4

[7]104

C1

C2C3う螵C10C11C12

=0.0

0.00.0う螵0.010.79·10-310.0

.

邊界條件為C﹛=0=0和C﹛=L=C璍，長度為L的區(qū)域上解析解為C=C璍e│(x-L)·玸in玥(βx)[]玸in玥(βL)，其中α=u[]2D和│=α2+K璻[]D1[]2.下表是由有限元解C環(huán)和解析解C制成的：

通過編程計(jì)算，給出n=100,200,400,800時(shí)的有限元解與真解的絕對(duì)誤差及L2范數(shù)А篇琻﹊=1(C-C環(huán))2

1[]2.通過數(shù)值實(shí)驗(yàn)可知，用有限元法作數(shù)值模擬時(shí)，當(dāng)網(wǎng)格剖分?jǐn)?shù)越來越大由100到800，絕對(duì)誤差的數(shù)量級(jí)可達(dá)10-6，L2范數(shù)的數(shù)量級(jí)可達(dá)10-5，其收斂階高于1階玶ate=log2E璑[]E2N.

通過研究一維固定長度的水槽化學(xué)污染物的穩(wěn)態(tài)分布情況，利用有限元水槽化學(xué)污染的模型作數(shù)值模擬，我們得到了各節(jié)點(diǎn)上高精度的污染物分布情況.

【參考文獻(xiàn)】オ

［1］Reddy J.N.有限元法概論［M］.鄒仲康，譯.長沙：湖南科學(xué)技術(shù)出版社，1988：173-189.

［2］李開泰，黃艾香，黃慶懷.有限元方法及其應(yīng)用［M］.北京：科學(xué)出版社，2006：20-45.

［3］Brenner S.C.,Scott L.R.The Mathematical Theory of Finite Element Methods［M］.北京：世界圖書出版公司，2008: 155-170.

［4］Babuska I.and Osborn E.Genera Lized finite element methods: Their performance and their relation to mixed methods［J］.SIAM J.Numer.Anal.,1983,20: 510-536.

［5］Melenk J.M.,Babuska I.The partition of unity finite element methods: basic theory and applications［J］.Comput.Methods Appl.Mech.Engrg.,1996,139(1/4): 289-314.

［6］Huang Yun-qing,Xu Jin-chao.A conforming finite element method for overlapping and nonmatching grids［J］.Math Comput.,2002,72(243):1057-1066.