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

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

doi:10.11858/gywlxb.20190836

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

doi: 10.11858/gywlxb.20190836

赵伟业^{1, 2, 3,},
赵聃^{1, 2, 3,*, ,},
吕品^{1, 2, 3},
金涛^{1, 2, 3},
马胜国^{1, 2, 3}

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

基金项目: 国家自然科学基金（11602158, 11802199, 11572214）；山西省自然科学青年基金（201601D021026）；山西省“1331工程”重点创新团队

详细信息

作者简介:
赵伟业（1992－），男，硕士研究生，主要从事多晶体有限变形数值计算研究.E-mail：939029385@qq.com

通讯作者:
赵　聃（1985－），男，博士，讲师，主要从事多晶体有限变形算法及多尺度本构理论研究.E-mail：zhaodan@tyut.edu.cn

中图分类号: O344.1
计量
- 文章访问数: 8225
- HTML全文浏览量: 3004
- PDF下载量: 30
出版历程
- 收稿日期: 2019-09-19
- 修回日期: 2019-10-25
- 刊出日期: 2020-02-25

Finite Element Calculation of Polycrystalline Shear-Compression Specimens with Static Loading

ZHAO Weiye^{1, 2, 3
,},
ZHAO Dan^{1, 2, 3,*
, ,},
LÜ Pin^{1, 2, 3},
JIN Tao^{1, 2, 3},
MA Shengguo^{1, 2, 3}

1.
Institute of Applied Mechanics, School of Mechanical and Vehicle Engineering, Taiyuan University of Technology, Taiyuan 030024, Shanxi, China
2.
Shanxi Key Laboratory of Material Strength & Structural Impact, Taiyuan 030024, Shanxi, China
3.
Taiyuan University of Technology Mechanics National Experimental Teaching Demonstration Center, Taiyuan 030024, Shanxi, China

摘要

摘要: 基于晶体塑性理论研究了晶体织构对数值计算结果的影响，建立了带有织构的多晶体压剪试样（SCS）模型。从材料和试样结构两方面研究了静态加载条件下微观晶粒在有限变形过程中对试样宏观力学性能的影响。由于模型几何结构的特殊性，重点对模型斜槽部分的应力、应变及变形特点进行了分析。考虑到试样在压缩过程中受摩擦的影响，数值分析了不同摩擦系数对变形过程的影响，在此基础上计算了相同摩擦系数下不同晶粒数目、不同单元数目以及单元类型对多晶体压剪模型力学性能的影响，并对试件关键部位不同取向晶粒的应力状态进行了分析。
- 多晶体 /
- 晶体织构 /
- 晶体塑性 /
- 压剪 /
- 摩擦系数
Abstract: The effect of crystal texture on the numerical results was studied based on the theory of crystal plasticity, and the polycrystalline compression shear sample (SCS) model with texture was established. The influence of micro-grain on the macroscopic mechanical properties in the process of finite deformation under static loading condition was studied in terms of material and sample structure. Because of the particularity of model geometry, the stress, strain and deformation characteristics of skewed slot were computed. Considering the effect of friction on the specimen during compression, the influence of friction coefficients on the deformation process was analyzed numerically. The influences of grain number, element number and element type on the mechanical properties of polycrystalline compression shear model under the same friction coefficient were calculated. The stress states of grain with different orientations in key parts of the specimen were also studied.
- polycrystalline /
- crystal texture /
- crystal plasticity /
- shear-compression /
- friction coefficient

HTML全文

多晶体是由许多不同取向的单晶体构成，晶粒的取向对金属材料的塑性变形特征以及损伤、断裂等都有很大影响^[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}} = {{{F}}^{{*}}}{{{F}}^{\rm{p}}}$

(1)

式中： ${{{F}}^{*}}$ 为弹性变形梯度，包含弹性变形和刚性转动； ${{{F}}^{\rm{p}}}$ 是体现晶体滑移的塑性变形梯度。速度梯度的表达式为

${{L}} = {\dot{ F}}{{{F}}^{{{ - 1}}}} = {{{L}}^{*}} + {{{L}}^{\rm{p}}} = {{\dot{ F}}^{*}} \cdot {{{F}}^{{{*}^{{{ - 1}}}}}} + {{{F}}^{*}} \cdot {{\dot{ F}}^{\rm{p}}} \cdot {{{F}}^{{{\rm{p}}^{{{ - 1}}}}}} \cdot {{{F}}^{{{*}^{{{ - 1}}}}}}$

