有限容積法對流項離散格式的穩定性與有界性研究

宇波,焦開拓,陳宇杰,李敬法,鄧雅軍,王鵬,孫東亮

(1. 北京石油化工學院機械工程學院,102617,北京;2. 西安交通大學動力工程多相流國家重點實驗室,710049,西安)

有限容積法是數值求解工程流動傳熱問題的主流方法。基于有限容積法進行控制方程離散時,若對流項處理不當,可能導致模擬結果出現無物理意義的振蕩或越界現象。因此,對流項離散格式的性質成為了有限容積法研究的重點[1-2]。對流項離散常用的低階格式有一階迎風格式(first-order upwind difference, FUD)、混合格式、指數格式和乘方格式[3]。對于強對流問題,采用低階格式計算穩定,但存在假擴散,計算精度較低。二階精度及以上的高階格式可以顯著減輕假擴散現象[4],常見的高階格式有二階迎風格式(second-order upwind difference, SUD)、中心差分格式(central difference, CD)和QUICK格式等。然而,采用高階格式可能由于不滿足對流項離散格式的穩定性得到非物理意義的振蕩解[5],對流項穩定性分析是有限容積法高階格式研究的關鍵點[6]。SUD格式只涉及到上游點,為絕對穩定格式,而CD格式、QUICK格式涉及下游點,在網格貝克萊數較大時會產生數值振蕩,即為條件穩定格式[7]。此外,組合SUD、CD格式的SCSD(stability-controllable second-order difference)格式[8],通過調節SUD、CD格式的權重實現穩定性可控。基于SCSD格式,Li等[9]提出了一種通過網格貝克萊數大小調節SUD、CD格式權重的SGSD(stability-guaranteed second-order difference)格式,該格式可根據當地網格貝克萊數調整權重系數,使其自動滿足對流項穩定性條件,兼顧了計算精度和對流項的穩定性。權重系數的計算方法并不唯一,不同權重計算方法對計算精度性存在影響,但對此鮮有研究。此外,不同對流項離散格式的穩定性條件由當地網格貝克萊數確定,對于與邊界相鄰的節點,無論是采用內節點法還是外節點法,對于均分網格由于界面與節點之間位置關系特殊,可能導致穩定性條件與內部節點存在差異。但是,目前對流項離散格式的穩定性分析均基于計算區域內部節點,推導得出的對流穩定性條件對于與邊界相鄰節點的適用性有待探討。

絕對穩定的對流項離散格式可以避免產生非物理意義的振蕩解,但部分絕對穩定的離散格式并不能避免越界現象(待求變量超出物理問題本身所規定的物理量的上下限),即無法保證有界性。對流項離散格式的有界性也是有限容積法研究的關鍵點[3]。為了克服越界現象,已發展了多類方法,例如通量密度修正法[10-12]、規正變量圖(normalized variable diagram,NVD)[13-14]和總變差減小法(total variation diminishing,TVD)[15-17]等。基于規正變量圖,Gaskell等[18]提出了對流有界性CBC(convection boundedness criterion)準則,滿足該準則的高階格式被稱為高階有界格式。經過幾十年的發展,國內外學者們提出了多種高階有界格式,常見的有SMART[18]、MINMOD[19]、HOAB[20]、SUPERBEE[21]、MUSCL[22]等。這些格式在計算精度、收斂速率和健壯性上表現各異[23],如何構造一種平衡各類性能的高階有界格式值得探究。

對于上述有限容積法對流項離散格式在穩定性和有界性方面存在的不足,本文采用理論分析和數值試驗相結合的方法開展研究:對于對流項穩定性,探究與邊界相鄰節點的穩定性條件和SGSD格式的權重系數計算方法;對于對流項有界性,分析常見的高階有界格式特性,總結有利于編程實施的通用表達形式,提出性能可控的高階有界格式,并進一步優化得到綜合計算性能較優的高階有界格式。

1 對流項離散格式的穩定性

若采用某一對流項離散格式得到的收斂解是不具有物理意義的振蕩解,則稱該對流項離散格式是不穩定的,否則是穩定的。本節主要針對與邊界相鄰節點的對流穩定性條件以及絕對穩定的SGSD格式權重系數計算方法進行說明。

