彌散核燃料等效輻照腫脹計算模擬

蔡 維，趙云妹，龔 辛，丁淑蓉，霍永忠

（復旦大學力學與工程科學系，上海 200433）

彌散核燃料等效輻照腫脹計算模擬

蔡 維，趙云妹，龔 辛，丁淑蓉＊，霍永忠

（復旦大學力學與工程科學系，上海 200433）

本文將彌散核燃料芯體看作一種特殊的顆粒復合材料，利用細觀計算力學的方法，假設燃料顆粒在芯體中周期性分布，建立了對芯體等效輻照腫脹進行計算模擬的有限元模型?？紤]顆粒的輻照腫脹和基體材料的輻照硬化效應，分別建立了燃料顆粒和基體材料的應力更新算法，編制了用戶材料子程序，在Abaqus軟件中實現(xiàn)了芯體等效輻照腫脹的有限元模擬。計算分析了顆粒大小和體積含量對芯體等效輻照腫脹的影響，并得到了等效輻照腫脹的擬合公式。研究結(jié)果表明，影響芯體等效輻照腫脹的主要因素是顆粒的輻照腫脹和體積含量。

彌散核燃料；有限元法；均勻化方法；等效輻照腫脹

相對于傳統(tǒng)的棒狀核燃料元件，彌散核燃料元件［1］具有高熱導，可達較高的燃耗，已廣泛應用于研究試驗堆中，在先進核反應堆和核廢料處理方面也具有良好的應用前景［1－2］。

彌散核燃料元件由彌散核燃料芯體和包殼組成，分棒狀和板型兩種形式，是組成彌散核燃料組件的基本構(gòu)件。彌散核燃料芯體在結(jié)構(gòu)上類似于顆粒復合材料，由一定體積含量的燃料顆粒彌散分布于金屬基體中構(gòu)成。為達到一定的功率密度，芯體中需裝載大量的燃料顆粒。

為對彌散核燃料組件進行優(yōu)化設計，需研究其在堆內(nèi)的輻照－熱力耦合行為。輻照試驗研究耗時、花費極高且無法在線觀測，建立元件或組件堆內(nèi)輻照熱力耦合行為的理論模型和數(shù)值模擬手段，對其堆內(nèi)行為的演化進行預測成為一種重要的輔助手段［3－5］。對彌散核燃料元件或組件的堆內(nèi)熱力耦合行為進行數(shù)值模擬，需將彌散核燃料芯體均勻化［6］，等效成一種均勻的介質(zhì)，不再考慮顆粒和基體的相互作用。

文獻［7－8］運用均勻化理論，通過計算模擬對彌散核燃料的等效熱傳導率和彈性性能進行了研究，而對其等效輻照腫脹的計算模擬研究尚未見報道。在核反應堆苛刻的環(huán)境中，核燃料顆粒所產(chǎn)生的裂變產(chǎn)物導致其發(fā)生輻照腫脹，并隨燃耗增長而增大［9］，成為影響燃料元件或組件堆內(nèi)安全的一重要的載荷。

本文假設燃料顆粒在芯體中周期性分布，根據(jù)燃料芯體的微觀結(jié)構(gòu)，選取代表性體元（RVE）和有限元模型，考慮金屬基體的輻照硬化效應，建立計算模擬燃料芯體等效輻照腫脹的數(shù)值方法。同時，計算分析顆粒體積含量和大小對芯體等效輻照腫脹的影響，并對解析公式進行擬合。

1 有限元模型

1．1 有限元幾何模型

考慮無限大的彌散核燃料，如圖1a所示，假設燃料顆粒在基體中呈簡單立方分布。此彌散核燃料芯體可看作圖1b所示的代表性體元在空間周期性排列而成。由于代表性體元包含了燃料芯體所有的微觀結(jié)構(gòu)信息，可用代表性體元代替整個燃料芯體，對其等效輻照腫脹性能進行研究。在代表性體元中，燃料顆粒隨燃耗的發(fā)展產(chǎn)生輻照腫脹，這導致體元隨燃耗的發(fā)展而產(chǎn)生相對的體積變化，稱為等效輻照腫脹。體元位于無限大復合燃料芯體中，根據(jù)對稱性，隨燃耗的發(fā)展體元邊界保持為平面?？紤]到體元結(jié)構(gòu)和載荷的對稱性，選取其1／8作為有限元模型進行研究，如圖1c所示。將其進行空間離散，如圖1d所示。