(2)

式中： ${{{L}}^{*}}$ 为速度梯度的弹性部分， ${{{L}}^{\rm{p}}}$ 为速度梯度的塑性部分。

${{\dot{ F}}^{\rm{p}}} \cdot {{{F}}^{{{\rm{p}}^{{{ - 1}}}}}} = \sum\limits_\alpha{{{\dot {\gamma} }^{({\alpha} )}}{{{s}}_{({\alpha} )}} \otimes{{{m}}_{({\alpha} )}}}$

(3)

式中： ${\dot \gamma} ^{({\alpha )}}$ 为第 $\alpha$ 个滑移系上的剪切应变率， ${{{s}}_{({\alpha} )}}$ 为初始构型下第 $\alpha$ 个滑移系滑移方向上的单位向量， ${{{m}}_{({\alpha} )}}$ 为初始构型下第 ${\alpha}$ 个滑移系滑移面法线方向的单位向量。晶格畸变后，滑移面的法向量与滑移方向上的向量可以表示为

${{s}}_{({\alpha} )}^{*} = {{{F}}^{*}} \cdot {{{s}}_{({\alpha} )}},\quad {{m}}_{({\alpha} )}^{*} = {{{m}}_{({\alpha} )}} \cdot {{{F}}^{{*}{{{ - 1}}}}}$

(4)

式中： ${{s}}_{({\alpha} )}^{*}$ 和 ${{m}}_{({\alpha} )}^{*}$ 为相互正交的非单位向量。在当前构型下，速度梯度可表达成其对称部分 ${{D}}$ （变形率张量）和反对称部分 ${{\varOmega }}$ （旋率张量）之和的形式

${{L}} = {\dot{ F}} \cdot {{{F}}^{ {{- 1}}}} = {{D}} + {{\varOmega}}$

(5)

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

${{D}} = \frac{1}{2}({{L}} + {{{L}}^{\rm{T}}}) = {{{D}}^*} + {{{D}}^{\rm{p}}},\quad {{\varOmega}} = \frac{1}{2}({{L}} - {{{L}}^{\rm{T}}}) = {{{\varOmega}} ^*} + {{{\varOmega}} ^{\rm{p}}}$

(6)

并且满足

${{{D}}^*} + {{{\varOmega}} ^*} = {{\dot{ F}}^*} \cdot {{{F}}^*}^{ - 1},\quad {{{D}}^{\rm{p}}} + {{{\varOmega}} ^{\rm{p}}} = \sum\limits_\alpha{{{\dot {\gamma} }^{(\alpha )}}{{s}}_{(\alpha )}^* \otimes {{m}}_{(\alpha )}^*}$

(7)

1.2 本构模型

在晶体变形过程中，弹性变形率 ${{{D}}^{{*}}}$ 和Cauchy应力张量的Jaumann率 $\mathop {{{\bf{\sigma }}^{\bf{*}}}}\limits^{\nabla}$ ^[20]表示为

$\mathop {{{{\sigma}} ^*}}\limits^{\nabla} + {{\sigma}} ({{I}}:{{{D}}^*}) = {{C}}:{{{D}}^*}$

(8)

式中： ${\bf{\sigma }}$ 为Cauchy应力，I为二阶单位张量， ${{C}}$ 为四阶对称弹性模量张量。

根据Cauchy应力的本构方程， $\mathop {{{\bf{\sigma }}^{\bf{*}}}}\limits^{\nabla}$ 可以表示为

$\mathop {{{{\sigma}} ^*}}\limits^{\nabla} = \mathop {{\sigma}} \limits^{\nabla} + ({{\varOmega}} - {{{\varOmega}} ^*}) \cdot {{\sigma}} - {{\sigma}} \cdot ({{\varOmega}} - {{{\varOmega}} ^*})$

(9)

式中： $\mathop {{\sigma }}\limits^{\nabla}$ 为共旋应力率， ${{{\varOmega }}^{\bf{*}}}$ 对应于旋率张量的弹性部分，其中

$\mathop {{\sigma}} \limits^{\nabla} = \dot {{\sigma}} - {{\varOmega}} \cdot {{\sigma}} + {{\sigma}} \cdot {{\varOmega}}$

(10)

分切应力 ${{\tau} ^{(\alpha )}}$ 与Cauchy应力之间有如下关系

