本發(fā)明針對(duì)核反應(yīng)堆堆芯中子學(xué)計(jì)算領(lǐng)域，提出了一種針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法。

背景技術(shù)：

核反應(yīng)堆中子學(xué)計(jì)算研究以核反應(yīng)堆堆芯為應(yīng)用對(duì)象，其堆芯由許多不同種類的組件構(gòu)成。根據(jù)堆型的不同，組件內(nèi)部的幾何結(jié)構(gòu)和材料布置復(fù)雜多變。因此，實(shí)際的反應(yīng)堆中子學(xué)問題是一個(gè)三維非均勻幾何的中子學(xué)問題。對(duì)核反應(yīng)堆進(jìn)行快速、精確的中子學(xué)計(jì)算，是反應(yīng)堆設(shè)計(jì)和校核的基本保障。目前在堆芯物理設(shè)計(jì)過程中，主要針對(duì)均勻幾何的進(jìn)行擴(kuò)散方程的求解。

變分節(jié)塊法是核反應(yīng)堆物理設(shè)計(jì)中常用的方法之一。它以二階偶宇稱形式的中子擴(kuò)散方程為出發(fā)點(diǎn)，方程呈現(xiàn)橢圓方程的形式，有利于garlerkin方法的應(yīng)用，更適合有限元方法的空間離散。變分節(jié)塊法的計(jì)算思想是：首先通過變分方法在非均勻求解區(qū)域建立包含二階中子輸運(yùn)方程和自然邊界條件的泛函；然后采用標(biāo)準(zhǔn)正交多項(xiàng)式進(jìn)行ritz離散，同時(shí)利用球諧函數(shù)實(shí)現(xiàn)角度展開，并構(gòu)造響應(yīng)矩陣；最后分別在三維堆芯的各個(gè)節(jié)塊內(nèi)分別求解響應(yīng)矩陣方程；節(jié)塊之間以流和其高階矩耦合，最終得到問題區(qū)域的中子通量密度分布。變分節(jié)塊法能夠達(dá)到較高的精度，但僅能處理均勻節(jié)塊的問題。在基于均勻節(jié)塊的中子學(xué)計(jì)算過程中，無法避免均勻化過程，以及隨之而來的誤差。隨著科學(xué)技術(shù)的發(fā)展和計(jì)算機(jī)水平的提高，人們已經(jīng)越來越開始重視減少近似和假設(shè)，追求更高精度的中子學(xué)計(jì)算方法，消除均勻化過程是中子學(xué)發(fā)展的必然趨勢(shì)。因此，研究具備非均勻幾何處理能力的計(jì)算方法對(duì)于中子學(xué)計(jì)算具有十分重要的意義。

基于變分節(jié)塊法，后續(xù)研究在變分節(jié)塊法的框架下提出了非均勻幾何變分節(jié)塊方法，通過引入有限元方法以精確描述柵元內(nèi)部的非均勻結(jié)構(gòu)。在此工作中，各柵元被單獨(dú)處理為節(jié)塊，節(jié)塊內(nèi)部采用有限元網(wǎng)格進(jìn)行細(xì)化，精確描述柵元和冷卻劑的材料、幾何分布。為了能夠?qū)嶋H描述柵元中的曲邊幾何結(jié)構(gòu)，且保證足夠的計(jì)算精度，往往需要大量的有限元網(wǎng)格剖分，因此計(jì)算量很大。以圖1中所示的一個(gè)壓水堆柵元的有限元剖分為例，柵元包含燃料和冷卻劑兩區(qū)，分別以紅色和藍(lán)色標(biāo)示，黑點(diǎn)代表有限元網(wǎng)格的節(jié)點(diǎn)。采用二階等參有限元對(duì)其進(jìn)行非結(jié)構(gòu)網(wǎng)格剖分，共包含32個(gè)有限元，對(duì)應(yīng)97個(gè)有限元節(jié)點(diǎn)。在實(shí)際計(jì)算過程中，需要存儲(chǔ)和求解反應(yīng)堆堆芯中每一個(gè)柵元中97個(gè)有限元節(jié)點(diǎn)上的中子通量密度展開矩，會(huì)花費(fèi)很大的計(jì)算代價(jià)。

技術(shù)實(shí)現(xiàn)要素：

為了有效地減小非均勻幾何變分節(jié)塊方法的計(jì)算代價(jià)，本發(fā)明提出了一種針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法，可以大幅度地減少計(jì)算時(shí)間和計(jì)算內(nèi)存。

為了實(shí)現(xiàn)上述目的，本發(fā)明采取了以下技術(shù)方案予以實(shí)施：

