多晶体压剪试样静态加载有限元计算

赵伟业，赵 聃，吕 品，金 涛，马胜国

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

多晶体是由许多不同取向的单晶体构成，晶粒的取向对金属材料的塑性变形特征以及损伤、断裂等都有很大影响[1]。从微观机制上看，晶体塑性变形方式有位错滑移、孪生、相变等多种方式，晶体塑性理论主要分析上述机制影响下材料的变形规律[2]。由于各晶粒的取向不同，塑性变形首先发生在位错滑移容易开动的晶粒上，且晶粒间的相互作用协调各自的变形，导致多晶体内部的应变场和应力场具有强烈的不均匀性。晶体塑性有限元法是一种重要的分析晶粒间相互作用的方法[3]。国内一些学者基于多晶体材料的特点建立了Voronoi有限元模型，针对不同工况并结合材料自身性质，对多晶体模型进行了数值分析[4-6]。

在材料成型过程中，经常出现复杂的受力状态，须明晰复杂应力状态下材料的强度和变形特性，其中压剪复合加载导致的复合应力状态是国内外学者的关注点。郑文等[7]通过添加具有不同角度倾斜端面的垫块，研究了静态和动态压剪复合加载下材料的力学响应，给出了相应的数据处理方法，并利用有限元模拟验证了实验的可行性。章超等[8]基于郑文的实验设计方法，利用材料试验系统（MTS）静态实验装置和分离式霍普金森压杆（SHPB)动态实验装置，分别进行了不同角度斜端面垫块下花岗岩的准静态与动态加载试验，分析了不同加载路径下岩石的力学性能以及破坏情况。李雪艳等[9]通过添加单斜端面垫块和长方体套筒，利用材料试验机对闭孔泡沫铝进行准静态条件下的压剪复合加载实验，并通过改变垫块的角度得到泡沫铝在不同加载路径下的力学性能。国外Rittel等[10-11]提出带有两个与竖直方向呈45°的对称矩形斜槽的圆柱体压剪试样，用于研究较宽应变率及应变范围内试样的力学响应。Dorogoy等[12-15]利用准静态和动态两种加载方式，结合有限元模拟探究了压剪试样在较宽应变率范围内发生塑性应变时的力学性能，测试了不同应变率下材料的剪切力学响应，并将压剪试样的矩形斜槽变为半圆形斜槽，改善了压剪试样矩形斜槽边角的应力集中现象。Vural等[16]对压剪试样进行了数值分析和实验测试，建立了外部施加载荷、位移与斜截面内的等效应力、等效塑性应变等之间的关系式，在此基础上，为了增强斜截面应力场的均匀性，改变斜槽的角度并进行分析。Zhao等[17]将金属压剪试样改变为斜槽部分由聚合物代替的全新压剪试样，通过动态及准静态实验结合数值模拟，进行应力、应变分析，并将这种新试样与原来的全金属压剪试样进行全面对比。

上述研究大都假设材料为各向同性，未能体现出微观晶粒对宏观各向异性行为的影响。本研究引入晶体塑性理论，试图揭示多晶体材料的微观变形机理，以期提高材料加工成型产品的质量。将多晶体材料属性赋予压剪试样，结合材料的宏、微观特性进行有限元分析，并对压剪复合应力状态下微观晶粒的变形演化进行分析。

1 晶体塑性理论

1.1 运动学

晶体塑性理论的基本框架基于Perice、Asaro和Needlman[18-19]的研究工作。材料有限变形过程中，总体变形梯度F可分解为

式中：F∗为弹性变形梯度，包含弹性变形和刚性转动；Fp是体现晶体滑移的塑性变形梯度。速度梯度的表达式为

式中：L∗为速度梯度的弹性部分，Lp为速度梯度的塑性部分。

其中张量D和 Ω 也可以分解为弾性变形部分和塑性变形部分（上标*和p分别表示弹性和塑性），即

并且满足

1.2 本构模型

在晶体变形过程中，弹性变形率D∗和Cauchy应力张量的Jaumann率表示为

式中：σ为Cauchy应力，I为二阶单位张量，C为四阶对称弹性模量张量。

分切应力 τ(α)与Cauchy应力之间有如下关系

本研究以面心立方（FCC）晶体为研究对象。面心立方晶体包含4个滑移面，每个滑移面包含3个滑移方向，因此FCC晶体共有{111}〈110〉 12个滑移系。

1.3 率相关硬化模型

由Pierce等[18]提出的以Schmid准则为基础的率相关剪切应变率模型为

应变硬化通过临界分切应力的增量关系演化表示为

式中：hαβ为硬化矩阵。hαα（α = β）表示自硬化模量，hαβ（α ≠ β）表示潜硬化模量。

