考虑位错密度的铝合金动态压缩晶体塑性有限元分析

赵英杰，胡 广，赵 聃，张团卫，马胜国，王志华

(太原理工大学 a.机械与运载工程学院，应用力学研究所，b.材料强度与结构冲击山西省重点实验室，c.力学国家级实验教学示范中心，太原 030024)

7075系铝合金因其具有较高的比强度和较低的密度，在航空工业、结构材料、航天飞行器和汽车部件等方面得到了广泛的应用，因此高应变率条件下对其力学性能的研究具有重要的研究意义[1]。分离式霍普金森压杆(SHPB)试验装置是一种用于表征材料在高应变率下力学特性常用的动态测试技术[2]。通过SHPB可以获得较宽应变率范围以及不同温度下材料的力学响应。KIM et al[3]利用SHPB研究了具有不同长径比的7075-T6铝合金在高应变率下的力学行为，观察到其应力、应变行为受长径比的影响。SASSO et al[4]使用SHPB测试了不同初始回火(T6和O回火)对7075铝合金的高应变率行为，发现T6回火表现为脆性，而O的回火行为是韧性的，在高应变率下，二者均表现出高应变率下正应变率敏感。CHUNG et al[5]在室温下通过SHPB对Fe27Co24Ni23Cr26合金的动态力学性能进行了研究，结果表明高应变率变形诱导微米孪晶、纳米孪晶和亚纳米孪晶组成分级结构，极大地提高了材料的强度和韧性。此外，SUZANE et al[6]采用SHPB对不同回火温度的HHA(High Hardness Armor)钢的高应变率变形行为及其微观结构演化进行了研究，发现不同回火温度的HHA钢均出现了表现率敏感的特性，并在变形过程中形成了绝热剪切带。

相比于传统宏观模型，晶体塑性有限元方法(CPFEM)通过考虑位错滑移、孪晶以及相变等微观机制对金属塑性变形和微观结构演化进行了准确地描述[7]。有限元模型中网格每个积分点上包含一个或多个单晶晶粒并与特定的晶体取向相关联，通过描述单晶变形本构关系在积分点上的运行，进而实现了用有限单元描述复杂的连续体变形[8]。ZHANG et al[9]将内热的形成理论耦合到动态CPFEM模型中，成功预测了HCP多晶金属材料微观绝热剪切带的形成，研究发现热扩散和塑性耗散相关常数之间的平衡以及微观织构对控制应变局部化至关重要。LI et al[10]通过建立具有不同初始织构多晶模型的动态单轴压缩和单轴剪切的晶体塑性有限元模型，研究了FCC金属在高应变率载荷下典型织构对塑性变形的影响，其预测的变形行为与实验结果大致相同。本文从宏观和微观尺度出发，发展了考虑位错密度的跨尺度晶体塑性模型，采用CPFEM与SHPB实验相结合的方式，研究了7075-T6铝合金动态加载下的宏观和微观力学响应，探讨了其动态变形行为以及微观结构的演化规律。

1 晶体塑性理论

1.1 运动学

晶体的变形由弹性变形和塑性变形组成，其中弹性变形主要为晶格的畸变，塑性变形则包括晶体的位错滑移、孪晶等。根据连续介质力学的理论[11]，假设材料初始构型为dX，当前构型为dx，则二者可以通过变形梯度张量F进行关联：

dX=FdX.

(1)

其中，变形梯度张量可以根据晶体变形的机制进一步分解为弹性变形梯度和塑性变形梯度：

F=FeFp.

(2)

式中：Fe为弹性变形梯度，包含弹性变形及刚体转动等可恢复的变形；Fp为塑性变形梯度，描述由位错滑移和变形孪晶等产生的非弹性剪切变形等不可恢复变形。根据晶体塑性理论可以认为材料首先由初始构型dX经塑性变形Fp转变至中间构型，随后通过弹性变形Fe转变至当前构型dx，该转变过程示意如图1所示。

(3)

图1 变形梯度的分解Fig.1 Multiplicative decomposition of deformation gradient

(4)

L=D+Ω.

(5)

式中：D为变形率张量；Ω为旋率张量。需要注意的是变形率张量为对称张量，旋率张量为反对称张量。与变形梯度张量相类似，变形率张量和旋率张量也可分解为弹性部分和塑性部分：

D=De+Dp,Ω=Ωe+Ωp.

(6)

1.2 本构模型

晶体内的弹性本构关系可根据文献[12]的研究得出：

