付元光1)2) 鄧力1) 李剛1)

1)(北京應用物理與計算數(shù)學研究所,北京 100088)2)(中物院高性能數(shù)值模擬軟件中心,北京 100088)(2017年12月14日收到;2018年4月13日收到修改稿)

1 引 言

核系統(tǒng)運行過程中,核材料不斷被輻照,其核素組成隨時間發(fā)生著復雜變化.在核反應堆物理計算中,通常用燃耗方程來描述核素含量隨時間的變化規(guī)律:

其中n是材料中核素的種類;Ni是第i種核素處于任意時刻t的含量;Ni0是第i種核素初始時刻的含量;λij是第j種核素反應(包括中子誘發(fā)、衰變等)產生第i種核素的轉換系數(shù);λi是第i種核素發(fā)生反應(包括中子誘發(fā)、衰變等)的消亡系數(shù);fi是第i種核素的遷移率.燃耗方程也可寫成矩陣-向量形式

其中N是n維核素含量向量,N0是N在t=0時刻的值,A是由系數(shù)λij和λi構成的n×n維燃耗矩陣.在實際問題中,燃耗方程包含數(shù)百至上千種核素,核素之間存在復雜的轉換關系,嚴格來講,描述核素產生、消亡的系數(shù)λij和λi是隨時間變化的,但變化比較緩慢,在適當?shù)臅r間步長內,可將λij和λi視為常量,使燃耗方程成為一階線性常微分方程組,再通過數(shù)值方法進一步求解.

求解燃耗方程有兩大技術路線,其一是線性子鏈方法(TTA)[1?3],通過解耦核素之間的轉換關系,生成一系列線性子鏈,再根據(jù)常微分方程理論逐一解析求解每條鏈;線性子鏈算法精度高,但鏈的構造過程十分耗時,同時為保證算法穩(wěn)定性,需要更費時的高精度浮點計算支持;其二是矩陣方法,通過數(shù)值方法求解矩陣tA的指數(shù)etA,主要有泰勒展開法[4]、Krylov子空間法[5]、切比雪夫有理近似法 (CRAM)[6?8]、Padé近似法[9]等,矩陣方法求解速度快很多,部分方法與TTA精度相當.目前,大部分燃耗數(shù)值計算程序主要借助以上算法求解齊次燃耗方程問題[10],即在(1)式中令fi=0,這是由于在常見的核反應堆中,不同空間區(qū)域之間核素遷移不顯著,燃耗區(qū)相對獨立,求解齊次方程已能夠滿足精度.對于一些新型核能系統(tǒng)(如熔鹽堆、液態(tài)散裂靶等),具有顯著的核素遷移,必須求解非齊次燃耗方程.在眾多程序中,ORIGEN[4]和CINDER90[11]是兩款可以求解非齊次燃耗方程的程序,但僅能處理fi為常量的情況,這種情況已能夠滿足大部分實際應用需求.對于CINDER90,雖然能給出正確解,但計算耗時很長;對于ORIGEN,如果指定fi＜0,可能導致結果發(fā)散.

本文研究了非齊次項fi具有特殊含時形式時方程的解法,在若干已有程序只能計算fi為常量的基礎上,擴展了fi適用范圍.設fi能夠在t=0附近(t=t0類似,只需要一個平移變換)通過有限階泰勒展開有效逼近,即

其中al是展開項向量,則非齊次燃耗方程變?yōu)?/p>

基于方程(4)的形式,本文使用了兩種求解方法:首先推導出了基于TTA方法的解析解形式,然后計算了基于矩陣方法的矩陣級數(shù)解的近最佳有理逼近表達式,通過不同方法的計算結果比對,驗證了兩種方法求解特定形式非齊次燃耗方程的可行性.

2 理論推導

2.1 非齊次燃耗方程基于TTA方法的解

TTA方法的基本思想是將核素之間的復雜轉換關系解耦成若干線性鏈,在一條線性鏈上,一個核素只能經由一個母核反應生成,且該核素發(fā)生反應只能生成一個子核,只有鏈首核素具有遷移率和初始核子密度,其余核素都為0.這種簡單的傳遞關系和初始條件使得描述一條線性子鏈的微分方程組可解析求解.逐一求解每一條線性子鏈,再把每條鏈的貢獻累加即得到結果.描述一條線性子鏈的微分方程組為

