單片式核燃料板軋制過程的數(shù)值模擬研究

萬繼波，孔祥喆，丁淑蓉，徐宏彬，霍永忠

（復(fù)旦大學(xué)力學(xué)與工程科學(xué)系，上海 200433）

單片式核燃料板軋制過程的數(shù)值模擬研究

萬繼波，孔祥喆，丁淑蓉＊，徐宏彬，霍永忠

（復(fù)旦大學(xué)力學(xué)與工程科學(xué)系，上海 200433）

針對(duì)UMo合金單片式核燃料板鋯合金包殼材料應(yīng)變率相關(guān)的力學(xué)本構(gòu)關(guān)系，推導(dǎo)出其三維應(yīng)力更新算法，相應(yīng)地編寫了定義其本構(gòu)關(guān)系的VUMAT子程序并驗(yàn)證了程序的正確性；建立了對(duì)UMo合金單片式板元件的框架軋制過程進(jìn)行計(jì)算模擬的有限元模型；利用顯式動(dòng)力有限元法，計(jì)算分析了復(fù)合坯內(nèi)部的變形以及接觸壓強(qiáng)在軋制過程中的演化規(guī)律。研究結(jié)果表明，利用VUMAT用戶材料子程序能方便正確地定義材料應(yīng)變率相關(guān)的本構(gòu)關(guān)系；燃料芯體與蓋板之間的軋制接觸壓力隨時(shí)間而演化，在靠近寬度方向的對(duì)稱面處具有最大的接觸壓力。本研究為優(yōu)化UMo合金單片式核燃料板的制造工藝參數(shù)提供了理論基礎(chǔ)和計(jì)算手段。

單片式核燃料板；軋制；數(shù)值模擬；大變形彈塑性；應(yīng)變率相關(guān)

UMo合金燃料［1］具有鈾密度高、輻照穩(wěn)定性好等優(yōu)點(diǎn)，UMo合金單片式核燃料板在高通量研究試驗(yàn)堆中具有良好的發(fā)展前景?？蚣苘堉茝?fù)合法是一種有效制備UMo合金單片式核燃料板的工藝方法［2］。結(jié)合工藝試驗(yàn)，利用有限元法模擬軋制過程，成為優(yōu)化工藝參數(shù)的發(fā)展趨勢(shì)［3］。

UMo合金單片式核燃料板的制造質(zhì)量，如芯體與其他部分的結(jié)合質(zhì)量，與其堆內(nèi)能密切相關(guān)。根據(jù)軋制結(jié)合的理論［4－12］，只有界面接觸壓力或等效應(yīng)變達(dá)到一定程度才能使不同金屬在界面處達(dá)到冶金結(jié)合，而結(jié)合強(qiáng)度受到多個(gè)工藝參數(shù)的影響［13－16］。為了提高燃料板的加工質(zhì)量，需優(yōu)化其加工工藝參數(shù)。對(duì)于軋制過程的有限元模擬，目前大部分的研究均針對(duì)單片板進(jìn)行［17－19］，雖然存在3層板軋制的二維解析分析［4］，但對(duì)于燃料板復(fù)合坯這種具有復(fù)雜結(jié)構(gòu)工件的軋制模擬研究尚未見報(bào)道。

本工作推導(dǎo)鋯合金材料應(yīng)變率相關(guān)的三維應(yīng)力更新算法，編寫定義材料力學(xué)本構(gòu)關(guān)系的VUMAT子程序，驗(yàn)證其正確性；針對(duì)UMo單片燃料板的框架軋制工藝，建立對(duì)單片燃料板的軋制過程進(jìn)行數(shù)值模擬的有限元模型；利用顯式動(dòng)力學(xué)模擬方法，計(jì)算分析燃料板復(fù)合坯內(nèi)部的變形以及接觸壓強(qiáng)在軋制過程中的演化規(guī)律，為板形核燃料元件的先進(jìn)制造提供理論基礎(chǔ)和計(jì)算手段。

1 理論模型

1．1 材料屬性

1）UMo合金材料性能

本文主要對(duì)復(fù)合坯在800℃下的軋制過程進(jìn)行模擬，此溫度下含Mo量為10%（質(zhì)量分?jǐn)?shù)）的材料性能參數(shù)［1，20］列于表1。

