, 刘相华1, 孙祥坤11. 东北大学 轧制技术及连轧自动化国家重点实验室， 辽宁 沈阳 110819;
2. 东北大学 材料科学与工程学院， 辽宁 沈阳 110819
收稿日期: 2015-03-27
基金项目: 国家自然科学基金资助项目(51374069, U1460107).
作者简介: 陈守东(1987-),男,安徽蚌埠人,东北大学博士研究生; 刘相华(1953-),男,黑龙江双鸭山人,东北大学教授,博士生导师.
摘要: 为了定量描述厚度尺寸对极薄带轧制微观变形不均匀性的影响，采用率相关晶体塑性理论和Voronoi图的多晶模型，考虑试样尺寸，晶体取向及其分布，模拟了不同厚度铜箔在相同压下率条件下的变形行为，得到了在细观尺度上铜箔的微观应力应变和滑移系分布.模拟获得的应力-应变曲线和实验测得的曲线基本一致.通过对20%压下率不同厚度铜箔的轧制变形的研究表明，无论是在晶粒内部还是在晶粒间，材料变形都非常不均匀，主要是由在铜箔厚度方向上只有一层晶粒时，晶粒尺寸、形貌和取向的不均匀分布，近邻晶粒取向差以及滑移系启动特性所引起的.
关键词: 极薄带轧制晶体塑性有限元非均匀变形滑移系启动
Crystal Plasticity Finite Element Simulation of Slip System Evolution in Pure Copper Foil Rolling
CHEN Shou-dong1, LIU Xiang-hua1
, LIU Li-zhong2, SUN Xiang-kun11. State Key Laboratory of Rolling and Automation, Northeastern University, Shenyang 110819, China;
2. School of Materials Science & Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: LIU Xiang-hua, professor, E-mail: liuxh@mail.neu.edu.cn
Abstract: The distributions of the stress, strain and slip systems in pure copper foil rolling with the same reduction were simulated by the rate-dependent crystal plasticity theory and Voronoi polycrystalline model with respect to specimen dimension, grain orientation and its distribution to quantitatively evaluate the influence of thickness on inhomogeneous deformation behavior of foil rolling at mesoscale. The simulated stress-strain curves agree well with the experimental results. The simulation results reveal that the deformation behavior in the polycrystalline aggregate is inhomogeneous not only in intracrystalline but also in intergranule with a 20% reduction in foil rolling. The foils are composed of only a single layer grain across thickness, the grains with different sizes, shapes and orientations are unevenly distributed in the foil, the misorientation of neighboring grains and the property of active slip systems.
Key words: foil rollingcrystal plasticity finite elementinhomogeneous deformationslip system activity 近年来微制造、微电子、计算机等高技术的快速发展，提出了一系列微材料和微成形问题，对金属箔材的需求量增加.当材料的本征尺寸减小到微尺度时，材料的各向异性增加，材料的局部性能有可能超过材料宏观结构的平均效应[1].当箔材尺寸降低到微尺度时，由于尺寸效应的出现，宏观连续的变形理论不能直接应用到微成形中[2].金属箔材厚度尺寸与其晶粒尺寸具有相同的数量级，就意味着在箔材厚度方向上只分布一个或几个少数晶粒.
晶体塑性有限元方法(CPFEM)对离散分布的位错进行均匀化处理，采用连续介质的场变量变形梯度张量来描述位错滑移的宏观效应，能够反映材料的微观结构特性[3]，被广泛地用来研究单晶体、多晶体的变形行为和织构演化等[4, 5, 6].Parasiz等[7]发现当晶粒尺寸减小到与试样尺寸同一数量级时，材料的变形行为主要受到单个晶粒性能的影响，导致非均匀应力和硬度的分布.Kim等[8]和Yun等[9]通过修正经典的霍尔-佩奇公式来描述薄带拉伸过程中应力-应变变化.尽管经过适当的修正得到的应力-应变关系与实验结果一致，但他们都没有考虑晶粒取向的影响.Lu等[10]建立了一种新的微观变形模型，发现保持试样尺寸不变的情况下，材料变形的分散性随着晶粒尺寸的增大而急剧增加.
变形区内不同晶粒分布、尺寸、取向等对材料变形具有重要影响，目前还没有考虑这些因素针对极薄带轧制建立相应的微观模型.本文采用晶体塑性有限元和Voronoi图的多晶模型，考虑试样尺寸、晶体取向分布，分析极薄带轧制变形中铜箔厚度对微观应力应变场和滑移系分布的影响.
1 理论和计算模型1.1 晶体塑性理论采用Peirce[11]提出的率相关晶体塑性本构模型.整个晶体变形分解为首先发生位错沿着滑移系滑移的塑性变形(Fp)，然后进行晶格畸变和旋转的弹性变形(Fe).总的变形梯度F分解为
率相关晶体塑性模型采用指数方程描述第α滑移系上滑移剪切率与分剪切应力之间的关系：
式中：
为参考剪切应变率；n为率敏感系数；τcα为临界分剪切应力；τα为滑移系上的分剪切应力.
采用Bassani和Wu[12]提出的适合FCC金属的硬化模型，其表达式为
式中:
为β滑移系的剪切率；hαα为自硬化系数；hαβ为潜硬化系数；γ0为参考剪切应变；γ为累积剪切应变；h0为初始硬化率；hs为易滑移阶段硬化模量；τ0为初始临界分剪切应力；τ1为临界分剪切应力饱和值；fαβ表示滑移系α与滑移系β间的相互作用系数，由5个常数a1~a5决定；N为滑移系个数.1.2 计算模型纯铜的弹性模量为C11= 168.4 GPa,C12= 121.4 GPa,C44=75.4 GPa.因纯铜为FCC晶体结构，位错滑移发生在12个滑移系{1 1 1} < 1 1 0>上，如表 1所示.模拟过程未知的材料硬化参数，通过模拟计算和实验相结合的方法确定.对多晶体铜箔拉伸变形进行模拟，将得到的应力-应变曲线与实验测定的曲线进行对比，选取合适的参数，如表 2所示.对厚度为125 μm的铜箔进行拉伸实验和模拟，拉伸速度为0.05 mm/min，模拟和实验得到的应力-应变曲线如图 1所示.由图 1可知模拟结果和实验结果相吻合.
图1(Fig. 1)
 | 图1 模拟和实验测得的应力-应变曲线的比较 Fig. 1 Comparison of the experimental and CPFEM simulated uniaxial tensile stress-strain curve |
