ZK61镁合金热轧板晶体塑性变形和织构演化机制

Crystal plastic deformation and texture evolution mechanism of ZK61 magnesium alloy hot-rolled plate

摘要

采用宏观有限元模拟与微观晶体塑性分析相结合的方法对六道次热轧ZK61镁合金中厚板的晶体内应力和织构演化进行研究。根据轧制实验与拉伸实验建立相应同步热轧模型并进行仿真计算，将模拟结果作为边界条件代入到基于Voronoi图的多晶体塑性模型中，再应用晶体塑性有限元方法对多晶体模型进行同步热轧模拟，得到各道次热轧态ZK61镁合金的塑性参数和织构极图，并分别与拉伸实验结果和电子背散射衍射实验导出的织构极图进行对比分析，归纳出不同道次下热轧ZK61镁合金晶体塑性变形与织构演化机制。结果表明：经多道次热轧后的ZK61板材存在大量孪晶与动态再结晶，晶粒均匀细化效果明显，其基面织构极密度与轧制道次整体呈正相关，不同道次下的峰值取向差角存在较大差异（8°、28°、88°），且合金力学性能得到改善，强度提升了7.55%，伸长率达到19.5%。

Abstract

The crystal internal stress and texture evolution of six-pass hot rolled ZK61 magnesium alloy plate are studied by means of combining macroscopic finite element simulation and microscopic crystal plasticity analysis. According to the rolling experiment and tensile test， the corresponding synchronous hot rolling model is established and simulated， and the simulation results are entered as boundary conditions in the polycrystalline plasticity model based on Voronoi diagram. Then the crystal plasticity finite element method is used to simulate the synchronous hot rolling of the polycrystalline model， and the plastic parameters and texture pole figures of the hot-rolled ZK61 magnesium alloy in each pass are obtained. Compared with the tensile test results and the texture pole figures derived from the electron backscatter diffraction experiment， the crystal plastic deformation and texture evolution mechanism of the hot-rolled ZK61 magnesium alloy under different passes are summarized. The results show that there are a large number of twins and dynamic recrystallization in the ZK61 plate after multi-pass hot rolling， the grain homogenization and refinement effect is obvious. The pole intensity of basal texture is correlated with the whole rolling passes， the peak misorientation angle of the alloy under different passes has a quite diversity （8°， 28°， and 88°）， and the mechanical properties of the alloy are improved. The alloy strength is increased by 7.55%， and the elongation is achieved by 19.5%. The results can provide reference for improving the plastic processing ability of magnesium alloy.

Graphical abstract

关键词

ZK61镁合金 / 热轧 / 数值模拟 / 晶体塑性有限元方法 / 微观结构 / 织构演化 / 力学性能

Key words

ZK61 magnesium alloy / hot rolling / numerical simulation / crystal plasticity finite element method / microstructure / texture evolution / mechanical property

作为“21世纪绿色工程材料”，镁合金在诸多热门及重要领域得到了广泛应用，但因其具有密排六方晶格（hexagonal close packed，HCP）的晶体结构，导致室温塑性变形能力较弱，通常需要采用热轧方式来启动锥面滑移和孪生变形^［1］。由于滑移和孪生的本质都是位错的剪切滑移，在研究HCP金属材料的塑性变形机制时，要明确在塑性变形过程中可能开动的孪晶，并将其当作一种滑移方式来处理，其中孪生产生的剪切应变所用本构方程与滑移的本构方程相同^［2-3］。另外热轧变形区内的微观晶粒状态参数对合金塑性变形能力也具有重要影响。

常见的研究塑性变形及相应晶体织构的非线性本构分析方法有黏塑性自洽（visco-plastic self-consistent，VPSC）模型^［4］与晶体塑性有限元方法^［5］（crystal plastic finite element method，CPFEM）等，其中，CPFEM可以很好地预测镁合金塑性变形过程的滑移、孪生等织构演化机制及晶间与晶粒内部的应力状态，以揭示材料的微观组织变化对力学性能的改善作用，且对织构的预测比VPSC模型更准确（VPSC模型通常过度预测织构^［6］），从而能够提供更有意义的结果。Pierce等^［7］在20世纪80年代首先提出了一种晶体塑性本构理论，将单晶塑性变形行为简化为两个特定的对称滑移系进行研究。该理论随后得到了发展完善并开始用于有限元模拟计算来预测多晶体的织构演化^［8-9］，被认为更能体现出金属材料的塑性变形本质^［10］。近年来对CPFEM的研究逐渐出现了许多新进展，如使用全场晶体塑性模型^［11］等改进方案进行多晶体有限元模拟；引入辐照效应^［12］、尺度效应^［13］等理论综合分析预测材料的微观演化机制；以及对于表面粗糙度^［14］、变形不均匀性^［15］与损伤失效机理^［16］等方面的研究应用。然而，现有研究大多仍处于在介观尺度上对于位错滑移与变形孪晶的激活与演变等金属塑性变形机制的分析与预测方面的优化补充，对于深入联系宏观、微纳尺度模拟分析的纵向结合上仍需进一步发展，以达到多尺度综合分析阐述演化机理与促进数学模型的创新开发的目的。

