基于生死单元法的双辊铸轧过程热-力耦合数值模拟

黄华贵　刘文文　王　巍　杜凤山

1.国家冷轧板带装备及工艺工程技术研究中心，秦皇岛,0660042.燕山大学，秦皇岛,066004

基于生死单元法的双辊铸轧过程热-力耦合数值模拟

黄华贵1,2刘文文1,2王巍2杜凤山1,2

1.国家冷轧板带装备及工艺工程技术研究中心，秦皇岛,0660042.燕山大学，秦皇岛,066004

以二辊φ160 mm×150 mm立式铝带实验铸轧机为对象，基于MSC.Marc商用有限元软件及其二次开发接口，引入纯铝固-液两相材料本构模型和界面压力热阻数学模型，建立了双辊铸轧过程热-力耦合非线性有限元模型。采用生死单元法模拟铝液的连续浇入，解决了辊套和铸轧区铝带材间的连续耦合传热问题。通过数值模拟，给出了铸轧速度、浇铸温度、熔池高度等因素对KISS点位置和辊面温度时间历程的影响规律，并对典型工况温度模拟结果进行了实验验证。

双辊铸轧；温度场；有限元法；生死单元法；数值模拟

0　引言

双辊薄带铸轧工艺是一种集金属快速凝固和热轧于一体的复合工艺，具有短流程和低能耗等优点，属于冶金及材料研究领域的前沿技术[1]。近年来，国内外学者围绕铸轧区金属熔体及带坯与铸轧辊间的换热行为、铸轧带材组织性能控制以及铸轧工艺，开展了大量的理论与实验研究工作[2-7]。经过近十年的快速发展和应用实践，该技术已逐渐被应用于铝合金、低碳钢以及不锈钢和硅钢等难变形金属带材的工业生产[8]。作为双辊铸轧工艺的重要基础，铸轧区内固-液相组成及其对轧制力和耦合传热行为的影响一直是国内外学者研究的热点，也是铸轧辊冷却结构的设计和铸轧工艺优化的重要依据。

数值模拟技术作为铸轧过程仿真的重要手段[9]，可实现铸轧过程热-力耦合建模和辊体热平衡的动态仿真分析，其核心技术包括：金属熔体固-液混合相(或半固态)变形抗力模型的建立、铸轧区带坯与铸轧辊表面换热系数的计算、熔体连续浇入模型的实现等。在实际铸轧生产过程中，铸轧区内金属形态主要由液相、固相和固-液两相混合体组成，且需要轧制多圈后辊体温度场才会达到平衡状态。一般认为，金属熔体具有牛顿流体的流变性能，固相金属的塑性流变性能目前已有较为成熟的理论模型，而固-液两相混合体铸轧变形过程的流变应力曲线则需通过半固态金属的热压缩实验测试获得[10-11]。湛利华等[5,12]在Gleeble上开展了铝合金连续铸轧过程流变行为物理模拟研究，并将其成功应用于铸轧工艺热-力耦合分析。铸轧区金属熔体、固相带坯与辊面之间换热行为方面，国内外学者主要采用接触热阻测试实验，得到不同温度、接触压力和表面粗糙度等条件下的接触换热系数，并应用于铸轧辊温度场模拟[13-15]。为获得铸轧辊辊体的平衡态温度场，对铸轧辊与带坯进行耦合传热有限元建模时，持续多圈的铸轧过程模拟需要大量的网格作为支持，这给模型求解带来较大困难。因此，解决金属熔体连续浇入建模问题，对实现铸轧过程热-力耦合和辊体热平衡的动态仿真分析具有重要意义。

本文以燕山大学二辊φ160 mm×150 mm铸轧机为对象，以铸轧区流变本构模型和接触换热数学模型为基础，通过MSC.Marc二次开发子程序接口，采用生死单元法建立了纯铝双辊铸轧过程的热-力耦合有限元模拟模型，解决了金属熔体的连续浇铸建模问题。目前，生死单元法在带钢卷取[16]、焊接[17]和工程施工[18]等工艺研究方面已有较多应用。

1　双辊铸轧热-力耦合有限元建模

1.1模型简化

