小方坯连铸过程温度场和流场的实时模拟

第36卷第12期 北京科技大学学报 Vol.36 No.12 2014年12月 Journal of University of Science and Technology Beijing Dec.2014 小方坯连铸过程温度场和流场的实时模拟 唐娜娜，李宏祥⑧，张济山，庄林忠 北京科技大学新金属材料国家重点实验室，北京100083 ☒通信作者，E-mail:hxli@skL.usth.cdu.cm 摘要为了寻求在瞬态下仍然有效的实时模拟模型，以便模拟连铸的全过程，基于任意拉格朗日一欧拉(ALE)方法和Lam一 Bremhorst低雷诺数ke方程，建立了考虑湍流的传热和流体流动的实时模拟模型.该模型有效地实现了移动计算域下的动态 网格扩展，可以模拟从拉坯到切割坯连铸全过程的温度场与流场的变化.作为实例，X70钢小方坯的模拟效果与实验验证吻 合很好 关键词连铸：小方坯：温度场：流场：有限元法：实时模拟 分类号TG249.7 Real-time simulation of temperature field and flow field for a billet in continuous casting TANG Na-na,LI Hong-xiang,ZHANG Ji-shan,ZHUANG Lin-zhong State Key Laboratory for Advanced Metals and Materials,University of Science and Technology Beijing.Beijing 100083,China Corresponding author,E-mail:hxli@skl.ustb.edu.cn ABSTRACT To develop a real-time mathematical model which can be valid ever at a transient state for simulating full continuous casting processes,a real-time simulation mathematical model of heat transfer and fluid flow was constructed considering turbulent flow on the basis of the arbitrary Lagrangian-Eulerian (ALE)algorithm and the Lam-Bremhorst low-Reynolds number formulation.By using this model,dynamic grids expansion under moving computational domains could be realized and the temperature field and fluid flow field could be simulated for the full process of steel continuous casting.As an example,the simulation results of X70 steel are in accordance with experimental measurements. KEY WORDS continuous casting:billets;temperature field;flow field:finite element method;real-time simulation 因为成本低、能耗低且效率高，连铸己经成为生 优势，因此国内外钢铁研究的学者已经开展了很多 产方坯和板坯主要的手段：但其技术复杂，合理控制 工作1.Liu等四比较了紧凑式带钢生产(CSP)和 方坯或者板坯的温度场和流场分布无疑是关注的重 可变薄板坯连铸TSC两种工艺下结晶器中的流场 点之一.不合理的温度分布与流体流动，不但影响 与温度场，发现水口结构对流动模式、熔体过热量的 钢水的凝固和夹杂物的去除，诱发鼓肚缺陷，而且与 耗散以及凝固坯壳的均匀生长影响比较大.Zhu 铸坯裂纹、偏析等表面及内部质量密切相关，严重时 等回研究了铸造速度、水口设计和水口插入深度对 还易导致漏钢等事故的发生.计算机数值模拟作为 漏斗形结晶器流场和温度场的影响，为连铸工艺优 经济、省力且高效的研究手段，在优化温度场和流场 化提供了依据.Zhang等同则通过建立二维有限点 分布，寻求最佳连铸工艺参数等方面展示出明显的 无孔模型模拟了连铸结晶器内的传热和凝固.除了 收稿日期：201309-17 基金项目：中央高校基本科研业务费专项资金资助项目(FRF-TD-12001):教有部博士学科点专项科研基金资助项目(20120006110019)：新 金属材料国家重点实验室开放课题(2012Z-13) DOI:10.13374/j.issnl1001-053x.2014.12.010:http://journals.ustb.edu.cn

第 36 卷 第 12 期 2014 年 12 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 36 No． 12 Dec． 2014 小方坯连铸过程温度场和流场的实时模拟 唐娜娜，李宏祥，张济山，庄林忠 北京科技大学新金属材料国家重点实验室，北京 100083  通信作者，E-mail: hxli@ skl． ustb． edu． cn 摘 要 为了寻求在瞬态下仍然有效的实时模拟模型，以便模拟连铸的全过程，基于任意拉格朗日--欧拉( ALE) 方法和Lam-- Bremhorst 低雷诺数 k-ε 方程，建立了考虑湍流的传热和流体流动的实时模拟模型． 该模型有效地实现了移动计算域下的动态 网格扩展，可以模拟从拉坯到切割坯连铸全过程的温度场与流场的变化． 作为实例，X70 钢小方坯的模拟效果与实验验证吻 合很好． 关键词 连铸; 小方坯; 温度场; 流场; 有限元法; 实时模拟 分类号 TG249. 7 Ｒeal-time simulation of temperature field and flow field for a billet in continuous casting TANG Na-na，LI Hong-xiang ，ZHANG Ji-shan，ZHUANG Lin-zhong State Key Laboratory for Advanced Metals and Materials，University of Science and Technology Beijing，Beijing 100083，China  Corresponding author，E-mail: hxli@ skl． ustb． edu． cn ABSTＲACT To develop a real-time mathematical model which can be valid ever at a transient state for simulating full continuous casting processes，a real-time simulation mathematical model of heat transfer and fluid flow was constructed considering turbulent flow on the basis of the arbitrary Lagrangian--Eulerian ( ALE) algorithm and the Lam--Bremhorst low-Ｒeynolds number k-ε formulation． By using this model，dynamic grids expansion under moving computational domains could be realized and the temperature field and fluid flow field could be simulated for the full process of steel continuous casting． As an example，the simulation results of X70 steel are in accordance with experimental measurements． KEY WOＲDS continuous casting; billets; temperature field; flow field; finite element method; real-time simulation 收稿日期: 2013--09--17 基金项目: 中央高校基本科研业务费专项资金资助项目( FＲF--TD--12--001) ; 教育部博士学科点专项科研基金资助项目( 20120006110019) ; 新 金属材料国家重点实验室开放课题( 2012Z--13) DOI: 10． 13374 /j． issn1001--053x． 2014． 12． 010; http: / /journals． ustb． edu． cn 因为成本低、能耗低且效率高，连铸已经成为生 产方坯和板坯主要的手段; 但其技术复杂，合理控制 方坯或者板坯的温度场和流场分布无疑是关注的重 点之一． 不合理的温度分布与流体流动，不但影响 钢水的凝固和夹杂物的去除，诱发鼓肚缺陷，而且与 铸坯裂纹、偏析等表面及内部质量密切相关，严重时 还易导致漏钢等事故的发生． 计算机数值模拟作为 经济、省力且高效的研究手段，在优化温度场和流场 分布，寻求最佳连铸工艺参数等方面展示出明显的 优势，因此国内外钢铁研究的学者已经开展了很多 工作［1--6］． Liu 等［1］比较了紧凑式带钢生产( CSP) 和 可变薄板坯连铸 FTSC 两种工艺下结晶器中的流场 与温度场，发现水口结构对流动模式、熔体过热量的 耗散以及凝固坯壳的均匀生长影响比较大． Zhu 等［2］研究了铸造速度、水口设计和水口插入深度对 漏斗形结晶器流场和温度场的影响，为连铸工艺优 化提供了依据． Zhang 等［3］则通过建立二维有限点 无孔模型模拟了连铸结晶器内的传热和凝固． 除了