1.1 與邊界相鄰節點的對流穩定性條件

上述穩定性條件僅適用于與邊界不相鄰的內點。對于與邊界相鄰的節點,無論是采用內節點法還是外節點法,邊界面上對流通量(CD格式、QUICK格式等)和擴散通量的表達式與其他內節點不同,由此得到的邊界相鄰節點的對流穩定性條件與其他內點不同。內節點法西邊界面網格示意圖如圖1所示。

圖1 內節點法西邊界面網格示意圖Fig.1 Schematic of the grid at the west boundary based on the cell centered scheme

(1)

式中:φ為待求變量;下標e、w代表東、西界面。

若對流項采用QUICK格式,擴散項一階導數采用界面相鄰兩點近似,采用圖1所示的均分網格,則式(1)可離散為

(2)

整理得

(3)

由正型系數法[3]可知,對流穩定性應滿足

(4)

式(2)中擴散項在邊界面處一階導數采用一階精度格式離散,而對流項中各界面處待求變量和東界面處一階導數采用二階精度格式離散。為使得各界面對流通量和擴散通量離散階數相等,邊界面處一階導數需采用三點近似的二階格式,此時離散表達式為

(5)

整理得

(6)

由正型系數法可知,對流穩定性應滿足

(7)

同理,推導了其他情況與邊界相鄰節點采用CD格式、QUICK格式時的對流項穩定性條件,如表1所示,表中上標1表示邊界處擴散通量中的一階導數采用一階精度,2表示采用二階精度。由此可知,與邊界相鄰節點的對流穩定性條件不僅與離散格式相關,還與區域離散方式(內、外節點法)、界面流速方向和邊界法線方向相關。

(a)內節點法

(b)外節點法

對流項穩定性條件是在理想和苛刻的條件下得到的,事實上對流離散格式的穩定性條件受到流動的維度、非線性、內熱源和網格類型等多種因素的影響,在流動與傳熱實際工程問題中使數值解發生振蕩的PeΔ值往往比理想條件下得到的結果大得多。

1.2 絕對穩定的SGSD格式權重計算

數值模擬計算中,為減輕假擴散的影響,對流項推薦采用高階格式進行離散。從物理意義上講,大PeΔ表示離散的對流項影響較強、擴散項影響較弱,采用絕對穩定的SUD格式離散較為合適;反之,小PeΔ表示離散的對流項影響較弱、擴散項影響較強,采用條件穩定的CD格式較為合適。據此,結合SUD格式、CD格式的優點,文獻[9]提出對強對流、弱對流問題均適用的絕對穩定的SGSD格式

(8)

為了減小假擴散并提高計算精度,王賢鋼等[24]提出了一種新的權重系數計算方法

(9)

式中:(PeΔ)max為計算區域中PeΔ的最大值。相較于原權重系數計算方法,該方法在PeΔ>6時增大了CD格式的權重,但修正后的權重系數計算方法并不滿足絕對穩定性條件。

CD格式、SUD格式的截差首項系數不同,當滿足對流項穩定性時CD格式、SUD格式造成的計算誤差也不同,因此SUD格式、CD格式的權重配比會影響SGSD格式整體的計算誤差。本文以滿足絕對穩定性條件為前提,提出了兩種權重系數計算方法:

方法Ⅰ

(10)

方法Ⅱ

(11)

采用SGSD格式不同權重系數計算方法得到的平均誤差隨PeΔ的變化如圖3所示。由圖3可知:當PeΔ為1～2時,按誤差從大到小排序為方法Ⅱ、原權重系數計算方法和文獻[24]方法、方法Ⅰ;當PeΔ為3～7時,排序為方法Ⅰ、方法Ⅱ、原權重系數計算方法和文獻[24]方法;當PeΔ為7～9時,排序為方法Ⅰ、原權重系數計算方法和文獻[24]方法、方法Ⅱ;當PeΔ>9時,排序為文獻[24]方法、原權重系數計算方法、方法Ⅰ、方法Ⅱ。綜合來看,SGSD格式原權重系數計算方法在較大PeΔ取值范圍內并不是最優的,可根據問題需要進行選擇。強對流問題可采用本文所提權重計算方法Ⅰ、方法Ⅱ。

