劉 聰，胡小利，張 斌,*，陳義學(xué)

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206；2.中國核電工程有限公司，北京 100840)

離散縱標(biāo)法(SN)因其計算效率高、魯棒性好、適合求解深穿透問題等特點，是核裝置屏蔽計算的主要數(shù)值方法之一。角度離散誤差和空間離散誤差共同影響著離散縱標(biāo)法的數(shù)值模擬精度，保證計算可靠性需要控制與問題相關(guān)、影響較大的離散誤差。在深穿透屏蔽問題中，粒子穿透若干自由程厚的屏蔽層，粒子通量密度的強衰減對空間離散格式的精度與效率提出挑戰(zhàn)?；谟邢薏罘址椒ǖ目臻g離散格式盡管其局部截斷誤差小，但對網(wǎng)格步長較敏感，對于強吸收、粒子平均自由程較短的介質(zhì)需劃分較細(xì)網(wǎng)格以保證計算精度。另一方面，差分格式未考慮粒子沿不同離散方向衰減的物理特性，在網(wǎng)格長寬比不佳情況下難以保證粒子傳遞方向的正確性，引起非物理震蕩問題[1]。對于大規(guī)模三維深穿透屏蔽問題，一般受計算資源的限制，用戶需合理設(shè)置計算條件以平衡計算時間和精度。

短特征線離散[2]是一類SN方程空間變量離散方法，該方法根據(jù)入射、出射間的貢獻(xiàn)關(guān)系將計算網(wǎng)格劃分為若干子網(wǎng)格，通過與特征線解的空間矩守恒構(gòu)造出射和網(wǎng)格角通量密度解。Larsen和Alcouffe對一階空間矩采用差商近似發(fā)展了二維線性短特征線離散格式(LC)[3]，Childs和Rhoades利用線性節(jié)塊SN方程輸出近似一階矩，發(fā)展了簡化處理的三維LC方法[4]，上述LC方法犧牲了一階空間矩守恒、忽略了分片物理特性以換取更高的計算速度。自1990年至今，Azmy和Mathews分別基于短特征線離散在二維和三維幾何下發(fā)展了一系列空間離散格式[5-7]。其中，指數(shù)短特征線[8]具有正定性以及出色的粗網(wǎng)格計算精度，對于包含大塊屏蔽體的模擬問題具有較高的計算效率。Brennan和Miller利用網(wǎng)格分割、子網(wǎng)格投影變換發(fā)展了四面體網(wǎng)格的指數(shù)短特征線(EC)格式[9]，Suriano將上述算法擴(kuò)展至三維笛卡爾網(wǎng)格[10]。子網(wǎng)格投影至標(biāo)準(zhǔn)單位網(wǎng)格的處理過程對于非結(jié)構(gòu)網(wǎng)格求解是不可避免的，但對直角六面體網(wǎng)格而言，相關(guān)的幾何參數(shù)計算可通過適當(dāng)變換省略簡化，以節(jié)省用于幾何參數(shù)計算的時間成本。本文基于三維六面體結(jié)構(gòu)網(wǎng)格推導(dǎo)建立新的EC格式，網(wǎng)格分割后直接計算出射通量矩，以避免將待計算子網(wǎng)格投影至標(biāo)準(zhǔn)單元后再投影回原網(wǎng)格的求解過程，實現(xiàn)對計算流程和參數(shù)計算的優(yōu)化。

1 指數(shù)短特征線離散方法

三維幾何下m方向的穩(wěn)態(tài)線性玻爾茲曼輸運方程如式(1)，即待求解的SN方程。

Σtψm(x,y,z)=Qm(x,y,z)

(1)

式中：ψ為中子角通量密度；Σt為宏觀總截面；Q為源強；Ω為方向向量，省略能群符號。

利用方向?qū)?shù)改寫式(1)中的泄漏項得到特征線方程，其解沿射線變量s具有以下形式：

(2)

在三維xyz幾何下s可表示為式(3)，利用該關(guān)系可將式(2)轉(zhuǎn)換至目標(biāo)網(wǎng)格坐標(biāo)系下。

(3)

式中：(x0,y0,z0)為射線起點；(x,y,z)為射線終點；μ、η和ξ為方向向量與坐標(biāo)軸夾角余弦。