第12期 唐娜娜等：小方坯连铸过程温度场和流场的实时模拟 ·1635· 结晶器，部分研究者也研究了二冷区或者连铸 u:=b1:+b,u (1) 全过程温度的变化.常运合等0通过开发异型坯连 式中，b,和b.分别是液相体积分数和固相体积分 铸动态二冷控制模型以铸坯温度场信息为基础对二 数，：和u分别是液相和固相的体积平均速度.固 冷水量进行了设定和优化.王迎春等因对连铸全过 相速度是恒定的，由下式计算得到： 程进行了热状态、凝固状态分析及工艺控制，并根据 u;=-V.6a (2) 连铸治金准则和目标温度控制进行二冷优化，获得 式中，V是拉坯速度，δ是表示方向的矢量.达西力 合理的温度分布和二冷水量分布.Ludwig等a模拟 计算公式如下： 了大方坯连铸全过程温度和流动行为.然而，综合 0 T>T, 分析国内外的研究状况可以看出，钢铁连铸温度场 和流场的模拟多采用商业化软件如Asys、CFX、 S= T<T<T, (3) Fluent、MSC.MARC或者自编程序等基于稳态连铸 M(u-u:） T<T 过程而展开，只是对固定计算域内的静态几何进行 式中，T、T,和T分别是温度、固相线温度和液相线温 研究，很少考虑到瞬态过程.稳态数值模拟的缺点 度，v和K分别是钢液的运动黏度和糊状区的渗透 是无法得知热和力的演化过程，不利于应力与变形 率，M是一个可使固体区域的混合流场达到固体速 的计算.为此人们不得不考虑建立瞬态下仍然有效 度的较大的数.渗透率K与糊状区的形状有关，由 的实时模拟模型，通过该模型可以模拟钢水从充满 Kozeny-Carman方程计算得到： 结晶器到拉坯再到板坯或者方坯切割全过程的温度 b 场和流场的变化.目前这方面的模拟研究开展得很 K=K1-b)2 (4) 少，明显的进展是ProCAST软件中采用Mile算法可 式中，K,是与二次枝晶臂间距（山2）有关的参数，其 以实现移动计算域下网格的移动扩展从而可以模拟 温度场与流场的变化P-),但Mle算法无法实现连 计算如下☒ 铸弧形段的瞬态计算回.基于任意拉格朗日一欧拉 d K=180 (5) (ALE)方法和Lam-Bremhorst低雷诺数ke(k为湍 浮力计算如下： 动能，ε为湍流耗散)方程，本文尝试建立考虑湍流 的传热和流体流动的实时模拟模型，并有效地实现 F=Ap Po 8. (6) 了移动计算域下的动态网格扩展.本文应用该模型 式中，△p是由于任意参考温度而引起的密度差，参 模拟了某钢厂X70钢小方坯直弧形连铸全过程温 考温度通常是指固相线温度或液相线温度.P和g 度场和流场的变化，模拟的凝固坯壳值与工厂实测 分别是密度常数和重力加速度.由此，质量守恒方 值比较接近，模拟效果良好. 程和动量守恒方程可写成下面的形式： 1数学模型 ai=0, (7) 8xi 传热与流体流动的模拟基于Bennon和Incrop-- ou; era推导的连续介质（混合物）模型而展开0，假设 dt 0=-12+2+) po ax ax) 0x 整个糊状区为枝晶网，忽略自由移动的固体晶粒的 (8) 影响.考虑到固相和液相之间速度不同而引起的界 g+8-号器 面摩擦，糊状区流场计算在混合动量方程中引入了 式中，P、山，和K分别是压力、湍流运动黏度和湍动 达西力S(i为x、y和z三个方向).模型忽略了凝 能.式(8)中最后一项为雷诺应力项.湍流由浇口 固收缩和由于合金成分差异引起的溶质传输（即宏 或者受浮力影响的流动所引起.湍流的程度也是非 观偏析)，热对流则采用Boussinesq近似，引入了浮 常重要的，不但影响流场的计算，而且对于由于湍流 力F.湍流采用Launder和Sharma开发的低雷诺 导致的热传输的计算也有影响.目前的模型中糊状 数(LRN)Ke模型. 区湍流的衰减通过在湍流方程中引入达西项来完 1.1流体流动 成.湍流参数K和E从LRNK一e传输方程中计算 液体溶池中的流动由浇口的强制性对流以及热 得到： 诱导的浮力（自由对流或自然对流）所引起.混合速 度的定义如下： +器品+)+

第 12 期 唐娜娜等: 小方坯连铸过程温度场和流场的实时模拟 结晶器，部分研究者［4--6］也研究了二冷区或者连铸 全过程温度的变化． 常运合等［4］通过开发异型坯连 铸动态二冷控制模型以铸坯温度场信息为基础对二 冷水量进行了设定和优化． 王迎春等［5］对连铸全过 程进行了热状态、凝固状态分析及工艺控制，并根据 连铸冶金准则和目标温度控制进行二冷优化，获得 合理的温度分布和二冷水量分布． Ludwig 等［6］模拟 了大方坯连铸全过程温度和流动行为． 然而，综合 分析国内外的研究状况可以看出，钢铁连铸温度场 和流场的模拟多采用商业化软件如 Ansys、CFX、 Fluent、MSC． MAＲC 或者自编程序等基于稳态连铸 过程而展开，只是对固定计算域内的静态几何进行 研究，很少考虑到瞬态过程． 稳态数值模拟的缺点 是无法得知热和力的演化过程，不利于应力与变形 的计算． 为此人们不得不考虑建立瞬态下仍然有效 的实时模拟模型，通过该模型可以模拟钢水从充满 结晶器到拉坯再到板坯或者方坯切割全过程的温度 场和流场的变化． 目前这方面的模拟研究开展得很 少，明显的进展是 ProCAST 软件中采用 Mile 算法可 以实现移动计算域下网格的移动扩展从而可以模拟 温度场与流场的变化［7--8］，但 Mile 算法无法实现连 铸弧形段的瞬态计算［9］． 基于任意拉格朗日--欧拉 ( ALE) 方法和 Lam-Bremhorst 低雷诺数 k-ε ( k 为湍 动能，ε 为湍流耗散) 方程，本文尝试建立考虑湍流 的传热和流体流动的实时模拟模型，并有效地实现 了移动计算域下的动态网格扩展． 本文应用该模型 模拟了某钢厂 X70 钢小方坯直弧形连铸全过程温 度场和流场的变化，模拟的凝固坯壳值与工厂实测 值比较接近，模拟效果良好． 1 数学模型 传热与流体流动的模拟基于 Bennon 和 Incrop￾era 推导的连续介质( 混合物) 模型而展开［10］，假设 整个糊状区为枝晶网，忽略自由移动的固体晶粒的 影响． 考虑到固相和液相之间速度不同而引起的界 面摩擦，糊状区流场计算在混合动量方程中引入了 达西力 Si ( i 为 x、y 和 z 三个方向) ． 模型忽略了凝 固收缩和由于合金成分差异引起的溶质传输( 即宏 观偏析) ，热对流则采用 Boussinesq 近似，引入了浮 力 Fi ． 湍流采用 Launder 和 Sharma［11］开发的低雷诺 数( LＲN) κ-ε 模型． 1. 1 流体流动 液体溶池中的流动由浇口的强制性对流以及热 诱导的浮力( 自由对流或自然对流) 所引起． 混合速 度的定义如下: ui = blul i + bsus i ． ( 1) 式中，bl 和 bs 分别是液相体积分数和固相体积分 数，ul i 和 us i 分别是液相和固相的体积平均速度． 固 相速度是恒定的，由下式计算得到: us i = － Vc δi3 ( 2) 式中，Vc是拉坯速度，δi3是表示方向的矢量． 达西力 计算公式如下: Si = 0 T ＞ Tl， ν K ( us i － ui ) Ts ＜ T ＜ Tl， M( us i － ui ) T ＜ Ts      ． ( 3) 式中，T、Ts和 Tl分别是温度、固相线温度和液相线温 度，ν 和 K 分别是钢液的运动黏度和糊状区的渗透 率，M 是一个可使固体区域的混合流场达到固体速 度的较大的数． 渗透率 K 与糊状区的形状有关，由 Kozeny--Carman 方程计算得到［12］: K = K0 b 3 l ( 1 － bl ) 2 ． ( 4) 式中，K0 是与二次枝晶臂间距( d2 ) 有关的参数，其 计算如下［12］: K0 = d2 2 180． ( 5) 浮力计算如下: Fi = Δρ ρ0 g． ( 6) 式中，Δρ 是由于任意参考温度而引起的密度差，参 考温度通常是指固相线温度或液相线温度． ρ0 和 g 分别是密度常数和重力加速度． 由此，质量守恒方 程和动量守恒方程可写成下面的形式: ui xi = 0， ( 7) ui t + uj ui xj = － 1 ρ0 p xi +  xj ( ν + νt ) ui xj + Fi + Si － 2 3 κ xi ． ( 8) 式中，p、νt 和 κ 分别是压力、湍流运动黏度和湍动 能． 式( 8) 中最后一项为雷诺应力项． 湍流由浇口 或者受浮力影响的流动所引起． 湍流的程度也是非 常重要的，不但影响流场的计算，而且对于由于湍流 导致的热传输的计算也有影响． 目前的模型中糊状 区湍流的衰减通过在湍流方程中引入达西项来完 成． 湍流参数 κ 和 ε 从 LＲN κ--ε 传输方程中计算 得到: κ t + uj κ xj =  x (j ν + νt C ) κ κ xj + · 5361 ·

