身管精锻过程跨尺度多晶体塑性有限元模拟与织构预测

徐笑，樊黎霞，王亚平，董晓彬

（南京理工大学机械工程学院，江苏南京210094）

身管精锻过程跨尺度多晶体塑性有限元模拟与织构预测

徐笑，樊黎霞，王亚平，董晓彬

（南京理工大学机械工程学院，江苏南京210094）

径向精锻是一种先进的身管制造方法。身管的多晶体材料经锻打后，晶粒会发生择优取向形成织构，在宏观上产生各向异性现象，这是普通有限元方法所无法模拟的。利用数学软件Matlab和程序设计语言Python联合编程并结合有限元软件ABAQUS，建立身管精锻宏观模型和细观多晶体模型；应用跨尺度方法，开发程序实现了宏观力学至微观力学边界条件的继承，使研究的晶粒尺寸达到了真实尺寸。分析微观多晶体模型的晶粒取向变化，使用Matlab编程画出锻后多晶体材料的极图和取向分布函数图，并用X射线衍射试验验证了晶粒取向模拟结果的正确性。建立了将身管径向精锻从宏观锻造到细观织构联系起来的一套研究方法，为预测锻后身管性能变化提供了思路。

兵器科学与技术；身管径向锻造；有限元分析；晶体塑性模型；织构

0 引言

身管是自动武器最重要的部件之一，在射击中承受着高温、高压的火药气体的冲刷和弹丸的剧烈摩擦。与传统身管制造工艺相比，径向精锻技术自动化程度高，提高了线膛和弹膛的同轴度，降低了成本并且提高了寿命，很大程度提高了自动武器的性能。有限元分析方法是身管径向精锻研究的一种常用方法。Tszeng等［1］最早建立了径向锻造的轴对称模型。Domblesky等［2］使用ABAQUS有限元软件对径向锻造过程进行了数值模拟，并且研究了锻件进给速度和锻件不同点的等效塑性应变对其温度的影响。在国内，樊黎霞等［3］使用ABAQSU有限元软件模拟了身管的径向锻造过程，对径向锻压力及应力进行了数值分析，并用试验验证了模拟结果。刘力力等［4］建立了身管径向锻造及回弹过程的轴对称模型，并分析了不同工艺参数对锻后残余应力分布的影响。

大量研究表明，在金属塑性变形中晶体取向对变形和断裂等行为是敏感的。变形后多晶体的择优取向称为织构。Bunge［5］提出织构影响材料性能的比例达到了20%～50%之高。这就为研究身管锻后性能提供了新思路，以微观织构的研究预测宏观精锻产品的锻后性能。在板材冲压领域已有类似先进的研究，马全仓等［6］根据取向、织构和简单的滑移理论预测了3104铝合金板的初始塑性应变比，并且计算结果与实测值变化趋势一致。陈攀宇等［7］利用Deform-3D软件模拟和试验研究AZ31镁合金管轧制过程，并使用金相显微镜对轧制前后织构进行分析，表明多道次轧制使得晶粒再结晶和细化。晶体塑性有限元能够模拟材料的各向异性行为，此方法更接近材料的本质特性。Balasubrmanian等［8］研究了多晶体塑性有限元在优化深冲工艺的可行性。Zhou等［9］用晶体塑性有限元研究了面心立方材料深冲制耳与织构的关系。在国内，李大永等［10］以轧制铝板为研究对象，利用塑性有限元方法对方盒冲压成形和织构演化进行了数值模拟。

跨尺度分析方法作为解决宏观和微观耦合问题的新方法在多晶体塑性领域已经得到广泛关注。郭一娜等［11］提出一种宏观和微观尺度结合的建模方法来研究铸态42CrMo环件的铸辗复合成形工艺。在国外，Bronkhorst等［12］在研究材料受大变形和不同热力的微观结构的变化时，建立了钽的基于状态变量的宏观内部连续模型，并提出了一种多尺度的多晶体塑性模型，对比两种模型在进行剪切实验时剪切区的VonMises应力，等效塑性应变等。在国内，杜凤山等［13］针对Cr5锻件锻压过程中晶体间的滑移进行了多尺度有限元仿真，得到了宏观滑移趋势和细观滑移趋势的对比。王国春［14］利用有限元重合网格方法开发程序，并对带圆孔的金属板等进行多尺度应力分析。