一種針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法，通過將有限元?jiǎng)澐譃槠皆磪^(qū)，以提高計(jì)算效率和減小計(jì)算內(nèi)存，步驟如下：

步驟1：在各平源區(qū)中，將各個(gè)有限元節(jié)點(diǎn)的源項(xiàng)近似為平源區(qū)內(nèi)的平均源項(xiàng)，通過此近似減小非均勻求解區(qū)域中的自由度數(shù)目；首先給出公式(1)和公式(2)響應(yīng)矩陣方程的表達(dá)式；根據(jù)公式(4)中的中子通量密度分布φ(x,y)的表達(dá)式，求得公式(5)中平源區(qū)的平均中子通量密度展開矩和中子通量密度分布表達(dá)式φ(x,y)的關(guān)系：

在非均勻幾何變分節(jié)塊方法中，節(jié)塊內(nèi)部中子通量密度展開矩的求解方程為

其中：

-1—矩陣的求逆；

φ—節(jié)塊內(nèi)部中子通量密度展開矩向量，其中展開矩代表展開系數(shù)的值；

j—節(jié)塊表面凈中子流密度展開矩向量；q—中子源項(xiàng)展開矩向量；

—響應(yīng)矩陣，二者表達(dá)形式不同，僅與節(jié)塊內(nèi)部的材料布置、幾何形狀有關(guān)；

響應(yīng)矩陣方程為：

式中：

j⁺—出射中子流密度展開矩向量；

j^-—入射中子流密度展開矩向量；

u—中子流源項(xiàng)展開矩向量；

—響應(yīng)矩陣，二者表達(dá)形式不同，僅與節(jié)塊內(nèi)部的材料布置、幾何形狀有關(guān)；

針對(duì)二維情況，已知中子源的空間分布表示為：

式中：

q(x,y)—中子源的空間分布；

σsgg′—第g群到第g’群的中子散射截面；

g,g′—兩個(gè)不同的能群標(biāo)識(shí)；

φg′(x,y)—第g^/群的中子通量密度分布；

keff—有效增殖因子；

νσfg′—第g’群的中子產(chǎn)生截面；

χg—第g群的中子裂變譜；

二維情況下，節(jié)塊內(nèi)部中子角通量密度的離散表達(dá)式與軸向無關(guān)：

φ(x,y)≈g^t(x,y)φ公式(4)

其中：

φ(x,y)—節(jié)塊內(nèi)部中子通量密度的分布；

g(x,y)—x‐y方向有限元形狀函數(shù)向量；

t—向量或矩陣的轉(zhuǎn)置；

φ—節(jié)塊內(nèi)部中子通量密度的展開矩向量；

在平源加速假設(shè)下，平源區(qū)的平均中子通量密度展開矩能夠表示為中子通量密度分布在本平源區(qū)內(nèi)部的積分平均值：

式中：

—各平源區(qū)平均中子通量的矩向量；

ht(x,y)—對(duì)應(yīng)于每個(gè)有限元網(wǎng)格的分片常量：ht(x,y)＝δe，e為每個(gè)有限元網(wǎng)格；

—對(duì)角線元素為每個(gè)有限元網(wǎng)格面積的對(duì)角矩陣；

步驟2：根據(jù)公式(5)中平源區(qū)的平均中子通量密度展開矩和中子通量密度分布φ(x,y)表達(dá)式之間的關(guān)系，代入中子通量密度展開矩的公式(1)中，就能夠得到公式(6)中平源區(qū)的平均中子通量密度展開矩和有限元節(jié)點(diǎn)上中子通量密度展開矩φ的關(guān)系；

將公式(5)代入公式(1)得：

其中：

—平源區(qū)的平均中子通量密度展開矩；

同時(shí)

步驟3：利用公式(6)中平源區(qū)的平均中子通量密度展開矩和有限元節(jié)點(diǎn)上中子通量密度展開矩φ之間的轉(zhuǎn)換關(guān)系，代入傳統(tǒng)變分節(jié)塊方法中響應(yīng)矩陣方程的表達(dá)式公式(1)和公式(2)，將求解對(duì)象從各個(gè)有限元節(jié)點(diǎn)上的未知量轉(zhuǎn)換為各個(gè)平源區(qū)的未知量，即公式(13)和公式(14)；

將公式(6)代入到公式(2)得：

此外，群內(nèi)源項(xiàng)能夠表示為：

假設(shè)每個(gè)平源區(qū)內(nèi)有限元節(jié)點(diǎn)上的源項(xiàng)與平源區(qū)的源項(xiàng)相等：

