牛進(jìn)林，鄔穎杰，郭 炯，李 富

(清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

隨著計(jì)算機(jī)計(jì)算能力的提升，采用全隱式耦合JFNK方法求解復(fù)雜核電廠(chǎng)系統(tǒng)已成為一個(gè)趨勢(shì)。本課題組[1-4]對(duì)JFNK方法應(yīng)用于高溫堆耦合計(jì)算的相關(guān)關(guān)鍵技術(shù)進(jìn)行了長(zhǎng)期研究，取得了一系列成果，主要是基于原有程序修改其迭代方法。國(guó)外很多JFNK的研究與應(yīng)用走基于平臺(tái)開(kāi)發(fā)的路線(xiàn)，以MOOSE(Multiphysics Object Oriented Simulation Environment)[5-6]平臺(tái)為代表。MOOSE是由美國(guó)愛(ài)達(dá)荷國(guó)家實(shí)驗(yàn)室針對(duì)大規(guī)模耦合計(jì)算問(wèn)題開(kāi)發(fā)的開(kāi)源通用平臺(tái)，目前國(guó)外基于MOOSE平臺(tái)開(kāi)發(fā)了Pronghorn、RattleSnake、Bison、Relap7[7-12]等中子物理、燃料性能分析和系統(tǒng)分析程序，但這些能用于反應(yīng)堆的程序?qū)χ袊?guó)是出口管制的，國(guó)內(nèi)目前無(wú)法獲得這些程序，這使得國(guó)內(nèi)基于MOOSE平臺(tái)開(kāi)發(fā)反應(yīng)堆多物理耦合程序必須從頭做起。本文基于MOOSE平臺(tái)開(kāi)發(fā)可用于計(jì)算中子本征值和動(dòng)力學(xué)的多群擴(kuò)散程序，并通過(guò)基準(zhǔn)題進(jìn)行驗(yàn)證。

1 MOOSE對(duì)本征值和瞬態(tài)問(wèn)題的實(shí)現(xiàn)

1.1 MOOSE程序開(kāi)發(fā)框架

MOOSE[5-6]是一面向?qū)ο蟮亩辔锢眈詈铣绦蜷_(kāi)發(fā)平臺(tái)，采用魯棒性的求解方法，具備以下能力：1) 與維度無(wú)關(guān)的用戶(hù)代碼，1套代碼對(duì)一、二、三維問(wèn)題都適用；2) 基于連續(xù)或非連續(xù)Galerkin有限元的網(wǎng)格離散；3) 全隱式耦合；4) 允許采用非結(jié)構(gòu)化網(wǎng)格，包含多種形狀函數(shù)；5) 很強(qiáng)的網(wǎng)格自適應(yīng)性；6) 內(nèi)置并行功能，用戶(hù)不用開(kāi)發(fā)額外的并行代碼；7) 內(nèi)置的后處理等。

MOOSE簡(jiǎn)化了開(kāi)發(fā)多物理程序的步驟，提供了一個(gè)面向?qū)ο?，包含?gòu)成模擬程序各方面的插入式系統(tǒng)。MOOSE程序開(kāi)發(fā)包含3個(gè)層級(jí)，如圖1所示，開(kāi)發(fā)者只需將精力集中于中間層和頂層的開(kāi)發(fā)，省去了大量重復(fù)性的底層工作，這種模塊化的編程方式能大幅加快開(kāi)發(fā)程序的速度，并有利于程序擴(kuò)展和維護(hù)[5-6]。

圖1 MOOSE層級(jí)Fig.1 MOOSE hierarchy

1.2 本征值

多群擴(kuò)散本征值問(wèn)題控制方程如下：

(1)

式中：D為擴(kuò)散系數(shù)，cm；g為當(dāng)前能群編號(hào)；g′為散射能群編號(hào)；G為最大能群數(shù)；φ為中子通量，cm-2·s-1；χ為裂變份額；Σr為移出截面，cm-1；Σs為散射截面，cm-1；Σf為裂變截面，cm-1；ν為平均有效裂變中子數(shù)；keff為有效增殖因數(shù)。

在MOOSE中，控制方程經(jīng)過(guò)Galerkin有限元離散，處理成弱解形式的殘差方程。為達(dá)到此目的，控制方程需要與由1組試函數(shù)構(gòu)成的矩陣在求解域內(nèi)進(jìn)行內(nèi)積。矩陣形式如下：

(2)

式中：B為基函數(shù)組；B為展開(kāi)基函數(shù)；L為節(jié)點(diǎn)數(shù)目。將式(2)中的矩陣與式(1)進(jìn)行內(nèi)積得到弱解形式的殘差方程(式(3))，方程的解由展開(kāi)系數(shù)和基函數(shù)的積組合而成，如式(4)所示。