(a)S=0

(b)S=800x

2 對流項離散格式的有界性

本節主要總結常見對流項高階有界格式的特點,提出通用表達式和兼顧計算精度、收斂速率、健壯性的高階有界格式。

2.1 CBC準則

關于對流項有界性的討論常基于規正變量的形式,規正變量和規正空間坐標的定義為

(12)

圖4 節點坐標及下標U、C、D和f的相對位置關系Fig.4 Node coordinate and relative position of subscripts U, C, D, f

對流項滿足有界性的充分必要條件為在規正變量圖上滿足CBC準則,具體表達式為

(13)

式(13)所示在規正變量圖中對應的區域如圖5中陰影部分所示。在規正變量圖中滿足CBC準則的格式稱為有界格式。

圖5 對流項有界性的CBC準則示意圖 Fig.5 Schematic of boundedness criteria for convection term

常見對流項離散格式在均分網格下規正變量圖中的圖線如圖6所示。格式圖線與縱軸的交點等于該格式對流穩定性條件的倒數,與原點相交的SUD、FUD格式為絕對穩定格式。圖6中位于CD、SUD格式之間的陰影區域可視作CD、SUD兩種二階格式的組合,因此具有二階截差精度。具有二階及以上截差精度并同時滿足CBC準則的格式稱為高階有界格式,可以證明高階有界格式是絕對穩定的格式。基于高階有界格式的這些優點,建議在工程實際問題中應盡可能采用高階有界格式離散對流項,以保證得到精度較高且具有物理意義的數值解。

圖6 均分網格下常見對流項離散格式的NVD圖Fig.6 NVD plot of typical discretized schemes for convection term in the uniform grid system

2.2 高階有界格式通用表達形式

(14)

表2 常見高階有界格式坐標點

圖7 通用高階有界格式示意圖Fig.7 Schematic of general high-order bounded scheme

圖8 性能可控高階有界格式示意圖Fig.8 Schematic of high-order bounded scheme with controllable performance

(15)

AHB格式的非均分網格表達式為

(16)

Nerr=max(|φn+1-φn|)

(17)

式中φn+1、φn為相鄰兩迭代步的待求變量向量。設置迭代收斂標準為Nerr≤1×10-8,最大外迭代次數為5 000。

圖9 二維單階梯突變純對流問題示意圖 Fig.9 Schematic of the two-dimensional convective problem with a sharp interface

為保證不同C、D取值情況下,絕大部分區域均可獲得收斂解,松弛因子為0.1。平均誤差、計算時間隨C、D取值的變化如圖10所示,其中平均誤差由x=0.5處中垂線待求變量結果獲得。

(a)平均誤差

(b)計算時間分布

(18)

2.3 案例分析

2.3.1 單階梯突變的純對流問題

該問題計算參數與2.2小節單階梯突變的純對流問題相同。松弛因子為0.2時各高階有界格式在x=0.5處中垂線的變量分布如圖11所示。由圖11可知,MINMOD格式假擴散最嚴重,待求變量在階梯處平緩變化與解析解差距最大,SMART、MUSCL格式與解析解的差距次之,SUPERBEE、SAHB、HOAB和ABC這4種格式均可較好地捕捉到階梯突變結果。表3給出了松弛因子為0.8、0.6、0.4、0.2時各格式平均誤差、迭代次數和計算時間,迭代次數為5 000表示該格式并未完全收斂。由表3可知,各格式之間計算精度、收斂速率和健壯性差別較大,且沒有一種格式能在多個方面明顯優于其他格式。單階梯突變純對流問題采用不同高階有界格式的平均誤差如圖12所示,單階梯突變純對流問題采用不同高階有界格式的最大松弛因子如圖13所示,可知SUPERBEE、SAHB、HOAB和ABC格式計算精度相對較高,但這些格式達到收斂標準對松弛因子的要求較高,尤其是SUPERBEE格式,在測試中并未收斂。在這4種計算精度相對較高的格式中,ABC格式可獲得的最大松弛因子為0.45,小于HOAB格式,但前者收斂過程更加穩定,所需的收斂時間明