二辊φ160 mm×150 mm立式实验铸轧机如图1所示，设计最大轧制力50 kN(5 t)，最大铸轧力矩1200 N·m，主电机功率为3 kW，具有振动铸轧功能[19]，以改善铸轧带材的内部质量。为便于有限元建模，进行如下假设：①铸轧辊为刚性辊，铝液入口速度和温度恒定，带坯与辊面不发生相对滑动；②铸轧区内熔体、铸坯与铸轧辊沿宽度方向传热均匀，轧制力轴向分布均匀，铸轧成形过程简化为平面应变问题；③辊套与辊内冷却水对流换热系数为定值；④铸轧辊与熔体间的换热系数取常数，辊面与带坯间的换热系数与接触压力、温度有关。

1.压下装置　2.平衡弹簧　3.铸轧辊　4.铸轧带材5.浇铸系统　6.铸轧辊　7.振动机构图1　二辊铸轧机结构简图

为减小计算规模，根据传热和变形的对称性，取二分之一结构进行建模，在MSC.Marc软件平台上建立纯铝双辊铸轧过程热-力耦合仿真模型，如图2所示。以铸轧厚度为3 mm的铝带产品为例进行建模分析，辊套厚30 mm，取铸轧区出口中心为坐标原点o，高度方向为y轴，辊面A点(图2)为验证实验温度检测点。图2中，ya为溶池液面高度，yk为单元“杀死”的临界坐标，h2为辊套和冷却水之间的对流换热系数。

图2　铸轧过程热-力耦合模型

1.2传热边界条件

(1)铸轧带坯与辊套之间接触换热边界表达为

qf=h1(tZP-tGTW)

(1)

其中,qf为轧辊与带坯间的热流密度；tZP、tGTW分别为带坯表面温度和辊套表面温度；h1为带坯与辊套之间的接触换热系数，与界面接触压力有关，文献[14]通过界面热阻测试实验，给出了其回归公式：

(2)

(2)辊套与内侧冷却水间的对流换热边界表达为

qw=h2(tGTN-tS)

(3)

式中，qw为辊套与冷却水之间的热流密度；tGTN和tS为辊套内表面温度和冷却水温度。

模型中辊体、冷却水初始温度和环境温度均为30 ℃。h2取14kW/(m2·K)。

(3)辊套与空气间的对流换热和辐射换热等价为综合换热系数，取25W/(m2·K)。

(4)铸轧带坯中心对称面取绝热边界条件。

1.3材料本构模型

工业纯铝的凝固温度区间为614～659 ℃，凝固潜热为3.9414J/g，其热物性参数如表1所示[20]。

表1　纯铝热物性参数

考虑固-液两相混合体和固相纯铝的热变形流变特性差异，引用文献[12]中纯铝的材料本构模型：

(4)

材料本构模型通过MSC.Marc二次开发接口子程序urpflo.f嵌入有限元模型[21]。铸轧辊材料选用42CrMo，相关热物性参数直接从MSC.Marc材料数据库中读取。

1.4生死单元法及连续浇铸建模

由于铸轧辊表面与铸坯间的传热具有间断周期性特点，生产实践表明，铸轧辊达到热平衡至少需要工作10圈以上，铸轧带材数十米。若采用传统带-辊耦合传热建模方法，则入口处需设置数量巨大的“金属熔体”网格，给模型求解带来困难。

如图2所示，本文利用MSC.Marc中的生死单元法(二次开发接口子程序uactive.f)，在铸轧入口浇铸区上方预先划分足够数量的网格，并将节点坐标y大于ya(即熔池液面高度)的单元全部“杀死”。随着铸轧辊的转动，铸轧区内金属网格会将“死单元”带入铸轧区，当其节点坐标y小于ya时，通过子程序将新进入铸轧区的“死单元”重新“激活”。同样，为避免铸轧出口铝带长度过大而导致模型网格数量过多，在出口方向以yk作为单元“杀死”的临界坐标，单元节点坐标y小于yk时将被“杀死”。因此，整个模拟过程中仅节点坐标在yk和ya之间的单元保持“激活”状态，参与模型求解的单元数量大大减少，解决了连续浇铸模拟和快速仿真问题。

2　模拟结果分析

2.1铝带单元“激活”与“杀死”处理

图3是铸轧速度为2.4m/min时，完成第一圈轧制后的辊面温度变化云图。从图3中可以看出，金属熔体入口液面始终保持不变，当铸轧出口带坯延伸至yk点后，单元即被“杀死”。在有限元法中，对于被“杀死”的单元只是将其单元刚度矩阵乘上一个较小的参数(如1.0×10-6)，并非将其从实际模型中删除。由于被“杀死”单元的单元载荷、质量和热边界条件等其他同类参数均为0，故可大大提高模型求解效率。铝液入口单元温度初值则可在模型中设置，不受单元生死操作的影响。

