目录结构

详情信息展示

参考文献

摘要：
1 Processing Map的发展▲
1.1 早期的加工图
1.2 基于动态材料模型的加工图
1.3 复杂本构方程材料加工图构建的数学基础
2 基于Matlab平台的数值计算过程▲
2.1 试验数据的获得及处理
2.2 Matlab下热加工图的数值计算步骤
2.3 Matlab平台的具体计算程序
3 结论▲
图表▲

表1不同温度、应变、应变速率下新锆合金的应力值 (MPa) Table 1Flow stress (in MPa) of new Zr alloy at various strains, strain rates and temperatures

图1 应变水平为0.5时的N18新锆合金加工图 (阴影为不稳定流变区) Fig.1 Processing map for N18 Zr alloy at a strain of 0.5. Shadou region corresponds to flow instability

基于Matlab的热加工图的数值构造方法

李中奎张建军田锋

西北工业大学材料学院,西北有色金属研究院难熔金属研究所,西北有色金属研究院难熔金属研究所,西北有色金属研究院难熔金属研究所陕西西安710072,西北有色金属研究院难熔金属研究所,陕西西安710016,陕西西安710016,陕西西安710016,陕西西安710016

摘要：

利用热加工图 (Processing Map) 可以分析材料在不同区域 (不同的变形温度和应变速率) 的高温变形机制, 从而获得热加工“安全区”和“不安全区”, 达到控制组织演变、避免缺陷的产生和优化工艺参数的目的。Murty提出热加工图计算方法适用于复杂本构方程的材料, 但其公式中含有积分项, 计算过程复杂。Matlab具有强大的矩阵计算功能, 本文给出了基于Matlab平台的计算过程和简易程序, 可快捷得出计算结果, 并绘制出材料的热加工图。

关键词：

热加工图;Matlab;数值法;

中图分类号： TG301

作者简介：李中奎, 通讯联系人 (E-mail:Lizhangyi@c-nin.com) ;

收稿日期：2007-05-08

基金：国家自然科学基金资助 (50571083);

A Numerical Computation Method for Hot Processing Map Based on Matlab

Abstract：

Processing map is an effective tool for optimization of materials workability, product property control and defect avoidance.Processing map is constructed by a superimposition of efficiency of power dissipation and the instability maps, and can be pided into safe domain and unsafe domain.The calculating theory for the hot processing map proposed by Murty is suitable for complex constitutive equation, but because of the integral in the computing formula, it is difficult to apply without computer program aid.A new numerical computation method and program for Murty hot processing map based on Matlab was presented, by using this method, the processing map can be easily figure out.

Keyword：

hot processing map;matlab;numerical computation method;

Received： 2007-05-08

金属塑性成形是指金属材料在一定的外力作用下, 利用其塑性发生塑性变形并获得一定力学性能的加工方法。热加工图主要目的是通过对塑性成形过程的分析与模拟, 掌握各种工艺参数对成形过程金属流动规律的影响, 以避免缺陷的产生, 达到控制组织和性能的目的。利用热加工图 (Processing Map) 可以分析材料在不同区域 (不同的变形温度和应变速率) 的高温变形机制, 从而获得热加工“安全区”和“不安全区”, 达到控制组织演变、避免缺陷的产生和优化工艺参数的目的。

1 Processing Map的发展

1.1 早期的加工图

Frost和Ashby首先以变形机制图的形式研究了材料的变形特征, 其重点主要放在低应变速率低蠕变机制应用研究方面 ^[1] 。

Raj发展了Ashby的思想构建了加工图, 从Raj的加工图中可得出硬质点在低温高应变速率下于软基体上形成的空洞区域以及在高温低应变速率下三角晶界处低楔形开裂产生区。 Raj研究构建了一些简单合金的加工图, 但对于复杂合金, Raj加工图还是存在局限性 ^[2] 。

Semiatin 和Lahoti将流变软化和材料性质相关联, 提出了归一化流变软化速率、应变速率敏感因子和局部流变判据α^*=-γ/m, 根据显微组织观察, Semiatin等 ^[3] 认为钛合金中α^*>5则会发生局部流变现象的经验值, 但是, 这完全根据实际经验取值, 没有严密的理论依据, 故参数α^*的取值范围会因材料不同而异。

1.2 基于动态材料模型的加工图

Prasad提出用动态材料模型 (DMM) 描述材料的变形行为 ^[4] 。根据动态材料学模型的观点, 材料的加工过程可看作为一个能量耗散系统, 能量的消耗取决于材料的加工流变行为, 它服从幂律方程:

