改進的節塊展開法求解對流擴散方程的穩定性和誤差分析

鄧志紅1，孫玉良1，李 富1，Rizwan-uddin2

（1.清華大學 核能與新能源技術研究院，北京 100084；2.伊利諾伊大學 香檳分校，伊利諾伊州 61801，美國）

改進的節塊展開法求解對流擴散方程的穩定性和誤差分析

鄧志紅1，孫玉良1，李 富1，Rizwan-uddin2

（1.清華大學 核能與新能源技術研究院，北京 100084；2.伊利諾伊大學 香檳分校，伊利諾伊州 61801，美國）

對改進的節塊展開法（MNEM）求解對流擴散方程進行了深入研究。利用符號不變原則理論上分析了MNEM的穩定性，分析結果顯示，MNEM是一種具有固有穩定性的求解方法。通過數值實驗對MNEM的計算誤差進行了分析。數值結果表明：對于一維問題，MNEM至少具有3階精度；對于多維問題，受橫向泄漏近似的影響，MNEM表現為2階精度。

改進的節塊展開法；對流擴散方程；穩定性；誤差分析

現代節塊法，因其在精度和效率方面的優勢，已廣泛應用于反應堆物理計算中。但在反應堆安全分析程序中，熱工問題的計算依然普遍采用有限差分或有限體積方法求解。為改進物理-熱工計算效率，在節塊法的構架下統一求解物理-熱工耦合問題是樸素而直接的想法。應用節塊法求解反應堆熱工問題的研究相對較少，其中，節塊積分法（NIM）在求解對流擴散方程和Navier-Stokes方程方面已取得不錯進展［12］。但其對虛擬源項通常采用0階近似，無法實現物理-熱工的高階耦合，損失了計算精度。研究表明，傳統的節塊展開法求解對流擴散方程只具有條件穩定性（網格臨界Peclet數小于4.644），其穩定性雖優于中心差分方法，但在對流占優的情況下采用粗網節塊將會帶來非物理的數值振蕩，精度很差［3］。為了保持節塊展開法對高階源項的處理能力，且改善節塊展開法的求解穩定性，本文提出改進的節塊展開法（MNEM）來求解對流擴散方程，并基于離散格式穩定性分析方法（符號不變原則）對MNEM的穩定性進行理論分析，根據數值實驗對其求解對流擴散方程的誤差進行分析。

1 基本方程

穩態三維對流擴散方程為：

式中：T為溫度；ρ為密度；cp為比定壓熱容；λ為導熱系數；u、v、w分別為3個方向的速度；S為源項。

與求解中子擴散方程類似，節塊展開法求解對流擴散方程的一般過程為：首先采用橫向積分技術將多維的偏微分方程轉化為多個一維的常微分方程；對節塊內的偏溫度采用基函數進行展開，利用Galerkin剩余權重法得到3個矩權重方程；利用節塊界面處的連續條件（偏溫度連續和偏溫度梯度連續）最終得到截面偏溫度的離散方程；利用平衡方程求解節塊的平均溫度。

圖1 節塊劃分示意圖Fig.1 Domain discretization for nodal expansion method

對流擴散方程與中子擴散方程不同，由于對流項的存在，其解可能發生劇烈變化，尤其當對流占優時。由于對流擴散方程的解劇烈振蕩，傳統差分方法在求解該方程時會遇到穩定性問題，通常需選取非常細的網格來避免數值振蕩，因此極大限制了其計算效率。為改善傳統節塊展開法（NEM）求解對流擴散方程的穩定性，本文中提出的MNEM采用了速度的指數函數代替原來第4階多項式作為展開函數，數值結果表明其極大地提高了傳統NEM的穩定性。

節塊的x方向的偏溫度展開如下：

在MNEM中，代替原來第4階展開多項式的指數函數為：

該指數函數滿足：

節塊邊界處的偏溫度和偏溫度梯度為：

在MNEM中，與傳統NEM類似，源項仍采用2階多項式展開。但其泄漏項的處理與中子擴散方程不同，由于對流擴散方程的解變化劇烈，傳統NEM中關于橫向泄漏的假設（3個相鄰節塊的橫向泄漏形狀相同，且為2階多項式）已不再成立。數值實驗表明，求解對流擴散方程時若仍采用原來的橫向泄漏近似方法，會由于簡單近似帶來的數值振蕩（尤其當對流占優時）影響方法的整體計算精度。因此，本文僅對橫向泄漏進行0階近似，即：

