２００６年ＡＮＳＹＳ年会论文集

yhlhzg

免费下载

yhlhzg

2006 年用户年会论文
钢板轧后冷却过程的数值模拟及变形分析
崔青玲 袁国 刘相华 王国栋
[东北大学轧制技术及连轧自动化国家重点实验室 沈阳 110004]
[ 摘 要 ]：利用热弹塑性有限元法对中厚板轧后冷却过程的温度场及应力应变场进行数值模拟。采用分
步循环加载的方法精确模拟不同换热边界条件，温度场模拟结果与实测值吻合较好。以热分析
为基础，将热分析得到的温度场作为载荷进行钢板冷却过程应力应变的模拟，讨论了钢板产生
横向变形的主要原因，为获得平直板形提供了理论基础。
[ 关键词 ]：轧后冷却，数值模拟，热弹塑性，冷却变形
中图分类号: TG 335.1 文献标识码: A
Numerical simulation of steel plate runout table
cooling and its deformation analysis
CUI Qing-ling, YUAN Guo, LIU Xiang-hua, WANG Guo-dong
[The state key laboratory of rolling and automation, northeastern university, shenyang110004]
[Abstract]: The temperature and stress-strain field of steel plate runout table cooling is simulated using
thermo-elasto-plastic finite element method. Different heat transfer boundary condition was
exactly simulated by stepped-circularly loading method. And the simulation results of
temperature keep good accordance with the experimental one. Based on the thermal analysis,
the stress-strain field is calculated loading with the temperature from result of thermal analysis.
The main reason of transverse deformation is discussed，and it can be seen as the theoretical
bases of getting flat plate.
[Key words]: runout table cooling, numerical simulation, thermo-elasto-plastic, cooling deformation
中厚板轧后控制冷却的工艺要求，除获得所需的组织性能外，平直的板形是冷却设备实用
化的关键。钢板冷却过程中引起翘曲的内在因素是相变、屈服应力和导热系数，其诱导因素则
是不均匀冷却。影响冷却均匀的原因很多，在采用同一种冷却方式的前提下，冷却区长度、上
国家自然科学基金资助项目：50474016
作者简介：崔青玲,1973,女，博士，讲师.Tel: 024-83681811, Email: Cuiql@mail.neu.edu.cn.
2006 年用户年会论文
下水管水量比、纵向集管分布间距和横向射流间距等都对钢板的变形产生不同程度的影响[1]。
本文采用数值模拟方法对钢板冷却过程的温度场及应力应变场进行模拟，分析钢板横向变形的
主要原因，为中厚板冷却获得平直板形提供理论基础。
1 热弹塑性有限元的基本原理
钢板冷却过程不存在外载荷，引起热应力和热变形的原因是温度分布不均匀，各点的膨胀
量不同，因此钢板控冷过程的应力应变属于热弹塑性问题。在塑性区域内，由于屈服条件依赖
于温度，则等向强化Mises 屈服准则为[2]：
− H(∫d ,T) = 0 p σ ε （1）
其微分形式为
{ } dT
T
d H d H p ∂
∂ ′
= ′ +
⎭ ⎬ ⎫
⎩ ⎨ ⎧
∂
∂
σ ε
σ
σ Τ
（2）
设温度增量为dT，则全应变增量可分解为
{ } { } { } { } { } 0 d ε d ε d ε d ε d ε e p T = + + + （3）
式中， { } e d ε 为弹性应变， { } p d ε 为塑性应变， { } T d ε 为温度应变， { } 0 d ε 为考虑温度对弹
性模量影响时的附加应变。
考虑基于Mises 屈服准则基础上的流动法则，则应力增量可写做
{ } [ ] { } { } { } { }⎟ ⎟⎠
⎞
⎜ ⎜⎝
⎛
− −
∂
∂
= − 0 ε ε ε
σ
σ
d σ D d ε d d d e p T （4）
将式（4）代入式（2）可得
[ ]( { } { } { })
[ ] {σ }
σ
σ
σ
ε ε ε
σ
σ
ε
∂
∂
⎭ ⎬ ⎫
⎩ ⎨ ⎧
∂
∂
′ +
∂
∂ ′
− − −
⎭ ⎬ ⎫
⎩ ⎨ ⎧ ∂
∂
=
e
e T
p
H D
dT
T
D d d d H
d Τ
0
Τ
（5）
将式（5）代入式（4）并令
2006 年用户年会论文
[ ]
[ ] { } { } [ ]
[ ] {σ }
σ
σ
σ
σ
σ
σ
σ
∂
∂
⎭ ⎬ ⎫
⎩ ⎨ ⎧
∂
∂
′ +
⎭ ⎬ ⎫
⎩ ⎨ ⎧
∂
∂
∂
∂
=
e
e
T
e
p
H D
D D
D Τ （6）
[ ] [ ] [ ] ep e p D = D − D （7）
则可得
{ } [ ]( { } { } { }) { } 0 0 d σ D d ε d ε d ε d σ ep T = − − + （8）
式中， { }
[ ] { }
[ ] {σ }
σ
σ
σ
σ
σ
σ
∂
∂
⎭ ⎬ ⎫
⎩ ⎨ ⎧
∂
∂
′ +
∂
∂ ′
∂
∂
=
e
e
H D
dT
T
D H
d 0 Τ 是由温度对塑性模量影响而引起的附加应力。
2 中厚板冷却的数值模拟
2.1 某厂钢板冷却装置组成简介
某中厚板厂采用的高密直集管层流冷却装置如图1 所示[3]，该控冷装置由上下集管、高压
水侧喷装置、压缩空气侧吹装置等基本设备组成，冷却区入口、出口处分别装有红外测温仪，
可以将开冷温度、终冷温度和返红温度在线输入计算机系统。针对钢板的厚度和所获得的组织
性能不同而采用不同的冷却策略（集管的开启组数、喷水流量和钢板速度）。
图1 层流冷却装置结构示意图
Fig.1 Scheme of the laminar cooling equipment
2.2 数值模拟的建模过程
采用热应力分析的间接耦合法对钢板冷却过程应力应变场进行模拟，即在热分析的基础上，
将求得的节点温度作为体载荷施加在结构分析中。
2006 年用户年会论文
由于钢板的长度远大于钢板的厚度和宽度，且钢板运行速度较快，可忽略钢板长度方向的
温度变化，将问题简化为二维问题。以钢板长度方向任一截面为研究对象，钢板上下表面受冷
却水射流冲击发生换热，表面温度梯度比心部大，因此表面网格划分比心部密，如图2 所示。
图2 钢板有限元网格（沿宽度截取部分）
Fig.2 mesh of the steel plate(partition along width)
模拟计算从钢板进入水冷区开始。为使数值模拟真实反映控冷过程，通过对钢板表面对流
换热边界条件的分析，对不同的换热方式定义不同的换热系数，并根据钢板通过冷却区的速度
分别计算不同冷却方式的作用时间，采用分步循环加载的方法加载。
2.3 换热系数的确定
在冷却过程中冷却水与热带钢表面之间存在着复杂的局部热交换过程，如图3 所示。根据
热交换方式，对带钢上表面划分为两个不同的局部冷却区域，即流冲击换热区和膜换热区；对
于带钢下表面，由于受重力影响，不存在稳定的汽膜换热区域，因此仅考虑射流冲击区和辐射
换热区。
此外，带钢还受到空气对流换热以及与输出辊道接触产生的热传导散热，这些换热方式对
带钢整体换热影响不大，可将其综合考虑到辐射换热过程中。
（a）上表面 （b）下表面
图3 钢板冷却过程局部换热方式示意图
Fig.3 The local heat transfer mode description of steel plate during cooling
2006 年用户年会论文
（1） 射流冲击区对流换热系数
冷却水冲击钢板，在约为集管直径1.5 倍大小的范围内形成射流冲击区，其换热方式为核
沸腾换热，根据实验研究[4]确定其对流换热系数为：
0.063 Re0.8 Pr1 3
r
w λ
α = （12）
式中,α 为对流换热系数，Re 为雷诺数；Pr 为普朗特数；λw 为冷却水的热传导率；r 为距冲
击点的距离。
（2）汽膜换热区
水流落在带钢上表面，冷却水沿冲击区径向流动，形成壁面射流区，其换热方式为膜沸腾
换热。Jocobi 等人[5]的实验研究认为，一旦带钢表面建立起稳定的汽膜沸腾换热，带钢表面的换
热系数将显著减小。在带钢温度大于350~400℃时，汽膜沸腾换热系数与带钢表面温度无关，
其换热系数范围[5,6]约为150~ 350 Wm-2K-1。
（3）辐射换热区
根据Stefan-Boltzmann 定律，不同温度的带钢发射率约在0.6 到0.95 之间。辐射换热区的
换热系数较小，一般为5~30 Wm-2K-1 之间[7]。
3 模拟结果及分析
（1）温度场模拟结果
图4 为规格为2.4m×0.02m×15m 的Q235B 钢板，集管水量为100(大)/80(小)m3/h，上下集
管水量比2.65(大)/2.15(小)，集管开启数为2(大)/ 7(小)组，边部遮蔽量在55mm 左右，钢板运行
速度为1m/s 的工况下，钢板刚出水冷区时即冷却时间t=20s 后的温度分布情况。可以看出在采
用边部遮蔽后，宽向温度基本均匀，边部温度过低现象得到改善。图5 为t=20s 时计算值和实
测值[3]的比较，可以看出实测温度和计算值相差为20℃左右，两者吻合较好。
2006 年用户年会论文
图4 冷却20s 后钢板温度分布图（截取边部部分）
Fig.4 temperature distribution of plate after cooling 20s (partition of width near side)
0 200 400 600 800 1000 1200
650
660
670
680
690
700
温度(°C)
板宽(mm)
实验值
计算值
图5 钢板宽度方向温度计算值和实测值比较
Fig.5 Comparison of calculated temperature and experimental ones along plate width
图6 为钢板厚度方向的温度时间历程曲线。可以看出，在水冷阶段，钢板上下表面的温度
变化曲线呈锯齿形，而心部温度的降低是一个连续平稳的过程。这是由于钢板表面受冲击换热
区、汽膜换热区和辐射换热区三种换热载荷的周期性作用，存在温度的不断降低和回升，而心
部温度变化是钢板内部导热的结果，变化相对平缓。水冷阶段结束后，钢板表面和心部的温度
趋于一致，正是钢板返红的过程。
从此温度时间历程曲线可以很清楚地看出钢板冷却的温度变化过程。曲线上的一个个温度
速降区域（锯齿的下降边），是钢板上下表面受冲击区换热区作用的结果，也就是说每一个温
度速降区域对应钢板长度方向的一组集管。
2006 年用户年会论文
0 5 10 15 20 25 30 35 40
600
650
700
750
800
温度(°C)
冷却时间(s)
上表面
中心
下表面
图6 钢板冷却过程温降曲线
Fig.6 Temperature drop curve of plate during cooling
（2）应力应变场结果分析
图7 和图8 分别为t=2.04s 时刻钢板的y 向位移和钢板厚度方向三个截面的等效应力沿板宽
方向的分布。钢板冷却过程中的热应力是由截面温度差造成的，冷却速度越大，截面温差越大，
产生的热应力越大。比较图7 和图8 看出，冷速最大的上表面等效应力最大，冷速最小的中心
层等效应力最小。由于钢板上表面冷速大于下表面冷速，上表面的热应力大于下表面，则上表
面比下表面收缩得快，因此造成钢板向快冷得上表面横向翘曲变形。
图7 钢板冷却2.04s 的变形
Fig.7 Deformation of plate at cooling time 2.04s
0.0 0.2 0.4 0.6 0.8 1.0 1.2
0
20
40
60
80
100
120
Mises等效应力(MPa)
板宽(m)
上表面
中心
下表面
图8 冷却2.04s 板宽方向等效应力分布
Fig.8 Distribution of equivalent stress along plate width at cooling time 2.04s
2006 年用户年会论文
图9 为上下表面和对应心部的应力和应变随时间变化的计算结果。因为残余翘曲的大小由
钢板冷却前的温度分布和冷却过程中板厚方向的塑性应变决定，假设冷却前钢板温度分布均匀，
则塑性应变决定残余翘曲的大小。由图9 可以看出钢板上下表面的塑性应变差出现在进入水冷
的开始阶段，即冷却初期就决定了钢板的残余翘曲。因此，为了获得平直的板形，必须调整上
下水量比，尽量减小水冷开始时上下表面的应力和塑性应变差。
0 10 20 30 40
0
30
60
90
120
150
180
冷却时间 (s)
Mises等效应力(MPa)
上表面
中心
下表面
(a)
0 10 20 30 40
-0.4
0.0
0.4
0.8
1.2
1.6
2.0
总应变(×10-3)
冷却时间(s)
上表面
中心
下表面
(b)
0 10 20 30 40
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8 塑
性
应
变
(
×10-
3
)
冷却时间(s)
上表面
中心
下表面
(c)
2006 年用户年会论文
图9 钢板应力和应变随冷却时间的变化
(a)等效应力 (b)等效全应变 (c)等效塑性应变
Fig.9 Stress and strain variation during cooling
(a) equivalent stress (b) total equivalent strain (c)plastic equivalent strain
4 结论
利用热弹塑性有限元方法对钢板轧后冷却过程进行数值模拟，温度场和应力应变场模拟结
果表明在上表面冷却速率大于下表面时，上表面的内应力大于下表面，则上表面比下表面收缩
得快，造成钢板向冷速较快的上表面翘曲。塑性应变差出现在水冷的开始阶段，表明冷却初期
就决定了钢板的残余翘曲。因此，为了获得平直的板形，必须调整上下水量比，尽量减小水冷
开始时上下表面的应力和塑性应变差。
参考文献
[1] 王峰丽,李谋渭,张少军等.中厚板控冷过程的数值模拟与冷却变形分析[J]. 机械工程学报，2003，
39(5):111-115.
Wang fengli, Li mouwei, Zhang shaojun, et,al. Study on numerical simulation of controlled cooling
medium steel plate and analysis on its deformation[J]. Journal of mechanical engineering, 2003，
39(5):111-115.
[2] 刘庄,吴肇基,吴景之等.热处理过程的数值模拟[M].北京：科学出版社, 1996.
Liu zhuang, Wu zhaoji, Wu jingzhi,et al. Numerical simulation of heat treatment processing[M].
Beijing: Science Press, 1996
[3] 于明.中厚板控制冷却过程数学模型及控制策略研究[D].东北大学，2004，12.
Yu ming. Study on cooling mathematical model and controlled strategy for medium and heavy plate
[D]. Northeastern University,2004,12.
[4] Natsuo Hatta, Jun-ichi Kokado, Hirohiko Takuda, et al. Predictable modeling foe cooling process
of a hot steel plate by a laminar water bar[J]. Arch. Eisenhuttenwes, 1984, 55(4): 143-148
[5] H.Jacobi,G.Kaestle et al.Heat Transfer in Cyclic Secondary Cooling during Solidification of
Steel[J]. Ironmaking and Steelmaking,1984，11（3）：132-144.
2006 年用户年会论文
[6] M.M.Prieto,L,S.Ruíz et al.Thermal performance of Numerical Model of Hot Strip Mill Runout
Table[J]. Ironmaking and Steelmaking,2001，28（6）：474-480.
[7] S.D.Cox,S.J.Hardy et al.Influence of Runout Table Operation Setup on Hot Strip Quality,Subject
to Initial Strip Condition:Heat Transfer Issues.Ironmaking and Steelmaking, 2001，28（5）：
363-372.