式中σ为流变应力, 为应变速率, m为应变速率敏感指数, K为常数。在应变速率一定的情况下, 加工过程中单位体积内材料所获得的能量P可分为两部分:

即外部输入总能量P分为两个部分: G是材料发生塑性变形所消耗的能量, 其中大部分能量转化成了热能, 小部分以晶体缺陷能的形式存储; J是材料变形过程中组织演化所耗的能量。这两种能量所占比例由材料在一定应力下的应变速率敏感指数m决定, 即:

并提出耗散效率系数η, 其物理意义为材料成形过程中显微组织演变所耗散的能量同线性耗散能量的比例关系, 其值为:

通常可绘出由和T所构成的平面上的能量耗散因子 (η) 等值线图, 从图中可直接分析不同区域的变形机制, 如动态再结晶、动态回复和超塑变形等加工安全区, 以及楔形开裂、空洞形成等缺陷区。

根据大应变塑性变形的极大值原理, 得到流变失稳的判据为:

1.3 复杂本构方程材料加工图构建的数学基础

Prasad的热加工图是基于公式 (1) 的, 即认为应变速率敏感指数m不随和T变化, 为一定值, Murty等 ^[5] 认为, 对于大多数工程合金, m值应该随和T变化, 流变应力是不符合幂率方程, 因此利用公式 (4) 计算得出η并不精确。 Murty提出了任意类型曲线的流变失稳准则。

在公式 (2) 中 , 是一个从0到σ的对应力所进行的积分, 但在实际试验中, 流动应力在变形过程会发生振荡, 因此应力值不是一个单调函数, 积分计算的过程中比较困难。如果能将流变应力σ表示为ε, 和T的函数, 由于应变的增加是单调的, 积分就能够容易精确的计算出。

由

即 (6)

G的积分转化为 (7)

式中 _min取一般试验中的最小应变速率值。

耗散效率

2 基于Matlab平台的数值计算过程

Murty提出的加工图的计算方法是基于任意的本构方程的, 应具有更普遍的适用性, 但是涉及到对应力或应变加载过程的积分, 计算过程较为复杂。一般可以使用两种方法计算, 一种为解析法。根据试验数据, 选择合适的本构方程, 将合金应力——应变关系用该本构方程表示出来, 然后依据公式 (7) ~ (8) 进行计算得出耗散效率η和不稳定判据ζ的解析解, 然后作出热加工图。但这种方法需要大量人工计算过程和良好的数学功底, 对计算人的要求较高, 其误差来源主要为本构方程的选择与构建过程所产生的偏差。

第二种计算方法为数值计算法, 该方法最大的优点在于不需要知道试验材料本构方程复杂的数学表达式, 直接对试验数据进行必要的拟合、运算, 最终得到数值解, 然后作图。数值计算法的误差来源于对试验数据的拟合过程, 但如果选择较好的拟合方法, 如三次样条法, 可以将这种误差尽量降低。

数值计算法需要运算大量的高阶矩阵, 如果用人工计算是几乎无法完成的任务, 但是借助于Matlab程序的强大的矩阵运算能力可轻松的完成。本文将对Matlab的计算过程作一介绍。

2.1 试验数据的获得及处理

为获得构建加工图的数据, 需进行热模拟机试验。试验的温度和应变速率范围的选择要根据需优化的加工工艺的范围而制定。为了便于理解, 本文以新锆合金N18 ^[6] 为例进行说明。试验温度为600, 660, 730, 775, 860, 980 ℃, 应变速率为0.001, 0.1 1, 10, 50 s^-1, 试验后的数据根据需要, 应转化为真应力-真应变曲线, 可得出同一应变水平下不同温度、不同应变速率时的流变应力值, 如表1。

2.2 Matlab下热加工图的数值计算步骤

(1) 试验数据的插值由于试验过程的温度和应变速率只能选择有限个点, 一般只能形成4

表1不同温度、应变、应变速率下新锆合金的应力值 (MPa) Table 1Flow stress (in MPa) of new Zr alloy at various strains, strain rates and temperatures

Strain	Strain rate/s^-1	Temperature/℃
Strain	Strain rate/s^-1	600	730	775	860	980
0.2	0.001	113.8	63.3	42.0	25.2	12.8
	0.01	155.5	78.1	55.1	42.7	13.6
	0.1	177.6	121.4	87.8	55.1	19.4
	1	218.1	164.7	108.6	73.6	38.9
	10	252.4	214.3	190.9	111.8	65.4
	50	270.8	234.4	209.9	121.9	95.4