晶体变形时，不同硬化阶段的硬化矩阵不同，准确表达的难度很大，所以通常需要对硬化矩阵进行简化。对自硬化模量可做简单假设[20]

式中：h0为初始硬化模量； τ0为临界切应力；τs为饱和切应力；γ为所有滑移系上切应变的Taylor累积积分，即总切应变，表示为

潜硬化模量表示为

式中：q为硬化矩阵系数。

2 有限元数值模拟

2.1 材料初始织构对数值计算的影响

以6061铝合金棒材为研究对象进行数值计算。考虑到棒材轴向拉拔后会带有丝织构[21]，将对数值结果产生影响，为此使用Abaqus有限元软件和Python语言建立了包含54个晶粒的Voronoi代表性体积单元（Representative volume element, RVE）模型进行压缩计算，模型分别包含了初始随机织构和初始丝织构，其他条件完全相同。

图1给出了多晶体的两种初始取向，即随机织构和丝织构在（111）面上的极图[22]，利用表1所示的6061铝合金的材料参数[23]进行数值计算，单向压缩时的应力-应变（σ33-ε33）曲线计算结果如图2所示。表1中，C11、C12、C44为弹性模量分量，m为率敏感指数，q为硬化矩阵系数。

图1 多晶体随机织构和丝织构的极图Fig.1 Polar graphs of polycrystalline with random texture and fiber texture

表1 6061铝合金材料参数[23]Table 1 Material parameters of 6061 aluminum alloy[23]

由图2可知：当塑性变形较小时，初始织构为丝织构的模型在压缩方向上的应力绝对值高于随机取向对应模型，最大相对偏差约为7%；当塑性变形增加到一定程度时，两种模型对应的曲线基本重合。这表明当塑性变形较小时，初始织构为丝织构的模型受影响较为明显；随着塑性变形的增大，多晶体的晶体取向发生了一定程度的旋转，使得塑性变形较大时两种模型的应力比较接近。在后续的有限元计算中，均以具有初始丝织构的材料为研究对象。

图2 具有不同初始织构的单向压缩应力-应变曲线Fig.2 Unidirectional compression stress-strain curves with different initial textures

2.2 创建多晶体压剪模型

使用Abaqus有限元软件和Python语言建立多晶体模型，其中1个单元代表1个晶粒，创建两个平板刚体模型，平板刚体与多晶体模型上下底面进行接触设置，多晶体压剪试样几何尺寸及模型如图3和图4所示。接触属性为面面硬接触，并且设置相应的摩擦系数。给予下板参考点全约束边界条件；上板参考点处设置一个2.5 mm的压缩位移载荷，其他边界条件全约束。

2.3 摩擦力对数值模拟结果的影响

图3 多晶体压剪试样几何尺寸Fig.3 Geometrical dimension of polycrystal shear-compression specimen

考虑到压剪模型与刚体接触板之间摩擦的影响，首先对摩擦系数µ进行分析。根据齐康等[24]的总结，铝合金在良好润滑情况下摩擦系数可达到0.02左右。于是设置3组摩擦系数进行数值模拟，分别为0.025、0.050、0.100。每组摩擦系数对应的Mises应力分布如图5所示，可以看到试件内应力分布不均匀，随着摩擦系数的增加，斜槽以外部分的应力明显增加，且整个模型中间向外膨胀的变形更明显。每组模型变形后上板都会与压剪模型发生相对滑动，如黑色箭头所示。相对滑动距离分别为1.282、0.688、0.214 mm，其中摩擦系数最小组的相对滑动距离几乎是摩擦系数最大组的6倍，对斜槽及斜槽以外部分的交界面的变形产生了明显影响。

图4 多晶体压剪模型Fig.4 Polycrystal shear-compression model

图6和图7分别为模型斜槽部分单元的平均切应力-平均切应变曲线和平均正应力-平均正应变曲线。由两图可知，摩擦系数对切应变和正应变有明显影响，且最大切应变（见图6插图）和最大正应变（见图7插图）都与摩擦系数呈负相关。切应变的最大相对偏差为13.7%，正应变的最大相对偏差为56.1%。图8为模型斜槽部分对应单元的平均Mises应力-平均等效应变曲线，其最大等效应变随摩擦系数的变化规律与图6、图7相似，由插图知不同摩擦系数下等效应变的最大相对偏差为34.5%。图9为模型顶面的位移-载荷曲线，当模型压缩2.5 mm时，不同摩擦系数下所需载荷的最大相对偏差为12.4%，即摩擦系数越大，压缩模型时需要的载荷越大。对不同摩擦系数情况下开槽部位的最大剪切应变、最大正应变以及等效应变进行总结。结果列于表2。

