本發(fā)明涉及高保真先進(jìn)核反應(yīng)堆多物理場耦合模擬和分析，具體涉及一種液態(tài)金屬冷卻快堆pin-by-pin級穩(wěn)態(tài)物理熱工耦合分析方法。

背景技術(shù)：

1、在液態(tài)金屬冷卻快堆運(yùn)行中，堆芯中子物理與熱工水力具有強(qiáng)烈耦合關(guān)系，對核反應(yīng)堆設(shè)計(jì)和安全運(yùn)行具有關(guān)鍵影響。在液態(tài)金屬冷卻快堆中，還具有與壓水堆不同的特點(diǎn)。第一，快中子能譜變化范圍廣，且液態(tài)金屬良好的熱工性能使反應(yīng)堆出入口溫差更大，導(dǎo)致核反應(yīng)截面變化更加明顯。第二，快中子具有較長自由程，這導(dǎo)致反應(yīng)堆局部功率變化可能影響多個組件功率分布。第三，液態(tài)金屬冷卻快堆采用池式或六邊形組件特殊結(jié)構(gòu)，存在熱分層、橫流、回流、倒流等三維熱工水力現(xiàn)象。正是由于以上原因，開發(fā)高保真三維pin-by-pin級物理熱工耦合分析程序?qū)τ谏钊肜斫馕锢頍峁は嗷プ饔?，從而?yōu)化反應(yīng)堆設(shè)計(jì)，提高核反應(yīng)堆安全性至關(guān)重要。

2、現(xiàn)有液態(tài)金屬冷卻快堆pin-by-pin級物理熱工耦合分析程序中子物理多采用點(diǎn)堆或擴(kuò)散方程描述，熱工水力多采用子通道方法，不能體現(xiàn)中子在空間中的散射，以及冷卻劑橫流、回流、倒流等關(guān)鍵三維現(xiàn)象。隨著計(jì)算機(jī)發(fā)展，未來研究將需要集成更高精度的中子物理模型、更復(fù)雜的熱工水力計(jì)算以及全面三維cfd模擬，以全面提升快堆設(shè)計(jì)精度和運(yùn)行安全性。

技術(shù)實(shí)現(xiàn)思路

1、為了克服上述現(xiàn)有技術(shù)存在問題，本發(fā)明的目的在于提供一種液態(tài)金屬冷卻快堆pin-by-pin級穩(wěn)態(tài)物理熱工耦合分析方法，本發(fā)明基于統(tǒng)一平臺的有限體積框架openfoam，實(shí)現(xiàn)中子物理計(jì)算和熱工水力計(jì)算在內(nèi)存中共享數(shù)據(jù)。本發(fā)明中子物理部分采用三維中子輸運(yùn)方程描述物理現(xiàn)象，熱工水力部分采用三維cfd控制方程描述物理現(xiàn)象，能夠高保真地體現(xiàn)物理作用行為。本發(fā)明整體求解中，少群常數(shù)采用openmc精細(xì)建模生成常數(shù)庫并插值調(diào)用；采用相同網(wǎng)格(即子通道網(wǎng)格劃分)以簡化計(jì)算，不同模塊之間通過物理場映射和邊界條件映射傳遞參數(shù)；采用picard迭代耦合。本發(fā)明可以為分析液態(tài)金屬冷卻快堆的物理熱工耦合特性提供一種快速高保真pin-by-pin級的耦合分析方法，為液態(tài)金屬冷卻快堆的設(shè)計(jì)與安全特性分析提供建議與指導(dǎo)。

2、為了達(dá)到上述目的，本發(fā)明采用如下技術(shù)方案：

3、一種液態(tài)金屬冷卻快堆pin-by-pin級穩(wěn)態(tài)物理熱工耦合分析方法，包括如下步驟：

4、步驟1：采用openmc對液態(tài)金屬冷卻快堆的結(jié)構(gòu)和材料進(jìn)行pin-by-pin精細(xì)建模，并按照堆芯實(shí)際材料和結(jié)構(gòu)劃分均勻化的群常數(shù)提取區(qū)域，在不同的全局溫度和冷卻劑密度條件下通過蒙特卡洛計(jì)算提取少群常數(shù)，包括中子份額χg、能量釋放裂變截面(κσf)g、總截面即輸運(yùn)修正截面有效裂變中子產(chǎn)額截面(νσf)g和散射截面σs,g'→g，通過后處理制作為openfoam可讀取的少群常數(shù)庫；