0.4	0.001	111.4	59.1	38.0	24.3	11.6
	0.01	155.5	72.8	49.7	37.9	12.0
	0.1	183.9	116.1	80.5	53.8	16.5
	1	238.9	156.3	108.1	69.3	36.3
	10	274.5	225.9	184.6	113.5	64.3
	50	287.1	251.3	215.4	117.7	93.9
0.5	0.001	108.3	57.1	36.1	24.3	10.9
	0.01	152.4	71.7	48.8	35.4	11.4
	0.1	184.9	112.9	77.8	52.6	15.2
	1	239.4	152.0	107.6	67.6	34.7
	10	269.8	221.7	178.3	110.1	61.8
	50	280.5	247.1	211.8	112.6	89.7

价或5阶矩阵, 无法进行运算, 必须通过插值计算获得中间温度、应变速率的应力值。选择插值的方法有多种, 可以为线性、多项式、最近插值、样条插值等, 采用三次多项式或三次样条插值法可以较好满足试验精度, 在Matlab中可用“interp2”命令实现插值计算。其基本命令格式为:

ZI=interp2 (X, Y, Z, XI, YI, method)

(其含义为: 用“method”指定的算法进行插值计算, 返回矩阵ZI, 其元素包含对应于参量XI与YI (可以是向量、或同型矩阵) 的元素, 即Zi (i, j) [Xi (i, j) , Yi (i, j) ]。用户可以输入行向量和列向量Xi与Yi, 此时, 输出向量Zi与矩阵meshgrid (Xi, Yi) 是同型的。同时取决于由输入矩阵X、 Y与Z确定的二维函数Z=f (X, Y) 。参量X与Y必须是单调的, 且相同的划分格式。

Method算法包含有:

“linear”: 双线性插值算法 (缺省算法) ;

“nearest”: 最临近插值;

“spline”: 三次样条插值;

“cubic”: 双三次插值。

通过插值计算, 得到50×50的二维矩阵就可足够满足计算精度需求了。

(2) 根据公式 (3) 计算每个试验数据点的m值。可利用Matlab中的“gradient” —— 梯度函数实现。

(3) 根据公式 (7) , (8) 计算耗散效率η值, 其难点在于的积分的运算。可利用“cumtrapz”函数——累计梯形积分函数实现。其基本命令格式为: Z=cumtrapz (X, Y) (其含义为用梯形法计算Y在X点上点积分)

(4) 根据公式 (5) 计算不稳定判据ζ。

(5) 由所计算出的η值、 ζ值矩阵, 利用绘图软件 (如Origin) , 在变形温度T和应变速率所构成的二维平面上绘制出η的等值线图-热加工图。图1为所绘制的新锆合金的热加工图。

图1 应变水平为0.5时的N18新锆合金加工图 (阴影为不稳定流变区) Fig.1 Processing map for N18 Zr alloy at a strain of 0.5. Shadou region corresponds to flow instability

2.3 Matlab平台的具体计算程序

X=[600, 660, 730, 775, 860, 980]% X轴, 温度

Y=[-3, -2, -1, 0, 1, 1.70]%Y轴 log应变速率

Z=[108.3, 78.5, 57.1, 36.1, 24.3, 10.9, 152.4, 112.9, 71.7, 48.8, 35.4, 11.4, 184.9, 154.4, 112.9, 77.8, 52.6, 15.2, 239.4, 215.4, 152.0, 107.6, 67.6, 34.7, 269.8, 249.1, 221.7, 178.3, 110.1, 61.8, 280.5, 264.1, 247.1, 211.8, 112.6, 89.7]

Zlog=log10 (Z) % log应力

i=50%设置插值矩阵的行数

j=50%设置插值矩阵的列数

Xi= (Ts: (Tf-Ts) / (j-1) :Tf) %赋值给X插值点

Ylogi= (-1: (1.70- (-1) ) / (i-1) :1.70) %赋值给Y插值点

Yi=10.^Ylogi %应变速率的插值向量

Zlogi=interp2 (X, Y, Zlog, Xi, Ylogi′, ‘spline’) % 用spline方法插值 (log应力) 矩阵

Zi=interp2 (X, Y, Z, Xi, Ylogi′, ‘spline’) % 用spline方法插值应力矩阵

[Fx, M]=gradient (Zlogi, 1, (1.70- (-1) ) / (i-1) ) % 计算log力/log速率=m%以下为通过递归循环计算积分

S=cumtrapz (Yi, Zi) %计算累积梯形积分

P=zeros (i, j) %设置全0矩阵

Es=-1%赋值应变速率最小值