孫美玲, 丁 曉, 江 山*

(1. 南通大學(xué)理學(xué)院, 江蘇 南通 226019; 2. 南通職業(yè)大學(xué)基礎(chǔ)部, 江蘇 南通 226007)

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

1 模型方程及其變分形式

考慮二維非穩(wěn)態(tài)對流擴散反應(yīng)方程[4-5]

(1)

(2)

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

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

(3)

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

(4)

2 空間尺度的有限元格式

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

(5)

在δ-準(zhǔn)則限定下, 可求得參考單元的6個局部基函數(shù)為

(6)

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

(7)

(8)

其中對應(yīng)的仿射變換為

(9)

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

(10)

令Tb為所有單元的全局節(jié)點對應(yīng)的數(shù)據(jù)索引, 則

(11)

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

(12)

進一步可得

(13)

(14)

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

(15)

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

(16)

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

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

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

(17)

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

(18)

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

(19)

定義

(20)

(21)

將式(20)(21)代入式(19), 即得全新的與時間有關(guān)的全離散線性代數(shù)系統(tǒng):

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

(22)

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

本文將數(shù)學(xué)與美學(xué)相結(jié)合, 對參數(shù)θ取黃金分割比例值0.618, 數(shù)值驗證黃金比例下全離散格式的數(shù)值精度優(yōu)勢．

4 數(shù)值算例

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

設(shè)置非穩(wěn)態(tài)問題(1)中系數(shù)α=γ=1,β=(1,1)T, 則其精確解為

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

(23)

(24)

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

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

(25)

(26)

(27)

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

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

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

圖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 結(jié)語

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