袁寶新，鄭杰，楊萬(wàn)奎，曾和榮

?

基于ICEM-CFD的中子學(xué)力學(xué)耦合程序研制

袁寶新，鄭杰，楊萬(wàn)奎，曾和榮

（中國(guó)工程物理研究院 核物理與化學(xué)研究所，四川 綿陽(yáng) 621900）

針對(duì)反應(yīng)堆工程分析中面臨的中子學(xué)力學(xué)耦合這一問(wèn)題，在同一軟件框架下實(shí)現(xiàn)中子學(xué)力學(xué)多物理場(chǎng)計(jì)算。利用商業(yè)CAD前處理軟件ICEM-CFD對(duì)堆芯幾何進(jìn)行網(wǎng)格剖分，開(kāi)發(fā)基于有限元方法的中子學(xué)力學(xué)耦合計(jì)算程序，對(duì)同一堆芯網(wǎng)格劃分進(jìn)行中子學(xué)和力學(xué)計(jì)算。初步實(shí)現(xiàn)了基于同一軟件框架下網(wǎng)格層面的中子學(xué)力學(xué)耦合計(jì)算，得到了中子學(xué)力學(xué)初步耦合的計(jì)算結(jié)果。通過(guò)基準(zhǔn)例題對(duì)程序進(jìn)行計(jì)算校驗(yàn)，該耦合程序計(jì)算能力是可信的。

ICEM-CFD；中子學(xué)；力學(xué)；耦合；程序

有限元方法因具有非結(jié)構(gòu)網(wǎng)格剖分能力可以適應(yīng)任意幾何形狀，因而從20世紀(jì)60年代起逐漸在中子學(xué)計(jì)算領(lǐng)域得到應(yīng)用和推廣。國(guó)內(nèi)已開(kāi)發(fā)了如3DFEMJS、TFEM-2D、FEM2D-655等程序[1-4]，國(guó)外則有FEM3DJAR、FENDER、FELICIT、FELTRAN等程序[5-6]。這些程序大多采用手工方式進(jìn)行幾何剖分，這就使得這一時(shí)代有限元程序的幾何處理模塊不可避免地具有網(wǎng)格處理繁復(fù)及無(wú)通用性等劣勢(shì)。

為開(kāi)展堆內(nèi)部件振動(dòng)條件下堆芯核噪聲頻譜特性研究，需要開(kāi)發(fā)一套復(fù)雜幾何下的中子學(xué)力學(xué)耦合計(jì)算程序。在這個(gè)背景下，利用商業(yè)CAD前處置軟件ICEM-CFD，對(duì)三維堆芯幾何做四面體或六面體網(wǎng)格剖分，對(duì)二維堆芯幾何做三角形或四邊形網(wǎng)格剖分，開(kāi)發(fā)了基于ICEM-CFD前處理器的中子學(xué)力學(xué)耦合計(jì)算程序FEMN_Couple，可以方便地獲得復(fù)雜堆芯幾何內(nèi)細(xì)致的中子噪聲頻譜。

1 基本方程

1.1 中子學(xué)有限元方程

對(duì)三維中子瞬態(tài)方程做一階微擾和傅里葉變換，并進(jìn)行有限元處理[7-13]。

離散快群實(shí)部方程：

離散快群虛部方程：

離散熱群實(shí)部方程：

離散熱群虛部方程：

1.2 振動(dòng)力學(xué)有限元方程

構(gòu)建結(jié)構(gòu)動(dòng)力學(xué)的有限元方程[14]：

對(duì)于一個(gè)懸臂梁，不計(jì)其軸向位移，其模態(tài)有限元方程為：

特征方程為：

2 ICEM-CFD接口技術(shù)

目前ICEM-CFD的網(wǎng)格輸出支持Abaqus/ANSYS/ AUTODYN/LS-Dyna等格式，通過(guò)對(duì)比這些輸出網(wǎng)格文件的格式，選用LS-Dyna網(wǎng)格文件格式開(kāi)發(fā)FEMN_Couple程序幾何模塊。