本工作建立了一种基于率相关晶体塑性本构方程的三维Voronoi多晶体应变硬化模型进行微观多晶体CPFEM计算，并结合宏观热轧ZK61镁合金中厚板有限元模拟进行多尺度分析，同时与轧件的拉伸试验和EBSD实验结果加以对比，得出热轧ZK61镁合金中厚板的塑性变形及织构演化机制。

1 实验材料与方法

选用规格为65 mm×200 mm×17 mm的ZK61镁合金铸态坯料，其具体成分见表1所示^［17］。

在对ZK61镁合金开展热轧实验前，为消除ZK61合金中晶界处的脆性Mg-Zn共晶相，需要对铸坯进行400 ℃×12 h的均匀化处理，其加热设备为LDK-A隧道炉，图1为合金经均匀化处理后的初始金相，其平均晶粒尺寸为66.67 μm。将铸坯由辽宁科技大学自研的中厚板实验轧机（图2）进行六道次同步热轧至4.5 mm，其中轧制温度为400 ℃，轧制速度为0.16 m/s，每道次间轧件需要保温10 min，板材轧前加热与道次间加热均使用RX3-45-12箱式电阻炉。在轧前和每道次轧制后均需取样。

合金拉伸试样如图3所示，分别在轧制方向（rolling direction，RD）、45°斜向与轧制横向（transeverse direction，TD）制备第三到第六道次轧件的1 mm厚标准拉伸试样，并在由800^#金相砂纸将表面打磨光亮后，使用WDW3100微机控制电子万能试验机，参照GB/T228—2002《金属材料室温拉伸试验方法》进行室温单向拉伸实验，其沿RD拉伸速率为1 mm/min。重复两次实验并分别取各相应位置拉伸数据的平均值以减小误差，随后将其导入到数据分析处理软件Origin中绘制试样的拉伸应力-应变曲线图。

将不同轧制状态下的EBSD试样在打磨和机械抛光后需要电解抛光，其电解液成分为20 mL高氯酸与180 mL无水乙醇，并倒入液氮降温至-30 ℃，在15 V电压下电解150 s，随后使用YM-060S超声波清洗机在丙酮中清洗2 min后吹干。试样采用ZEISS Sigma-300热场发射扫描电子显微镜进行EBSD实验，并运用Channel 5与MATLAB软件处理实验数据，得到晶粒反极图（inverse pole figure，IPF）、取向差角分布统计和织构极图。

2 有限元模型的建立

2.1 宏观轧制模型

采用商用软件ABAQUS对ZK61镁合金板材多道次同步热轧进行有限元建模仿真分析。为兼顾模拟的精确性和可操作性，需要对轧制模型进行如下假定：①模型为均匀对称轧制，上下两轧辊辊径、辊速和表面摩擦系数均相同，转动方向相反；②轧制过程中板件仅受轧制力作用，不存在其它外力影响；③仅对轧件划分较细的网格，轧辊为解析刚性体薄壁圆柱筒，不需要划分网格，以提高运算效率。

据此可建立基于弹塑性变形的宏观同步热轧三维模型如图4所示，其中轧辊直径为270 mm，轧制温度为400 ℃，具体轧制方案见表2，材料受温度影响的物性参数由CALPHAD计算相图方法得到，见表3所示，而塑性参数由铸坯试样的拉伸实验得到（由于数据量过大在此不一一列出）。

有限元模拟后可输出轧件等效塑性应变场和应力场，以及代表单元的位移差-时间曲线等计算结果。其中前者可经Origin软件处理生成应力-应变曲线并将其与各道次拉伸曲线进行对比拟合；后者可作为边界条件代入到多晶体塑性模型中，运用CPFEM模拟ZK61镁合金沿RD单向拉伸过程，将模拟得到的应力-应变曲线和晶体织构分别与拉伸曲线与EBSD织构图进行对比分析。

2.2 微观CPFEM模型

2.2.1 晶体塑性本构方程的建立

假设单晶体的非弹性变形完全是由晶体滑移引起，并可将晶体孪生认为是一种“伪滑移”^［18-19］，从而以处理位错剪切滑移的方式来统一解决晶体滑移和孪生的变形机制，而在高温塑性变形下，镁合金晶界处会存在黏滞性^［20］等不利于建立热轧态晶体塑性本构模型的复杂情况，故将其放在后续的研究方向之中。由于晶体在塑性变形时相继经历弹性变形和旋转，以及塑性剪切作用^［21］，因此总变形梯度F由式（1）表示：

