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

丁 曉，陳 璐，江 山

（南通大學 理學院，江蘇 南通 226019）

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

1 Stokes 方程與泛函空間

考慮二維矢量型Stokes 方程

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

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

再進行分部積分，有

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

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

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

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

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

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

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

使得對任意i，j=1，…，6 都有

比如，對于?1，將三角形單元的六個點代入，有

解六元一次方程組，可得

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

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

3 數(shù)值算例

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

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

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

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

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

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

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

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

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