非穩(wěn)態(tài)二維擴(kuò)散方程的高次有限元黃金比例有限差分格式求解

孫美玲, 丁 曉, 江 山*

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

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

1 模型方程及其變分形式

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

(1)

(2)

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

令檢驗(yàn)函數(shù)v(x,y)在?Ω上為零．于是, 問(wèn)題(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方向進(jìn)行一致剖分, 并在剖分后的單元上構(gòu)造合適的有限元基函數(shù)φj(j=1,2,…,Nb), 其中Nb為全局基函數(shù)的總數(shù)目．以基函數(shù)φj為基底張成有限維逼近空間Uh, 且Uh?H1(Ω)．分別選用二維Lagrange型線性元和二次元作為有限維空間的基函數(shù)進(jìn)行比較．為了闡述清晰, 以更復(fù)雜的二維二次元基函數(shù)為例, 進(jìn)行空間尺度下的有限元空間構(gòu)造．

(5)

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

(6)

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

(7)

(8)

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

(9)

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

(10)

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

(11)

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

(12)

進(jìn)一步可得

(13)

(14)

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

(15)

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

(16)

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

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

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

(17)

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

(18)

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

(19)

定義

(20)

(21)

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

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

(22)

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

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

4 數(shù)值算例

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

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

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

(23)

(24)

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

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

(25)

(26)

(27)

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

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

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

圖1 θ=0.618, t=0.25 s時(shí), 線性有限元(a)與二次有限元(b)在N=24下的絕對(duì)誤差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é)語(yǔ)

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