${\tau ^{(\alpha )}} = {{\sigma }}:({ {{s}}_{({\alpha} )}^{*}} \otimes { {{m}}_{({\alpha} )}^{*}})$

(11)

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

1.3 率相关硬化模型

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

${\dot \gamma ^{(\alpha )}} = {\dot \gamma _0}{\rm{sgn}}\;{\tau ^{(\alpha )}}{\left| {\frac{{{\tau ^{(\alpha )}}}}{{{{{g}}^{(\alpha )}}}}} \right|^{1/m}}$

(12)

式中： ${\dot \gamma _0}$ 为参考切应变率， $\tau ^{\left( \alpha \right)}$ 为第 $\alpha$ 滑移系上的分解切应力，m为率敏感系数， $g^{(\alpha)}$ 为临界分切应力。当m=0时，式（12）对应于黏弹性材料；当 $m \to \infty$ 时，式（12）对应于率无关材料。

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

${\dot {g}^{(\alpha )}} = \sum\limits_\beta {{{h}_{\alpha \beta}} {\dot \gamma ^{(\beta )}} }$

(13)

式中： ${{h}_{\alpha \beta }}$ 为硬化矩阵。 ${{h}_{\alpha \alpha }}$ （ $\alpha = \beta$ ）表示自硬化模量， ${{h}_{\alpha \beta}}$ （ $\alpha \ne \beta$ ）表示潜硬化模量。

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

${{h}_{\alpha \alpha }} = {h}(\gamma ) = {{h}_0}\sec {{h}^2}\left| {\frac{{{{h}_0}\gamma }}{{{\tau _{\rm{s}}} - {\tau _0}}}} \right|$

(14)

式中： ${{h}_0}$ 为初始硬化模量； ${{\tau} _0}$ 为临界切应力； ${{\tau} _{\rm{s}}}$ 为饱和切应力； $\gamma$ 为所有滑移系上切应变的Taylor累积积分，即总切应变，表示为

$\gamma = \sum\limits_\alpha {\int\limits_0^{t} {\left| {\dot \gamma ^{(\alpha )}} \right|{\rm{d}}t} }$

(15)

潜硬化模量表示为

${{h}_{\alpha \beta}} = {qh}(\gamma )\;\;\;\;(\alpha \ne \beta)$

(16)

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

2. 有限元数值模拟

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

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

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

图 1 多晶体随机织构和丝织构的极图

Figure 1. Polar graphs of polycrystalline with random texture and fiber texture

下载: 全尺寸图片幻灯片

表 1 6061铝合金材料参数^[23]

Table 1. Material parameters of 6061 aluminum alloy^[23]

C₁₁/GPa	C₁₂/GPa	C₄₄/GPa	m	q	${ { {\dot {\gamma } } }_{{0}}}$ /s^–1	h₀/MPa	τ₀/MPa	τ_s/MPa
108.200	61.300	28.500	20	1.1	1	96	23.5	54

下载: 导出CSV

| 显示表格

图 2 具有不同初始织构的单向压缩应力-应变曲线

Figure 2. Unidirectional compression stress-strain curves with different initial textures

下载: 全尺寸图片幻灯片

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

2.2 创建多晶体压剪模型

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

图 3 多晶体压剪试样几何尺寸

Figure 3. Geometrical dimension of polycrystal shear-compression specimen

下载: 全尺寸图片幻灯片

图 4 多晶体压剪模型

Figure 4. Polycrystal shear-compression model

下载: 全尺寸图片幻灯片

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

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

图 5 不同摩擦系数情况下的应力云图

Figure 5. Stress nephogram with different friction coefficients

下载: 全尺寸图片幻灯片

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

表 2 不同摩擦系数对应模型的数值结果比较

Table 2. Numerical results of models corresponding to different friction coefficients

Friction coefficient	Maximum shear strain	Maximum normal strain	Maximum equivalent strain	Maximum force/N
0.025	0.740	−0.331	0.741	8 899.73
0.050	0.697	−0.263	0.625	9 566.97
0.100	0.651	−0.212	0.553	10 001.50
Relative maximum difference	13.7%	56.1%	34.0%	12.4%

下载: 导出CSV

| 显示表格

图 6 不同摩擦系数情况下模型斜槽单元的平均切应力-平均切应变曲线