本文将以材料的本质属性入手研究身管径向精锻，结合多晶体塑性有限元方法和跨尺度方法，建立一种将宏观身管精锻过程和微观织构联系起来的建模分析方法，为研究身管径向精锻的织构与宏观性能的定量关系打下基础。结合前人的理论框架和经验，采用非线性分析功能强大的ABAQUS有限元软件和Matlab程序语言、Python脚本编程联合建模，分别建立身管精锻过程的轴对称宏观模型和微观多晶体有限元模型，利用线性插值理论实现宏观至微观边界条件的继承。最后，编译程序提取所需晶粒取向绘制极图和取向分布函数（ODF）图并与X射线衍射（XRD）试验所测的实际图形作对比，验证本方法的可行性。

1 精锻过程宏观仿真模拟

身管的径向精锻加工是一个较为复杂的金属材料塑性流动过程。毛坯的底部由精锻机的夹头加持，并给予其一定的夹持力以实现轴向进给。毛坯的另一端与挡块配合。毛坯四周对称分布了4个锤头，加工时4个锤头对毛坯进行高频率的锻打。毛坯内部置有芯棒，芯棒具有身管线膛和弹膛的特征。为了保证圆形的横截面和较高的表面光滑度，毛坯与芯棒在两次锻打间隙一起做旋转进给。

锤头的速度对材料应变率及再结晶有较大影响，在精锻加工中，此速度为一谐函数，如（1）式所示:

式中:v为锤头速度（mm/s）；n为锻打频率（次/min）；h为距上止点距离（mm）；s为锤头全行程（mm）。

1.1 轴对称有限元模型建立

运用ABAQUS对身管的径向锻造过程建立有限元模型。身管径向锻造在忽略内部膛线、身管旋转轴向进给和锤头的闭合间隙的情况下，是一个轴对称的加工过程建立身管径向锻造的轴对称模型。重要部件几何尺寸如表1所示，身管材料力学性能参数如表2所示。

表1 身管毛坯几何尺寸Tab.1 Geometric sizes of barrel workpiece

表2 材料力学性能Tab.2 Mechanical properties of material

身管的整个锻造过程为冷精锻，挡块、锤头和芯棒视为刚体。身管毛坯材料为30SiMn2MoVA，对其作网格划分，采用4节点轴对称缩减积分单元CAX4R，力学性能使用应变硬化模型σ= 1 225ε0.065，σ为流变应力，ε为应变，所建立宏观模型图1所示。

图1 宏观轴对称有限元模型Fig.1 Macroscopic axisymmetric FEA model

1.2 仿真结果

在ABAQUS/CAE中设定上述参数后，于Job模块中创建分析并求解，得到宏观模型的Mises等效应力分布云图（见图2）。锻压力是精锻实际生产中非常重要的一个工艺参数，加工时不能超过机器的最大锻压力值，所以在图2中选取内表面应力较大的一点作为研究对象。

图2 Mises应力布云图Fig.2 Mises stress nephogram

2 微观尺度上多晶体材料建模仿真

微观尺度多晶体建模采用ABAQUS/CAE建模与INP文件建模两种方法结合的方式。INP文件包含了对整个模型的完整描述，在前处理器和求解器之间建立了一座数据桥梁。具体建模方法分为4步:1）建立与宏观模型所选单元尺寸相同的模型，划分网格，设置分析步等重要信息；2）使用Voronoi方法进行晶粒划分；3）将宏观模型中得到的节点位移信息通过程序进行插值计算，得到新的微观模型上4条边界上所有节点的位移作为新的边界条件；4）修改UMAT子程序信息，于求解器求解时提交。

图3 微观有限元模型建立流程Fig.3 Microscopic FEA modeling process

2.1 微观基础模型的建立