為計算分析燃料顆粒體積含量變化對芯體等效腫脹的影響，本文分別選取顆粒體積含量為10．657 4%、15．514 0%、19．081 6%、23．832 4%、30．300 9%的幾種情況進行建模。假設顆粒的半徑r為0．05mm，有限元幾何模型的邊長列于表1。對于顆粒尺寸不同的燃料芯體，可類似地建立有限元模型。

1．2 邊界條件

根據(jù)前面的分析可知，彌散核燃料芯體的等效輻照腫脹需通過計算有限元模型的體積變化得到。有限元模型中的主要載荷來自于顆粒的輻照腫脹，其使得顆粒和基體之間產(chǎn)生強烈的力學相互作用，導致顆粒和基體內(nèi)部產(chǎn)生大應變、有限元模型的棱邊產(chǎn)生變化。針對圖1c所施加的位移邊界條件如下：

圖1 燃料顆粒在芯體中的分布形式（a）、代表性體元（b）、有限元幾何模型（c）及有限元網(wǎng)格（d）Fig．1 Distribution pattern of fuel particles in fuel meat（a），RVE（b），finite element geometric model（c）and finite element grid（d）

表1 顆粒體積含量與有限元幾何模型邊長的關(guān)系Table 1 Relationship between particle volume fraction and size of finite element geometric model

1）分別對x＝L、y＝L和z＝L的這3個面，施加節(jié)點位移約束方程，使這3個面始終保持為平面；

2）對x＝0、y＝0和z＝0的這3個面，根據(jù)結(jié)構(gòu)和載荷對稱性，施加對稱邊界條件。

1．3 材料參數(shù)

1）燃料顆粒的材料參數(shù)［10－12］

（1）彈性模量和泊松比

式中：E為彈性模量，Pa；T為溫度，K；D為顆粒的相對理論密度，取95%。

泊松比ν取0．316。

（2）塑性模型

（3）腫脹經(jīng)驗公式

燃料顆粒的體積腫脹SWfuel包括裂變氣體腫脹（gs）、裂變固體腫脹（ss）和裂變收縮（ds）3部分。

式中，V0為燃料顆粒的原始體積。裂變氣體腫脹、裂變固體腫脹導致的每產(chǎn)生1020cm－3的裂變密度下的體積變化分數(shù)［13－14］分別為：

裂變密實所導致的體積收縮［15］為：

式中：BU為燃耗，MW·d／kgUO2；ADST為調(diào)整因子，取0．6。

2）金屬基體的材料參數(shù)

受快中子的沖擊，基體材料將產(chǎn)生輻照硬化效應，下面給出考慮輻照硬化效應的基體材料性能，與文獻［11］中所用模型相同。

（1）彈性模量和泊松比［16］

式中：φ為快中子注量率，m－2·s－1；t為時間，s。

（2）輻照硬化塑性模型［17］

式中：σ為真應力，Pa；ε為真應變；˙ε為真應變率，若˙ε＜10－5s－1，設˙ε＝10－5s－1；參數(shù)K、n和m分別為強度系數(shù)、應變硬化指數(shù)和應變率敏感指數(shù)。本文中快中子注量率φ取2× 1018m－2·s－1，T取300～730K。

2 數(shù)值模擬方法

為對彌散核燃料芯體的等效輻照腫脹進行計算，需將所考察的燃耗區(qū)間分為多個時間步。針對前面給出的特殊的材料性能，考慮大變形效應，需在旋轉(zhuǎn)坐標系下［12］分別建立燃料顆粒和基體材料的三維應力更新算法。

2．1 燃料顆粒的應力更新算法

