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

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

改進的節(jié)塊展開法求解對流擴散方程的穩(wěn)定性和誤差分析

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

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

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

改進的節(jié)塊展開法；對流擴散方程；穩(wěn)定性；誤差分析

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

1 基本方程

穩(wěn)態(tài)三維對流擴散方程為：

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

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

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

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

節(jié)塊的x方向的偏溫度展開如下：

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

該指數(shù)函數(shù)滿足：

節(jié)塊邊界處的偏溫度和偏溫度梯度為：

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

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

其中，wn（x）為權重函數(shù)。

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

2 MNEM穩(wěn)定性分析

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

考慮一維穩(wěn)態(tài)線性對流擴散方程為：

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

其中，R＝ρcpu／λ。

根據(jù)符號不變原則，穩(wěn)定性保持的條件為：

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

3 數(shù)值結果與分析

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

3.1 一維問題

1）一維無源

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

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

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

2）一維有源

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

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

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

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

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

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

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

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

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

3.2 二維無源問題

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

其中：

其解析解［5］為：

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

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

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

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

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

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

4 結論

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

［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—），男，湖南衡陽人，博士研究生，核科學與技術專業(yè)