微观有限元模型的建立流程如图3所示。选取一个在图3（a）中编号为3 366的单元进行跨尺度分析，即为图3（b）的标识出的单元。使用ABAQUS/ CAE建立与宏观模型中与此单元尺寸完全一致的模型。所建模型单元尺寸为250 μm×249 μm，共划分为15 625个单元，即为图3（c）。为了达到完全模拟在宏观模型中的受力情况，单元采用与宏观模型一样的CAX4R单元。分析步等建模信息也与宏观模型设置同样的参数。

假设在锻造前，模型内的晶粒大小相似、取向随机并且没有取向择优的织构。这与Voronoi图形生成方法随机、均匀的思想相符合。所以采用该方法划分晶粒。使用Matlab软件编译程序，读取微观模型的基本信息，如网格划分和节点坐标等。Voronoi图形的结构受种子点的数量及坐标的影响非常大。为了使种子点的信息受均匀性的控制，在程序中编写控制函数，并调用Voronoi函数生成晶粒的拓扑信息，再使用inpolygon函数循环判断出每个单元中心点属于哪个晶粒，最后生成晶粒分组信息，并修改入模型的INPUT文件中。模型中的晶粒取向用米勒指数法来表示。

图4为微观多晶体模型建模的过程示意图，右图已经成功划分晶粒，其含有258个晶粒，平均晶粒尺寸为15.5 μm，符合金属晶粒的真实大小。

图4 微观多晶体有限元模型Fig.4 Microscopic polycrystalline FEA model

图5 单元顶点位移Fig.5 Displacement of element's vertex

图6 矩形单元广义坐标Fig.6 Rectangular element generalized coordinates

2.2 跨尺度边界条件的继承

为使宏观模型计算过程中该单元的边界条件继承至微观模型中，使用有限元线性插值理论，并结合材料内部的位移连续理论，进行宏观至微观跨尺度边界条件的继承。

从宏观模型计算结果中可以得到该单元的4个顶点在整个精锻过程中随时间变化的位移曲线。图5为单元第一顶点的x、y方向位移图。x轴的正方向为轴向进给方向，y轴的正方向为毛坯的径向。由于冷锻工艺成形是利用金属的塑性流动来实现的，当所取单元位置受到锻打后，金属会向进给方向的反方向流动，开始阶段流动速度大于进给速度，所以在开始阶段y轴会出现负值，随着所取单元通过轴向进给离锻打区域越来越远，其位移主要受轴向进给控制，所以位移值斜率会趋近于进给速度。

微观建模型共1 000个边界节点。使用二维4节点线性单元形函数和位移函数联合计算得出任意节点的位移曲线。首先提取出1 000个边界点的坐标值，如图6所示，先将其化成广义坐标，再代入（2）式～（6）式中，利用得到的广义坐标值计算每个点的4个原单插值函数Ni，再将4个顶点的位移值ui和4个单元插值函数Ni代入（7）式位移函数即可得到所需边界点的相应位移值｛r｝，具体公式为

式中:i（ξi，ηi）分别为1（1，1）、2（-1，1）、3（-1，-1）、4（1，-1）；x0和y0为自然坐标系的原点坐标；a、b为图6中矩形单元的长和宽。

采用Python脚本语言计算插值结果，程序生成INPUT格式的文本文件，将其修改入模型文件中，最后达到添加边界条件的目的。首先，读取微观模型的OBD文件，于根装配中提取模型边界的节点信息，包括节点编号和坐标，并将每个节点定义为一个单独的集合，最后以INPUT格式输出。具体格式为

读取宏观模型的OBD结果文件，使用（1）式的计算方法得到1 000个节点的随时间变化的位移曲线。在ABAQUS中添加约束，为上一步中定义的单个节点集合设置位移载荷约束，位移大小设置为1，再将求得的任意节点的时间位移数据以Amplitude方式加于约束中。具体格式为

如上例所示，约束类型为位移/旋转型，施加目标为set-1集合，即前面所创建的边界的单节点集合，后面的1代表约束方向为x方向，最后的一个1代表位移值为1.此约束受Amplitude控制，上例所示的Amplitude模块中左边一列为时间，单位为s，右边一列为位移，单位为mm.将上述生成的集合信息、约束信息和振幅信息添加入微观模型的INPUT文件中，完整的微观材料模型就建成了。