yhlhzg

2006 年用户年会论文
激光快速成形激光、粉末与熔池交互数值模拟
贾文鹏1,2，陈静2，林鑫2，钟诚文1，黄卫东2
1.西北工业大学高性能计算研究与发展中心 710065
2. 西北工业大学凝固技术国家重点实验室 710065
[ 摘 要 ] 本文针对激光快速成形过程中激光、粉末与熔池交互，建立了气/粉两相流场、激光快速成形
件及熔池温度场模型。采用ANSYS 单元生死技术模拟了激光加热及熔覆层的生长及熔池及自
由界面形状，采用CFX 基于拉格朗日方法的粒子跟踪（particle tracking model）模型追踪粉
末颗粒并根据界面温度设置粒子平性（parallel）和垂直(perpendicular)动量损失模拟熔池对粒
子的捕捉以及工件对粉末的反射。计算了316L 不锈钢粉末、激光与熔池的交互，计算结果与
实验相吻合。
[ 关键词 ] 激光快速成形, 气/粉两相流场, 熔池，交互
The numerical simulation of interaction between melting
pool and metal powder during laser rapid forming
W.P.Jia1,2, Jing Chen2, Lin xin2, Zhong wencheng1,weidong Huang2
1.Northwestern Polytechnical University, Center of High Performance Computing
Xi’an, China, 710072, E-mail:hpc@nwpu.edu.cn
2. Northwestern Polytechnical University, State Key Laboratory of Solidification
Processing Xi’an, China, 710072,E-mail:cgs@nwpu.edu.cn
[ Abstract ] The Laser Rapid Forming (LRF), an advanced solid freedom fabrication technology, has a
bright future in manufacturing of high performance dense metal components with complex structure in
aero, spaceflight and many other fields. During LRF, the melting pool free surface serves as a moving
window of laser power and powder mass input and a dynamic boundary of the LRF continuous growth
domain, so the interaction between laser, powder and melting pool free surface is the basic problem that
should be coped with in LRF process. In this paper, the melting pool transient temperature field and
gas/powder two phase fluid field numerical model are established. Using ANSYS`s finite element birth
and death technology, the growth of cladding layers and shape of melting pool free surface are simulated,
and using CFX`s particle tracking model and momentum loss (parallel and perpendicular)option, the
interaction between melting pool free surface and powder such as capture and reflect are simulated. At
last, the interaction between 316L stainless steel powder, laser and melting pool free surface is
calculated, and the numerical result is in agreement with that of the experiment.
[ Keyword ]. Laser Rapid Forming, gas/powder two phase fluid, melting pool, interaction
2006 年用户年会论文
1 前言
激光快速成形技术是一项新的先进的制造技术,能够实现高性能复杂结构致密金属零
件的快速、无模具、近终形制造,在航空、航天、汽车等高技术领域具有光明的应用前景
[1-2] 。随着激光快速成形技术研究的深入开展,迫切要求发展建立能够准确描述激光熔覆
过程的理论模型以准确把握其内在机理。
激光快速成形的核心是激光熔覆——激光熔化粉末并逐层堆积的过程（如图1），在
此过程中熔池自由表面是激光能量和粉末质量进入熔池的自由界面，同时也是熔覆层生长
的动态边界，所以粉末与熔池交互是激光快速成形过程不可回避的基本问题，而要实现高
性能复杂结构致密金属零件的整体精确制造则必须建立可靠的激光、粉末与熔池交互，从
而在此基础之上实现激光快速成形过程的模拟。
2 激光、粉末与熔池交互
在激光快速成形过程中，由于激光加热和粉末流的输入，导致熔覆层的生长，而熔覆
层的不断生长及逐层堆积又导致了成形件连续体域的不断变化，而成形件连续体域的这种
变化又会影响动态温度场、流场的形成及发展，反过来，动态温度场、流场又会影响熔池
对激光能量和粉末吸收，导致成形件连续体域的不断变化。上述过程如图2 所示。
图1 激光熔覆粉末熔化与逐层堆积过程
2006 年用户年会论文
图2 激光、粉末与熔池交互过程
由图可以看出，熔覆层生长是粉末与熔池自由表面作用、粉末与熔池熔液混合流动及
熔化最后凝固堆积的结果，粉末之所以能进入熔池是因为激光加热成形件表面形成温度场
分布，而在熔池处的温度超过熔点，形成了粉末可进入熔池的窗口——液态的熔池自由表
面。熔池自由表面是成形件熔覆层生长的动态边界，是能量和质量进入的自由界面，所以
激光、粉末与熔池交互伴随着熔池自由表面界面移动及熔覆层生长，是和熔池温度、流场
及气/粉两相流场的形成及发展密不可分的。总的来说，激光、粉末与熔池交互主要包括：
a. 激光与粉末交互，包括粉末对激光的遮蔽、粉末温度的升高及激光能量的空间分布；b. 气
粉两相流的形成及粉末在工件及熔池表面的空间质量分布; c. 粉末与熔池表面以及工件表
面交互，包括粉末的被捕捉及反弹；d. 熔池自由表面界面移动及熔覆层生长，包括温度及
流场及自由表面界面演化等几大问题。而以上问题又相互有机的紧密联系着，期望采用统
一模型同时解决难度相当大。
2.1 基本假设及方法步骤：
鉴于激光、粉末与熔池交互问题的复杂性，本文在其他学者实验及理论计算的基础上
作如下假设：
1) 激光能量在光束横截面上服从高斯分布，粉末对激光的遮蔽是粉末体积分数的线性函
数[3]；
2) 粉末在气体中的体积分数很低，只考虑气体流动和重力对粉末的影响，不考虑粉末颗
粒之间的作用[4]；
2006 年用户年会论文
3) 固体颗粒撞到固体表面就会反弹而撞入液态熔池就会被吸收[5]；
4) 粉末进入熔池后在熔化前已经完全混合或已经完全分布在液态金属中[6][7]；
5) 对于液态金属对流换热假设材料高于熔点的热传导系数是平均的2.5 倍[8]；
这样主要考虑气粉两相流场、粉末在工件及熔池表面的空间质量分布、熔池温度场及
熔化自由表面形状、粉末的被捕捉及反弹和熔覆层生长。采用ANSYS 单元生死技术[9]模
拟激光加热及熔覆层的生长及熔池及自由界面形状，采用CFX 基于拉格朗日方法的粒子跟
踪（particle tracking model）模型追踪粉末颗粒，具体方法步骤如下：
1) 在ANSYS 中建立喷嘴及工件几何模型并划分网格；
2) 将ANSYS 几何模型导入CFX 设置气/粉两相流场与工件表面交互规则，根据界面温度
设置粒子平性（parallel）和垂直(perpendicular)动量损失模拟熔池对粒子的捕捉以及工
件对粉末的反射，并由CFX 计算喷嘴内外气/粉两相流场分布并与工件表面耦合得到
粉斑的大小、粉末固相分数的分布以及在气/粉两相流冷却下的工件表面温度；
3) 将由CFX 计算得到粉斑分布、得粉率以及工件表面温度导入ANSYS，与激光能量分
布耦合计算激光快速成形温度场，设定生长规则：如果熔化表面相邻单元得粉大于单
元体积则由ANSYS 单元生死激活生长单元，并对生长界面进行界面重构得到新的成
形件几何模型；
4) 将由ANSYS 界面重构得到新的成形件几何模型经网格划分后连同激光快速成形温度
场导入CFX 设置气/粉两相流场与工件及熔池表面交互规则，根据界面温度设置粒子
平性（parallel）和垂直(perpendicular)动量损失模拟熔池对粒子的捕捉以及工件对粉末
的反射，由CFX 计算生长后模型的粉斑粉末固相分数的分布、熔池得粉率以及在气/
粉两相流冷却下的工件及熔池表面温度；
5) 随着激光的移动重复3 和4 步直到激光成形进入准稳态；
6) 采用准稳态CFX 计算得到粉斑分布、得粉率结果在ANSYS 中计算激光成形温度场直
到结束。
2.2 计算及结果分析：
2.2.1 物理模型：
在LRF855 激光快速成形机上采用316L 不锈钢粉末在30 mm×10 mm×7 mm（长×宽
×高）45 钢块状试样上沿宽度中心线从一端向另一端进行激光快速成形，成形参数见表1，
表1 激光快速成形工艺参数
为简化计算，取对称1/2 模型，激光垂直照射工件表面，喷嘴与工件倾角Φ=46°，其
轴线与激光光轴相交于工件表面，如图3：
激光功率（W） 光斑直径（mm） 扫描速度（mm/s） 送粉率（g/min）
2300~3400 3~4 5~10 15~20
2006 年用户年会论文
图3 计算模型
整个计算域包括气粉两相外流场100mm×25 mm×50 mm（长×宽×高）及成形试样
固体两个domain,它们之间通过interface 自动连接，外流场边界条件见图4 a，其中入口气
体流速0.01m/s, 出口相对平均静压力0 Pa,开口（opening 边界）相对开口压力0 Pa，约束
气体流速12m/s, 粉末为球形粉末，最小直径0.01~0.1mm，按质量正态分布，初始速度
0.01m/s；整个计算域采用非结构四面体网格划分如图4b：
图4 网格模型及边界条件
激光快速成形温度场及熔覆层生长在ANSYS 中采用单元生死模拟，整个计算域采用
结构六面体Solid70 单元网格划分如图5 a，首先杀死所有成形单元，并加载换热边界条件
如图5 b：
a 边界条件 b 网格模型
2006 年用户年会论文
a ANSYS Solid70 单元网格模型 b 杀死所有成形单元并加载换热边界条件
图5 ANSYS Solid70 单元网格模型及边界条件
2.2.2 计算结果：
开始时工件表面未熔化，气/粉两相流场及工件表面上平均体积分数分布如图6：
a 气/粉两相流场 b 工件表面上平均体积分数分布
图6 未熔化表面气/粉两相流场及平均体积分数分布
由图可见：首先粉末随载粉气体到达喷嘴并在重力作用下从出粉口下落，在重力及约束气
流的带动下加速下落形成气/粉两相流，并在工件表面上形成近似椭圆质量分布，由于工件
表面未熔化，此时粉末颗粒全部反弹，没有熔覆层形成。
随激光的加热及温度场的演化，熔化表面形成，熔池开始捕捉粉末颗粒，混合熔化后
形成新的熔覆表面形状，如图所示7：
2006 年用户年会论文
a 单元生死模拟熔覆表面 b 界面重构
图7 粉末熔覆表面形状及界面重构
基于固体颗粒撞到固体表面就会反弹而撞入液态熔池就会被吸收的假设和该处粉末的体积
分数分布，可得在时间步内是否激活单元使其成为熔覆层，图7 a 为单元生死模拟熔覆表
面，而图7 b 为通过界面重构及重新温度插值，将ANSYS 模型转换到CFX 中。
在CFX 中计算气粉两相流外流场，熔覆层表面气/粉两相流场及平均体积分数分布如
图8 所示：
a 气/粉两相流场 b 熔覆层表面上方平均体积分数分布
图8 熔覆层表面气/粉两相流场及平均体积分数分布
比较图8 a, b 和图6 a, b 可以看出：形成熔覆层后由于熔化自由表面对粉末的捕捉反弹颗粒
大为减少，最大粉末平均体积份数分布也集中在熔覆层表面前沿小斜面。增加追踪粒子样
本数至200，粉末进入熔池和反弹如图9 所示：
2006 年用户年会论文
a 粉末进入熔池 b 粉末反弹
图9 粉末进入熔池和反弹
由图可知：200 个样本中有166 个进入熔池其余反弹，根据计算粉末利用率为：83%。
由ANSYS 温度场计算可知：熔池自由表面温度高达1800℃左右而约束气体温度只有
25℃，所以在熔池表面存在换热，熔池自由表面温度和表面气体温度如图10 所示：
a 熔池自由表面温度 b 熔池表面气体温度
图10 熔池表面换热
其中图10 a 为熔池自由表面温度，图10 b 为熔池表面气体温度，由图可知：由于约
束气体的冷却,熔池自由表面最高温度为1781℃（图10 a），同时熔池表面气体温度也从
25℃最高上升至159℃，但仅限于熔池表面的气体，距熔池表面稍远部分气体温度即不受
影响（图10 b）。
2006 年用户年会论文
根据激光成形进入准稳态时的粉斑分布、得粉率及换热边界条件，在ANSYS 中计算
激光成形温度场及熔覆层生长如图11 所示：
图11 激光成形温度场及熔覆层生长
其中图11 a 为进入准稳态时的温度场及熔覆层，图11 b~c 熔覆和结束时的温度场分
布及熔覆层形状，图11 d 为冷却时的温度场分布。对比图11 a 和11 b，可以看到进入准稳
态后熔池温度基本保持在1800℃左右，而熔池大小熔覆层平均高度为0.4mm 和平均宽度
3.5mm 基本保持恒定。
2.2.3 验证分析
通过近距离连续拍摄系统对316L 不锈钢粉末激光快速成形过程中激光、粉末与熔池交
互行为的实时观察，成形参数见表1，熔覆层、熔池自由表面及粉末的捕捉及反弹如图12
所示，图中橘红亮部为熔池自由表面，下部黑色部分为成形试样，左半部高出部分为熔覆
层，经测量其高度为0.43mm，宽度3.6mm，计算结果熔覆层平均高度为0.4mm 和平均宽
度3.5mm 与实验吻合较好。
a 进入准稳态 b. 熔覆中
c. 熔覆结束 d 冷却
2006 年用户年会论文
图12 近距离连续拍摄激光、粉末与熔池交互行为
图12 中熔池自由表面上部为气粉两相流场域，粉末流在约束气体和重力的作用下从右上方
喷向熔池表面，白色箭头1 所指为将被熔池捕捉粉末颗粒，箭头2 所指为反弹粉末颗粒，
而箭头3 所指为动量较大粉末颗粒被熔池捕捉后又反弹，表面还附着液态熔池金属。由于
本文作了固体颗粒撞到固体表面反弹而撞入液态熔池被吸收的假设，对于前两种情况作了
比较好的描述，如图13 所示：
图13 粉末与熔池交互模拟
图中白色箭头1 和2 分别与图12 图中白色箭头1 和2 相对应，所指为1 将被熔池捕作粉末
颗粒，2 所指为反弹粉末颗粒，对于粉末颗粒被熔池捕作后又反弹的情况，在实际激光快
1 2
2006 年用户年会论文
速成形中是应该尽力避免的，只有当使用大颗粒粉末或大流量约束气时才会着重考虑，所
以一般情况下使用固体颗粒撞入液态熔池被完全吸收的假设是正确的。
3 结论
激光、粉末与熔池交互是激光快速成形不可回避的基本问题，通过合理假设简化问题，
采用ANSYS 单元生死技术模拟激光加热及熔覆层的生长及熔池及自由界面形状，采用CFX
基于拉格朗日方法的粒子跟踪（particle tracking model）模型追踪粉末颗粒，根据界面温度
设置粒子平性（parallel）和垂直(perpendicular)动量损失模拟熔池对粒子的捕捉以及工件对
粉末的反射，计算表明对于316L 不锈钢粉末、激光与熔池的交互，计算结果与实验相吻合。
[参考文献]
[1] 陈 静 激光快速成形过程熔池行为的实时观察研究 应用激光 2005 年4 月
[2] J . Mazumder et al , Optics and L asers in Engineering 2000 (34) : 397
[3] Yunchang Fu et al, Journal of material processing technology 128(2002)106-112
[4] Jehnming Lin, Optics & Laser Technology 31 (1999) 233±238
[5] 刘振侠, 中国激光 0258-7025(2003)06-0567-04
[6] J.Powell, Ph.D.dissertation, Imperial College of Science and Technology, 1983.
[7] C.L.Atwood, M.L.Griffith, L.D.Harwell, et al, Sandia Report ,1999
[8] Ehsan Toyserkani, Optics and Lasers in Engineering 41 (2004) 849–867
[9] Ansys Inc. ANSYS User Manual, Ansys Inc., USA,1998