F=F^{\mathrm{*}}\cdot F^{\mathrm{p}}（1）

式中：F^{\mathrm{*}}为晶格的拉伸变形和旋转变形梯度，其晶格方向和间距与原始参考位形相同；F^{\mathrm{p}}为材料中间参考位形的塑性剪切变形梯度。则此时的速度梯度L如式（2）：

L\equiv \dot{F}\cdot F^{-\mathrm{1}}=D+\Omega（2）

式中：D为拉伸对称速率；\Omega为反对称自旋张量；L是一种具有完整对称集的张量。

在研究晶体塑性变形过程中，需要假设晶体滑移符合施密特定律，即忽略弹性晶格畸变，剪应力仅由施密特应力决定。而弹性势能的存在保证了L与晶格的对称拉伸速率D^{\mathrm{*}} 的关系如式（3）：

\stackrel{{\nabla }^{\mathrm{*}}}{\sigma }+\sigma \left(I:D^{\mathrm{*}}\right)=L:D（3）

式中：\sigma为柯西应力；\stackrel{{\nabla }^{\mathrm{*}}}{\sigma }为焦曼应力率，与晶格旋转轴的关联应力率有关。而在有限变形的应力应变尺度上，可以较好地保持小变形理论架构，因此可定义分解切应力{\tau }^{\left(\mathrm{\alpha }\right)}如式（4）：

{\tau }^{\left(\mathrm{\alpha }\right)}=m^{\mathrm{*}\left(\mathrm{\alpha }\right)}\cdot \frac{{\rho }_{\mathrm{0}}}{\rho }\sigma \cdot s^{\mathrm{*}\left(\mathrm{\alpha }\right)}（4）

式中：s^{\mathrm{*}\left(\mathrm{\alpha }\right)}为滑移系\alpha的滑移向量；m^{\mathrm{*}\left(\mathrm{\alpha }\right)}为相应滑移面法向量；{\rho }_{\mathrm{0}}和\rho分别为参考状态和瞬时状态下的密度。

根据施密特定律，速率相关的晶体中，滑移系\mathrm{\alpha }的滑移速率{\dot{\gamma }}^{\left(\mathrm{\alpha }\right)}如式（5）：

{\dot{\gamma }}^{\left(\mathrm{\alpha }\right)}={\dot{a}}^{\left(\mathrm{\alpha }\right)}{f}^{\left(\mathrm{\alpha }\right)}\cdot \left({\tau }^{\left(\mathrm{\alpha }\right)}/g^{\left(\mathrm{\alpha }\right)}\right)（5）

式中：{\dot{a}}^{\left(\mathrm{\alpha }\right)}为滑移系\mathrm{\alpha }上的参考应变率；f^{\left(\mathrm{\alpha }\right)}\left(x\right)=x{\left|x\right|}^{n-\mathrm{1}}是一个描述应变率与应力关系的初等函数，其中n为速率敏感度指数；g^{\left(\mathrm{\alpha }\right)}为滑移系\mathrm{\alpha }的当前应变硬化强度，而其增量关系的演化能够体现出应变硬化的特征如式（6）：

{\dot{g}}^{\left(\mathrm{\alpha }\right)}=\sum _{\mathrm{\beta }}{h}_{\mathrm{\alpha }\mathrm{\beta }}{\dot{\gamma }}^{\left(\mathrm{\beta }\right)}（6）

式中：h_{\mathrm{\alpha }\mathrm{\beta }}为滑移硬化模量，即一切已激活滑移系\mathrm{\alpha }\mathrm{、}\mathrm{\beta }的求和区间，其中h_{\mathrm{\alpha }\mathrm{\alpha }}（不对\mathrm{\alpha }求和）为自硬化模量，而h_{\mathrm{\alpha }\mathrm{\beta }}（\mathrm{\alpha }\ne \mathrm{\beta }）为潜在硬化模量。而它们的表达式取决于一切滑移系的剪切应变{\gamma }^{\left(\mathrm{\alpha }\right)}如式（7）：

