罩式退火爐熱過程數學模型與實驗驗證

楊培培，溫治,2，豆瑞鋒

?

罩式退火爐熱過程數學模型與實驗驗證

楊培培1，溫治1,2，豆瑞鋒1

(1. 北京科技大學機械工程學院，北京，100083；2. 北京科技大學冶金工業節能減排北京市重點實驗室，北京，100083)

考慮對流板復雜的幾何結構，建立三維模型，采用正交試驗法模擬獲得流量、外徑和爐內溫度對爐內各處努塞爾數分布的影響，得到流量在考察的3種因素內流量的影響作用最強的結論。接著進一步分析隨流量的變化，并擬合得到各處與流量的關系式；為了解決爐內復雜幾何形狀條件下的輻射換熱角系數求解問題，采用Monte Carlo方法建立爐內輻射換熱角系數求解模型，并最終采用熱輻射網絡圖法求解爐內輻射換熱。通過現場進行熱電偶插片實驗驗證本文所建立模型的準確性。研究結果表明：模擬計算得到的溫度與現場實測值基本一致，證明了模型的準確可靠性。

正交試驗法；蒙特卡洛法；熱電偶插片實驗；罩式退火爐；數學模型

帶鋼是以鋼卷的形式在全氫罩式退火爐中退火的，影響傳熱的主要參數是徑向等效導熱系數、對流換熱系數和輻射熱流密度。目前國內外許多學者[1?2]針對接觸傳熱問題展開了深入廣泛的研究，得到了較為準確的徑向等效導熱系數的理論計算公式。相比而言，對內罩空間內各處對流換熱系數和輻射換熱的研究，則進行了過多的簡化假設，使結果的計算精度降低。保護性循環氣體在爐內的流量和流速分配直接影響表面對流換熱，因此，準確計算爐內保護氣體流量分配是準確計算爐內對流換熱的關鍵。李衛杰等[3?4]采用流動阻力分析和預測校正方法，將導流板結構簡化為一定厚度的中間層和部分界面流通的上下兩層鍥形槽層，對爐內循環氣體流量分配進行了模型化并對影響流動分配的相關影響因素進行了敏感性測試。林林等[5]通過對氫氣流動阻力分析，采用預測校正法，在爐型尺寸、對流板結構尺寸、風機循環風量一定的前提下，將對流板結構進行簡化，進而確定了不同鋼卷堆垛數、不同鋼卷尺寸時爐內各處的流量分配。張西和等[6]利用動網格技術和流體力學基本原理對罩式爐內循環氣體的冷態流場進行計算，結果顯示沿罩式爐高度方向無論在鋼卷的芯部還是外側，循環氣體的流速沿高度方向均是逐漸衰減的?？偠灾?，由于對流板結構復雜，鋼卷尺寸規格各異，內罩內的氫氣流動情況及其復雜，準確確定爐內流動情況特別困難。本文作者通過CFD商業軟件建立三維模型，充分考慮對流板的復雜結構，對罩式退火爐內罩內保護氣體的流場分布進行了數值模擬，得到保護氣體流量的分配規律及對各處努塞爾數()的影響。采用Monte Carlo法計算了爐內輻射換熱角系數，最終采用輻射換熱組網絡圖法求解得到了爐內各處的輻射熱流密度，為鋼卷溫度場的求解提供準確的邊界條件。采用熱電偶插片的實驗方法，對鋼卷退火過程內部溫度進行現場測試，為傳熱數學模型的調試和驗證提供可靠的參數。目前對于加熱罩傳熱過程的處理，采用一維徑向導熱方程，而內罩則是視為溫度均勻一致，采用集總參數法進行處理[7?10]。

1 傳熱過程數學模型

1.1 數學模型的建立

考慮到鋼卷和罩式退火爐結構的軸對稱性，根據退火過程中鋼卷的受熱特點，進行如下的簡化。

1) 鋼卷內部無熱源；

2) 鋼卷邊界與外界換熱符合第三類邊界條件；

3) 忽略鋼卷沿環向的導熱，即只考慮鋼卷沿徑向和軸向的導熱；