·1636 北京科技大学学报 第36卷 p+G--2-2 K 数，表明温度边界层与流动边界层的关系，根据经验 (9) 取为1.0.从式(13)中的第二项可以看出，潜热以 "1e+ E+山=证"+C/远 固相速度释放，此处忽略了糊状区自由凝固的钢液， 假设整个固相分数与以拉坯速度运动的树枝晶组织 C,P+(1-C)G月s-CfE2 有关 K 1.3网格扩展几何 子4：12 2e+2m(,） (10) 模拟连铸的拉坯过程，要求能够实现移动计算 域下网格随着拉坯的进行而动态进行扩展.单纯的 式中，C1、C2、C3、C.和C,是常数，壁面边界条件是 Eulerian方法，网格保持固定，表达自由边界比较困 ￡=K=0.式(9)中的最后一项来自于对湍流耗散参 难，而且无法适应计算域是移动边界条件下的计算 数的重新定义.壁面处耗散是零，额外项将逐渐接 单纯的Lagrangian算法尽管相比单纯的Eulerian算 近真实耗散的非零值.P和G是湍流产生项. 法由于缺少对流项，控制方程简单，但是由于网格点 与物质点重合，物质的扭曲容易造成计算网格的畸 形，从而导致计算失败，同时单纯的Lagrangian算法 v.dfi G二pi,ax: ”=Cf (11) 无法处理一些接触边界条件，如拐角或者某些尖锐 的边.在连铸过程中部分计算域是固定的（如浸入 式中，C是常数，P,是湍流普朗特数.LRN阻尼函 式浇口)，部分计算域则是以拉坯速度运动的（如板 数定义如下： 坯随着引锭头而移动)，所以使用单纯的Eulerian算 34 f=1.0-03e=ie, 法和单纯的Lagrangian算法都无法准确表达钢坯的 es2 拉坯行为.例如浇口等在实际连铸过程中位置固 (12) vE 定，如采用单纯的Lagrangian算法处理移动网格，则 需要注意的是式(11)中的G,当冷的液体上部 容易引起整个计算域的速度波动.因此，引入了任 为热的液体时，需要取负值，因为这种情况下抑制了 意拉格朗日-欧拉(ALE)方法很好地解决了连铸过 动量的垂直湍流流动.在式(I2)中f,为Launder和 程中网格固定以及移动扩展的问题.该方法计算网 Sharma定义的函数，表示了与局部湍流雷诺数有关 格可以在空间中以任意的形式运动，即独立于物质 的单项流的阻尼机制.f可以用函数修改，表示 坐标系和空间坐标系运动，通过规定合适的网格运 了糊状区液相分数变化的影响，这与Shy等的研究 动形式可以准确地描述物体的移动界面，并维持合 相似，取α=1/2国.式(9)和(10)的达西阻尼项只 理的单元形状.纯拉格朗日和纯欧拉方法实际上是 有在计算糊状区时才予以考虑，这时LRN阻尼函数 ALE描述的两个特例，当网格的运动速度等于物体 f将取Launder和Sharma原始定义的数值（即不考 的运动速度时就退化为拉格朗日描述，当网格固定 虑b“),而Launder等计算使用的经验常数值己经成 于空间不动时就退化为欧拉描述.有关任意拉格朗 为K传输方程的标准值，这些经验值具体可参考 日-欧拉方法方法的详尽描述，参考文献14].ALE 文献01]. 方法其基本原理是：引入一个参考构型，该构型由一 1.2热传输 组无约束网格点所组成，这些网格点可以在空间任 考虑到由于瑞流混合而引起的热扩散的增加， 意运动.参考构型中的每一点都由三个独立坐标X: 混合的能量守恒方程也需要加入湍流混合项来进行 表示，这样参考构型中每点的运动都可以用空间坐 修改，改后的方程如下： 标x来表达，其是K和时间t的函数，即 aH (14) 识+p叫+p%-g xi=x;(,t). ax; 参考构型中的坐标X、空间坐标x:和物质坐标X:之 品a+院)品 间的关系可以通过物质点和空间点在参考构型中的 (13) ax; 匹配建立起来，如下所示： 式中，P、t、H、T、c。、入和Pr,分别是密度、时间、混合 xi=f(xj,t)=f (X,t) (15) 焓、温度、比热容、热导率和湍流普朗特数.能量方 式中，f和：是单值连续一阶导数.由此，在ALE描 程中忽略了黏滞热.湍流普朗特数表征流体流动中 述中，守恒方程用网格计算节点的速度（此处等于 动量交换与热交换相对重要性的一个量纲一的参 固体速度u)来修改.Eulerian-agrangian形式的质

北 京 科 技 大 学 学 报 第 36 卷 P + G － ε － 2 ν K κ － 2ν ( 槡κ x ) j 2 ， ( 9) ε t + uj ε xj =  x (j ν + νt C ) ε ε xj + C1［P + ( 1 － C3 ) G］ε － C2 fε ε2 κ － 2 ν K ε + 2ννt (  2 ui xj x ) k 2 ． ( 10) 式中，C1、C2、C3、Cκ 和 Cε 是常数，壁面边界条件是 ε = κ = 0． 式( 9) 中的最后一项来自于对湍流耗散参 数的重新定义． 壁面处耗散是零，额外项将逐渐接 近真实耗散的非零值． P 和 G 是湍流产生项． P = νt ( ui xj + uj x ) i ui xj ， G = νt Prt fi xi ，νt = Cvfv κ2 ε ． ( 11) 式中，Cv是常数，Prt是湍流普朗特数． LＲN 阻尼函 数定义如下: fε = 1. 0 － 0. 3e － Ｒe2 ，fv = gα l e ( － 3. 4 1 + Ｒet 50 ) 2 ， Ｒe = κ2 νε ． ( 12) 需要注意的是式( 11) 中的 G，当冷的液体上部 为热的液体时，需要取负值，因为这种情况下抑制了 动量的垂直湍流流动． 在式( 12) 中 fv 为 Launder 和 Sharma 定义的函数，表示了与局部湍流雷诺数有关 的单项流的阻尼机制． fv 可以用函数 bα l 修改，表示 了糊状区液相分数变化的影响，这与 Shyy 等的研究 相似，取 α = 1 /2 ［13］． 式( 9) 和( 10) 的达西阻尼项只 有在计算糊状区时才予以考虑，这时 LＲN 阻尼函数 fv 将取 Launder 和 Sharma 原始定义的数值( 即不考 虑 bα l ) ，而 Launder 等计算使用的经验常数值已经成 为 κ-ε 传输方程的标准值，这些经验值具体可参考 文献［11］． 1. 2 热传输 考虑到由于湍流混合而引起的热扩散的增加， 混合的能量守恒方程也需要加入湍流混合项来进行 修改，改后的方程如下: ρ H t + ρus j H xj + ρ( uj － us j ) cpT xj =  x (j λ + ρcpνt Pr ) t T xj ． ( 13) 式中，ρ、t、H、T、cp、λ 和 Prt分别是密度、时间、混合 焓、温度、比热容、热导率和湍流普朗特数． 能量方 程中忽略了黏滞热． 湍流普朗特数表征流体流动中 动量交换与热交换相对重要性的一个量纲一的参 数，表明温度边界层与流动边界层的关系，根据经验 取为 1. 0． 从式( 13) 中的第二项可以看出，潜热以 固相速度释放，此处忽略了糊状区自由凝固的钢液， 假设整个固相分数与以拉坯速度运动的树枝晶组织 有关． 1. 3 网格扩展几何 模拟连铸的拉坯过程，要求能够实现移动计算 域下网格随着拉坯的进行而动态进行扩展． 单纯的 Eulerian 方法，网格保持固定，表达自由边界比较困 难，而且无法适应计算域是移动边界条件下的计算． 单纯的 Lagrangian 算法尽管相比单纯的 Eulerian 算 法由于缺少对流项，控制方程简单，但是由于网格点 与物质点重合，物质的扭曲容易造成计算网格的畸 形，从而导致计算失败，同时单纯的 Lagrangian 算法 无法处理一些接触边界条件，如拐角或者某些尖锐 的边． 在连铸过程中部分计算域是固定的( 如浸入 式浇口) ，部分计算域则是以拉坯速度运动的( 如板 坯随着引锭头而移动) ，所以使用单纯的 Eulerian 算 法和单纯的 Lagrangian 算法都无法准确表达钢坯的 拉坯行为． 例如浇口等在实际连铸过程中位置固 定，如采用单纯的 Lagrangian 算法处理移动网格，则 容易引起整个计算域的速度波动． 因此，引入了任 意拉格朗日--欧拉( ALE) 方法很好地解决了连铸过 程中网格固定以及移动扩展的问题． 该方法计算网 格可以在空间中以任意的形式运动，即独立于物质 坐标系和空间坐标系运动，通过规定合适的网格运 动形式可以准确地描述物体的移动界面，并维持合 理的单元形状． 纯拉格朗日和纯欧拉方法实际上是 ALE 描述的两个特例，当网格的运动速度等于物体 的运动速度时就退化为拉格朗日描述，当网格固定 于空间不动时就退化为欧拉描述． 有关任意拉格朗 日--欧拉方法方法的详尽描述，参考文献［14］． ALE 方法其基本原理是: 引入一个参考构型，该构型由一 组无约束网格点所组成，这些网格点可以在空间任 意运动． 参考构型中的每一点都由三个独立坐标 χi 表示，这样参考构型中每点的运动都可以用空间坐 标 xi来表达，其是 χi和时间 t 的函数，即 xi = xi ( χj ，t) ． ( 14) 参考构型中的坐标 χi、空间坐标 xi和物质坐标 Xi之 间的关系可以通过物质点和空间点在参考构型中的 匹配建立起来，如下所示: χi = fi ( xj ，t) = ^ fi ( Xk，t) ( 15) 式中，fi 和 ^ fi 是单值连续一阶导数． 由此，在 ALE 描 述中，守恒方程用网格计算节点的速度( 此处等于 固体速度 us j ) 来修改． Eulerian-Lagrangian 形式的质 · 6361 ·