1）三維計(jì)算中，LS-Dyna網(wǎng)格文件將網(wǎng)格分為固體網(wǎng)格（相當(dāng)于體網(wǎng)格）和殼網(wǎng)格（相當(dāng)于面網(wǎng)格）。中子學(xué)力學(xué)耦合三維計(jì)算的體網(wǎng)格需求為三棱柱和四棱柱，面網(wǎng)格需求為三角形和四邊形。二維計(jì)算中，LS-Dyna網(wǎng)格文件將網(wǎng)格分為殼網(wǎng)格（相當(dāng)于面網(wǎng)格）和梁網(wǎng)格（相當(dāng)于線網(wǎng)格），中子學(xué)力學(xué)耦合二維計(jì)算的面網(wǎng)格需求為三角形和四邊形，線網(wǎng)格需求為直線。LS-Dyna網(wǎng)格文件的網(wǎng)格類(lèi)型可支持方形燃料組件、六角形燃料組件、棒形燃料組件等類(lèi)型的三維/二維中子學(xué)力學(xué)耦合計(jì)算。

2）中子學(xué)計(jì)算中面中子流為連續(xù)的，在三維計(jì)算中，LS-Dyna網(wǎng)格文件可以不輸出幾何體的內(nèi)部面網(wǎng)格。在二維計(jì)算中，LS-Dyna網(wǎng)格文件可以不輸出幾何體的內(nèi)部線網(wǎng)格。該幾何處理策略有利于系數(shù)矩陣的組裝。

3）LS-Dyna按part來(lái)對(duì)網(wǎng)格賦予材料屬性，每個(gè)組件定義為一個(gè)part，賦予該組件一個(gè)與part編號(hào)對(duì)應(yīng)的材料號(hào)。在幾何前處理時(shí)，按材料號(hào)準(zhǔn)備多群截面文件，程序在單元分析計(jì)算時(shí)，查詢單元所屬的part編號(hào)，根據(jù)該part編號(hào)查詢對(duì)應(yīng)的材料號(hào)，根據(jù)材料號(hào)查詢多群截面文件，結(jié)合單元當(dāng)前的節(jié)點(diǎn)參數(shù)完成矩陣系數(shù)的組裝。

3 程序結(jié)構(gòu)

FEMN_Couple程序的結(jié)構(gòu)如圖1所示，在前處理軟件ICEM-CFD中，對(duì)堆芯建?；?qū)隒AD堆芯幾何模型，手工處理邊界條件，并劃分網(wǎng)格，輸出網(wǎng)格文件。程序讀入ICEM-CFD的網(wǎng)格輸出文件，建立單元、節(jié)點(diǎn)的領(lǐng)域表和稀疏矩陣結(jié)構(gòu)表，手工輸入或從材料文件讀取各區(qū)域的材料信息，讀入多群截面和力學(xué)參數(shù)，填充并裝配系數(shù)矩陣。中子學(xué)穩(wěn)態(tài)參數(shù)分析部分對(duì)堆芯群中子通量和群中子伴隨通量進(jìn)行計(jì)算，為中子學(xué)力學(xué)耦合瞬態(tài)計(jì)算提供穩(wěn)態(tài)中子學(xué)輸入?yún)?shù)。模態(tài)分析部分對(duì)指定結(jié)構(gòu)件進(jìn)行模態(tài)計(jì)算，為中子學(xué)力學(xué)耦合瞬態(tài)計(jì)算提供力學(xué)輸入?yún)?shù)。中子學(xué)力學(xué)耦合瞬態(tài)分析部分，讀入中子學(xué)穩(wěn)態(tài)參數(shù)分析和模態(tài)分析的計(jì)算結(jié)果，計(jì)算中子學(xué)力學(xué)耦合方程的源項(xiàng)，處理邊界條件，求解中子學(xué)力學(xué)耦合方程，得到群中子噪聲頻譜，遍歷各群直至完成，輸出計(jì)算結(jié)果并結(jié)束。

圖1 程序結(jié)構(gòu)

