高次有限元格式數(shù)值求解Stokes 方程

丁 曉，陳 璐，江 山

（南通大學(xué) 理學(xué)院，江蘇 南通 226019）

20 世紀(jì)50 年代末，受大型水壩計(jì)算問題的啟發(fā)，馮康院士創(chuàng)建了系統(tǒng)化、理論化的有限元法。隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，有限元法已廣泛應(yīng)用于求解熱傳導(dǎo)、電磁場(chǎng)、流體力學(xué)等連續(xù)介質(zhì)問題。Stokes 方程可以描述雷諾數(shù)很小時(shí)液滴在黏性流體中的運(yùn)動(dòng)問題，近年來(lái)被廣泛關(guān)注。利用有限元法求解Stokes 方程一直是計(jì)算數(shù)學(xué)和流體力學(xué)領(lǐng)域的熱點(diǎn)，其研究取得了諸多成果[1-7]。比如：文獻(xiàn)[1]考慮非定常Stokes 方程，將流函數(shù)方程和渦度方程分開討論，再用混合有限元法進(jìn)行誤差估計(jì)；王小軍和祝家麟[2]針對(duì)含開邊界的二維Stokes 問題研究Galerkin 邊界元解法，數(shù)值模擬了不可壓黏性流體的繞流；吳妍等[3]使用自適應(yīng)變分多尺度方法對(duì)Stokes 方程進(jìn)行優(yōu)化，將細(xì)尺度問題分解為若干互不影響的子問題再進(jìn)行求解；段獻(xiàn)葆等[4]基于Stokes 問題，提出了一種與優(yōu)化方法相耦合的自適應(yīng)網(wǎng)格方法；文獻(xiàn)[7]借助差商法研究了非線性穩(wěn)態(tài)Stokes 系統(tǒng)的弱解對(duì)高階分?jǐn)?shù)的可微性；等等。此外，作者課題組應(yīng)用有限元法與多尺度計(jì)算，針對(duì)二維奇異攝動(dòng)邊界層問題[8]、非穩(wěn)態(tài)擴(kuò)散方程[9]、彈性力學(xué)方程[10]等進(jìn)行了一系列富有成效的理論分析和數(shù)值驗(yàn)證。但上述成果均未涉及Stokes 方程求解，且所得精度或效率仍有提升空間。本文主要針對(duì)Dirichlet 邊界下的二維矢量型Stokes 方程，基于有限維逼近無(wú)限維的數(shù)值計(jì)算思想和變分虛功原理，將原方程區(qū)域離散化分割成有限個(gè)單元，再對(duì)每個(gè)局部單元選定二次基函數(shù)，形成并組裝信息矩陣和載荷向量，進(jìn)而針對(duì)矢量型方程求解代數(shù)方程組，并進(jìn)行圖像繪制和誤差分析，即采用高次有限元格式處理Stokes 方程，以期得到更好的數(shù)值精度、更高的穩(wěn)定收斂階數(shù)。

1 Stokes 方程與泛函空間

考慮二維矢量型Stokes 方程

類似地，應(yīng)力張量也可寫為

二維區(qū)域Ω=[0，1] × [-0.25，0]，?Ω 是其邊界。

再進(jìn)行分部積分，有

2 有限元格式與二次基函數(shù)構(gòu)造

有限元法的基本思想是用有限維空間去逼近無(wú)限維空間，在保證數(shù)值解正確性的前提下，進(jìn)一步考慮數(shù)值解的穩(wěn)定性與收斂性。

下面構(gòu)造所需的有限元空間。針對(duì)流速與壓力，分別考慮有限元空間Uh?H1（Ω）和Wh?L2（Ω），記Uh0是由Uh在邊界上值為零的所有函數(shù)張成的空間。有限元法的變分形式為：

接著，對(duì)區(qū)域Ω 進(jìn)行網(wǎng)格剖分，用三角形單元將其離散化。對(duì)每個(gè)局部單元選取二次基函數(shù)，形成并組裝信息矩陣A 和載荷向量，求解代數(shù)方程組，得到流速的有限元解。

為了取得更優(yōu)越的模擬效果，在局部單元上采用二次元基函數(shù)。設(shè)初始的三角形單元E=ΔA1A2A3，再取三條邊的中點(diǎn)A4、A5、A6，見圖1。

圖1 線性單元3 節(jié)點(diǎn)與二次單元6 節(jié)點(diǎn)的對(duì)應(yīng)關(guān)系

定義參考單元上的6 個(gè)二次元基函數(shù)

使得對(duì)任意i，j=1，…，6 都有

比如，對(duì)于?1，將三角形單元的六個(gè)點(diǎn)代入，有

解六元一次方程組，可得

所以，?1=2x2+2y2+4xy-3x-3y+1。類似地，經(jīng)上述步驟可以求得全部6 個(gè)基函數(shù)，即

由此可使總剛度矩陣的稠密程度更高，每個(gè)單元與節(jié)點(diǎn)間的聯(lián)系相對(duì)線性基函數(shù)而言也更緊密，二次基函數(shù)的有限元解能更好地反映Stokes方程的精確解分量。

3 數(shù)值算例

為驗(yàn)證有限元法求解二維Stokes 方程的可行性與有效性，通過實(shí)際算例和MATLAB 編程計(jì)算相應(yīng)的數(shù)值精度和收斂階數(shù)。

令黏度系數(shù)v=1，在區(qū)域Ω=[0，1]×[-0.25，0]，邊界?Ω 上有則右端為

為了直觀展現(xiàn)有限元解與精確解之間的近似程度，用MATLAB 畫出單方向剖分?jǐn)?shù)N=32 時(shí)的精確解的兩個(gè)分量u1，u2如圖2，二次有限元解的兩個(gè)近似分量u1h，u2h如圖3，以及其絕對(duì)誤差的三維圖如圖4。由圖可見，有限元解與精確解非常接近，且u1的最大誤差數(shù)量級(jí)為3×10-6，u2的最大誤差數(shù)量級(jí)為6×10-7，驗(yàn)證了有限元解逼近精確解的可靠性。

圖2 精確解的分量u1，u2 圖像

圖3 有限元解的分量u1h，u2h 圖像

圖4 精確解與有限元解之間的絕對(duì)誤差三維圖像

通過計(jì)算得到三種不同范數(shù)意義下有限元解與精確解之間的數(shù)值誤差和收斂階，見表1。由表1 可見：二次元基函數(shù)格式下該算例的L∞范數(shù)和L2范數(shù)都達(dá)到了三階收斂；H1半范數(shù)達(dá)到了二階收斂，并且最高精度可以實(shí)現(xiàn)10-8數(shù)量級(jí)。數(shù)值實(shí)踐結(jié)果充分表明，高次有限元格式對(duì)于數(shù)值求解Stokes 方程十分有效，得到的結(jié)果也很理想。

表1 二次有限元格式的誤差范數(shù)及其收斂階數(shù)

綜上可知，通過高次有限元格式的理論構(gòu)造和數(shù)值模擬，可直觀系統(tǒng)地展示二次元基函數(shù)下有限元法求解二維矢量型Stokes 方程的有效性和精確性，得到了高精度、高收斂的一致穩(wěn)定化數(shù)值模擬結(jié)果。