第12期 唐娜娜等：小方坯连铸过程温度场和流场的实时模拟 ·1637· 量、动量和能量守恒方程转换为： 隐式的（仅仅发散项）.守恒方程(7)和(8)的半隐 a=0, 式方法如下.考虑到在时间步长数n的时候速度场 (16) ax: 为，压力场为p”、u+1和p+1满足： 测+(y,- ai u* =0， (19) dx; ax; -L2+是(+w)驰+F,+S-?s,(1) Po ox:dxj 3x: 山1=-△1心 uiap p dx; p兴+p%-)=（A+当g aH dx;dxi △日(u+w)3 aug dx; +AF7+A2k 30x: (18) (20) K和ε从式(9)和(10)计算得到，其用网格计算节 定义初步或中间可压缩速度场如下： 点的速度修改与上述动量方程类似.在引锭头和板 坯的已经凝固的部分，Eulerian-Lagrangian方法简化 u;=u-△t"u +a…(+n) ou 0x dx; 为单纯的Lagrangian方法(u=),因此在式(18) 中没有潜热的对流释放.在浸入式浇口部分， 4F-42k (21) 3 ax Eulerian-Lagrangian方法简化为单纯的Eulerian方 通过Hodge分解投射到自由发散场的公式如下（达 法(u=0).由此，方坯或者板坯的计算域划分为三 西法则的源项S与多孔媒介的标量梯度类似)： 个部分，包括钢水入流的浇口在内的方坯或者板坯 .1dp"+ u+1=u4-△r -+△t…Sg+1. (22) 的顶部固定为Euler层，Euler层下边分别为 p dx: Expansion层和Lagrange层.在Expansion层，一层单 在Perot D后，这个方法中的速度场u*1有一 元不断地垂直增长，当它的厚度达到某一临界值时， 个误差项等于△2(u+,)2p"+'p.对于连铸中 单元被分割，下边的分割单元并入Lagrange层并随 的黏度，这个误差项可以忽略，但是如果黏度的数量 引锭头以拉坯速度向下移动.当Lagrange层移动进 级是非常高的（例如金属的挤压），迭代过程是必须 入弧形段时，网格计算节点的速度将用相应的角速 的.压力刚度方程可由质量守恒方程(19)、映射方 度表示，角速度等于拉坯速度与弧形段对应半径的 程(22)得到： 比值.通过上述ALE算法，将能够实现包括弧形段 1 ap""1 aui as:" (23) 在内的动态网格扩展，这使得连铸全过程的实时计 p时da+ 算成为可能 求解质量和动量守恒方程的分离数值方案为在 1.4数值求解 每个时间步长内依次求解方程(21)、(23)和(22). 上述不同的守恒方程的计算均采用有限元法，动 方程(23)中采用压力边界条件（入口流和出口流）， 态Navier-Stokes的计算使用分数步长法，混合对流和 方程(22)中采用速度狄利克雷边界条件.标准有限 扩散方程的求解采用标准Galerkin法和Streamline Up- 元近似是基于加权残差法的伽辽金法，但是这个方 wind Petrov-Galerkin Method (SUPG). 法不适合处理实际关心的多流问题.产生的中心差 1.4.1流体流动 分近似是数值求解困难的原因.在SUPG公式中， 求解NS方程采用了投影法或速度场和压力等 对标准伽辽金加权函数Φ使用流动方向上的流线 阶插值的分离速度压力分步法.为处理不可压缩流 迎风加以修正如下所示： 体，Chorin的于1968年引入了投影法，将不可压缩 性作为一个约束.在Zhu和Sethian的研究工作 =D+△rU中 (24) 后，求解NS方程的中心思想变为先忽略不可压缩 式中，Φ是标准伽辽金加权函数，玉为标准伽辽金 性来更新速度，计算一个时间步，然后投射到不可压 权函数修改后的近似值，△r是迎风参数，U.,是速度 缩流体的空间.这个中心思想是受霍奇分解的启 分量的节点未知解(i定义速度分量，q定义节点 发.霍奇分解描述了定义在一个区域2上的速度场 数)，④，为节点数为q时的标准伽辽金加权函数 山：可被分解为自由发散的部分u,这满足了边界条 在一维中，存在一个最优逆向参数的选择，这个选择 件um:=0和一些标量函数平的梯度.现在模型中 来自于一维对流扩散方程的解析，例如文献8]. 应用的这个方法，压力是隐式的，速度是显式的或半 在Brooks和Hughes9、Zienkiewicz和Taylor D0之

第 12 期 唐娜娜等: 小方坯连铸过程温度场和流场的实时模拟 量、动量和能量守恒方程转换为: ui xi = 0， ( 16) ui t + ( uj － us j ) ui xj = － 1 ρ0 p xi +  xj ( ν + νt ) ui xj + Fi + Si － 2 3 κ xi ，( 17) ρ H t + ρ( uj － us j ) cpT xj =  x (j λ + ρcpνt Pr ) t T xj ． ( 18) κ 和 ε 从式( 9) 和( 10) 计算得到，其用网格计算节 点的速度修改与上述动量方程类似． 在引锭头和板 坯的已经凝固的部分，Eulerian--Lagrangian 方法简化 为单纯的 Lagrangian 方法( uj = us j ) ，因此在式( 18) 中没有 潜 热 的 对 流 释 放． 在浸入式浇口部分， Eulerian--Lagrangian 方法简化为单纯的 Eulerian 方 法( us j = 0) ． 由此，方坯或者板坯的计算域划分为三 个部分，包括钢水入流的浇口在内的方坯或者板坯 的 顶 部 固 定 为 Euler 层，Euler 层 下 边 分 别 为 Expansion层和 Lagrange 层． 在 Expansion 层，一层单 元不断地垂直增长，当它的厚度达到某一临界值时， 单元被分割，下边的分割单元并入 Lagrange 层并随 引锭头以拉坯速度向下移动． 当 Lagrange 层移动进 入弧形段时，网格计算节点的速度将用相应的角速 度表示，角速度等于拉坯速度与弧形段对应半径的 比值． 通过上述 ALE 算法，将能够实现包括弧形段 在内的动态网格扩展，这使得连铸全过程的实时计 算成为可能． 1. 4 数值求解 上述不同的守恒方程的计算均采用有限元法，动 态 Navier-Stokes 的计算使用分数步长法，混合对流和 扩散方程的求解采用标准 Galerkin 法和Streamline Up￾wind Petrov--Galerkin Method ( SUPG) 法． 1. 4. 1 流体流动 求解 N-S 方程采用了投影法或速度场和压力等 阶插值的分离速度压力分步法． 为处理不可压缩流 体，Chorin［15］于 1968 年引入了投影法，将不可压缩 性作为一个约束． 在 Zhu 和 Sethian［16］的研究工作 后，求解 N-S 方程的中心思想变为先忽略不可压缩 性来更新速度，计算一个时间步，然后投射到不可压 缩流体的空间． 这个中心思想是受霍奇分解的启 发． 霍奇分解描述了定义在一个区域 Ω 上的速度场 u* i 可被分解为自由发散的部分 ui，这满足了边界条 件 uimi = 0 和一些标量函数 Ψ 的梯度． 现在模型中 应用的这个方法，压力是隐式的，速度是显式的或半 隐式的( 仅仅发散项) ． 守恒方程( 7) 和( 8) 的半隐 式方法如下． 考虑到在时间步长数 n 的时候速度场 为 un i ，压力场为 pn 、un + 1 i 和 pn + 1满足: un + 1 i xi = 0， ( 19) un + 1 i = un i － Δt·un j un i xj － Δt·1 ρ pn + 1 xi + Δt·  xj ( ν + νt ) un + 1 i xj + Δt·Fn i + Δt·Sn + 1 i － 2 3 kn xi ． ( 20) 定义初步或中间可压缩速度场如下: u* i = un i － Δt·un j un i xj + Δt·  xj ( ν + νt ) u* i xj + Δt·Fn i － Δt·2 3 kn xi ． ( 21) 通过 Hodge 分解投射到自由发散场的公式如下( 达 西法则的源项 Si与多孔媒介的标量梯度类似) : un + 1 i = u* i － Δt·1 ρ pn + 1 xi + Δt·Sn + 1 i ． ( 22) 在 Perot［17］后，这个方法中的速度场 un + 1 i 有一 个误差项等于 Δt 2 ( v + vt ) 2 Δ Δ pn + 1 /ρ． 对于连铸中 的黏度，这个误差项可以忽略，但是如果黏度的数量 级是非常高的( 例如金属的挤压) ，迭代过程是必须 的． 压力刚度方程可由质量守恒方程( 19) 、映射方 程( 22) 得到: 1 ρ  2 pn + 1 x 2 i = 1 Δt u* i xi + Sn + 1 i xi ． ( 23) 求解质量和动量守恒方程的分离数值方案为在 每个时间步长内依次求解方程( 21) 、( 23) 和( 22) ． 方程( 23) 中采用压力边界条件( 入口流和出口流) ， 方程( 22) 中采用速度狄利克雷边界条件． 标准有限 元近似是基于加权残差法的伽辽金法，但是这个方 法不适合处理实际关心的多流问题． 产生的中心差 分近似是数值求解困难的原因． 在 SUPG 公式中， 对标准伽辽金加权函数 Φ 使用流动方向上的流线 迎风加以修正如下所示: Φ 槇 = Φ + ΔτUi，qΦq Φ xi ( 24) 式中，Φ 是标准伽辽金加权函数，Φ 槇 为标准伽辽金 权函数修改后的近似值，Δτ 是迎风参数，Ui，q是速度 分量的节点未知解( i 定义速度分量，q 定义节点 数) ，Φq 为节点数为 q 时的标准伽辽金加权函数． 在一维中，存在一个最优逆向参数的选择，这个选择 来自于一维对流扩散方程的解析，例如文献［18］． 在 Brooks 和 Hughes［19］、Zienkiewicz 和 Taylor［20］ 之 · 7361 ·

