本發(fā)明涉及核工程中燃耗數(shù)值模擬領(lǐng)域，具體涉及一種基于線積分有理近似方法求解微擾燃耗靈敏度系數(shù)的方法。
背景技術(shù)：
：反應(yīng)堆中能量的釋放來自于核素與中子發(fā)生反應(yīng)和放射性核素自發(fā)衰變反應(yīng)。而正是由于這些反應(yīng)，使得反應(yīng)堆的核素組成在不斷的發(fā)生變化。燃耗計(jì)算的目標(biāo)即在給定核素之間轉(zhuǎn)換關(guān)系以及相應(yīng)系數(shù)(包括中子反應(yīng)截面，中子通量，衰變常數(shù)，分支比等)，和初始核素原子核密度向量的情況下，得到核素原子核密度向量隨時間變化的特性。描述放射性核素自發(fā)衰變性質(zhì)的參數(shù)直接來源于核評價數(shù)據(jù)庫；核素中子反應(yīng)截面參數(shù)一般是首先由核評價庫加工得到點(diǎn)截面或者多群截面，再經(jīng)中子輸運(yùn)計(jì)算后按照通量能譜歸并為問題相關(guān)的單群截面；中子通量密度則是根據(jù)輸運(yùn)計(jì)算得到的各中子反應(yīng)單群截面和中子反應(yīng)的反應(yīng)熱，按照給定功率水平進(jìn)行歸一而得到。核評價數(shù)據(jù)庫中存儲的核素衰變性質(zhì)參數(shù)，中子反應(yīng)截面以及反應(yīng)熱來自于實(shí)驗(yàn)測量或者理論預(yù)測，均存在一定的不確定性，從而會導(dǎo)致燃耗計(jì)算結(jié)果的不確定性。定義響應(yīng)為依賴于計(jì)算結(jié)果函數(shù)式，對該響應(yīng)不確定性的評價可以通過計(jì)算不同輸入?yún)?shù)的靈敏度系數(shù)的方法而實(shí)現(xiàn)。經(jīng)過一系列理論推導(dǎo)，在僅考慮燃耗計(jì)算自身時，靈敏度系數(shù)的計(jì)算可歸結(jié)為關(guān)于燃耗方程、擾動參數(shù)和響應(yīng)的一個積分運(yùn)算；在考慮燃耗-輸運(yùn)耦合計(jì)算時，該積分計(jì)算也是靈敏度系數(shù)表達(dá)式中的一部分。燃耗方程的求解算法有泰勒展開與截?cái)?、線性子鏈法(TTA)、切比雪夫有理近似(CRAM)和線積分有理近似(QRAM)等。泰勒展開與截?cái)嗍禽^為經(jīng)典的方法，具有效率較高的優(yōu)勢，然而該方法應(yīng)用的前提是短壽期核素單獨(dú)予以處理以降低燃耗矩陣的范數(shù)，導(dǎo)致其求解精度較低。切比雪夫有理近似和線積分有理近似適于帶通量的燃耗方程求解，兩種方法的有效性基于燃耗矩陣特征值分布于負(fù)實(shí)軸附近的特點(diǎn)，具有優(yōu)異的求解精度和效率。線性子鏈法在衰變?nèi)己姆匠痰那蠼馍蟽?yōu)于其他方法。靈敏度系數(shù)計(jì)算所涉及的積分運(yùn)算，目前通常采用數(shù)值積分法，劃分計(jì)算時間區(qū)間為若干小區(qū)間，每個小區(qū)間進(jìn)行兩次燃耗計(jì)算，并存儲計(jì)算結(jié)果。為了達(dá)到計(jì)算精度需求，往往需要劃分上百個小區(qū)間，因而產(chǎn)生了巨大的時間和存儲開銷。技術(shù)實(shí)現(xiàn)要素：為了克服上述現(xiàn)有技術(shù)存在的問題，本發(fā)明的目的在于提供一種基于線積分有理近似方法求解微擾燃耗靈敏度系數(shù)的方法，線積分有理近似算法屬于燃耗方程求解的矩陣數(shù)值逼近類算法，具有較高的精度和效率；雖然其效率略低于切比雪夫有理近似方法，但是其求得的解精度隨時間變化的穩(wěn)定性更高；使用線積分有理近似方法可以展開原子核密度向量以及共軛原子核密度向量為求解時間區(qū)間上的向量函數(shù)式，進(jìn)而代入積分項(xiàng)中求得微擾燃耗靈敏度系數(shù)，預(yù)期在達(dá)到同樣計(jì)算精度的條件下該方法的效率和存儲需求均優(yōu)于數(shù)值積分方法。為了實(shí)現(xiàn)以上目的，本發(fā)明采取如下的技術(shù)方案：基于線積分有理近似方法求解微擾燃耗靈敏度系數(shù)的方法，包括如下步驟：步驟1：確定線積分有理近似的離散化參數(shù)：選取積分曲線Γ為公式(1)確定的拋物型曲線，根據(jù)燃耗矩陣特征值分布于負(fù)實(shí)軸附近的特性，確保所有特征值被該積分曲線圍繞一次；z＝φ(u)＝μ(iu+1)2,u∈R公式(1)式中：z——復(fù)數(shù)變量u——實(shí)數(shù)變量φ(u)——映射實(shí)軸為復(fù)平面上拋物型積分曲線的函數(shù)，簡稱映射函數(shù)μ——控制曲線尺寸大小的因子選取離散積分點(diǎn)以及權(quán)重因子如下：zk=φ^(uk)=μ^(i·(k-n+12)h^+1)2k=1,...,n]]>式中：zk，ck——離散積分點(diǎn)以及權(quán)重因子——優(yōu)化后的映射函數(shù)及其導(dǎo)函數(shù)——優(yōu)化后的控制離散點(diǎn)選取的量i——虛數(shù)單位n——線積分有理近似階數(shù)uk——離散積分點(diǎn)對應(yīng)的實(shí)數(shù)μ^=π·n4t18Λ+1]]>h^=8Λ+1n]]>式中：t1——控制解穩(wěn)定區(qū)間的量Λ——控制解穩(wěn)定區(qū)間的量關(guān)于控制解穩(wěn)定區(qū)間的量t1和Λ，其值的選取在后面的步驟3有進(jìn)一步討論；步驟2：選取有理近似階數(shù)n為偶數(shù)，計(jì)算原子核密度向量函數(shù)以及共軛原子核密度向量函數(shù)，其中，原子核密度向量函數(shù)如下式：式中：N——原子核密度向量t——燃耗時間tf——燃耗步末時間I——單位矩陣A——燃耗方程系數(shù)矩陣，簡稱燃耗矩陣共軛原子核密度向量函數(shù)如下式：式中：N*——共軛原子核密度向量A*——燃耗矩陣的共軛矩陣其中共軛原子核密度向量N*(tf)依賴于響應(yīng)R的定義。公式(4)可寫做：式中：Pk——Re((zkI-tfA)-1N(0))，即式(4)第k個解分量的實(shí)部Qk——Im((zkI-tfA)-1N(0))，即式(4)第k個解分量的虛部xk，yk——積分點(diǎn)zk的實(shí)部與虛部ak，bk——權(quán)重因子ck的實(shí)部與虛部τ——t/tf，稱為相對時間公式(5)可寫做：式中：Pk*——Re((zkI-tfA*)-1N*(tf))，式(5)第k個解分量的實(shí)部Qk*——Im((zkI-tfA*)-1N*(tf))，式(5)第k個解分量的虛部步驟3：將原子核密度向量函數(shù)以及共軛原子核密度向量函數(shù)代入求解微擾燃耗靈敏度系數(shù)的積分公式中，進(jìn)行解析積分，得到微擾燃耗靈敏度系數(shù)；微擾燃耗靈敏度系數(shù)由下面的積分公式計(jì)算得到：式中：R——響應(yīng)——微擾燃耗靈敏度系數(shù)B——燃耗矩陣的擾動方向矩陣將公式(6)和公式(7)代入上式，使用解析的方式得到微擾燃耗靈敏度系數(shù)；當(dāng)相對時間τ接近0時，原子核密度向量函數(shù)的精度較差；而相對時間τ接近1時，共軛原子核密度向量函數(shù)的精度較差；為了降低靈敏度系數(shù)的求解誤差，采用相對時間分段的方法。以相對時間分三段為例，分別得到在相對時間區(qū)間[0.0,0.01]，[0.01,0.1]和[0.1,1.0]上的原子核密度向量函數(shù)，以及相對時間區(qū)間[0.99,1.0]，[0.9,0.99]和[0,0.9]上的共軛原子核密度向量函數(shù)。在三段相對時間區(qū)間上，控制解穩(wěn)定區(qū)間的量Λ取為3，t1分別取為0.01，0.1以及1.0。公式(8)分五段積分，即[0.0,0.01]，[0.01,0.1]，[0.1,0.9]，[0.9,0.99]和[0.99,1.0]。上面的分段策略可以保證靈敏度系數(shù)計(jì)算的相對誤差在千分位。對于不同的求解精度和效率要求，可以采取不同的分段數(shù)目以及相應(yīng)的調(diào)整控制穩(wěn)定區(qū)間的量的取值。和現(xiàn)有技術(shù)相比較，本發(fā)明具備如下優(yōu)點(diǎn)：對于求解微擾燃耗靈敏度系數(shù)，本發(fā)明所提出的方法能夠在保證精度需求的前提條件下，相比數(shù)值積分方法獲得更高的效率和節(jié)省內(nèi)存使用。以上述相對時間三段劃分為例，本發(fā)明提出的方法的計(jì)算量相當(dāng)于12次燃耗計(jì)算，存儲量相當(dāng)于300個原子核密度向量(按照有理近似階數(shù)為24估計(jì))。在達(dá)到相同精度的條件下，預(yù)計(jì)數(shù)值積分方法至少需要上千次燃耗計(jì)算的計(jì)算量以及相對應(yīng)數(shù)目的原子核密度向量的存儲量。因此本發(fā)明提出的方法可以作為微擾燃耗敏感性分析的有力工具。附圖說明圖1QRAM和CRAM對指數(shù)函數(shù)的逼近誤差比較。具體實(shí)施方式下面結(jié)合附圖對本發(fā)明做進(jìn)一步詳細(xì)說明：本發(fā)明使用線積分有理近似方法展開原子核密度向量函數(shù)與共軛原子核密度向量函數(shù)，以此為基礎(chǔ)得到微擾燃耗靈敏度系數(shù)。如圖1所示，線積分有理近似方法得到的解(QRAM)，其精度隨時間變化的穩(wěn)定性高于切比雪夫有理近似方法(CRAM)。具體包括如下步驟：步驟1：確定線積分有理近似的離散化參數(shù)。選取積分曲線Γ為公式(1)確定的拋物型曲線，根據(jù)燃耗矩陣特征值分布于負(fù)實(shí)軸附近的特性，可以確保所有特征值被該積分曲線圍繞一次。因此可以應(yīng)用矩陣指數(shù)基于復(fù)平面上線積分的定義式。z＝φ(u)＝μ(iu+1)2,u∈R公式(1)式中：z——復(fù)數(shù)變量u——實(shí)數(shù)變量φ(u)——映射實(shí)軸為復(fù)平面上拋物型積分曲線的函數(shù)，簡稱映射函數(shù)μ——控制曲線尺寸大小的因子線積分的計(jì)算借助于數(shù)值積分。根據(jù)平衡截?cái)嗾`差與數(shù)值積分誤差的原則，選取離散積分點(diǎn)以及權(quán)重因子如下：zk=φ^(uk)=μ^(i·(k-n+12)h^+1)2k=1,...,n]]>式中：zk，ck——離散積分點(diǎn)以及權(quán)重因子——優(yōu)化后的映射函數(shù)及其導(dǎo)函數(shù)——控制離散點(diǎn)選取的量，后面有表達(dá)式i——虛數(shù)單位n——線積分有理近似階數(shù)uk——離散積分點(diǎn)對應(yīng)的實(shí)數(shù)μ^=π·n4t18Λ+1]]>h^=8Λ+1n]]>式中：t1——控制解穩(wěn)定區(qū)間的量Λ——控制解穩(wěn)定性的量關(guān)于控制解穩(wěn)定區(qū)間的量t1和Λ，其值的選取在后面的步驟3有進(jìn)一步討論。步驟2：計(jì)算原子核密度向量函數(shù)以及共軛原子核密度向量函數(shù)。由于采用了數(shù)值積分方法計(jì)算線積分，存在近似，所以使用約等號。其中，原子核密度向量函數(shù)如下式：式中：N——原子核密度向量t——燃耗時間tf——燃耗步末時間I——單位矩陣A——燃耗方程系數(shù)矩陣，簡稱燃耗矩陣有理近似階數(shù)n一般取為偶數(shù)，此時離散積分點(diǎn)以及權(quán)重因子呈共軛復(fù)數(shù)對的形式出現(xiàn)，因此上式可以簡化為：類似的，共軛原子核密度向量函數(shù)如下式：式中：N*——共軛原子核密度向量A*——共軛燃耗矩陣其中共軛原子核密度向量N*(tf)依賴于響應(yīng)R的定義。通過求解一系列復(fù)數(shù)線性方程組，公式(4)可寫做：式中：Pk——Re((zkI-tfA)-1N(0))，即式(4)第k個解分量的實(shí)部Qk——Im((zkI-tfA)-1N(0))，即式(4)第k個解分量的虛部xk，yk——積分點(diǎn)zk的實(shí)部與虛部ak，bk——權(quán)重因子ck的實(shí)部與虛部τ——t/tf，稱為相對時間公式(5)可寫做：式中：Pk*——Re((zkI-tfA*)-1N*(tf))，式(5)第k個解分量的實(shí)部Qk*——Im((zkI-tfA*)-1N*(tf))，式(5)第k個解分量的虛部步驟3：將原子核密度向量函數(shù)以及共軛原子核密度向量函數(shù)代入求解微擾燃耗靈敏度系數(shù)的積分公式中，進(jìn)行解析積分，得到微擾燃耗靈敏度系數(shù)。微擾燃耗靈敏度系數(shù)可由下面的積分公式計(jì)算得到：式中：R——響應(yīng)——微擾燃耗靈敏度系數(shù)B——燃耗矩陣的擾動矩陣將公式(6)和公式(7)代入上式，即可使用解析的方式得到微擾燃耗靈敏度系數(shù)。當(dāng)相對時間τ接近0時，原子核密度向量函數(shù)的精度較差；而相對時間τ接近1時，共軛原子核密度向量函數(shù)的精度較差。為了降低靈敏度系數(shù)的求解誤差，采用相對時間分段的方法。以相對時間分三段為例，分別得到在相對時間區(qū)間[0.0,0.01]，[0.01,0.1]和[0.1,1.0]上的原子核密度向量函數(shù)，以及相對時間區(qū)間[0.99,1.0]，[0.9,0.99]和[0,0.9]上的共軛原子核密度向量函數(shù)，在三段相對時間區(qū)間上，控制解穩(wěn)定區(qū)間的量Λ取為3，t1分別取為0.01，0.1以及1.0。公式(8)分五段積分，即[0.0,0.01]，[0.01,0.1]，[0.1,0.9]，[0.9,0.99]和[0.99,1.0]。如下式：式中：na(τ)，nb(τ)，nc(τ)——t1分別取0.01，0.1和1.0時得到的以相對時間為因變量的原子核密度向量函數(shù)na*(τ)，nb*(τ)，nc*(τ)——t1分別取0.01，0.1和1.0時得到的以相對時間為因變量的共軛原子核密度向量函數(shù)對于式(9)的計(jì)算，以相對時間區(qū)間[0,0.01]上的積分為例：式中：B(i,j)——矩陣B的第i行第j列元素計(jì)算矩陣Ca：Ca(2i-1,2j-1)=∫00.01exc,i(1-τ)+xa,jτ(ac,icos(yc,i(1-τ))-bc,isin(yc,i(1-τ)))×(aa,jcos(ya,jτ)-ba,jsin(ya,jτ))dτ]]>Ca(2i-1,2j)=-∫00.01exc,i(1-τ)+xa,jτ(ac,icos(yc,i(1-τ))-bc,isin(yc,i(1-τ)))×(ba,jcos(ya,jτ)+aa,jsin(ya,jτ))dτ]]>Ca(2i,2j-1)=-∫00.01exc,i(1-τ)+xa,jτ(bc,icos(yc,i(1-τ))+ac,isin(yc,i(1-τ)))×(aa,jcos(ya,jτ)-ba,jsin(ya,jτ))dτ]]>式(10)可轉(zhuǎn)化為：式中：M*c(i,:)——矩陣M*c的第i行Ma(j,:)——矩陣Ma的第j行上面的分段策略可以保證靈敏度系數(shù)計(jì)算的相對誤差在千分位量級。對于不同的求解精度和效率要求，可以采取不同的分段數(shù)目以及相應(yīng)的調(diào)整控制穩(wěn)定區(qū)間的量Λ和t1的取值。當(dāng)前第1頁1 2 3