考察主象限內(nèi)某一離散方向，計算網(wǎng)格按入射來源可分為3個子網(wǎng)格(圖1a)，每個子網(wǎng)格根據(jù)出射面可被細(xì)分(圖1b)，如此根據(jù)出射與入射角通量密度間的貢獻(xiàn)關(guān)系，計算網(wǎng)格至多被分為7個子網(wǎng)格，每個子網(wǎng)格僅含1個入射面和1個出射面。

a——根據(jù)入射面分割；b——根據(jù)入射出射貢獻(xiàn)關(guān)系分割圖1 六面體網(wǎng)格分割示意圖Fig.1 Hexahedral mesh splitting

1.1 基函數(shù)和空間矩定義

在SN求解過程中網(wǎng)格入射條件和源強分布為已知值，由給定初值或上風(fēng)向網(wǎng)格出射值確定，EC格式假設(shè)入射角通量密度和源強為指數(shù)函數(shù)分布，即：

(4)

式中：i、j和k分別表示x、y、z方向相關(guān)參數(shù)；a和q為分布系數(shù)；b、c、d為對應(yīng)x、y、z方向的指數(shù)因子。

選取變形的勒讓德多項式作為基函數(shù)，基函數(shù)與角通量密度或源強的空間分布進(jìn)行內(nèi)積可得到空間矩，見式(5)和(6)。利用空間矩將角通量密度分布在網(wǎng)格之間進(jìn)行傳遞，零階矩為平均角通量密度，一階矩度量其分布形狀。

(5)

(6)

1.2 指數(shù)因子求解

目標(biāo)網(wǎng)格接收3個上風(fēng)向網(wǎng)格的出射面通量空間矩和本網(wǎng)格源空間矩作為已知條件，需解出其對應(yīng)分布的指數(shù)因子。以網(wǎng)格源矩為例，一階x矩與3倍零階矩相除可求解沿x方向指數(shù)因子b的非線性方程，即：

(7)

選取其他方向一階矩可得對應(yīng)指數(shù)因子的待求解方程，求解入射面角通量密度的指數(shù)分布與此相同。為方便初值選取，將待求解方程映射至第1象限，求解式(8)得到某方向的指數(shù)因子α，根據(jù)對稱性確定該因子正負(fù)。

(8)

式(8)為超越方程，其函數(shù)曲線如圖2所示，需進(jìn)行非線性數(shù)值求解，指數(shù)因子的求解算法需具有優(yōu)良的精度、效率和魯棒性。

當(dāng)α接近0或數(shù)值足夠大時，可對式(8)進(jìn)行泰勒展開或分解，選取主導(dǎo)項替代原方程以漸近解直接作為指數(shù)因子的值，節(jié)省計算量。

(9)

(10)

圖2 指數(shù)因子求解方程函數(shù)圖Fig.2 Function of exponential factor solving equation

1.3 角通量空間矩求解

由已知的入射和源強條件，需求解出射面和網(wǎng)格角通量密度分布，即零階和一階通量空間矩。從圖2可見，使方程f(α)-β=0有解的條件是β<1。通量空間矩的求解涉及若干包含多個指數(shù)項疊乘的多重積分計算，積分計算精度不足將導(dǎo)致指數(shù)因子求解方程無解，輸運計算無法進(jìn)行。由于指數(shù)因子的存在范圍廣、可正可負(fù)，多重積分存在可去奇點等特點，常規(guī)計算算法難以滿足EC格式要求。本文選用Mathews提出的多重指數(shù)矩[9]解決上述問題，其中指數(shù)矩定義如下：

(11)

因指數(shù)矩M積分上下限特點，在求解出射通量矩時需對圖1b的分割方案進(jìn)行圖3所示的細(xì)分。

圖3 EC出射面通量矩計算分割方案Fig.3 Splitting scheme for EC out-going flux spatial moment calculation

出射面通量矩分別源于邊界入射和網(wǎng)格源強貢獻(xiàn)，以圖1中x方向出射面為例說明推導(dǎo)過程。首先確定子網(wǎng)格范圍和粒子來源，利用式(2)和(3)表示出射面上的通量密度分布函數(shù)，來源于邊界入射和網(wǎng)格源強的零階空間矩即式(12)和(13)。

(12)

(13)

式中，角標(biāo)sub-1/2/3、in和Q分別表示子網(wǎng)格編號、入射貢獻(xiàn)和源強貢獻(xiàn)，總截面角標(biāo)省略。