5、步驟2：采用幾何建模軟件與網(wǎng)格劃分軟件，對液態(tài)金屬冷卻快堆進(jìn)行幾何建模與網(wǎng)格劃分，其中組件區(qū)域劃分子通道形式的網(wǎng)格，結(jié)構(gòu)材料區(qū)域根據(jù)需求自由劃分網(wǎng)格；根據(jù)openmc計(jì)算中均勻化的群常數(shù)提取區(qū)域進(jìn)行命名；

6、步驟3：設(shè)置初始物理場和邊界條件：初始物理場包括中子通量場分布和中子通量密度二階矩φ20場分布，冷卻劑溫度t0場分布、冷卻劑密度ρ0場分布和冷卻劑速度u0場分布；物理計(jì)算外邊界設(shè)置真空邊界，熱工計(jì)算入口邊界設(shè)置冷卻劑質(zhì)量流量，其余均按照默認(rèn)設(shè)定；

7、步驟4：進(jìn)入物理熱工耦合外迭代，根據(jù)少群常數(shù)的區(qū)域和網(wǎng)格中初始或上次迭代熱工計(jì)算的作為全局溫度的冷卻劑溫度場分布和冷卻劑密度場分布，進(jìn)行插值計(jì)算，更新少群常數(shù)；

8、步驟5：進(jìn)入中子物理部分內(nèi)迭代，采用簡化球諧函數(shù)sp3方法求解三維穩(wěn)態(tài)中子輸運(yùn)方程；具體如下：

9、步驟5-1：計(jì)算上一迭代步散射源項(xiàng)sscatter和裂變源項(xiàng)sfission；

10、步驟5-2：求解采用簡化球諧函數(shù)sp3方法近似后的三維穩(wěn)態(tài)中子輸運(yùn)方程；

11、步驟5-3：計(jì)算有效增殖系數(shù)keff；

12、步驟5-4：判斷中子通量密度和有效增殖系數(shù)keff是否收斂；若收斂則結(jié)束中子物理內(nèi)迭代，進(jìn)入步驟6；若不收斂，則返回步驟5-1；

13、步驟6：根據(jù)物理計(jì)算得到的中子通量場分布和液態(tài)金屬冷卻快堆設(shè)計(jì)基準(zhǔn)功率，計(jì)算得到液態(tài)金屬冷卻快堆實(shí)際功率場分布，作為熱工水力計(jì)算的熱量源項(xiàng)進(jìn)行更新；

14、步驟7：根據(jù)上次計(jì)算得到的冷卻劑速度場分布，冷卻劑密度場分布更新子通道網(wǎng)格信息；

15、步驟7-1：首先通過微元體相鄰面數(shù)量區(qū)分內(nèi)通道，然后通過微元體所有面法線向量相加的計(jì)算結(jié)果值區(qū)分邊通道和角通道，避免標(biāo)記通道帶來的錯誤，提高計(jì)算效率；

16、步驟7-2：根據(jù)上次熱工水力計(jì)算得到的速度場，更新信息包括流通面積a、潤濕周長p、水力直徑dh、雷諾數(shù)re；這些信息將通過與阻力模型、液態(tài)金屬傳熱模型、子通道形狀因子變換模型和湍流交混模型計(jì)算得到阻力二階張量k，湍流粘性系數(shù)μt，湍流導(dǎo)熱系數(shù)κt，子通道網(wǎng)格形狀因子a；

17、步驟8：進(jìn)入熱工水力部分內(nèi)迭代，采用simple算法求解三維cfd控制方程；

18、三維cfd控制方程形式如下：

19、質(zhì)量守恒方程：

20、

21、動量守恒方程：

22、

23、能量守恒方程：

24、

25、其中：ρ是冷卻劑密度，u是冷卻劑速度，p是冷卻劑壓力，g是重力系數(shù)，h是冷卻劑比焓，q是熱量源項(xiàng)；

26、simple算法將速度場與壓力場解耦進(jìn)行迭代求解，加入能量守恒方程后，內(nèi)迭代方程求解的順序依次是，能量守恒方程，速度預(yù)測方程，壓力泊松方程，速度修正方程；當(dāng)速度場，壓力場和溫度場收斂時結(jié)束熱工部分內(nèi)迭代，得到冷卻劑溫度場分布、冷卻劑速度場分布、冷卻劑密度場分布、冷卻劑壓力場場分布；

27、步驟9：判斷物理熱工耦合外迭代是否收斂，若兩次迭代計(jì)算溫度場最大值小于0.01k，有效增殖系數(shù)keff小于1e-5，則結(jié)束物理-熱工耦合計(jì)算；若不收斂，則返回步驟4；

