目录结构

The specific coefficients of thermal viscoplastic constitutive equation for TB2 were obtained by use of the least square method and C coded programs.The method of estimate adiabatic temperature proposed was used to estimate the temperature rising in adiabatic shear band(ASB).The estimate value is in agreement with the results obtained by FEM simulation. The estimation result of adiabatic temperatures rising indicates that the temperature in ASB is between the recrystallization temperature and phase transition point, and the result is also in agreement with the characteristic of the microstructure in ASB.

Keyword：

adiabatic shearing; thermal viscoplastic constitutive relation; TB2 alloy; adiabatic temperatures rising; FEM;

Received： 2002-12-02

绝热剪切现象普遍存在于爆炸复合、高速撞击、侵彻、冲孔、切削、高速成型、冲蚀等涉及冲击载荷的高速变形过程中 ^[1,2,3,4] 。由于整个变形过程是在很短的时间内完成的, 90%塑性功转化为热能并且来不及散失, 所以近似认为在高应变速率(>10³/s)下的变形过程为绝热过程。材料在高应变速率条件下的变形特点由3个互相竞争的因素来决定: 1) 由于变形程度的增加而导致的应变强化效应; 2) 应变率强化效应; 3) 由于绝热温升而导致的热软化效应。当热软化效应占优时, 材料就会发生热粘塑性本构失稳 ^[3] , 而形成所谓绝热剪切带(adiabatic shear band, ASB) ^[3,4,5] 。

因而, 绝热温升是决定材料是否发生热粘塑性本构失稳的重要因素。但事实上, 温度是影响一切物理冶金过程的重要因素。材料在绝热剪切变形中, 由于存在绝热温升, ASB内的温度可能达到或者超过其再结晶温度(0.3~0.5 T_m)。研究 ^{[2,6,7,8,9,10]} 发现, ASB内的晶粒由原始-12 μm减小至-0.1 μm, 并且晶粒呈等轴状, 晶粒内位错密度低。 ASB内的组织结构特征表明, ASB内己发生了动态再结晶, 然而在如此短的时间内, 晶粒急剧细化而发生动态再结晶的机制尚不清楚。绝热温升在其中扮演着重要角色, 然而目前的实验手段尚无法直接测量其分布规律和量比数据。

TB2合金是一种对绝热剪切相当敏感的合金, 其流变应力对应变速率有较大的敏感性, 在高应变速率形变条件下表现出粘塑性。本文作者以TB2合金为研究对象, 采用Johnson-Cook热粘塑性本构模型 ^[2,7] , 利用最小二乘法 ^[11,12,13] 拟合本构方程中的待定参数, 根据热功转换原理估算绝热温升, 在DEFORM-2D中建立有限元计算模型 ^[14,15] , 模拟TB2合金在绝热剪切条件下的应力、应变、绝热温升等, 为进一步模拟ASB内组织演化奠定基础。

1TB2合金的热粘塑性本构模型

本文采用的Johnson-Cook热粘塑性本构关系 ^[2] 形式为: $σ = (σ_{0} + E_{1} ε) (1 + g l n \frac{\dot{ε}}{\dot{ε}_{0}}) (1 - α \frac{Τ - Τ_{0}}{Τ_{0}})$ (1)

式中 σ₀、 $\dot{ε}_{0}$ 、 T₀分别是由室温准静态试验确定的TB2合金的特征应力、特征应变率和参考温度; E₁、 g、 α分别是表征TB2合金的应变强化、应变率强化和热软化的参数。

通过分离式Hopkinson压杆(SHPB)试验得到的数据 ^[2] 如图1所示。

图1 TB2在不同应变速率、不同温度下的真实应力—应变关系 Fig.1 Relationship between stress and strain under different strain rates and temperatures

结合实验数据, 利用最小二乘法拟合这3个参数的值。式(1)可表示为式(2)所示的一般形式 ^[13] 。

y=f(ζ₁, ζ₂, …, ζ_k, θ₁, θ₂, …, θ_p)+ε

=f(ζ, θ)+ε (2)

若获得了n组的独立观测值: (ζ_i1, ζ_i2, …, ζ_ik, y_i) i=1, 2, …, n

则有模型