表1 UMo合金材料性能參數(shù)Table 1 Material performance parameter of UMo alloy

2）鋯合金材料的材料屬性

鋯合金的彈性模量和泊松比采用如下模型［21］：

式中：E為未受輻照時(shí)鋯合金的彈性模量，Pa；T為溫度，K。

鋯合金未受輻照時(shí)的應(yīng)力應(yīng)變曲線為：

1．2 應(yīng)變率相關(guān)塑性材料的三維應(yīng)力更新算法

如式（3）所示，鋯合金材料具有與應(yīng)變率相關(guān)的硬化曲線，要實(shí)現(xiàn)燃料板復(fù)合坯軋制過程的數(shù)值模擬，需導(dǎo)出其對(duì)應(yīng)的三維應(yīng)力更新算法，編制用戶材料子程序，方便引入軋制的有限元分析中。

首先，假設(shè)該時(shí)間步無新的塑性應(yīng)變產(chǎn)生，得到試驗(yàn)應(yīng)力：

其中，G和λ為材料的拉梅常數(shù)。

求解試驗(yàn)應(yīng)力σprij所對(duì)應(yīng)的Mises應(yīng)力為：

對(duì)于有新塑性應(yīng)變產(chǎn)生的情況，需考慮塑性應(yīng)變?cè)隽繉?duì)應(yīng)力的影響。根據(jù)塑性本構(gòu)理論，對(duì)微小塑性應(yīng)變?cè)隽縟εpij［22］，有：

對(duì)式（8）右邊進(jìn)行后歐拉積分，得到該時(shí)間步的塑性應(yīng)變?cè)隽浚?/p>

則t＋Δt時(shí)刻的柯西應(yīng)力為：

其偏斜部分為：

則：

在小學(xué)數(shù)學(xué)教學(xué)的過程中，習(xí)題練習(xí)能夠有效地幫助學(xué)生鞏固知識(shí)記憶，同時(shí)也能夠幫助教師了解學(xué)生的學(xué)習(xí)難點(diǎn)，從而能夠有針對(duì)性地進(jìn)行講解。而在此過程中，教師也可以結(jié)合微課教學(xué)視頻開展復(fù)習(xí)工作，提升學(xué)生的學(xué)習(xí)效率。

對(duì)式（12）兩端自乘，得：

可得：

即：

將式（17）代入式（16），得到Δˉεp所滿足的非線性方程：

收斂條件為：

其中：

按照上述算法，編制了用戶材料子程序VUMAT，并引入軋制過程的有限元計(jì)算中。

2 有限元模型

2．1 幾何模型

復(fù)合坯件（圖1a）由鋯合金上下蓋板、框架和UMo單片式燃料芯體裝配成整體，并在四周焊接成多材質(zhì)復(fù)合坯件［2］，經(jīng)過多道次的熱軋和冷軋工藝，最后實(shí)現(xiàn)芯體與包殼之間的冶金結(jié)合，達(dá)到預(yù)定的形狀和尺寸。軋制過程中軋輥與復(fù)合坯件之間產(chǎn)生相互作用，其具體參數(shù)列于表2［3］。由于對(duì)稱性，可取整個(gè)工件的1／4和一個(gè)軋輥的1／2建立有限元模型（圖1b）。軋輥與可產(chǎn)生較大彈塑性變形的軋板相比，其變形量極小，在此用解析剛體殼來描述。軋件采用線性減縮積分六面體單元（C3D8R）來離散（圖2），共離散為68 544個(gè)單元。

圖1 復(fù)合坯件（a）及有限元模型（b）Fig．1 Composite slab（a）and finite element model（b）

表2 復(fù)合坯件各部分參數(shù)Table 2 Parameter of composite slab

圖2 軋件對(duì)稱邊界及初始條件Fig．2 Symmetric boundaries and initial condition of work piece

2．2 邊界條件與初始條件

