DOI: 10.11817/j.issn.1672-7207.2015.04.045
罩式退火炉热过程数学模型与实验验证
杨培培1，温治1,2，豆瑞锋1
(1. 北京科技大学 机械工程学院，北京，100083；
2. 北京科技大学 冶金工业节能减排北京市重点实验室，北京，100083)
摘要：考虑对流板复杂的几何结构，建立三维模型，采用正交试验法模拟获得流量、外径、和炉内温度对炉内各处努塞尔数Nu分布的影响，得到流量在考察的3种因素内流量的影响作用最强的结论。接着进一步分析Nu数随流量的变化，并拟合得到各处Nu与流量的关系式；为了解决炉内复杂几何形状条件下的辐射换热角系数求解问题，采用Monte Carlo方法建立炉内辐射换热角系数求解模型，并最终采用热辐射网络图法求解炉内辐射换热。通过现场进行热电偶插片实验验证本文所建立模型的准确性。研究结果表明：模拟计算得到的温度值与现场实测值基本一致，证明了模型的准确可靠性。
关键词：正交试验法；蒙特卡洛法；热电偶插片实验；罩式退火炉；数学模型
中图分类号：TF062 文献标志码：A 文章编号：1672-7207(2015)04-1518-09
Thermal process mathematical model and experimental verification in Bell-type annealing furnace
YANG Peipei1, WEN Zhi1,2, DOU Ruifeng1
(1. School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China;
2. Beijing Key Laboratory for Energy Saving and Emission Reduction of Metallurgical Industry,
University of Science and Technology Beijing, Beijing 100083, China)
Abstract: Considering the complex structure characteristics of convection plate, three dimensional model and effect of flow rate, steel coil outer diameter and furnace temperature on Nu distribution were built using orthogonal test. The results show that flow rate has strongest influence on the range of consideration. Moreover, Nu with the flow rate change was analyzed, and the relationship between flow rate and Nu was obtained. In order to solve the problem of solving the angle factor under the condition of complex geometry, using Monte Carlo method the model of radiative heat transfer angle factor in furnace was set, and the heat radiation network diagram for radiation heat transfer was used in the end. In addition, thermocouple insert experiment was also used to verify the accuracy of the model. The results show that the simulation values are consistent with experimental data.
Key words: orthogonal test method; Monte Carlo method; thermocouple insert experiment; batch furnace; mathematical model
带钢是以钢卷的形式在全氢罩式退火炉中退火的，影响传热的主要参数是径向等效导热系数、对流换热系数和辐射热流密度。目前国内外许多学者[1-2]针对接触传热问题展开了深入广泛的研究，得到了较为准确的径向等效导热系数的理论计算公式。相比而言，对内罩空间内各处对流换热系数和辐射换热的研究，则进行了过多的简化假设，使结果的计算精度降低。保护性循环气体在炉内的流量和流速分配直接影响表面对流换热，因此，准确计算炉内保护气体流量分配是准确计算炉内对流换热的关键。李卫杰等[3-4]采用流动阻力分析和预测校正方法，将导流板结构简化为一定厚度的中间层和部分界面流通的上下两层锲形槽层，对炉内循环气体流量分配进行了模型化并对影响流动分配的相关影响因素进行了敏感性测试。林林等[5]通过对氢气流动阻力分析，采用预测校正法，在炉型尺寸、对流板结构尺寸、风机循环风量一定的前提下，将对流板结构进行简化，进而确定了不同钢卷堆垛数、不同钢卷尺寸时炉内各处的流量分配。张西和等[6]利用动网格技术和流体力学基本原理对罩式炉内循环气体的冷态流场进行计算，结果显示沿罩式炉高度方向无论在钢卷的芯部还是外侧，循环气体的流速沿高度方向均是逐渐衰减的。总而言之，由于对流板结构复杂，钢卷尺寸规格各异，内罩内的氢气流动情况及其复杂，准确确定炉内流动情况特别困难。本文作者通过CFD商业软件建立三维模型，充分考虑对流板的复杂结构，对罩式退火炉内罩内保护气体的流场分布进行了数值模拟，得到保护气体流量的分配规律及对各处努塞尔数(Nu)数的影响。采用Monte Carlo法计算了炉内辐射换热角系数，最终采用辐射换热组网络图法求解得到了炉内各处的辐射热流密度，为钢卷温度场的求解提供准确的边界条件。采用热电偶插片的实验方法，对钢卷退火过程内部温度进行现场测试，为传热数学模型的调试和验证提供可靠的参数。目前对于加热罩传热过程的处理，采用一维径向导热方程，而内罩则是视为温度均匀一致，采用集总参数法进行处理[7-10]。
1 传热过程数学模型
1.1 数学模型的建立
考虑到钢卷和罩式退火炉结构的轴对称性，根据退火过程中钢卷的受热特点，进行如下的简化。
1) 钢卷内部无热源；
2) 钢卷边界与外界换热符合第三类边界条件；
3) 忽略钢卷沿环向的导热，即只考虑钢卷沿径向和轴向的导热；
根据钢卷的受热分析情况和物理模型情况，建立钢卷二维轴对称非稳态导热微分方程。物理模型如图1所示。由于钢卷径向由多层带钢组成，带钢层间存在气隙层，使径向热导率小于轴向热导率，且随温度有一定的变化，因此描述钢卷温度分布规律的方程应为二维轴对称、各向异性的导热微分方程，如式(1)所示。
图1 钢卷导热的物理模型
Fig. 1 Physical model of steel coil conduction
(1)
钢卷导热微分方程的定解条件主要包括：几何条件(如图1所示)、初始条件和边界条件。
初始条件为：
T(r，z)=T0(t=0，Ri≤r≤Ro，0≤z≤H) (2)
边界条件为：
(r=Ri，0≤z≤H) (3)
(r=Ro，0≤z≤H) (4)
(z=0，Ri≤r≤Ro) (5)
(z=H，Ri≤r≤Ro) (6)
式中：T为钢卷温度，K；r和z分别为径向坐标和轴向坐标，m；hi，ho，hb和ht分别为钢卷内壁、外壁、上侧壁和下侧壁的换热系数，W/(m-2·K)；Tg为保护气体温度，K；qRi，qRo，qRb和qRt分别为钢卷内壁、外壁、上侧壁和下侧壁的辐射换热热流密度，W/m-2；λr和λz分别为钢卷径向和轴向的导热系数，W/(m·K)；Ri和Ro分别为钢卷内径和外径，m；H为钢卷高度(带钢宽度)，m。
影响钢卷温度场的重要参数包括钢卷径向等效导热系数、对流换热系数和辐射热流密度，因此准确确定这几个参数是求解钢卷温度场的关键。
1.2 径向等效导热系数
径向等效导热系数采用如下公式计算[11]。
(7)
式中：s为带钢厚度，m；b为带钢层与层之间的间隙，m；λs为带钢的导热系数，W/(m·K)；λF为保护气体的导热系数，W/(m·K)；Tm为相邻两层带钢的温度，K；ε为带钢黑度；σ为斯蒂芬.玻尔兹曼常数，W/(m2.K4)；σP为带钢表面的粗糙度，μm；tanθ为带钢表面形状的平均斜度；P为钢卷打卷张应力，MPa；H为接触固体中材料较软者的硬度，MPa；φ=P/(H+P)。
1.3 对流换热系数
钢卷位于炉台上方，炉台下方循环风机处的流场十分复杂，如图2所示，但对钢卷区域的气流分布影响不大(因为进气口和出气口保护气体流量相同)，所以只模拟炉台以上的区域。由于内罩和内部钢卷的布置是轴对称的，对流板沿周向具有周期性，由20个相同的结构组成，如图3所示。因此为了减少网格数，提高计算效率，取整体结构的1/20，采用周期性边界条件进行模拟计算，结构如图4所示。
影响炉内对流换热系数的因素很多，包括保护气体的流速，炉内几何结构和材质的热物性等，每种因素又有多种变化，试验量非常大，因此采用正交试验法进行研究，按照三因素三水平的设计思想，通过改变保护气体流量、钢卷的外径及炉内温度设计9个物理模型进行模拟，正交表头设计方案如表1所示。
1.3.1 物理模型的选择
计算中选取Realizable k-ε湍流模型作为流动模型，差分格式采用一阶迎风格式，假设内罩内氢气处于定温，氢气的热物理性质为其温度下的物理性质。
图2 罩式退火炉结构
Fig. 2 Structure of bell-type annealing furnace
图3 钢卷卷间对流板的几何结构
Fig. 3 Geometric structure of convection plate between steel coil
图4 内罩内保护气体流通区域几何模型
Fig. 4 Geometric model of protection gas circulation region inside inner cover
表1 正交表头设计方案
Table 1 Design scheme of orthogonal table
1.3.2 边界条件的设定
进口边界：给定氢气的质量流量(kg/s)和温度(K)；出口边界：压力出口，相对压力为0 Pa；壁面条件为无滑移绝热壁面。
1.4 辐射换热模型
辐射换热角系数是辐射换热求解的关键，本文采用Monte Carlo法计算各个表面之间的角系数。Monte Carlo法是20世纪50年代发展起来，在求解辐射换热问题时(包括：直接求解辐射换热方程、求解辐射传递系数、表面之间的角系数等)，不需要进行太多的简化假设，不仅适合于多维几何问题，也适合于非灰气体介质[12-14]。
为了获得各个表面上辐射热流密度的分布规律，如图5所示，将钢卷内外表面、内罩内表面、顶层钢卷上表面等划分为多个单元面并进行编号。采用Monte Carlo方法计算各单元面之间的角系数，得到图6所示的结果。通过热辐射网络图法，求解有效辐射方程组，并进一步得出各个单元面的辐射换热热流密度，进而代入边界条件方程式(3)~(6)。
图5 内罩截面图
Fig. 5 Inner cover cross section
图6 三维空间角系数fi-j计算结果
Fig. 6 Calculated Result of three dimensional space angle factor fi-j
2 方程的数值求解
钢卷传热的计算区域网格划分采用外节点法，如图7所示，建立圆柱轴对称坐标系，取一个弧度的中心角所包含的范围作为研究对象[15]，推导得到的离散方程式为
(8)
系数简化整理得：
apTp=aETE+aWTW+aNTN+aSTS+b (9)
采用交替方向块迭代法进行扫描计算，则可以加快收敛速度，即先逐行(或逐列)进行一次扫描，再逐列(或逐行)进行一次扫描，两次全场扫描组成一轮迭代。具体公式为：
apTp(n+1/2)=aETE(n+1/2)+aWTW(n+1/2)+ aNTN(n)+aSTS(n)+b (10)
apTp(n+1)=aNTN(n+1)+aSTS(n+1)+ aETE(n+1/2)+aWTW(n+1/2)+b (11)
图7 圆柱轴对称坐标网格系统
Fig. 7 Axisymmetric coordinate grid system
3 模拟结果分析
3.1 对流换热系数
正交试验的极差分析可以确定各因素对试验指标综合影响的次序，极差值最大的为主要因素，极差值较小的为次要因素。表2列出了试验结果的极差分析，从表2可以明显看出流量对炉内各处的对流换热系数影响最大，因此进一步考察流量对于图4所示的编号为①~④各处的流量分配的影响，①~③依次为从上到下每个对流板的编号，④为最顶部空间的编号。结果如表3所示。
从表3可以看出：大部分循环气体从罩内最上部空间(通道④)流过，从对流板处流过的气体较少。随着循环风量的增加，从通道④流过的流量所占比例变化不大。对于气体工况，如钢卷的卷数和尺寸发生变化时，也得到了相近的结果。
表2 正交试验结果的极差分析
Table 2 Extreme difference analysis of orthogonal experimental results
表3 不同风量条件下各通道内气体流量分配比例
Table 3 Gas flow distribution ratio of each channel under different flow conditions
显然，保护气体流量的变化对炉内各处的流量分配影响较小，可以忽略，在进行退火过程温度场模拟计算时，仅需计算一次得到各处的流量分配比例，就可获得任意入炉气体流量下各处的流量。
对流换热系数h=Nu·λ/l，其中λ为流体的导热系数，l为几何特征长度。
显然，在一定温度和一定结构下，求得Nu数就可获得对应的对流换热系数。
对于强制对流换热，Nu=f(Re,Pr)=aRemPrn，当温度一定时，普朗特数Pr一定，因此Nu只与雷诺数Re有关。而Re=v·l/υ，其中v为流体速度，等于Q/A, Q为流量，A为面积，l为几何特征长度，υ为运动黏度。因此，在一定温度和一定几何结构下，Nu与流体速度成正比，而流体速度与保护气体流量成正比，因此，确定Nu与保护气体流量的关系，就可以准确求得炉内结构一定的情况下各处对流换热系数与流量的关系。
图8所示为对流板上、下表面各处Nu随气体流量的变化，图8中虚线为第2个对流板上、下表面Nu±15%的误差限，显然，第1个和第3个对流板各处Nu基本上落在误差限内，因此拟合得到的第2个对流板各处Nu与流量Q的关系式，同样可以适用于第1个和第3个对流板。b与A，l，υ和Pr有关，炉内结构一定和温度一定时为常数。
对流板下表面：
Nu=f(Q，A，l，υ，Pr)=bQ0.83 (12)
对流板上表面：
Nu=f(Q，A，l，υ，Pr)=bQ0.78 (13)
图9所示为钢卷芯部和外部各处Nu随气体流量的变化，图9中虚线为第3个钢卷芯部和外部Nu±15%的误差限。图9(a)所示为钢卷芯部，第1个和第2个钢卷各处Nu基本上落在误差限内，因此拟合得到的第3个钢卷各处Nu与流量Q的关系式，同样可以适用于第1个和第2个钢卷。由图4可知：第4个钢卷芯部为气体入口，由于入口效应Nu较大，对于第4个钢卷芯部需要进行入口效应的修正。图9(b)所示为钢卷外部，保护气体从第1个钢卷外部进入炉内，但是由于保护气体经过炉台，因此此处入口效应较弱，仅在低保护气体流量下作用明显，其余基本都落在第3个钢卷的±15%误差限内，因此拟合得到的公式对于其余3个钢卷外部同样适用。c与A，l，υ和Pr有关，炉内结构一定和温度一定时为常数。
钢卷芯部：
Nu=f(Q，A，l，υ，Pr)=cQ0.87 (14)
钢卷外部：
图8 对流板表面各处Nu分布
Fig. 8 Distribution of Nu throughout convection plate surface
图9 钢卷内外表面各处Nu分布
Fig. 9 Distribution of Nu inside and outside steel coil surface
Nu=f(Q，A，l，υ，Pr)=cQ0.70 (15)
3.2 对流换热热流和辐射换热热流
图10所示为辐射换热量、对流换热量和氢气温度随时间的变化曲线。由图10可知：加热初始阶段，对流换热量急剧增加，一直持续到400 ℃左右，此时对流换热最大值约为700 kW，这是因为刚开始加热时，气体温度上升比较快，气体与钢卷温差大。随着加热过程的进行，对流换热量又迅速下降，这是因为钢卷温度上升显著，而气体温度上升相对平缓，使二者温差减小。随着温度的升高，辐射换热量也逐渐升高，到保温阶段开始，内罩温度不变，钢卷温度继续增加，使得辐射换热量开始下降。总的来说，对流换热量始终大于辐射换热量。图11所示为对流换热量与辐射换热量的比值直观地描述了这2个量相差的量级。这个最大值出现在加热的初始阶段，其比值的最大值为65。随着温度的上升，辐射热流密度和对流热流密度逐渐接近，在加热后期和保温阶段基本上处于同一个数量级。这也说明，炉内对流换热系数的准确性关系到钢卷温度预测的准确性，全氢罩式退火炉中，对流换热是主要换热方式。
图10 换热量和氢气温度随时间的变化曲线
Fig. 10 Heat transfer and hydrogen temperature as a function of time
图11 对流换热量与辐射换热量的比
Fig. 11 Ratio of convection heat transfer and radiation heat transfer
4 数学模型验证
4.1 钢卷温度场
用第1卷和第4卷的实验实测钢卷温度曲线与模型计算的钢卷温度曲线比较，对模型进行了验证，对比结果如图12所示。由图12可见：卷1的芯部温度预测中，相对误差以90%以上的概率落在3%以内，外部温度预测中，相对误差以85%的概率落在3%以内，以90%以上的概率落在10%以内；卷4的芯部温度预测中，相对误差以70%以上的概率落在3%以内，外部温度预测以80%以上的概率落在3%以内，以90%以上的概率落在10%以内。
卷1芯部在加热阶段末期误差较大的原因是测试时热电偶与钢卷接触不紧密，导致测量值无法真实反映当前位置的测点温度，总而而言模型计算曲线和实验数据曲线基本上吻合较好，说明此数学模型是准确、可靠地，可以用于预测钢卷的温度场分布，指导实际生产。
图12 模型计算和实测退火曲线比较
Fig. 12 Comparison of calculated data and experimental data
4.2 氢气温度
将模型计算得到的氢气温度值与实验测得值进行比较，如图13所示。从图13可以看出：模拟值和实验值基本吻合；由图13(b)的相对误差分析可以明显看出：95%的相对误差在±1%以内，只在加热最开始和末期相对误差值较大，这是因为在加热最开始，内罩温度视为均匀一致存在误差，造成了氢气温度的不准确，而末期氢气温度的测量不准确。
图13 氢气温度的模拟值与实验值比较
Fig. 13 Comparison of hydrogen temperature calculated data and experimental data
5 结论
1) 通过CFD建立三维模型模拟炉内退火过程，通过正交试验法得到在所考察的因素内流量对各处的换热影响作用的结论，并通过进一步分析得到各处Nu数与流量的关系式。
2) 采用Monte Carlo法求解炉内参与辐射换热的各面角系数，通过热辐射网络图获得各面的辐射热流密度，为温度场的求解提供精准的边界条件。
3) 通过钢卷及氢气的模拟计算和实测温度曲线的对比，说明该计算模型正确，并具有一定的精度，可以用于钢卷罩式退火过程的仿真计算，对优化罩式退火工艺具有一定的指导意义。
参考文献：
[1] Wu M, Chen Y, Pangal K, et al. High-performance polysilicon thin film transistor on steel substrates[J]. J Non-Crystalline Solids, 1999, 266: 1284-1288.
[2] Lochner H, Schweiger G. Annealing cold rolled strip in HICON/H2 bell annealers[J]. Iron and Steel Engineers, 1998, 4: 45-51.
[3] 李卫杰, 陈继东, 张向, 等. 罩式炉内循环保护气体流动特性的模型化研究[J]. 华中科技大学学报(自然科学版), 2009, 37(5): 116-119.
LI Weijie, CHEN Jidong, ZHANG Xiang, et al. Modeling study for flow characteristic of circulation protective gas in bell annealers[J]. Huazhong University of Science and Technology (Natural Science Edition), 2009, 37(5): 116-119.
[4] 黄文文. 全氢罩式炉料室空间流动特性研究[D]. 武汉: 华中科技大学能源与动力工程学院, 2011: 27-34.
HUANG Wenwen. Study on flow characteristic of the work base for high performance hydrogen bell-type annealers[D]. Wuhan: Huazhong University of Science and Technology. School of Energy and Power Engineering, 2011: 27-34.
[5] 林林. 全氢罩式退火炉退火过程的研究[D]. 北京: 北京科技大学机械工程学院, 2003: 41-50.
LIN Lin. Study of annealing thermal process in high performance hydrogen bell-type annealing furnace[D]. Beijing: University of Science and Technology Beijing. School of Mechanical Engineering, 2003: 41-50.
[6] 张西和, 陈光, 杨进. 动网格技术在罩式炉流场模拟中的应用[J]. 工业加热, 2008, 37(2): 37-40.
ZHANG Xihe, CHEN Guang, YANG Jin. Numerical simulation of gas flow in bell-type furnace using moving grids[J]. Industrial Heating, 2008, 37(2): 37-40.
[7] 莫春立, 詹志东, 李殿中, 等. 冷轧钢卷在罩式炉中退火过程中的温度场变化模拟计算[J]. 金属热处理, 2005, 30(5): 67-71.
MO Chunli, ZHAN Zhidong, LI Dianzhong, et al. Simulating the temperature evolution of the cold-rolled steel coil during annealing in bell-type furnace[J]. Heat Treatment of Metals, 2005, 30(5): 67-71.
[8] 孙晨, 姜泽毅, 刘志成, 等．全氢罩式炉钢卷退火过程在线数值仿真[J]. 工业加热, 2007, 36(2): 22-24.
SUN Chen, JIANG Zeyi, LIU Zhicheng, et al. On-line numerical simulation of steel coils’ annealing process in bell-type furnace with high performance hydrogen[J]. Industrial Heating, 2007, 36(2): 22-24.
[9] 高伟. 全氢炉退火过程钢卷温度场和应力场的耦合研究[D]. 武汉: 华中科技大学能源与动力工程学院, 2008: 21.
GAO wei. Simulate research of the coupling relationship between temperature field and thermal stress field of annealing process in full-hydrogen bell-type annealing furnace[D]. Wuhan: Huazhong University of Science and Technology. School of Energy and Power Engineering, 2008: 21.
[10] 方顺利. 全氢罩式炉退火工艺设备的仿真与分析[D]. 武汉: 华中科技大学能源与动力工程学院, 2012: 45-48.
FANG Shunli. The simulation and optimization of anneal technique of the full-hydrogen bell-type annealer furnace[D]. Wuhan: Huazhong University of Science and Technology. School of Energy and Power Engineering, 2012: 45-48.
[11] Park S J, Hong B H, Baik S C, et al. Finite element analysis of hot rolled coil cooling[J]. ISIJ Int, 1998, 38(11): 1262-1269.
[12] 卞伯绘. 辐射换热的分析与计算[M]. 北京: 清华大学出版社, 1988: 79-82.
BIAN Bohei. Analysis and calculation of radiation heat transfer[M]. Beijing: Tsinghua University press, 1988: 79-82.
[13] Zhou H C, Chen D L, Cheng Q. A new way to calculate radiative intensity and solve radiative transfer equation through using the Monte Carlo method[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2004, 83(3/4): 459-481.
[14] Mirhosseini M, Saboonchi A, View factor calculation using the Monte Carlo method for a 3D strip element to circular cylinder[J]. International Communications in Heat and Mass Transfer 2011, 38: 821-826.
[15] 陶文铨. 数值传热学[M]. 西安: 西安交通大学出版社, 2009: 86-88.
TAO Wenquan. Numerical heat transfer[M]. Xi’an: Xi’an Jiaotong University Press, 2009: 86-88.
(编辑 杨幼平)
收稿日期：2014-04-03；修回日期：2014-06-11
基金项目(Foundation item)：中央高校基本科研业务费专项资金资助项目(FRF-TP-14-021A1)(Project (FRF-TP-14-021A1) supported by the Fundamental Research Funds for the Central Universities)
通信作者：豆瑞锋，讲师，从事工业热设备传热、传质过程的模拟与实验等研究；：douruifeng@126.com