其中n是鏈包含的核素數(shù)目;Ni是鏈上第i個核素在任意t時刻的含量;γi是第i個核素生成第i+1個核素的轉換系數(shù)(γi≤λi).若(5)式中非齊次項滿足則方程組中第一個方程的非齊次項變成m+1項之和,根據(jù)線性常微分方程理論[12],方程組的解為1個齊次方程組的通解和另m+1個非齊次方程組的特解之和,第l個非齊次方程的非齊次項正是tlal.若通解為Nbase,特解分別為N(l)(l=0,···,m),則方程組的解為對第l個非齊次方程組進行Laplace變換,設變換后的向量為n(l)(p),其第K個分量滿足

再對(7)式進行Laplace逆變換,可得如下遞推關系:

(8)式說明第l個和第l?1個非齊次方程組的解存在遞推關系.綜合上述推導得到解方程組的流程:

1)求解齊次方程組,先對方程組的第一個方程Laplace變換,得到第一個變換分量,并利用遞推關系(6)式得到所有變換分量,再對所有變換分量進行Laplace逆變換,得到Nbase;

2)求解非齊次方程組,先解l=0的方程組,類似1)的流程,得到N(0),根據(jù)遞推關系(8)式,再對N(0)進行多次積分,可算出全部N(l);

3)將以上所有解相加,得到最終解.

經過一系列推導,可得到解的普通形式,它給出了一條線性鏈上任意一個核素在任意t＞0時刻的含量:

在以上推導中做了λi/=λj的假設,現(xiàn)實情況中,線性鏈上可能存在多對相同的核素,使鏈出現(xiàn)局部閉環(huán),導致假設不滿足.為防止計算式發(fā)散,必須進行處理.一種方式是考慮λi=λj的情況運用極限進行重新推導,得到一組新的解析表達式[2];另一種方式是對相等的λi做微小的擾動Δλ,人為確保λi互不相等.有研究表明[13],前者對計算精度略有提升,但會額外引入計算量,后者對計算精度影響不大.本文選擇相對更容易實現(xiàn)的后者.

2.2 非齊次燃耗方程基于有理近似方法的求解

有理近似方法的基本思想是用一個有理分式在rk(x)某個區(qū)域上逼近目標函數(shù)f(x),如

其中Pk,Qk是k階實多項式,無公因子,一般情況下令Qk常數(shù)項為1.(10)式也可以寫成如下形式:

矩陣函數(shù)同樣可以用關于矩陣的有理分式逼近,根據(jù)矩陣函數(shù)的定義可以將(11)式推廣到矩陣[15],即

其中I是單位矩陣.如果k是偶數(shù),則{θi},{αi}分別是k/2個共軛數(shù)對,(12)式的運算可以減半,變?yōu)?/p>

對于燃耗方程(4)式,若方程的初始條件為N(0)=N0,將其代入如下遞推式:

可求出N1,再用N1求出N2,一直進行下去,直到m→∞,得矩陣級數(shù)形式的解[16]

在實際燃耗問題中,有些核素的衰變非?？?致使消亡系數(shù)很大,有些問題中時間步長t也會取很大,使得‖tA‖＞1020,直接求解矩陣tA的上述級數(shù)展開不可承受.如果找到相應的有理逼近式,就能避開直接求解的困難.定義和函數(shù)Hl(x)滿足:如能找到Hl(x)的有理逼近式,根據(jù)(12)式,可將級數(shù)求和轉換成有限項計算:

圖1 Δ(t)=Hl?r14的圖像 (a)H?1的Δ(t);(b)H0的Δ(t);(c)H1的Δ(t);(d)H2的Δ(t)Fig.1.Curves of Δ(t)=Hl?r14:(a)Δ(t)of H?1;(b)Δ(t)of H0;(c)Δ(t)of H1;(d)Δ(t)of H2.