Figure 6. Average shear stress-average shear strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 7 不同摩擦系数情况下模型斜槽单元的平均正应力-平均正应变曲线

Figure 7. Average normal stress-average normal strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 8 摩擦系数不同时模型斜槽单元的平均Mises应力-平均等效应变曲线

Figure 8. Average Mises stress-average equivalent strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 9 摩擦系数不同时模型顶面的力-位移曲线

Figure 9. Force-displacement curve of the top surface of the model corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

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

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

表 3 特征晶粒欧拉角

Table 3. Euler angle of characteristic grains

Grain number	φ₁/(°)	$\psi$ /(°)	φ₂/(°)	Grain number	φ₁/(°)	$\psi$ /(°)	φ₂/(°)
1	69.24	167.34	8.61	9	334.93	38.19	29.33
2	44.05	92.39	45.56	10	30.45	127.11	356.19
3	356.63	154.92	76.61	11	59.93	105.73	313.11
4	305.74	23.07	333.37	12	88.45	166.05	33.32
5	45.14	101.42	64.90	13	74.86	154.71	302.25
6	282.88	86.40	282.11	14	76.35	51.48	307.37
7	350.37	146.78	286.51	15	79.65	77.10	351.57
8	49.77	3.76	329.80	16	357.28	56.00	72.22

下载: 导出CSV

| 显示表格

图 10 模型斜槽部分一条棱上所选取的16个单元

Figure 10. 16 selected elements on an edge of the chute part of the model

下载: 全尺寸图片幻灯片

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

图 11 不同摩擦系数下特征晶粒在压缩方向的应变

Figure 11. Strain of characteristic grain in compression direction under different friction coefficients

下载: 全尺寸图片幻灯片

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

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

图 12 包含不同晶粒数目的压剪有限元模型

Figure 12. Finite element models of compression shear with different grain numbers

下载: 全尺寸图片幻灯片

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

图 13 不同晶粒数目对应模型的载荷-位移曲线

Figure 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。

图 14 单元数目不同的模型

Figure 14. Models with different number of elements

下载: 全尺寸图片幻灯片

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

图 15 单元数目不同时模型的Mises变形云图

Figure 15. Mises deformation nephograms of models with different numbers of elements

下载: 全尺寸图片幻灯片

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

表 4 单元数目不同时模型在不同压缩距离下的载荷及其最大相对偏差

Table 4. Loads of the model with different numbers of elements at different compression distance and their maximum relative differences

Number of elements	Load/N
Number of elements	1.0 mm	1.4 mm	1.8 mm	2.5 mm
8 026	−7 790.82	−8 180.42	−8 469.96	−8 726.07
14 604	−7 487.38	−7 824.37	−8 008.19	−8 157.61
23 835	−7 669.64	−8 038.55	−8 270.62	−8 436.06
Maximum relative difference	4.05%	4.55%	5.77%	6.97%

下载: 导出CSV

| 显示表格

图 16 单元数目不同时模型的载荷-位移曲线

Figure 16. Force-displacement curves of models with different numbers of elements

下载: 全尺寸图片幻灯片

图 17 斜槽相同位置处晶粒的Mises应力-等效应变曲线

Figure 17. Mises stress-equivalent strain curve of grains at the same position in the chute

下载: 全尺寸图片幻灯片

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

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

图 18 两种单元类型对应的模型斜槽部分Mises变形云图

Figure 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 两种单元类型对应的模型顶面力-位移曲线

Figure 19. Top surface force-displacement curve of model with two element types

下载: 全尺寸图片幻灯片

图 20 两种单元类型对应的模型斜槽部分单元的平均Mises应力-平均等效应变曲线

Figure 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个晶粒），数值结果并未出现明显差异；而固定晶粒数目、改变单元数目时，由于单元尺寸的变化导致不同组模型、相同位置的晶粒形状有所差异，晶粒间的相互作用出现差异，从而对数值结果产生了一些影响。不同单元类型的计算结果表明，六面体单元的变形程度大于四面体单元，计算精度更高。在相同压缩量下六面体单元模型所需压力比四面体单元小，但计算时间较长。若考虑计算成本，可选择单元数目相对较少的四面体单元对多晶体宏观力学响应进行预估。

图 1 多晶体随机织构和丝织构的极图

Figure 1. Polar graphs of polycrystalline with random texture and fiber texture

下载: 全尺寸图片幻灯片