圖11 單階梯突變純對流問題中垂線計算結果 Fig.11 Variable profile at the mid-vertical line for the convective problem with a sharp interface

圖12 單階梯突變純對流問題采用不同高階有界格式的平均誤差Fig.12 Mean error of high-order bounded schemes for the convective problem with a sharp interface

顯小于后者,并且計算效率接近可取高松弛因子的MINMOD、MUSCL格式,單階梯突變純對流問題采用不同高階有界格式的計算時間如圖14所示。

圖13 單階梯突變純對流問題采用不同高階有界格式的最大松弛因子Fig.13 Maximum relaxing factor of high-order bounded schemes for the convective problem with a sharp interface

因此,ABC格式具有良好的收斂速率和健壯性,且計算精度較高。計算時間為不同松弛因子下計算收斂所需的最短時間,后文如無特殊說明,均與此相同。

圖14 單階梯突變純對流問題采用不同高階有界格式的計算時間Fig.14 Computation time of high-order bounded schemes for the convective problem with a sharp interface

表3 單階梯突變純對流問題中各高階有界格式計算性能

2.3.2 Smith-Hutton純對流問題

第一個案例基于均勻流場進行分析,為對比不均勻流場中不同高階有界格式的計算性能,采用Smith-Hutton純對流問題進行研究。如圖15所示,幾何區域大小為2 m×1 m,區域中存在旋轉流場,通過邊界條件控制:下邊界被分成左、右兩部分,左半邊入口為第一類邊界條件,變量分布滿足雙曲正切函數φ=1+tanh[α(1+2x)],其中α取40;右半邊出口為第二類邊界條件,界面法向導數為0。該問題的解析解為φ=1+tanh[α(1-2x)]。采用均分網格內節點法對計算區域進行離散,網格數為120×60,求解方法與案例1相同。

圖15 Smith-Hutton純對流問題示意圖 Fig.15 Schematic of the Smith-Hutton convection problem

松弛因子為0.2時各高階有界格式在右半部分出口處的變量分布如圖16所示。由于其解析解分布為雙曲正切函數,因此變量分布并非完全的階梯突變。由圖16可知,該問題中各格式計算精度差異與前兩個案例大致相同,SUPERBEE、SAHB、HOAB和ABC格式的計算結果更接近解析解。表4和圖17～19給出了該問題中各格式的計算精度、收斂速率和健壯性對比結果。不同于第一個算例,ABC的計算精度、收斂速率、健壯性均優于HOAB、MUSCL格式。MUSCL、SUPERBEE和SAHB格式未完全收斂。

圖16 Smith-Hutton純對流問題出口變量分布Fig.16 Variable profile of the Smith-Hutton convection problem at the outlet

表4 Smith-Hutton純對流問題中各高階有界格式計算性能

圖17 Smith-Hutton純對流問題采用不同高階有界格式的平均誤差Fig.17 Mean error of high-order bounded schemes for the Smith-Hutton convection problem

圖18 Smith-Hutton純對流問題采用不同有界格式的最大松弛因子 Fig.18 Maximum relaxing factor of high-order bounded schemes for the Smith-Hutton convection problem

圖19 Smith-Hutton純對流問題采用不同高階有界格式的計算時間Fig.19 Computation time of high-order bounded schemes for the Smith-Hutton convection problem

2.3.3 頂蓋驅動問題

上述兩個案例均為二維純對流問題,為了分析本文所提出的高階有界格式在三維對流擴散問題以及非線性流動問題中的計算性能,本小節采用三維方腔頂蓋驅動問題對此進行研究。雷諾數為1 000,所有壁面均為無滑移,采用50×50×50的結構化均分網格對計算區域進行離散。采用SIMPLE算法對該問題進行求解,收斂標準為動量方程、連續性方程余量的二范數分別小于6.0×10-6、1.0×10-10,對于收斂困難的SUPERBEE格式,動量方程收斂標準采用1.0×10-4。