式中：

σx(x,y)—平源區(qū)內(nèi)部不同位置處的x反應(yīng)的截面，x可以為吸收，總截面，散射；

—各平源區(qū)x反應(yīng)截面所組成的對(duì)角矩陣；

因此，將公式(6)代入公式(9)，得：

其中：

將公式(11)代入公式(8)得：

同理，公式(2)變?yōu)椋?/p>

步驟4：利用通用的迭代方法，對(duì)公式(13)和公式(14)進(jìn)行求解，從而得到整個(gè)非均勻求解區(qū)域的中子通量密度分布和中子流密度分布，實(shí)現(xiàn)非均勻幾何變分節(jié)塊方法的求解。

本發(fā)明通過將計(jì)算對(duì)象從各個(gè)有限元節(jié)點(diǎn)上中子通量密度展開矩向量φ，轉(zhuǎn)變?yōu)榍蠼飧鱾€(gè)平源區(qū)內(nèi)平均中子通量密度展開矩向量展開矩向量，能夠降低存儲(chǔ)內(nèi)存和計(jì)算時(shí)間，實(shí)現(xiàn)針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法。

與現(xiàn)有技術(shù)相比，本發(fā)明具有如下突出優(yōu)點(diǎn)：

1.本發(fā)明通過采用平源加速方法劃分每個(gè)有限元為一個(gè)平源區(qū)，能夠顯著減少柵元中自由度的數(shù)目。以圖1為例，若采用本發(fā)明中的方法，柵元中的自由度可以從97減少到32，從而顯著降低中子通量密度展開矩的存儲(chǔ)計(jì)算內(nèi)存。

2.本發(fā)明通過利用公式(14)中的表達(dá)式，能夠?qū)⒂?jì)算對(duì)象從各個(gè)有限元節(jié)點(diǎn)上中子通量密度展開矩向量φ，轉(zhuǎn)變?yōu)榍蠼飧鱾€(gè)平源區(qū)內(nèi)平均中子通量密度展開矩向量減小解向量的長(zhǎng)度，降低響應(yīng)矩陣的規(guī)模，從而節(jié)省求解過程中的浮點(diǎn)數(shù)操作數(shù)目，顯著提高計(jì)算效率。

附圖說明

圖1等參有限元描述下的壓水堆柵元非均勻幾何。

具體實(shí)施方式

下面結(jié)合具體實(shí)施方式對(duì)本發(fā)明作進(jìn)一步詳細(xì)說明：

本發(fā)明一種針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法，通過將有限元?jiǎng)澐譃槠皆磪^(qū)，以提高計(jì)算效率和減小計(jì)算內(nèi)存，步驟如下：

步驟1：在各平源區(qū)中，將各個(gè)有限元節(jié)點(diǎn)的源項(xiàng)近似為平源區(qū)內(nèi)的平均源項(xiàng)，通過此近似減小如圖1所示的非均勻求解區(qū)域中的自由度數(shù)目；首先給出公式(1)和公式(2)中的響應(yīng)矩陣方程的表達(dá)式；根據(jù)公式(4)中的中子通量密度分布φ(x,y)的表達(dá)式，求得公式(5)中平源區(qū)的平均中子通量密度展開矩和中子通量密度分布表達(dá)式φ(x,y)的關(guān)系：

在非均勻幾何變分節(jié)塊方法中，節(jié)塊內(nèi)部中子通量密度展開矩的求解方程為

其中：

-1—矩陣的求逆；

φ—節(jié)塊內(nèi)部中子通量密度展開矩向量，其中展開矩代表展開系數(shù)的值；

j—節(jié)塊表面凈中子流密度展開矩向量；

q—中子源項(xiàng)展開矩向量；

—響應(yīng)矩陣，二者表達(dá)形式不同，僅與節(jié)塊內(nèi)部的材料布置、幾何形狀有關(guān)；

響應(yīng)矩陣方程為：

式中：

j⁺—出射中子流密度展開矩向量；

j^-—入射中子流密度展開矩向量；

u—中子流源項(xiàng)展開矩向量；

—響應(yīng)矩陣，二者表達(dá)形式不同，僅與節(jié)塊內(nèi)部的材料布置、幾何形狀有關(guān)。

以二維情況為例，已知中子源的空間分布可以表示為：

式中：

q(x,y)—中子源的空間分布；

σsgg′—第g群到第g’群的中子散射截面；

g,g′—兩個(gè)不同的能群標(biāo)識(shí)；