\left\{\begin{array}{l}{h}_{\mathrm{\alpha }\mathrm{\alpha }}=\left\{\left(h_{\mathrm{0}}-h_{\mathrm{s}}\right)\mathrm{s}\mathrm{e}\mathrm{c}{\mathrm{h}}^{\mathrm{2}}\left[\frac{\left(h_{\mathrm{0}}-h_s\right){\gamma }^{\left(\mathrm{\alpha }\right)}}{{\tau }_{\mathrm{s}}-{\tau }_{\mathrm{0}}}\right]+h_{\mathrm{s}}\right\}G\left[{\gamma }^{\left(\mathrm{\beta }\right)};\mathrm{\beta }\ne \mathrm{\alpha }\right]\stackrel{}{}(\mathrm{不对}\mathrm{\alpha }\mathrm{求和})\\ h_{\mathrm{\alpha }\mathrm{\beta }}=q{h}_{\mathrm{\alpha }\mathrm{\alpha }}\stackrel{}{}\left(\alpha \ne \beta \right)\end{array}\right.（7）

式中：h_{\mathrm{0}}，h_{\mathrm{s}}分别为材料的初始硬化模量和第一阶段的硬化模量；{\tau }_{\mathrm{0}}，{\tau }_{\mathrm{s}}分别为与当前强度值g^{\left(\mathrm{\alpha }\right)}\left(\mathrm{0}\right)等价的初始屈服应力值和第I阶段应力（或者是剧烈塑性流动时产生的穿透应力）值；q为两种滑移硬化模量之间的关系常数；而交叉硬化函数G\left[{\gamma }^{\left(\mathrm{\beta }\right)};\mathrm{\beta }\ne \mathrm{\alpha }\right]=\mathrm{1}+\sum _{\mathrm{\beta }\ne \mathrm{\alpha }}{f}_{\mathrm{\alpha }\mathrm{\beta }}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left({\gamma }^{\left(\mathrm{\beta }\right)}/{\gamma }_{\mathrm{0}}\right)，其中{\gamma }_{\mathrm{0}}为滑移系间相互作用达到峰值强度后的滑移量，以及每个组成部分f_{\mathrm{\alpha }\mathrm{\beta }}表示一种特定滑移的相互作用的强度大小。

而孪生的剪切变形机制与滑移剪切不同，由于高能势垒的存在，孪生过程的剪切是不对称的^［22］，可据此采用与滑移相似的速率相关晶体塑性流动方程来表示孪生变形时的剪切应变如式（8）：

（8）

式中：{\dot{\gamma }}^{\mathrm{0}}为参考剪切应变率；{\tau }_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\left(\alpha \right)}为背应力的等效分切应力。这里考虑到孪生的极性，在式（8）中定义孪生变形只有在其作用位置的等效剪应力为正时才开始进入计算。

以上是本次研究所建立的晶体塑性本构方程及相应描述了三个强化阶段的率相关晶体材料的应变硬化模型公式推导过程，可以将其与有限元分析方法相结合，进行CPFEM模拟运算。

2.2.2 CPFEM模拟的实现

相应晶体塑性本构方程调用及CPFEM模拟流程如图5所示。在晶体塑性模拟过程中，需要将材料参数输入到材料常数矩阵，主要包括HCP晶体弹性模量、潜在激活的滑移系参数、参考状态下全局坐标系中晶体的初始取向、晶体塑性本构方程相关参数、自硬化与潜在硬化模量、前向梯度时间积分参数以及迭代法参数。将矩阵代入到子例程UMAT执行函数运算，使应力和解相关状态变量更新到增量结束时的值，并为本构模型提供了材料的雅可比矩阵；而产生的应力、应变和状态变量由ABAQUS软件进行基于隐式积分方法的增量求解。当调用子例程时，将提供增量开始时的状态（应力、解相关状态变量）以及（估计的）应变增量和时间增量。其中所需晶体参数详见表4，激活的滑移系与孪晶系详见表5，滑移系硬化参数详见表6。

对试样大变形区取一个代表单元，以初始晶粒尺寸为66.67 μm、变形温度为400 ℃、变形量为25%的实验为例，输出代表单元的轴向位移与时间关系，并设置为多晶体模型位移边界条件。在试样大变形区域选取具有代表性的单元，选取宏观模拟中代表单元的位移差-时间曲线作为多晶体模型的边界条件，利用代表单元的位移边界条件作为微观尺度下多晶体变形模型的边界条件。多晶体模型的建立采用基于Voronoi图的三维晶粒划分来重构晶粒形貌，以近似反映真实的微观晶粒分布，其中X轴对应RD，Y轴对应TD，Z轴对应ND。分别输出第三、五道次的轧制多晶模型，其中第三道次（图6（a），（b））晶粒尺寸为50×50×50，网格划分为10×10×10，晶粒数为27个；第五道次（图6（c），（d））晶粒尺寸为30×30×30，网格划分为20×20×20，晶粒数为30个。两道次多晶网格划分结果如图6所示。

3 结果与分析

3.1 实验结果分析

3.1.1 拉伸实验结果分析

对第三到第六道次轧制板材取样后进行拉伸实验，数据分三个方向进行处理比较：0°（RD）方向，45°方向，90°（TD）方向。各道次试样真应力-真应变曲线如图7所示，相应拉伸性能变化如图8所示。