·1638+ 北京科技大学学报 第36卷 后，这个方法可以延伸到瞬态和多维案例中： Ar=-B h 2模拟实例与讨论 如前所述，基于ALE算法与Lam-Bremhorst低 Pe= (25) 雷诺数湍流方程，笔者开发了非稳态与时间相关的 温度场和流场计算的实时模型。 式中，Pe是单元佩克莱数，Pe和单元尺寸h定义 2.1几何模型、热物性参数及边界条件设置 如下： 以某钢厂断面为100mm×100mm的连铸小方 =风mm(告舒告） (26) 坯为研究对象，为节省计算时间，本研究只模拟了小 方坯的二维情况，模拟计算域选取小方坯的纵向中 式中，△x△y和△z是单元在x、y和z方向的长度， 心断面，同时考虑到对称性，只选取纵向断面的1/2 山，山，和山，是在x、y和z方向的速度分量.每个时间 开展建模.结晶器长度700mm,铸坯在结晶器中冷 步长的每个单元都要计算逆向参数.SUPG方法用 却形成一定厚度的凝壳以后，离开结晶器进入垂直 于带有混合扩散和对流的所有方程中.在其他方程 段，长度2m,然后进入扇形段，半径3m,扇形段弧 中采用标准伽辽金法.质量矩阵是集中的 度为π2，经过扇形段以后进入水平段，水平段长度 在方程(21)的离散后，求解通过雅克比迭代法 7m,然后切割小方坯.为表征连铸过程铸坯的移 进行迭代，方程(23)的求解通过不完整的柯列斯基 动，将其计算域划分为Euler、Expansion和Lagrange 共轭梯度法进行迭代，这源于Kershaw的研究u. 三层（图1)，随着引锭头的向下移动，Expansion层 方程(22)的求解采用替代法，因为在质量矩阵集中 不断地向下扩展，并且以5mm高度为一个新网格连 之后，这一步中不涉及矩阵.在上面的流体流动方 续不断地加入到Lagrange层中，而结晶器和Euler 程的瞬态求解中，单元库朗数对于数值方案的精确 层的位置保持不动 和稳定是非常重要的.模型通过增大或减小与来自 最后时间步长的有关的时间步长来控制计算中最高 ■ 单元库朗数 1.4.2能量方程 在能量方程中，内部的能源项和扩散项必须是 隐式的：否则，就需要特别小的时间步长（或每个时 Eler层 间步长的迭代数)，计算时间较长.在动态起步阶 段，能量方程中的时间偏导数项的精确是非常重要 2 Expansion层 的.焓是一个独立的变量，温度作为焓的函数通过 下面数学关系的使用被线性化四： Lagrange层 aT ah dT(x）_ax;ar (27) dh (x)ahah ox;ax 这个表达式在每个单元的中心点计算.在时间离散 引锭头 后，能量式(13)转换为： 8ha+1 +p则+p,(g1-4)C- s oh" ox; 图1计算模型 aT"ah" Fig.1 Computational model 日a+凸)，ash (28) ax;( Pr.ah"ah"ax; 连铸小方坯材质为X70钢，部分热物性参数如 dxi dxi 比热容、热导率、运动黏度等都是与温度有关的函 这里n为时间步长（温度和焓已知），n+1代表了新 数，具体数值详见表1.密度是与温度无关的常数， 的时间步长（温度和焓未知）.空间通过SUPG方法 为7800kg°m-3.固相线温度、液相线温度和浇铸温 进行离散，产生的代数方程通过ICCG(O)方法进行 度分别是1461、1520和1540℃，熔化潜热为256.4 求解. kJ.kg-1

北 京 科 技 大 学 学 报 第 36 卷 后，这个方法可以延伸到瞬态和多维案例中: Δτ = β 槡15 h V ，β = coth ( Pe ) 2 － 2 Pe， Pe = u2 槡i ·h ν ． ( 25) 式中，Pe 是单元佩克莱数，Pe 和单元尺寸 h 定义 如下: h = u2 槡i ·min ( Δx | ux | ，Δy | uy | ，Δz | uz ) | ． ( 26) 式中，Δx、Δy 和 Δz 是单元在 x、y 和 z 方向的长度， ux、uy和 uz是在 x、y 和 z 方向的速度分量． 每个时间 步长的每个单元都要计算逆向参数． SUPG 方法用 于带有混合扩散和对流的所有方程中． 在其他方程 中采用标准伽辽金法． 质量矩阵是集中的． 在方程( 21) 的离散后，求解通过雅克比迭代法 进行迭代，方程( 23) 的求解通过不完整的柯列斯基 共轭梯度法进行迭代，这源于 Kershaw 的研究［21］． 方程( 22) 的求解采用替代法，因为在质量矩阵集中 之后，这一步中不涉及矩阵． 在上面的流体流动方 程的瞬态求解中，单元库朗数对于数值方案的精确 和稳定是非常重要的． 模型通过增大或减小与来自 最后时间步长的有关的时间步长来控制计算中最高 单元库朗数． 1. 4. 2 能量方程 在能量方程中，内部的能源项和扩散项必须是 隐式的; 否则，就需要特别小的时间步长( 或每个时 间步长的迭代数) ，计算时间较长． 在动态起步阶 段，能量方程中的时间偏导数项的精确是非常重要 的． 焓是一个独立的变量，温度作为焓的函数通过 下面数学关系的使用被线性化［22］: dT( xj ) dh( xj ) = T xj h xj h xj h xj ． ( 27) 这个表达式在每个单元的中心点计算． 在时间离散 后，能量式( 13) 转换为: ρ hn + 1 t + ρus j hn xj + ρcp ( un + 1 j － us j ) Tn xj =  x (j λ + ρcpνt Pr ) t Tn xj hn xj hn xj hn xj hn + 1 xj ． ( 28) 这里 n 为时间步长( 温度和焓已知) ，n + 1 代表了新 的时间步长( 温度和焓未知) ． 空间通过 SUPG 方法 进行离散，产生的代数方程通过 ICCG( 0) 方法进行 求解． 2 模拟实例与讨论 如前所述，基于 ALE 算法与 Lam-Bremhorst 低 雷诺数湍流方程，笔者开发了非稳态与时间相关的 温度场和流场计算的实时模型． 2. 1 几何模型、热物性参数及边界条件设置 以某钢厂断面为 100 mm × 100 mm 的连铸小方 坯为研究对象，为节省计算时间，本研究只模拟了小 方坯的二维情况，模拟计算域选取小方坯的纵向中 心断面，同时考虑到对称性，只选取纵向断面的 1 /2 开展建模． 结晶器长度 700 mm，铸坯在结晶器中冷 却形成一定厚度的凝壳以后，离开结晶器进入垂直 段，长度 2 m，然后进入扇形段，半径 3 m，扇形段弧 度为 π/2，经过扇形段以后进入水平段，水平段长度 7 m，然后切割小方坯． 为表征连铸过程铸坯的移 动，将其计算域划分为 Euler、Expansion 和 Lagrange 三层( 图 1) ，随着引锭头的向下移动，Expansion 层 不断地向下扩展，并且以 5 mm 高度为一个新网格连 续不断地加入到 Lagrange 层中，而结晶器和 Euler 层的位置保持不动． 图 1 计算模型 Fig． 1 Computational model 连铸小方坯材质为 X70 钢，部分热物性参数如 比热容、热导率、运动黏度等都是与温度有关的函 数，具体数值详见表 1． 密度是与温度无关的常数， 为 7800 kg·m － 3 ． 固相线温度、液相线温度和浇铸温 度分别是 1461、1520 和 1540 ℃，熔化潜热为 256. 4 kJ·kg － 1 ． · 8361 ·

第12期 唐娜娜等：小方还连铸过程温度场和流场的实时模拟 ·1639· 表1X70钢的部分热物理性能 Table 1 Some thermophysical properties of X70 steel 温度/ 比热容/ 热导率/ 运动黏度/ 温度/ 比热容/ 热导率/ 运动黏度/ ℃ (Jkg-1.K-1) (W.m-1.K-1) (m2.s1) ℃ (小kg1k) (Wm1.K-1) (m2s1) 0 388.0 40.0 700 904.5 32.3 7.71×1011 50 410.7 41.5 750 904.5 30.6 7.72×101 100 434.8 43.0 800 904.5 25.4 7.70×101 200 487.3 42.0 860 904.5 26.0 7.72×1011 300 546.2 41.0 1450 1081.5 31.4 8.06×1011 400 612.1 40.0 1470 1067.5 40.0 1.60×10-5 500 686.1 39.0 1520 1067.5 200 8.03×10-7 600 768.9 35.5 2000 1067.5 200 3.27×10- 模拟过程中没有考虑结晶器与铸坯之间的空气 后将被切割，该切割面设为绝热边界条件.连铸全 层，铸坯与结晶器之间的换热、铸坯与冷却水之间的 过程铸造速度取为l.2m·min,壁面设为无滑移边 换热、铸坯与空气之间的换热选择来自某钢厂的经 界条件，随着钢液的凝固，速度边界条件从壁面的零 验值，该换热系数是距离弯液面距离的函数，部分换 速度逐渐转变为铸造速度，固液界面考虑了由于速 热系数经验值列于表2.二冷水温度取为定值 度不同而引起的界面摩擦，关于界面摩擦参见前面 100℃，对称轴处的热流设为零.铸坯在经过水平段 的数学模型部分的阐述 表2模拟所用部分换热系数 Table 2 Partial heat transfer coefficient used in this study 距弯液面 换热系数/ 距弯液面 换热系数/ 距弯液面 换热系数/ 距弯液面 换热系数/ 距离/m (Wm2·K1) 距离/m (Wm2·K-l) 距离/m (Wm2·K1) 距离/m (Wm2·Kl) 0.100 5800 1.171 7800 1.309 7800 1.466 6550 0.300 4800 1.172 700 1.310 700 1.467 700 0.400 3900 1.202 700 1.354 700 1.509 700 0.500 3100 1.203 60 1.355 60 1.510 60 1.100 2400 1.224 60 1.379 60 1.534 60 1.101 700 1.225 700 1.380 700 1.147 700 1.284 700 1.438 700 1.148 7800 1.285 7800 1.439 6550 2.2模拟结果与讨论 放使得温度下降缓慢，在降低到固相线以下的温 2.2.1温度场模拟 度后，温度下降较快 图2为连铸全过程中铸坯表面及心部两个监 铸坯坯壳的形成过程可由铸坯的固相分数分 控点1、2处温度随时间的变化曲线，两个监控点 布图表示.图3为铸坯垂直段不同时刻的固相分 的具体位置如图1计算模型所示.从图中可以看 数，最下端蓝色区域为引锭头.由图可以看出，随 出，结晶器内（曲线ab段）的冷却强度最大，表面 着拉锭的进行，铸坯表面逐渐形成凝壳，而且随着 温度在结晶器内下降最快，在结晶器出口（点b) 时间的增加，凝固壳厚度逐渐增大，但是铸坯在 处，由于冷却强度的降低出现了明显的温度回升.120s和150s时的固相分数分布几乎没有区别，这 之后，由于钢水凝固结晶潜热的释放量小于喷水表明垂直段铸坯的凝固行为在120s以后己经达 冷却的强度，温度又出现了下降.在二冷区（曲线 到了稳态.图4为120s和150s时铸坯垂直段不 bc段)结束后，冷却方式变为空冷，冷却强度低，所 同位置处的固相分数分布.坯壳的厚度也就是固 以温度又出现了回升.心部温度开始时由于处在 相分数等于1所对应的x方向的长度.由图可以 固相线以上的温度范围，钢水凝固结晶潜热的释 看出，坯壳厚度从y=0m到y=-2m是逐渐增大

