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

第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 ·