单片式核燃料板轧制过程的数值模拟研究

万继波，孔祥喆，丁淑蓉，徐宏彬，霍永忠

（复旦大学力学与工程科学系，上海 200433）

单片式核燃料板轧制过程的数值模拟研究

万继波，孔祥喆，丁淑蓉＊，徐宏彬，霍永忠

（复旦大学力学与工程科学系，上海 200433）

针对UMo合金单片式核燃料板锆合金包壳材料应变率相关的力学本构关系，推导出其三维应力更新算法，相应地编写了定义其本构关系的VUMAT子程序并验证了程序的正确性；建立了对UMo合金单片式板元件的框架轧制过程进行计算模拟的有限元模型；利用显式动力有限元法，计算分析了复合坯内部的变形以及接触压强在轧制过程中的演化规律。研究结果表明，利用VUMAT用户材料子程序能方便正确地定义材料应变率相关的本构关系；燃料芯体与盖板之间的轧制接触压力随时间而演化，在靠近宽度方向的对称面处具有最大的接触压力。本研究为优化UMo合金单片式核燃料板的制造工艺参数提供了理论基础和计算手段。

单片式核燃料板；轧制；数值模拟；大变形弹塑性；应变率相关

UMo合金燃料［1］具有铀密度高、辐照稳定性好等优点，UMo合金单片式核燃料板在高通量研究试验堆中具有良好的发展前景。框架轧制复合法是一种有效制备UMo合金单片式核燃料板的工艺方法［2］。结合工艺试验，利用有限元法模拟轧制过程，成为优化工艺参数的发展趋势［3］。

UMo合金单片式核燃料板的制造质量，如芯体与其他部分的结合质量，与其堆内能密切相关。根据轧制结合的理论［4－12］，只有界面接触压力或等效应变达到一定程度才能使不同金属在界面处达到冶金结合，而结合强度受到多个工艺参数的影响［13－16］。为了提高燃料板的加工质量，需优化其加工工艺参数。对于轧制过程的有限元模拟，目前大部分的研究均针对单片板进行［17－19］，虽然存在3层板轧制的二维解析分析［4］，但对于燃料板复合坯这种具有复杂结构工件的轧制模拟研究尚未见报道。

本工作推导锆合金材料应变率相关的三维应力更新算法，编写定义材料力学本构关系的VUMAT子程序，验证其正确性；针对UMo单片燃料板的框架轧制工艺，建立对单片燃料板的轧制过程进行数值模拟的有限元模型；利用显式动力学模拟方法，计算分析燃料板复合坯内部的变形以及接触压强在轧制过程中的演化规律，为板形核燃料元件的先进制造提供理论基础和计算手段。

1 理论模型

1．1 材料属性

1）UMo合金材料性能

本文主要对复合坯在800℃下的轧制过程进行模拟，此温度下含Mo量为10%（质量分数）的材料性能参数［1，20］列于表1。

表1 UMo合金材料性能参数Table 1 Material performance parameter of UMo alloy

2）锆合金材料的材料属性

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

式中：E为未受辐照时锆合金的弹性模量，Pa；T为温度，K。

锆合金未受辐照时的应力应变曲线为：

1．2 应变率相关塑性材料的三维应力更新算法

如式（3）所示，锆合金材料具有与应变率相关的硬化曲线，要实现燃料板复合坯轧制过程的数值模拟，需导出其对应的三维应力更新算法，编制用户材料子程序，方便引入轧制的有限元分析中。

首先，假设该时间步无新的塑性应变产生，得到试验应力：

其中，G和λ为材料的拉梅常数。

求解试验应力σprij所对应的Mises应力为：

对于有新塑性应变产生的情况，需考虑塑性应变增量对应力的影响。根据塑性本构理论，对微小塑性应变增量dεpij［22］，有：

对式（8）右边进行后欧拉积分，得到该时间步的塑性应变增量：

则t＋Δt时刻的柯西应力为：

其偏斜部分为：

则：

在小学数学教学的过程中，习题练习能够有效地帮助学生巩固知识记忆，同时也能够帮助教师了解学生的学习难点，从而能够有针对性地进行讲解。而在此过程中，教师也可以结合微课教学视频开展复习工作，提升学生的学习效率。

对式（12）两端自乘，得：

可得：

即：

将式（17）代入式（16），得到Δ¯εp所满足的非线性方程：

收敛条件为：

其中：

按照上述算法，编制了用户材料子程序VUMAT，并引入轧制过程的有限元计算中。

2 有限元模型

2．1 几何模型

复合坯件（图1a）由锆合金上下盖板、框架和UMo单片式燃料芯体装配成整体，并在四周焊接成多材质复合坯件［2］，经过多道次的热轧和冷轧工艺，最后实现芯体与包壳之间的冶金结合，达到预定的形状和尺寸。轧制过程中轧辊与复合坯件之间产生相互作用，其具体参数列于表2［3］。由于对称性，可取整个工件的1／4和一个轧辊的1／2建立有限元模型（图1b）。轧辊与可产生较大弹塑性变形的轧板相比，其变形量极小，在此用解析刚体壳来描述。轧件采用线性减缩积分六面体单元（C3D8R）来离散（图2），共离散为68 544个单元。

图1 复合坯件（a）及有限元模型（b）Fig．1 Composite slab（a）and finite element model（b）

表2 复合坯件各部分参数Table 2 Parameter of composite slab

图2 轧件对称边界及初始条件Fig．2 Symmetric boundaries and initial condition of work piece

2．2 边界条件与初始条件

