崔 冬,何小燕,楊海燕

(1.上海市水利工程設(shè)計研究院有限公司,上海 200061;2.上海圍海工程技術(shù)研究中心,上海 200061)

對流擴(kuò)散方程數(shù)值耗散的定量研究方法

崔 冬1,2,何小燕1,2,楊海燕1,2

(1.上海市水利工程設(shè)計研究院有限公司,上海 200061;2.上海圍海工程技術(shù)研究中心,上海 200061)

針對對流擴(kuò)散方程數(shù)值耗散問題普遍缺乏明確定量指標(biāo)的現(xiàn)狀,提出了度量數(shù)值耗散程度的定量指標(biāo)及其研究方法,即借助廣義“假擴(kuò)散”系數(shù)這一定量指標(biāo)來度量某一離散格式的數(shù)值耗散程度,其數(shù)值上近似等于具有相同擴(kuò)散效應(yīng)的真擴(kuò)散所對應(yīng)的擴(kuò)散系數(shù),其中真擴(kuò)散效應(yīng)可通過求解純擴(kuò)散問題的高精度數(shù)值解近似得到。采用該方法研究了三階QUICKEST離散格式的數(shù)值耗散規(guī)律,分析了時間步長、流速大小及空間步長對“假擴(kuò)散”系數(shù)的敏感性以及相關(guān)性,得到了“假擴(kuò)散”系數(shù)的具體表達(dá)式,進(jìn)而驗(yàn)證了該定量研究方法的可行性。

對流擴(kuò)散方程;數(shù)值耗散;QUICKEST離散格式;流體力學(xué)

保守性物質(zhì)的對流擴(kuò)散是流體力學(xué)、環(huán)境科學(xué)領(lǐng)域中經(jīng)常遇到的一類經(jīng)典問題,可用對流擴(kuò)散方程加以描述,如式(1)所示。該方程由非恒定項(xiàng)對流項(xiàng)u擴(kuò)散項(xiàng)組成,一般采用數(shù)值離散的方法求解,但求解過程中不可避免地會帶來一定的離散誤差而影響數(shù)值解的準(zhǔn)確性,其中影響較大的離散誤差是由一階導(dǎo)數(shù)項(xiàng)(非恒定項(xiàng)、對流項(xiàng))離散帶來的數(shù)值耗散[1]。數(shù)值耗散,也稱“假擴(kuò)散”,是截斷誤差的一種,會放大擴(kuò)散作用的效應(yīng),是對流擴(kuò)散方程數(shù)值求解中普遍存在的一個問題,且數(shù)值耗散程度越大,給數(shù)值解帶來的計算誤差越大。在應(yīng)用某種研究手段研究實(shí)際對流擴(kuò)散問題之前,有必要對研究手段的數(shù)值耗散程度有較為準(zhǔn)確的把握,以便判斷研究手段的適用性,并盡可能地減小這種“假擴(kuò)散”效應(yīng)[2]。

式中:c為物質(zhì)的質(zhì)量濃度;t為時間;x、y為直角坐標(biāo)系坐標(biāo);u、v分別為沿x、y方向的流速;Dx、Dy分別為物質(zhì)沿x、y方向的擴(kuò)散系數(shù)。當(dāng)u=v=0時,式(1)可簡化為純擴(kuò)散方程;當(dāng)Dx=Dy=0時,式(1)可簡化為純對流方程。

對流擴(kuò)散方程求解存在的數(shù)值耗散問題,一直受到研究人員的關(guān)注,并已取得了諸多研究成果[1-9],如提出了改進(jìn)離散格式、減小空間步長等減少“假擴(kuò)散”效應(yīng)的措施。但針對不同離散格式以及同一離散格式不同參數(shù)下的數(shù)值耗散程度尚缺乏統(tǒng)一的度量指標(biāo),已有的研究往往是通過比較經(jīng)典算例數(shù)值解與解析解的差異程度、純對流問題數(shù)值解最大濃度隨時間的降低程度等間接方式來判斷數(shù)值耗散程度的大小,這給實(shí)際應(yīng)用帶來不便。針對當(dāng)前研究現(xiàn)狀,本文提出度量數(shù)值耗散程度的定量指標(biāo)及研究方法,并將該方法應(yīng)用到三階QUICKEST離散格式的數(shù)值耗散規(guī)律的研究中,以進(jìn)一步驗(yàn)證方法的可行性。