變量代換改變積分變量使積分上下限與M指數(shù)矩相同，式(12)經(jīng)變換得式(14)。

(14)

將指數(shù)分布假設(shè)(式(4))代入改寫后的公式，經(jīng)過適當(dāng)?shù)暮喜⒒?，可將多重積分替換為指數(shù)矩函數(shù)M，式(12)和(13)最終分別被改寫為式(15)和(16)。

M0[ck(1-rzy)]

(15)

(16)

(17)

(18)

參照上述處理過程對y、z方向出射面空間矩進(jìn)行推導(dǎo)，當(dāng)某些子網(wǎng)格面積占比極小時取消該子網(wǎng)格計算提高算法穩(wěn)定性。EC格式是一種非線性空間離散方式，指數(shù)矩函數(shù)需要在輸運掃描中在線計算。指數(shù)矩函數(shù)計算占據(jù)了輸運計算大部分的時間，在Mathews和Miller的算法基礎(chǔ)上優(yōu)化了部分參數(shù)的計算，合理設(shè)置收斂準(zhǔn)則提高算法效率。對于網(wǎng)格通量空間矩，采用網(wǎng)格平衡方程(19)和一階矩平衡方程(20)進(jìn)行計算。網(wǎng)格空間矩用于更新網(wǎng)格源強，出射面空間矩傳遞給下風(fēng)向網(wǎng)格。

(19)

(20)

相比Suriano等的方法，本文建立了新的EC格式，通過對計算網(wǎng)格細(xì)分避免進(jìn)行單元投影映射，無需進(jìn)行雅可比矩陣和相關(guān)向量計算，指數(shù)因子求解由耦合非線性方程組問題降為單個非線性方程問題，網(wǎng)格空間矩由矩平衡方程直接給出，降低了計算成本，并對求解算法進(jìn)行了優(yōu)化。EC格式具有恒正性同時要求保證非負(fù)源，其他求解步驟與標(biāo)準(zhǔn)SN求解相同。

2 基準(zhǔn)驗證

2.1 空間離散格式計算速度比較

為衡量不同空間離散格式的計算速度，選取“box-in-box”簡單幾何固定源問題進(jìn)行模擬。網(wǎng)格劃分為25×25×25，采用S8階求積組。粒子散射的各向異性在實際屏蔽計算中不可忽略，散射源計算的相關(guān)開銷同樣對計算時間有較大影響，因此選取P1、P3和P5階散射展開進(jìn)行測試。以P1展開、帶置零修正菱形差分(DZ)計算時間為參考，各算例的歸一化計算時間列于表1，表中括號內(nèi)是以相同展開階數(shù)DZ用時為參考的歸一化時間。指數(shù)定向差分(EDW)[11]是一種基于預(yù)估指數(shù)分布的差分格式，提高了粗網(wǎng)格精度，但具有差分格式無法正確考慮射線傳遞過程的共同缺點。ARES輸運程序中二維LC格式[12]采用差商近似計算一階矩，同時引入斜率旋轉(zhuǎn)置零修正，三維LC精確考慮各子網(wǎng)格零階和一階矩平衡，根據(jù)網(wǎng)格類型預(yù)先計算所需的多重積分矩。由表1可見，EC格式的單網(wǎng)格計算成本遠(yuǎn)高于其他幾種格式的。實際屏蔽模擬多采用P3或P5散射展開，隨著通量角度矩、散射源計算占用的時間增加，不同空間離散格式用時差距將縮小。EC格式主要用于大塊強吸收屏蔽體區(qū)域，其較高的單網(wǎng)格計算成本可由其出色的粗網(wǎng)格精度補償，具體數(shù)值分析結(jié)果見后文。

表1 簡單固定源問題歸一化計算時間Table 1 Normalized calculation time of simple fixed-source problem

2.2 Kobayashi基準(zhǔn)題