(7)

α滑移系的塑性剪切变形率可以根据率相关晶体塑性模型[13-14]表示为如下形式：

(8)

1.3 硬化模型

图2 滑移面上的林位错和运动位错Fig.2 Forest dislocation and moving dislocation on slip system

基于位错密度演化的硬化模型假设位错运动的阻碍来源于林位错，林位错的空间分布密度ρα与滑移阻力gα的关系为：

(9)

(10)

式中：Gαβ为位错密度相互作用矩阵G中的系数，其中存在G0～G5六个强度因子，它们之间相互独立，代表不同滑移系之间的相互作用。其中，G0为滑移系中位错的自硬化作用，G1为伯氏矢量不同的滑移面平行位错间的相互作用，二者均描述了滑移面平行(包括共面)的位错间相互作用强度，其余4个系数则描述了非共面位错间的相互作用强度。ARSENLIS et al[15]通过对FCC金属的研究得出了FCC金属中6个强度因子分别为：G0=0.10，G1=0.22，G2=0.30，G3=0.38，G4=0.16，G5=0.45.

位错密度不仅存在增殖，同时由于滑移系的相互作用还会产生位错密度湮灭的现象，据此KOCKS et al[16]提出了Kocks-Mecking(K-M)模型用以描述位错密度与剪切应变率的依赖关系：

(11)

式中：ρα为第α个滑移系的位错密度；k1表示位错的增殖系数；k2表示位错的湮灭系数。将其与式(10)联立可以得到通过位错密度表示的当前强度变化率与剪切应变率之间的关系为：

(12)

上式为以位错密度为中间变量的滑移阻力变化率与剪切应变率的关系，即基于位错密度的硬化模型。

2 结果与分析

2.1 动态压缩实验

通过分离式霍普金森压杆对7075-T6铝合金制成的直径6 mm，长12 mm的圆柱形试件进行动态压缩实验，实验设备示意如图3所示。试件放置于入射杆和透射杆之间，通过控制气压以调控子弹的速度从而达到对试件压缩应变率的控制，入射杆和透射杆上的中心处对称位置分别粘贴两个电阻应变片，组成半桥电路以消除两杆本身弯曲的影响，经过动态应变仪采集到应变片电阻的变化，并通过示波器将电阻的变化转换成电压的变化并进行显示，通过计算机的进一步处理可以得到应变的变化。

图3 分离式霍普金森压杆示意Fig.3 Diagram of split Hopkinson pressure bar

虽然本文所研究的7075-T6铝合金在材料内部晶粒存在不同取向，但从宏观角度而言其各部分力学性能基本一致，符合材料的均匀性假设，所以本文采用二波法[17]对实验数据进行处理。

(13)

图4 压缩后试件的工程应力-应变曲线Fig.4 Engineering stress-strain curve of sample after compression

为了对本文所用7075-T6铝合金材料的初始微观结构及变形后的微观结构进行表征，获取其织构信息，在初始圆柱形试件及实验后试件上分别按照图5所示方式取出厚度为1 mm的试样，阴影面为其扫描面，样品坐标系如下图所示，选取其压缩方向为RD方向。

图5 EBSD样品示意图Fig.5 Diagram of EBSD simple

对取出后样品做如下处理：1) 将样品表面依次使用200目至2 000目的SiC砂纸进行初步打磨，随后依次使用2 000目至4 000目的水磨砂纸对其进行进一步打磨，直至表面光亮且不存在肉眼可见的划痕；2) 将通过砂纸打磨好的样品依次使用W2.5至W0.5的金刚石抛光膏在机械抛光机上进行进一步的抛光处理，抛光至在光学显微镜下观察不到明显划痕为止；3) 将机械抛光后的样品在体积比为9∶1的乙醇高氯酸溶液中进行电解抛光，电解在25 V的直流电压下进行，时间为15 s，并将电解后样品在乙醇溶液中清洗，取出后使其自然风干。使用热场发射扫描电子显微镜(JSM-7100F)配备的EBSD探针(HKL NordlysNano)对处理后样品进行扫描，加速电压和扫描步长分别设置为20 kV和0.5 μm.图6(a)为实验后得到铝合金初始状态的反极图(inverse pole figure，IPF)，图6(b)为实验前后试件的极图(pole figure，PF)对比。本文所采用铝合金由拉拔成型工艺制作，由图6(a)可知，本文所采用的铝合金材料内部晶粒尺寸沿拉拔方向为300 μm左右，沿横向为100 μm左右，其晶粒分布颜色可知，其晶粒的初始取向基本一致，以<001>方向和<111>方向为主，同时包含少量其他取向晶粒，大多数晶粒沿拉拔方向变形为纤维状。当其经过动态压缩后，材料内部晶粒发生旋转，由初始方向向α-fiber方向发生聚集。