在整個(gè)軋制過程中，軋輥沿著其對(duì)稱軸做勻速定軸轉(zhuǎn)動(dòng)，軋輥的轉(zhuǎn)動(dòng)角速度設(shè)為10rad／s，則其邊緣線速度為1．0m／s；由于軋件建立的是1／4模型，底面與內(nèi)側(cè)面均為對(duì)稱邊界條件（圖2）；利用預(yù)定義場(chǎng)賦予軋件初始條件，定義的初速度為0．8m／s，其略小于軋輥的邊緣線速度，可保證軋件能順利咬入輥縫中，軋件進(jìn)入輥縫后靠軋輥與軋件之間的摩擦來帶動(dòng)軋件運(yùn)動(dòng)。

2．3 層間接觸屬性

本問題使用接觸對(duì)算法。軋件上表面與軋輥表面、芯體表面與框架及上蓋板下表面、框架上表面與上蓋板下表面之間均要設(shè)置接觸關(guān)系。接觸類型為顯式面面接觸，采用經(jīng)典的庫倫摩擦模型，其摩擦行為用摩擦系數(shù)μ來表征，臨界剪應(yīng)力τcrit取決于法向接觸壓應(yīng)力p：

只有在接觸面之間的剪應(yīng)力等于極限摩擦剪應(yīng)力μp時(shí)，接觸面之間才會(huì)發(fā)生相對(duì)滑動(dòng)，接觸面的摩擦系數(shù)設(shè)置為0．3。

3 計(jì)算結(jié)果與討論

3．1 鋯合金用戶材料子程序正確性的驗(yàn)證

圖3 鋯合金材料VUMAT子程序正確性的驗(yàn)證Fig．3 Verification of VUMAT subroutine for zircaloy

圖4 蓋板寬度方向位移云圖（a）、外側(cè)面3條路徑（b）和3條路徑上的位移分布（c）Fig．4 Contour plot of through－width displacements in cover（a），three paths at outer surface（b）and displacement distributions along three paths（c）

圖5 蓋板位于寬度對(duì)稱面上的兩條路徑（a）和軋制后兩個(gè)路徑上各點(diǎn)沿厚度方向的位移（b）Fig．5 Two paths on symmetric plane along width direction in cover（a）and through－thickness displacements on two paths after rolling（b）

3．2 工件變形及受力分析

通過數(shù)值模擬能得到復(fù)合坯件各部分軋制之后的形狀和尺寸，下面主要給出蓋板在寬度方向和厚度方向的變形情況。

1）蓋板寬度方向的變形

圖4a示出了有限元模型中上蓋板在寬度方向的位移云圖，可看出，蓋板外側(cè)面具有非均勻分布的寬展。取圖4b所示的3條路徑，繪出其軋制后沿寬度的位移（圖4c）。由圖4c可知，3條路徑上對(duì)應(yīng)的位移沿路徑具有類似的變化規(guī)律，最大的寬展位于路徑3的兩端，即靠近復(fù)合坯中面部分的兩端處，其寬展均經(jīng)歷了由大變小再變大的情況。從具體數(shù)值看，蓋板各處的寬展并不大，最大的寬度變化為0．65mm左右。

2）工件厚度方向的變形

復(fù)合坯各部分軋制后的厚度變化備受關(guān)注，如包殼的厚度直接關(guān)系到單片式板元件在堆內(nèi)使用的安全性。在有限元模型的上蓋板上選取如圖5a所示的兩條路徑，繪出路徑上各點(diǎn)沿厚度方向的位移（圖5b）。根據(jù)路徑5上的對(duì)應(yīng)位移，可看出其下面芯體厚度的變化；而根據(jù)路徑4和5的位移差值，可計(jì)算得到蓋板位于寬度對(duì)稱面處厚度的變化。

由圖5b可看出，芯體和蓋板的厚度變化沿路徑呈非均勻分布；由路徑4、5的結(jié)果可知，蓋板上表面軋制后相對(duì)較平整，而芯體與蓋板接觸的位置表面不夠平整，這是由于端部在軋制過程中蓋板與芯體之間產(chǎn)生的復(fù)雜力學(xué)相互作用；相對(duì)接近芯體的部分，蓋板上表面產(chǎn)生了很大的位移量?？傻贸?，軋制之后，蓋板的厚度變化明顯高于芯體的厚度變化。

3．3 接觸面的接觸壓力