對于一典型的時間步［t，t＋Δt］，t時刻的柯西應力σtij和彈性對數(shù)應變εiej，t滿足如下關(guān)系：

式中，λ和G為拉梅常數(shù)，由材料的彈性模量和泊松比計算得到。則t＋Δt時刻的柯西應力σt＋Δtij為：

式中：Δσij為柯西應力增量；Δεeij為彈性對數(shù)應變的增量。對于燃料顆粒內(nèi)部的每個積分點，有：

其中：Δεtotalij為總的應變增量；Δεswij和Δεpij分別為輻照腫脹和塑性應變增量。

先假設無塑性應變產(chǎn)生，即Δεpij＝0，計算得到試驗應力σprij，并由此得到Mises等效應力，其中，sprij為試驗應力的偏斜部分，sprij＝σprij－σprkkδij／3。若為根據(jù)式（2）計算所得的t時刻的等效塑性應變），則說明此增量步有新的塑性應變產(chǎn)生，否則，試驗應力σprij為t＋Δt時刻的真實應力σt＋Δtij，應力更新完成。

對于有新塑性應變產(chǎn)生的情況，由塑性流動法則［12］：

對式（17）進行后歐拉積分，得到：

則

2．2 金屬基體的應力更新算法

考慮輻照硬化效應，基體材料的彈性參數(shù)和塑性模型隨時間而變化。對于時間步［t，t＋Δt］，積分點的彈性應變增量為：

t＋Δt時刻的柯西應力σt＋Δtij為：

式中：G（t＋Δt）和λ（t＋Δt）分別為積分點t＋Δt時刻的拉梅常數(shù)；ΔG和Δλ為因時間變化所導致的拉梅常數(shù)的變化。

針對2．1和2．2節(jié)的應力更新算法，分別編制了定義材料力學本構(gòu)關(guān)系的子程序，并引入ABAQUS軟件的有限元計算中，實現(xiàn)彌散核燃料等效輻照腫脹的計算模擬。

3 數(shù)值模擬方法有效性的驗證

針對體積含量為10．657 4%、溫度為600K的彌散核燃料，計算其隨燃耗增長的微觀的位移、應力和應變場，并根據(jù)有限元計算的結(jié)果來驗證所建立的數(shù)值模擬方法的有效性。

3．1 輻照硬化塑性的驗證

3．2 顆粒輻照腫脹的驗證

圖2 基體中的節(jié)點A（a）及節(jié)點A計算結(jié)果的驗證（b）Fig．2 Node A in matrix（a）and validation of computation results of node A（b）

圖3 燃料顆粒體積變化的有限元計算結(jié)果與輻照腫脹理論值的比較Fig．3 Comparison of volumetric variation result obtained by FEM with theoretical irradiation swelling result

一方面，通過有限元計算可得到燃料顆粒邊緣處節(jié)點在不同燃耗下的位移ΔR，進而可計算燃料顆粒在不同燃耗下的實際體積變化ΔV／V0＝（（R0＋ΔR）3－R30）／R30，其中，R0為燃料顆粒的原始半徑。另一方面，根據(jù)式（3）可計算出燃料顆粒不同燃耗下的理論輻照腫脹結(jié)果（需對不同的燃耗單位進行換算，參見文獻［10］）。通過上述兩種途徑得到的結(jié)果示于圖3。由圖3a可明顯看出，在很低的燃耗下，燃料顆粒的體積變化主要是輻照所致裂變收縮引起的，實際體積變化為負，在FIMA為0．3%左右出現(xiàn)了最低點。且在很低的燃耗下，有限元計算所得燃料顆粒的實際體積縮小量（絕對值）低于輻照腫脹的理論計算結(jié)果（絕對值）。由于燃料顆粒的實際體積變化包括輻照腫脹和彈性變化兩部分。彈性變形所導致的體積變化主要與燃料顆粒內(nèi)各點的靜水壓力－σkk／3有關(guān)（圖4）。