图 2 具有不同初始织构的单向压缩应力-应变曲线

Figure 2. Unidirectional compression stress-strain curves with different initial textures

下载: 全尺寸图片幻灯片

图 3 多晶体压剪试样几何尺寸

Figure 3. Geometrical dimension of polycrystal shear-compression specimen

下载: 全尺寸图片幻灯片

图 4 多晶体压剪模型

Figure 4. Polycrystal shear-compression model

下载: 全尺寸图片幻灯片

图 5 不同摩擦系数情况下的应力云图

Figure 5. Stress nephogram with different friction coefficients

下载: 全尺寸图片幻灯片

图 6 不同摩擦系数情况下模型斜槽单元的平均切应力-平均切应变曲线

Figure 6. Average shear stress-average shear strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 7 不同摩擦系数情况下模型斜槽单元的平均正应力-平均正应变曲线

Figure 7. Average normal stress-average normal strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 8 摩擦系数不同时模型斜槽单元的平均Mises应力-平均等效应变曲线

Figure 8. Average Mises stress-average equivalent strain curve of the model’s chute element corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 9 摩擦系数不同时模型顶面的力-位移曲线

Figure 9. Force-displacement curve of the top surface of the model corresponding to different friction coefficients

下载: 全尺寸图片幻灯片

图 10 模型斜槽部分一条棱上所选取的16个单元

Figure 10. 16 selected elements on an edge of the chute part of the model

下载: 全尺寸图片幻灯片

图 11 不同摩擦系数下特征晶粒在压缩方向的应变

Figure 11. Strain of characteristic grain in compression direction under different friction coefficients

下载: 全尺寸图片幻灯片

图 12 包含不同晶粒数目的压剪有限元模型

Figure 12. Finite element models of compression shear with different grain numbers

下载: 全尺寸图片幻灯片

图 13 不同晶粒数目对应模型的载荷-位移曲线

Figure 13. Force-displacement curves of models corresponding to different numbers of grains

下载: 全尺寸图片幻灯片

图 14 单元数目不同的模型

Figure 14. Models with different number of elements

下载: 全尺寸图片幻灯片

图 15 单元数目不同时模型的Mises变形云图

Figure 15. Mises deformation nephograms of models with different numbers of elements

下载: 全尺寸图片幻灯片

图 16 单元数目不同时模型的载荷-位移曲线

Figure 16. Force-displacement curves of models with different numbers of elements

下载: 全尺寸图片幻灯片

图 17 斜槽相同位置处晶粒的Mises应力-等效应变曲线

Figure 17. Mises stress-equivalent strain curve of grains at the same position in the chute

下载: 全尺寸图片幻灯片

图 18 两种单元类型对应的模型斜槽部分Mises变形云图

Figure 18. Mises deformation nephogram of the model with two element types

下载: 全尺寸图片幻灯片

图 19 两种单元类型对应的模型顶面力-位移曲线

Figure 19. Top surface force-displacement curve of model with two element types

下载: 全尺寸图片幻灯片

图 20 两种单元类型对应的模型斜槽部分单元的平均Mises应力-平均等效应变曲线

Figure 20. Average Mises stress-strain curves of the sloped part of the model with two element types

下载: 全尺寸图片幻灯片

表 1 6061铝合金材料参数^[23]

Table 1. Material parameters of 6061 aluminum alloy^[23]

C₁₁/GPa	C₁₂/GPa	C₄₄/GPa	m	q	${ { {\dot {\gamma } } }_{{0}}}$ /s^–1	h₀/MPa	τ₀/MPa	τ_s/MPa
108.200	61.300	28.500	20	1.1	1	96	23.5	54

下载: 导出CSV

表 2 不同摩擦系数对应模型的数值结果比较

Table 2. Numerical results of models corresponding to different friction coefficients

Friction coefficient	Maximum shear strain	Maximum normal strain	Maximum equivalent strain	Maximum force/N
0.025	0.740	−0.331	0.741	8 899.73
0.050	0.697	−0.263	0.625	9 566.97
0.100	0.651	−0.212	0.553	10 001.50
Relative maximum difference	13.7%	56.1%	34.0%	12.4%

下载: 导出CSV

表 3 特征晶粒欧拉角

Table 3. Euler angle of characteristic grains