4 程序分步計(jì)算測(cè)試

4.1 中子學(xué)計(jì)算基準(zhǔn)例題測(cè)試

選取文獻(xiàn)[15]發(fā)布的BIBLIS例題進(jìn)行中子學(xué)計(jì)算校驗(yàn)，該堆芯布局如圖2所示，各區(qū)域材料可參考文獻(xiàn)[15]。根據(jù)文獻(xiàn)[16]，在坐標(biāo)(195,195)、(250,195)和(310,195)處分別設(shè)置一個(gè)熱群中子噪聲點(diǎn)源，使用FEMN_Couple程序的中子學(xué)模塊對(duì)該問(wèn)題進(jìn)行計(jì)算，得到的熱群中子噪聲分布如圖3所示?？梢钥闯?，隨著點(diǎn)源位置從堆芯正中心往外圈移動(dòng)，熱群中子噪聲的波峰也隨之向外圈移動(dòng)，熱群中子噪聲波形與文獻(xiàn)[16]計(jì)算結(jié)果基本一致。

圖2 BIBLIS堆芯布局

4.2 力學(xué)計(jì)算例題測(cè)試

對(duì)于固定組件，其在堆芯內(nèi)部的振動(dòng)可以等效為一個(gè)懸臂梁振動(dòng)問(wèn)題，不計(jì)其軸向位移有：

圖3 源位置在不同坐標(biāo)引起的熱群中子噪聲幅值分布

5 中子學(xué)力學(xué)耦合計(jì)算

振動(dòng)情況下的吸收截面可以表示為：

部件振動(dòng)條件下的噪聲表示為：

選取文獻(xiàn)[15]發(fā)布的BIBLIS例題，將該中心組件的一階振動(dòng)頻率代入上述方程，使用FEMN_ Couple程序進(jìn)行中子學(xué)力學(xué)耦合計(jì)算，得到熱群中子噪聲實(shí)部波形如圖4所示。該結(jié)果與文獻(xiàn)[16]給出的控制棒振動(dòng)情況下堆芯熱群中子的實(shí)部波形分布基本一致。

圖4 中心組件的一階振動(dòng)頻率下的熱群中子噪聲實(shí)部

6 結(jié)語(yǔ)

針對(duì)反應(yīng)堆工程分析中面臨的中子學(xué)力學(xué)耦合這一問(wèn)題，文中利用商業(yè)CAD前處理軟件ICEM-CFD對(duì)堆芯幾何進(jìn)行了網(wǎng)格剖分，開(kāi)發(fā)了基于有限元方法的中子學(xué)力學(xué)耦合計(jì)算程序FEMN_Couple，對(duì)同一堆芯網(wǎng)格劃分進(jìn)行中子學(xué)和力學(xué)計(jì)算，初步實(shí)現(xiàn)了基于同一軟件框架下網(wǎng)格層面的中子學(xué)力學(xué)耦合，給出了中子學(xué)力學(xué)耦合的初步計(jì)算結(jié)果。下一步將在程序支持部件受迫振動(dòng)力學(xué)計(jì)算及部件受迫振動(dòng)中子學(xué)力學(xué)耦合計(jì)算兩個(gè)方面進(jìn)一步開(kāi)展工作。

[1] 屠柱國(guó), 湯裕仁, 黃艾香, 等. 用有限元法求解三維中子擴(kuò)散方程的研究—三維中子擴(kuò)散有限元程序3DFEMJS[J]. 核科學(xué)與工程, 1983, 3(2): 115-126.

[2] 顧麗珍, 王永慶, 于素花, 等. 有限元法在計(jì)算二維反應(yīng)堆中子擴(kuò)散方程中的應(yīng)用[J]. 清華大學(xué)學(xué)報(bào), 1980, 20(4): 97-110.

[3] 郭海兵, 劉永康, 施工. 基于CAD建模和有限元方法的三維中子擴(kuò)散程序[J]. 清華大學(xué)學(xué)報(bào), 2012, 52(7): 901-905.