通過軋制模擬還可得到軋制過程中芯體與復(fù)合坯件其他部分相互接觸表面上的接觸壓力，接觸壓力決定了各部分之間是否能達(dá)到冶金結(jié)合。下面主要給出芯體上表面與蓋板下表面之間的接觸壓力在軋制過程中的演化情況。

選取芯體表面上沿寬度的1條路徑，如圖6a所示，繪出該路徑上各點(diǎn)在不同軋制時(shí)間下的接觸壓力分布（圖6b）。幾個(gè)不同時(shí)刻分別代表幾個(gè)典型階段，即復(fù)合坯件入輥階段，軋輥接近和遠(yuǎn)離所考察路徑階段。由圖6b可看出，在復(fù)合坯件剛?cè)胼侂A段，路徑6附近的芯體與蓋板上表面尚未產(chǎn)生力學(xué)相互作用，接觸壓力為零；隨著軋輥的接近，芯體上表面與蓋板下表面之間產(chǎn)生了較為強(qiáng)烈的力學(xué)相互作用，最大接觸壓力產(chǎn)生于路徑上接近于寬度方向上的對(duì)稱面一端，沿著路徑逐漸減少，在靠近寬度方向的外側(cè)處最低；隨著軋輥遠(yuǎn)離該路徑，接觸壓力隨之降低。如果知道達(dá)到冶金結(jié)合所需達(dá)到的接觸壓力，可根據(jù)接觸表面上的接觸壓力情況來判斷實(shí)現(xiàn)冶金結(jié)合的程度和情況。

圖6 芯體上表面的輸出路徑6（a）和路徑6上各點(diǎn)在不同軋制時(shí)間下的接觸壓力分布（b）Fig．6 Path 6at upper surface of fuel meat（a）and distribution of contact pressure along path 6at different rolling time（b）

4 結(jié)論

本研究針對(duì)UMo單片式核燃料板復(fù)合坯件，建立了對(duì)其進(jìn)行軋制模擬的有限元模型；針對(duì)包殼材料應(yīng)變率相關(guān)的本構(gòu)關(guān)系，推導(dǎo)出適于顯式動(dòng)力學(xué)模擬的三維應(yīng)力更新算法，并編制了用戶材料子程序VUMAT，實(shí)現(xiàn)了軋制過程的有限元模擬，對(duì)軋制結(jié)果進(jìn)行了初步分析。主要結(jié)論如下：

1）用戶材料子程序VUMAT能正確地定義應(yīng)變率相關(guān)的材料本構(gòu)關(guān)系，為正確有效地模擬軋制過程，以及優(yōu)化工藝參數(shù)奠定了基礎(chǔ)。

2）蓋板外側(cè)面的寬展呈現(xiàn)非均勻分布，在沿長(zhǎng)度方向的兩端處具有最大的寬展量，從具體數(shù)值上來看非常小。

3）軋制后蓋板的相對(duì)厚度變化低于芯體，這是由不同的材料特性和受力情況決定的。

4）芯體上表面與蓋板下表面之間的接觸壓力隨時(shí)間變化，在靠近寬度方向的對(duì)稱面處具有最大的接觸壓力，向外側(cè)遞減。

［1］ OZALTUN H，SHEN M H，MEDVEDEV P．Assessment of residual stresses on U10Mo alloy based monolithic mini－plates during hot isostatic pressing［J］．Journal of Nuclear Materials，2011，419（1－3）：76－84．

［2］ KEISER D，RICE F．CNEA fresh fuel plate characterization summary report［R］．USA：Idaho National Laboratory，2012．

［3］ 劉建生，陳慧琴，郭曉霞．金屬塑性加工有限元模擬技術(shù)與應(yīng)用［M］．北京：冶金工業(yè)出版社，2003．

［4］ BAY N，CLEMENSEN C，JUELSTORP O．Bond strength in cold roll bonding［J］．CIRP Annals－Manufacturing Technology，1985，34（1）：221－224．

［5］ VAIDYANATH L R，NICHOLAS M G，MILNER D R．Pressure welding by rolling［J］．British Welding Journal，1959，6：13－28．

［6］ CAVE J A，WILLIAMS J D．Mechanism of cold pressure welding by rolling［J］．Journal of the Institute of Metals，1973，101（7）：203－207．