28、步驟10：進(jìn)行物理熱工耦合特性分析，包括分析確定物理熱工耦合研究的必要性，優(yōu)化液態(tài)金屬冷卻快堆設(shè)計(jì)，提高液態(tài)金屬冷卻快堆安全性；分析液態(tài)金屬冷卻快堆的燃料多普勒溫度反饋和冷卻劑密度反饋現(xiàn)象；熱工參數(shù)如何通過影響堆芯及結(jié)構(gòu)材料的截面，從而影響中子行為的三維中子物理效應(yīng)；中子物理如何影響熱工參數(shù)的分布，并分析物理反饋?zhàn)饔孟乱簯B(tài)金屬冷卻快堆中橫流、倒流、回流三維熱工水力現(xiàn)象。

29、所述步驟1具體如下：

30、步驟1-1：首先，在openmc中進(jìn)行液態(tài)金屬冷卻快堆幾何結(jié)構(gòu)與材料的pin-by-pin精細(xì)建模，包括燃料、氣隙、包殼、冷卻劑、反射層、堆芯結(jié)構(gòu)支撐材料；

31、步驟1-2：將低富集度燃料區(qū)域、高富集度燃料區(qū)域、軸向反射層、徑向反射層、控制組件區(qū)域、外部支撐結(jié)構(gòu)材料區(qū)域、冷卻劑下降段區(qū)域進(jìn)行劃分和均勻化；

32、步驟1-3：設(shè)定能群結(jié)構(gòu)劃分和數(shù)量，設(shè)定蒙特卡洛計(jì)算粒子迭代次數(shù)，設(shè)定不同的全局溫度和冷卻劑密度進(jìn)行蒙特卡洛計(jì)算和均勻化的少群常數(shù)提??；

33、步驟1-4：openmc生成少群常數(shù)庫后，根據(jù)不同的全局溫度和冷卻劑密度后處理為openfoam可讀的少群常數(shù)庫。

34、所述步驟2具體如下：

35、步驟2-1：采用幾何建模軟件建立幾何模型；

36、步驟2-2：組件區(qū)域采用子通道棱柱的網(wǎng)格劃分方式，燃料棒中心位于棱柱交點(diǎn)處；徑向反射層，結(jié)構(gòu)材料及冷卻劑下降段采用任意自由的網(wǎng)格劃分方式；

37、步驟2-3：按照openmc均勻化提取少群常數(shù)的區(qū)域，為幾何體劃分同樣區(qū)域并命名，步驟4中將根據(jù)群常數(shù)的區(qū)域和網(wǎng)格中的全局溫度與冷卻劑密度進(jìn)行群常數(shù)插值計(jì)算。

38、所述步驟4中

39、插值調(diào)用的公式如下：

40、

41、其中，影響因素包含全局溫度和冷卻劑密度，基準(zhǔn)點(diǎn)為全局溫度t0和冷卻劑密度ρ0，溫度變化點(diǎn)為全局溫度t1和冷卻劑密度ρ0，密度變化點(diǎn)為全局溫度t0和冷卻劑密度ρ1。

42、所述步驟5具體如下：

43、步驟5-1：裂變源項(xiàng)sfission與散射源項(xiàng)sscatter的計(jì)算公式為：

44、

45、sscatter＝σs,g'→g,iφ0,g',i

46、式中χg,i是網(wǎng)格i中能群g的中子份額，(νσf)g',i是網(wǎng)格i中能群g’的有效裂變中子產(chǎn)額截面，φ0,g',i是網(wǎng)格i中能群g’的中子通量密度，∑s,g'→g,i是網(wǎng)格i中由能群g’散射到能群g的散射截面；

47、步驟5-2：采用簡化球諧函數(shù)sp3方法簡化后，穩(wěn)態(tài)的三維中子輸運(yùn)方程寫為：

48、

49、其中，dg,i＝1/(3σt,g,i)，σr,g,i＝σt,g,i-σs,g→g,i

50、式中，dg,i是網(wǎng)格i中能群g的擴(kuò)散系數(shù)，σt,g,i是網(wǎng)格i中能群g的輸運(yùn)截面，φ0,g,i是網(wǎng)格i中能群g中子通量密度，φ2,g,i是網(wǎng)格i中能群g中子通量密度二階矩，σr,g,i是網(wǎng)格i中能群g的移出截面，σs,g→g,i是網(wǎng)格i中由能群g散射到能群g的散射截面；

51、邊界條件采用近似簡化mashark邊界條件：