Azmy和Mathews等的研究表明，當(dāng)短特征線離散用于空腔介質(zhì)計算時，舍入誤差問題可能引起數(shù)值精度退化、迭代不穩(wěn)定。本文EC格式采用矩平衡方程計算網(wǎng)格矩，節(jié)省了大量多重積分的計算量，但由于總截面處于分母項，計算空腔介質(zhì)對求解涉及的參數(shù)精度要求較高。經(jīng)測試發(fā)現(xiàn)，采用指數(shù)矩函數(shù)改寫的公式及優(yōu)化后的參數(shù)計算算法具有較高的魯棒性。選取Kobayashi基準(zhǔn)題[13]問題一半散射算例驗證EC格式對于含空腔介質(zhì)問題的精度和穩(wěn)定性，基準(zhǔn)題幾何模型示于圖4。x、y、z負(fù)方向邊界為反射邊界，其余為真空邊界，圖中10～50 cm區(qū)域材料截面為1.0×10-4cm-1，用于模擬空腔介質(zhì)，其余區(qū)域材料截面為0.1 cm-1，中心0～10 cm區(qū)域立方體為源區(qū)。

圖4 Kobayashi基準(zhǔn)題問題1幾何模型Fig.4 Sketch of Kobayashi benchmark problem 1

考慮大塊真空區(qū)域和弱散射介質(zhì)的存在，該基準(zhǔn)題主要用于考驗SN程序控制角度離散誤差、消除射線效應(yīng)的能力。為盡可能降低角度離散誤差，選取S40階勒讓德切比雪夫求積組，EC格式的關(guān)鍵點計算結(jié)果和蒙特卡羅參考基準(zhǔn)[14]示于圖5，其中1A區(qū)域沿直線x=z=5 cm，1B區(qū)域沿對角線x=y=z，1C沿直線y=55 cm、z=5 cm。計算結(jié)果相對偏差小于±4%，本文的EC格式可正確處理粒子在空腔介質(zhì)中的輸運過程，保證迭代穩(wěn)定收斂和結(jié)果正確。

圖5 Kobayashi基準(zhǔn)題問題1關(guān)鍵點計算結(jié)果Fig.5 Results of key points for Kobayashi benchmark problem 1

2.3 MMS基準(zhǔn)問題

(21)

式中:wm為求積組離散方向權(quán)重;M為方向總數(shù)。

該方法通過調(diào)整各邊界入射和結(jié)合源分布以構(gòu)造出不同光滑程度的通量密度解。本文設(shè)置結(jié)合源為階梯遞減分布，構(gòu)造具有函數(shù)連續(xù)、指數(shù)衰減特征的通量解分布。x、y和z負(fù)方向為反射邊界，其余為真空邊界，模型邊長為20 cm，材料總截面為1 cm-1，散射比分別設(shè)為0.01、0.1、0.3、0.5和0.7，網(wǎng)格步長設(shè)為1 cm和2.5 cm，采用S16階EON求積組，參考基準(zhǔn)分布示于圖6。

根據(jù)式(22)統(tǒng)計系統(tǒng)標(biāo)通量密度L2范數(shù)誤差因子，該因子衡量系統(tǒng)相對誤差總和。由于模擬采用的網(wǎng)格步長已超過粒子平均自由程，該基準(zhǔn)問題主要考驗空間離散在粗網(wǎng)格條件下的計算精度。

(22)

式中：NIJK為網(wǎng)格總數(shù)；ΔVijk為網(wǎng)格體積。

誤差因子的統(tǒng)計結(jié)果示于圖7，對于該問題EC和LC格式精度相當(dāng)且優(yōu)于其他方法。隨著散射比的提高，網(wǎng)格內(nèi)通量密度衰減減弱，EC和LC的精度優(yōu)勢隨之減弱。網(wǎng)格步長從1 cm變?yōu)?.5 cm，除SC格式外，其他空間離散格式誤差因子均有所上升。由于結(jié)合源按階梯狀分布，區(qū)域交界處通量密度下降較快，而各區(qū)域內(nèi)部為平源且通量解分布相對平坦，較符合常數(shù)特征線SC的基礎(chǔ)假設(shè)。全局解的非物理振蕩問題造成DZ、EDW和θ權(quán)重差分TW 3種差分格式的誤差因子較大。

圖6 MMS基準(zhǔn)問題Fig.6 MMS benchmark problem

2.4 平板入射問題

MMS基準(zhǔn)問題構(gòu)造了全局網(wǎng)格的基準(zhǔn)解，但由于結(jié)合源分布的設(shè)置，通量密度衰減速度慢。為模擬屏蔽問題中的深穿透、強衰減特征，通過平板入射問題比較空間離散格式的精度和效率，模型左邊界給定常數(shù)入射條件，右邊界為真空邊界，其余邊界為反射邊界。模型x方向長50 cm，材料總截面為1 cm-1，散射比為0.5，粒子通量密度在系統(tǒng)中將經(jīng)歷超過20個量級的衰減。