Grain number	φ₁/(°)	$\psi$ /(°)	φ₂/(°)	Grain number	φ₁/(°)	$\psi$ /(°)	φ₂/(°)
1	69.24	167.34	8.61	9	334.93	38.19	29.33
2	44.05	92.39	45.56	10	30.45	127.11	356.19
3	356.63	154.92	76.61	11	59.93	105.73	313.11
4	305.74	23.07	333.37	12	88.45	166.05	33.32
5	45.14	101.42	64.90	13	74.86	154.71	302.25
6	282.88	86.40	282.11	14	76.35	51.48	307.37
7	350.37	146.78	286.51	15	79.65	77.10	351.57
8	49.77	3.76	329.80	16	357.28	56.00	72.22

下载: 导出CSV

表 4 单元数目不同时模型在不同压缩距离下的载荷及其最大相对偏差

Table 4. Loads of the model with different numbers of elements at different compression distance and their maximum relative differences

Number of elements	Load/N
Number of elements	1.0 mm	1.4 mm	1.8 mm	2.5 mm
8 026	−7 790.82	−8 180.42	−8 469.96	−8 726.07
14 604	−7 487.38	−7 824.37	−8 008.19	−8 157.61
23 835	−7 669.64	−8 038.55	−8 270.62	−8 436.06
Maximum relative difference	4.05%	4.55%	5.77%	6.97%

下载: 导出CSV

参考文献(24)