当然，需要指出的是，目前没有任何迹象表明有火星微生物在那里“游泳”。事实上，这样环境也不适宜生命的存在。火星南极冰层下的水温预计在华氏零下90度，远低于水的冰点。

2.3 运用UMAT实现多晶体自定义材料属性

微观多晶体模型的塑性本构关系使用ABAQUS中的子程序UMAT来实现。本文采用Huang［15］的UMAT程序框架，硬化模型采用Asaro模型。模型中的某些参数如初始屈服时的硬化模量h0、应变率敏感系数n、参考剪切应变率γ·0等材料系数很难准确测量到，即使个别方法如潜硬化实验等也只能得到单晶的个别材料参数。现存的文献中也没有人做过本文所用材料30SiMn2MoVA的这些材料属性的研究。所以本文采用拟合实验所得应力-应变曲线的方法来获得这些材料系数。

3 数值模拟结果与试验验证

3.1 微观模型模拟结果

进行划分晶粒、材料本构关系赋予和边界条件的继承后就能得到微观多晶体塑性有限元模型。提交解算器计算后得到微观模型的真实应变分布如图7所示。

图7 应变分布云图Fig.7 Strain nephogram

从图7中可以看出，模型内部真实应变分布严重的不均匀，在应变带之间交替分布着大小不同的小应变区域和应变集中区域，这是由于不同晶粒的晶粒取向和滑移系的差别在塑性变形过程中引发应变的累积，从而使得相邻的晶粒的应变也可能会有很大的差别，最后导致应变及其不均匀的分布。有些晶粒内部应变大小差别在5倍以上，这是传统的有限元模拟达不到的效果。应变的集中不仅仅出现在晶界附近，也能出现于晶粒的内部。由于晶粒本身的各向异性、几何形状和晶粒间的变形协调导致晶粒的形状变得更加的不规则。

3.2 极图绘制与试验验证

晶体学一般采用极图和ODF图来描述织构。身管的材料为30SiMn2MoVA，属于体心立方（BCC）晶体，具有48个滑移系，研究发现在冷塑性变形中只有一部分滑移系有大量的开动，主要开动的滑移系分别为12个｛110｝〈111〉滑移系和12个｛112｝〈111〉滑移系。编译Python脚本程序提取出结果文件中上述24个滑移系的密勒指数信息，并将其转换成欧拉角（φ1，φ，φ2）格式。使用Matlab的插件MTEX编译程序，读取上述欧拉角信息，绘制可以较好表达变形中晶粒取向变化的｛100｝方向极图和以φ1为视角的ODF图。

为了验证模拟结果的正确性，在5.8 mm口径锻后身管上切割了表面积为10 mm×10 mm正方形的试样进行了XRD试验。该试验从材料、身管尺寸及精锻加工参数等采用与本文一致的工艺参数。XRD分析是利用晶体形成的XRD，对物质进行内部原子在空间分布状况的结构分析方法。本试验采用德国Bruker公司产D8 Advance XRD仪。试验样品采用在锻后身管中采用线切割方法得到边长10 mm的块状材料，并在干燥环境下进行多次打磨，如图8所示。试验完成后得到XRD图谱和3个方向上的不完全极图两种数据文件。

图8 XRD试验试样Fig.8 XRD experimental samples

使用Jade5.0分析衍射图谱得到物相信息，再配合不完全极图3个不完全极图导入Multex软件处理得到相同方向的极图和ODF图。

图9 BCC金属取向线Fig.9 Orientations of BCC metal

对比观察模拟结果和试验结果所得到的极图（见图10）和ODF图（见图11），其强度较大的织构分布位置和对称性都较高的一致性。BCC金属变形纤维织构主要的织构类型有｛111｝〈UVW〉和｛HKL｝〈110〉两类，而这些织构大都处在α取向线和γ取向线。从DOF图中可以看出，在BCC中重要取向线α取向线上（0°，54.7°，45°）附近以及γ取向线上（90°，54.7°，45°）处有织构集中，如黑色十字标记所示。所以判定这两个位置有织构生成，将这两个位置的织构转化为密勒指数法表示，分别为｛111｝〈110〉和｛111｝〈112〉。