(a)时间τ=0(b)时间τ=4 s

(c)时间τ=12 s(d)时间τ=18 s图3　铸轧过程辊面温度分布变化

模拟获得了铸轧区稳态温度场(图4a)和带坯等效应力分布(图4b)，从图4中可以看出铸轧区内液相区、固-液两相区、固相区和KISS点高度，以及铝带等效应力分布。结合铸轧辊与铝带间接触压力模拟结果(图5)和式(2)，可计算出铸轧区内铝带和铸轧辊间的接触换热系数分布曲线(图5)，接触压力从入口到出口先增大后减小，接触换热系数与接触压力的变化趋势相一致。

(a)温度场(b)等效应力图4　稳态铸轧区(τ=150 s)

图5　铝带与铸轧辊间接触压力和接触换热系数

2.2铸轧速度对KISS点高度和出口温度的影响

取熔池液面高度ya=30 mm，铸轧速度v分别为1.6 m/min、2.4 m/min、3.2 m/min，浇铸温度t0=700 ℃，KISS点高度和铝带出口温度变化规律如图6所示。从图6中可以看出，当v=3.2 m/min时，铝带出口温度处于固液两相区，会发生漏钢事故；当v=1.6 m/min时，出口温度较低，KISS点高度过高，易造成卡壳事故。模拟结果表明，当v=2.4 m/min时，KISS点高度由开始的12.74 mm随时间逐渐降低并稳定于8 mm(图6b)，铸轧过程稳定。

(a)出口温度

(b)KISS点高度(v=2.4 m/min)图6　铸轧速度对KISS点高度和铝带出口温度的影响

2.3浇铸温度对KISS点高度和出口温度的影响

(a)出口温度

(b)KISS点高度图7　浇铸温度对KISS点高度和出口温度的影响

取熔池液面高度ya=30 mm，铸轧速度v=2.4 m/min，浇铸温度t0分别为680 ℃、700 ℃、720 ℃，KISS点高度和铝带出口温度的变化规律如图7所示。从图7中可以看出，三种工况下KISS点高度随铸轧圈数的变化趋势基本相似，浇铸温度每降低20 ℃，铝带出口温度约降低10 ℃(图7a)，相同时刻KISS点高度升高约1 mm(图7b)。

2.4熔池高度对KISS点高度和出口温度的影响

取浇铸温度t0=700 ℃，铸轧速度v=2.4 m/min，熔池高度ya分别为20 mm、30 mm和40 mm，KISS点高度和铝带出口温度变化规律如图8所示。当熔池高度为20 mm时，铝带出口温度超过625 ℃，会出现漏钢事故(图8a)；当熔池高度由30 mm升高至40 mm时，稳态KISS点高度由8 mm提高至19.3 mm，铝带出口温度由570 ℃降低至498 ℃。

2.5辊套热平衡分析及实验验证

根据上述模拟结果，确定适于本实验铸轧机的合理工艺参数为：熔池高度30 mm、浇铸温度700 ℃、铸轧速度2.4 m/min，并以以上参数为基础进行辊套热平衡分析，模拟获得的辊面及辊套径向温度随时间变化曲线如图9所示，图中，h′为测点距辊面的距离。

(a)出口温度

(b)KISS点高度图8　熔池高度对KISS点高度和铝带出口温度的影响

(a)辊面温度

1.h′=5 mm　2.h′=10 mm　3.h′=15 mm4.h′=20 mm　5.h′=25 mm　6.h′=30 mm(b)辊套径向温度图9　辊套热平衡模拟结果

铸轧过程辊套内各点温度随时间呈周期性变化，距辊面越近，温度波动周期性越明显(图9b)。轧制9圈后辊套温度变化达到稳定状态，辊面最高温度为298 ℃，最低温度为125 ℃(图9a)，辊套内壁温度为45 ℃(图9b)。

为验证模型精度，利用Raytek红外测温仪对相同工况下(图10a)辊面温度进行采集。由于工作过程轧辊处于转动状态，无法对辊面进行定点测温，本文通过事先在辊面中心A点(图2)对应的辊身端部贴签做标识，分别连续铸轧1圈、2圈、3圈后停机对A点进行温度测试，结果表明，实测值与有限元模拟结果基本吻合(图10b)。

(a)实验现场