表1(Table 1)
 表1 面心立方金属的滑移系Table 1 The slip system of face centered cubic metal
| 表1 面心立方金属的滑移系 Table 1 The slip system of face centered cubic metal |
表2(Table 2)
表2 计算中采用的本构模型参数Table 2 Parameters in the constitutive model considered by the simulation
| 表2 计算中采用的本构模型参数 Table 2 Parameters in the constitutive model considered by the simulation |
/ s-1采用基于Voronoi图的多晶模型重构晶粒的形貌，来近似反映真实的材料内部微结构.构造的多晶模型的平均晶粒尺寸d=100 μm.将极薄带轧制变形简化为二维平面应变变形，其中X轴对应轧制方向 (RD)，Y轴对应法向 (ND)，Z轴对应横向 (TD).上轧辊角速度设置为1.04 rad/s，异速比设置为1.1.对每一厚度的铜箔进行一道次20%压下率的轧制变形，轧辊与轧件之间采用库伦摩擦，摩擦系数μ=0.1.影响极薄带轧制微观不均匀变形的另一个重要因素是晶体取向分布.
退火后铜箔沿法向的极图如图 2所示.由图 可知，经过完全退火后，铜箔几乎没有织构，因此 对建立的多晶模型赋予随机晶体取向.建立多晶 模型时，保持多晶模型长度L=5 mm和晶粒尺寸d=100 μm不变，只改变铜箔厚度尺寸.当铜箔厚度(t)/晶粒尺寸(d)=1时，建立的多晶模型中晶粒形状为竹节晶；t/d<1时，建立的多晶模型中晶粒形状为饼状晶.模拟过程中，轧辊设定为刚体，铜箔为变形体，采用CPE4R单元离散多晶模型，每个晶粒包含100个单元.长度方向上有50个晶粒，厚度方向上有1个晶粒的不同厚度铜箔多晶模型如图 3所示，其中不同灰度代表不同晶粒取向.
图2(Fig. 2)
 | 图2 退火铜箔法向极图 Fig. 2 Normal direction of the grain orientations on the experimental foil |
图3(Fig. 3)
 | 图3 不同厚度铜箔初始计算模型 Fig. 3 Initial mesh and orientation of simulation for different thickness of copper foils(a)—t=20 μm； (b)—t=50 μm; (c)—t=100 μm. |
