節塊展開法求解對流擴散方程的穩定性和數值耗散特性分析

周夏峰，李 富

（清華大學 核能與新能源技術研究院 先進反應堆工程與安全教育部重點實驗室，北京 100084）

節塊展開法求解對流擴散方程的穩定性和數值耗散特性分析

周夏峰，李 富

（清華大學 核能與新能源技術研究院 先進反應堆工程與安全教育部重點實驗室，北京 100084）

本文研究了節塊展開法求解對流擴散方程的穩定性和數值耗散特性。通過離散方程精確解和數值實驗方法分析不同階節塊展開法的穩定性和數值耗散特性，并將其與有限體積法中的中心差分和一階迎風格式進行對比。結果表明：偶數階節塊展開法的穩定是有條件的，即Peclet數（Pe）小于限值，且Pe限值會隨展開階數的增大而增大，其穩定性范圍和精度均優于中心差分格式；奇數階節塊展開法是無條件穩定的，但隨Pe的增大，數值耗散增大、計算誤差增大，且當Pe大于一定值后，產生的數值耗散大于一階迎風格式。

節塊展開法；對流擴散方程；穩定性；數值耗散

節塊法由于其計算速度快、計算精度高等優點，已廣泛用于反應堆物理計算中［1］，典型節塊法有節塊展開法（NEM）［1］、節塊格林函數法［2］、解析節塊法［3］（包括節塊積分法）及半解析節塊法［4］等，NEM已被應用于高溫氣冷堆物理問題的高效求解［5］。目前，節塊法應用到熱工計算的研究相對較少，主要以節塊積分法為主［6－8］；而針對節塊展開法求解熱工問題開展的研究則更少。文獻［9］對節塊展開法進行了初步改進，針對4階基函數展開的節塊展開法的穩定性和誤差進行了分析，并將改進后的節塊展開法推廣到圓柱幾何下對流擴散方程的求解［10］；然而針對其他階基函數展開的節塊展開法的穩定性目前還未進行具體分析，且各階節塊展開法的數值耗散特性目前亦未進行分析。

對流擴散方程數值計算結果誤差的主要來源在于對流項的離散格式。當對流項的離散格式不合理時，會人為地增大或削弱擴散項的作用，使計算結果產生較大誤差，即數值耗散特性；在一定條件下（如流速較高或網格劃分較粗時）甚至會導致數值解發生振蕩，產生不穩定［11］。在計算流體力學和計算傳熱學中離散格式的穩定性和數值耗散特性是一基本問題。因此，本文研究不同階基函數展開的節塊展開法求解對流擴散方程的穩定性和數值耗散特性，并與有限體積法中的中心差分和一階迎風格式進行對比和分析。

1 基本方程和求解過程

采用求解中子擴散方程的類似思路［1］，節塊展開法求解對流擴散方程也可利用橫向積分法將m維對流擴散方程轉化為m個一維的對流擴散方程（即分方向的橫向積分方程），各一維對流擴散通過橫向泄漏項相互耦合。因此，研究一維節塊展開法離散形式的穩定性和數值耗散特性具有重要意義。此外，在進行離散格式的特性分析時，大家普遍針對一維穩態模型方程進行研究［9］，這是由于采用的模型相對簡單，可更方便地突出所研究的離散格式存在的問題。因此，基于以上原因，本文也選取一維穩態的對流擴散方程作為研究分析起點。直角坐標系下一維穩態對流擴散方程可寫為：

其中：ρ（x）為流體密度；u（x）為x方向的速度函數；Q（x）為源項。

選擇不同的φ（x）、Г（x）和Q（x），可代表不同的物理方程：當φ（x）代表x方向的速度函數，Г（x）代表動力黏度，Q（x）代表壓力梯度項時，式（1）代表Navier-Stokes方程；當φ（x）代表溫度，Г（x）代表導熱系數時，式（1）代表能量守恒方程。用NEM求解該方程時，計算區域被分為K個節塊，節塊k的尺寸為［－Δxk，Δxk］，k＝1，…，K。局部坐標原點在節塊中心位置，Δxk為節塊k在x方向的半寬度，如圖1所示。

圖1 節塊劃分示意圖Fig.1 Nodal discretization for NEM

