非定常對(duì)流占優(yōu)擴(kuò)散方程的龍格庫塔伽遼金有限元方法

馮立偉, 席 偉

(1.沈陽化工大學(xué) 數(shù)理系，遼寧 沈陽 110142； 2.沈陽化工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，遼寧 沈陽 110142)

大量的自然現(xiàn)象的理論模型是對(duì)流擴(kuò)散方程，如海水中熱量和鹽分的變化、大氣和水體中污染物的流動(dòng)分布、藥物在人體中的流動(dòng)和稀釋代謝等[1-3].所以，對(duì)流擴(kuò)散方程的求解方法研究早已成為重要的研究領(lǐng)域.有限差分法、譜元法、有限元法是求解對(duì)流擴(kuò)散方程的典型方法[4-7].但是，在對(duì)流占優(yōu)擴(kuò)散方程中，由于對(duì)流現(xiàn)象起主要作用，擴(kuò)散現(xiàn)象比較微弱，導(dǎo)致傳統(tǒng)有限元方法產(chǎn)生偽數(shù)值振蕩[8-10]，使得求解失敗.流線擴(kuò)散有限元(streamline diffusion,SD)[8,10]通過在檢驗(yàn)函數(shù)上沿流場(chǎng)方向添加人工黏性項(xiàng)，提高其解的穩(wěn)定性，但其未考慮時(shí)間效應(yīng)，只能求解穩(wěn)態(tài)對(duì)流擴(kuò)散問題.對(duì)于發(fā)展型對(duì)流擴(kuò)散問題，由于物質(zhì)在空間的分布隨時(shí)間發(fā)生變化，經(jīng)典求解方法為時(shí)空有限元，但其將時(shí)間項(xiàng)與空間項(xiàng)同等對(duì)待，導(dǎo)致求解維數(shù)較高，顯著增加了算法的計(jì)算復(fù)雜度.另一求解方法為差分流線擴(kuò)散有限元[11-13](finite difference streamline diffusion,FDSD)，其在時(shí)間方向上采用差分離散，空間上使用流線擴(kuò)散有限元，成功抑制了數(shù)值振蕩的產(chǎn)生.但是，由于其在檢驗(yàn)函數(shù)上添加流線項(xiàng)，導(dǎo)致變分方程的每一項(xiàng)都對(duì)應(yīng)添加了流線項(xiàng)，使得求解格式復(fù)雜，計(jì)算量顯著增大.本文對(duì)非定常對(duì)流占優(yōu)擴(kuò)散方程給出一種新的求解方法：首先，在擴(kuò)散項(xiàng)上引入擬合因子來抑制偽數(shù)值振蕩；其次，使用伽遼金有限元方法將方程轉(zhuǎn)化為半離散化常微分方程；最后，使用四階四級(jí)龍格庫塔方法將常微分方程轉(zhuǎn)變?yōu)榫€性方程組，從而獲得其數(shù)值解.對(duì)求解方法進(jìn)行了理論分析，并給出了詳細(xì)的計(jì)算過程.使用一個(gè)數(shù)值例子和天然氣管道泄漏實(shí)例進(jìn)行了實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果驗(yàn)證了方法的有效性.

1 計(jì)算格式描述

非定常對(duì)流占優(yōu)擴(kuò)散方程的初邊值問題:

(1)

(f,v)，?v∈V.

定義兩個(gè)泛函：

此變分問題的Galerkin有限元離散為求uh(·;t)∈Uh，使得

(2)

成立.將式(2)改寫為半離散化常微分方程組

(3)

對(duì)此方程組再采用四階四級(jí)龍格庫塔法離散時(shí)間導(dǎo)數(shù)項(xiàng)，得到全離散化線性方程組

按時(shí)間層，逐層推進(jìn)可求得各時(shí)間層上每個(gè)網(wǎng)格節(jié)點(diǎn)上數(shù)值解.

2 詳細(xì)計(jì)算過程

2.1 空間Galerkin有限元半離散格式

對(duì)空間區(qū)域進(jìn)行三角形網(wǎng)格剖分，獲得網(wǎng)格節(jié)點(diǎn)坐標(biāo)、單元信息以及邊界節(jié)點(diǎn)和邊界單元信息，并對(duì)各單元和節(jié)點(diǎn)分別進(jìn)行編號(hào).在第k個(gè)剖分單元ek上使用一次線性形狀函數(shù)Ni、Nj、Nm，則

uh(x,y,t)=Ni(x,y)Ui(t)+Nj(x,y)Uj(t)+

Nm(x,y)Um(t),(x,y)∈e.

上式簡(jiǎn)記為uh=NUe(t),同理vh=NVe，記uh=BUe,vh=BVe，代入近似變分方程(2)，方程(2)中的各項(xiàng)計(jì)算如下：

cNTN+εBTB)dxdy

)Ue+

其中：Ue=CU,Ve=CV，Me、Ke、Fe分別為單元e質(zhì)量矩陣、單元?jiǎng)偠染仃嚭蛦卧奢d向量，M、K、F分別為總質(zhì)量矩陣、總剛度矩陣和總荷載向量.此時(shí)，變分方程(2)變?yōu)?/p>

因?yàn)閂是任意的，所以上式變?yōu)?/p>

2.2 單元計(jì)算及總體合成

Ke=?e(NTb·B+cNTN+εBTB)dxdy+