(3)

(4)

式中：F為非線(xiàn)性方程組殘差；U為解向量[φg]；U為有限元展開(kāi)系數(shù)；(·)表示體積分；〈·〉表示面積分。式(3)可進(jìn)一步改寫(xiě)成式(5)形式，其中中子不考慮向上散射。

Fφ(U)=Aφ-Mφ/keff

(5)

(6)

(7)

(8)

對(duì)于本征值問(wèn)題，傳統(tǒng)是采取冪迭代法進(jìn)行求解，Knoll等[13]于2009年實(shí)現(xiàn)了用JFNK求解本征值問(wèn)題。其基本思路是在式(5)的基礎(chǔ)上補(bǔ)充一個(gè)關(guān)于有效增殖因數(shù)的殘差方程，如式(9)所示。

(9)

式中，Ω為整個(gè)求解域。

聯(lián)立式(5)、(9)，并采用JFNK方法[14]對(duì)殘差方程進(jìn)行求解，從而得到本征值方程的解{φg，keff}。

在MOOSE中，物理場(chǎng)都是以kernel，即算子(以下統(tǒng)稱(chēng)為算子)的形式存在的。每個(gè)殘差方程包含多個(gè)算子，每個(gè)算子皆對(duì)應(yīng)1個(gè)C++類(lèi)，只有把殘差方程的所有算子都開(kāi)發(fā)出來(lái)，控制方程才是完備的。式(3)包括5個(gè)算子類(lèi)型，每個(gè)圓括號(hào)代表1個(gè)全求解域的算子，主要包括擴(kuò)散項(xiàng)、碰撞移出項(xiàng)、群間散射源項(xiàng)、裂變?cè)错?xiàng)；尖括號(hào)代表邊界條件算子。雖然每個(gè)能群殘差方程的算子類(lèi)型會(huì)有重復(fù)，但由于群截面和物性的不同，都需重新定義本能群的算子。而相關(guān)的物性需在額外的物性模塊定義，需有相關(guān)的C++類(lèi)相匹配。目前已實(shí)現(xiàn)了多群擴(kuò)散方程各算子的開(kāi)發(fā)，完成了常物性開(kāi)發(fā)，也實(shí)現(xiàn)了與變量耦合的變物性。式(9)的殘差方程是MOOSE提供的，不需要額外開(kāi)發(fā)。

陳大勇睜圓了雙眼，雙手握緊鬼子的槍刺，一聲虎嘯，刺刀被他帶著噴血拔了出來(lái)，小鬼子驚呆了，手一撒，就那么兩秒的時(shí)間，陳大勇反手把刺刀扎進(jìn)了鬼子的胸膛。

1.3 瞬態(tài)計(jì)算

考慮多組緩發(fā)中子的多群擴(kuò)散瞬態(tài)方程如下：

g=1,2，…，G

(10)

k=1,2，…，K

(11)

式中：C為緩發(fā)中子先驅(qū)核濃度，cm-3；keff,0為系統(tǒng)初始時(shí)刻的有效增殖因數(shù)；v為當(dāng)前群中子平均速度，cm/s；k為緩發(fā)中子先驅(qū)核組標(biāo)；K為緩發(fā)中子先驅(qū)核組數(shù)；λ為緩發(fā)中子先驅(qū)核衰變常量，s-1；β為緩發(fā)中子總份額；βk為第k組緩發(fā)中子先驅(qū)核份額。

將式(10)、(11)在求解域分別與式(2)進(jìn)行內(nèi)積，得到弱解形式的殘差方程：

(12)

(13)

聯(lián)立式(12)、(13)，采用JFNK方法對(duì)殘差方程進(jìn)行求解，得到瞬態(tài)方程的解{φ，C}。

在本征值程序的基礎(chǔ)上，殘差方程(12)主要是增加了瞬態(tài)項(xiàng)和緩發(fā)中子源項(xiàng)，裂變?cè)错?xiàng)因緩發(fā)中子的緣故需重新調(diào)整；開(kāi)發(fā)了緩發(fā)中子殘差方程(13)的相關(guān)算子和物性。

1.4 預(yù)處理矩陣

(14)

式中：J為雅可比矩陣；N為預(yù)處理陣；v為子空間基向量；ε為差分步長(zhǎng)。

目前程序中預(yù)處理陣采用塊對(duì)角矩陣形式，不考慮不同變量間的耦合關(guān)系。本征值問(wèn)題采用的預(yù)處理陣為式(15)，瞬態(tài)問(wèn)題采用的預(yù)處理陣為式(16)。