图5 不同摩擦系数情况下的应力云图Fig.5 Stress nephogram with different friction coefficients

图6 不同摩擦系数情况下模型斜槽单元的平均切应力-平均切应变曲线Fig.6 Average shear stress-average shear strain curve of the model's chute element corresponding to different friction coefficients

图7 不同摩擦系数情况下模型斜槽单元的平均正应力-平均正应变曲线Fig.7 Average normal stress-average normal strain curve of the model's chute element corresponding to different friction coefficients

图8 摩擦系数不同时模型斜槽单元的平均Mises应力-平均等效应变曲线Fig.8 Average Mises stress-average equivalent strain curve of the model's chute element corresponding to different friction coefficients

图9 摩擦系数不同时模型顶面的力-位移曲线Fig.9 Force-displacement curve of the top surface of the model corresponding to different friction coefficients

表2 不同摩擦系数对应模型的数值结果比较Table 2 Numerical results of models corresponding to different friction coefficients

以上结果表明：随着摩擦系数的增加，模型与刚体接触板之间的滑动阻力增大；当3组模型都压缩2.5 mm时，摩擦系数的增加使模型与刚体接触板的相对滑动减小，斜槽部分沿切向和正向的变形都减小。因此，斜槽部分沿切向和正向的变形在竖直方向的分量之和减小，而斜槽以外部分在竖直方向上的压缩量增加。此分析结果与Dorogoy等[13]关于摩擦力的分析相吻合。因此，在模型压缩2.5 mm时，为了减少所需压力，增加斜槽的压剪变形程度，应当尽可能地减小摩擦系数。

为了进一步分析压剪模型关键区域的微观变形特征，在斜槽部分按照图10箭头所示方向选取了16个单元作为特征晶粒。各晶粒的初始欧拉角 (φ1,ψ,φ2)如表3 所示。

图10 模型斜槽部分一条棱上所选取的16个单元Fig.10 16 selected elements on an edge of the chute part of the model

表3 特征晶粒欧拉角Table 3 Euler angle of characteristic grains

图11显示了不同摩擦系数下16个晶粒在z方向上的应变变化。可以看出，大多数晶粒即使位置接近，但其应变值差异也较为明显，这是因为初始晶粒取向差较明显。图11中不同摩擦系数下晶粒1和晶粒16的应变绝对值较其他晶粒偏小，而且随着摩擦系数的增加而减小，原因是摩擦系数增加使模型中间部分产生更明显的涨型，斜槽两端在z方向上的变形减小，1号和16号晶粒位于斜槽端部，变形过程中会受到端末效应的影响，影响其变形。

图11 不同摩擦系数下特征晶粒在压缩方向的应变Fig.11 Strain of characteristic grain in compression direction under different friction coefficients

2.4 晶粒尺寸对数值结果的影响

多晶体有限元计算时，为了考虑晶粒尺寸对数值结果造成的影响，在单元数目不变的情况下建立了如图12所示的包含不同晶粒数目的多晶体计算模型，通过晶粒数目的差异体现晶粒尺寸的影响，模型上下刚体接触板与模型之间的摩擦系数设置为0.025，其他条件不变。

图13为不同晶粒数目模型的位移-载荷曲线，同时与1个单元表示1个晶粒的模型计算结果（Number of grains 20 000）进行对比。由图13可知，虽然晶粒数目不同，但其体现的初始拉伸织构基本一致，且几种模型计算的载荷结果在变化趋势上基本一致，各曲线间相差不足1%。

图12 包含不同晶粒数目的压剪有限元模型Fig.12 Finite element models of compression shear with different grain numbers

图13 不同晶粒数目对应模型的载荷-位移曲线Fig.13 Force-displacement curves of models corresponding to different numbers of grains

2.5 单元数目对数值结果的影响

为了探究单元数对多晶体压剪模型斜槽区域力学性能的影响，控制模型的晶粒数不变（固定为500），建立了3组单元数分别为8 026、14 604、23 835的模型进行数值分析，摩擦系数设置为0.025，单元类型为C3D8，其他条件不变，见图14。

图15为不同单元数目的3组模型的变形Mises应力云图，3组模型的变形情况基本一致。由于单元数目的变化，单元尺寸也随之改变，导致建模过程中，3组模型处于同一位置上的晶粒根据临近原则所包含的单元数目存在一定的差异，所以变形云图也产生了轻微的差异。

图14 单元数目不同的模型Fig.14 Models with different number of elements