yhlhzg

2006 ANSYS 中国用户论文集
板坯连铸结晶器铜板－坯壳接触模型及其应用
王卫华，李永林，王凤琴，薄凤华
[首钢技术研究院，北京 100041]
[摘 要]：本文利用ANSYS有限元软件，建立了板坯结晶器铜板－坯壳接触有限元模型，并利用此模型计算分析了不同的足辊位置对板坯结晶器内坯壳受力的影响。结果表明，足辊向外侧移动将造成结晶器下沿附近节点形成明显的应力集中，对结晶器造成强烈的摩擦，将使结晶器发生局部磨损。
[ 关键词 ]：板坯连铸，结晶器铜板，有限元法，接触模型，磨损
The FEM Contact Model Between Mould Copper Plate and Slab Shell of ContinuousCasting and Its Application
Weihua Wang, Yonglin Li, Fengqin Wang, Fenghua Bo
(Shougang Research Institute of Technology, Beijing100041)
Abstract: By using ANSYS FEM software, the authors build a mathematics model and investigate the influence of foot roller position on the abrasion of mould copper plates during slab continuous casting. If the foot rollers are moved away from slab shell, the pressure in the lower part of mould copper plate will increase. This will cause high friction force between the slab shell and the copper plate. The abrasion of mould copper plate will increase in some area.
Key words: Slab Continuous Casting, Mould Copper Plate, FEM, Contact model, Abrasion
1. 前言
有限单元法是20世纪中叶在电子计算机诞生之后，在计算数学、计算力学和计算工程科学领域里诞生的最有效的工程计算方法。经过40多年的发展不仅使各种不同的有限元方法形态相当丰富，理论基础相当完善，还发展了一大批先进的有限元分析软件，其中美国ANSYS公司开发的计算机模拟工程结构有限元分析软件ANSYS就是一个目前被广泛应用的软件。
本文利用ANSYS通用有限元软件建立板坯结晶器铜板与坯壳接触模型。并利用此模型对足辊对中情况引起的钢连铸结晶器铜板与坯壳的接触应力进行有限元模拟，分析各种足辊位置对结晶器铜板磨损受力的影响，从而为减少因机械应力而造成的结晶器不均匀磨损提供理论依据。
造成结晶器失效的方式主要包括结晶器变形和结晶器磨损两种，它们与结晶器冷却系统设计、结晶器传热状况、结晶器铜板和镀层的材质及结构、结晶器支撑系统设计、足辊对中等因素密切相关。首钢板坯连铸机在使用中出现的主要问题有：铜板磨损，即过钢量达到一定值后，结晶器铜板镀层磨损达到规程上限；铜板变形，即结晶器宽面出现鼓肚；窄面和足辊都出现了不同程度的位移。
2. ANSYS接触分析概述
在有限元结构力学模拟中，非线性包括三种情况：材料非线性、状态非线性和几何非线性。在状态非线性分析中，接触分析是一种很普遍的非线性行为。接触是状态变化非线性类型中一个特殊而重要的子集，是一种高度非线性行为，需要很大的计算资源。为了进行有效的计算，理解问题的特性和建立合理的模型是很重要的。接触问题存在两个较大的难点：其一在求解问题之前接触区域表面之间的接触或分开是未知的，它随载荷材料边界条件和其他因素而定；其二，大多数的接触问题需要计算摩擦，摩擦使问题的收敛性变得困难。
2006 ANSYS 中国用户论文集
接触问题分为两种基本类型：刚一柔接触和柔一柔接触。为了给接触问题建模，首先必须认识到模型中的哪些部分可能会相互接触，以通过指定的接触单元来识别可能的接触对。
3.模型的建立
3.1模型的基本假设
根据研究重点，主要对铜板－坯壳接触模型作以下假设：
（1）结晶器铜板和足辊均为刚性，不产生变形；
（2）初始阶段结晶器与坯壳紧密接触，即不考虑气隙；
（3）只考虑钢水静压力对凝固坯壳的作用，而不考虑坯壳的热应力。坯壳温度采用结晶器内板坯凝固、流动耦合模型的计算结果。
（4）坯壳的材料本构关系为弹－塑性。
以上假设能突出坯壳与结晶器下沿的接触情况，忽略了次要的因素，使模型易于调整。
3.2模型的描述
利用ANSYS有限元软件作接触非线性分析主要分三个阶段进行：前处理，求解及后处理。前处理阶段选择单元类型，并建立窄面铜板、坯壳和足辊模型并对其进行网格划分，定义计算所需实常数、材料特性和工艺参数。前处理阶段完成后，就进入求解阶段进行分析计算，该阶段可执行如下操作：定义瞬态或稳态分析类型，施加初始条件和边界条件载荷，进入接触向导定义接触对及接触类型，定义时间步长并选择求解器，进行有限元求解。后处理过程紧接在前处理过程和求解过程之后，通过友好的用户界面，可很容易获得铜板及坯壳的应力分布，位移量及变形程度等结果。
考虑到建立模型和划分网格的工作量，以及涉及计算工作所需的大量数据。根据模型尺寸和模拟的实际自然状态，建立二维模型分析坯壳窄面的情况，考虑了铜板的锥度，将坯壳－铜板作为第一接触对，采用刚-柔接触类型；分别将三个足辊与坯壳作为第二、三、四接触对，采用柔-柔接触类型来进行受力分析。四个接触对均采用面-面接触方式。
计算采用四边形八节点的183结构单元，其中接触对采用169、172二维接触分析单元。为保证模型精度，模型共划分2940个单元。
计算中坯壳沿拉坯方向的厚度分布由流动、传热和凝固耦合模型分析给出，考虑了沿高度方向分布的结晶器瞬时热流密度，并与漏钢坯壳进行过对照。模型计算中钢水静压力对坯壳的影响，采用带液芯坯壳应力计算中常用的“剔除液芯法”[1-2]，即根据铸坯温度场计算结果剔除掉液芯单元，把沿高度分布的钢水静压力直接施加到坯壳的凝固前沿上。由于坯壳受力情况比较复杂，为突出关心的重点，只考虑凝固坯壳所受机械应力的影响。几何模型以弯月面与铜板的交点为原点，以水平方向中弯月面从外到内方向为x正方向，竖直方向为y方向，拉坯方向与y轴正向相同。考虑到结构、边界条件的对称性，只取沿窄面中心线的1/2坯壳、结晶器窄面铜板和三个足辊。为了得到结构化网格，建模时将坯壳在纵向1/2处的剖开。
3.3模型的主要参数
根据凝固模型计算，得到结晶器下沿附近凝固坯壳的温度。
板坯连铸机与模型相关的主要参数为：铸坯断面尺寸1800×220mm；结晶器高度900mm，直型；铜板厚度45mm；窄面锥度7mm；足辊宽、窄面各三排；弯月面距结晶器上口100mm。
参考相关文献中的物性参数，采用钢在高温下的物理性能，取弹性模量为15GPa，泊松比为0.37，屈服强度为20MPa[3-4]。
2006 ANSYS 中国用户论文集
3.4.模型的计算及结果分析
将模型处理为稳态进行计算，描述首钢板坯连铸机坯壳分别在以下三种情况下结晶器铜板－坯壳的接触和受力情况：
（1）足辊工作面与结晶器铜板下沿内侧平行；
（2）足辊沿x轴负向移动；
（3）足辊沿x轴正向移动。
图1左上为二维几何模型，其余三幅图为坯壳在三种足辊位置时结晶器下沿到第一足辊局部的x方向应力云图。由图1可见，在钢水静压力的作用下，如果足辊向－x方向移动，则坯壳不与第一足辊接触，在结晶器下沿处形成一个应力集中区；如果足辊不移动，则坯壳与第一足辊接触处形成应力集中区；在足辊沿＋x方向移动时，坯壳与第一足辊接触处也形成应力集中区，且x方向应力较足辊不移动时大，同时坯壳与结晶器下沿脱离接触。
图1 坯壳x方向应力云图
（左上：几何模型图；右上：足辊沿x轴负向移动；
左下：足辊不移动，右下：足辊沿x轴正向移动）
left top: geometry model; right top: foot roller moving away form slab shell
left below: foot roller has the same position as mould copper
2006 ANSYS 中国用户论文集
right below: foot roller moving close to slab shell
Figure 1, X direction stress nephogram of slab shell
为了进一步考察坯壳与结晶器的接触情况，只讨论与结晶器接触的一段坯壳。总共考虑六种足辊位置对坯壳受力的影响，即-0.3mm，-0.2mm，-0.1mm，0mm，0.1mm，0.2mm，0.3mm。首先观察各种足辊位置时坯壳的受力情况，见图2。 400500600700800900-1.0x105-8.0x104-6.0x104-4.0x104-2.0x1040.0 x方向应力，Pa到弯月面距离，mm足辊移动 0.1mm 0.2mm 0.3mm 0mm
a, 足辊向偏离坯壳方向移动 b足辊向靠近坯壳方向移动
图2 坯壳在各种足辊位置时的受力情况
a. foot roller moving away from slab shell b. foot roller moving close to slab shell
Figure 2 Stress distribution on the slab shell with different foot roller position
由图2可见，当足辊向－x方向移动时，结晶器下沿节点受力随移动量增大而迅速增大，最大应力都出现在结晶器下沿（800mm），最大应力都大于足辊不移动的情况。当足辊向＋x方向移动时，最大应力随移动量增加而逐步增大，最大应力出现的位置也发生向上移动的趋势，但最大应力都小于足辊不移动的情况。
板坯浇铸过程中，坯壳在钢水静压力作用下，易被压向铜板工作面，这种情况在结晶器下部尤其明显。当铜板和坯壳接触产生摩擦，易使得镀层脱落和铜板磨损。同时如果结晶器背板支撑系统不合理，导致浇铸过程中在钢水静压力作用下，足辊区开口度变化，使得坯壳发生扭曲，加大铜板磨损。以上情况分析表明，即使足辊向外侧（－x方向）移动一个很小的距离，都将引起坯壳受力的急速增大。由于坯壳受力的同时也将力反作用于铜板，对铜板法向的压力在拉坯作用下将产生巨大的摩擦力，造成铜板的局部磨损。
5.结论
1）利用ANSYS有限元软件建立结晶器铜板－坯壳接触二维模型，并只考虑钢水静压力对凝固坯壳的作用，而不考虑坯壳的热应力，计算分析了不同的足辊位置对板坯结晶器内坯壳受力的影响，得出结晶器局部磨损为结晶器失效的主要方式。
2）计算分析结果表明，足辊的外移会造成结晶器下沿附近节点形成明显的应力集中，坯壳的反作用力对结晶器造成强烈的摩擦，加剧结晶器的局部磨损。
参考文献：
[1] Kelly J E, Michalek K P, O’Connor T G, et al. Initial development of thermal and stress fields in continuously cast steel billets. Metallurgical transaction A, 1988, 9A(10): 2589~2602
[2] A.Grill,J.K. Brimacombe and F.Weinberg, Mathematical analysis of stresses in continuous casting of steel, Ironmaking and Steelmaking,1976,Vol.3,No.1,38~47
2006 ANSYS 中国用户论文集
[3] Mizukami H, Murakami K, Miyashita Y. Mechanical properties of continuously cast steel at high temperatures, Tetsu-to-hagane (Iron and Steel, in Japanese), 1977, 63: 146, s652
[4] Uehara M, Samarasekera I V, Brimacombe J K. Mathematical modeling of unbending of continuously cast steel slabs. Ironmaking and Steelmaking, 1986, 13(3): 138~153