(15)

(16)

(17)

Nkeff,keff=1

(18)

NCk,i,Ck,j=(Bi,λkBj)

(19)

式中，i、j分別為試函數(shù)和展開(kāi)基函數(shù)的組標(biāo)。

1.5 邊界條件

MOOSE提供了第1類(lèi)和第2類(lèi)邊界條件的接口，但缺乏真空邊界條件和反照率邊界條件。因此開(kāi)發(fā)這兩個(gè)邊界條件的接口對(duì)中子計(jì)算是很有必要的，現(xiàn)已在程序中完成了真空邊界條件和反照率邊界條件的接口開(kāi)發(fā)。對(duì)于反照率邊界條件，在邊界處滿(mǎn)足式(20)。

(20)

式中，α為反照率。當(dāng)反照率取0時(shí)，得到真空邊界條件，如式(21)所示。

(21)

1.6 MOOSE實(shí)現(xiàn)

基于MOOSE開(kāi)發(fā)的中子程序基本框架如圖2所示，中間淡藍(lán)色框圖是MOOSE的基本框架，外圍灰色部分是1.2～1.5節(jié)的內(nèi)容在MOOSE基本框架下的實(shí)現(xiàn)，也是程序開(kāi)發(fā)的核心，其中中子擴(kuò)散方程、緩發(fā)中子方程、邊界條件、物性參數(shù)、預(yù)處理器是構(gòu)成中子程序開(kāi)發(fā)的關(guān)鍵技術(shù)難點(diǎn)。目前尚無(wú)專(zhuān)門(mén)的截面庫(kù)與本程序?qū)?，現(xiàn)階段是根據(jù)特定問(wèn)題的截面信息定制一套物性參數(shù)。MOOSE提供了多種預(yù)處理的接口，除塊對(duì)角預(yù)處理外還有有限差分預(yù)處理、基于物理的預(yù)處理等，目前在本程序中主要是實(shí)現(xiàn)塊對(duì)角的預(yù)處理方式。式(3)、(12)、(13)的各物理場(chǎng)算子構(gòu)成了程序的基本模塊，將各基本模塊按照問(wèn)題需求在輸入卡中像搭積木一樣“拼接”起來(lái)便構(gòu)成了基本程序，如圖3所示，本征值問(wèn)題和瞬態(tài)問(wèn)題所需的基本模塊是不完全一致的。中子程序流程圖如圖4所示。

圖2 中子程序框架Fig.2 Neutron program diagram

2 數(shù)值驗(yàn)證

通過(guò)2D-TWIGL基準(zhǔn)題[15]對(duì)程序進(jìn)行驗(yàn)證。本基準(zhǔn)題包括1個(gè)穩(wěn)態(tài)算例和2個(gè)瞬態(tài)算例。

圖3 模塊化程序架構(gòu)Fig.3 Modularized program architecture diagram

圖4 中子程序流程圖 Fig.4 Neutron program flow diagram

2.1 基準(zhǔn)題描述

TWIGL基準(zhǔn)題是一簡(jiǎn)化的二維雙群擴(kuò)散加1組緩發(fā)中子的動(dòng)力學(xué)問(wèn)題。反應(yīng)堆是正方形的，邊長(zhǎng)160 cm，包含3個(gè)物質(zhì)區(qū)域，其1/4堆芯幾何模型和邊界條件如圖5所示。其截面參數(shù)列于表1。緩發(fā)中子先驅(qū)核的份額為0.007 5，衰變常量為0.08 s-1。

圖5 TWIGL基準(zhǔn)題幾何模型Fig.5 TWIGL benchmark geometry

2.2 本征值計(jì)算

采用結(jié)構(gòu)化網(wǎng)格，共6 561個(gè)網(wǎng)格點(diǎn)，6 400個(gè)有限元，如圖6所示。本征值計(jì)算結(jié)果為：數(shù)值解0.913 20，參考解0.913 21[15]，兩者間相對(duì)誤差為0.001%?？梢?jiàn)計(jì)算結(jié)果與參考解符合得很好，精度很高?？烊汉蜔崛褐凶油糠植既鐖D7所示。MOOSE提供冪迭代(PI)求解本征值和用切比雪夫加速的接口，1次冪迭代也可看作1次非線(xiàn)性迭代。JFNK求解本征值和經(jīng)過(guò)契比雪夫加速的冪迭代進(jìn)行對(duì)比，求解的本征值是相同的，計(jì)算性能結(jié)果對(duì)比列于表2，殘差收斂結(jié)果對(duì)比示于圖8，發(fā)現(xiàn)在同樣的收斂精度下，JFNK能大幅減少非線(xiàn)性迭代的步數(shù)，從而大幅減少計(jì)算時(shí)間，計(jì)算效率提高接近1個(gè)數(shù)量級(jí)。另外可看出，即使經(jīng)過(guò)契比雪夫加速，冪迭代的殘差在非線(xiàn)性步也不是單調(diào)下降的，在小范圍會(huì)有波動(dòng)上升。