為便于計算和分析，每個節塊內ρ（x）、u（x）和Г（x）分別取節塊平均值ρk、uk和Гk，將φ（x）在節塊k內用N階勒讓德正交基函數展開，Q（x）用M階勒讓德正交基函數展開，則：

上述勒讓德正交基函數滿足以下遞推和正交關系式：

Q（x）的展開系數已知，而φ（x）的展開系數未知，即每個節塊有N＋1個未知變量，所以需找到N＋1個約束條件。

1）相鄰兩節塊邊界處滿足φ（x）連續性

2）相鄰兩節塊邊界處滿足流連續

利用Γk連續性條件，可得到以下關系式：

3）矩方程條件

由于每個節塊需N＋1個約束條件，而上述僅有兩個約束條件，這里采用剩余權重法得到0～N－2矩方程來補充剩余所需的N－1個約束條件。

根據式（5）、（6）、（9）得到的N＋1個等式，可得到φ（x）展開系數的表達式為：

其中：Ek、Fk和Gk分別為的離散系數；Hk為源項經過上述處理形成的值。

2 穩定性分析

本文通過給出離散方程精確解，來研究節塊展開法求解對流擴散方程穩定性問題［12］。即在一定條件下，通過一定的方法，求出離散方程對應的精確解，根據離散方程精確解的形式，來分析其穩定性所滿足的條件。此外，對于離散格式的穩定性分析，通常都是針對一維穩態無內熱源的問題進行的。因此，本文選取經典的一維穩態無內熱源的對流擴散方程作為研究對象。

其中，ρu和Г分別為對流和擴散項系數，ρu＞0，Г＞0。

為方便計算，空間網格采用均勻網格，網格數為K，網格間距為。各物性參數皆假定為常數，設Pe＝ρu·2Δxk／Γ，代表對流和擴散的強度之比，并分別采用3～6階勒讓德正交基函數展開來離散方程（12），最終可得到式（11）的形式，并進行進一步的整理和變換，可得到如下離散形式：

選取不同階數的勒讓德正交基函數展開，對應的R不同，表1為不同階基函數展開的節塊展開法離散形式對應的R。

表1 不同階基函數展開的節塊展開法離散形式對應的RTable 1 Rfor discretization scheme of NEM with different order basis functions

此外，中心差分格式對應的離散關系式為：

一階迎風格式對應的離散關系式經過變換，可整理為式（13）的形式為：

由式（14）、（15）及表1可知，不同階基函數展開的節塊展開法的離散格式以及有限體積法中的中心差分格式、一階迎風格式均可表示為相同的形式，區別僅是R不同。式（13）表示的離散形式經簡單地變化，即可寫成類似式（11）的三點關系式：

其中：A＝－（1＋R）；B＝2；C＝－（1－R）。式（16）的精確解為：

其中：λ1、λ2分別為式（16）對應的特征等式A＋Bλ＋Cλ2＝0的根；f1、f2的值由兩個邊界條件決定。

計算得到：

今年6月，李克強總理考察銀川經開區育成中心。面對圍攏過來的創業者、海外歸來的留學生、返鄉創客和自主擇業的軍轉干部，總理勉勵他們說：“這里有荒涼的沙漠，這里卻沒有荒涼的人生。希望你們用自己的汗水和智慧，把這里澆灌成創業創新的綠洲?！?/p>

由式（18）和（19）可得到滿足穩定性條件時R的取值范圍為：

根據R的取值范圍，即可得到不同階基函數展開的節塊展開法離散格式以及有限體積法中的中心差分格式、一階迎風格式的穩定性條件，表2列出上述所有離散格式滿足穩定性條件所對應的Pe取值范圍。其中，3NEM、4NEM、5NEM、6NEM分別代表3～6階勒讓德正交基函數展開對應的節塊展開法。

表2 不同階基函數展開的節塊展開法、中心差分及一階迎風格式的穩定性范圍Table 2 Stability of NEM with different order basis functions，center difference scheme and first order upwind difference scheme

由表2可知，有限體積法中的中心差分和一階迎風格式穩定性分析結果與文獻［11，13］均是一致的；而對于節塊展開法，當采用奇數階勒讓德正交基函數展開時，節塊展開法是無條件穩定的（如3階和5階基函數展開的節塊展開法）；當采用偶數階勒讓德正交基函數展開時，節塊展開法是條件穩定的，并隨展開階數的增大，其穩定性區域增大（如4階基函數展開的節塊展開法，Pe≤4.644 4；6階基函數展開的節塊展開法，Pe≤7.293 5），當Pe超過一定值后，就可能發生數值振蕩，導致不穩定。