由圖4a所示的FIMA為0．3%的靜水壓力云圖可看出，燃料顆粒內(nèi)各點的靜水壓力均為負，即靜水拉力有使燃料顆粒體積增大的貢獻，抵消了一部分輻照所致裂變收縮的貢獻，所以，燃料顆粒在很低燃耗下的實際體積縮小量低于理論的輻照腫脹結(jié)果。與圖3a不同，由圖3b可看出，在較高燃耗下，燃料顆粒的實際體積變化略低于理論輻照腫脹結(jié)果。同樣，從圖4b所示的FIMA為30%的靜水壓力云圖看，F(xiàn)IMA為30%的燃料顆粒各點具有很大的靜水壓力，其驅(qū)使燃料顆粒的體積縮小，抵消了一部分燃料顆粒的輻照腫脹，所以，有限元計算所得實際體積變化低于輻照腫脹理論結(jié)果是合理的。通過上述分析，可看出燃料顆粒的輻照腫脹得到了正確的模擬。

圖4 不同燃耗下的靜水壓力云圖Fig．4 Contour plots of hydrostatic pressure with different burnups

4 結(jié)果與討論

本節(jié)針對溫度為600K的彌散核燃料，計算分析顆粒體積含量和大小對其等效輻照腫脹的影響，并得到等效輻照腫脹的擬合公式。

4．1 網(wǎng)格收斂性考察

為保證有限元計算結(jié)果的精度和計算的經(jīng)濟性，對不同顆粒體積含量情況下的有限元模型，由粗到細逐步加密網(wǎng)格，并考察結(jié)果的收斂性，力求得到一種經(jīng)濟高效的有限元網(wǎng)格。網(wǎng)格的逐步加密圖如圖5所示。針對顆粒體積含量為20%的情況，給出了FIMA為30%下有限元模型Path 1上的位移和Mises應力隨網(wǎng)格加密而變化的情況，如圖6所示。從圖6可看出，位移很易收斂；而對于Mises應力，隨著網(wǎng)格的進一步加密才得到收斂?？紤]到經(jīng)濟性，可采用具有3 024個單元的網(wǎng)格進行計算。對于其他體積含量的有限元模型，也采用類似的方法獲得經(jīng)濟性和計算精度兼具的網(wǎng)格。

圖5 網(wǎng)格逐步加密圖Fig．5 Refinement of finite element mesh

圖6 不同網(wǎng)格下Path 1上各點的位移和Mises應力Fig．6 Displacements and Mises stresses on Path 1for different cases of mesh grid

4．2 顆粒體積含量對芯體等效輻照腫脹的影響

對不同體積含量的有限元模型進行計算，得到的芯體的等效輻照腫脹如圖7所示。從圖7可看到，芯體的等效輻照腫脹隨顆粒體積含量的增加而增加。在較低燃耗下，高體積含量的情況體積收縮也較大；對于不同體積含量的情況，由收縮轉(zhuǎn)為膨脹的轉(zhuǎn)變?nèi)己膸缀跸嗤?，與燃料顆粒的轉(zhuǎn)變?nèi)己南嘟晦D(zhuǎn)變?nèi)己暮?，等效輻照腫脹與燃耗呈線性關(guān)系；對于顆粒體積含量大的芯體，其等效輻照腫脹的變化率也較高。由于受到金屬基體的約束，在FIMA為30%以下時，等效輻照腫脹均小于6%。

4．3 顆粒大小對芯體等效輻照腫脹的影響

對于體積含量分別為10．657 4%和30．300 9%兩種情況，分別針對0．02、0．05、0．08mm的顆粒建立有限元模型，對芯體的等效腫脹進行計算。計算的不同顆粒大小下的等效輻照腫脹如圖8所示?？煽闯鲱w粒大小對芯體等效腫脹幾乎無影響。在一定燃耗下，在不同顆粒大小的有限元模型中，顆粒和基體之間的力學相互作用存在差異，而這種差異對芯體等效輻照腫脹的總體貢獻幾乎為零，故造成圖4所示的結(jié)果。

圖7 不同體積含量下的等效輻照腫脹Fig．7 Equivalent irradiation swelling for different particle volume fractions