第 12 期 唐娜娜等: 小方坯连铸过程温度场和流场的实时模拟 表 1 X70 钢的部分热物理性能 Table 1 Some thermophysical properties of X70 steel 温度/ ℃ 比热容/ ( J·kg － 1·K － 1 ) 热导率/ ( W·m － 1·K － 1 ) 运动黏度/ ( m2 ·s － 1 ) 0 388. 0 40. 0 — 50 410. 7 41. 5 — 100 434. 8 43. 0 — 200 487. 3 42. 0 — 300 546. 2 41. 0 — 400 612. 1 40. 0 — 500 686. 1 39. 0 — 600 768. 9 35. 5 — 温度/ ℃ 比热容/ ( J·kg － 1·K － 1 ) 热导率/ ( W·m － 1·K － 1 ) 运动黏度/ ( m2 ·s － 1 ) 700 904. 5 32. 3 7. 71 × 1011 750 904. 5 30. 6 7. 72 × 1011 800 904. 5 25. 4 7. 70 × 1011 860 904. 5 26. 0 7. 72 × 1011 1450 1081. 5 31. 4 8. 06 × 1011 1470 1067. 5 40. 0 1. 60 × 10 － 5 1520 1067. 5 200 8. 03 × 10 － 7 2000 1067. 5 200 3. 27 × 10 － 7 模拟过程中没有考虑结晶器与铸坯之间的空气 层，铸坯与结晶器之间的换热、铸坯与冷却水之间的 换热、铸坯与空气之间的换热选择来自某钢厂的经 验值，该换热系数是距离弯液面距离的函数，部分换 热系数 经 验 值 列 于 表 2． 二冷水温度取为定值 100 ℃，对称轴处的热流设为零． 铸坯在经过水平段 后将被切割，该切割面设为绝热边界条件． 连铸全 过程铸造速度取为 1. 2 m·min － 1，壁面设为无滑移边 界条件，随着钢液的凝固，速度边界条件从壁面的零 速度逐渐转变为铸造速度，固液界面考虑了由于速 度不同而引起的界面摩擦，关于界面摩擦参见前面 的数学模型部分的阐述． 表 2 模拟所用部分换热系数 Table 2 Partial heat transfer coefficient used in this study 距弯液面 距离/m 换热系数/ ( W·m － 2·K － 1 ) 0. 100 5800 0. 300 4800 0. 400 3900 0. 500 3100 1. 100 2400 1. 101 700 1. 147 700 1. 148 7800 距弯液面 距离/m 换热系数/ ( W·m － 2·K － 1 ) 1. 171 7800 1. 172 700 1. 202 700 1. 203 60 1. 224 60 1. 225 700 1. 284 700 1. 285 7800 距弯液面 距离/m 换热系数/ ( W·m － 2·K － 1 ) 1. 309 7800 1. 310 700 1. 354 700 1. 355 60 1. 379 60 1. 380 700 1. 438 700 1. 439 6550 距弯液面 距离/m 换热系数/ ( W·m － 2·K － 1 ) 1. 466 6550 1. 467 700 1. 509 700 1. 510 60 1. 534 60 2. 2 模拟结果与讨论 2. 2. 1 温度场模拟 图 2 为连铸全过程中铸坯表面及心部两个监 控点 1、2 处温度随时间的变化曲线，两个监控点 的具体位置如图 1 计算模型所示． 从图中可以看 出，结晶器内( 曲线 ab 段) 的冷却强度最大，表面 温度在结晶器内下降最快，在结晶器出口( 点 b) 处，由于冷却强度的降低出现了明显的温度回升． 之后，由于钢水凝固结晶潜热的释放量小于喷水 冷却的强度，温度又出现了下降． 在二冷区( 曲线 bc 段) 结束后，冷却方式变为空冷，冷却强度低，所 以温度又出现了回升． 心部温度开始时由于处在 固相线以上的温度范围，钢水凝固结晶潜热的释 放使得温度下降缓慢，在降低到固相线以下的温 度后，温度下降较快． 铸坯坯壳的形成过程可由铸坯的固相分数分 布图表示． 图 3 为铸坯垂直段不同时刻的固相分 数，最下端蓝色区域为引锭头． 由图可以看出，随 着拉锭的进行，铸坯表面逐渐形成凝壳，而且随着 时间的 增 加，凝 固 壳 厚 度 逐 渐 增 大，但 是 铸 坯 在 120 s 和 150 s 时的固相分数分布几乎没有区别，这 表明垂直段铸坯的凝固行为在 120 s 以后已经达 到了稳态． 图 4 为 120 s 和 150 s 时铸坯垂直段不 同位置处的固相分数分布． 坯壳的厚度也就是固 相分数等于 1 所对应的 x 方向的长度． 由图可以 看出，坯壳厚度从y = 0 m到 y = － 2 m 是逐渐增大 · 9361 ·

·1640 北京科技大学学报 第36卷 1600 1.0 ■—1208 y=0m -。—150s 1500 0.5 心部 1400 0 0.01 002003 0.040.05 1300 1.0 1200 表面 0.5 )y-0.5m 1100 0.01 0.02 0.030.04 0.05 1000 0 1.0- 900 =-1m 0.5 800 1000 100200300 400500600 700 0 时间/s 0 0.01 0.020.03 0.040.05 图2铸坯表面及心部的温度变化曲线 1.0 y=-15m Fig.2 Temperature change of the slab surface and center 0.5 的，并且在这些位置处120s和150s时的固相分 0 0 0.01 0.02 0.030.040.05 数分布曲线是重合的，这也更准确地证实了上文 提到的垂直段铸坯的凝固行为在120s以后达到 =-2m 0.5 稳态的结论 0 0.01 0.02 0.03 0.040.05 xim 图4120：和150s时铸坯垂直段不同位置处的固相分数 Fig.4 Solid fraction for different locations of the billet vertical part at 固相分数 120s and 150s 1.0 0.9 全凝固的坯壳部分.根据此图，可以算出铸坯在出 .8 结晶器时的坯壳厚度大约为12.5mm,这与工厂实 07 测值13mm非常接近，说明本研究所开发的有限 0.6 元分析模型用于钢铁连铸数值模拟的可行性 0.5 2.2.2流场模拟 图6为计算的湍动能和速度场在不同时刻的分 布轮廓图.由图可以看出，在30s和60s时的湍动 能和速度场分布几乎看不出什么区别，说明流场在 0 30s以后就基本趋于稳定. 5s10s20850s80=120=150s 图7为30s和60s时不同位置处的流场y方 图3铸坯垂直段不同时刻的固相分数 向的速度山，拉坯方向为y的负方向.由图可知，y Fig.3 Slab solid fraction at different time for the vertical part of a billet 方向的速度沿x方向存在波动，不同位置处的流场 y方向的速度在x方向的变化趋势是相同的，都是 图5(a)为连铸全过程不同时刻铸坯的温度分 在处于液体时(x=0m附近)速度最大.在y= 布，铸坯最低端的灰色区域为引锭头.由图可以看 0.095m和y=0.035m处大约为0.11m·s-1,在 出，400s时铸坯已经历了垂直段与扇形段，并且拉 y=-0.035m处大约为0.07ms1.随着x增大， 到了水平段，观察发现在扇形段的结束处铸坯心 y方向的速度逐渐变小至零，然后速度反向，速度 部的温度己经降到1461℃以下，即固相线温度以 先变大后减小，减小至零后又开始变大，在固体坯 下，也就是说，铸坯在达到水平段前（接近水平段 壳处速度达到0.02m·s,也就是以拉坯速度运 时)心部己经完全凝固.图5(b)为结晶器内的温 动.在30s和60s时，y方向的速度在x方向的变 度分布，图中灰色区域为结晶器内温度在固相线 化曲线几乎重合在一起，进一步说明流场在30s 温度1461℃以下的铸坯表面，也就是铸坯己经完 以后趋于稳定