表1 TWIGL截面參數(shù)Table 1 TWIGL cross section parameter

2.3 瞬態(tài)計(jì)算

采用一階向后歐拉差分處理時(shí)間項(xiàng)，以穩(wěn)態(tài)為初始時(shí)刻?？紤]兩個(gè)瞬態(tài)工況：1) 物質(zhì)區(qū)1的熱群吸收截面在0時(shí)刻突然減少0.003 5 cm-1，而后物性保持不變，模擬0.5 s(工況1)；2) 物質(zhì)區(qū)1的熱群吸收截面從0時(shí)刻開(kāi)始，在0.2 s時(shí)間內(nèi)線(xiàn)性減少0.003 5 cm-1，而后物性保持不變，模擬0.5 s(工況2)。兩個(gè)工況的時(shí)間步長(zhǎng)均分別選取0.01、0.05、0.1 s，參考解[15]時(shí)間步長(zhǎng)為0.001 s，模擬結(jié)果與參考解的對(duì)比示于圖9?？煽闯?，無(wú)論是急速瞬變還是緩慢瞬變，程序的跟蹤能力都很好，數(shù)值結(jié)果和模擬結(jié)果相一致；時(shí)間步長(zhǎng)取0.01 s和0.05 s時(shí)，結(jié)果的精度也是有保障的；時(shí)間步長(zhǎng)取0.1 s(參考解時(shí)間步長(zhǎng)的100倍)時(shí)，結(jié)果會(huì)出現(xiàn)一些偏差，這是由于時(shí)間步長(zhǎng)如果取得過(guò)大，時(shí)間離散帶來(lái)的誤差就會(huì)變大，因此需要綜合考慮時(shí)間離散帶來(lái)的誤差和計(jì)算效率來(lái)選取時(shí)間步長(zhǎng)。

圖6 TWIGL基準(zhǔn)題網(wǎng)格Fig.6 TWIGL benchmark mesh

圖7 快群和熱群中子通量分布Fig.7 Fast neutron and thermal neutron flux distribution

表2 JFNK和冪迭代計(jì)算性能對(duì)比Table 2 Performance comparison of JFNK and PI

3 結(jié)論

本文基于MOOSE平臺(tái)開(kāi)發(fā)了可用于中子擴(kuò)散的求解器，采用JFNK方法計(jì)算本征值問(wèn)題和瞬態(tài)問(wèn)題，開(kāi)發(fā)了反照率邊界條件和真空邊界條件接口。通過(guò)基準(zhǔn)題TWIGL驗(yàn)證，發(fā)現(xiàn)對(duì)于本征值問(wèn)題，JFNK方法精度高、收斂快，與傳統(tǒng)的采用切比雪夫加速的冪迭代方法相比，能大幅減少非線(xiàn)性步的數(shù)量，極大提高計(jì)算效率。對(duì)于瞬態(tài)問(wèn)題，瞬變的跟蹤能力強(qiáng)，在較大的時(shí)間步長(zhǎng)下，也能保證一定的精度。

圖8 殘差收斂性對(duì)比Fig.8 Residual convergence comparison

圖9 兩個(gè)工況的結(jié)果對(duì)比 Fig.9 Result comparison of two operation conditions

本程序具有基于MOOSE平臺(tái)開(kāi)發(fā)程序的優(yōu)勢(shì)，1套代碼能解決多維問(wèn)題，對(duì)規(guī)則幾何或不規(guī)則幾何都有強(qiáng)大的模擬能力，能很方便地進(jìn)行并行計(jì)算等。之后還可對(duì)程序進(jìn)行進(jìn)一步的優(yōu)化，如選取不同的預(yù)處理矩陣、采用高階有限元、高階時(shí)間步長(zhǎng)、并行計(jì)算等，從而不斷提高求解精度和計(jì)算效率。在未來(lái)的研究中，中子程序可與熱工計(jì)算在同一框架下進(jìn)行全隱式耦合，為實(shí)現(xiàn)反應(yīng)堆的多物理耦合計(jì)算提供支持。但現(xiàn)階段還沒(méi)有可與中子程序?qū)拥暮藬?shù)據(jù)庫(kù)接口，限制了程序的模擬能力，這也是需要重點(diǎn)解決的問(wèn)題。