yhlhzg

2006 年用户年会论文
1
箱-椭孔连轧过程三维塑性变形有限元模拟
管志杰 洪慧平 康永林 朱国明
（北京科技大学材料科学与工程学院，北京100083）
摘 要 借助有限元分析软件ANSYS/LS-DYNA，对方坯在箱-椭孔连轧生产过程进行了模拟仿真。分
析了轧制区等效应力场和等效应变场，并分析了在方坯横截面上等效应变和等效应力的分布以
及横截面上两处的等效应变随时间的变化曲线。改变孔型参数后重新模拟了轧制过程，把得到
的模拟结果与原始孔型下的结果进行比较，得出了增大箱形孔型内圆角半径和减小椭圆孔辊缝
有利于热轧过程中裂纹压合的观点。
关键词 箱形孔型 椭圆孔型 连轧 有限元法 ANSYS/LS-DYNA
Simulation of Continuous Rolling with Box-Oval
passes by Finite Element Method
Guan Zhijie, Hong Huiping, Kang Yonglin, Zhu Guoming
(School of Materials Science and Engineering, University of Science and Technology,
Beijing 100083)
Abstract The numerical simulation of square billet during continuous dynamic hot rolling with box-oval passes
has been conducted by finite element method with software ANSYS/LS-DYNA, and the equivalent
stress field and the equivalent strain field of rolled area are analyzed in this paper. Besides,
distribution of the equivalent stress field and the equivalent strain field in the cross-section as well as
the curves of the equivalent strain correlated to time of the two points in the cross-section have been
calculated. Another similar simulation of rolling process has been achieved after some parameters of
passes are changed. The results show that it is beneficial to compact the inner cracks by enlarging
fillet radius of box passes and reducing gap of oval passes.
Keyword Box pass Oval pass Continuous Rolling FEM ANSYS/LS-DYNA
1 前言
型钢轧制过程是一个复杂的塑性变形过程，孔型前后金属断面差异大、形状复杂，且
在轧制过程中金属各部分受力不均、变形不均，轧件尺寸难以精确计算。轧件在孔型中发
生复杂的三维变形，对该过程进行理论解析非常困难。随着计算机三维模拟仿真技术的飞
2006 年用户年会论文
2
速发展，特别是大型非线性CAE 模拟仿真软件和硬件平台具有的超高性能指标，都为准确
地模拟材料加工过程金属三维大变形从而全面提高产品质量和开发新品种提供了可能[1]。
近年来，国内外一些学者利用大型有限元模拟软件模拟了型钢的轧制过程，取得了一
定的成果。R.Iankov[2]利用有限元模拟仿真软件MSC.MARC 模拟了圆坯料在扁平孔型中轧
制生产线材的三维过程，得出了最终产品的等效应变、应变率、等效应力和残余应力的分
布；并且利用MSC.MARC 中的等效平面应变法模拟出横截面上等效塑性应变、全应变和
应力的分布。C.S.Li 等[3]利用有限元法模拟了50CrV4 汽车齿轮圆钢在连轧过程中的表面、
心部、1/2 半径和平均温度的分布情况。洪慧平等[1][4] 采用有限元模拟仿真软件
MSC.MARC/Auto- forge 对大圆钢两道次连轧过程进行了三维热力耦合模拟，得到了大规格
合金芯棒圆钢热连轧时各道次轧辊孔型中金属的三维变形，并准确地计算了轧件的应力场、
应变场、温度场和相关力能参数的动态分布值(如轧制力、轧制力矩等)。李建超等[5][6]借助
显式动力学分析有限元软件ANSYS/LS-DYNA 对方坯分别在箱型孔和椭圆孔中轧制情况
进行了模拟，分析了轧制过程中轧件的变形情况、轧制区应力场和力学参数，并用实验数
据和模拟数据作比较，从而证明了数值模拟的可行性。
本文是基于ANSYS/LS-DYNA 对方坯在箱-椭圆孔型中连轧过程进行有限元模拟，同
时在改变两个轧辊的孔型参数的基础上重新模拟。通过对比孔型改变前后的模拟结果，找
出孔型改变对钢坯内裂纹愈合的影响。
2 有限元模型建立
2.1 原始孔型有限元模型建立
本次模拟采用的钢种为45 钢，成品尺寸为Φ48mm，坯料尺寸为75mm×75mm×1000mm，
设定坯料有内圆角，半径为10mm。原始孔型粗轧前两道次的孔型结构如图1 所示，这两
个孔型只是连轧孔型系统的一部分。相关轧制和设备工艺参数见表1。
(a) 箱形孔型 (b) 椭圆孔型
图1 方坯连轧前两道次孔型
Fig .1.1st pass and 2nd pass of continuous rolling square billet
2006 年用户年会论文
3
表1 两道次连轧初始工艺设备参数
Tab.1. The initial rolling parameters of two-pass continuous rolling
孔型 轧制速度（m/s） 轧辊转速（r/min） 最大辊径（mm） 辊缝（mm） 延伸系数
箱形 0.244 7.882 685 18 1.2427
椭圆 0.325 9.86 685 20 1.3319
首先在ANSYS 的前处理器中定义单元和材料特性，坯料和轧辊均采用SOLID164 单
元，假设轧制过程温度恒为一定值，然后确定材料特性，其中包括密度、杨氏模量、屈服
强度、泊松比和切线模数等。
由于坯料具有对称性，为了减少计算时间和节省存储空间，取1/4 轧件建立有限元模
型，考虑到坯料在孔型中发生连轧，给定坯料长度为1000mm，同时轧辊间距取500mm。
轧辊做为刚性接触体定义，轧辊与坯料之间的摩擦采用库仑摩擦模型，坯料与箱形孔和椭
圆孔的动态摩擦系数均取0.3，坯料与两个轧辊的静态摩擦系数都设定为0.35。给定坯料一
定的初始速度，轧辊与坯料接触后，摩擦力使坯料完成整个轧制过程。
在ANSYS 前处理器中建立几何模型并划分网格。坯料在长度方向取100 等份，在横
截面上划分64 个单元，这样轧件的有限元模型共采用了6400 个单元和8181 个节点。箱形
孔在圆周上划分160 份，椭圆孔在圆周上划分120 份，两个轧辊在厚度上都划分为一份，
在孔型线上都划分为15 份。建立的有限元模型如图2 所示。
2.2 改变两个轧辊孔型参数后建立有限元模型
将箱形孔孔型内圆角半径由原来的30mm 改为40mm，同时改变第二个孔型辊缝值，
把辊缝值由原来的20mm 减小到18mm，此时孔型图如图3 所示。在不改变连轧常数的情
况下，坯料过箱形孔时轧制速度增加为0.246m/s，同时轧辊转速增加到0.8274rad/s；相应
的，坯料过椭圆孔时轧制速度增加为0.335m/s，轧辊转速增到1.061rad/s。其他参数保持不
变，建立有限元模型。
图2 两道次连轧方坯有限元模型
Fig.2. Finite element model of continuous rolling square billet by two passes
2006 年用户年会论文
4
3 模拟结果的分析
3.1 原始孔型下的结果
图4 为轧件在原始阶段、第一机架咬入阶段、第二机架咬入阶段和稳定轧制阶段时网
格变形情况。轧件刚被箱形孔型咬入时，轧件的角部首先和轧辊进行接触而发生变形，轧
件的中心最后和轧辊接触。坯料在第二机架咬入时与第一机架咬入时相似，只是角部变形
更剧烈，部分网格变形加剧。在稳定轧制阶段，轧件沿轧制方向方向网格变形比较均匀。
在后三个阶段中，变形后的网格与变形前相比在高度方向被压缩，在轧制方向被拉伸，从
而使轧件变长。而且在两个机架咬入时，轧件角部表面金属在轧制方向上流动速度大于心
部的节点流动速度，则轧件角部在轧制方向延伸较大，轧件中心在轧制方向延伸较小，从
而使轧件头部中心向内凹。模拟结果中的这种轧件在孔型中的流动特征与实际轧制情况下
轧件在孔型中的流动特性相吻合，说明利用ANSYS/LS-DYNA 对金属连轧过程进行模拟仿
真分析是可行的。
45 钢在给定温度下的屈服强度约为4×107Pa。图5 为在稳定轧制阶段不同视角的等效
应力和等效应变的分布云图。由图5（a）和图5（b）可见，红色区表示的等效应力值与钢
(a) 箱形孔型 (b) 椭圆孔型
图3 两个孔型参数均改变时的孔型图
Fig.3. Two passes of changed Parameters
(a) (b) (c) (d)
a)为原始阶段；b)为第一机架咬入阶段；c)为第二机架咬入阶段；d)为稳定轧制阶段
图4 两机架连轧不同阶段的变形情况
Fig.4. Deformation of square billet during continuous rolling in two passes
2006 年用户年会论文
5
的屈服强度相当，说明钢在轧制过程中处于塑性变形状态。在咬入过程中，等效应力逐渐
变大，一直达到塑性状态时的屈服应力值，说明轧件的咬入过程就是塑性变形发生的过程。
由图5（c）和图5（d）可见，与轧辊首先接触的坯料角部在两辊轧制后积累的等效应变值
较坯料心部大，坯料经过第二道次后的等效应变比经过第一道次后的等效应变明显增大，
与孔型内圆角接触的钢的等效应变比中间部分的钢的等效应变值大。
为了更好地观察钢在轧制过程中的等效应力和等效应变分布，取变形过程中的横截面
做分析，在这里取两个轧辊正下方的面即第一道次下坯料横截面和第二道次下坯料横截面。
横截面上等效应力和等效应变的分布分别如图6 和图7 所示。
由图6（a）可知，在第一道次下坯料横截面上左边部的金属的等效应变值明显比其它
位置的值大；图6（b）表明处于横截面右上方区域的金属还处于塑性状态。图7 为第二道
次下坯料横截面上等效应变和等效应力的分布，由图7（a）可见，横截面上的等效应变值
由左边部向右依次减小。图7（b）表明，横截面上左下角的等效应力值最小，而且从左到
右等效应力值依次增大，横截面右半部处于塑性状态。
(a) (b)
(c) (d)
a)为从模型的左向观察到的等效应力图；b) 从模型的右向观察到的等效应力图
c)为与a)对应的等效应变图；d)为与b)对应的等效应变图
图5 轧件在轧制稳定阶段下增量子步为201 时的等效应力与等效应变
Fig .5. Distribution of equivalent stress and equivalent strain in increment 201 during steadily rolling
2006 年用户年会论文
6
(a) 等效应变图 (b) 等效应力图
图6 第一道次下坯料横截面上等效应力和等效应变
Fig.6. Equivalent stress field and equivalent strain field in billet’s cross-section of the 1st pass
(a) 等效应变图 (b) 等效应力图
图7 第二道次下坯料横截面上等效应力和等效应变
Fig .7. Equivalent stress field and equivalent strain field in billet’s cross-section of the 2nd pass
图8 时间后处理采用的节点及对应单元在模型横截面上的位置
Fig.8. Location of nodes and correlative element on the cross-section of model
2006 年用户年会论文
7
为了了解在整个轧制过程某一点的等效应变随时间的变化，利用ANSYS 的时间后处
理器得出了有限元模型中单元和节点的等效应变随时间的变化曲线。所取的节点及其对应
的单元所在的位置见图8。后处理采用的节点的标号为6899 和6305，节点6899 对应的单
元是1406 单元，节点6305 对应的单元是2106 单元。为了使所取的单元和节点更具有普遍
性，所取单元和节点所在的模型横截面距离坯料模型端面为50mm。单元及其节点确定以
后，就可利用时间后处理器得出等效应变随时间的变化曲线，如图9 所示。
选用这两个单元及节点是因为成品中裂纹的分布位置与这两个位置相对应。观察这两
个位置的等效应变随时间的变化曲线，有利于分析变形过程中坯料在这两个位置的等效应
变的变化，从而找出压合裂纹的措施。由图9（a）和图9（b）两图可知，这两个位置处的
等效应变的变化趋势是一致的；所不同的是，1406 单元上6899 节点的等效应变值较2106
单元上6305 节点大，这说明1406 单元上6899 节点处的钢变形程度较大。
3.2 二个轧辊孔型参数同时改变后的结果与原始孔型下的结果对比
增大箱形孔型内圆角半径，同时减小椭圆孔型辊缝值后重新模拟，并得出了模拟结果。
图10 为孔型参数改变前后在第一道次下坯料横截面上等效应变和等效应力的对比图。由图
可见，截面上的等效应变和等效应力的分布在孔型参数改变前后的情况发生变化。
(a) (b)
a)为单元1406 上6899 节点的等效应力-时间曲线；b)为单元2106 上6305 节点的等效应变-时间曲线
图9 等效应变随时间的变化曲线
Fig.9. Curves of equivalent strain in response to time
(a) 改变孔型参数后的等效应变 (b) 原始孔型下的等效应变
2006 年用户年会论文
8
在金属发生热塑性变形的前提下，大的应变量有利于减少金属内部的缺陷，有利于金
属内的裂纹或孔洞的压合。由图10（a）和图10（b）可见，改变箱形孔型内圆角后，坯料
左边部的等效应变值变大，而且范围扩大。这说明改变内圆角后，坯料左边部应变进一步
深入，有利于左边部缺陷（如孔洞或裂纹）的压合。改变箱形孔型内圆角后，在轧辊的正
下方的坯料的等效应力的分布如图10（c），从图中可见改变孔型参数后，坯料处于塑性
变形状态的区域增大。
图11 为孔型参数改变前后在第二道次下坯料横截面上等效应力和等效应变的对比图。
由图11（a）图11（b）可知，改变孔型参数后模拟的横截面右上边部的等效应变变大，说
明此处变形进一步深入，这对压合于右上边部缺陷（如孔洞或裂纹）是有利的。
因此，增大箱形孔内圆角半径和减小椭圆孔辊缝值有利于压合边部的缺陷（如孔洞或
裂纹等）。
(c) 改变孔型参数后的等效应力 (d) 原始孔型下的等效应力
图10 孔型参数改变前后第一道次下坯料横截面上等效应力和等效应变对比
Fig .10. Comparison of equivalent stress field and equivalent strain field in billet’s cross-section
of the 1st pass after changing parameters of passes
(a) 改变孔型参数后的等效应变 (b) 原始孔型下的等效应变
2006 年用户年会论文
9
图12 为单元1406 上6899 节点处的等效应变随时间的变化曲线在孔型改变前后的对比
图。由图可知，两个轧辊孔型参数同时改变后，在第一道次轧制时单元1406 上6899 节点
处的等效应变值比原始孔型下增大0.06 左右，而且在轧制过程中该处发生稳定的塑性变形，
由于在稳定塑性状态下，较大的应变是有利于裂纹和孔洞的压合，所以增大箱形孔型内圆
角半径和减小椭圆孔型辊缝值变孔型参数有利于此处孔洞或裂纹的压合。
图13 为单元2106 上6305 节点处的等效应变随时间的变化曲线在孔型改变前后的对
比。由图可见，在改变参数后, 单元2106 上6305 节点处在经过第一道次轧制时的等效应
变值比在原始孔型下的等效应变大0.03 左右，这说明增大箱形孔型内圆角半径和减小椭圆
孔型辊缝值对抑制或压合坯料中此处的裂纹或孔洞等缺陷有积极的作用。
(c) 改变孔型参数后的等效应力 (d) 原始孔型下的等效应力
图11 孔型参数改变前后第二道次下坯料横截面上等效应力和等效应变对比
Fig .11. Comparison of equivalent stress field and equivalent strain field in billet’s cross-section
of the 2nd pass after changing parameters of passes
(a) 改变孔型参数后等效应变-时间曲线 (b) 原始孔型下等效应变-时间曲线
图12 等效应变随时间的变化曲线在孔型改变前后的对比（1406 单元上6899 节点）
Fig.12. Comparison of curves of equivalent strain in response to time
by changing parameters of passes(1406 element and 6899 node
2006 年用户年会论文
10
这说明：与原始孔型相比，改变孔型参数后，实际坯料中与后处理采用的单元1406 上
6899 节点和单元2106 上6305 节点对应的裂纹或孔洞将易于被压合。
4 结论
1) 用ANSYSY/LS-DYNA 模拟出了方坯在箱-椭圆孔型中的连轧过程，模拟结果中的
轧件在孔型中的流动特征与实际轧制情况下轧件在孔型中的流动特性相吻合，说明利用
ANSYS/LS-DYNA 对金属连轧过程进行模拟仿真分析是可行的。
2) 增大箱形孔内圆角半径和减小椭圆孔辊缝值后，第一道次下坯料横截面左边部应变
进一步深入，第二道次下坯料横截面上右上边部变形也进一步深入，这对边部缺陷（如孔
洞或裂纹等）的压合是有利的。
3) 增大箱形孔内圆角半径和减小椭圆孔辊缝值后，与后处理采用的单元和节点处对应
的裂纹或孔洞易于被压合。
[参考文献]
[1] 洪慧平,康永林,冯长桃等.连轧大规格合金芯棒钢三维热力耦合模拟仿真.钢铁,2002,37(10):23～42
[2] R.Iankov.Finite element simulation of profile rolling of wire.Journal of Materials Processing Technology,2003,
142:355～361
[3] C.S.Li,X.H.Liu,G.D.Wang.Simulation on temperature field of 50CrV4 automobile gear bar steel in continuous
rolling by FEM.Journal of Material Processing Technology,2002,120:26～29
[4] 洪慧平,康永林,冯长桃等.椭-圆孔型连轧大圆钢三维热力耦合弹塑性有限元分析.特殊钢,2002,23(5):5～8
[5] 李建超,吴迪,赵宪明.箱型孔轧制方坯的显式动力学有限元模拟.特殊钢,2002,23(2):23～25
(a) 改变孔型参数后等效应变-时间曲线 (b) 原始孔型下等效应变-时间曲线
图13 等效应变随时间的变化曲线在孔型改变前后的对比（2106 单元上6305 节点）
Fig.13.Comparison of curves of equivalent strain in response to time by changing
parameters of passes(2106 element and 6305 node)
2006 年用户年会论文
11
[6] 李建超,崔建忠,吴迪等.椭圆孔轧制方坯的显式动力学模拟.塑性工程学报,2003,10(1):57～59