图6 EBSD结果Fig.6 Results of EBSD

2.2 晶体塑性模型参数分析

由于本文所采用的动态隐式算法需要较小的分析步和较长的计算时间，综合考虑计算精度和时间成本，本文将动态压缩实验简化为刚性面板对圆柱形试样的压缩，试件尺寸与前文实验中相同，直径为6 mm，长12 mm.网格类型设置为六面体八节点线性减缩积分单元(C3D8R)，模型中一个网格代表一个晶粒，通过Python程序将EBSD得到的7075-T6铝合金的初始织构输入三维宏观多晶模型中，通过实验中得到的工程应力-应变曲线对不同参数的影响进行标定。由于本节主要分析某个变量对宏观应力-应变曲线的影响，所以只对分析时的参数做控制变量的分析，即仅保证其他参数在合理取值范围内且保持不变，并未使用确实符合实际的参数进行分析。

由于实验中试件与入射杆和透射杆之间均存在接触，所以在CPFEM模拟中应对其接触时摩擦力的影响进行分析比较。文献[18]的研究表明，当接触面较为光滑时，铝合金的摩擦系数为0.02左右，本文选取摩擦系数分别为0.01、0.02和0.05的3种模型进行分析，其工程应力-应变曲线如图7所示。由图可以看出，对于0.02左右的摩擦系数而言，相同模型的工程应力-应变曲线变化极小，说明摩擦系数对模拟结果力学性能的影响可以忽略不计。此外摩擦系数还会对存在接触位置的变形产生影响，由图7中的小图接触面对数应变分布云图可以看出，不同摩擦系数接触面的变形基本一致，三者直径变化的相对误差为1.5% 左右。综合上述分析，在后续CPFEM的模拟中将不再考虑摩擦系数对结果的影响。

图8(a)和图8(b)分别为不同位错增殖系数k1及不同位错湮灭系数k2的模型中工程应力-应变曲线的变化情况，可以看出在其他参数不变的情况下，k1和k2对材料的屈服强度并不产生影响，仅与其屈服后的硬化行为有关。随着k1的增加，材料的加工硬化行为更加显著，就微观结构而言，这是由于k1的增加导致晶粒内部位错增殖速度上升，在材料发生塑性变形后位错大量增殖，由于位错间的相互作用导致材料内晶体位错的运动更加困难，宏观上即表现为材料的硬化行为加剧且极限强度提高。与k1相反，k2的增加则会导致材料极限强度下降且降低其塑性阶段的硬化行为，这是由于当k2的值增加时，材料塑性变形过程中的位错密度湮灭速度也随之增加，晶体内部位错运动的阻力减小，宏观上则表现为极限强度和硬化行为的减弱。

图7 不同摩擦系数模型工程应力-应变曲线Fig.7 Engineering stress-strain curves of different models with different friction coefficients

图8 不同模型工程应-应变曲线Fig.8 Engineering stress-strain curves of different models

2.3 参数标定

根据上文对相关参数的分析，在对动态压缩实验进行CPFEM模拟时，本文选用网格尺寸为0.3 mm的模型进行计算，摩擦系数设置为0.02.通过对模拟中得到的工程应力-应变曲线与实验中曲线比较以调整晶体塑性模型中的相关参数，当两条曲线可以较好地相符时，则认为此时的参数能够准确描述实际实验中试件的宏观力学行为及微观结构变化。图9为使用最终标定的参数计算时CPFEM模拟与实验的工程应力-应变曲线对比，由于分离式霍普金森压杆的二维弥散效应[19]，动态实验中得到的工程应力-应变曲线存在较大波动，一般取其震荡的中心值作为实验结果。本文进行了3次相同应变率的动态压缩实验，通过3次实验所得到铝合金的工程应力-应变曲线具有较好的重复性，说明通过SHPB对铝合金进行动态压缩实验是十分可靠的。需要注意的是，CPFEM模拟结果与实验结果在弹性段有较大的误差，这是由于动态压缩实验中材料迅速进入塑性段，处于弹性段的时间较短，由于数据采集和环境影响，导致实验中的弹性段数据并不完全准确，本文所采用晶体塑性模型中，通过立方晶系中3个相互独立的弹性常数C11、C12和C44[20]决定材料的弹性性质。由此可以看出使用CPFEM模拟7075-T6铝合金动态压缩实验得出的工程应力-应变曲线在塑性阶段能够与实验中得到的曲线基本相符，说明拟合后得到的参数与实际材料参数能够很好地吻合，表1为7075-T6铝合金CPFEM模型中的相关参数标定值，其中ρint为7075-T6铝合金初始位错密度[21]，μ为其剪切模量[22]。