1 定量指標(biāo)的提出

數(shù)值耗散最早定義為對流擴(kuò)散方程中一階導(dǎo)數(shù)項(xiàng)(非恒定項(xiàng)、對流項(xiàng))離散格式的截斷誤差小于二階而引起較大數(shù)值計算誤差的現(xiàn)象。以二維純對流方程為例,當(dāng)采用一階迎風(fēng)格式對非恒定項(xiàng)和對流項(xiàng)進(jìn)行數(shù)值離散時,略去三階及三階以上的項(xiàng),其等價表達(dá)式[10]為

式中:等號后兩項(xiàng)即為數(shù)值離散引入的數(shù)值耗散項(xiàng); Δx、Δy分別為沿x、y方向的空間步長;Δt為時間步長。將分別定義為沿x、y方向的“假擴(kuò)散”系數(shù),形式上起擴(kuò)散系數(shù)的作用。

事實(shí)上“假擴(kuò)散”效應(yīng)不僅存在于一階離散格式中,其他離散格式中也存在,如二階迎風(fēng)格式、三階迎風(fēng)格式、QUICK格式等[11]。參照式(2),設(shè)想其他離散格式下二維純對流方程也可寫成類似的等價表達(dá)式:

將Dfalse定義為廣義的“假擴(kuò)散”系數(shù),用來定量表征數(shù)值耗散程度,這樣就可將“假擴(kuò)散”系數(shù)這一重要參數(shù)的應(yīng)用范圍由一階迎風(fēng)離散格式推廣到更多離散格式中。

2 定量指標(biāo)的研究方法

2.1 研究方法的提出

由于數(shù)值離散的復(fù)雜性,很多離散格式難以直接寫出“假擴(kuò)散”系數(shù)的表達(dá)式?？紤]到“假擴(kuò)散”對數(shù)值解的影響表現(xiàn)為放大了擴(kuò)散效應(yīng),“假擴(kuò)散”應(yīng)與真擴(kuò)散具有類似的擴(kuò)散效應(yīng),若能找到與“假擴(kuò)散”效應(yīng)相當(dāng)?shù)恼鏀U(kuò)散,則可以認(rèn)為二者具有相同的擴(kuò)散系數(shù)。即“假擴(kuò)散”系數(shù)取值上近似等于具有相同擴(kuò)散效應(yīng)的真擴(kuò)散對應(yīng)的擴(kuò)散系數(shù),其中不同擴(kuò)散系數(shù)下的真擴(kuò)散效應(yīng),可以通過求解純擴(kuò)散方程的高精度數(shù)值解近似得到。

以下擬采用沿單一方向(x方向)的純對流問題來研究“假擴(kuò)散”效應(yīng),這樣既消除了物質(zhì)自身擴(kuò)散帶來的影響[6],又將二維問題簡化為一維,式(3)可簡化為式(4);相應(yīng)地,采用單一方向(x方向)的純擴(kuò)散問題來模擬真擴(kuò)散效應(yīng),如式(5)所示,可消除對流項(xiàng)離散帶來的計算誤差。

計算模型設(shè)置如下:計算區(qū)域?yàn)橐粋€長100km、寬8 km、水深10 m的河道,計算網(wǎng)格采用單一矩形網(wǎng)格,網(wǎng)格走向與河道走向一致。當(dāng)采用純對流問題研究“假擴(kuò)散”效應(yīng)時,初始時刻設(shè)置一個1 km× 1 km范圍、初始質(zhì)量濃度為20 mg/L的保守性物質(zhì)團(tuán)位于河道上游,其他計算區(qū)域初始濃度為零,水體沿河道縱向均勻流動,如圖1所示。當(dāng)采用純擴(kuò)散問題研究真擴(kuò)散效應(yīng)時,計算網(wǎng)格空間步長經(jīng)試算取25 m可達(dá)到足夠精度,水體靜止,初始時刻設(shè)置1 km×1 km范圍、初始質(zhì)量濃度為20 mg/L的保守性物質(zhì)團(tuán)位于河道中游,如圖2所示。