北 京 科 技 大 学 学 报 第 36 卷 图 2 铸坯表面及心部的温度变化曲线 Fig． 2 Temperature change of the slab surface and center 的，并且在这些位置处 120 s 和 150 s 时的固相分 数分布曲线是重合的，这也更准确地证实了上文 提到的垂直段铸坯的凝固行为在 120 s 以后达到 稳态的结论． 图 3 铸坯垂直段不同时刻的固相分数 Fig． 3 Slab solid fraction at different time for the vertical part of a billet 图 5( a) 为连铸全过程不同时刻铸坯的温度分 布，铸坯最低端的灰色区域为引锭头． 由图可以看 出，400 s 时铸坯已经历了垂直段与扇形段，并且拉 到了水平段，观察发现在扇形段的结束处铸坯心 部的温度已经降到 1461 ℃ 以下，即固相线温度以 下，也就是说，铸坯在达到水平段前( 接近水平段 时) 心部已经完全凝固． 图 5( b) 为结晶器内的温 度分布，图中灰色区域为结晶器内温度在固相线 温度 1461 ℃ 以下的铸坯表面，也就是铸坯已经完 图 4 120 s 和 150 s 时铸坯垂直段不同位置处的固相分数 Fig． 4 Solid fraction for different locations of the billet vertical part at 120 s and 150 s 全凝固的坯壳部分． 根据此图，可以算出铸坯在出 结晶器时的坯壳厚度大约为 12. 5 mm，这与工厂实 测值 13 mm 非常接近，说明本研究所开发的有限 元分析模型用于钢铁连铸数值模拟的可行性． 2. 2. 2 流场模拟 图 6 为计算的湍动能和速度场在不同时刻的分 布轮廓图． 由图可以看出，在 30 s 和 60 s 时的湍动 能和速度场分布几乎看不出什么区别，说明流场在 30 s 以后就基本趋于稳定． 图 7 为 30 s 和 60 s 时不同位置处的流场 y 方 向的速度 uy，拉坯方向为 y 的负方向． 由图可知，y 方向的速度沿 x 方向存在波动，不同位置处的流场 y 方向的速度在 x 方向的变化趋势是相同的，都是 在处于液 体 时( x = 0 m 附 近) 速 度 最 大． 在 y = 0. 095 m 和 y = 0. 035 m 处大约为 0. 11 m·s － 1，在 y = － 0. 035 m 处大约为 0. 07 m·s － 1 ． 随着 x 增大， y 方向的速度逐渐变小至零，然后速度反向，速度 先变大后减小，减小至零后又开始变大，在固体坯 壳处速度达到 0. 02 m·s － 1，也就是以拉坯速度运 动． 在 30 s 和 60 s 时，y 方向的速度在 x 方向的变 化曲线几乎重合在一起，进一步说明流场在 30 s 以后趋于稳定． · 0461 ·

第12期 唐娜娜等：小方还连铸过程温度场和流场的实时模拟 ·1641· 温度℃ 温度℃ 1540.0 1540 1532.1 1476 1524.2 1412 1516.3 30 -1348 1508.4 1284 1500.5 1220 1492.6 110s 1156 1484.7 180s -1092 300s 1476.8 -1028 1468.9 964 1461.0 -900 400s 图5铸坯的凝固.(：)连铸全过程的温度分布：(b)结品器内的温度分布 Fig.5 Billet solidification:(a)temperature distribution during the entire continuous casting process:(b)temperature distribution in the mould region 速度m*。 0.120 湍动能刀 0.108 0.00070(a 0.096 0.00063 -0.084 000056 0.072 0.00049 0.060 0.0042 0.00035 -0.048 0.00028 0.036 0.00021 0.024 0.0014 0.012 0.00007 15830s 60s 15s 30g 60 图6不同时刻计算的湍动能(a）和速度场(b) Fig.6 Turbulent energy (a)and speed field (b)calculated at the different time 3结论 和流场的变化.出结晶器时铸坯坯壳厚度12.5mm, 这与实测值I3mm比较接近，模拟效果良好.根据 (1)针对钢铁连铸非稳态的铸造过程，建立了 不同时刻的固相分数分布得出垂直段的铸坯凝固行 模拟传热和流体流动的实时模拟模型.模型开发过 为在120s以后达到稳态.模拟的流场趋于稳定所 程充分考虑了湍流的影响，其计算采用了Lam- 需时间较短，时间大约30s.y方向的速度沿x方向 Bremhorst低雷诺数k一e方程，并对糊状区使用与液 存在波动，液体处的速度最大，固体坯壳处以拉坯速 相分数相关的达西源项进行处理.为了模拟连续拉 度0.02ms1向下运动. 坯的过程，采用任意拉格朗日一欧拉(ALE)方法实 (3)所开发的模型因为其可以实现包括连铸弧 现了移动计算域下的动态网格扩展. 形段在内的连铸全过程温度与流场的模拟，较 (2)利用建立的实时模型模拟了X70钢小方 ProCAST软件采用MLE算法等无法实现连铸弧形 坯连铸全过程（垂直段+弧形段+水平段）温度场 段的非稳态计算而言，展示出明显的技术优势.连

第 12 期 唐娜娜等: 小方坯连铸过程温度场和流场的实时模拟 图 5 铸坯的凝固． ( a) 连铸全过程的温度分布; ( b) 结晶器内的温度分布 Fig． 5 Billet solidification: ( a) temperature distribution during the entire continuous casting process; ( b) temperature distribution in the mould region 图 6 不同时刻计算的湍动能( a) 和速度场( b) Fig． 6 Turbulent energy ( a) and speed field ( b) calculated at the different time 3 结论 ( 1) 针对钢铁连铸非稳态的铸造过程，建立了 模拟传热和流体流动的实时模拟模型． 模型开发过 程充分考虑了湍流的影响，其 计 算 采 用 了 Lam￾Bremhorst 低雷诺数 k--ε 方程，并对糊状区使用与液 相分数相关的达西源项进行处理． 为了模拟连续拉 坯的过程，采用任意拉格朗日--欧拉( ALE) 方法实 现了移动计算域下的动态网格扩展． ( 2) 利用建立的实时模型模拟了 X70 钢小方 坯连铸全过程( 垂直段 + 弧形段 + 水平段) 温度场 和流场的变化． 出结晶器时铸坯坯壳厚度 12. 5 mm， 这与实测值 13 mm 比较接近，模拟效果良好． 根据 不同时刻的固相分数分布得出垂直段的铸坯凝固行 为在 120 s 以后达到稳态． 模拟的流场趋于稳定所 需时间较短，时间大约 30 s． y 方向的速度沿 x 方向 存在波动，液体处的速度最大，固体坯壳处以拉坯速 度 0. 02 m·s － 1向下运动． ( 3) 所开发的模型因为其可以实现包括连铸弧 形段在内的连铸全过程温度与流场的模拟，较 ProCAST 软件采用 MILE 算法等无法实现连铸弧形 段的非稳态计算而言，展示出明显的技术优势． 连 · 1461 ·

·1642 北京科技大学学报 第36卷 0.04 flow and surface wave in continuous casting moulds /Proceedings 0y=0.095m of China Iron Steel Annual Meeting.Beijing,2009:1561 -0.04 (张焕鑫，刘瑞，沈厚发，等.薄板坯连铸结品器流场与液位 -0.08 -0.12 ●ees 波动不稳定性的模拟研究/中国钢铁年会论文集.北京， 0.010.020.030.040.05 2009:1561) 0.04 8] Xue J G,Jin X W,Yao G T.Numerical simulation of solidifica- 0 3=0.035m tion process in slab casting mould based on3D mile algorithm mod- -0.04 -5-55-38656880 -0.08 el.J Aeronaut Mater,2010,30(3):10 -0.12 (薛建国，金学伟，姚耕耘.三维MLE算法模型模拟连铸结 0 0.01 0.02 0.030.04 0.05 品器内的板坯凝固过程.航空材料学报，2010,30(3)：10) 0.02 y=-0.035m 9] Li R.Casting Process Simulation ProCAST from Entry to the Mas- -0.02 ter.Beijing:China Water Power Press,2010 -0.04 -0.06 (李日.铸造工艺仿真ProCAST从入门到精通.北京：中国水 -0.08 0 0.010.020.030.040.05 利水电出版社，2010) x/m [10]Bennon W D,Incropera F P.A continuum model for momentum, 图730s和60s时不同位置处的流场y方向的速度 heat and species transport in binary solid-iquid phase change sys- Fig.7 Fluid velocity in the y direction at different locations at 30s tems:I.Model formulation.Int J Heat Mass Transfer,1987, and 60s 30(10):2161 [11]Launder B E,Sharma B I.Application of the energy-dissipation 铸全过程温度与流动模拟的实现为模拟该过程的应 model of turbulence to the calculation of flow near a spinning 力和应变场提供了前提条件 disc.Lett Heat Mass Transfer,1974,1(2):131 [12]Fjaer H G,Mortensen D,Hakonsen A,et al.Coupled stress, thermal and fluid flow modelling of the start-up phase of alumin- 参考文献 ium sheet ingot casting /Light Metals 1999.Warrendale:The [1]Liu H P,Yang C Z,Zhang H,et al.Numerical simulation of flu- Minerals,Metals Materials Society (TMS),1999:743 id flow and thermal characteristics of thin slab in the funnel-type [13]Shyy W,Pang Y,Hunter G B,et al.Effect of turbulent heat molds of two casters.ISI/Int,2011,51(3):392 transfer on continuous ingot solidification.J Eng Mater Technol, Zhu M Y,Wang J,Zhang Y.Numerical simulation of fluid flow 1993,115:8 and heat transfer in funnel shaped mold of thin slab continuous 04] Ghosh S.Finite element simulation of some extrusion processes caster.J Iron Steel Res,2005,12(6)14 using arbitrary Lagrangian-Eulerian description.Mater Shaping B]Zhang L,Kang J W,Shen H F,et al.Finite point method for the 7 echnol,1990,8(1):53 simulation of solidification and heat transfer in continuous casting 5]Chorin A J.Numerical solution of the Navier-Stokes equations. mold.Tsinghua Sci Technol,2004,9(5):570 Math Comput,1968,22(104):745 4]Chang Y H,Zhang JQ,Qian H Z,et al.Development and appli- [16]Zhu J,Sethian J.Projection methods coupled to level set inter- cation of a dynamic secondary cooling control model for beam face techniques.J Comput Phys,1992,102 (1):128 blanks based on finite element method.I Unir Sci Technol Bei- [17]Perot J B.An analysis of the fractional step method.Comput img,2011,33(4):418 Phs,1993,108(1):51 (常运合，张家泉，钱宏智，等.基于有限元法异型坯动态二 18] Donea J.Generalized Galerkin methods for convection dominated 冷控制模型开发与应用.北京科技大学学报，2011,33(4)： transport phenomena.Appl Mech Rer,1991,44:205 418) [19]Brooks N,Hughes J R.Streamline upwind/Petrov-Galerkin for- [5]Wang Y C,Zhu L G.Li D Y,et al.The numerical simulation mulations for convection dominated flows with particular emphasis and control of the solidification and heat transfer in the entire on the incompressible Navier-Stokes equations.Comput Methods process of continuous casting.Spec Cast Nonferrous Alloys,2005 Appl Mech Eng,1982,32(1):199 25(9):531 20 Zienkiewicz O C,Taylor R L.The finite element method (王迎春，朱立光，李大永，等.连铸凝固传热全过程数值模 McGran-Hill Comp,1991,2 :320 拟与控制.特种铸造及有色合金，2005,25(9)：531) [21]Kershaw D S.The incomplete Cholesky-Conjugate gradient ]Ludwig O,Aloe M,Thevoz P.State of the art in modelling of con- method for the iterative solution of systems of linear equations.J tinuous casting /Proceedings of the 6th European Conference on Comp P%s,1978,26(1):43 Continuous Casting.Switzerland,2008 22]Amiez G,Gremaud P A.On a numerical approach to Stefan-ike Zhang H X,Liu R,Shen H F,et al.Modeling of unsteady fluid problems.Numer Math,1991,59(1):71