由图7、8可看出，ZK61镁合金热轧板在不同道次、不同方向上的弹性和塑性能力均有较大差异。在不同道次上，合金各方向的弹性模量均逐渐减小（图8（a）），0°与45°方向的屈服强度（图8（b））和抗拉强度（图8（c））也不断下降，而90°方向的塑性指标波动较大。第五道次的应变率（图7（c））和伸长率（图8（d））均处于高位，说明轧件在第五道次塑性最强；而在不同方向上，合金在第三、五道次的弹性模量（图8（a））与屈服强度（图8（b））由0°到90°逐渐下降，而抗拉强度呈上升趋势（图8（c）），应变率和伸长率的最大值分别出现在第三、五道次的45°方向，分别为0.232（图7（a））和22%（图8（d）），而第六道次板材强度提升了7.55%，伸长率也达到了19.5%（图8（d））。实验结果体现出热轧板材在不同道次及不同方向的弹性和塑性均有较大差异，说明热轧变形后的ZK61镁合金板材表现出了较强的各向异性，即热轧道次会对ZK61的微观结构演变产生极大影响。

3.1.2 EBSD数据处理与织构分析

镁合金板材在不同轧制道次下的微观组织结构会发生一定改变，进而影响其力学性能。分别对初始铸态板材与经各道次热轧后的板材试样进行EBSD检测分析，得出ZK61镁合金板材不同轧制道次的微观组织IPF图如图9所示。通过观察发现，与初始板坯（图9（a））对比，经过不同道次轧制后的板材晶粒明显细化，前四道次（图9（b）~（e））轧制形成了许多孪晶，而后两道次（图9（f），（g））基本都是细小均匀的晶粒，还可以看到很多晶粒联结在一起，形成了项链状组织形貌。对其中各状态的晶粒尺寸进行统计，初始态晶粒平均直径为66.67 μm，轧制第一道次晶粒平均直径便降低至16.81 μm，第二道次为15.32 μm，第三道次略升至16.42 μm，第四道次又降到11.99 μm，第五道次继续降至7.49 μm，而第六道次又回升为9.13 μm，表明第一、四、五道次对板材晶粒起到非常明显的破碎和细化作用，而在第三道次和第六道次则是温度对变形晶粒的回复和长大机制发挥了一定影响。这是由于镁合金在热轧时内部能量较高，在经过一定道次的轧制和道次间退火后晶粒产生回复再结晶，使微观组织发生局部变形，晶粒的畸变消失，能量释放，出现新的等轴晶组织。

进一步分析发现，初始状态下\left(\mathrm{0001}\right)、\left(\mathrm{01}\overline{\mathrm{1}}\mathrm{0}\right)和\left(\overline{\mathrm{1}}\mathrm{2}\overline{\mathrm{1}}\mathrm{0}\right)织构占比都很均匀，随着轧制道次的增加，\left(\mathrm{0001}\right)和\left(\overline{\mathrm{1}}\mathrm{2}\overline{\mathrm{1}}\mathrm{0}\right)织构占比慢慢的开始减少，\left(\mathrm{0001}\right)织构占比增加，第六道次\left(\mathrm{01}\overline{\mathrm{1}}\mathrm{0}\right)和\left(\overline{\mathrm{1}}\mathrm{2}\overline{\mathrm{1}}\mathrm{0}\right)织构占比很小，基本只有\left(\mathrm{0001}\right)织构。这是由于同步轧制过程使部分晶粒取向发生改变，而在这些晶粒发生再结晶时，周围刚开始形核的小晶粒取向没有任何规律，使晶体内部c轴取向偏离板材ND的晶粒逐渐变多，对基面织构产生一定影响。

观察轧制不同道次试样的取向差分布结果如图10所示。初始态（图10（a））和第一道次（图10（b））状态下样品在88°的大角度晶界呈现出峰值取向差角，第二、三、四道次（图10（c）~（e））状态下样品在8°的小角度晶界呈现出峰值取向差角，第五和第六道次（图10（f），（g））状态下样品在28°的大角度晶界呈现出峰值取向差角。而出现该情况的主要原因由位错滑移与孪生两种机制相互影响与共同作用导致的^［23］：轧制过程初期激活了较多取向差角为86.3°的\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}拉伸孪晶，而随着轧制的进行材料内部积累了大量位错，使位错滑移成为了主要的轧制变形机制，在位错滑移与缠结的过程中不断产生细小再结晶晶粒，最后两道次则是由于变形晶粒内部普遍激活了取向差角为22°的\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}-\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{3}\right\}二次孪晶主导了共同作用。

典型孪晶类型以及其和基体的取向角度差见表7，根据表中还可以观察到，随着轧制道次的增加，拉伸孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}、压缩孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{1}\right\}和\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{3}\right\}的总体趋势也逐渐减小，二次孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{1}\right\}-\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}的总体趋势基本不变，而二次孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}-\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{3}\right\}的总体趋势逐渐增加。

3.2 模拟结果分析