φg′(x,y)—第g’群的中子通量密度分布；

keff—有效增殖因子；

νσfg′—第g’群的中子產(chǎn)生截面；

χg—第g群的中子裂變譜。

二維情況下，節(jié)塊內(nèi)部中子角通量密度的離散表達(dá)式與軸向無關(guān)：

φ(x,y)≈g^t(x,y)φ公式(4)

其中：

φ(x,y)—節(jié)塊內(nèi)部中子通量密度的分布；

g(x,y)—x‐y方向有限元形狀函數(shù)向量；

t—向量或矩陣的轉(zhuǎn)置；

φ—節(jié)塊內(nèi)部中子通量密度的展開矩向量；

在平源加速假設(shè)下，平源區(qū)的平均中子通量密度展開矩能夠表示為中子通量密度分布在本平源區(qū)內(nèi)部的積分平均值：

式中：

—各平源區(qū)平均中子通量的矩向量；

ht(x,y)—對(duì)應(yīng)于每個(gè)有限元網(wǎng)格的分片常量：ht(x,y)＝δe，e為每個(gè)有限元網(wǎng)格；

—對(duì)角線元素為每個(gè)有限元網(wǎng)格面積的對(duì)角矩陣。

步驟2：根據(jù)公式(5)中平源區(qū)的平均中子通量密度展開矩和中子通量密度分布φ(x,y)表達(dá)式之間的關(guān)系，代入中子通量密度展開矩的公式(1)中，就能夠得到公式(6)中平源區(qū)的平均中子通量密度展開矩和有限元節(jié)點(diǎn)上中子通量密度展開矩φ的關(guān)系；

將公式(5)代入公式(1)得：

其中：

—平源區(qū)的平均中子通量密度展開矩；

同時(shí)

步驟3：利用公式(6)中平源區(qū)的平均中子通量密度展開矩和有限元節(jié)點(diǎn)上中子通量密度展開矩φ之間的轉(zhuǎn)換關(guān)系，代入傳統(tǒng)變分節(jié)塊方法中響應(yīng)矩陣方程的表達(dá)式公式(1)和公式(2)，將求解對(duì)象從各個(gè)有限元節(jié)點(diǎn)上的未知量轉(zhuǎn)換為各個(gè)平源區(qū)的未知量，即公式(13)和公式(14)；

將公式(6)代入到公式(2)得：

此外，群內(nèi)源項(xiàng)能夠表示為：

假設(shè)每個(gè)平源區(qū)內(nèi)有限元節(jié)點(diǎn)上的源項(xiàng)與平源區(qū)的源項(xiàng)相等：

式中：

σx(x,y)—平源區(qū)內(nèi)部不同位置處的x反應(yīng)的截面，x可以為吸收，總截面，散射等；

—各平源區(qū)x反應(yīng)截面所組成的對(duì)角矩陣；

因此，將公式(6)代入公式(9)，得：

其中：

將公式(11)代入公式(8)得：

同理，公式(2)變?yōu)椋?/p>

步驟4：對(duì)于非均勻幾何求解區(qū)域內(nèi)不同種類的節(jié)塊，分別計(jì)算各節(jié)塊內(nèi)部的空間相關(guān)響應(yīng)矩陣；

步驟5：針對(duì)特定能群，利用紅‐黑掃描的方式對(duì)響應(yīng)矩陣方程公式(14)進(jìn)行迭代求解，得到出、入射偏中子流密度展開矩j⁺、j^-；

步驟6：在求得出、入射偏中子流密度展開矩j⁺、j^-后，代入公式(13)中，即可解得節(jié)塊內(nèi)部對(duì)應(yīng)各個(gè)平源區(qū)的平均中子通量密度展開矩向量從而得知節(jié)塊內(nèi)部的中子通量密度分布；

步驟7：進(jìn)入下一能群的計(jì)算，最終求解整個(gè)非均勻幾何求解區(qū)域的中子通量密度分布和中子流密度分布。通過將計(jì)算對(duì)象從各個(gè)有限元節(jié)點(diǎn)上中子通量密度展開矩向量φ，轉(zhuǎn)變?yōu)榍蠼飧鱾€(gè)平源區(qū)內(nèi)平均中子通量密度展開矩向量能夠降低存儲(chǔ)內(nèi)存和計(jì)算時(shí)間，實(shí)現(xiàn)針對(duì)非均勻幾何變分節(jié)塊方法的平源加速方法。

驗(yàn)證結(jié)果顯示，在保證精度的前提下，平源加速方法可將計(jì)算速度提高4倍以上，內(nèi)存減少70％。