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

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

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

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

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

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

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

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

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

1 基本方程

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

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

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

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

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

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

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

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

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

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

在MNEM中，Galerkin剩余權(quán)重法被用來產(chǎn)生3個(gè)矩權(quán)重方程，其過程為：

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

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

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

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

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

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

其中，R＝ρcpu／λ。

根據(jù)符號(hào)不變?cè)瓌t，穩(wěn)定性保持的條件為：

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

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

基于文中的推導(dǎo)過程，開發(fā)了求解對(duì)流擴(kuò)散方程的程序。

3.1 一維問題

1）一維無源

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

圖2 不同C值時(shí)MNEM的計(jì)算結(jié)果與解析解的對(duì)比Fig.2 Comparison between MNEM and analytical solution for different C

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

2）一維有源

為進(jìn)一步驗(yàn)證MNEM求解帶源問題，本文將計(jì)算帶分布源的一維對(duì)流擴(kuò)散問題：

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

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

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

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

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

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

二維無源問題的計(jì)算公式為：

其中：

其解析解［5］為：

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

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

為說明MNEM計(jì)算多維對(duì)流問題的計(jì)算精度，表2列出C分別為100和1 000、不同節(jié)塊劃分時(shí)MNEM與NIM的RMS誤差。由表2可見，在各種節(jié)塊劃分情況下，MNEM和NIM的誤差基本相當(dāng)，但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時(shí)，MNEM和NIM的RMS誤差隨節(jié)塊尺寸的變化。由圖6可見，二者無論是誤差值還是誤差的變化率均基本一致，表明對(duì)于二維無源問題，二者的計(jì)算精度相當(dāng)，且均為2階。本算例的結(jié)果與一維有源問題的結(jié)果存在差異的原因是：本算例中源項(xiàng)為零，因而橫向泄漏近似是本問題中唯一的近似。由于MNEM和NIM均采用橫向泄漏項(xiàng)的0階近似，因此二者的計(jì)算結(jié)果與精度相當(dāng)。

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

4 結(jié)論

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

國家重大科技專項(xiàng)經(jīng)費(fèi)資助項(xiàng)目（ZX06901）

鄧志紅（1984—），男，湖南衡陽人，博士研究生，核科學(xué)與技術(shù)專業(yè)