图10 极图比较Fig.10 Comparison of pole figures

图11 ODF图比较Fig.11 Comparison of ODF figures

将模拟结果通过Horta加权法计算可以得到各｛HKL｝的轴密度，同样XRD试验测得各｛HKL｝的衍射强度也可计算得到轴密度。将轴密度使用地质等高绘图软件Winsurf绘制反极图后可以得到模拟和试验中较强的织构强度如表3和表4所示。

表3 模拟所得较强织构的织构强度Tab.3 Simulated texture intensity of strong texture

表4 模拟所得较强织构的织构强度Tab.4 Experimental texture intensity of strong texture

可以得到30SiMn2MoVA材料身管径向精锻锻后织构成分为｛111｝〈112〉和｛111｝〈110〉，其中前者为主要织构成分。

4 结论

利用多晶体塑性有限元方法建立了身管精锻的宏观至微观多尺度有限元模型，使用线性插值理论编译程序实现了宏观和微观两个模型边界条件的统一，使研究的晶粒尺寸达到了实际大小。绘制极图和ODF图并与试验结果对比，将宏观锻造和微观织构联系起来。为以后研究宏观性能和微观织构的关系提供了一套研究方法。

1）基于多晶体塑性有限元方法和跨尺度方法，开发程序，建立从宏观到微观的跨尺度身管多晶体塑性有限元模型，可以更好地反映身管精锻中的各向异性现象，使身管精锻的数值模拟更加符合实际。

2）基于所建立的跨尺度模型分析微观织构，提取分析数据绘制极图和ODF图并与XRD试验所得图形对比具有较好的一致性。

3）通过多晶体塑性有限元法和跨尺度方法的结合，不仅可以更加准确地模拟金属宏观变形，而且能够预测塑性变形中的织构变化。

（References）

［1］ Tszeng T C，Kobayashi S.Determination of residual stresses in radial forging，manufacturing processes simulation，PED-20［R］. Berkeley，US:Department of Mechanical Engineering，University of California，1986:31-45.

［2］ Domblesky J P，Shivpuri R，Painter B.Application of the finite-element method to the radial forging of large diameter tubes［J］. Journal of Materials Processing Technology，1995，49（1）:57-74.

［3］ 樊黎霞，刘力力，刘庆东，等.身管线膛精锻加工过程的数值分析［J］.兵工学报，2009，30（8）:1098-1102. FAN Li-xia，LIU Li-1i，LIU Qing-dong，et al.The numerical simulation of the precision radial forging processing of the gun barrel［J］.Acta Armamentarii，2009，30（8）:1098-1102.（in Chinese）

［4］ 刘力力，樊黎霞，董雪花.盲孔法测定径向锻造身管表面残余应力方法研究［J］.兵工学报，2012，33（6）:712-717. LIU Li-li，FAN Li-xia，DONG Xue-hua.Study on measurement of residual stresses in forged barrel surface［J］.Acta Armamentarii，2012，33（6）:712-717.（in Chinese）

［5］ Bunge H J.Texture analysis in materials science mathematical methods［M］.London，UK:Butterworth Heinemann Ltd，1982.

［6］ 马全仓，毛卫民，冯惠平.3104深冲铝板织构对初始R值的影响［J］.北京科技大学学报，2004，26（1）:79-81. MA Quan-cang，MAO Wei-min，FENG Hui-ping.Influence of texture on initial R-value of deep drawing 3104 aluminum sheet［J］.Journal of University of Science and Technology Beijing，2004，26（1）:79-81.（in Chinese）

［7］ 陈攀宇，李彩霞，李立州，等.基于Deform-3D的镁合金管轧制成形数值模拟与试验研究［J］.轻合金加工技术，2015，43（7）:33-37. CHEN Pan-yu，LI Cai-xia，LI Li-zhou，et al.Numerical simulation and experimental study of roll forming of magnesium alloy tubes based on Deform-3D［J］.Light Alloy Fabrication Technology，2015，43（7）:33-37.（in Chinese）