在MNEM中，Galerkin剩余權重法被用來產生3個矩權重方程，其過程為：

其中，wn（x）為權重函數。

選取w0（x）＝1、w1（x）＝fx1（x）、w2（x）＝fx2（x），將得到3個矩權重方程，分別為：

2 MNEM穩定性分析

為了說明MNEM較傳統NEM求解對流擴散方程的優勢，本文將從理論上分析MNEM的穩定性。穩定性分析通常基于簡單的一維問題［4］，且在NEM中，多維的偏微分方程將通過橫向積分技術變為多個一維的常微分方程，因此，基于一維問題的分析向多維問題拓展是直接的。

考慮一維穩態線性對流擴散方程為：

為便于穩定性分析，假設上述參數在計算區域內恒定，且計算區域被劃分為相同尺寸節塊，節塊半寬度為a。由于流速為0時，對流擴散方程蛻化為擴散方程，通常無需考慮穩定性問題。因此，推導過程中假設流速不為0。

其中，R＝ρcpu／λ。

根據符號不變原則，穩定性保持的條件為：

由式（17）、（18）可知，無論速度的大小和方向，MNEM的離散形式始終滿足穩定性條件（式（19）），即網格Peclet數可為任意值。本文給出的穩定性分析雖是基于簡單的一維、常系數、定溫邊界條件下得到的，但研究表明，此時得到的網格臨界Peclet數是較為苛刻條件下的值，對于原假設條件的偏離均會使得網格臨界Peclet數有所增加［5］。因此，最終可得到穩定性分析的結論：MNEM是一種具有固有穩定性的方法。

3 數值結果與分析

基于文中的推導過程，開發了求解對流擴散方程的程序。

3.1 一維問題

1）一維無源

其中，C＝ρcpu／λ。為了考察MNEM的穩定性和計算精度，圖2示出了C＝-100、-10、0.5、10、100時MNEM的計算結果。節塊的尺寸劃分均勻，均為0.2。為了更好地描述MNEM的計算精度，圖2也示出了MNEM在節塊中的詳細溫度分布。

圖2 不同C值時MNEM的計算結果與解析解的對比Fig.2 Comparison between MNEM and analytical solution for different C

由圖2可知，對于不同的C值，MNEM的計算結果均與解析解符合得非常好，即使當C＝±100、溫度梯度很大時，MNEM依然能很好地跟蹤溫度的變化。且無論流速的大小和方向，MNEM的計算結果均保持數值穩定性，表明MNEM具有固有的迎風特性。

2）一維有源

為進一步驗證MNEM求解帶源問題，本文將計算帶分布源的一維對流擴散問題：

圖3示出C＝1 000時，不同節塊數目下MNEM的數值解同解析解的比較。

由圖3可知，節塊數N＝2、5、10時均能得到很好的結果，即計算得到的節塊平均值同解析解平均值符合得非常好。且計算區域中溫度的變化劇烈，即便如此，MNEM采用非常少的節塊（如N＝2）依然能得到穩定、高精度的結果，表征了MNEM計算對流擴散問題的能力。

圖3 MNEM數值解與解析解的對比Fig.3 Comparison between MNEM and analytical solution

為進一步分析MNEM的特點，表1列出C取不同值時MNEM和NIM計算本算例的均方根（RMS）誤差隨節塊尺寸的變化情況。RMS誤差定義如下：

其中，Ti為數值解，Ti，ref為與之對應的參考解。

由表1可知，在各種C和節塊劃分的情況下，MNEM的RMS誤差均比NIM的小，精度優勢非常明顯。為比較二者的精度階次，圖4示出RMS誤差隨節塊尺寸變化的對數關系。由圖4可見，MNEM的誤差在對數圖中的斜率最小為3.10，表征計算該一維帶源問題，MNEM的最低收斂精度為3階。NIM的誤差斜率則全部為2.0，符合文獻［1］中采用0階源項展開時NIM是2階精度的論斷。

表1 MNEM和NIM的RMS誤差隨節塊尺寸的變化Table 1 RMS errors of MNEM and NIM versus node-size

圖4 MNEM和NIM求解一維問題的RMS誤差隨節塊尺寸的變化Fig.4 RMS errors of MNEM and NIM versus node-size for 1Dproblem

3.2 二維無源問題

二維無源問題的計算公式為：

其中：

其解析解［5］為：