[1]	董湘怀, 仲町英治. 晶体塑性模型在板材成形计算机模拟中的应用 [J]. 中国机械工程, 1997, 8(4): 27–30, 118. DONG X H, NAKAMACHI E. Application of crystal plastic model in computer simulation of sheet metal forming [J]. China Mechanical Engineering, 1997, 8(4): 27–30, 118.
[2]	刘海军, 方刚, 曾攀. 基于晶体塑性理论的大变形数值模拟技术 [J]. 塑性工程学报, 2006, 13(2): 1–8, 28. doi: 10.3969/j.issn.1007-2012.2006.02.001 LIU H J, FANG G, ZENG P. Numerical simulation of large deformation based on the theory of crystal plasticity [J]. Journal of Plastic Engineering, 2006, 13(2): 1–8, 28. doi: 10.3969/j.issn.1007-2012.2006.02.001
[3]	皮华春, 韩静涛, 薛永栋. 金属塑性成形的晶体塑性学有限元模拟研究进展 [J]. 机械工程学报, 2006, 42(3): 15–21. doi: 10.3321/j.issn:0577-6686.2006.03.003 PI H C, HAN J T, XUE Y D. Progress in finite element simulation of metal plasticity forming with crystal plasticity [J]. Journal of Mechanical Engineering, 2006, 42(3): 15–21. doi: 10.3321/j.issn:0577-6686.2006.03.003
[4]	司良英. FCC金属冷加工织构演变的晶体塑性有限元模拟 [D]. 沈阳: 东北大学, 2009. SI L Y. Finite element simulation of crystal plasticity in FCC metal cold working texture evolution [D]. Shenyang: Northeastern University, 2009.
[5]	司良英, 邓关宇, 吕程. 基于Voronoi图的晶体塑性有限元多晶几何建模 [J]. 材料与冶金学报, 2009, 8(3): 193–197, 216. SI L Y, DENG G Y, LÜ C. Crystal plastic finite element polycrystalline geometry modeling based on Voronoi diagram [J]. Journal of Materials and Metallurgy, 2009, 8(3): 193–197, 216.
[6]	张丰果, 董湘怀. 微塑性成形模拟材料细观建模 [J]. 模具技术, 2011(3): 16–19. ZHANG F G, DONG X H. Microplastic forming simulates material microscopic modeling [J]. Mould Technology, 2011(3): 16–19.
[7]	郑文, 徐松林, 蔡超. 基于Hopkinson压杆的动态压剪复合加载实验研究 [J]. 力学学报, 2012, 44(1): 124–131. ZHENG W, XU S L, CAI C. Experimental study on dynamic compression shear composite loading based on Hopkinson compression bar [J]. Journal of Mechanics, 2012, 44(1): 124–131.
[8]	章超, 徐松林, 王道荣. 花岗岩动静态压剪复合加载实验研究 [J]. 固体力学学报, 2014, 35(2): 115–123. ZHANG C, XU S L, WANG D R. Experimental study on dynamic and static compressor-shear composite loading of granite [J]. Journal of Solid Mechanics, 2014, 35(2): 115–123.
[9]	李雪艳, 李志斌, 张舵. 闭孔泡沫铝准静态压剪性能研究 [J]. 高压物理学报, 2018, 32(3): 52–59. doi: 10.11858/gywlxb.20170655 LI X Y, LI Z B, ZHANG D. Study on quasi-static compressive shear properties of aluminum foam with closed-cell [J]. Journal of High Pressure Physics, 2018, 32(3): 52–59. doi: 10.11858/gywlxb.20170655
[10]	RITTEL D, LEE S, RAVICHANDRAN G. A Shear-compression specimen for large strain testing [J]. Experimental Mechanics, 2002, 42(1): 58–64. doi: 10.1007/BF02411052
[11]	RITTEL D, RAVICHANDRAN G, LEE S. Large strain constitutive behavior of OFHC copper over a wide range of strain rates using the shear compression specimen [J]. Mechanics of Materials, 2002, 34(10): 627–642. doi: 10.1016/S0167-6636(02)00164-3
[12]	DOROGOY A, RITTEL D. A numerical study of the applicability of the shear compression specimen to parabolic hardening materials [J]. Experimental Mechanics, 2006, 46(3): 355–366. doi: 10.1007/s11340-006-6414-8
[13]	DOROGOY A, RITTEL D. Numerical validation of the shear compression specimen. Part Ⅱ: dynamic large strain testing [J]. Experimental Mechanics, 2005, 45(2): 178–185. doi: 10.1007/BF02428191
[14]	DOROGOY A, RITTEL D. Numerical validation of the shear compression specimen. Part I: quasi-static large strain testing [J]. Experimental Mechanics, 2005, 45(2): 167–177. doi: 10.1007/BF02428190
[15]	DOROGOY A, RITTEL D, GODINGER A. Modification of the shear-compression specimen for large strain testing [J]. Experimental Mechanics, 2015, 55(9): 1627–1639. doi: 10.1007/s11340-015-0057-6
[16]	VURAL M, MOLINARI A, BHATTACHARYYA N. Analysis of slot orientation in shear-compression specimen (SCS) [J]. Experimental Mechanics, 2010, 51(3): 263–273.
[17]	ZHAO J, KNAUSS W G, RAVICHANDRAN G. A new shear-compression-specimen for determining quasistatic and dynamic polymer properties [J]. Experimental Mechanics, 2008, 49(3): 427–436.
[18]	PIERCE D, ASARO R J, NEEDLEMAN A. Material rate sensitivity and localized deformation in crystalline solids [J]. Acta Metall, 1983, 31(12): 1951–1976. doi: 10.1016/0001-6160(83)90014-7
[19]	ASARO R J, NEEDLEMAN A. Overview No. 42 texture development and strain hardening in rate dependent polycrystals [J]. Acta Metallurgica, 1985, 33(6): 923–953.
[20]	HUANG Y. A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program: MECH 178 [R]. Harvard University, 1991.
[21]	王国军, 孙强. 4032合金热挤压棒材变形行为及形变织构 [J]. 黑龙江冶金, 2013, 33(2): 1–5, 8. WANG G J, SUN Q. Deformation behavior and deformation texture of hot extruded 4032 alloy bar [J]. Heilongjiang Metallurgy, 2013, 33(2): 1–5, 8.
[22]	HUANG S Y, ZHANG S R, LI D Y. Simulation of texture evolution during plastic deformation of FCC, BCC and HCP structured crystals with crystal plasticity based finite element method [J]. Transactions of Nonferrous Metals Society of China, 2011, 21(8): 1817–1825. doi: 10.1016/S1003-6326(11)60936-9
[23]	辛存, 赵聃, 闫晓鹏. 材料三维微结构表征及其晶体塑性有限元模拟 [J]. 计算力学学报, 2018, 36(2): 233–239. XIN C, ZHAO D, YAN X P. 3D modeling microstructure and crystal plasticity finite element simulation [J]. Journal of Computational Mechanics, 2018, 36(2): 233–239.
[24]	齐康, 闫昊, 陈祥瑶. 利用ABAQUS模拟不同模态下的金属切削过程 [J]. 机械工程与自动化, 2018(2): 93–94. QI K, YAN H, CHEN X Y. Metal cutting processes in different modes are simulated by ABAQUS [J]. Mechanical Engineering and Automation, 2018(2): 93–94.