52、

53、其中φ*i＝φ0,i+2φ2,i，φ*i是網(wǎng)格i中的伴隨中子通量密度，φ0,i是網(wǎng)格i中的中子通量密度，φ2,i是網(wǎng)格i中的中子通量密度二階矩，αi是進(jìn)入控制體中子流與離開控制體中子流之比，當(dāng)αi＝0時為真空邊界條件，當(dāng)αi＝1時，為全反射邊界條件；

54、步驟5-3：計(jì)算有效增殖系數(shù)keff，為當(dāng)前迭代步與上一迭代步的中子通量密度比；

55、步驟5-4：判斷中子通量密度和有效增殖系數(shù)keff是否收斂；若兩次迭代中子通量密度小于1e-4且兩次迭代有效增殖系數(shù)keff值小于1e-5，則進(jìn)入步驟6；若不收斂，則返回步驟5-1；

56、所述步驟6具體如下：

57、步驟6-1：液態(tài)金屬冷卻快堆實(shí)際功率分布計(jì)算公式如下：

58、

59、式中，pi是液態(tài)金屬冷卻快堆實(shí)際功率分布，p0是液態(tài)金屬冷卻快堆設(shè)計(jì)基準(zhǔn)功率，(νσf)g',i是網(wǎng)格i中能群g的有效裂變中子產(chǎn)額截面，vi是網(wǎng)格i的體積；

60、步驟6-2：熱工水力計(jì)算的熱量源項(xiàng)計(jì)算公式如下：

61、

62、式中，qv是熱量源項(xiàng)，vi,fluid是網(wǎng)格i中除去固體燃料后實(shí)際冷卻劑流體體積。

63、所述步驟7具體如下：

64、步驟7-1：首先通過微元體相鄰面數(shù)量區(qū)分內(nèi)通道，內(nèi)通道控制體的相鄰面?zhèn)€數(shù)為5，而邊通道和角通道的控制體相鄰面?zhèn)€數(shù)為6；然后通過微元體所有面法線向量相加的計(jì)算結(jié)果值區(qū)分邊通道和角通道，邊通道的面法線向量和為0，角通道的面法線向量和大于0，在正六邊形組件中，角通道的面法線向量和為0.732，能夠避免標(biāo)記通道帶來的錯誤，提高計(jì)算效率；

65、步驟7-2：根據(jù)上次熱工水力計(jì)算得到的速度場，更新信息包括流通面積a、潤濕周長p、水力直徑dh、雷諾數(shù)re；這些信息將通過與阻力模型、液態(tài)金屬傳熱模型、子通道形狀因子變換模型和湍流交混模型計(jì)算得到阻力二階張量k，湍流粘性系數(shù)μt，湍流導(dǎo)熱系數(shù)κt，子通道網(wǎng)格形狀因子a；

66、由幾何關(guān)系計(jì)算不同類型通道中，流通面積的計(jì)算公式如下：

67、

68、ab＝n1·abi+n2·abe+n3·abc

69、式中，ppin是燃料棒心距，dpin是燃料棒直徑，pedge是燃料棒直徑與邊通道寬度之和，abi是內(nèi)通道的流通面積，abe是邊通道的流通面積，abc是角通道的流通面積，ab是總流通面積，n1是內(nèi)通道個數(shù)，n2是邊通道個數(shù)，n3是角通道個數(shù)；

70、不同類型通道中，潤濕周長的計(jì)算公式如下：

71、

72、pb＝n1×pbi+n2×pbe+n3×pbc

73、式中，pbi是內(nèi)通道的潤濕周長，pbe是邊通道的潤濕周長，pbc是角通道的潤濕周長，pb是總潤濕周長；

74、不同類型通道中，水力直徑的計(jì)算公式如下：

75、

76、式中，dhbi是內(nèi)通道的水力直徑，dhbe是邊通道的水力直徑，dhbc是角通道的水力直徑，dhb是平均水力直徑；

77、不同類型通道中，雷諾數(shù)的計(jì)算公式如下：

78、

79、式中，rebi是內(nèi)通道的雷諾數(shù)，rebe是邊通道的雷諾數(shù)，rebc是角通道的雷諾數(shù)，reb是平均雷諾數(shù)，μ是冷卻劑的粘度系數(shù)；

80、通過修改網(wǎng)格中以上的特征量信息能夠?qū)崿F(xiàn)子通道網(wǎng)格拓?fù)?，從而保證了三維cfd控制方程計(jì)算的正確性。

81、所述步驟8具體如下：