北 京 科 技 大 学 学 报 第 36 卷 图 7 30 s 和 60 s 时不同位置处的流场 y 方向的速度 Fig． 7 Fluid velocity in the y direction at different locations at 30 s and 60 s 铸全过程温度与流动模拟的实现为模拟该过程的应 力和应变场提供了前提条件． 参 考 文 献 ［1］ Liu H P，Yang C Z，Zhang H，et al． Numerical simulation of flu￾id flow and thermal characteristics of thin slab in the funnel-type molds of two casters． ISIJ Int，2011，51( 3) : 392 ［2］ Zhu M Y，Wang J，Zhang Y． Numerical simulation of fluid flow and heat transfer in funnel shaped mold of thin slab continuous caster． J Iron Steel Ｒes，2005，12( 6) : 14 ［3］ Zhang L，Kang J W，Shen H F，et al． Finite point method for the simulation of solidification and heat transfer in continuous casting mold． Tsinghua Sci Technol，2004，9( 5) : 570 ［4］ Chang Y H，Zhang J Q，Qian H Z，et al． Development and appli￾cation of a dynamic secondary cooling control model for beam blanks based on finite element method． J Univ Sci Technol Bei￾jing，2011，33( 4) : 418 ( 常运合，张家泉，钱宏智，等． 基于有限元法异型坯动态二 冷控制模型开发与应用． 北京科技大学学报，2011，33( 4) : 418) ［5］ Wang Y C，Zhu L G，Li D Y，et al． The numerical simulation and control of the solidification and heat transfer in the entire process of continuous casting． Spec Cast Nonferrous Alloys，2005， 25( 9) : 531 ( 王迎春，朱立光，李大永，等． 连铸凝固传热全过程数值模 拟与控制． 特种铸造及有色合金，2005，25( 9) : 531) ［6］ Ludwig O，Aloe M，Thevoz P． State of the art in modelling of con￾tinuous casting / / Proceedings of the 6th European Conference on Continuous Casting． Switzerland，2008 ［7］ Zhang H X，Liu Ｒ，Shen H F，et al． Modeling of unsteady fluid flow and surface wave in continuous casting moulds / / Proceedings of China Iron ＆ Steel Annual Meeting． Beijing，2009: 1561 ( 张焕鑫，刘瑞，沈厚发，等． 薄板坯连铸结晶器流场与液位 波动不稳定性的模拟研究 / / 中国钢铁年会论文集． 北京， 2009: 1561) ［8］ Xue J G，Jin X W，Yao G T． Numerical simulation of solidifica￾tion process in slab casting mould based on 3D mile algorithm mod￾el． J Aeronaut Mater，2010，30( 3) : 10 ( 薛建国，金学伟，姚耕耘． 三维 MiLE 算法模型模拟连铸结 晶器内的板坯凝固过程． 航空材料学报，2010，30( 3) : 10) ［9］ Li Ｒ． Casting Process Simulation ProCAST from Entry to the Mas￾ter． Beijing: China Water Power Press，2010 ( 李日． 铸造工艺仿真 ProCAST 从入门到精通． 北京: 中国水 利水电出版社，2010) ［10］ Bennon W D，Incropera F P． A continuum model for momentum， heat and species transport in binary solid-liquid phase change sys￾tems: Ⅰ． Model formulation． Int J Heat Mass Transfer，1987， 30( 10) : 2161 ［11］ Launder B E，Sharma B I ． Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc． Lett Heat Mass Transfer，1974，1( 2) : 131 ［12］ Fjr H G，Mortensen D，Hakonsen A，et al． Coupled stress， thermal and fluid flow modelling of the start-up phase of alumin￾ium sheet ingot casting / / Light Metals 1999． Warrendale: The Minerals，Metals ＆ Materials Society ( TMS) ，1999: 743 ［13］ Shyy W，Pang Y，Hunter G B，et al． Effect of turbulent heat transfer on continuous ingot solidification． J Eng Mater Technol， 1993，115: 8 ［14］ Ghosh S． Finite element simulation of some extrusion processes using arbitrary Lagrangian-Eulerian description． J Mater Shaping Technol，1990，8( 1) : 53 ［15］ Chorin A J． Numerical solution of the Navier--Stokes equations． Math Comput，1968，22( 104) : 745 ［16］ Zhu J，Sethian J． Projection methods coupled to level set inter￾face techniques． J Comput Phys，1992，102( 1) : 128 ［17］ Perot J B． An analysis of the fractional step method． J Comput Phys，1993，108( 1) : 51 ［18］ Donea J ． Generalized Galerkin methods for convection dominated transport phenomena． Appl Mech Ｒev，1991，44: 205 ［19］ Brooks N，Hughes J Ｒ． Streamline upwind /Petrov--Galerkin for￾mulations for convection dominated flows with particular emphasis on the incompressible Navier--Stokes equations． Comput Methods Appl Mech Eng，1982，32( 1) : 199 ［20］ Zienkiewicz O C，Taylor Ｒ L． The finite element method， McGraw-Hill Comp，1991，2 : 320 ［21］ Kershaw D S． The incomplete Cholesky-- Conjugate gradient method for the iterative solution of systems of linear equations． J Comp Phys，1978，26( 1) : 43 ［22］ Amiez G，Gremaud P A． On a numerical approach to Stefan-like problems． Numer Math，1991，59( 1) : 71 · 2461 ·