根據鋼卷的受熱分析情況和物理模型情況，建立鋼卷二維軸對稱非穩態導熱微分方程。物理模型如圖1所示。由于鋼卷徑向由多層帶鋼組成，帶鋼層間存在氣隙層，使徑向熱導率小于軸向熱導率，且隨溫度有一定的變化，因此描述鋼卷溫度分布規律的方程應為二維軸對稱、各向異性的導熱微分方程，如式(1) 所示。

鋼卷導熱微分方程的定解條件主要包括：幾何條件(如圖1所示)、初始條件和邊界條件。

初始條件為：

(，)=0(=0，i≤≤o，0≤≤) (2)

邊界條件為：

(=i，0≤≤) (3)

(=o，0≤≤) (4)

(=0，i≤≤o) (5)

(=，i≤≤o) (6)

式中：為鋼卷溫度，K；和分別為徑向坐標和軸向坐標，m；i，o，b和t分別為鋼卷內壁、外壁、上側壁和下側壁的換熱系數，W/(m2?K)；g為保護氣體溫度，K；qi，qo，qb和qt分別為鋼卷內壁、外壁、上側壁和下側壁的輻射換熱熱流密度，W/m2；λ和λ分別為鋼卷徑向和軸向的導熱系數，W/(m?K)；i和o分別為鋼卷內徑和外徑，m；為鋼卷高度(帶鋼寬度)，m。

圖1 鋼卷導熱的物理模型

影響鋼卷溫度場的重要參數包括鋼卷徑向等效導熱系數、對流換熱系數和輻射熱流密度，因此準確確定這幾個參數是求解鋼卷溫度場的關鍵。

1.2 徑向等效導熱系數

徑向等效導熱系數采用如下公式計算[11]。

式中：為帶鋼厚度，m；為帶鋼層與層之間的間隙，m；λ為帶鋼的導熱系數，W/(m?K)；λ為保護氣體的導熱系數，W/(m?K)；T為相鄰兩層帶鋼的溫度，K；為帶鋼黑度；為斯蒂芬.玻爾茲曼常數，W/(m2.K4)；σ為帶鋼表面的粗糙度，μm；tan為帶鋼表面形狀的平均斜度；為鋼卷打卷張應力，MPa；為接觸固體中材料較軟者的硬度，MPa；=/(+)。

1.3 對流換熱系數

鋼卷位于爐臺上方，爐臺下方循環風機處的流場十分復雜，如圖2所示，但對鋼卷區域的氣流分布影響不大(因為進氣口和出氣口保護氣體流量相同)，所以只模擬爐臺以上的區域。由于內罩和內部鋼卷的布置是軸對稱的，對流板沿周向具有周期性，由20個相同的結構組成，如圖3所示。因此為了減少網格數，提高計算效率，取整體結構的1/20，采用周期性邊界條件進行模擬計算，結構如圖4所示。

圖2 罩式退火爐結構

圖3 鋼卷卷間對流板的幾何結構

圖4 內罩內保護氣體流通區域幾何模型

影響爐內對流換熱系數的因素很多，包括保護氣體的流速，爐內幾何結構和材質的熱物性等，每種因素又有多種變化，試驗量非常大，因此采用正交試驗法進行研究，按照三因素三水平的設計思想，通過改變保護氣體流量、鋼卷的外徑及爐內溫度設計9個物理模型進行模擬，正交表頭設計方案如表1所示。

表1 正交表頭設計方案

1.3.1 物理模型的選擇

計算中選取Realizable?湍流模型作為流動模型，差分格式采用一階迎風格式，假設內罩內氫氣處于定溫，氫氣的熱物理性質為其溫度下的物理性質。

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 內罩截面圖

圖6 三維空間角系數fi?j計算結果

2 方程的數值求解

鋼卷傳熱的計算區域網格劃分采用外節點法，如圖7所示，建立圓柱軸對稱坐標系，取一個弧度的中心角所包含的范圍作為研究對象[15]，推導得到的離散方程式為

系數簡化整理得：

aT=aT+aT+aT+aT+b(9)

圖7 圓柱軸對稱坐標網格系統