圖1 純對流問題下計算模型及初始濃度設(shè)置

圖2 純擴(kuò)散問題下計算模型及初始濃度設(shè)置

2.2 基本假定的驗(yàn)證

本文提出的“假擴(kuò)散”系數(shù)的研究方法有一個基本假定,即認(rèn)為當(dāng)真擴(kuò)散與“假擴(kuò)散”具有相同擴(kuò)散效應(yīng)時,二者應(yīng)當(dāng)具有相同的擴(kuò)散系數(shù);反過來說,當(dāng)二者具有相同擴(kuò)散系數(shù)時,二者應(yīng)當(dāng)具有相同的擴(kuò)散效應(yīng)。為了驗(yàn)證這一基本假定的正確性,以下采用“假擴(kuò)散”系數(shù)具有理論解的一階迎風(fēng)離散格式來模擬“假擴(kuò)散”效應(yīng),與具有相同擴(kuò)散系數(shù)的真擴(kuò)散效應(yīng)進(jìn)行對比,觀察二者是否相近。參數(shù)設(shè)置及一階迎風(fēng)離散格式的“假擴(kuò)散”系數(shù)理論解(根據(jù)式(2)計算得到)如表1所示。

表1 一階迎風(fēng)離散格式不同參數(shù)下的“假擴(kuò)散”系數(shù)理論解

采用具有一階迎風(fēng)格式的MIKE 21對流擴(kuò)散模型[12]求解如圖1所示的純對流問題,可分別得到2組參數(shù)下在時間t=6 h時的“假擴(kuò)散”效應(yīng)如圖3所示(截圖范圍以濃度團(tuán)為中心,沿河道縱向與橫向分別取10 km與4 km)。

圖3 一階迎風(fēng)離散格式下的“假擴(kuò)散”效應(yīng)(t=6 h,單位:mg/L)

真擴(kuò)散效應(yīng)通過求解圖2所示的純擴(kuò)散問題的高精度數(shù)值解近似得到,其中非恒定項(xiàng)采用一階迎風(fēng)格式離散,擴(kuò)散項(xiàng)采用二階中心差分格式離散,對大多數(shù)有實(shí)際意義的問題,二者的離散精度已能滿足需要[1]。表1所示的2組擴(kuò)散系數(shù)對應(yīng)的t=6 h時真擴(kuò)散效應(yīng)如圖4所示。對比圖3與圖4可以看出,“假擴(kuò)散”與真擴(kuò)散呈現(xiàn)出幾乎完全一致的擴(kuò)散效應(yīng)。這說明當(dāng)真擴(kuò)散與“假擴(kuò)散”具有相同擴(kuò)散系數(shù)時,二者具有相同擴(kuò)散效應(yīng)的基本假定是成立的;反過來說,當(dāng)真擴(kuò)散與“假擴(kuò)散”具有相同擴(kuò)散效應(yīng)時,二者具有相同擴(kuò)散系數(shù)的基本假定也是成立的。

圖4 純擴(kuò)散問題下的真擴(kuò)散效應(yīng)(t=6 h,單位:mg/L)

2.3 純擴(kuò)散問題的高精度數(shù)值解