圖8 顆粒大小對等效輻照腫脹的影響Fig．8 Effect of particle size on equivalent irradiation swelling

4．4 芯體等效輻照腫脹的擬合公式

通過以上計算分析可得出，芯體的等效輻照腫脹主要決定于顆粒的輻照腫脹和體積含量以及燃料顆粒和基體中彈性變形的貢獻，而彈性變形的貢獻相對較小。對于所建立的有限元模型，將其體積的絕對變化首先近似為：

其中：ΔVfuel為核燃料顆粒部分的絕對體積變化；Vfuel0為核燃料顆粒部分的原始體積。

據(jù)此，可將芯體的等效輻照腫脹初步近似為：

其中，Vf為顆粒的原始體積含量。SW為有限元計算結(jié)果，Vf×SWfuel為初步近似值。根據(jù)有限元計算所得的位移結(jié)果可得到不同燃耗下的實際棱邊長度a＋Δa（a為原始的棱邊長度，Δa為棱邊的絕對變形量），進而可得到SW＝（（a＋Δa）3－a3）／a3。同時，由式（3）所得燃料顆粒的SWfuel可計算得到不同燃耗下的Vf×SWfuel。

將芯體等效輻照腫脹的有限元計算值與通過式（25）得到的初步近似結(jié)果進行比較，結(jié)果如圖9所示。從圖9可看出，二者在較低燃耗下仍存在較大誤差，還需進一步擬合。

為達到更精確的近似值，進行分段考慮。

當FIMA≤3%時，有：

當FIMA＞3%時，有：

有限元計算結(jié)果與利用式（26）、（27）得到的擬合結(jié)果的比較示于圖10?？煽闯觯檬剑?6）、（27）可得到更好的擬合結(jié)果，可將式（26）、（27）作為彌散核燃料等效輻照腫脹的擬合公式。

圖9 芯體等效腫脹有限元計算結(jié)果與近似值的比較Fig．9 Comparison of finite element calculation result with approximate result

圖10 芯體等效腫脹有限元計算結(jié)果與分段考慮近似值的比較Fig．10 Comparison of finite element calculation result with segment approximate result

5 結(jié)論

本研究建立了對彌散核燃料等效輻照腫脹進行計算模擬的有限元模型，并考慮基體材料的輻照硬化效應，建立了對等效輻照腫脹進行計算模擬的數(shù)值方法，分別考察了燃料顆粒大小、體積含量對芯體等效輻照腫脹的影響，得到了等效輻照腫脹的擬合公式。在所考察的模型和參數(shù)范圍內(nèi)，得到的主要結(jié)論如下：

1）芯體等效輻照腫脹主要受顆粒輻照腫脹和體積含量的影響。在較低燃耗下，燃料顆粒的裂變收縮起主要作用，等效輻照腫脹為負值；在較高燃耗下，等效輻照腫脹主要來自于顆粒固體裂變產(chǎn)物的貢獻，與燃耗呈線性關(guān)系。

2）本文所建立的有限元模型和數(shù)值模擬方法是有效的。

［1］ XU Z．Design strategies for optimizing high burnup fuel in pressurized water reactors［D］．US：Massachu－setts Institute of Technology，2003．

［2］ DUYN L．Evaluation of the mechanical behavior of a metal－matrix dispersion fuel for plutonium burning［D］．US：Georgia Institute of Technology，2003．

［3］ WANG Q，YAN X，DING S，et al．Research on the interfacial behaviors of plate－type dispersion nuclear fuel elements［J］．Journal of Nuclear Materials，2010，399（1）：41－54．

［4］ WANG Q，YAN X，DING S，et al．Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effectⅠ：Effects of variations of the fuel particle volume fractions［J］．Journal of Nuclear Materials，2010，400（2）：157－174．

［5］ DING S，WANG Q，HUO Y．Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effectⅡ：Effects of variations of the fuel particle diameters［J］．Journal of Nuclear Materials，2010，397（1－3）：80－91．

［6］ 喬生儒．復合材料細觀力學性能［M］．西安：西北工業(yè)大學出版社，1997．