采用交替方向塊迭代法進行掃描計算，則可以加快收斂速度，即先逐行(或逐列)進行一次掃描，再逐列(或逐行)進行一次掃描，兩次全場掃描組成一輪迭代。具體公式為：

aT(n+1/2)=aT(n+1/2)+aT(n+1/2)

aT(n)+aT(n)(10)

aT(n+1)=aT(n+1)+aT(n+1)

aT(n+1/2)+aT(n+1/2)(11)

3 模擬結果分析

3.1 對流換熱系數

正交試驗的極差分析可以確定各因素對試驗指標綜合影響的次序，極差值最大的為主要因素，極差值較小的為次要因素。表2列出了試驗結果的極差分析，從表2可以明顯看出流量對爐內各處的對流換熱系數影響最大，因此進一步考察流量對于圖4所示的編號為①~④各處的流量分配的影響，①~③依次為從上到下每個對流板的編號，④為最頂部空間的編號。結果如表3所示。

真的，什么革命，什么轉折點，什么歷史的劇變。想想陸地開墾那緩慢而艱難的過程，那無休無止而意義不明的過程——人類淤積的過程[2]9。

表2 正交試驗結果的極差分析

表3 不同風量條件下各通道內氣體流量分配比例

從表3可以看出：大部分循環氣體從罩內最上部空間(通道④)流過，從對流板處流過的氣體較少。隨著循環風量的增加，從通道④流過的流量所占比例變化不大。對于氣體工況，如鋼卷的卷數和尺寸發生變化時，也得到了相近的結果。

顯然，保護氣體流量的變化對爐內各處的流量分配影響較小，可以忽略，在進行退火過程溫度場模擬計算時，僅需計算一次得到各處的流量分配比例，就可獲得任意入爐氣體流量下各處的流量。

對流換熱系數=·/，其中為流體的導熱系數，為幾何特征長度。

顯然，在一定溫度和一定結構下，求得就可獲得對應的對流換熱系數。

對于強制對流換熱，=(,)=aRePr，當溫度一定時，普朗特數一定，因此只與雷諾數有關。而=·/，其中為流體速度，等于/,為流量，為面積，為幾何特征長度，為運動黏度。因此，在一定溫度和一定幾何結構下，與流體速度成正比，而流體速度與保護氣體流量成正比，因此，確定與保護氣體流量的關系，就可以準確求得爐內結構一定的情況下各處對流換熱系數與流量的 關系。

圖8所示為對流板上、下表面各處隨氣體流量的變化，圖8中虛線為第2個對流板上、下表面±15%的誤差限，顯然，第1個和第3個對流板各處基本上落在誤差限內，因此擬合得到的第2個對流板各處與流量的關系式，同樣可以適用于第1個和第3個對流板。與，，和有關，爐內結構一定和溫度一定時為常數。

(a) 對流板下表面；(b) 對流板上表面

對流板下表面：

=(，，，，)=0.83(12)

對流板上表面：

=(，，，，)=0.78(13)

圖9所示為鋼卷芯部和外部各處隨氣體流量的變化，圖9中虛線為第3個鋼卷芯部和外部±15%的誤差限。圖9(a)所示為鋼卷芯部，第1個和第2個鋼卷各處基本上落在誤差限內，因此擬合得到的第3個鋼卷各處與流量的關系式，同樣可以適用于第1個和第2個鋼卷。由圖4可知：第4個鋼卷芯部為氣體入口，由于入口效應較大，對于第4個鋼卷芯部需要進行入口效應的修正。圖9(b)所示為鋼卷外部，保護氣體從第1個鋼卷外部進入爐內，但是由于保護氣體經過爐臺，因此此處入口效應較弱，僅在低保護氣體流量下作用明顯，其余基本都落在第3個鋼卷的±15%誤差限內，因此擬合得到的公式對于其余3個鋼卷外部同樣適用。與，，和有關，爐內結構一定和溫度一定時為常數。

(a) 鋼卷芯部；(b)鋼卷外部

鋼卷芯部：

=(，，，，)=0.87(14)

鋼卷外部：

=(，，，，)=0.70(15)

3.2 對流換熱熱流和輻射換熱熱流