為說明MNEM的計算精度，選定計算區域中y方向的中線位置（y＝0.5）處的區域。圖5示出C分別為100和1 000時，不同節塊尺寸下MNEM和NIM的數值解與解析解的結果比較。由圖5可見，MNEM的數值解同解析解均符合得非常好，且同NIM的結果一致。需要說明的是，由于圖5給出的均是節塊平均溫度，因此不同節塊劃分情況下的溫度曲線存在一定的差異。

圖5 MNEM和NIM的數值解同解析解的對比Fig.5 Numerical solution of MNEM and NIM versus analytical solution

為說明MNEM計算多維對流問題的計算精度，表2列出C分別為100和1 000、不同節塊劃分時MNEM與NIM的RMS誤差。由表2可見，在各種節塊劃分情況下，MNEM和NIM的誤差基本相當，但NIM的誤差稍微優于MNEM的。

表2 二維無源問題MNEM和NIM的RMS誤差隨節塊尺寸的變化Table 2 RMS errors of MNEM and NIM versus node-size for 2Dproblem without source

圖6示出C為100和1 000時，MNEM和NIM的RMS誤差隨節塊尺寸的變化。由圖6可見，二者無論是誤差值還是誤差的變化率均基本一致，表明對于二維無源問題，二者的計算精度相當，且均為2階。本算例的結果與一維有源問題的結果存在差異的原因是：本算例中源項為零，因而橫向泄漏近似是本問題中唯一的近似。由于MNEM和NIM均采用橫向泄漏項的0階近似，因此二者的計算結果與精度相當。

圖6 MNEM和NIM求解二維問題的RMS誤差隨節塊尺寸的變化Fig.6 RMS errors of MNEM and NIM versus node-size for 2Dproblem

4 結論

本文深入研究了MNEM求解穩態對流擴散方程的特性。基于一維線性無源對流擴散方程，利用符號不變原則從理論上分析了MNEM的穩定性，分析表明，該方法具有固有的穩定性，對網格Peclet數無限制。基于一系列數值實驗研究了MNEM求解對流擴散方程的精度，計算結果表明：對于一維問題，MNEM至少具有3階精度；對于多維問題，由于MNEM目前仍采用簡單的橫向泄漏近似，其精度表現為2階。后續工作將進一步研究有效的橫向泄漏近似方法，使MNEM計算對流擴散問題的精度整體提高到3階，同時將研究節塊法求解流動方程的可行性，從而最終實現在節塊法的框架下求解物理-熱工問題。

［1］ Rizwan-uddin.A second-order space and time nodal method for the one-dimensional convectiondiffusion equation［J］.Computers ＆Fluids，1997，26（3）：233-247.

［2］ TOREJA A J.A nodal approach to arbitrary geometries，and adaptive mesh refinement for the nodal method［D］.USA：University of Illinois，Urbana，IL，2002.

［3］ DENG Z H，Rizwan-uddin，LI F，et al.Stability and error analysis of nodal expansion method for convection-diffusion equation［C］∥Proceedings of ICAPP2012.Chicago：［s.n.］，2012.

［4］ TAO W Q，SPARROW E M.The transportive property and convective numerical stability of the steady-state convection-diffusion finite-difference equation［J］.Numerical Heat Transfer，1987，11（4）：491-497.

［5］ GUPTA M M，MANOHAR R P，STEPHENSON J W.A single cell high order scheme for the convection-diffusion equation with variable coefficients［J］.International Journal for Numerical Methods in Fluids，1984，7（4）：641-651.

Stability and Error Analysis on Modified Nodal Expansion Method for Transient Convection-diffusion Equation

DENG Zhi-hong1，SUN Yu-liang1，LI Fu1，Rizwan-uddin2
（1.Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing100084，China；2.Nuclear，Plasma，Radiological Engineering，University of Illinois at Urbana-Champaign，Urbana 61801，USA）

To further investigate the features of modified nodal expansion method（MNEM）for solving the convection-diffusion equation，the stability and error analysis were carried out.Based on sign preservation principle，the stability analysis reveals that the MNEM has inherent stability.The error analysis was implemented through a series of numerical experiments，and the results show that the MNEM is 3rd order scheme for one dimensional problem，while as 2nd order scheme for multi-dimensional problem because of using simple transverse leakage approximation.

modified nodal expansion method；convection-diffusion equation；stability；error analysis

O24

A

1000-6931（2014）02-0298-07

10.7538／yzk.2014.48.02.0298

2012-10-16；

2013-01-15

國家重大科技專項經費資助項目（ZX06901）

鄧志紅（1984—），男，湖南衡陽人，博士研究生，核科學與技術專業