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

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

1 引 言

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

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

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

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

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

其中al是展開(kāi)項(xiàng)向量,則非齊次燃耗方程變?yōu)?/p>

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

2 理論推導(dǎo)

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

TTA方法的基本思想是將核素之間的復(fù)雜轉(zhuǎn)換關(guān)系解耦成若干線性鏈,在一條線性鏈上,一個(gè)核素只能經(jīng)由一個(gè)母核反應(yīng)生成,且該核素發(fā)生反應(yīng)只能生成一個(gè)子核,只有鏈?zhǔn)缀怂鼐哂羞w移率和初始核子密度,其余核素都為0.這種簡(jiǎn)單的傳遞關(guān)系和初始條件使得描述一條線性子鏈的微分方程組可解析求解.逐一求解每一條線性子鏈,再把每條鏈的貢獻(xiàn)累加即得到結(jié)果.描述一條線性子鏈的微分方程組為

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

再對(duì)(7)式進(jìn)行Laplace逆變換,可得如下遞推關(guān)系:

(8)式說(shuō)明第l個(gè)和第l?1個(gè)非齊次方程組的解存在遞推關(guān)系.綜合上述推導(dǎo)得到解方程組的流程:

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

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

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

經(jīng)過(guò)一系列推導(dǎo),可得到解的普通形式,它給出了一條線性鏈上任意一個(gè)核素在任意t＞0時(shí)刻的含量:

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 方法適用范圍

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

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

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

3 數(shù)值實(shí)驗(yàn)

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

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

表2 小型矩陣問(wèn)題在不同時(shí)刻的解Table 2.Solution of small matrix problem at different time.

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

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

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

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

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

4 結(jié) 論

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