(b)辊面温度对比图10　温度模拟结果验证

3　结束语

本文较为系统地分析了双辊铸轧过程热-力耦合有限元建模的核心数学模型，并采用生死单元方法，解决了金属熔体连续浇铸建模问题，为铸轧辊热平衡分析和工艺优化提供了新途径。以二辊φ160 mm×150 mm立式实验铸轧机为对象，模拟分析了铸轧速度、熔池高度、浇铸温度对铸轧区KISS点高度和带材出口温度的影响规律，给出了合理的铸轧工艺规程及相应工况下的辊套热平衡演变过程，并通过铸轧实验对温度模拟结果进行了验证。

[1]赵红阳, 胡林, 李娜. 双辊薄带铸轧技术的进展及热点问题评述[J]. 鞍钢技术, 2007(6): 1-5.

Zhao Hongyang, Hu Lin, Li Na. Development of Two-high Thin Strip Casting and Rolling Technology and Key Point Review[J].Angang Technology,2007(6):1-5.

[2]Zhao Hu, Li Peijie, He Liangju. Coupled Analysis of Temperature and Flow During Twin-roll Casting of Magnesium Alloy Strip[J]. Journal of Materials Processing Technology, 2011, 211(6): 1197-1202.

[3]Wang Bo, Zhang Jieyu, Li Xiangmei, et al. Simulation of Solidification Microstructure in Twin-roll Casting Strip[J]. Computational Materials Science, 2010, 49(1):S135-S139.

[4]陈守东, 陈敬超. 双辊薄带连铸凝固组织模拟微观模型的验证[J]. 材料热处理学报, 2012, 33(7): 154-158.

Chen Shoudong, Chen Jingchao. Validation of Micro Model for Solidification Process of Twin-roll Continuous Casting Thin Strip[J]. Transactions of Materials and Heat Treatment, 2012,33(7):154-158.

[5]湛利华, 李晓谦, 唐朝阳, 等. 铝带坯连续铸轧过程热力耦合有限元分析[J]. 中国机械工程, 2005, 16(11): 979-984, 996.

Zhan Lihua, Li Xiaoqian, Tang Zhaoyang,et al. Thermo-mechanical Coupling Analysis of Aluminum Continuous Roll Casting Processes[J]. China Mechanical Engineering,2005, 16(11): 979-984, 996.

[6]Parka C M, Choia J T.Thermal Crown Analysis of the Roll in the Strip Casting Process[J]. Journal of Materials Processing Technology， 2009, 209(8):3714-3723.

[7]朱志华, 肖文锋, 李晓谦. 双辊连续铸轧过程中轧制界面接触压力切块法建模与数值模拟[J]. 中国机械工程, 2005, 16(13): 1091-1095.

Zhu Zhihua, Xiao Wenfeng, Li Xiaoqian. Modeling of Rolling Pressure Using Slab-method and the Numerical Simulation during Twin-roll Continuous Roll Casting Process[J]. China Mechanical Engineering,2005, 16(13): 1091-1095.

[8]Jerry W S, Robert J, Comstock J R. Method of Continuous Casting Non-oriented Electrical Steel Strip:US, 2006015z142A1[P].2006-07-13.

[9]许志强, 孟哲儒, 杜凤山, 等. 双辊薄带铸轧数值模拟研究现状及展望[J]. 燕山大学学报, 2014,39(2):95-101.

Xu Zhiqiang, Meng Zheru, Du Fengshan, et al. Current Situation and Prospect of Twin-roll Strip Casting Process Numerical Simulation[J]. Journal of Yanshan University, 2014,39(2):95-101.

[10]谢金乐, 刘允中, 吴汇江, 等. 半固态7050铝合金热压缩变形行为[J]. 特种铸造及有色合金, 2011, 31(9): 816-819.

Xie Jinle, Liu Yunzhong, Wu Huijiang, et al. Hot Compression Behavior of Semi-solid 7050 Aluminum Alloy[J]. Special Casting & Nonferrous Alloys, 2011, 31(9): 816-819.

[11]湛利华, 钟掘, 李晓谦, 等. 连续铸轧流变行为的物理模拟及其应力-应变关系的演变[J]. 中国有色金属学报,2004,14(12):1995-2002.

Zhan Lihua, Zhong Jue, Li Xiaoqian, et al. Physical Simulation of Rheological Behavior and Stress-strain Relation Evolvement in Continuous Roll Casting Process[J]. The Chinese Journal of Nonferrous Metals, 2004,14(12): 1995-2002.

