TVD格式求解對流擴散方程的最優(yōu)限制器研究

戚園春，劉 昉，侯慶志，2

(1. 天津大學建筑工程學院，天津 300354；2. 天津大學智能與計算學部，天津 300354)

1 引言

對流擴散方程作為描述自然界流動、傳熱及物質(zhì)輸運過程的主要控制方程，其可充分反映流體在流動過程中所攜帶的物理量的變化規(guī)律，在流體力學、能源開發(fā)及環(huán)境科學等領域發(fā)揮著重要的作用[1]。

對流擴散方程的求解方法主要有有限差分法、有限元法和有限體積法[3-6]。有限體積法物理意義明確、計算效率較高，基于該方法求解對流擴散方程時，常用的離散格式有一階迎風、二階迎風、中心差分、Lax-Friedrichs、Lax-Wendroff、Beam-Warming、Fromm等，然而這些格式會導致計算結(jié)果出現(xiàn)不同程度的數(shù)值耗散或者數(shù)值振蕩[7]，精度過低會抹平數(shù)值解中的間斷，精度過高又會導致計算結(jié)果在間斷附近出現(xiàn)偽振蕩。近年來，高分辨率格式的發(fā)展有效緩和了這一問題，其可敏銳地捕捉間斷，提高激波的分辨率并消除偽振蕩。常用的高分辨率離散方法有FCT (Flux Corrected Transport)方法和TVD (Total Variation Diminishing)方法。

FCT方法是幾乎所有保單調(diào)和非振蕩流體輸運算法的起源，該方法的本質(zhì)是在輸運格式里面施加一個擴散和反擴散項，使之能平滑激波區(qū)和流場陡變區(qū)，從而避免偽振蕩。而TVD方法是通過保證變量的變化不隨時間增加，從而消除偽振蕩，其在對流擴散方程求解中表現(xiàn)出了良好的性能。早期，Harten[8]依據(jù)TVD性質(zhì)構(gòu)造了一個通量限制函數(shù)，并與高階格式相結(jié)合，首次提出了“高分辨率方法”這一概念。Tamamidis等[9]采用3種高階單調(diào)格式分別求解了4類二維對流擴散問題，并與無通量限制器格式的計算結(jié)果進行對比，表明帶通量限制器的高階格式可顯著提高計算結(jié)果的準確性。Hutter等[10]采用多種數(shù)值格式求解一維對流擴散問題，表明帶有Superbee限制器的MTVDLF方法是求解對流占優(yōu)問題的最佳方法。郭建紅等[11]基于有限體積通量修正，應用高階迎風對流格式和帶TVD限制器的高分辨率格式研究了自然空泡湍流流動的數(shù)值計算方法，發(fā)現(xiàn)后者能更好地捕捉空泡界面附近物理量的階躍特性。馮定華等[12]應用8種差分格式和5種限制器計算分析了激波管Riemanna問題，并對比了各類限制器的壓縮性和耗散性。Zhang等[13]對現(xiàn)有的TVD格式進行歸類，表明合理的通量限制器是TVD格式的關鍵要素，并在此基礎上提出了一種改進的用于穩(wěn)態(tài)計算的SS-TVD(Steady-state TVD)限制器TCDF。

因此，為進一步優(yōu)選求解不同類型對流擴散問題的最優(yōu)算法，本文共采用了3種典型的TVD格式：MUSCL、TVDLF、MTVDLF格式，10種常用的通量限制器：Superbee[14]、Minmod[15]、Woodward[10]、Van Leer[9]、UMIST[16]、WACEB[17]、Albada[18]、OSPRE[19]、TCDF[13]、Koren[13]，求解了階躍型純對流問題、階躍型對流擴散問題、線性高斯型對流擴散問題、擬線性高斯型對流擴散問題、Burgers方程，共5類典型的對流擴散問題。同時，為了促進該問題的數(shù)值模擬研究發(fā)展，本文所有格式的相關代碼可在以下網(wǎng)站自行下載：http：∥114.55.218.152/wp-blog/，以期為國產(chǎn)自主軟件開發(fā)提供一定的技術(shù)支持。

2 數(shù)學模型

對流擴散方程的一維形式和守恒形式為

(1)

(2)

3 TVD方法與通量限制器

有限體積法守恒性較好，物理意義明確，以下主要采用該方法離散方程(2)中的導數(shù)，通過在[xj-1/2，xj+1/2]×[tn，tn+1]上積分，可以得到該方程的離散形式為

(3)

(4)

(5)

(6)

(7)

(8)

(9)

φj=φ(θj)

(10)

(11)

式中，σj為斜率限制器，θj為連續(xù)梯度比，φ(θ)函數(shù)的選用見表1。

表1 TVD格式中的常用通量限制器

本文主要對比帶有以上通量限制器的三種TVD格式：MUSCL、TVDLF和MTVDLF格式(見表2)，并通過數(shù)值實驗找出適應于不同類型問題的最優(yōu)格式。由于在對流擴散問題中，擴散項對不同的計算格式并不敏感，一般采用中心差分格式處理就可以滿足實際需要。因此，以下僅采用TVD格式構(gòu)造對流通量，三種格式中對流通量的構(gòu)造方法如下

(12)

表2 三種TVD格式

可以看出，當對流速度為常數(shù)時，MUSCL格式和MTVDLF格式本質(zhì)上是相同的。

4 數(shù)值實驗

針對以下五個典型的對流擴散問題，分別開展數(shù)值實驗。為了定量地比較這些數(shù)值格式的數(shù)值精度，進行了誤差分析，誤差定義為

(13)

4.1 線性對流擴散問題

(14)

(15)

(16)

