本發(fā)明涉及核反應(yīng)堆堆芯核設(shè)計(jì)及安全領(lǐng)域，具體涉及一種壓水堆穩(wěn)態(tài)計(jì)算的多物理耦合方法。

背景技術(shù)：

由于核電行業(yè)的特殊性，要進(jìn)行反應(yīng)堆的安全測(cè)試或是燃料優(yōu)化管理的實(shí)驗(yàn)成本很高且不易操作，所以現(xiàn)階段我們通常采用軟件仿真的手段來(lái)模擬計(jì)算反應(yīng)堆的真實(shí)運(yùn)行。核電站是由不同功能的系統(tǒng)組成的，具有多元化物理過(guò)程的特點(diǎn)，各個(gè)物理過(guò)程相互影響，十分復(fù)雜，尤其是在核電站的反應(yīng)堆堆芯部分，更是集中地存在多種物理過(guò)程。例如中子會(huì)引發(fā)重核的裂變，這一過(guò)程釋放的能量使得堆芯內(nèi)冷卻劑溫度升高，溫度場(chǎng)的變化又會(huì)引起流場(chǎng)以及核數(shù)據(jù)的變化，核數(shù)據(jù)的改變又會(huì)反過(guò)來(lái)作用于中子的碰撞，影響中子裂變率，這就是各個(gè)物理場(chǎng)之間相互影響。為了正確的描述這一多物理過(guò)程在壓水堆核設(shè)計(jì)軟件中針對(duì)反應(yīng)堆多物理過(guò)程需要進(jìn)行耦合，這樣復(fù)雜的多物理過(guò)程也為壓水堆堆芯穩(wěn)態(tài)計(jì)算軟件中的多物理耦合方法提出了很高的要求。由于在反應(yīng)堆系統(tǒng)中每增加一個(gè)物理場(chǎng)就會(huì)使計(jì)算難度增加，所以這個(gè)多物理的耦合方法對(duì)于壓水堆堆芯軟件最后計(jì)算得到結(jié)果的精確性以及正確性是十分重要的，一種優(yōu)秀的壓水堆堆芯軟件的多物理耦合方法是可以提高計(jì)算精度以及軟件穩(wěn)定度的。

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

為了克服上述現(xiàn)有技術(shù)存在的問(wèn)題，本發(fā)明的目的在于提供一種壓水堆穩(wěn)態(tài)計(jì)算的多物理耦合方法，在原有耦合方法上做出進(jìn)一步改進(jìn)以提高軟件計(jì)算的穩(wěn)定性，并使得計(jì)算結(jié)果更為精確。

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

一種壓水堆穩(wěn)態(tài)計(jì)算的多物理該耦合方法，該耦合方法包括了截面計(jì)算、三維中子擴(kuò)散計(jì)算、熱工水力計(jì)算以及燃耗計(jì)算，具體步驟如下：

步驟1：根據(jù)工況點(diǎn)插值擬合得到節(jié)塊的群常數(shù)：根據(jù)組件計(jì)算給出的群常數(shù)-狀態(tài)參數(shù)的離散關(guān)系，利用最小二乘擬合方法進(jìn)行函數(shù)化和回代獲得堆芯計(jì)算所需的某一特定工況點(diǎn)下的群常數(shù)；采用壓水堆常用的組合方式：

其中，bu表示燃耗，cb為硼濃度，tf為冷卻劑溫度，tm為慢化劑溫度，δ(cr)表示控制棒對(duì)截面影響；fbase為關(guān)于燃耗和硼濃度的基態(tài)函數(shù)，fcrt為關(guān)于冷卻劑、慢化劑以及硼濃度的重啟項(xiàng)函數(shù)；燃耗深度取四階擬合；

需要通過(guò)擬合獲得的群常數(shù)如下：

g＝1,…,g--------------------不包括毒物的影響的宏觀總截面

g＝1,…,g-------------------不包括毒物的影響宏觀吸收截面

g＝1,…,g-----------宏觀散射矩陣

g＝1,…,g-----------------------------------宏觀裂變截面

g＝1,…,g-------------------------------宏觀中子產(chǎn)生截面

g＝1,…,g--------------------------宏觀中子能量產(chǎn)生截面

χ^g,g＝1,…,g-----------------------------宏觀中子裂變譜

yi,yxe,ypm,ysm----------------------毒物及其先驅(qū)核的裂變產(chǎn)額

----------考慮毒物及其先驅(qū)核的微觀吸收截面

步驟2：由上述得到的群常數(shù)進(jìn)行全堆三維中子擴(kuò)散計(jì)算：采用變分節(jié)塊法進(jìn)行擴(kuò)散計(jì)算，變分節(jié)塊法從二階偶對(duì)稱形式的中子輸運(yùn)方程出發(fā)，利用伽俐金變分技術(shù)在整個(gè)求解域上建立一個(gè)包含節(jié)塊中子平衡方程的泛函，再利用拉格朗日乘子法將邊界條件包含在該泛函內(nèi)，再通過(guò)里茲法以空間上的蒸餃多項(xiàng)式函數(shù)和角度上的球諧函數(shù)為基函數(shù)將這個(gè)泛函展開(kāi)，得到耦合了節(jié)塊內(nèi)中子通量密度與結(jié)塊邊界上的分中子流密度的節(jié)塊響應(yīng)矩陣，從而求解中子擴(kuò)散方程，得到三維全堆中子通量分布；