图9 CPFEM模型与动态压缩实验的工程应力-应变曲线对比Fig.9 Comparison of the engineering stress-strain curves in dynamic compression test between CPFEM model and experiment

2.4 织构演化分析

本文通过ABAQUS中的UMAT子程序对晶体塑性有限元模型进行计算时，可以获得计算过程中每一步中的状态变量(solution-dependent state variables，SDV)，包括试件内部每个晶粒的取向及其位错密度等，据此，可以通过Python程序包将所需晶体取向的变化信息进行收集和处理，并将相应晶体的欧拉角信息输入JTEX软件中进行处理，便可将CPFEM模拟得到的晶粒取向变化结果与实验结果加以比较。动态压缩实验后得到的试件极图与CPFEM模拟中试件极图的对比(见图10)，可以看出二者可以较好地相符，但织构强度存在一定差异，这是由于实际实验中对材料的表征以及计算中边界条件的设置均存在一定的误差，但总体而言本文所研究的基于位错密度的晶体塑性有限元模型可以较为准确地预测材料的宏观力学行为以及微观结构的演化趋势。为了进一步分析7075-T6铝合金材料动态压缩后织构的演化情况，并验证CPFEM模拟计算结果的准确性，绘制了数值计算结果和实验结果的取向分布函数图(orientation distribution function，ODF)进行对比。

表1 7075-T6铝合金CPFEM模型相关参数Table 1 Related parameters of CPFEM model of 7075-T6 aluminum alloy

图10 CPFEM数值计算与实验试样极图Fig.10 Pole figure of sample in CPFEM numerical calculation and experiment

通过JTEX软件导出初始状态、实验后以及数值计算结果的ODF如图11所示，对应体积分数柱状图如图12所示，导出时误差角设置为15°[23]。由图可以看出，7075-T6铝合金在初始状态时织构类型主要为Cube {001}<100>织构以及一定量的Goss {110}<001>织构、A {110}<111>织构和Cu {112} <111>织构，其中Cube织构和Goss织构为铝合金中典型的再结晶织构。本文所采用的拉拔成型的7075-T6铝合金在初始状态时产生了<111>方向和<001>方向的聚集。动态压缩后由于试件内部晶粒发生了取向的择优转动，晶粒的织构类型发生了相应的变化，Cube{001}<100>织构、A {110}<111>织构、P{011}<211>织构和Cu {110}<001>织构的体积分数相较初始状态材料下降，同时生成了较多的Brass{110}<112>织构和Goss {110}<001>织构，这些织构的体积变化趋势与吴晓冬[24]对7003铝合金动态压缩后的微观结构变化相符。同时根据CPFEM数值计算结果与实验结果的对比可知，除S{123}<634>织构和F{111}<112>织构外，CPFEM模拟可以对动态压缩过程中的织构体积分数变化趋势作出较为准确地预测。

图11 实验、CPFEM计算结果及初始状态材料的ODFFig.11 ODF figure of experiment, CPFEM numerical calculation, and initial stage material

图12 实验、CPFEM数值计算结果及初始状态材料中的典型织构体积分数Fig.12 Volume fraction of ideal texture in experiment, CPFEM numerical calculation, and initial stage material

3 结论

1) 摩擦系数对试件宏观力学性能的影响较小，只对其最终变形模态的显示有一定影响。

2) 位错增殖系数k1增大会导致材料的硬化行为加剧且极限强度提高，位错湮灭系数k2增大则会导致材料的极限强度下降，且弱化其硬化行为。

3) 动态压缩后由于试件内部晶粒发生了取向的择优转动，晶粒的织构类型发生了相应的变化，生成了较多的Brass{110}<112>织构和Goss {110}<001>织构。

4) CPFEM模型可以较好地预测动态压缩过程中试样内部织构的变化趋势，但对织构体积分数的预测与实验存在一定差异。