[4] 薛友義. 有限元法在解中子擴(kuò)散方程中的應(yīng)用[J]. 核動(dòng)力工程, 1981, 2(2): 10-20.

[5] YOSHITAKA N RYUJI K, YUKIO T, et al. A Computer Program for Solving the Three-dimensional Multigroup Diffusion Equation by Finite Element Method with 20-node Isoparametric Element (FEM3DJAR)[J]. Progress in Nuclear Energy, 1986, 18(1-2): 207-214.

[6] SHUTTLEWORTH E. Fender: A Finite Element Code for the Solution of the Diffusion Equation in Shield Design Applications[J]. Annals of Nuclear Energy, 1981, 8(11-12): 597-607.

[7] 劉金匯, 谷芳毓. 反應(yīng)堆堆內(nèi)部件振動(dòng)的中子噪聲物理模型[J]. 核動(dòng)力工程, 1999, 20(1): 36-41.

[8] 彭鋼. 反應(yīng)堆堆內(nèi)部件振動(dòng)的中子噪聲物理模型[J]. 核科學(xué)與工程, 2001, 21(3): 264-270.

[9] 彭鋼. 反應(yīng)堆冷卻劑沸騰中子噪聲物理模型研究[J]. 核動(dòng)力工程, 2000, 21(5): 385-388.

[10] CHRISTOPHE D, IMRE P. Numerical Tools Applied to Power Reactor Noise analysis[J]. Progress in Nuclear Energy, 2009, 51: 67-81.

[11] VIKTOR L, CHRISTOPHE D, IMRE P et al. Neutron Noise Calculations Using the Analytical Nodal Method and Comparisons with Analytical Solutions[J]. Annals of Nuclear Energy, 2011, 38:808-816.

[12] CHRISTOPHE D. Development of a 2-D 2-group Neutron Noise Simulator[J]. Annals of Nuclear Energy, 2004, 31: 647-680.

[13] CHRISTOPHE D. Comparison of the Calculated Neutron Noise Using Finite Differences and the Analytical Nodal Method[J]. Annals of Nuclear Energy, 2012, 43: 176-182.

[14] 張相庭, 王志培, 黃本才, 等. 結(jié)構(gòu)振動(dòng)力學(xué)[M]. 上海: 同濟(jì)大學(xué)出版社, 2005.

[15] SEYED A H, FARAHNAZ S D. Galerkin and Generalized Least Squares Finite Element: A Comparativestudy for Multi-group Diffusion Solvers[J]. Progress in Nuclear Energy, 2015, 85: 473-490.

[16] SEYED A H, NASER V. Neutron Noise Simulation by GFEM and Unstructured Triangle Elements[J]. Nuclear Engineering and Design, 2012, 253: 238-258.

Research on Coupling Code of ICEM-CFD Based Neutronic and Mechanics

YUAN Bao-xin, ZHENG Jie, YANG Wan-kui, ZENG He-rong

(Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China)

In view of the neutronic and mechanical coupling problems in reactor engineering analysis, to realize the multi physics field calculation under the same software framework.The CAD pre-processing software ICEM-CFD was used to carry out mesh subdivision of core geometry to develop neutron computation module and the mechanical calculation module based on the finite element method, and execute neutronic and mechanical calculations for the same mesh subdivision of core geometry.It preliminary realized mechanical coupled neutronics calculation of mesh based on the software framework and obtained the calculation results of preliminary neutronics mechanics coupling.Calculation and verification of code based on benchmark examples shows that the code for calculating the capacity can be trusted.

ICEM-CFD; neutronic; mechanics; coupling; code

10.7643/ issn.1672-9242.2019.02.007

TL329

A

1672-9242(2019)02-0032-05

2018-11-22；

2018-12-18

國(guó)家自然科學(xué)基金面上項(xiàng)目（11475150）；四川省科技計(jì)劃項(xiàng)目（2018JY0622）。

袁寶新（1988—），男，湖北黃岡人，碩士研究生，主要研究方向?yàn)榉磻?yīng)堆物理。