[12]湛利华. 铝合金连续铸轧过程流变行为研究及热-力耦合分析[D]. 长沙：中南大学, 2005.

[13]刘菊. 固体界面接触热阻及热导率测量的实验研究[D]. 武汉：华中科技大学, 2011.

[14]张云湘. 瞬态法接触热阻实验研究及其在快速铸轧工艺参数仿真中的应用[D]. 长沙：中南大学, 2004.

[15]Wang De, Zhou Cheng, Xu Guojin, et al. Heat Transfer Behavior of Top Side-pouring Twin-roll Casting[J]. Journal of Materials Processing Technology, 2014, 214(6): 1275-1284.

[16]王晓晨, 杨荃, 刘瑞军, 等. 基于ANSYS有限元法的热卷箱内中间坯温度场分析[J]. 北京科技大学学报, 2013(4):454-458.

Wang Xiaochen, Yang Quan, Liu Ruijun, et al. Temperature Field Analysis of Intermediate Slabs in the Hot Coil Box Based on ANSYS Finite Element Method[J]. Journal of University of Science and Technology Beijing, 2013(4): 454-458.

[17]赵洪运, 舒凤远, 张洪涛, 等. 基于生死单元的激光熔覆温度场数值模拟[J]. 焊接学报, 2010, 31(5): 81-84.

Zhao Hongyun, Shu Fengyuan, Zhang Hongtao, et al. Numerical Simulationon Temperature Field of Laser Cladding Based on Birth-death Element Method[J]. Transactions of the China Welding Institution,2010, 31(5): 81-84.

[18]郑江, 葛鸿鹏, 王先铁, 等. 局部位形约束生死单元法及其在施工力学分析中的应用[J]. 建筑结构学报,2012, 33(8): 101-108.

Zheng Jiang, Ge Hongpeng, Wang Xiantie, et al. A Method of Element Birth and Death of Local Configuration Con-straint and Its Application in Construction Mechanics[J]. Journal of Building Structures,2012, 33(8): 101-108.

[19]李尧. 双辊薄带振动铸轧过程仿真模拟及实验研究[D]. 秦皇岛：燕山大学, 2011.

[20]王祝堂. 铝合金及其加工手册[M]．3版. 长沙：中南大学出版社，2005．

[21]陈火红, 尹伟奇, 薛小香. MSC.Marc二次开发指南[M]. 北京:科学出版社, 2004.

(编辑袁兴玲)

Thermal-mechanical Coupled Modelling and Numerical Simulation for Twin-roller Casting Process with Technique of Deactivate and Reactivate Element

Huang Huagui1,2Liu Wenwen1,2Wang Wei2Du Fengshan1,2

1.National Engineering Research Center for Equipment and Technology of Cold Strip Rolling,Qinhuangdao,Hebei,066004 2.Yanshan University,Qinhuangdao,Hebei,066004

Aiming at the φ160 mm×150 mm vertical type experimental twin-rollers continuous caster, a thermal-mechanical coupled FEM model of twin-roller casting process was established with MSC.Marc and its’ secondary development interface based on the solid-liquid two phases constitutive model of pure aluminum and the thermal resistance model considering the contact pressure. The technique of deactivate and reactivate element method was used to simulate the continuous casting of aluminums liquid, and the problem of coupled heat transmission between the cast roll shell and aluminum strip in the cast rolling area was therefore resolved. From the simulation results, the influences of casting speed, casting temperature, the height of cast rolling area on the position of KISS point and the time history of temperature on the roll surface were presented. An twin-roller casting experiment was carried out, and the numerical simulation results are in good agreements with the experimental data.

twin-roller casting; temperature field; finite element method(FEM); deactivate and reactivate element method; numerical simulation

2014-06-30

国家自然科学基金资助项目(51474189)；河北省自然科学基金资助项目(E2013203377)

TG233.6;TP391.9DOI：10.3969/j.issn.1004-132X.2015.11.013

黄华贵，男，1978年生。燕山大学机械工程学院教授、博士。主要研究方向为材料加工工艺及装备。获国家发明专利6项，发表论文30余篇。刘文文，男，1990年生。燕山大学机械工程学院硕士研究生。王巍，女，1978年生。燕山大学机械工程学院讲师。杜凤山，男，1960年生。燕山大学机械工程学院教授、博士研究生导师。