本文采用文獻[25]中z=0.5平面與x=0.5平面相交線上的u方向速度的計算結果作為基準解,進一步計算了相同網格節點位置處、不同格式計算結果與基準解之間的相對差異,如圖20所示。對比結果表明,ABC格式平均相對偏差最小,HOAB格式、SAHB格式次之。不同于上述給定流場問題,SUPERBEE格式無法收斂且計算精度最低。

圖20 三維頂蓋驅動問題采用不同高階有界格式的相對偏差Fig.20 Relative error of different discrete schemes of convection term for the three-dimensional lid-driven flow problem

計算時間對比結果如圖21所示。由圖21可以看出,ABC與MINMOD、MUSCL格式的收斂所需時間均在1 000 s左右,且明顯低于其他高階有界格式。此外,采用經典的QUICK格式得到的計算時間為922 s,與ABC格式差異不大。

圖21 三維頂蓋驅動問題采用不同對流項離散格式的計算時間Fig.21 Computation time of different discrete schemes of convection term for the three-dimensional lid-driven flow problem

采用不同對流項離散格式的最大松弛因子如圖22所示。由圖22可以看出,采用ABC、MINMOD和MUSCL格式時,亞松弛因子均可取到0.85,且大于其他高階有界格式。此外,采用QUICK格式的最大松弛因子同為0.85。這表明ABC格式具有較高的求解健壯性。

圖22 三維頂蓋驅動問題采用不同對流項離散格式的最大松弛因子Fig.22 Maximum relaxing factor of different discrete schemes of convection term for the three-dimensional lid-driven flow problem

綜上,對于階梯突變、Smith-Hutton純對流和頂蓋驅動問題,所測試的現有格式均無法兼顧計算精度、收斂速率和健壯性,且這些格式在不同問題下的表現不同。為了進一步對比說明ABC格式在計算精度、收斂速率和健壯性方面的綜合表現,以上述3個算例為對象,對不同格式的計算精度、計算時間和最大松弛因子進行分別排名,并以排名大小為分數進行加權平均得到總分,總分越小意味著綜合表現越佳。為了提高對比的可靠性,分別采用兩種不同的排名方式和權重因子。兩種排名方式為:①考慮計算精度和計算時間因編程風格和計算機運行狀態不同而導致的差異,當兩種格式差異較小時,認為兩者一致,并列排名;②不考慮編程風格和計算機運行狀態差異帶來的偏差,直接根據計算精度和計算時間的大小進行排名。兩種權重因子為:①計算精度、計算時間、最大松弛因子的權重取值分別為0.5、0.25、0.25;②計算精度、計算時間、最大松弛因子的權重取值均為1/3。不同排名方式下不同高階有界格式在3種經典算例中的計算性能排名如表5～8所示,可知ABC格式雖然在不同問題中的計算精度、收斂速率和健壯性的排名不同,但均整體表現最優,具有求解精度高、收斂快、健壯性好的特點。

表5 各高階有界格式計算性能對比(排名方式①,權重系數①)

表6 各高階有界格式計算性能對比(排名方式①,權重系數②)

表7 各高階有界格式計算性能對比(排名方式②,權重系數①)

表8 各高階有界格式計算性能對比(排名方式②,權重系數②)

3 結 論

本文對有限容積法對流項離散格式的穩定性和有界性進行了較為深入的研究,主要結論如下。

(1)通過對流項離散格式的穩定性分析,發現內、外節點法邊界相鄰節點的對流項穩定性條件與內點的存在較大差異;對于SGSD格式,存在多種權重系數計算方法。

(2)根據常見高階有界格式線型為分段線性函數的特點,推導得到通用的分段線型高階有界格式GHB,并給出了均分、非均分網格下常見格式的折點取值。

(3)提出了通過調節C、D兩個控制因子實現性能可控的高階有界格式AHB格式,并通過參數優化得到了綜合性能較優的ABC格式。采用不同的高階有界格式計算了3個經典算例,并對計算誤差、計算時間和最大松弛因子進行綜合評估,結果表明在不同評估方法下,ABC格式在計算精度、收斂速率和健壯性的整體性能方面均表現最優。