3 數值耗散特性分析

在數值計算過程中，當對流項的離散格式不合理時，會人為地增大或削弱擴散項的作用，這就是所謂的數值耗散特性，也叫做假擴散［11］。當增大擴散項的作用時，稱其為正數值耗散，正數值耗散越大，數值穩定性越好，但卻增加了數值誤差；當削弱擴散項的作用時，稱其為負數值耗散，負數值耗散可能會引起數值的不穩定性，產生數值振蕩。

根據式（17）、（18）可得到各離散格式對應的精確解為：

即離散格式所計算的問題實際上是Peeff對應的問題，而并不是原始Pe對應的問題，Peeff越接近Pe，說明該離散格式所計算的問題越接近真實問題。根據計算，得到有限體積法中的中心差分、一階迎風及不同階基函數展開的節塊展開法的Pe與Peeff關系曲線示于圖2，其相對誤差（Peeff－Pe）／Pe示于圖3，不同離散格式的｜Peeff－Pe｜／Pe示于圖4。對于條件穩定的格式，僅示出穩定區域的Peeff。

圖2 不同離散格式Pe與Peeff的關系Fig.2 Relationship between Pe and Peeff

圖3 不同離散格式的（Peeff－Pe）／PeFig.3 （Peeff－Pe）／Pe for different discretization schemes

圖4 不同離散格式的｜Peeff－Pe｜／PeFig.4 ｜Peeff－Pe｜／Pe for different discretization schemes

由圖2～4可得出以下結論：

1）Pe較小時，隨節塊展開法的基函數展開階數的增大，Peeff越接近真實的Pe，產生的誤差越小，節塊展開法明顯優于有限體積法中的中心差分和一階迎風格式。

2）有限體積法中的中心差分格式及偶數階基函數展開的節塊展開法的相對誤差（Peeff－Pe）／Pe＞0，說明Peeff＞Pe，相當于增加對流項的作用，而削弱擴散項的作用，即產生負數值耗散，隨Pe的增大負數值耗散增大，最終引起不穩定性，產生數值振蕩。這解釋了有限體積法中的中心差分及偶數階基函數展開的節塊展開法是條件穩定的原因。

3）有限體積法中的一階迎風格式及奇數階基函數展開的節塊展開法的相對誤差（Peeff－Pe）／Pe＜0，說明Peeff＜Pe，相當于削弱對流項的作用，而增大擴散項的作用，即產生正數值耗散，隨Pe的增大，正數值耗散及誤差增大，當Pe＞PeD（PeD為節塊展開法的數值耗散等于有限體積法中的一階迎風格式的數值耗散時對應的Pe）時，奇數階基函數展開的節塊展開法的數值耗散極其明顯，甚至大于有限體積法中的一階迎風格式產生的數值耗散，使誤差大幅增加。對于3階基函數展開的節塊展開法，PeD＝6；對于5階基函數展開的節塊展開法，PeD＝14.319 8。

4）隨展開階數的增大，節塊展開法的精度越來越高，針對一維穩態無源問題，3～6階基函數展開的節塊展開法的精度分別接近4、6、8、10階。

4 數值計算結果與分析

為驗證上述穩定性分析、數值耗散特性分析及精度分析的正確性，分別選取一維無源基準問題，進行數值實驗驗證，其中基準問題的解析解均為已知。其中，ρu和Г為常數，空間網格采用均勻網格，設Pe＝2ρuh／Г（h為網格的半寬度）。

此外，為驗證不同離散格式的數值誤差，此處將采用均方根誤差（RMS）來表征各種離散格式對應的數值解與參考解的偏離程度，RMS定義為：

4.1 各離散格式的數值誤差及精度數值驗證

圖5 Pe＝1.6時不同離散格式對應的數值解與精確解分布情況Fig.5 Numerical solutions of different discretization schemes and exact solution for Pe＝1.6

圖5為Pe＝1.6時不同離散格式對應的數值解與原微分方程的精確解分布情況，圖6為不同離散格式的RMS隨網格尺寸的變化情況。