實函數(shù)一般具有惟一的最佳有理逼近[14],稱為最佳Chebyshev有理逼近,尋找方法主要有Remez方法[17],Carathéodory-Fejér方法[18?20](CF方法)等.Remez方法的基本思路是在區(qū)間[?1,1]上選定2k+2個Chebyshev多項式零點{cn},n=1,···,2k+2, 在這些點上 目標函數(shù)和有理逼近式的殘差以δ水平交替變換,即由于Pk有k+1個未知系數(shù),Qk有k個,連同δ,共2k+2個未知數(shù),正好和方程數(shù)相等,于是可求出所有系數(shù)和δ,如果δ過大則需要調整{cn}的位置,并繼續(xù)求解方程,直至δ收斂到期望水平.Remez方法是一種迭代過程,需要求解非線性方程組,計算量大,難收斂.CF方法的基本思路是在區(qū)間[?1,1]上用Chebyshev多項式展開目標函數(shù),用展開系數(shù)構造Hankel矩陣,Carathéodory-Fejér定理具體給出了目標函數(shù)和最佳有理逼近式的殘差的具體形式,它可以用Hankel矩陣的特征值和特征向量表示,進而通過殘差求出最佳有理逼近式.事實上,CF方法求出的僅是最佳逼近式的一個近似逼近,由于近似逼近向最佳逼近的收斂速度很快,可認為二者基本等效[14].CF方法的實現(xiàn)流程較為復雜,但不需要迭代,數(shù)值穩(wěn)定性好.CF方法一般在[?1,1]區(qū)間上求解目標函數(shù)的最佳有理逼近式,對于燃耗矩陣tA,有研究表明[6]一般情況下其特征值分布在復平面負x軸上或附近,因而需要得到目標函數(shù)在(?∞,0)上的近似式,用以計算tA的矩陣函數(shù),為此可進行變換x=進而對目標函數(shù)使用CF方法,再代回x即求出目標函數(shù)f(x)在(?∞,0)上的最佳逼近式.

表1 Hl(x)的14階有理逼近式系數(shù)Table 1.Rational approximation coefficients of order 14 of Hl(x).

本工作借鑒文獻[18]所述流程,用C語言和高精度浮點計算庫第三方庫實現(xiàn)了CF方法.目前只計算了Hl(x)在l=?1,0,1,2情況下的最佳有理逼近式,每個和函數(shù)的逼近階數(shù)k從2算到20.最佳有理逼近的精度按O(9.29?k)收斂[18],k取值過小無法保證精度,過大則增加計算量,k取14是相對折中的做法.表1列出了k=14時,Hl(x)對應的最佳有理逼近系數(shù),取11位有效數(shù)字(實際計算中算到了30位有效數(shù)字);圖1是和14階最佳有理逼近式在t∈[?1,1]內的絕對誤差Δ(t)的圖像,可看出誤差以α0水平振蕩;圖2是Hl(x)與最佳逼近式最大誤差隨逼近階數(shù)k的變化情況,可看出當l越大,誤差越小,說明越大的l對應的Hl(x)越容易收斂.

圖2 Hl-rk最大絕對誤差隨k的變化情況Fig.2.Maximum absolute error of Hl－rkvarying with order k.

2.3 方法適用范圍

以上方法都是基于非齊次項滿足f(t)=推導出的,在許多情況下,非齊次項需要用非常高階的泰勒展開逼近,甚至不能被有限階泰勒展開很好地逼近,這時使用上述方法計算量非常大,可能還會失效.不過,以上方法的推導過程為更一般情況的處理提供了一種思路,沿用上述方法流程,有可能求解非齊次項f(t)為其他形式時的情況.

對于有理近似方法,可以先尋找和燃耗方程(2)具有相同形式的代數(shù)方程y′=?λy+f(t)的一個特解h(t),只要h(t)在問題所關心的t范圍內有界、無奇異,對于任意t能明確給出h(t)的值,無論h(t)是何種形式,即便h(t)可能存在不可積的情形,都可以通過CF方法求出h(t)的一個最佳有理逼近式,從這一點講有理近似方法的通用性比TTA方法要好.但需要說明,對于不同h(t),達到期望的逼近精度所需要的有理逼近式階數(shù)也不同,有些h(t)需要非常高的逼近階數(shù),這會增大計算量.本文求出了f(t)=tl,l=0,1,···,m情形下的h(t),用14階逼近就可以達到很高的逼近精度,但對其他形式的f(t)未必夠用.

本文只從操作實現(xiàn)上論述方程求解的可能性,對于更多形式的f(t)是否可解,限于作者水平無法系統(tǒng)展開.對一些實際工程問題,如液態(tài)散裂靶的運行,由于束流轟擊散裂靶的強度隨時間變化并不劇烈,散裂產物的累積速率相對恒定,可近似認為非齊次項為常數(shù)項,本文所述方法能夠勝任求解.

3 數(shù)值實驗