圖10所示為輻射換熱量、對流換熱量和氫氣溫度隨時間的變化曲線。由圖10可知：加熱初始階段，對流換熱量急劇增加，一直持續到400 ℃左右，此時對流換熱最大值約為700 kW，這是因為剛開始加熱時，氣體溫度上升比較快，氣體與鋼卷溫差大。隨著加熱過程的進行，對流換熱量又迅速下降，這是因為鋼卷溫度上升顯著，而氣體溫度上升相對平緩，使二者溫差減小。隨著溫度的升高，輻射換熱量也逐漸升高，到保溫階段開始，內罩溫度不變，鋼卷溫度繼續增加，使得輻射換熱量開始下降?？偟膩碚f，對流換熱量始終大于輻射換熱量。圖11所示為對流換熱量與輻射換熱量的比值直觀地描述了這2個量相差的量級。這個最大值出現在加熱的初始階段，其比值的最大值為65。隨著溫度的上升，輻射熱流密度和對流熱流密度逐漸接近，在加熱后期和保溫階段基本上處于同一個數量級。這也說明，爐內對流換熱系數的準確性關系到鋼卷溫度預測的準確性，全氫罩式退火爐中，對流換熱是主要換熱方式。

1—輻射換熱量；2—對流熱換量；3—氫氣溫度

圖11 對流換熱量與輻射換熱量的比

4 數學模型驗證

4.1 鋼卷溫度場

用第1卷和第4卷的實驗實測鋼卷溫度曲線與模型計算的鋼卷溫度曲線比較，對模型進行了驗證，對比結果如圖12所示。由圖12可見：卷1的芯部溫度預測中，相對誤差以90%以上的概率落在3%以內，外部溫度預測中，相對誤差以85%的概率落在3%以內，以90%以上的概率落在10%以內；卷4的芯部溫度預測中，相對誤差以70%以上的概率落在3%以內，外部溫度預測以80%以上的概率落在3%以內，以90%以上的概率落在10%以內。

(a) 卷1；(b) 卷1相對誤差；(c) 卷4；(d) 卷4相對誤差

卷1芯部在加熱階段末期誤差較大的原因是測試時熱電偶與鋼卷接觸不緊密，導致測量值無法真實反映當前位置的測點溫度，總而而言模型計算曲線和實驗數據曲線基本上吻合較好，說明此數學模型是準確、可靠地，可以用于預測鋼卷的溫度場分布，指導實際生產。

4.2 氫氣溫度

將模型計算得到的氫氣溫度值與實驗測得值進行比較，如圖13所示。從圖13可以看出：模擬值和實驗值基本吻合；由圖13(b)的相對誤差分析可以明顯看出：95%的相對誤差在±1%以內，只在加熱最開始和末期相對誤差值較大，這是因為在加熱最開始，內罩溫度視為均勻一致存在誤差，造成了氫氣溫度的不準確，而末期氫氣溫度的測量不準確。

(a) 氫氣溫度的實驗值和模擬值；(b) 相對誤差

5 結論

1) 通過CFD建立三維模型模擬爐內退火過程，通過正交試驗法得到在所考察的因素內流量對各處的換熱影響作用的結論，并通過進一步分析得到各處數與流量的關系式。

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/H2bell 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.

(編輯 楊幼平)

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)

Considering the complex structure characteristics of convection plate, three dimensional model and effect of flow rate, steel coil outer diameter and furnace temperature ondistribution were built using orthogonal test. The results show that flow rate has strongest influence on the range of consideration. Moreover,with the flow rate change was analyzed, and the relationship between flow rate andwas 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.

orthogonal test method; Monte Carlo method; thermocouple insert experiment; batch furnace; mathematical model

10.11817/j.issn.1672-7207.2015.04.045

TF062

A

1672?7207(2015)04?1518?09

2014?04?03；

2014?06?11

中央高?；究蒲袠I務費專項資金資助項目(FRF-TP-14-021A1)(Project (FRF-TP-14-021A1) supported by the Fundamental Research Funds for the Central Universities)

豆瑞鋒，講師，從事工業熱設備傳熱、傳質過程的模擬與實驗等研究；E-mail：douruifeng@126.com