步驟3：由步驟2得到中子通量分布進(jìn)行熱工水力計(jì)算：已知堆芯內(nèi)中子通量分布能夠得到堆芯的功率分布；堆芯功率和中子通量密度按下式進(jìn)行歸一：

其中：vl為第l個(gè)節(jié)塊的體積；

vcore為堆芯的總裂變體積；

pav為堆芯的歸一化因子；

為第l節(jié)塊第g群的通量；

l為總節(jié)塊數(shù)；

g為總能群數(shù)；

為g能群第l節(jié)塊的宏觀中子能量產(chǎn)生截面；

熱工水力計(jì)算模型采用單通道模型；

步驟4：判斷核熱耦合是否收斂，若收斂，進(jìn)行下一步；否則返回步驟1；核熱耦合的收斂條件有兩個(gè)，分別是：

其中，是第m次迭代的有效增殖系數(shù)；

pi^m是第m次迭代第i個(gè)節(jié)塊內(nèi)的功率；

是兩次迭代計(jì)算的有效增殖系數(shù)的相對(duì)誤差；

是兩次迭代計(jì)算的節(jié)塊功率相對(duì)誤差最大值；

和是收斂誤差線值；

步驟5：燃耗計(jì)算：由步驟4可知在某工況點(diǎn)下堆芯的有效增殖系數(shù)以及堆芯功率，此時(shí)進(jìn)行燃耗計(jì)算，堆芯內(nèi)某處的燃耗深度由下式表示：

式中，表示空間i處在第k步時(shí)間點(diǎn)上的燃耗深度，mw·d/tu

pi^k表示空間i處在第k步時(shí)間點(diǎn)上功率，mw；

δt^k表示第k步燃耗段的時(shí)間步長(zhǎng)，d；

ρi表示空間i處的對(duì)應(yīng)于初始鈾裝載量，t；

步驟6：更新節(jié)塊宏觀截面：由步驟5中更新了燃耗，其他熱工參數(shù)不變，燃耗會(huì)改變截面值所以需要重新插值擬合截面，同步驟1；

步驟7：三維中子擴(kuò)散計(jì)算：由步驟6得到新的截面后再進(jìn)行三維中子擴(kuò)散計(jì)算，得到中子通量分布；

步驟8：根據(jù)用戶的選擇判斷是否進(jìn)行臨界硼搜索，若不進(jìn)行臨界硼搜索，進(jìn)行下一計(jì)算模塊；否則進(jìn)行臨界硼搜索并調(diào)整硼濃度至臨界并返回步驟1；由步驟7得到的中子通量計(jì)算有效增值因子判斷反應(yīng)堆是否臨界，若有效增殖因子值為1則反應(yīng)堆達(dá)到臨界狀態(tài)，進(jìn)入下一非多物理耦合的計(jì)算模塊；若有效增殖因子值不為1，則反應(yīng)堆尚未臨界，根據(jù)當(dāng)前的有效增殖因子調(diào)整硼濃度，將調(diào)整后的硼濃度帶回步驟1作為工況輸入?yún)?shù)之一去插值擬合截面。

本發(fā)明耦合方法與現(xiàn)有技術(shù)相比有如下突出特點(diǎn)：

將截面計(jì)算與燃耗計(jì)算脫耦，使得核熱耦合過(guò)程更易收斂，使軟件計(jì)算穩(wěn)定，減少迭代次數(shù)降低運(yùn)行時(shí)間。

附圖說(shuō)明

圖1為本發(fā)明多物理耦合流程圖。

具體實(shí)施方式

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

如圖1所示，本發(fā)明一種壓水堆穩(wěn)態(tài)計(jì)算的多物理該耦合方法，該耦合方法包括了截面計(jì)算、三維中子擴(kuò)散計(jì)算、熱工水力計(jì)算以及燃耗計(jì)算：具體步驟如下：

步驟1：根據(jù)工況點(diǎn)插值擬合得到節(jié)塊的群常數(shù)：根據(jù)組件計(jì)算給出的群常數(shù)-狀態(tài)參數(shù)的離散關(guān)系，利用最小二乘擬合方法進(jìn)行函數(shù)化和回代獲得堆芯計(jì)算所需的某一特定工況點(diǎn)下的群常數(shù)；采用壓水堆常用的組合方式：

其中，bu表示燃耗，cb為硼濃度，tf為冷卻劑溫度，tm為慢化劑溫度，δ(cr)表示控制棒對(duì)截面影響；fbase為關(guān)于燃耗和硼濃度的基態(tài)函數(shù)，fcrt為關(guān)于冷卻劑、慢化劑以及硼濃度的重啟項(xiàng)函數(shù)。燃耗深度取四階擬合；