［7］ BAY N．Mechanisms producing metallic bonds in cold welding［J］．Welding Research Supplement，1983，62（5）：137－142．

［8］ ERSHOV A A，MYL’NIKOV A S，SYCHEVA T A．Conditions of joint formation in dissimilar metals during combined cold rolling［J］．Mechanics of Composite Materials，1985，21（1）：84－87．

［9］ CLEMENSEN C，JUELSTORP O，BAY N．Cold welding，Part 3：Influence of surface preparation on bond strength［J］．Metal Construction，1986，18（6）：625－629．

［10］ZHANG W，BAY N．Cold welding：Experimental investigation of the surface preparation methods［J］．Welding Research Supplement，1997，76（8）：326－330．

［11］SHTERTSER A A．Welding wave on the contact spot of solids［J］．Tribology International，1998，31（4）：169－174．

［12］STEINLAGE G A，BOWMAN K J，TRUMBLE K P．A comparison of three mechanical models for cold roll bonding metal laminates［J］．Journal of Adhesion Science and Technology，1996，10（3）：199－229．

［13］MANESH H D，SHAHABI H S．Effective parameters on bonding strength of roll bonded Al／St／Al multilayer strips［J］．Journal of Alloys and Compounds，2009，476（1－2）：292－299．

［14］ABBASI M，TOROGHINEJAD M R．Effects of processing parameters on the bond strength of Cu／Cu roll－bonded strips［J］．Journal of Materials Processing Technology，2010，210（3）：560－563．

［15］JAMAATI R，TOROGHINEJAD M R．Investigation of the parameters of the cold roll bonding（CRB）process［J］．Materials Science and Engineering A，2010，527（9）：2 320－2 326．

［16］YAN H，LENARD J G．A study of warm and cold roll－bonding of an aluminium alloy［J］．Materials Science and Engineering A，2004，385（1－2）：419－428．

［17］YU H，LIU X，BI H，et al．Deformation behavior of inclusions in stainless steel strips during multi－pass cold rolling［J］．Journal of Materials Processing Technology，2009，209（1）：455－461．

［18］CHAUDHARI G P，ACOFF V．Cold roll bonding of multi－layered Bi－metal laminate composites［J］．Composites Science and Technology，2009，69（10）：1 667－1 675．

［19］ANGELOV T，NEDEV A．Finite－element analysis of a hot－rolling problem with nonlinear friction［J］．Advances in Engineering Software，1997，28（9）：555－561．

［20］李冠興，武勝．核材料科學(xué)與工程：核燃料［M］．北京：化學(xué)工業(yè)出版社，2007．

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

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

Numerical Simulation Research on Rolling Process of Monolithic Nuclear Fuel Plate

WAN Ji－bo，KONG Xiang－zhe，DING Shu－rong＊，XU Hong－bin，HUO Yong－zhong
（Department of Mechanics and Engineering Science，F(xiàn)udan University，Shanghai 200433，China）

For the strain－rate－dependent constitutive relation of zircaloy cladding in UMo monolithic nuclear fuel plates，the three－dimensional stress updating algorithm was derived out，and the corresponding VUMAT subroutine to define its constitutive relation was developed and validated；the finite element model was built to simulate the frame rolling process of UMo monolithic nuclear fuel plates；with the explicit dynamic finite element method，the evolution rules of the deformation and contact pressure during the rolling process within the composite slab were obtained and analyzed．The research results indicate that it is convenient and efficient to define the strain－ratedependent constitutive relations of materials with the user－defined material subroutine VUMAT；the rolling－induced contact pressure between the fuel meat and the covers varies with time，and the maximum pressure exits at the symmetric plane along the plate width direction．This study supplies a foundation and a computation method for optimi－zing the processing parameters to manufacture UMo monolithic nuclear fuel plates．

monolithic nuclear fuel plate；rolling；numerical simulation；large deformation elasto－plasticity；strain－rate－dependent

TB125

：A

：1000－6931（2015）03－0511－07

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

2014－01－08；

2014－02－27

國家自然科學(xué)基金資助項(xiàng)目（11172068，91226101，11272092）

萬繼波（1987—），男，湖北洪湖人，碩士研究生，工程力學(xué)專業(yè)

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