2 模拟结果与讨论图 4~图 5为经过20%压下率轧制变形后3种厚度铜箔切应力和真应变计算结果.由图可知，铜箔的微观应力分布极不均匀，存在应力集中区，具有明显的应力梯度特点，不仅在整个厚度方向上分布不均匀，而且在晶粒内部应力分布也不均匀，尤其在晶界处和100 μm厚铜箔的上下表面处.晶界不一定是应力最大处，晶粒内部可能出现应力集中处，如图 4所示，H和N的晶粒内部和D,G,L晶界处都出现了应力集中.随着铜箔厚度的增加，应力集中区逐渐向箔材中心和表面转移.图 5所示的真应变也表现出同样的特性.在当前计算条件下，在晶粒内部和晶粒之间出现了剪切变形带，随着铜箔厚度的减小，变形带基本贯穿了整个晶粒.这一方面由于晶粒具有不同的取向，滑移系启动先后不同造成晶粒的变形顺序不同，晶粒的变形量也不同；另一方面由于多晶体变形的协调性，晶粒的变形受到相邻晶粒变形的影响，从而形成应变集中区和变形带.当铜箔厚度较小时，晶粒尺寸远远超过厚度尺寸，在晶粒内部形成了很多方向不同的小变形带，尤其集中在铜箔的上下表面处.相邻晶粒具有不同的应力应变和集中区分布，可见晶粒取向分布在很大程度上影响极薄带轧制微观变形行为.
图4(Fig. 4)
 | 图4 3种厚度铜箔变形后切应力分布 Fig. 4 Contour plots of shear stress in roll bite with three thickness(a)—t=20 μm； (b)—t=50 μm; (c)—t=100 μm. |
图5(Fig. 5)
 | 图5 3种厚度铜箔变形后对数应变分布 Fig. 5 Contour plots of logarithmic strain in roll bite with three thickness(a)—t=20 μm； (b)—t=50 μm; (c)—t=100 μm. |
为了更好地理解极薄带轧制变形的位错滑移机制，对轧制过程中滑移系的启动进行了分析.晶体的[1 0 0]方向与X轴一致，[0 0 1]和[0 1 0]分别与Y轴和Z轴相对应.由晶体塑性理论可知，滑移系是否启动由该滑移系上剪切应变率控制.轧制变形中，FCC金属12个滑移系的启动过程如图 6所示.极薄带轧制过程可以简化为二维平面应变变形，有利于变形沿轧制方向延伸的滑移系启动，4个主要滑移系满足条件：a3 (1 1 1)[-1 1 0],b2 (-1 1 1)[1 1 0],c2 (1 -1 1)[1 1 0] 和 d3 (1 1 -1)[-1 1 0].
图6(Fig. 6)
 | 图6 轧制过程 (1 1 1) 滑移系的启动 Fig. 6 Illustration of (1 1 1) slip plane for initial rolling progress |
图 7，图 8分别是厚度t=20 μm和t=100 μm铜箔20%压下率变形特定位置处滑移系启动和运动过程.剪切应变率的负值表示滑移系沿主滑移系和次滑移系的反方向进行滑移.A,B,C和J,K,L点分别位于t=20 μm和t=100 μm铜箔的上表面、中心和晶界.由图可知，各个位置上激活的滑移系个数和运动状态都是不一样的.对于t=20 μm的铜箔，A,B,C位置上4个主滑移系都被激活运动且一直持续到变形结束，c2滑移系的剪切应变率一直保持正值，沿正方向进行滑移，而a3,b2,d3滑移系的剪切应变率一直保持负值，表示一直沿反方向进行滑移.在A处，经过一段时间变形后，可以观察到新的滑移系a1和d1被激活，这两个新激活的滑移系具有相反的滑移方向.在B处，除4个主滑移系被激活外，4个次滑移系也被激活，a1,c3沿其正方向滑移，b1,d1沿其反方向滑移.在C处，a3被激活后保持较为稳定的阶段，经过变形后，主滑移系的运动引起了次滑移系a1和d1的被激活.表面和晶界启动滑移系的剪切应变率大于晶粒中心，表明滑移集中发生在这些区域.对于t=100 μm的铜箔，J,K,L位置上4个次滑移系a1,b1,c3和d1被激活运动，a1和c3沿反方向滑移，b1和d1沿正方向滑移.4个次滑移系的剪切应变率在变形过程具有明显的波动性.在J处，4个主滑移系和4个次滑移系的剪切应变率呈对称分布，4个滑移系沿正方向滑移，另外4个滑移系沿其反方向滑移.这种相反的滑移引起在铜箔上表面J处形成了应力和应变的集中.变形一段时间后，次滑移系的运动引起了4个主滑移系的被激活.在K处，滑移系的剪切应变率不再呈对称分布.4个次滑移系a1,b1,c3和d1激活后迅速增强，而后迅速减弱，保持稳定.在变形的中间阶段，可以明显观察到新的次滑移系a2和d2被激活.在L处，4个次滑移系的剪切应变率呈对称分布，但4个主滑移系的运动则没有对称性.变形一段时间后，次滑移系的运动引起了4个主滑移系的运动，a3,b2,d3沿正方向滑移，而c2沿反方向滑移.极薄带轧制变形的不均匀性和力学性能与滑移系的运动有密切关系，不同位置处被激活的主滑移系和次滑移系差异很大，同时滑移系激活后的运动状态也不同，正是这种差异性造成了变形的不均匀性和不同的力学响应.
图7(Fig. 7)
 | 图7 图 4a中A, B, C点上激活滑移系的剪切应变率 Fig. 7 The shear strain rate evolution of active slip systems at point A, B and C in Fig. 4a (a)—A点； (b)—B点; (c)—C点. |