在整个轧制过程中，轧辊沿着其对称轴做匀速定轴转动，轧辊的转动角速度设为10rad／s，则其边缘线速度为1．0m／s；由于轧件建立的是1／4模型，底面与内侧面均为对称边界条件（图2）；利用预定义场赋予轧件初始条件，定义的初速度为0．8m／s，其略小于轧辊的边缘线速度，可保证轧件能顺利咬入辊缝中，轧件进入辊缝后靠轧辊与轧件之间的摩擦来带动轧件运动。

2．3 层间接触属性

本问题使用接触对算法。轧件上表面与轧辊表面、芯体表面与框架及上盖板下表面、框架上表面与上盖板下表面之间均要设置接触关系。接触类型为显式面面接触，采用经典的库伦摩擦模型，其摩擦行为用摩擦系数μ来表征，临界剪应力τcrit取决于法向接触压应力p：

只有在接触面之间的剪应力等于极限摩擦剪应力μp时，接触面之间才会发生相对滑动，接触面的摩擦系数设置为0．3。

3 计算结果与讨论

3．1 锆合金用户材料子程序正确性的验证

图3 锆合金材料VUMAT子程序正确性的验证Fig．3 Verification of VUMAT subroutine for zircaloy

图4 盖板宽度方向位移云图（a）、外侧面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 盖板位于宽度对称面上的两条路径（a）和轧制后两个路径上各点沿厚度方向的位移（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 工件变形及受力分析

通过数值模拟能得到复合坯件各部分轧制之后的形状和尺寸，下面主要给出盖板在宽度方向和厚度方向的变形情况。

1）盖板宽度方向的变形

图4a示出了有限元模型中上盖板在宽度方向的位移云图，可看出，盖板外侧面具有非均匀分布的宽展。取图4b所示的3条路径，绘出其轧制后沿宽度的位移（图4c）。由图4c可知，3条路径上对应的位移沿路径具有类似的变化规律，最大的宽展位于路径3的两端，即靠近复合坯中面部分的两端处，其宽展均经历了由大变小再变大的情况。从具体数值看，盖板各处的宽展并不大，最大的宽度变化为0．65mm左右。

2）工件厚度方向的变形

复合坯各部分轧制后的厚度变化备受关注，如包壳的厚度直接关系到单片式板元件在堆内使用的安全性。在有限元模型的上盖板上选取如图5a所示的两条路径，绘出路径上各点沿厚度方向的位移（图5b）。根据路径5上的对应位移，可看出其下面芯体厚度的变化；而根据路径4和5的位移差值，可计算得到盖板位于宽度对称面处厚度的变化。

由图5b可看出，芯体和盖板的厚度变化沿路径呈非均匀分布；由路径4、5的结果可知，盖板上表面轧制后相对较平整，而芯体与盖板接触的位置表面不够平整，这是由于端部在轧制过程中盖板与芯体之间产生的复杂力学相互作用；相对接近芯体的部分，盖板上表面产生了很大的位移量。可得出，轧制之后，盖板的厚度变化明显高于芯体的厚度变化。

3．3 接触面的接触压力

通过轧制模拟还可得到轧制过程中芯体与复合坯件其他部分相互接触表面上的接触压力，接触压力决定了各部分之间是否能达到冶金结合。下面主要给出芯体上表面与盖板下表面之间的接触压力在轧制过程中的演化情况。

选取芯体表面上沿宽度的1条路径，如图6a所示，绘出该路径上各点在不同轧制时间下的接触压力分布（图6b）。几个不同时刻分别代表几个典型阶段，即复合坯件入辊阶段，轧辊接近和远离所考察路径阶段。由图6b可看出，在复合坯件刚入辊阶段，路径6附近的芯体与盖板上表面尚未产生力学相互作用，接触压力为零；随着轧辊的接近，芯体上表面与盖板下表面之间产生了较为强烈的力学相互作用，最大接触压力产生于路径上接近于宽度方向上的对称面一端，沿着路径逐渐减少，在靠近宽度方向的外侧处最低；随着轧辊远离该路径，接触压力随之降低。如果知道达到冶金结合所需达到的接触压力，可根据接触表面上的接触压力情况来判断实现冶金结合的程度和情况。

图6 芯体上表面的输出路径6（a）和路径6上各点在不同轧制时间下的接触压力分布（b）Fig．6 Path 6at upper surface of fuel meat（a）and distribution of contact pressure along path 6at different rolling time（b）

4 结论

本研究针对UMo单片式核燃料板复合坯件，建立了对其进行轧制模拟的有限元模型；针对包壳材料应变率相关的本构关系，推导出适于显式动力学模拟的三维应力更新算法，并编制了用户材料子程序VUMAT，实现了轧制过程的有限元模拟，对轧制结果进行了初步分析。主要结论如下：

1）用户材料子程序VUMAT能正确地定义应变率相关的材料本构关系，为正确有效地模拟轧制过程，以及优化工艺参数奠定了基础。

2）盖板外侧面的宽展呈现非均匀分布，在沿长度方向的两端处具有最大的宽展量，从具体数值上来看非常小。

3）轧制后盖板的相对厚度变化低于芯体，这是由不同的材料特性和受力情况决定的。

4）芯体上表面与盖板下表面之间的接触压力随时间变化，在靠近宽度方向的对称面处具有最大的接触压力，向外侧递减。

［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］ 刘建生，陈慧琴，郭晓霞．金属塑性加工有限元模拟技术与应用［M］．北京：冶金工业出版社，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］李冠兴，武胜．核材料科学与工程：核燃料［M］．北京：化学工业出版社，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，Fudan 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

国家自然科学基金资助项目（11172068，91226101，11272092）

万继波（1987—），男，湖北洪湖人，硕士研究生，工程力学专业

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