通過求解圖2所示的純擴(kuò)散問題的高精度數(shù)值解,近似得到單一方向(x方向)擴(kuò)散系數(shù)取值為0～50 m2/s之間的濃度場分布情況,其中部分取值下的濃度場縱斷面分布如圖5所示。為方便實(shí)際應(yīng)用,本文根據(jù)圖5呈現(xiàn)的濃度場縱斷面分布特征,并參考以往研究成果,選用物質(zhì)團(tuán)在靜水中擴(kuò)散6 h與12h后對應(yīng)的擴(kuò)散范圍(以0.5mg/L包絡(luò)范圍來統(tǒng)計擴(kuò)散范圍)與最大濃度值來表征不同擴(kuò)散系數(shù)下的擴(kuò)散效應(yīng),如圖6所示。當(dāng)研究某種離散格式數(shù)值耗散程度時,可采用相應(yīng)的離散格式求解圖1所示的純對流問題,統(tǒng)計得到t=6 h與t=12 h的物質(zhì)團(tuán)擴(kuò)散范圍與最大濃度值,通過查圖6,找到相對應(yīng)的“假擴(kuò)散”系數(shù)。

圖5 部分?jǐn)U散系數(shù)下的濃度場縱斷面分布

圖6 擴(kuò)散系數(shù)與物質(zhì)團(tuán)擴(kuò)散范圍及最大濃度值的關(guān)系

此外,通過圖5、圖6還可以看出,擴(kuò)散作用的存在使保守性物質(zhì)團(tuán)最大濃度值降低程度明顯,當(dāng)擴(kuò)散系數(shù)分別取5 m2/s與30 m2/s時,經(jīng)過6 h單一方向的擴(kuò)散,物質(zhì)團(tuán)最大濃度分別由20 mg/L降低至14.6mg/L與6.95mg/L,降低幅度分別為27%與65%。

3 研究方法的初步應(yīng)用

以下以三階QUICKEST離散格式為例介紹度量數(shù)值耗散程度的定量方法在研究某一離散格式下數(shù)值耗散規(guī)律中的應(yīng)用,三階QUICKEST格式由廣泛應(yīng)用于一維數(shù)值模擬的二次迎風(fēng)插值格式——QUICK格式演變而來,三階QUICKEST格式在QUICK格式的基礎(chǔ)上加以改進(jìn),將應(yīng)用范圍擴(kuò)展至多維[3]。利用當(dāng)前節(jié)點(diǎn)值、上游最近的2個節(jié)點(diǎn)值、下游最近的1個節(jié)點(diǎn)值等10個節(jié)點(diǎn)共同構(gòu)造差分[3](式(6)),表達(dá)式復(fù)雜,難以直接確定“假擴(kuò)散”系數(shù)。

式中:n為時間步長的編號;j為節(jié)點(diǎn)沿x方向的編號;k為節(jié)點(diǎn)沿y方向的編號;A1～A10為系數(shù)。

由式(2)可知,“假擴(kuò)散”系數(shù)一般與時間步長、流速大小、空間步長有關(guān)。為了得到三階QUICKEST離散格式的數(shù)值耗散規(guī)律,分別研究了時間步長、流速大小、空間步長對“假擴(kuò)散”系數(shù)的敏感性以及相關(guān)性。

3.1 時間步長的敏感性及相關(guān)性

采用三階QUICKEST離散格式分別研究了Δx =100 m、200 m兩組空間步長(流速以0.5 m/s為例)下的“假擴(kuò)散”效應(yīng),通過查圖6可以得到相對應(yīng)的“假擴(kuò)散”系數(shù),如圖7所示。可以看出,時間步長對“假擴(kuò)散”系數(shù)不敏感,改變時間步長,“假擴(kuò)散”系數(shù)幾乎不受影響(理論上改變時間步長會對舍入誤差略有影響,因而造成計算結(jié)果之間的微小差異)。

圖7 三階QUICKEST離散格式下時間步長對“假擴(kuò)散”系數(shù)的敏感性

3.2 流速大小的敏感性及相關(guān)性

采用三階QUICKEST離散格式分別研究了100m、200m兩組空間步長下流速大小取0.1～2 m/s時對應(yīng)的“假擴(kuò)散”效應(yīng),通過查圖6可以得到相對應(yīng)的“假擴(kuò)散”系數(shù)如圖8所示?？梢钥闯?流速大小對“假擴(kuò)散”系數(shù)敏感,改變流速大小,“假擴(kuò)散”系數(shù)會隨之變化。采用趨勢線擬合發(fā)現(xiàn),三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)與流速的0.5次方呈相當(dāng)好的擬合關(guān)系。