3.2.1 宏观轧制模拟结果分析

为保证宏观轧制模拟的准确性和代表性，需要选取板材的RD-TD面，并分别从等效应力场和等效塑性应变场两个方面的变化结果进行分析，其中轧制板材的等效应力场变化结果如图11所示。

其中第一道次轧制的最大等效应力值为204.0 MPa（图11（a）），第二道次轧制的最大等效应力值为70.3 MPa（图11（b）），第三道次轧制的最大等效应力值为218.7 MPa（图11（c）），第四道次轧制的最大等效应力值为218.8 MPa（图11（d）），第五道次轧制的最大等效应力值为219.2 MPa（图11（e）），第六道次轧制的最大等效应力值为219.4 MPa（图11（f））。由图可以看出，从第一道次到第二道次轧制，最大等效应力值迅速降低，之后一直到第六道次轧制的最大等效应力值逐步有极小幅增加。因为第一道次轧制变形量很小，且一开始有温度的影响；而第二道次变形量与第一道次变形量差别很大，导致其最大等效应力值降低。从第二道次到第六道次的变形量差别不大，是由于热轧板材在受到变形抗力和温度效应的共同作用下趋近于稳定状态^［24-25］，较符合轧制使其最大等效应力值增加的规律，即等效应力在变形量作用的情况下发生了微小变化。

观察图12中各道次轧制等效塑性应变云图可知，第一道次轧制的最大等效塑性应变值为0.2335（图12（a）），第二道次轧制的最大等效塑性应变值为0.98（图12（b）），第三道次轧制的最大等效塑性应变值为0.8883（图12（c）），第四道次轧制的最大等效塑性应变值为0.633（图12（d）），第五道次轧制的最大等效塑性应变值为0.632（图12（e）），第六道次轧制的最大等效塑性应变值为0.498（图12（f））。由图可以看出，从第一道次到第二道次轧制过程中，最大等效塑性应变值迅速增加，之后一直到第六道次轧制过程中的最大等效塑性应变值逐步减小；沿厚度方向观察，上下两表面等效塑性应变值比心部的大，因为第一道次轧制变形量很小，且一开始有温度的影响，而第二道次变形量与第一道次变形量差别很大，进而使其最大等效塑性应变值增加。后续各轧制道次变形量差别不大，符合轧制变形使板材最大等效塑性应变减小的规律。

3.2.2 微观晶粒状态变化分析

基于拉伸试验结果，选取力学性能最突出的第三与第五道次相关数据引入多晶体塑性模型中，进行CPFEM数值计算以模拟ZK61镁合金同步热轧时微观晶粒状态变化情况。

图13为RD-TD面第三道次的多晶模型晶粒演化模拟结果，其中晶粒尺寸为50×50×50，网格划分为10×10×10，晶粒个数为27个。初始时晶粒1区域的应力较大，晶粒2、3、5、6、9区域的应力居中（图13（a））；轧制过程中晶粒1、6、10区域的应力居中，晶粒2、3、4、9区域的应力较小（图13（b））；道次结束后晶粒2、4、7、9区域的应力较小，晶粒1、6区域的应力居中（图13（c））。由图可知该道次轧制过程板材受力整体分布尚不均匀，塑性加工强化作用不足。

图14为RD-TD面第五道次的多晶模型晶粒演化模拟结果，其中晶粒尺寸为30×30×30，网格划分为20×20×20，晶粒个数为30个。初始时晶粒6区域的应力较大，晶粒1、3、7、8、9区域的应力居中，晶粒2、3、4、5区域的应力较小（图14（a））；轧制过程中晶粒1、5区域的应力较大，晶粒2、4、6、7、8区域的应力居中（图14（b））；道次结束后晶粒1、3、4区域的应力居中，晶粒2、6、7、8区域的应力较小（图14（c））。而整体上随着轧制过程的进行，晶粒所受应力值均获得较大增长，但涨幅与比例已不及第三道次，说明随着轧制的进行，板材受到一定的加工硬化与回复再结晶的共同作用，使微观组织的变化逐渐趋于稳定。

3.2.3 织构演化过程分析

ZK61镁合金的初始实验极图和第三、五、六道次的模拟极图与实验极图如图15所示。ZK61的初始状态极密度为15.84（图15（a）），第三道次为8.57（图15（c）），第五道次为15.28（图15（e）），第六道次为17.87（图15（g））。将CPFEM的模拟织构与实验织构对比分析，可以观察到随着轧制道次的增加，基面织构极密度的总体趋势是增加的。基面织构\left\{\mathrm{0001}\right\}随着轧制道次的增加越来越强烈，最后两道次\left\{\mathrm{0001}\right\}织构在中心时最大，而其他两个织构\left\{\mathrm{11}\overline{\mathrm{2}}\mathrm{0}\right\}和\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{0}\right\}一直变化不大，因此基面织构是镁合金最主要的织构，由轧制道次来看，最大织构强度的总体趋势随着轧制道次的增加而增加。通过实验与模拟结果对比，可观察到随着镁合金轧制道次的增加，大量基面织构\left\{\mathrm{0001}\right\}的晶粒靠近c轴中心，因此各向异性较强，基面织构较强烈。