图16为3组模型同时压缩2.5 mm时刚体接触板上参考点的位移-载荷曲线。可见，3组曲线轨迹基本相似，但是在相同位移下，载荷大小关于单元数目没有明显的规律，当模型压缩至2.5 mm时，3组模型载荷的最大相对偏差约6.97%，如表4所示。图17为3组模型斜槽部分相同位置处晶粒的Mises应力-等效应变曲线。由图17插图可见，单元尺寸及数目的变化导致3个模型的晶粒形状有一定差异，Mises应力云图的差异比较明显；3条曲线的弹性部分基本一致，塑性部分的差异较明显。以上表明，当单元数目或尺寸发生变化时，3组模型相同位置处晶粒的形状受到一定程度的影响，原因是指定位置的晶粒与周围晶粒间的相互作用发生了改变，周围晶粒的取向和形貌发生变化，进而影响该位置晶粒的应力和应变。

图15 单元数目不同时模型的Mises变形云图Fig.15 Mises deformation nephograms of models with different numbers of elements

图16 单元数目不同时模型的载荷-位移曲线Fig.16 Force-displacement curves of models with different numbers of elements

图17 斜槽相同位置处晶粒的Mises应力-等效应变曲线Fig.17 Mises stress-equivalent strain curve of grains at the same position in the chute

表4 单元数目不同时模型在不同压缩距离下的载荷及其最大相对偏差Table 4 Loads of the model with different numbers of elements at different compression distance and their maximum relative differences

2.6 单元类型对模拟结果的影响

为了探究单元类型的影响并提高模型计算效率，将斜槽区域单元数为2 093的四面体单元模型与斜槽区域单元数为2 040的六面体单元模型进行对比分析，两种模型除了单元类型外其他条件与前述相同。斜槽局部变形如图18所示，两种网格模型的斜槽整体变形相似，应力分布都不均匀，从微观的角度看六面体单元变形程度更大，且由云图颜色分布可知六面体单元的计算精度更高，表明六面体单元比四面体单元的刚度更小，更适合分析模型局部变形情况。

图18 两种单元类型对应的模型斜槽部分Mises变形云图Fig.18 Mises deformation nephogram of the model with two element types

图19、图20分别是两种不同类型单元的斜槽的位移-载荷曲线和Mises应力-平均等效应变曲线。可见，当模型被压缩到一定程度时，四面体单元的应力更大，位移为2.5 mm时，两条曲线的应力偏差为3.4%。六面体单元模型的最大平均等效应变比四面体单元模型大4.2%。这表明当整个压剪模型压缩2.5 mm时，六面体单元模型在承受更小压力的同时，模型斜槽部分单元发生了更大的变形。

图19 两种单元类型对应的模型顶面力-位移曲线Fig.19 Top surface force-displacement curve of model with two element types

图20 两种单元类型对应的模型斜槽部分单元的平均Mises应力-平均等效应变曲线Fig.20 Average Mises stress-strain curves of the sloped part of the model with two element types

3 结 论

(1)塑性变形较小时，初始织构为丝织构的RVE单向压缩模型相对于晶粒为随机取向的RVE单向压缩模型的应力值更高，当塑性变形增加到一定程度时，两种模型的应力值接近。这表明模型变形时晶粒取向会发生一定程度的转动，从而对数值结果造成影响。

(2)数值计算了3种不同摩擦系数对多晶压剪模型变形的影响。结果表明，随着摩擦系数的增加，模型与刚体接触板相对滑动减小，模型中间部分发生明显涨型，模型斜槽部分的压缩与剪切变形都减小，使得在压缩相同高度时，摩擦系数大的模型的斜槽部分沿竖直方向的压缩量减小，而斜槽以外部分沿竖直方向的压缩量增加，同时压缩模型所需要的力也增加。数值模拟结果显示，选取摩擦系数为0.025进行计算能够得到较好的压剪应力状态。从微观角度分析，由于晶粒初始取向不同导致Mises应力分布不均匀，且对于位置接近的晶粒，晶粒初始取向差的存在会导致其应变值差异明显。

(3)在相同摩擦条件下，分析固定单元数目而改变晶粒数目（尺寸）、固定晶粒数目而改变单元数目两种情况，并讨论了单元类型对数值结果的影响。结果表明：在单元数目确定的情况下，晶粒数目增加到一定程度时（1个单元代表1个晶粒），数值结果并未出现明显差异；而固定晶粒数目、改变单元数目时，由于单元尺寸的变化导致不同组模型、相同位置的晶粒形状有所差异，晶粒间的相互作用出现差异，从而对数值结果产生了一些影响。不同单元类型的计算结果表明，六面体单元的变形程度大于四面体单元，计算精度更高。在相同压缩量下六面体单元模型所需压力比四面体单元小，但计算时间较长。若考虑计算成本，可选择单元数目相对较少的四面体单元对多晶体宏观力学响应进行预估。