yhlhzg

2006 年用户年会论文
1
用有限元法对铸铁件铸造缺陷的预测
陈玲1，郑清春1，安莉珠2，殷飞1
1 天津理工大学，天津 300191 2 天津市第二维修齿轮厂，天津 300113
[ 摘 要 ] 以有限元法为理论基础，ANSYS 软件为平台，建立铸件凝固过程的数学模型，对铸造凝固过程
进行了温度场、温度梯度场和应力应变场的数值模拟。模拟时考虑相变中结晶潜热、材料随温
度变化的非线性因素。得到了铸件在凝固过程的各种瞬态场的分布和变化规律。预测出了铸件
可能出现缩孔、缩松等铸造缺陷的位置；同时探讨了判断热裂纹的方法，结果表明模拟预测与
生产验证结果基本一致。
[ 关键词 ] 铸造凝固过程, 缩孔缩松，热裂，应力应变场
Forecast of Casting Drawback with
the Finite Element Analysis
CHEN Ling1, ZHENG Qing-chun1,AN Li-Zhu2 YIN Fei1
1 Tianjin University of Technology, Tianjin 300191
2 The Second Gear Maintenance Factory of Tianjin Tianjin 300113
[ Abstract ] The mathematic model in the solidification process of the cast is built up with finite element
method and ANSYS software. Considering the material nonlinear caused by the
temperature and the crystal latent heat in the process of the phase transformation, the
process of the cast solidification is simulated, then get the transient temperature field, the
gradient field, the stress and strain field, the distribution and the change regulation of in
solidification process. The positions where the shrinkage hole and shrinkage porosity may
be occurred will be forecasted by analysis. Meanwhile, the method with which the heat crack
is judged is discussed. The result shows that the simulated result is correspond to the
experimental result.
[ Keyword ] Cast solidification process, Shrinkage hole and shrinkage, Heat crack ,Stress and strain
field
作者简介：陈玲 （1964-），女，教授
天津市高等学校科技发展基金项目（20051010）
2006 年用户年会论文
2
1 前言
铸造行业是一个古老而又非常重要的传统行业，无论是在发达国家还是在发展中国家，
铸造行业都有着不可替代的地位。铸造过程的计算机模拟仿真技术已成为当今国际公认的
制造科学与材料科学的主要前沿领域，也是改造传统铸造产业的必由之路。由于铸件在形
成过程中会产生缩孔缩松，热裂纹等等，这些缺陷对产品的质量和使用性能具有潜在的危
害性，但长期以来铸造的质量很难得到保证。在实际生产之前对铸件在浇注、凝固过程中
可能产生的缺陷采用计算机数值模拟的方法进行有效的预测，从而在浇注前就采取有效对
策，以确保铸件的质量，是可行的途径之一。本文对铸造凝固过程的瞬态温度场和应力场
进行了数值模拟，预测可能产生缩孔缩松以及热裂纹的位置，并且与实际铸件进行了验证
对比；同时探讨了铸件热裂纹的判断方法，为降低铸件的废品率提供了参考。
2 热传导微分方程及边界条件
在一般三维问题里，瞬态温度场的场变量T（x,y,z,t）在直角坐标中应满足的微分方
程是
2 2 2
2 2 2
c T k T T T g
t x y z
ρ
∂ ⎛∂ ∂ ∂ ⎞
∂ − ⎜⎝∂ +∂ +∂ ⎟⎠=
（在Ω内） (1)
边界条件是：第一类边界条件：T f (x, y, z, t )
1
= Γ
（2）
第二类边界条件： q
n
k T =
∂
∂
Γ 2
（3）
第三类边界条件： ( ) Γ
Γ
= −
∂
∂ h T T
n
k T
3
（4）
式中： Ｔ—温度 ℃，ｔ—时间ｓ，ｃ—比热Ｊ/(ｋｇ·℃)，ρ—密度ｋｇ/ｍ3，k—热
传导系数 Ｗ/(ｍ·℃)， g—物体内部热源强度 Ｊ/(ｓ·ｍ3)
式（1）为热量平衡方程，表明微体升温所需的热量应与传入微体的热量以及微体内热
源产生的热量相平衡。求解瞬态温度场问题是求解在初始条件下，即在T=T0 (当t=0 )条
件下满足瞬态热传导方程及边界条件的场函数T ，T 应是坐标和时间的函数。
3 温度场模拟分析
铸造系统温度场和应力应变场的数值模拟过程中，各冷却阶段由于材料的性能取决于
温度，因此是典型的非线性瞬态热分析，模拟中对结晶潜热和边界条件的处理尤为重要。
铸件在凝固过程中存在金属由液相转变为固相的相变过程，因此需要考虑潜热的释放,对潜
热的处理有各种方法，本文通过定义材料随温度变化的热焓方法来考虑潜热[1]，由于 ANSYS
软件中有热焓的功能，单位是J/m3，故在模拟中直接输入用公式（5）计算的数据即可。
2006 年用户年会论文
3
H = ∫ρc(T )dT （5）
在铸型与空气之间存在着热辐射和热对流，考虑到铸型外表面在浇筑后温度不是很高，
将热辐射的影响换算成砂型的对流系数。根据相关资料，在冷却过程中，铸型边界的温度
变化不是很大，可取为常数，在分析时砂型顶部、底部和侧面采用了不同的数值。
3．1 建立铸件有限元模型
铸件材料为灰铸铁，图1 示出铸件的几何模型，铸件和铸型材料特性均随温度变化。
本文在ANSYS 软件环境下，采用热分析单元solid70 建立有限元模型。铸造系统的初始条
件，即ｔ=0 时系统各部分所具有的温度：砂型为30℃，铸件为1500℃(浇注温度)；铸型与
空气之间采用对流边界条件，砂型顶部对流系数可选为7，底部为6，四个侧面为9。
首先对支架铸件进行温度场分析，得到铸件在凝固过程中温度场以及温度梯度分布，
从而预测缩孔缩松等缺陷位置；同时为求解铸件应力应变场打基础。
图1 铸件的几何模型 图2 3500 秒时温度分布图
Fig.1 Model of Cast Fig. 2 3500s Temperature Distribution
3．2 缩孔缩松的预测结果
缩孔、缩松产生的主要原因是体积收缩引起的。预测铸件中缩孔及缩松的发生部位及
其程度是实现铸件质量预测和控制的一个重要内容。文中采用新山判据G R [2]，即温度梯
度和凝固速度平方根的比值判断缩孔的产生，判据的临界值的取值范围为0.1～1.2，数值
越大，缩松倾向越大，一般取0.7 作为临界值时，通常与实际情况符合得较好。即：
( )( ) 2 1 2 1 T T t t
G
R
G
− −
= （6）
其中：G—某一节点在某时间的温度梯度值 ℃/m，R—温度变化速度 ℃/s，（T2-T1）—某
节点在某一段时间的温度差℃，t2-t1—该节点的时间差s。当G R <0.7 时，存在缺陷；当
G R ≥0.7 时就没有缺陷。
图2 和图3 采集了时间为3500 秒和44 秒时支架温度和温度梯度分布云图，表1 列举
了计算结果。其中1682，1683，1684，2311 四个节点的G R < 0.7，是可能产生缩孔缩松
的危险的位置，具体的位置请见图4。经工程实际验证，该位置恰是产生缩孔缩松的地方。
2006 年用户年会论文
4
图3 44 秒梯度分布云图 图4 铸件缺陷位置图
Fig. 3 44s Gradient Distribution Fig.4 Drawback Position of Cast
表1 计算结果
Table1 Results of Calculation
NODE TEMP1（℃） GRA（℃/m） TEMP2（℃） ΔT（℃） Δt（s） R1/2（℃/s） G/ R1/2
1682 152.85 6.91 144.71 8.14 500 0.13 0.53
1683 198.24 8.00 184.10 14.14 500 0.17 0.47
1684 235.71 9.61 214.89 20.82 500 0.20 0.48
2311 287.39 19.33 265.32 22.07 250 0.30 0.64
486 196.83 12.17 195.13 1.7 250 0.08 1.52
1663 328.42 98.77 320.60 7.82 250 0.18 5.49
145 189.23 72.76 186.80 2.43 250 0.10 7.28
552 236.39 525.68 213.18 23.21 134.86 0.41 12.82
4 热应力计算
当物体各部分温度发生变化时，物体由于热变形而产生线应变。如果物体各部分的热
变形不受任何约束，则物体上有变形而不引起应力。但是物体由于约束或各部分温度变化
不均匀，热变形不能自由进行时，则在物体产生应力。
热裂纹是铸件在凝固过程收缩受到阻碍产生于晶界的一种缺陷，要实现对该缺陷的预
测，就需分析铸造凝固过程应力场。由于产生铸造裂纹的机理非常复杂，矶部、刘驰、P N
Hansen、张家泉等人提出若干种判断方法[3]，但是，至今也未找出一种通用的热裂纹判据。
其中张家泉提出以第四强度理论中的等效应力作为判据的可操作性较强，他认为热裂属于
表观脆性实为塑性的断裂，具体判据如下：
b
W K
σ
σ
= （7）
2006 年用户年会论文
5
式中：σ—等效应力MPa，σb—材料凝固过程中抗拉强度极限值Mpa。
为不失一般性，首先分析与本文研究铸件相同材料的标准应力框的应力场，找出温度
与抗拉强度极限σb 的关系，再以此为基础，分析支架铸件，预测缺陷发生的位置。图5 和
图6 是在1400 秒时捕捉的标准应力框温度场和应力场的云图；表2 列出了600 秒和1400
秒时部分节点的应力值。在考察了600s、1400s、2200s、3000s、3800s 等时间段的等效应
力值σ 后，将温度与凝固过程中的抗拉强度极限的关系拟合，发现呈线性关系：请见公式
（8）及图7 所示。其中X 轴是节点的温度，Y 轴是强度极限，依据图7 可知相同铸件材料
在任一温度下的抗拉强度极限。
Y (1500 X) 16
= 1 &#8722; （8）
]
图5 1400s 温度分布 图6 1400s 应力分布 图7 抗拉强度极限随温度变化
Fig. 5 1400s Temp. Distribution Fig.6 1400s Stress Distribution Fig.7 Resistance Intension with Temp.
表2 600 秒和1400 秒应力值
Table 2 Stress at 600s and 1400s
600s 1400s
节点 主应力和等效应力 （MPa） 主应力和等效应力 （MPa）
σ1 σ2 σ3 σ σ1 σ2 σ3 σ
8 -3.1 -3.9 -11 7.56 -7.1 -9 -26 18
1263 -0.99 -4.5 -7.8 6 -2.3 -10 -18 13.6
972 -2.4 -2.9 -9.5 6.9 -5.7 -6.8 -22 15.8
862 -1.4 -5.3 -5.7 4.1 -3.4 -12 -13 9
为分析出铸件的热裂缺陷，本文在500～1400℃区间对该铸件进行了模拟计算，得到
等效应力值σ。再用凝固过程中的抗拉强度极限σb 与节点温度之间的函数关系式
Y (1500 X) 16
= 1 &#8722; 求出抗拉强度极限，找出铸件哪些节点的应力值偏大，再运用判据（7）计
算Kw，判断该值是否大于1，从而预测发生热裂的部位，具体数据请见表3。
2006 年用户年会论文
6
表3 铸件应力与抗拉强度极限比较表
Table3 Comparison of Stress and Intension for cast
时间 （s） 温度(℃) 应力值σ （MPa） 抗拉强度极限σb （MPa） Kw=σ/σb
92 1106.9 0.117 24.6 0.005
117 1075.7 0.211 25 0.008
242 920 0.284 34 0.0084
442 720 0.548 46 0.012
542 635 0.610 51 0.012
642 580 0.726 57.5 0.013
742 538 0.900 60.6 0.015
通过上表的计算可以看出铸件的应力值比凝固过程中的抗拉极限强度值要小，所以铸
件不会产生热裂缺陷；还可看出无论是抗拉强度极限值还是模拟的应力值都是随着铸件的
节点温度的降低而增加，这与符合实际情况相符。
5 结论
本文以真实铸件作为分析对象，采用数值模拟的方法，得出铸造凝固过程温度场分布
规律，结合新山判据，得到了可能产生缩孔缩松的位置，该位置与生产实际吻合，说明用
该法判据缩孔缩松缺陷的准确度较高；在应力应变场分析中，本文提出了铸件在铸造凝固
过程中的抗拉强度极限随温度变化的关系式，并且用此关系式分析了铸件的应力场分布情
况，结合第四强度理论得到了判断铸件发生热裂纹的方法。得出的结论，旨在为工程设计
人员改进铸造工艺提供依据。
[参考文献]
[1] 陈玲，王鹏林，张敬宇等. 金属型铸造温度场的有限元数值模拟及确定热物理参数的“逆方法”. 机
械设计[J]，2000，20（12）：29-31
[2] [日]大中逸雄.计算机传热凝固解析入门[Ｍ].许云祥译，北京：机械工业出版社，1985
[3] 王旭明，程军，党惊知，铸件热裂纹预测的判据[J]，华北工学院学报，1996，17（1）：37-41
[4] 程军、柳百成，有限元在铸件凝固过程中应用[J]，大型铸锻件，1989，第二期，14-22
2006 年用户年会论文
7
[5] 陈玲，王鹏林，杨秀萍编写.计算计辅助工程分析——ANSYS 导学[M].天津：天津理工学院机械教研
室与振动研究室，1999
[6] 陈玲，钱灵，王鹏林.基于ANSYS 的三维铸铁件缩孔缩松形成的数值模拟[J].机械设计，2004，21
（7）：12-14
[7]赵宇宏，程军，徐宏等.铸钢件凝固过程热应力场数值模拟研究[J].华北工学院学报，1996，19（4）：
306-308
[8] 王勖成，邵敏.有限单元法基本原理和数值方法[M].北京：清华大学出版社，1997 年3 月，第二版。
[9] Hou-Cheng Huang, Asif S.Usmani.Finite Element Analysis for Heat Transfer[M],Spriger-verlag,
London,1994

sonwy

大哥，怎么下阿？