由于單元質(zhì)量矩陣中被積函數(shù)為多項(xiàng)式函數(shù)，可采用歐拉積分計(jì)算.由于源項(xiàng)形式不唯一，依賴于具體求解問題，在單元?jiǎng)偠染仃嚭蛦卧奢d向量中與其有關(guān)的積分，可以采用高斯型積分公式

在求得單元質(zhì)量矩陣、單元?jiǎng)偠染仃嚭蛦卧奢d向量后，并按本單元所包含的節(jié)點(diǎn)的編號(hào)將其疊加獲得總質(zhì)量矩陣M、總剛度矩陣K和總荷載向量F，從而得到常微分方程組式(3).

2.3 常微分方程組離散求解

將常微分方程組(3)改寫為

為了獲得高精度解，對(duì)此常微分方程采用四階四級(jí)龍格庫塔有限差分法離散為

(4)

其中：k1=(Fn-1-KUn-1)，

全離散化線性方程組(4)可簡(jiǎn)記為

AUn=bn-1.

(5)

2.4 方程初邊界條件的加載

2.5 計(jì)算程序的編寫

根據(jù)上面的詳細(xì)計(jì)算過程描述，給出了求解方法的算法流程(圖1)，通過使用MATLAB編程實(shí)現(xiàn).根據(jù)具體問題，編寫求解區(qū)域的幾何描述矩陣g，將求解區(qū)域進(jìn)行三角形剖分，獲得三角形單元矩陣t、節(jié)點(diǎn)坐標(biāo)矩陣p和邊界單元矩陣e，可使用MATLAB函數(shù)initmesh(g,h)實(shí)現(xiàn).編寫主函數(shù)MassStiffnessLoadCalculate計(jì)算各個(gè)時(shí)刻的總質(zhì)量矩陣、總剛度矩陣和總荷載向量，并對(duì)它們進(jìn)行總體組裝，獲得求解常微分方程的龍格庫塔法中系數(shù)及右端，再按照式(4)獲得線性方程組式(5)，最后線性方程組采用雙共軛梯度法進(jìn)行求解.其中主函數(shù)通過調(diào)用子函數(shù)ElementMassStiffnessLoadCalculate實(shí)現(xiàn)單元上的質(zhì)量矩陣、剛度矩陣和荷載向量的計(jì)算.另外，單元計(jì)算上述兩類積分，分別編寫函數(shù)EulerIntegral、GaussIntegral2D和GaussIntegral1D實(shí)現(xiàn)標(biāo)準(zhǔn)三角形單元上的高斯重積分和標(biāo)準(zhǔn)邊界單元上的高斯線積分.編寫函數(shù)DealFirstBounCond采用置1法實(shí)現(xiàn)第一類邊界條件的加載.先劃掉系數(shù)矩陣A的該節(jié)點(diǎn)對(duì)應(yīng)行和列中的元素，并按上述修正右端列向量b的其余行，得到最終的線性方程組.使用pdemesh和pdecont函數(shù)分別繪制最終時(shí)刻的解及解的等值線圖像.

圖1 算法流程

3 數(shù)值實(shí)驗(yàn)

3.1 數(shù)值算例

Ω=[0,1]×[0,1].

其中：b=(1 1),c=1,f=e-ty(1-y)(1-2x+2ε)+e-tx(1-x)(1-2y+2ε)，ε=0.001.其解析解為u=e-txy(1-x)(1-y).采用三角形網(wǎng)格，空間網(wǎng)格最大尺寸為0.05，時(shí)間步長(zhǎng)為0.01，計(jì)算到最終時(shí)刻的結(jié)果見圖2，誤差見圖3.在不同空間和時(shí)間剖分尺寸下，最終時(shí)刻的總體相對(duì)誤差見表1.

圖2 數(shù)值解與解析解

圖3 最終時(shí)刻的誤差

表1 相對(duì)誤差

3.2 管道泄漏模擬

對(duì)半徑為1 m的天然氣管道的泄漏進(jìn)行數(shù)值模擬，天然氣泄漏的水平速度和豎直速度b=[1,0.7]m/s，擴(kuò)散系數(shù)ε=0.5，在管道45°角處有一半徑為0.1 m的缺口，以f=300 m/s的速度向外泄漏天然氣，使用時(shí)間步長(zhǎng)τ=0.001 s，空間最大網(wǎng)格尺寸h=0.1 m，圖4為四分之一管道周圍空間的網(wǎng)格剖分結(jié)果，求得的各個(gè)時(shí)刻的天然氣濃度等值線如圖5所示.從圖5可看出：計(jì)算方法能夠較好地模擬天然氣的泄漏.

圖4 天然氣泄漏模型網(wǎng)格剖分

圖5 各個(gè)時(shí)刻的天然氣濃度等值線

4 結(jié) 論

通過在對(duì)流項(xiàng)上添加擬合因子消除對(duì)流占優(yōu)的不利影響，并聯(lián)合使用四階四級(jí)龍格庫塔方法和伽遼金有限元方法，得到了對(duì)流占優(yōu)擴(kuò)散方程的一種求解方法.對(duì)方法進(jìn)行了詳細(xì)闡述，并進(jìn)行了數(shù)值實(shí)驗(yàn)和天然氣泄漏實(shí)驗(yàn).擬合因子抑制了偽數(shù)值振蕩，成功實(shí)現(xiàn)了對(duì)流占優(yōu)擴(kuò)散方程的求解.