需要通過(guò)擬合獲得的群常數(shù)如下：

g＝1,…,g--------------------不包括毒物的影響的宏觀總截面

g＝1,…,g-------------------不包括毒物的影響宏觀吸收截面

g＝1,…,g-----------宏觀散射矩陣

g＝1,…,g-----------------------------------宏觀裂變截面

g＝1,…,g-------------------------------宏觀中子產(chǎn)生截面

g＝1,…,g--------------------------宏觀中子能量產(chǎn)生截面

χ^g,g＝1,…,g-----------------------------宏觀中子裂變譜

yi,yxe,ypm,ysm----------------------毒物及其先驅(qū)核的裂變產(chǎn)額

----------考慮毒物及其先驅(qū)核的微觀吸收截面

步驟2：由上述得到的群常數(shù)進(jìn)行全堆三維中子擴(kuò)散計(jì)算：采用變分節(jié)塊法進(jìn)行擴(kuò)散計(jì)算，變分節(jié)塊法從二階偶對(duì)稱形式的中子輸運(yùn)方程出發(fā)，利用伽俐金變分技術(shù)在整個(gè)求解域上建立一個(gè)包含節(jié)塊中子平衡方程的泛函，再利用拉格朗日乘子法將邊界條件包含在該泛函內(nèi)，再通過(guò)里茲法以空間上的蒸餃多項(xiàng)式函數(shù)和角度上的球諧函數(shù)為基函數(shù)將這個(gè)泛函展開(kāi)，得到耦合了節(jié)塊內(nèi)中子通量密度與結(jié)塊邊界上的分中子流密度的節(jié)塊響應(yīng)矩陣，從而求解中子擴(kuò)散方程，得到三維全堆中子通量分布；

步驟3：由步驟2得到中子通量分布進(jìn)行熱工水力計(jì)算：已知堆芯內(nèi)中子通量分布能夠得到堆芯的功率分布；堆芯功率和中子通量密度按下式進(jìn)行歸一：

其中：vl為第l個(gè)節(jié)塊的體積；

vcore為堆芯的總裂變體積；

pav為堆芯的歸一化因子；

為第l節(jié)塊第g群的通量；

l為總節(jié)塊數(shù)；

g為總能群數(shù)；

為g能群第l節(jié)塊的宏觀中子能量產(chǎn)生截面；

熱工水力計(jì)算模型采用單通道模型；

步驟4：判斷核熱耦合是否收斂，若收斂，進(jìn)行下一步；否則返回步驟1；核熱耦合的收斂條件有兩個(gè)，分別是：

其中，是第m次迭代的有效增殖系數(shù)；

pi^m是第m次迭代第i個(gè)節(jié)塊內(nèi)的功率；

是兩次迭代計(jì)算的有效增殖系數(shù)的相對(duì)誤差；

是兩次迭代計(jì)算的節(jié)塊功率相對(duì)誤差最大值；

和是收斂誤差線值；

步驟5：燃耗計(jì)算：由步驟4可知在某工況點(diǎn)下堆芯的有效增殖系數(shù)以及堆芯功率，此時(shí)進(jìn)行燃耗計(jì)算，堆芯內(nèi)某處的燃耗深度由下式表示：

式中，表示空間i處在第k步時(shí)間點(diǎn)上的燃耗深度，mw·d/tu

pi^k表示空間i處在第k步時(shí)間點(diǎn)上功率，mw；

δt^k表示第k步燃耗段的時(shí)間步長(zhǎng)，d；

ρi表示空間i處的對(duì)應(yīng)于初始鈾裝載量，t；

步驟6：更新節(jié)塊宏觀截面：由步驟5中更新了燃耗，其他熱工參數(shù)不變，燃耗會(huì)改變截面值所以需要重新插值擬合截面，同步驟1；

步驟7：三維中子擴(kuò)散計(jì)算：由步驟6得到新的截面后再進(jìn)行三維中子擴(kuò)散計(jì)算，得到中子通量分布；

步驟8：根據(jù)用戶的選擇判斷是否進(jìn)行臨界硼搜索，若不進(jìn)行臨界硼搜索，進(jìn)行下一計(jì)算模塊；否則進(jìn)行臨界硼搜索并調(diào)整硼濃度至臨界并返回步驟1；由步驟7得到的中子通量計(jì)算有效增值因子判斷反應(yīng)堆是否臨界，若有效增殖因子值為1則反應(yīng)堆達(dá)到臨界狀態(tài)，進(jìn)入下一非多物理耦合的計(jì)算模塊；若有效增殖因子值不為1，則反應(yīng)堆尚未臨界，根據(jù)當(dāng)前的有效增殖因子調(diào)整硼濃度，將調(diào)整后的硼濃度帶回步驟1作為工況輸入?yún)?shù)之一去插值擬合截面。