图8(Fig. 8)
 | 图8 图 4c中J, K, L点上激活滑移系的剪切应变率 Fig. 8 The shear strain rate evolution of active slip systems at point J , K and L in Fig. 4c(a)—J点； (b)—K点; (c)—L点. |
图 9为两种厚度铜箔20%压下率变形区内晶粒a3和c2的剪切应变率演化云图.由图可知，滑移在晶界处受阻，晶粒内部和自由表面的剪切应变率的绝对值较大，说明在这些区域发生了显著地滑移.滑移系优先在晶界处启动，然后引起晶粒内部和自由表面滑移系的启动，滑移在晶界处受阻，形成滑移塞积，在垂直晶界处形成硬化区；由于自由表面对滑移系运动的阻碍比晶界小得多，滑移系最终在上下表面终止运动，形成软化区.对于t=20 μm的铜箔，a3和c2在变形开始时被激活，分别沿a3和-c2的方向滑移，由于t/d较小，形成了贯穿整个厚度方向的滑移剪切变形带.在厚度方向存在剪切应变率分布不均匀的现象，滑移系在铜箔的上下自由表面停止运动.对于t=100 μm的铜箔，剪切应变率分布更加不均匀，不仅相邻晶粒间出现方向相反的滑移，而且在单个晶粒内部也出现了方向相反的滑移，这种滑移机制造成了晶间和晶内剪切变形带的形成.晶内出现方向相反的滑移，可促进亚晶的形成.滑移系在晶界处被激活，而后引起晶内滑移系的启动，最终在自由表面停止运动.可见，晶粒取向对滑移系启动和运动状态具有重要影响，晶粒取向不同，造成滑移系激活顺序和滑移方向不同，从而引起相邻晶粒发生方向不同的滑移；保持晶粒尺寸不变，随铜箔厚度的增加，应变率分布的不均匀性增大，滑移系在晶内和自由表面都被激活滑移.这种滑移系启动位置和方向的差异性造成了变形的不均匀分布和不同区域的力学响应.
图9(Fig. 9)
 | 图9 主要滑移系a3和c2的剪切应变率的演化 Fig. 9 The evolution of shear strain rate in the two principal activated slip systems of a3 and c2 (a)—t=20 μm； (b)—t=100 μm. |