${\begin{cases} y = f (ζ_{i 1}, ζ_{i 2}, ?, ζ_{i k}, θ_{1}, θ_{2}, ?, θ_{p}) + ε_{i} \\ ? ? = f_{i} (θ) + ε_{i} ? i = 1, 2, ?, n \\ ε_{i} ? Ν (0, σ^{2}) \end{cases} ? ? ? (3)$

现在问题是要求参数θ的最小二乘法估计 $\overset{\land}{θ}$ 的求法, 式(3)记为

$S (θ) = \sum_{i} [y_{i} - f_{i} (θ)]^{2} ? ? ? (4)$

而θ的最小二乘法估计 $\overset{\land}{θ}$ 应满足

$S (\overset{\land}{θ}) = \min_{θ} S (θ)$

如果f对θ的每一分量均可微, 则求 $\overset{\land}{θ}$ 相当于解下述正规方程组:

$\frac{? S (θ)}{? θ_{j}} |_{θ = \overset{\land}{θ}} = 0 ? j = 1, 2, ?, p$

即:

$\begin{array}{l} \sum_{i} [y_{i} - f_{i} (θ)] \frac{? S (θ)}{? θ_{j}} |_{θ = \overset{\land}{θ}} = 0 \\ j = 1, 2, ?, p ? ? ? (5) \end{array}$

式中 y_i为y变量的第i组实测值; f_i(θ)为y变量的第i组估计值。

用高斯-牛顿迭代法 ^[13] 解方程组(5)。迭代公式的一般形式为

θ^(l+1)=θ^(l)+β^(l)l=0, 1, 2, …

求β^(l)先将f_i(θ)在θ^(l)处作泰勒展开, 近似有

$f_{i} (θ) = f_{i} (θ^{(l)}) + \sum_{j = 1}^{p} \frac{? f_{i} (θ)}{? θ_{j}} |_{θ = θ^{(l)}} (θ_{j} - θ_{j}^{(l)})$

式中 p为参数的个数。

若

$\begin{array}{l} S S (θ) = \sum [y_{i} - f_{i} (θ^{(l)} - \sum_{j = 1}^{p} \frac{? f_{i} (θ)}{? θ_{j}} |_{θ = θ^{(l)}} \\ ? (θ_{j} - θ_{j}^{(l)})]^{2} \\ = \sum_{i} [y_{i} - f_{i} (θ^{(l)} - \sum_{j = 1}^{p} x_{i j}^{(l)} β_{j}^{(l)}]^{2} ? ? ? (6) \end{array}$

达到最小, 则求出了参数的最小二乘估计。

将式(6)写为矩阵形式, 则

β^(l)=(X^(l)′X^(l))^-1X^(l)′Y

l=0, 1, 2, … (7)

式中 $Y^{(l)} = [\begin{array}{l} y_{1} - f_{1} (θ^{(l)}) \\ y_{2} - f_{2} (θ^{(l)}) \\ ? ? \\ ? ? \\ y_{n} - f_{n} (θ^{(l)}) \end{array}], X^{(l)} = [\begin{array}{l} x_{11}^{(l)} x_{12}^{(l)} ? ? x_{1 p}^{(l)} \\ x_{21}^{(l)} x_{22}^{(l)} ? ? x_{2 p}^{(l)} \\ ? ? \\ x_{n 1}^{(l)} x_{n 2}^{(l)} ? ? x_{n p}^{(l)} \end{array}], β^{(l)} = [\begin{array}{l} β_{1}^{(l)} \\ β_{2}^{(l)} \\ ? \\ ? \\ β_{p}^{(l)} \end{array}] = [\begin{array}{l} θ_{1} - θ_{1}^{(l)} \\ θ_{2} - θ_{2}^{(l)} \\ ? \\ ? \\ θ_{p} - θ_{p}^{(l)} \end{array}], β_{j}^{(l)} = θ_{j} - θ_{j}^{(l)}, x_{i j}^{(l)} = \sum_{j = 1}^{p} \frac{? f_{i} (θ)}{? θ_{j}} |_{θ = θ^{(l)}}$ 。

对任意给定的δ>0, 当|S(θ^(l+1))-S(θ^(l))|<δ时, 可取 $\overset{\land}{θ} = θ^{(l + 1)}$ 。用C语言实现上述求解过程, 通过程序读入实验数据和精度控制条件, 运用程序即可得到拟合参数的结果, 从而得到热粘塑性本构方程的具体表达式。

2 有限元模拟

为了模拟TB2合金在绝热剪切条件下的内部应力、应变及温度分布情况, 采用如图2所示的“帽形”试样 ^[7] 作为在DEFORM-2D中进行模拟计算模型。

图2 用于有限元模拟的试样 Fig.2 Specimen for FEM simulation

由于试样为轴对称, 模拟时只取一半进行分析计算, 如图3所示, 帽形试样放置于上下模之间。上模以不同的速度向下压缩试样, 在试样中产生不同的应变速率, 下模保持静止, 试样与上模之间就保持从主运动关系。在图3所示的模型中, 上面的小圆柱与下面的大圆环相交部位产生强烈的剪切变形而形成变形局域化区域, 因而在进行模拟计算之前, 必须先定义好材料属性、生成网格、设定变形边界条件与温度边界条件等。表1, 2 ^[2] 分别列出了TB2合金的材料属性和材料模拟参数。

材料属性中, 除了表1, 2中已列出的参数外,

图3 试样与模子关系示意图 Fig.3 Schematic diagram of relationship between specimen and dies

还要设定其本构模型。 DEFORM-2D中提供了几种常用的本构关系通式, 不能用这些通式表示的, 可以在一定范围内直接输入若干组数据, 这样系统会根据这些数据自动进行插值。本模拟中, 选定温度范围是20~500 ℃、应变速率范围是(1.0~2.0)×10⁴/s、应变范围是0.02~0.13, 将这些条件代入通过拟合得到的热粘塑性本构方程中, 得到对应的应力值, 并将得到的每一组应力、应变、应变速率、温度数据以表格的形式输入到DEFORM-2D中。本模拟采用四边形网格, 其中试样网格数设置为2000, 上模、下模网格数均设置为200, 并设上、下模为刚体。图3中的箭头所指的边为模型的对称轴, 其在x轴上没有位移, 故设其速度为0作为变形边界条件, 温度边界条件是试样初始温度为20

表1 TB2的材料属性 Table 1 Material properties of TB2

Elastic modulus/ MPa	Poisson ratio		Density/ (g·cm^-3)		Specific heat/ (J·g^-1·℃^-1)
1.15×10⁵	0.21		4.81		0.542
Thermal conductivity/ (J·cm^-1·s^-1·℃^-1)		Linear expansion coefficient/ ℃^-1		Melting point/ ℃		Phase transition point/℃
0.126		8.53×10^-6		1 610		750

表2 TB2的材料模拟参数 Table 2 Material simulation parameters of TB2

Environment temperature/℃	Velocity of top die/(m·s^-1)	Coefficient of friction
20	2	0.001

℃, 与空气接触的外表面温度保持初始温度不变。假设试样与模子之间的摩擦很小, 摩擦系数取为0.001。由于变形局域化区域高为2.03 mm, 所以设定上模向下的位移为1 mm时终止变形。

3 模拟结果与讨论

利用编制的C语言程序对式(1)中的待定参数进行拟合, 结果为: E₁=1 007.168 MPa, α=0.428 609, g=0.016 688。将参数值代回式(1), 即得到TB2合金的热粘塑性本构方程。将实测值与本构方程的模拟值进行比较, 发现该本构方程较好地表达了TB2合金的流变应力与应变、应变速率及温度的关系, 具有较高的精度。图4所示是在 $Τ = 293 Κ ? \dot{ε} = 1 500 / s$ 时, 模拟值与实测值的比较。

图4 流变应力的模拟值与实测值的比较 Fig.4 Comparison of simulated with measured flow stresses

DEFORM-2D的模拟计算结果如图5~8所示。图5表明, 最大应力分布于如图3中所示的变形局域化区域内。图6结果表明, 应变都集中在变形局域化区域内, 而基体基本上没有什么变形, 对于整个试祥来说, 在y轴方向上被压缩了1 mm, 所以其整体应变为|ln(14/15)|=0.07, 但其变形集中在上圆柱与下圆环相交部位, 即整个试样的单向压缩变形转化为此相交部位的剪切变形, 若剪切带的宽度为200 μm, 则可以这样来估计出其内的剪应

图5 帽形试样上的等效应力分布图 Fig.5 Distributions of equal stress in hat-shaped specimen

图6 帽形试样上的等效应变分布图 Fig.6 Distributions of equal strain in hat-shaped specimen

变 ^[9] , γ=1 mm/200 μm=5。图7表明, 在上模下行速度为2 m/s时, 其等效应变速率约1.5×10⁴/s, 同样可以通过用上模下行速度除以剪切带宽度来估计。图8表明, 基体内温度梯度很小, 与初始温度相差不大, 而剪切变形局域化区域内的温度明显高于基体温度, 且由外围到中心的温度不断升高, 因而试样两拐角处温度最高, 变形也最为严重, 但与整个变形局域化区域相比, 具体所占的体积百分比很小, 可以忽略不计。估算剪切变形局域化区域内的平均温度(从图8中的等温线B算到H)为501 ℃, 这是温度的有限元模拟结果。

图7 帽形试样上的等效应变速率分布图 Fig.7 Distributions of equal strain rate in hat-shaped specimen

图8 帽形试样上的温度分布图 Fig.8 Distributions of temperature in hat-shaped specimen

通过理论估算的方法来求绝热温升, 且与有限元模拟的结果进行比较, 验证其可靠性。因为上模对试样所作的功大部分转化为热量, 同时由于变形集中在局部区域内, 因而由变形转化而来的热量也将集中在这个区域, 并且导致该区域的温度升高。因此, 确定了应力应变关系, 就可以求出塑性功; 确定了变形集中的区域, 就可以求出热影响区的体积; 确定了密度、比热容等材料常数, 就可以由得到的热量算出其温升。从而得出

$Δ Τ = \frac{β W}{C_{v} ρ V} = F ? S = A \int_{0}^{δ_{e n d}} σ d δ ? ? ? (8)$

式中 W为外力所作的功; β为热功转化系数, 一般塑性功转化为热量的系数取0.9, 但为计算方便, 忽略弹性功的影响, 而从0到δ_end积分, 故取此系数为0.85; C_p为TB2合金的比热容; ρ为TB2合金的密度; V为剪切变形局域化区域的体积; A为试样与上模接触面的面积, 即为受力面积; σ为应力; δ为上模的位移。

若真应变值取正, 则位移与真应变之间的关系为

$ε = | \ln \frac{L - δ}{L} | ? e^{ε} = \frac{L}{L - δ} ? δ = L (1 - \frac{1}{e^{δ}})$

式中 L为试样的高度。

将式(8)中的位移替换为应变, 再带入本构方程积分得:

$\begin{array}{l} W = A \int_{0}^{ε_{e n d}} σ d (L - L e^{- ε}) \\ = A L (1 + g l n \frac{\dot{ε}}{\dot{ε}_{0}}) (1 - α \frac{Τ - Τ_{0}}{Τ_{0}}) \\ ? (σ_{0} - σ_{0} e^{- ε_{e n d}} - E_{1} ε_{e n d} e^{- ε_{e n d}} + \\ ? E_{1} - E_{1} e^{- ε_{e n d}}) ? ? ? (9) \end{array}$

式中 ε_end为位移δ_end所对应的真应变值。剪切变形局域化区域的体积可以根据图8估算, 将其简化为图9所示的阴影部分体积(一环状), 尺寸已经标出。

根据图6、 7所示的应变与应变速率场分别设定ε_end为0或1, 应变速率为1.5×10⁴/s, 环境温度设为20 ℃, 将体积、材料常数等代入式(8), 得出ΔT=477 ℃, 则剪切变形局域化区域内的温度T=477+20=497 ℃。这与有限元模拟的结果501 ℃非常接近, 表明本文作者根据Johnson-Cook模型建立的热粘塑性本构方程具有一定的实用价值, 用此方法估算绝热温升具有较高的可靠性。

图9 剪切变形局域化区域示意图 Fig.9 Schematic diagram of shearing localized region

4 结论

1) 采用Johnson-Cook热粘塑性本构模型, 用最小二乘法拟合出本构方程中的待定参数, 经验证得到的热粘塑性本构方程具有较高的精度, 反映了TB2合金对应变速率和温度敏感的流变特性。

2) 在DEFORM-2D中建立了模拟TB2合金绝热剪切行为的模型, 得到了TB2合金在高应变速率下的应力、应变、应变速率和温度的分布, 发现其最大应力位于上圆柱与下圆环相交的区域内, 并且应变、应变速率和温度分布也在这个区域内集中。估算了在应变速率为 1.5×10⁴/s、环境温度为20 ℃时, 剪切变形局域化区域内的平均温度为501 ℃。利用本文作者提出的方法, 估算了在此形变条件下, 剪切变形局域化区域内的温度为497 ℃, 说明建立的热粘塑性本构方程具有一定的实用价值, 估算的绝热温升具有较高的可信度。

3) 计算出的剪切变形局域化区域内的平均温度(501 ℃)达到了TB2合金的再结晶温度(483 ℃), 而低于其相转变温度(750 ℃), 这与作者前期所作的研究结果 ^[6,8] 相一致, 即其微观组织发生了动态再结晶而未发现相转变。