圖6 不同離散格式的RMSFig.6 RMS of different discretization schemes

由圖5可知，3～6階基函數展開的節塊展開法的數值計算結果與精確解吻合較好，而有限體積法中的一階迎風和中心差分格式數值解與精確解有較大的誤差。由圖6可知，隨展開階數的增大，節塊展開法產生的數值誤差減小、精度增加，各離散格式的精度與上述分析的基本一致。因此，節塊展開法明顯優于中心差分和一階迎風格式。

4.2 偶數階基函數展開的節塊展開法的穩定性數值驗證

為驗證偶數階基函數展開的節塊展開法的穩定性分析結果，表3列出不同Pe對應的數值計算結果。由表3可知，對于4階基函數展開的節塊展開法，當Pe＝4時，數值解未發生數值振蕩；當Pe＝5時，數值解發生數值振蕩，產生不穩定性。而對于6階基函數展開的節塊展開法，當Pe＝7時，數值解未發生數值振蕩；當Pe＝8時，數值解發生數值振蕩，產生不穩定性。以上結果與4階基函數展開的節塊展開法穩定性條件Pe≤4.644 4，6階基函數展開的節塊展開法穩定性條件Pe≤7.293 5分析結果吻合。

4.3 奇數階基函數展開的節塊展開法的數值耗散特性數值驗證

為驗證奇數階基函數展開的節塊展開法的數值耗散特性分析結果，選取不同Pe的數值計算結果，將各離散點數值解的絕對誤差示于圖7。由圖5、7可知，在Pe＝1.6時，3階和5階基函數展開的節塊展開法的數值耗散誤差均小于有限體積法中的一階迎風格式的；當Pe＝10時，3階基函數展開的節塊展開法的數值耗散誤差大于有限體積法中的一階迎風格式的，5階基函數展開的節塊展開法的數值耗散誤差仍然小于有限體積法中的一階迎風格式的；但當Pe＝20時，3階和5階基函數展開的節塊展開法的數值耗散誤差均大于有限體積法中的一階迎風格式的。因此，隨Pe的不斷增大，奇數階基函數展開的節塊展開法的數值耗散誤差增大，最終大于有限體積法中的一階迎風格式的，使計算結果誤差過大。

此外，將Pe＝50時各離散格式的數值計算結果及RMS示于圖8。由圖8可知：偶數階基函數展開的節塊展開法、有限體積法中的中心差分格式皆發生數值振蕩，奇數階基函數展開的節塊展開法的數值耗散遠超過有限體積法中的一階迎風；有限體積法中的一階迎風格式計算效果最好。由此可見，節塊展開法求解強對流占優問題仍存在一定的問題，需改進和提高。

表3 偶數階基函數展開的節塊展開法數值解Table 3 Numerical solutions of NEM with even order basis function

圖7 不同Pe時奇數階基函數展開的節塊展開法的絕對誤差Fig.7absolute error of NEM with odd order basis function for different Peclet numbers

圖8 Pe＝50時各離散格式的數值計算結果Fig.8 Numerical solutions of different discretization schemes for Pe＝50

5 總結

本文通過精確解和數值實驗方法，研究了不同階基函數展開的節塊展開法求解對流擴散方程的穩定性、數值耗散特性及誤差分析，并與有限體積法中的中心差分和一階迎風格式進行對比，可得出以下結論：偶數階基函數展開的節塊展開法是條件穩定的，即要求Pe＝2ρuh／Г小于Pe限值，且穩定性范圍和精度均優于有限體積法中的中心差分格式的，Pe限值會隨展開階數的增大而增大；而奇數階基函數展開的節塊展開法是無條件穩定的，但隨Pe的增大，數值耗散增大，計算誤差增大。當Pe大于一定值后，產生的數值耗散大于有限體積法中的一階迎風格式的。

節塊展開法求解強對流占優問題或網格劃分較粗時，會存在不穩定性或數值耗散過大等問題，這一定程度上影響了節塊展開法求解對流擴散方程的使用范圍和效率，需進行改進和提高。本文雖僅研究了一維對流擴散方程，但由于節塊展開法是利用橫向積分方法將多維對流擴散方程轉化為多個一維的對流擴散方程，因此研究一維節塊展開法離散形式對于節塊展開法的發展和改進具有重要的指導意義。