4 结论

（1）对ZK61镁合金板材进行轧制实验，实验结果表明：经轧制后的镁合金板材在不同轧制道次条件下，其微观组织结构也会发生改变，进而影响其力学性能。经过不同道次轧制后的板材形成了许多孪晶，晶粒明显细化，分布较为均匀。

（2）由四道次轧制板材三个方向的室温单向拉伸实验可得出，不同方向不同道次的力学性能各不相同，这说明在室温下，ZK61镁合金板材表现出较强的塑性各向异性，而其力学性能及塑性变形能力均随各道次轧制有着不同的变化趋势，说明热轧道次对ZK61的微观结构演变产生了较大影响，经六道次热轧后板材的强度提升了7.55%，伸长率达到19.5%。

（3）观察轧制不同道次试样的取向差分布结果，可以观察到，初始状态和第一道次状态下样品在大角度晶界范围（88°）时呈现出峰值取向差角，第二道次状态、第三道次状态和第四道次状态下样品在小角度晶界范围（8°）时呈现出峰值取向差角，第五道次状态和第六道次状态下样品在大角度晶界范围（28°）时呈现出峰值取向差角，主要因为轧制试样中存在很多再结晶晶粒。根据取向差分布图还可以观察到，随着轧制道次的增加，拉伸孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}的总体趋势逐渐减小，压缩孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{1}\right\}和\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{3}\right\}的总体趋势也逐渐减小，二次孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{1}\right\}-\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}的总体趋势基本不变，而二次孪晶\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{2}\right\}-\left\{\mathrm{10}\overline{\mathrm{1}}\mathrm{3}\right\}的总体趋势逐渐增加。

（4）在CPFEM模拟中观察到轧制板材晶粒应力的变化，可得出网格划分为10×10×10的第三道次板材内侧和外侧的晶粒应力都减小；网格划分为20×20×20的第五道次板材的内侧晶粒应力不变，外侧晶粒应力减小。将CPFEM的模拟织构与实验织构对比分析，可观察到随着轧制道次的增加，基面织构极密度的总体趋势是增加的。

参考文献

[1] 李骆宾，江鹏，任广猛，轧制路径对ZK60镁合金板材组织和性能的影响［J］.金属热处理，2020，45（2）：36-41.

[2] LI L B， JIANG P， REN G M，et al.Effect of rolling route on microstructure and properties of ZK60 magnesium alloy sheets［J］.Heat Treatment of Metals，2020，45（2）：36-41.

[3] 蓝永庭，陈渊，任一方，基于滑移与孪生镁合金晶体塑性本构模型及微结构关联分析［J］.中国有色金属学报，2019，29（8）：1660-1674.

[4] LAN Y T， CHEN Y， REN Y F，et al.Constitutive model based on slip and twinning of AZ31 Mg alloy and analysis of microstructural relatedness［J］.The Chinese Journal of Nonferrous Metals，2019，29（8）：1660-1674.

[5] 陈渊.含孪生机制的镁合金本构关系及微结构关联性研究［D］.南宁：广西大学，2018.

[6] CHEN Y.Study on constitutive modeling including the mechanism of twinning and microstructure relevance in magnesium alloy［D］.Nanning：Guangxi University，2018.

[7] CHEN S F， SONG H W， ZHANG S H，et al.Effect of shear deformation on plasticity，recrystallization mechanism and texture evolution of Mg-3Al-1Zn alloy sheet：experiment and coupled finite element-VPSC simulation［J］.Journal of Alloys and Compounds，2019，805：138-152.

[8] 李文雪.利用晶体塑性有限元法研究晶粒尺度下的Mg-Y合金变形机制［D］.上海：上海交通大学，2019.

[9] LI W X.The study on grain-scale deformation mechanism of Mg-Y alloys by crystal plastic finite element method［D］.Shanghai：Shanghai Jiao Tong University，2019.

[10] SCHMIDT C， KAWALLA R， WALDE T，et al.Experimental and numerical investigation of texture development during hot rolling of magnesium alloy AZ31［J］.Materials Science Forum，2007，539-543，3448-3453.

[11] PEIRCE D， ASARO R J， NEEDLEMAN A.Material rate dependence and localized deformation in crystalline solids［J］.Acta Metallurgica，1983，31（12）：1951-1976.

[12] MATHUR K K， DAWSON P R， KOCKS U F.On modeling anisotropy in deformation processes involving textured polycrystals with distorted grain shape［J］.Mechanics of Materials，1990，10（3）：183-202.