圖8 三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)與流速大小的相關(guān)關(guān)系

3.3 空間步長的敏感性及相關(guān)性

采用三階QUICKEST離散格式分別研究了0.5 m/s、2m/s兩組流速下空間步長取25～500m時的“假擴(kuò)散”效應(yīng),通過查圖6可以得到相對應(yīng)的“假擴(kuò)散”系數(shù)如圖9所示?？梢钥闯?空間步長對“假擴(kuò)散”系數(shù)相當(dāng)敏感,改變空間步長,“假擴(kuò)散”系數(shù)會發(fā)生較大變化。采用趨勢線擬合發(fā)現(xiàn),三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)與空間步長的2次方呈相當(dāng)好的擬合關(guān)系。

圖9 三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)與空間步長的相關(guān)關(guān)系

3.4 “假擴(kuò)散”系數(shù)的表達(dá)式

由圖8、圖9可以看出,三階QUICKEST離散格式下的“假擴(kuò)散”系數(shù)與流速的0.5次方及空間步長的2次方呈正相關(guān)關(guān)系。因此,認(rèn)為三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)可寫成

式中η為系數(shù),可由圖8及圖9的擬合關(guān)系推求得到。以圖8為例,當(dāng)空間步長為100m時,由Dfalse(x)= ηu1/2Δx2=1002u1/2η=0.75u1/2,可 推 求 出η=7.5×10-5。由圖9亦可得到η=7.5×10-5。因此三階QUICKEST離散格式下“假擴(kuò)散”系數(shù)表達(dá)式可寫為

4 結(jié) 論

a.可借助廣義“假擴(kuò)散”系數(shù)這一定量指標(biāo)來度量某一離散格式的數(shù)值耗散程度,“假擴(kuò)散”系數(shù)數(shù)值上近似等于具有相同擴(kuò)散效應(yīng)的真擴(kuò)散所對應(yīng)的擴(kuò)散系數(shù),其中真擴(kuò)散效應(yīng)可通過求解純擴(kuò)散問題的高精度數(shù)值解近似得到。

b.采用該定量研究方法,成功得到了三階QUICKEST離散格式的數(shù)值耗散規(guī)律,驗(yàn)證了該定量研究方法的可行性。

[1]陶文銓.計算傳熱學(xué)的近代進(jìn)展[M].北京:科學(xué)出版社,2000:105-154.

[2]楊屹松,李煒.混合有限分析格式的數(shù)值耗散及頻散[J].水利學(xué)報,1992,22(10):26-31.(YANG Yisong,LI Wei.Numerical dissipation and dispersion of hybrid finite analytic method[J].Journal of Hydraulic Engineering, 1992,22(10):26-31.(in Chinese))

[3]EKEBJARG L,JUSTSEN P.An explicit scheme for advection-diffusion modeling in two dimensions[J]. ComputerMethodsinAppliedMechanicsand Engineering,1991,88(3):287-297.

[4]陳炳睿,朱建榮.物質(zhì)輸運(yùn)方程中平流項(xiàng)數(shù)值格式的改進(jìn)[J].華東師范大學(xué)學(xué)報:自然科學(xué)版,2006(4):64-70.(CHEN Bingrui,ZHU Jianrong.Numerical scheme improvement of the advection term in mass transport equation[J].Journal of East China Normal University: Natural Science,2006(4):64-70.(in Chinese))

[5]謝志華,林建國,由曉丹.污染物對流擴(kuò)散方程的幾種新的高階QUICK組合顯格式比較研究[J].水動力學(xué)研究與進(jìn)展:A輯,2005,20(3):346-356.(XIE Zhihua, LIN Jianguo,YOU Xiaodan.Comparative investigation of some high-order explicit schemes combined with QUICK for the convection-diffusion equation of pollutants[J]. Chinese Journal of Hydrodynamics,2005,20(3):346-356.(in Chinese))