［1］ LAWRENCE R D.Program in nodal methods for the solution of the neutron diffusion and transport equations［J］.Progress in Nuclear Energy，1986，17（3）：271-301.

［2］ LAWRENCE R D，DORNING J J.A nodal Green’s function method for multidimensional neutron diffusion calculations［J］.Nuclear Science and Engineering，1980，76（2）：218-231.

［3］ SMITH K.An analytical nodal method for solving the two-group，multi-dimensional，static and transient neutron diffusion equations［D］.America：Massachusetts Institute of Technology，1979.

［4］ KIM Y I，KIM Y J.A semi-analytic multigroup nodal method［J］.Annals of Nuclear Energy，1999，26（8）：699-708.

［5］ WANG D Y，LI F，GUO J，et al.Improved nodal expansion method for solving neutron diffusion equation in cylindrical geometry［J］.Nuclear Engineering and Design，2010，240（8）：1 997-2 004.

［6］ MICHAEL E P E，DORNING J J，RIZWAN U.Studies on nodal integral methods for the convection-diffusion heat equation［J］.Nuclear Science and Engineering，2001，137（3）：380-399.

［7］ WANG F，RIZWAN U.Modified nodal integral method for the three-dimensional time-dependent incompressible Navier-Stokes equations［J］.Nuclear Science and Engineering，2005，149（1）：107-114.

［8］ SINGH S.Simulation of turbulent flows using nodal integral method［D］.America：University of Illinois，2008.

［9］ 鄧志紅，孫玉良，李富，等.改進的節塊展開法求解對流擴散方程的穩定性和誤差分析［J］.原子能科學技術，2014，58（2）：298-304.

DENG Zhihong，SUN Yuliang，LI Fu，et al.Stability and error analysis on modified nodal expansion method for transient convection-diffusion equation［J］.Atomic Energy Science and Technology，2014，58（2）：298-304（in Chinese）.

［10］鄧志紅，孫玉良，李富，等.改進的節塊展開法求解圓柱幾何對流擴散方程［J］.計算物理，2013，30（4）：475-482.

DENG Zhihong，SUN Yuliang，LI Fu，et al.Modified nodal expansion method for convectiondiffusion equation in cylindrical geometry［J］.Chinese Journal of Computational Physics，2013，30（4）：475-482（in Chinese）.

［11］陶文銓.數值傳熱學［M］.2版.西安：西安交通大學出版社，2001.

［12］GRESHO P M，LEE R L.Don’t suppress the wiggles：They are telling you something［J］.Computers ＆Fluids，1981，9（2）：223-251.

［13］PATANKAR S V.Numerical heat transfer and fluid flow［M］.America：Hemisphere Publishing Corporation，1980.

［14］KANG Y H，MICHAEL W.A method for reduction of numerical diffusion in the donor cell treatment of convection［J］.Journal of Computational Physics，1986，63（1）：201-221.

Stability and Numerical Diffusion Analysis of Nodal Expansion Method for Convection-diffusion Equation

ZHOU Xia-feng，LI Fu
（Key Laboratory of Advanced Reactor Engineering and Safety of Ministry of Education，Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing100084，China）

The stability and numerical diffusion analysis of nodal expansion method（NEM）for convection-diffusion equation was studied.The stabilities and numerical diffusion analyses for NEM with different order basis functions by exact solution of the discretization equation and numerical experiment results were done，and the results from NEM were compared with the results from both the center difference scheme and the first order upwind difference scheme of the finite volume method.The results show that the even order NEM is conditionally stable with Peclet number（Pe）being less than limit value，Pe limit value also increases with the order of expansion functions，and stability range and accuracy are better than those of center difference scheme.The odd order of NEM is unconditionally stable，however，with the increase of Pe，the numerical diffusion becomes lager and the calculation error becomes bigger.In addition，whenPe reaches certain value，the numerical diffusion even exceeds that of the first order upwind difference scheme.

nodal expansion method；convection-diffusion equation；stability；numerical diffusion

TL329

：A

：1000-6931（2015）04-0705-08

10.7538／yzk.2015.49.04.0705

2013-12-20；

2014-04-30

國家科技重大專項資助項目（ZX06901）；國家自然科學基金資助項目（11375099）

周夏峰（1989—），男，河南宜陽人，博士研究生，核科學與技術專業