82、根據(jù)simple算法原理，將質(zhì)量守恒方程與動量守恒方程求解問題轉(zhuǎn)化為求解速度預(yù)測方程、壓力泊松方程和速度修正方程；

83、步驟8-1：求解能量守恒方程；

84、步驟8-2：求解速度預(yù)測方程；

85、步驟8-3：求解壓力泊松方程；

86、步驟8-4：求解速度修正方程；

87、步驟8-5：若兩次迭代溫度小于1e-5k，速度小于1e-8ms-1且壓力小于1e-8pa達(dá)到收斂，進(jìn)入步驟9；若不收斂，則返回步驟8-1

88、所述步驟10具體如下：

89、步驟10-1：通過外迭代過程中溫度場最大值和有效增殖系數(shù)keff的變化情況分析確定物理熱工耦合研究的必要性，優(yōu)化液態(tài)金屬冷卻快堆設(shè)計(jì)，提高液態(tài)金屬冷卻快堆安全性；

90、步驟10-2：通過耦合前后的中子通量分布變化和堆芯裂變功率分布的變化，分析液態(tài)金屬冷卻快堆的燃料多普勒溫度反饋和冷卻劑密度反饋現(xiàn)象；

91、步驟10-3：通過耦合前后的堆芯裂變功率分布變化和少群常數(shù)分布變化，分析熱工參數(shù)如何通過影響堆芯及結(jié)構(gòu)材料的截面，從而影響中子行為的三維中子物理效應(yīng)；

92、步驟10-4：通過耦合前后堆芯裂變功率分布變化和冷卻劑溫度場分布變化與流場分布變化，得到中子物理如何影響熱工參數(shù)的分布，并分析物理反饋?zhàn)饔孟乱簯B(tài)金屬冷卻快堆中入口效應(yīng)和流量流經(jīng)堆芯的分配情況，包括橫流、倒流、回流三維熱工水力現(xiàn)象。

93、和現(xiàn)有技術(shù)比較，本發(fā)明具有如下優(yōu)點(diǎn)：

94、1.高保真，計(jì)算精度高。物理和熱工均采用具有能夠描述三維現(xiàn)象特性的物理方程，其中中子物理部分使用三維穩(wěn)態(tài)中子輸運(yùn)方程描述并采用sp3方法求解，熱工水力部分使用三維cfd控制方程描述并采用simple算法將質(zhì)量守恒方程和動量守恒方程變?yōu)樗俣阮A(yù)測方程、壓力泊松方程和速度修正方程，與能量守恒方程進(jìn)行迭代求解，能夠高保真地解析三維物理現(xiàn)象。相比傳統(tǒng)中子擴(kuò)散方法能夠解析空間散射，相比傳統(tǒng)子通道方法能夠解析回流倒流等三維物理現(xiàn)象，相比多孔介質(zhì)方法考慮了燃料棒包殼和冷卻劑的相對位置，更接近實(shí)際物理模型，提高計(jì)算精度。

95、2.計(jì)算速度快。本發(fā)明方法中子物理部分采用sp3方法求解三維穩(wěn)態(tài)中子輸運(yùn)方程，在openfoam中通過線性插值算法，根據(jù)不同的區(qū)域、溫度、密度、能群調(diào)用，獲得熱工反饋?zhàn)饔孟碌纳偃撼?shù)，相比蒙特卡洛方法具有更快的求解速度。熱工水力部分采用子通道簡化網(wǎng)格，相比精細(xì)建模大幅度減少網(wǎng)格數(shù)量，提高計(jì)算速度。物理、熱工部分分別進(jìn)行內(nèi)迭代，提高picard外迭代收斂速度。此外，本發(fā)明在物理熱工耦合計(jì)算時基于統(tǒng)一平臺openfoam進(jìn)行數(shù)據(jù)傳遞，共享內(nèi)存，避免讀寫接口程序傳輸速度慢的問題。

96、3.基于開源平臺，代碼靈活性高，適用范圍廣。本發(fā)明方法完全基于開源平臺開發(fā)，包括openmc少群常數(shù)生成部分，和openfoam堆芯pin-by-pin級物理-熱工耦合求解部分，其中，物理熱工耦合部分采用模塊化方式開發(fā)，代碼靈活性高，便于后續(xù)對方法進(jìn)行改進(jìn)。此外，本發(fā)明方法中，物理熱工均采用相同網(wǎng)格，避免跨尺度計(jì)算帶來的復(fù)雜性和誤差，邊界條件無需特殊化處理，適用范圍廣。