[6]程愛杰,趙衛(wèi)東.對流擴(kuò)散方程的經(jīng)濟(jì)差分格式[J].計算數(shù)學(xué),2000,22(3):309-318.(CHENG Aijie,ZHAO Weidong.An economical difference scheme for convectiondiffusion equations[J].Mathematica Numerica Sinica, 2000,22(3):309-318.(in Chinese))

[7]林建國,殷佩海.污染濃度對流擴(kuò)散的三節(jié)點(diǎn)高精度格式[J].大連海事大學(xué)學(xué)報,2002,28(4):64-67.(LIN Jianguo,YINPeihai.Three-pointschemewithhigh accuracy for the convection-diffusion equation of the pollutant concentration[J].Journal of Dalian Maritime University,2002,28(4):64-67.(in Chinese))

[8]由同順.非線性對流-擴(kuò)散方程的高階特征:差分格式及其誤差估計[J].數(shù)值計算與計算機(jī)應(yīng)用,1994(4): 312-317.(YOUTongshun.High-ordercharacteristicdifference schemes for the nonlinear convection:diffusion equation andtheirerrorestimates[J].Journalon Numerical Methods and Computer Applications,1994(4): 312-317.(in Chinese))

[9]盧洪波,謝鴻罡,楊沫,等.一種改進(jìn)的QUICK差分方案[J].東北電力學(xué)院學(xué)報,1998,18(4):1-6.(LU Hongbo,XIE Honggang,YANG Mo,et al.A modified numerical method with QUICK scheme[J].Journal of Northeast China Institute of Electric Power Engineering, 1998,18(4):1-6.(in Chinese))

[10]汪德爟.計算水力學(xué)理論與應(yīng)用[M].南京:河海大學(xué)出版社,1989:117-118.

[11]胡洪,黃虎,宋倩倩,等.不同對項(xiàng)離散格式對流場計算的影響[J].數(shù)值計算與計算機(jī)應(yīng)用,2010,31(2):153-160.(HU Hong,HUANG Hu,SONG Qianqian,et al.The influence of different discrete scheme of the convective term on term on fluid field calculation[J].Journal on Numerical Methods and Computer Applications,2010,31 (2):153-160.(in Chinese))

[12]許婷.丹麥MIKE21模型概述及應(yīng)用實(shí)例[J].水利科技與經(jīng)濟(jì),2010,16(4):867-869.(XU Ting.Calculation principle and application example of a two-dimensional flow model:MIKE21 HD[J].Water Conservancy Science and Technology and Economy,2010,16(4):867-869.(in Chinese))

Quantitative research method of numerical dissipation levels of convective-diffusive equation//

CUI Dong1,2,HE Xiaoyan1,2,YANG Haiyan1,2
(1.Shanghai Water Engineering Design and Research Institute Co.,Ltd.,Shanghai 200061, China;2.Shanghai Engineering Research Center of Reclamation,Shanghai 200061,China)

Considering the fact that quantitative indicator on numerical dissipation is lacking,in this paper the quantitative indicator and research method were given.The method is,namely,using generalized“false diffusive”coefficient of the certain indicators to measure the degree of numerical dissipation of a discrete format.The numerical value of the above method approximately equals to the same effect on the true diffusion of the diffusion coefficient.As to the true diffusion effect,the highly precise numerical value can be approximately obtained by the solution of the pure diffusion formulas.This method was successfully applied to third-order QUICKEST discrete scheme.The overall results show that false diffusion coefficient could be introduced to characterize the numerical dissipation level,whose value can be obtained by comparing the true diffusion effect.

convective-diffusive;numerical dissipation;QUICKEST discrete scheme;fluid mechanics

TV131.2;O241.82

:A

:1006-7647(2014)05-0008-04

10.3880/j.issn.1006-7647.2014.05.002

201309-18 編輯:駱超)

上海市科學(xué)技術(shù)委員會資助課題(13DZ2251500)

崔冬(1982—),女,江蘇徐州人,工程師,碩士,主要從事河口海岸工程研究。E-mail:njcuidong.student@sina.com