［8］ Balasubramanian S，Anand L.Single crystal and polycrystal elasto-viscoplasticity:application to earing in cup drawing of FCC materials［J］.Computational Mechanics，1996，17（4）:209-225.［9］ Zhou Y，Jonas J J，Savoie J，et al.Effect of texture on earing in FCC metals:finite element simulations［J］.International Journal of Plasticity，1998，14（1/2/3）:113-138.

［10］ 李大永，张少睿，彭颖红，等.板材冲压成形的晶体塑性有限元模拟［J］.机械工程学报，2008，44（1）:190-194. LI Da-yong，ZHANG Shao-rui，PENG Ying-hong，et al.Finite element simulation of sheet metal stamping with polycrystalline plasticity［J］.Journal of Mechanical Engineering，2008，44（1）: 190-194.（in Chinese）

［11］ 郭一娜，丁双凤，李永堂，等.42CrMo环件铸辗复合成形的跨尺度建模［J］.机械工程学报，2014，50（14）:81-88. GUO Yi-na，DING Shuang-feng，LI Yong-tang，et al.Multiscale modeling for 42CrMo ring during blank-casting and rolling compound forming process［J］.Journal of Mechanical Engineering，2014，50（14）:81-88.（in Chinese）

［12］ Bronkhorst C A，Maudlin P J，Gray G T，et al.Accounting for microstructure in large deformation models of polycrystalline metallic materials［M］∥Ghosh S，Dimiduk D.Computational Methods for Microstructure Property Relationships.NY，US: Springer，2010:239-276.

［13］ 杜凤山，李源，王珉，等.锻压过程宏-细观跨尺度仿真研究［J］.计算力学学报，2014，31（6）:799-810. DU Feng-shan，LI Yuan，WANG Min，et al.Macro-meso crossscale simulation of forging process［J］.Chinese Journal of Computational Mechanics.2014，31（6）:799-810.（in Chinese）

［14］ 王国春.基于有限元重合网格法的结构多尺度分析研究［D］.新疆:新疆大学，2006. WANG Guo-chun.Study on multi-scale analysis of structure based on the finite element mesh superposition method［D］. Xinjiang:Xinjiang University，2006.（in Chinese）

［15］ Huang Y.A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program［R］.Cambridge，MA，US:Division of Engineering and Applied Sciences，Harvard University，1991:2-9.

Trans-scale Polycrystalline Finite Element Simulation of Radial Forging Process for Barrel and Prediction of Texture

XU Xiao，FAN Li-xia，WANG Ya-ping，DONG Xiao-bin
（School of Mechanical Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

Radial forging is an advanced manufacturing method for barrel.After the polycrystalline material of barrel is forged，the preferred orientation of the materials grain forms texture which presents the anisotropy phenomenon macroscopically.This phenomenon can not be simulated using the common finite element analysis（FEA）method.A macroscopic barrel radial forging FEA model and a microscopic polycrystalline plasticity FEA model are established based on ABAQUS code，the mathematical software Matlab and the programming language Python.The multi-scale method is used to program the boundary conditions inheriting from macro-mechanics state to meso-mechanics state which make the size of grain as reality.The pole figures of forged polycrystalline material and the orientation distribution function（ODF）figures are drawn by Matlab，which are used to analyze the orientation change of grain in the microscopic polycrystalline plasticity FEA model.The XRD experiment is used to validate the accuracy of the simulated results.A research method to connect macroscopic forging process and microscopic texture is presented，which provides a new way to predict the mechanical property of barrel after forging.

ordnance science and technology；radial forging of barrel；finite element analysis；crystalline plasticity model；texture

TJ05

A

1000-1093（2016）07-1180-07

10.3969/j.issn.1000-1093.2016.07.004

2015-08-06

国家自然科学基金项目（51575279）

徐笑（1989—），男，硕士研究生。E-mail:xuxiaoam@sina.com；

樊黎霞（1965—），女，教授，硕士生导师。E-mail:fanlixia151@njust.edu.cn