將以上三式代入方程(1)，可得線性對流擴散方程的無量綱形式

(17)

4.1.1 階躍型純對流問題

該問題的初始濃度分布為

(18)

式中x0=0.1，為初始間斷位置，精確解為

(19)

模擬時間T=0.5，Δx=0.005，Δt=0.000005。分別采用MUSCL格式和TVDLF格式進行計算，數(shù)值結(jié)果如表3所示。

表3 算例1計算誤差

由上表可以看出，在網(wǎng)格數(shù)相同的情況下，選用同種通量限制器時，TVDLF格式的計算誤差遠大于MUSCL格式。且在兩種TVD格式中，Superbee限制器的誤差最小，而Minmod誤差最大。兩種格式的計算結(jié)果如圖1和圖2。

圖1 帶有Superbee和Minmod限制器的MUSCL格式計算結(jié)果

圖2 帶有Superbee和Woodward限制器的TVDLF格式

可以看出，基于TVD格式的數(shù)值解中沒有出現(xiàn)明顯的震蕩，只存在小部分激波被抹掉的現(xiàn)象，而基于TVDLF格式的計算結(jié)果出現(xiàn)了嚴重的數(shù)值耗散。因此，對于階躍型純對流問題，應優(yōu)先選用帶有Superbee限制器的MUSCL格式進行求解。

4.1.2 階躍型對流擴散問題

該問題的初始濃度分布及精確解為

c(x，0)=0， 0

(20)

(21)

該問題中模擬時間T=1，Δx=0.01，Δt=0.0001，Pe=100。分別采用MUSCL和TVDLF格式計算，發(fā)現(xiàn)TVDLF格式的計算結(jié)果仍會出現(xiàn)不可接受的數(shù)值耗散，故不再考慮該格式。MUSCL格式的計算結(jié)果見表4。

表4 算例2計算誤差

可以看出，基于MUSCL格式的計算結(jié)果與精確解高度吻合，計算誤差均在0.8%以下，該問題中限制器的類型基本不影響計算結(jié)果，Minmod限制器的計算結(jié)果如圖3。

圖3 帶有Minmod限制器的MUSCL格式計算結(jié)果

4.1.3 高斯型對流擴散問題

該問題中高斯分布的初始位置為x0=1/4，σ=1/30，Pe=200，初始濃度分布及精確解為：

(22)

(23)

該問題中模擬時間T=0.4，Δx=0.01，Δt=0.00001。分別采用MUSCL格式和TVDLF格式計算，TVDLF格式的計算結(jié)果仍會出現(xiàn)不可接受的數(shù)值耗散，故不再考慮。MUSCL格式的計算結(jié)果見表5。

表5 算例3計算誤差

可以看出，基于MUSCL格式的計算誤差介于1.2%～3.6%之間，并且 Woodward限制器的計算誤差最小，Minmod限制器的誤差最大，基于該限制器的計算結(jié)果如圖4。

圖4 帶有Woodward限制器的MUSCL格式計算結(jié)果

4.2 擬線性對流擴散問題

為了比選出TVD方法求解擬線性對流擴散問題時性能最優(yōu)的通量限制器，求解如下方程：

(24)

該問題中對流速度a=0.02，擴散系數(shù)Γ=0.0001，Pe=200，高斯分布的初始位置x0=1/4，σ=1/30，初始濃度分布及精確解為

(25)

(26)

模擬時間T=10s，Δx=0.01m，Δt=0.0005s。控制體界面兩側(cè)對應同一個a(1+x)，故MUSCL和MTVDLF格式相同?；赥VDLF格式的計算結(jié)果仍表現(xiàn)出了嚴重的耗散，故不再考慮。MUSCL格式的數(shù)值計算結(jié)果見表6。

表6 算例4計算誤差

可以看出，選用MUSCL格式計算時，Minmod限制器誤差最大，而Woodward限制器誤差最小。因此，對于擬線性高斯型對流擴散問題，應優(yōu)先選用帶有Woodward限制器的MUSCL或MTVDLF格式求解。選用Woodward和Minmod限制器的計算結(jié)果如圖5。

圖5 帶有Woodward和Minmod限制器的MUSCL格式計算結(jié)果

4.3 非線性對流擴散問題-Burgers 方程

Burgers方程是模擬沖擊波傳播和反射的非線性偏微分方程，為找出TVD方法求解該方程時的最優(yōu)限制器，設置了本算例。該方程的守恒形式為

(27)

(28)

該問題中計算域為[-1，2]，初始流速分布及精確解為

(29)

(30)

模擬時間T=2s，Δx=0.01m，Δt=0.001s，三種格式的數(shù)值計算結(jié)果見表7。

表7 算例5計算誤差

可以看出，MUSCL格式和MTVDLF格式的計算誤差均很小，而TVDLF格式的數(shù)值計算誤差較大，尤其是采用Albada和OSPRE限制器的情況。因此，對于Burgers方程，應優(yōu)先選用MUSCL或MTVDLF格式求解，且對于該問題限制器的類型基本不影響計算結(jié)果。帶有Superbee限制器的三種格式的計算結(jié)果如圖6。

圖6 帶有Superbee限制器的MUSCL、TVDLF、MTVDLF三種格式計算結(jié)果

5 結(jié)論

基于有限體積法，本文首先推導了求解對流擴散方程的通量格式，為提高計算效率，進一步對線性對流擴散方程的無量綱形式進行了推導，接下來將10種常用的通量限制器分別應用于3種典型的TVD格式進行了數(shù)值實驗，并通過比較以上格式的數(shù)值特性，得出了5類典型對流擴散問題對應的最優(yōu)通量限制器。