非穩態二維擴散方程的高次有限元黃金比例有限差分格式求解

孫美玲, 丁 曉, 江 山*

(1. 南通大學理學院, 江蘇 南通 226019; 2. 南通職業大學基礎部, 江蘇 南通 226007)

非穩態對流擴散反應問題廣泛存在于流體力學、空氣動力學、化工制造、生態環境及運籌控制等領域, 其模型方程含有時間導數項、空間二階/一階導數項和右端非穩態項, 其中的非對稱項增加了理論分析與數值求解的難度．由于非穩態對流擴散反應方程的解析解通常難以得到, 所以有效且實用的數值求解方法便顯得尤為重要．近年來, 非穩態擴散問題常用的數值求解方案主要有有限差分法[1-3]、有限元法[4-5]、間斷有限元法[6-7]、有限體積法[8-9]、時空有限元法[10]和虛擬元方法[11]等．例如, Gupta等[1]僅構造了單一的有限差分格式求解含參數的拋物型對流擴散問題; Lin等[4]提出對流擴散反應方程的一種穩定化有限元逼近方法, 但計算量偏大．本文擬構造高次有限元并結合有限差分的數值隱格式求解非穩態對流擴散問題．在選擇有限差分參數θ時不再簡單使用Euler向后差分或Crank-Nicolson六點對稱等格式, 而是選擇具備美學特性的黃金分割比例差分格式進行時空間的全離散, 以期實現數值計算的高階精度與穩定收斂．

1 模型方程及其變分形式

考慮二維非穩態對流擴散反應方程[4-5]

(1)

(2)

其中n為外法線向量, ds表示曲線積分.

令檢驗函數v(x,y)在?Ω上為零．于是, 問題(1)轉化為: 尋找試探函數u(x,y,t),須滿足

(3)

(ut,v)+a(u,v)=(f,v)

(4)

2 空間尺度的有限元格式

在有限元格式下完成變分形式(4)的與空間尺度有關的半離散化格式[4-5]．在空間區域Ω的x、y方向進行一致剖分, 并在剖分后的單元上構造合適的有限元基函數φj(j=1,2,…,Nb), 其中Nb為全局基函數的總數目．以基函數φj為基底張成有限維逼近空間Uh, 且Uh?H1(Ω)．分別選用二維Lagrange型線性元和二次元作為有限維空間的基函數進行比較．為了闡述清晰, 以更復雜的二維二次元基函數為例, 進行空間尺度下的有限元空間構造．

(5)

在δ-準則限定下, 可求得參考單元的6個局部基函數為

(6)

實際單元的頂點記為(xi,yj)(i=1,2,3;j=1,2,3), 通過仿射變換建立參考單元與實際單元之間的對應關系:

(7)

(8)

其中對應的仿射變換為

(9)

以式(8)所示的局部基函數ψni(x,y)為基底張成局部有限元空間Sh(En)=span{ψn1,ψn2,ψn3,ψn4,ψn5,ψn6}．在剖分形成的每一個全局節點xj(j=1,2,…,Nb)對應定義全局基函數φj, 使得i,j=1,2,…,Nb時有φj|En∈Sh(En)且滿足

(10)

令Tb為所有單元的全局節點對應的數據索引, 則

(11)

(uht,vh)+a(uh,vh)=(f,vh),

(12)

進一步可得

(13)

(14)

求解得到待定系數uj(t)．令

(15)

MX′(t)+A(t)X(t)=b(t).

(16)

由于基函數φj與檢驗函數φi只與空間變量x,y有關,A(t)及b(t)隨時間t的變化關系取決于系數α、β、γ和源項f, 且在某個確定時刻A(t)和b(t)退化為僅依賴于空間．

3 時間尺度的黃金比例有限差分格式

采用有限差分格式[12]處理時間尺度．有限差分格式的類型由參數θ的取值確定, 如θ=0時, 為Euler向前差分顯格式;θ=1時, 為Euler向后差分隱格式;θ=0.5時, 為Crank-Nicolson六點對稱隱格式．本文將突破固有的差分格式, 通過構造與θ取值有關的全離散格式, 研究不同取值時差分格式下的優化數值精度結果．

(17)

當式(15)中的系數α,β,γ與時間無關時,A(t)不依賴于時間, 簡記為A．定義Xm+θ=θXm+1+(1-θ)Xm, 有

(18)

將式(18)代入全離散格式(17), 化簡可得

(19)

定義

(20)

(21)

將式(20)(21)代入式(19), 即得全新的與時間有關的全離散線性代數系統:

BXm+θ=d,m=0,1,…,Nm-1.

(22)

求解式(22)可得未知向量組Xm+θ．

本文將數學與美學相結合, 對參數θ取黃金分割比例值0.618, 數值驗證黃金比例下全離散格式的數值精度優勢．

4 數值算例

運用MATLAB 2020b軟件編寫數值算例程序, 硬件環境為Intel Core i9-10900K CPU 3.70 GHz, 32 GB RAM．

設置非穩態問題(1)中系數α=γ=1,β=(1,1)T, 則其精確解為

u(x,y,t)=-cos(πy)ex+t,

(23)

(24)

對于具有精確解的二維非穩態對流擴散反應問題(24), 在空間尺度上使用有限元格式的線性基函數與二次基函數分別進行半離散化, 在時間尺度上取θ=0.618的黃金比例有限差分格式進行全離散化．為了衡量新的混合格式下數值解uh對真解u的逼近效果, 計算兩者誤差的L∞范數、L2范數與H1半范數:

‖u-uh‖∞=supx,y∈Ω|u(x,y,t)-uh(x,y,t)|,

(25)

(26)

(27)

表1 θ=0.618時, 線性有限元誤差的L∞范數、L2范數、H1半范數和收斂階

表2 θ=0.618時, 二次有限元誤差的L∞范數、L2范數、H1半范數和收斂階

為更直觀地對比, 圖1給出了θ=0.618, 空間剖分數N=24,t=0.25 s時方程(1)的真解與線性有限元解、二次有限元解之間的離散節點絕對誤差．由圖1可知: 在t=0.25 s,N=24時, 線性有限元的最大誤差數量級為10-4, 而二次有限元的最大誤差數量級為10-6且誤差遠低于線性有限元的,表明相同條件下二次有限元格式的精度優于線性有限元．

圖1 θ=0.618, t=0.25 s時, 線性有限元(a)與二次有限元(b)在N=24下的絕對誤差Fig.1 When θ=0.618, t=0.25 s, absolute error of linear finite elements (a) and quadratic finite elements (b) on partition number N=24

5 結語

本文提出空間高次有限元法結合時間全離散的黃金比例有限差分隱格式, 精確且高效求解了非穩態二維對流擴散反應問題．構造半離散化時的二次有限元,并系統建立全離散化的時間參數θ格式的代數系統, 通過3種范數的誤差估計數值驗證了黃金比例差分格式下二次有限元具備更高的計算精度和收斂階數．本文方法為后續研究提供了一種新思路,可以通過合理選取θ值改善全離散格式對精度的影響．