比較3種具有較好粗網(wǎng)格精度的空間離散格式EDW、LC和EC，沿x方向25～50 cm間的逐網(wǎng)格通量密度計算結(jié)果示于圖8，參考解為網(wǎng)格步長為0.05 cm的LC解，圖例注明了網(wǎng)格尺寸。當(dāng)網(wǎng)格步長足夠大時需考慮網(wǎng)格內(nèi)的通量密度變化，EC基于指數(shù)分布假設(shè)可重構(gòu)出通量密度分布，圖8d展示了EC粗網(wǎng)格下的重構(gòu)解。

a——網(wǎng)格步長1 cm；b——網(wǎng)格步長2.5 cm圖7 MMS基準(zhǔn)問題誤差因子結(jié)果Fig.7 Error factor result of MMS benchmark problem

a——EDW結(jié)果；b——LC結(jié)果；c——EC結(jié)果；d——EC重構(gòu)解結(jié)果圖8 平板入射問題不同空間離散格式的計算結(jié)果Fig.8 Results of plane incident problems under different spatial discretization schemes

對于該問題，LC和EDW在網(wǎng)格2 cm以下時精度相當(dāng)，隨著網(wǎng)格步長繼續(xù)增大，EDW解趨于保守，LC由于負(fù)解的產(chǎn)生將低估粒子通量密度。在測試的全部網(wǎng)格條件下，EC保證了較高的計算精度，當(dāng)網(wǎng)格步長為20倍自由程時也正確模擬了粒子衰減過程。由細(xì)網(wǎng)LC解獲得重構(gòu)位置處的參考值，按式(23)統(tǒng)計了EC重構(gòu)解的均方根相對偏差列于表2。

(23)

表2 25～50 cm間的EC重構(gòu)解均方根相對偏差Table 2 RMS relative deviations of EC reconstructed solutions from 25 cm to 50 cm

該平板入射問題由于無分布源存在，粒子通量密度呈現(xiàn)強衰減分布，由表2結(jié)果可見，EC在極端粗網(wǎng)條件下也可保證較高計算精度，達(dá)到同一精度要求下所需網(wǎng)格數(shù)量極大降低。根據(jù)本例題EC重構(gòu)解結(jié)果和網(wǎng)格量，達(dá)到均方根偏差小于10%的精度水平，考慮在三維幾何下，EC格式所需網(wǎng)格數(shù)量相比EDW和LC格式將減少約(5/2)3至(20/2)3倍，EC格式較高的單網(wǎng)格計算成本問題可由其所需網(wǎng)格量減少解決，相比帶權(quán)重差分等粗網(wǎng)格精度不佳的空間離散方法，EC效率優(yōu)勢將更為顯著。對于大塊厚屏蔽體，EC格式可在較少網(wǎng)格數(shù)量條件下達(dá)到所需精度。

3 結(jié)論

本文基于離散縱標(biāo)方法提出了新的指數(shù)特征線空間離散格式，通過適當(dāng)?shù)淖泳W(wǎng)格分割和變量代換，避免進(jìn)行單元投影映射，節(jié)省了幾何參數(shù)計算量，借助指數(shù)矩函數(shù)保證了算法穩(wěn)定性，指數(shù)因子求解和指數(shù)矩函數(shù)計算相關(guān)算法在已有研究基礎(chǔ)上進(jìn)行了優(yōu)化改進(jìn)。通過4個簡單幾何固定源問題的模擬分析，展現(xiàn)了EC格式出色的粗網(wǎng)格計算精度。綜合考慮單網(wǎng)格計算成本和粗網(wǎng)格下的計算精度，EC格式對于實際屏蔽問題的計算成本可接受。對于含有大塊均勻屏蔽體的三維深穿透屏蔽問題，相比其他空間離散方法，EC格式達(dá)到相同計算精度所需網(wǎng)格量可減少約1～3個量級，在未來自適應(yīng)空間離散策略研究中EC將作為有力補充。合適的網(wǎng)格劃分是發(fā)揮EC方法計算能力和效率的前提，網(wǎng)格步長設(shè)置過小將極大增加其計算成本，未來將利用間斷非匹配網(wǎng)格增強EC對實際復(fù)雜幾何模型的應(yīng)用能力。