[13] KALIDINDI S R.Incorporation of deformation twinning in crystal plasticity models［J］.Journal of the Mechanics and Physics of Solids，1998，46（2）：267-290.

[14] 刘相华.塑性有限元在金属轧制过程中应用的进展［J］.金属学报，2010，46（9）：1025-1033.

[15] LIU X H.Progress and application of plastic finite element method in metals rolling process［J］.Acta Metallurgica Sinica，2010，46（9）：1025-1033.

[16] LI Y Q， ZHANG H M， SHANG X Q，et al.A multiscale investigation on the preferential deformation mechanism of coarse grains in the mixed-grain structure of 316LN steel［J］.International Journal of Plasticity，2022，152：103244.

[17] FU J Q， LIU W B， SUI H N，et al.Simulations of the localized necking and Luders band in irradiated metals by crystal plasticity［J］.Acta Mechanica Sinica，2023，39：122339.

[18] ZHANG H M， LIU J， SUI D S，et al.Study of microstructural grain and geometric size effects on plastic heterogeneities at grainlevel by using crystal plasticity modeling with high-fidelity representative microstructures［J］.International Journal of Plasticity，2018，100：69-89.

[19] SIGNORELLI J W， BERTINETTI M A， ROATTA A.A review of recent investigations using the Marciniak-Kuczynski technique in conjunction with crystal plasticity models［J］.Journal of Materials Processing Technology，2021，287：116517.

[20] TANG X F， PENG L F， WANG L Y，et al.Influence of crystal structure on size dependent deformation behavior and strain heterogeneity in micro-scale deformation［J］.International Journal of Plasticity，2019，118：147-172.

[21] SHANG X Q， ZHANG H M， SHANG X Q，et al.The effect of stress state and strain partition mode on the damage behavior of a Mg-Ca alloy［J］.International Journal of Plasticity，2021，144：103040.

[22] 王浩雷.ZK61镁合金中厚板轧制过程的数值模拟及晶体塑性变形行为研究［D］.鞍山：辽宁科技大学，2022.

[23] WANG H L.Numerical simulation of ZK61 magnesium alloy rolling process and study of crystal plastic deformation behavior［D］.Anshan：University of Science and Technology Liaoning，2022.

[24] 徐帅.AZ31镁合金薄板尺寸效应及其全场晶体塑性模拟研究［D］.上海：上海交通大学，2019.

[25] XU S.Size effects of AZ31 Mg foils and the full-field crystal plasticity modeling［D］.Shanghai：Shanghai Jiao Tong University，2019.

[26] 宋畅.多晶铜中形变孪晶的晶体塑性有限元模拟［D］.大连：大连理工大学，2022.

[27] SONG C.Crystal plasticity finite element simulation of deformation twins in polycrystalline copper［D］.Dalian：Dalian Universityof Technology，2022.

[28] 高建中.Zn添加调控Mg-9Gd-2Nd-0.5Zr合金析出相演变及对合金高温性能的影响［D］.太原：太原理工大学，2022.

[29] GAO J Z.The effect of Zn addition on the evolution of precipitated phase and high temperature properties of Mg-9Gd-2Nd-0.5Z ralloy was investigated［D］.Taiyuan：Taiyuan University of Technology，2022.

[30] PAUDEL Y R， INDECK J， HAZELI K，et al.Characterization and modeling of 10 1 ¯ 2 twin banding in magnesium ［J］.Acta Materialia，2020，183：438-451.

[31] MUHAMMAD W，ALI U， BRAHME A P，et al.Experimental analyses and numerical modeling of texture evolution and the development of surface roughness during bending of an extruded aluminum alloy using a multiscale modeling framework［J］.International Journal of Plasticity，2019，117：93-121.

[32] 陈扬.高应变速率下 AZ31 镁合金的各向异性研究［D］.沈阳：沈阳工业大学，2022.

[33] CHEN Y.Research on anisotropy of AZ31 magnesium alloy under high strain rates［D］.Shenyang：Shenyang University of Technology，2022.

[34] 郑玄玄.AZ31镁合金板材轧制过程有限元模拟与实验研究［D］.镇江：江苏大学，2018.

[35] ZHENG X X.Finite element simulation and experiment research on rolling process of AZ31 magnesium alloy plate［D］.Zhenjiang：Jiangsu University，2018.

[36] 周永强.304/Q345R不锈钢复合板轧制数值模拟及轧后组织性能分析［D］.秦皇岛：燕山大学，2016.

[37] ZHOU Y Q.304/Q 345R stainless steel composite plate rolling numerical microstructure performance analysis plate［D］.Qinhuangdao：Yanshan University，2016.

基金资助

辽宁省教育厅重点基金(LJKZ0281)