將以上算法集成到了北京應用物理與計算數(shù)學研究所開發(fā)的燃耗計算程序JBURN中,通過兩個數(shù)值算例驗證程序功能的正確性.首先是一個小型矩陣問題,矩陣規(guī)模為6×6,非齊次項展開到2階.各系數(shù)和初始條件詳細如下:其中A的特征值全為負,最大特征值和最小特征值差距11個量級,且A能夠對角化. 計算t=105,106,107時刻的解,此情況下tA的范數(shù)約為109—1011.若A的特征值對角陣和特征向量陣為Λ和P,即A=PΛP?1,代入(15)式易得Hl(tA)=PHl(tΛ)P?1,根據(jù)矩陣函數(shù)定義[15],Hl(tΛ)也是對角陣,它的每個元素可由tΛ的每個元素代入Hl(x)中直接求出,以這種方法得到的結果作為基準,分別使用TTA和有理近似方法求解,結果列于表2.可以看出,兩種數(shù)值解法與參考解所得結果相差很小,有些完全一致,驗證了方法的正確性.

第二個算例基于真實的核素反應數(shù)據(jù)庫,計算MOX核材料經5.3075×1015cm?2·s?1大小的中子通量密度輻照20天后各種核素的含量,材料的具體成分和各種展開系數(shù)列于表3,其中正展開系數(shù)表示核素不斷地由外界補給,負展開系數(shù)表示核素向外界流出.

表2 小型矩陣問題在不同時刻的解Table 2.Solution of small matrix problem at different time.

表3 MOX燃耗算例計算初始條件Table 3.Initial conditions of MOX burnup example.

核反應數(shù)據(jù)庫中包含3400種核素,因此矩陣A的規(guī)模為3400×3400,tA的范數(shù)約為1030.用數(shù)值方法求解tA的特征值和特征向量耗時且不穩(wěn)定,因此無法像算例一一樣給出參考解,僅給出TTA方法和有理近似方法的計算結果,其中核素含量小于10?24時強制設為0.圖3給出了3400種核素的計算結果,并給出計算相對偏差,從圖中可以看出,絕大部分核素含量的計算結果符合得很好,偏差小于10%.進一步統(tǒng)計,3400種核素中有1251種核素含量不為0,其中26個核素的偏差大于10%.303個核素偏差在0.1%—10%之間,481個核素偏差在10?3%—0.1%之間,441個核素的偏差小于10?3%.表4列出了工程計算中幾種重要重金屬和裂變產物核素的含量對比情況,可以看出,對于重核而言,計算偏差非常小,偏差最大的243Am也只有5.91×10?2%;對于裂變產物而言偏差稍大,但偏差最大的151Sm也不超過1%,說明兩種方法的精度足夠滿足工程應用需求.

圖3 兩種方法計算MOX材料燃耗結果與相對偏差Fig.3.Results and relative error of MOX burnup calculated by two methods.

表4 重要重金屬與裂變產物核素含量Table 4.Quantity of important actinides and fission products.

在計算耗時方面,TTA方法耗時約85 s,有理近似方法僅耗時0.8 s,顯著快于TTA方法,主要原因是TTA方法需要通過回溯算法搜索構造線性子鏈,本問題中構造出300多萬條核素鏈,成為最耗時的因素.一般而言,對于考慮重金屬裂變的問題,線性子鏈個數(shù)會相當龐大,TTA方法耗時會顯著高于有理近似方法;對于純粹的核素衰變問題,只需要若干個線性子鏈就能描述,這時TTA方法則會顯著快于有理近似方法.有理近似方法耗時主要在解線性方程環(huán)節(jié),根據(jù)(12)式,方程的個數(shù)和有理逼近式階數(shù)成正比,對不同問題,在相同逼近階數(shù)下,有理近似方法求解速度沒有顯著波動.

4 結 論

本文初步研究了非齊次燃耗方程的兩種數(shù)值解法,用以求解方程非齊次項含時的情況.在非齊次項能夠被有限階泰勒展開逼近的條件下,首先基于Laplace變換推導出了方程基于線性子鏈方法的解析表達式,然后使用CF方法計算了矩陣級數(shù)解的近最佳有理逼近式.使用不同方法計算兩個數(shù)值算例,結果符合很好,驗證了方法的正確性與精度.本文的方法在計算精度和效率上能夠滿足工程計算要求.方法的實現(xiàn)流程也為非齊次燃耗方程的求解提供了一種思路,在后續(xù)工作中,可借鑒上述方法,求解含有其他形式非齊次項的燃耗方程.