［7］ 姜馨，丁淑蓉，霍永忠．彌散型燃料等效彈性性質(zhì)的有限元模擬［J］．核動力工程，2011，32（3）：43－47．

JIANG Xin，DING Shurong，HUO Yongzhong．Finite element simulation for equivalent elastic properties of dispersion fuel elements［J］．Nuclear Power Engineering，2011，32（3）：43－47（in Chinese）．

［8］ 姜馨，丁淑蓉，霍永忠．彌散型燃料元件等效熱傳導系數(shù)的有限元模擬［J］．核動力工程，2010，31（3）：45－49．

JIANG Xin，DING Shurong，HUO Yongzhong．Finite element simulation for effective thermal conductivity coefficient of dispersion fuel element［J］．Nuclear Power Engineering，2010，31（3）：45－49（in Chinese）．

［9］ MATPRO－09：A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior［M］．US：NRC，1976．

［10］SUZUKI M，SAITOU H．Light water reactor fuel analysis code FEMAⅪ－6［R］．Japan：Atomic Energy Agency，2005．

［11］GONG X，DING S，ZHAO Y，et al．Effects of irradiation hardening and creep on the thermomechanical behaviors in inert matrix fuel elements［J］．Mechanics of Materials，2013，65：110－123．

［12］BELYTSCHKO T，LIU W K，MORAN B．Nonlinear finite elements for continua and structures［M］．UK：John Wiley ＆Sons Ltd．，2000．

［13］CHUBB W，STORHOK V W，KELLER D L．Factors affecting the swelling of nuclear fuel at high temperatures［J］．Nuclear Technology，1973，18：231－255．

［14］NAKAJIMA T，ICHIKAWA M，IWANO Y，et al．FEMAⅪ－Ⅲ：A computer code for the analysis of thermal and mechanical behavior of fuel rods［R］．Japan：JAERI，1985．

［15］WIESENACK W，VANKEERGHEN M，THANKAPPAN R．Assessment of UO2conductivity degradation based on in－pile temperature data［R］．Norway：HWR，1996．

［16］FISHER E F，RENKEN C J．Single－crystal elastic moduli and the hcp－transformation in Ti，Zr and Hf［J］．Physical Review A，1954，135（2）：482－494．

［17］HAGRMAN D L，REYMAN G A．MATPRO－11：A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior［M］．US：NRC，1979．

Calculation Simulation of Equivalent Irradiation Swelling for Dispersion Nuclear Fuel

CAI Wei，ZHAO Yun－mei，GONG Xin，DING Shu－rong＊，HUO Yong－zhong
（Department of Mechanics and Engineering Science，F(xiàn)udan University，Shanghai 200433，China）

The dispersion nuclear fuel was regarded as a kind of special particle composites．Assuming that the fuel particles are periodically distributed in the dispersion nuclear fuel meat，the finite element model to calculate its equivalent irradiation swelling was developed with the method of computational micro－mechanics．Considering irradiation swelling in the fuel particles and the irradiation hardening effect in the metal matrix，the stress update algorithms were established respectively for the fuel particles and metal matrix．The corresponding user subroutines were programmed，and the finite element simulation of equivalent irradiation swelling for the fuel meat was performed in Abaqus．The effects of the particle size and volume fraction on the equivalent irradiation swelling were investigated，and the fitting formula of equivalent irradiation swelling was obtained．The results indicate that the main factors to influence equivalent irradiation swelling of the fuel meat are the irradiation swelling and volume fraction of fuel particles．

dispersion nuclear fuel；finite element method；homogenization method；equivalent irradiation swelling

TL341

：A

：1000－6931（2015）03－0502－09

10．7538／yzk．2015．49．03．0502

2013－12－20；

2014－03－08

國家自然科學基金資助項目（11172068，91226101，11272092）；教育部博士點基金資助項目（20110071110013）

蔡 維（1987—），男，安徽六安人，碩士研究生，從事彌散燃料的等效輻照性能研究

＊通信作者：丁淑蓉，E－mail：dsr1971＠163．com