3 结 论1) 采用晶体塑性有限元方法模拟了厚度变化对铜箔轧制不均匀变形的影响，将模拟得到的应力-应变曲线与实验测定的曲线进行对比，两者基本一致.计算了晶粒尺度下铜箔微观应力应变及滑移系分布.
2) 由于各个晶粒的取向不同以及相邻晶粒的取向差不同，晶粒在变形过程中不同位置处滑移系被激活的条件不同，滑移系的启动和运动也不同.
3) 材料内部变形很不均匀，应力应变集中区可能出现在晶界，也可能出现在晶粒内部.这种变形的不均匀性主要是由初始晶粒取向不同，相邻晶粒取向差以及变形时滑移系激活和运动特性不同引起的.
参考文献
| [1] | Klusemann B,Svendsen B,Vehoff H.Investigation of the deformation behavior of Fe-3% Si sheet metal with large grains via crystal plasticity and finite-element modeling[J].Computational Materials Science,2012,52(1):25-32.( 1) |
| [2] | Fan H D,Li Z H,Huang M S.Size effect on the compressive strength of hollow micropillars governed by wall thickness[J].Scripta Materialia,2012,67(3):225-228.(1) |
| [3] | Hill R,Rice J R.Constitutive analysis of elastic-plastic crystals at arbitrary strain[J].Journal of the Mechanics and Physics of Solids,1972,20(6):401-413.(1) |
| [4] | Zhang P,Karimpour M,Balint D,et al.A controlled Poisson Voronoi tessellation for grain and cohesive boundary generation applied to crystal plasticity analysis[J].Computational Materials Science,2012,64:84-89.(1) |
| [5] | Kim H K,Oh S K.An interfacial energy incorporated couple stress crystal plasticity and the finite element simulation of grain subdivision[J].Journal of the Mechanics and Physics of Solids,2012,60(11):1815-1841.(1) |
| [6] | Rawat S,Chandra S,Chavan V M,et al.Integrated experimental and computational studies of deformation of single crystal copper at high strain rates[J].Journal of Applied Physics,2014,116(21):213507.(1) |
| [7] | Parasiz S A,Kinsey B,Krishnan N,et al.Investigation of deformation size effects during microextrusion[J].Journal of Manufacturing Science and Engineering,2007,129(4):690-697.(1) |
| [8] | Kim G Y,Ni J,Koc M.Modeling of the size effects on the behavior of metals in microscale deformation processes[J].Journal of Manufacturing Science and Engineering,2007,129(3):470-477.(1) |
| [9] | Wang Y,Dong P L,Xu Z Y,et al.A constitutive model for thin sheet metal in micro-forming considering first order size effects[J].Materials & Design,2010,31(2):1010-1014.(1) |
| [10] | Lu H N,Wei D B,Jiang Z Y,et al.Modelling of size effects in microforming process with consideration of grained heterogeneity[J].Computational Materials Science,2013,77:44-52.(1) |
| [11] | Peirce D,Asaro R J,Needleman A.Material rate dependence and localized deformation in crystalline solids[J].Acta Metallurgical,1983,31(12):1951-1976.(1) |
| [12] | Bassani J L,Wu T Y.Latent hardening in single crystals Ⅱ.analytical characterization and predictions[J].Proceedings of the Royal Society A,1991,435(1893):21-41.(1) |
本文献在全文中的定位：
... 的各向异性增加，材料的局部性能有可能超过材料宏观结构的平均效应[1].当箔材尺寸降低到微尺度时，由于尺寸效应的出现，宏观连续的变形 ...
1
本文献在全文中的定位：
... ，由于尺寸效应的出现，宏观连续的变形理论不能直接应用到微成形中[2].金属箔材厚度尺寸与其晶粒尺寸具有相同的数量级，就意味着在箔材 ...
1
本文献在全文中的定位：
... 形梯度张量来描述位错滑移的宏观效应，能够反映材料的微观结构特性[3]，被广泛地用来研究单晶体、多晶体的变形行为和织构演化等 ...
1
本文献在全文中的定位：
... ，被广泛地用来研究单晶体、多晶体的变形行为和织构演化等4, 5, 6 ...
1
本文献在全文中的定位：
... 、多晶体的变形行为和织构演化等[4, 5, 6].Parasiz等 ...
1
本文献在全文中的定位：
... [4, 5, 6].Parasiz等[7]发现当晶粒尺寸减 ...
1
本文献在全文中的定位：
... .Parasiz等[7]发现当晶粒尺寸减小到与试样尺寸同一数量级时，材料的变形行为主要 ...
1
本文献在全文中的定位：
... 行为主要受到单个晶粒性能的影响，导致非均匀应力和硬度的分布.Kim等[8]和Yun等[9]通过修正经典的霍 ...
1
本文献在全文中的定位：
... 和硬度的分布.Kim等[8]和Yun等[9]通过修正经典的霍尔-佩奇公式来描述薄带拉伸过程中应力-应变变化.尽 ...
1
本文献在全文中的定位：
... 应力-应变关系与实验结果一致，但他们都没有考虑晶粒取向的影响.Lu等[10]建立了一种新的微观变形模型，发现保持试样尺寸不变的情况下，材料 ...
1
本文献在全文中的定位：
... 采用Peirce[11]提出的率相关晶体塑性本构模型.整个晶体变形分解为首先发生位错沿 ...
1
本文献在全文中的定位：
... 采用Bassani和Wu[12]提出的适合FCC金属的硬化模型，其表达式为 ...
