熱成形工藝中界面換熱系數求解方法對比分析

常穎，唐行輝，王斌，盈亮，王存宇，趙坤民

?

熱成形工藝中界面換熱系數求解方法對比分析

常穎1，唐行輝1，王斌1，盈亮1，王存宇2，趙坤民1

(1. 大連理工大學汽車工程學院，工業裝備結構分析重點實驗室，遼寧大連，116024；2. 鋼鐵研究總院特殊鋼研究所，北京，100081)

根據高強度鋼熱成形特點，通過圓臺模具實驗，分別利用熱平衡法、Beck非線性估算法以及基于優化的有限元仿真法(FEM優化法)求解熱成形用22MnB5鋼在模具淬火下的實時界面換熱系數(IHTC)，并對3種求解方法的結果和適用性進行比較和仿真驗證。研究結果表明：Beck非線性估算法、熱平衡法與FEM優化法3種方法求解IHTC的平均誤差率分別為3.7%，7.5%和10.3%左右，Beck非線性估算法最準確；同時，得出馬氏體相變的發生會對樣件與模具間的IHTC有正的增益效果，壓強與IHTC間存在著近似冪函數關系。

熱成形；22MnB5鋼；界面換熱系數(IHTC)；熱平衡法；Beck非線性估算法；FEM優化

高強鋼熱成形技術是將高強度鋼板(如22MnB5)加熱至均勻奧氏體化后，在高溫下沖壓成形并在模具內以大于27℃/S的降溫速率連續冷卻淬火從而獲得較高強度和機械性能零部件的工藝過程。運用此技術生產的汽車零部件在車身領域得到了應用[1]。零件和模具間的界面換熱系數表征了高強度鋼熱成形沖壓淬火階段熱鋼板與水冷模具間熱量交換的劇烈程度，將直接決定了熱鋼板的降溫速率和零件溫度場的變化情況，進而間接影響到成形后零件的強度與硬度。故高強鋼熱成形過程中零件與模具間的界面換熱系數的精確求取將對成形零件溫度場和零件性能的預測起到至關重要的作用。界面換熱系數取決于接觸表面載荷、溫度(模具和零件)、表面粗糙程度和材料熱物理性能等諸多復雜因素[2]。針對熱成形過程的界面換熱系數，Merklei等[3?4]通過試驗對給定壓強和溫差的接觸界面換熱系數進行了研究，并對不同間隙值的間隙傳熱系數進行了比較。Abdulhay等[5?6]通過實驗得到U型件實沖過程中不同部位估計2~30 MPa壓力范圍的接觸熱阻，進而得出接觸熱阻和壓強的近似指數函數關系。目前為止，用于瞬態界面換熱系數主要有3種獲取方法：一是基于恒定的模具溫度和Biot數較小的假設，引入牛頓冷卻定律，求取界面換熱系數，可稱為熱平衡法；二是根據實驗獲取的模具或零件內部溫度場的變化，逆向求解界面的熱流密度和界面溫度，從而求出界面換熱系數[2?9]，即反向熱傳導問題處理方法(IHCP)；三是基于優化的有限元仿真法，即FEM優化法，是工業中常用的方法，可以通過優化軟件指導在FEM 軟件中建立好的仿真模型，通過調整IHTC，使仿真溫度場與實驗值逼近吻合[10?11]。上述理論方法在熱成形過程的應用上主要局限有三點：

1) 由于熱成形淬火過程中模面溫度是實時變化的，所以牛頓冷卻定律無法直接應用，需要加以改進；

2) 反向熱傳導處理方法大都依賴于內部點的溫度采集，為了盡量減小熱傳導本身的延遲和滯后性，需要對溫度的采集精度和位置標定提出較高的要求；

3) FEM優化法一般只能得到淬火過程中的一個有效IHTC[12]，無法表征整個淬火過程中IHTC的實時變化情況。

因此，基于上述問題，本文作者主要研究：1) 建立圓臺模型，通過表面熔焊熱電偶的方式直接采集的模具表面溫度歷程，并利用自編的以Beck非線性估計法[13]為核心的求解程序，根據模面以下測點溫度反求得到的模具表面溫度進行對比，探討Beck的非線性估計法對熱成形過程中換熱系數估計的適用性與準確性；2) 將熱平衡法、Beck非線性估算法與FEM優化法進行對比，比較三者計算IHTC的準確性，以期獲取熱成形工藝中最佳的IHTC求解方式；3) 分析沖壓壓強因素對界面換熱系數的影響。

1 實驗

采用的實驗裝置如圖1(a)所示。可調壓機的公稱壓力為400 kN，模具材質為45號模具鋼，壓力平衡支座的作用是保證零件表面所受壓力均勻。下模具設有三根單線細絲K型熱電偶，分別采集距模面中心處2，4和6 mm處溫度(如圖1(b)所示)，在上模表面增加了表面熔焊的熱電偶，因熱電偶已與模具表面熔焊，磨平處理后，可以認為采集的溫度為有效模面溫度。試件采用22MnB5熱成形鋼,切割成直徑為75 mm、高為2 mm的圓柱樣件，熱電偶固定在板料中心附近。試件加熱至900 ℃后保溫約3 min，待完全均勻奧氏體化后迅速轉移，上模下行并保壓30 s后取出樣件。熱電偶采集的溫度信號通過采溫裝置MX100實時生成數據文檔，試驗過程如圖1(c)所示。

(a) 實驗裝置簡圖；(b) 熱電偶布置位置下垂量測試裝置示意圖；(c) 實驗沖壓模具

2 實驗結果與分析

2.1 界面換熱系數的3種估算方法

在界面傳熱的研究過程中，系統與外界的熱交換包括沖壓件和模具與空氣的熱輻射和熱對流換熱,但因其量級較小,可假設系統封閉；對于試件而言，其Biot數小于0.1[14]，可采用簡化的集中參數法進行分析，認為其厚度方向溫度不變。基于以上假設可將實驗模型可簡化成一維熱傳導問題。

2.1.1 熱平衡法

在忽略零件和模具與空氣的輻射和對流換熱情況下，根據能量守恒定律，零件損失的熱量等于模具吸收的熱量，故可推得變形后的牛頓冷卻公式：

式中：為零件比熱容；w為模具溫度；為接觸面積；為對應時刻的換熱系數；(t)為時刻的板料溫度，(t?1)為?1時刻的板料溫度。

事實上?時間內w是連續變化的，即使和保持不變，也難以求出微分方程的解析解。所以采用分段數值積分對連續變化的溫度歷程求解Δ對應的換熱系數，其求解算式為

式中：eq為Δ時段的等效比熱容。

22MnB5鋼馬氏體相變過程的eq表示為：

式中：為單位質量熱焓(潛熱)，J/kg；為馬氏體體積分數，可根據Koistinen?Marburger定律[15]計算：

式中：s為馬氏體轉變溫度；為相變動力因子。s和與材料相關，對于熱成形過程，研究表明s還與壓力相關，但針對本文的壓力研究范圍對應的實驗觀察結果，s取為400 ℃。對于普通鋼材，取為0.011。求得的等效熱容曲線如圖2所示。

1—比熱容；2—等效比熱容

2.1.2 Beck非線性估計法

對于圓形板料的熱沖壓實驗來說，不存在內熱源，板料與模具間的傳熱問題可以近似等效為一維傳熱問題，故對與常物性、無內熱源的一維非穩態導熱問題，可表示為

式中：為導熱系數；為材料密度；為材料的比定壓熱容；為時間；為圓臺下模沿軸線方向任意位置的溫度；為能量傳遞方向。

若給出初始時刻溫度分布的初始條件(initial condition)，以及規定了邊界上的熱流密度，則可以對溫度場進行求解，此類求解方式稱為第二類邊界條件求解，也稱為Neumann條件求解[16]。

Beck非線性估算法處理反向熱傳導問題是最為經典和應用最為廣泛的方法之一，其實質為一種優化方法與利用第二類邊界條件求解溫度場的經典組合，其優化目標函數是：

式中：為分段內計算點數量；為計算中所需測點數量。

求解思路為：給定0時刻熱流密度0，求出后“1，2，…，”時刻測點溫度，求出溫度對熱流密度的靈敏系數，修正熱流密度直至滿足收斂條件。通過調用Matlab中的Pdepe函數求解工具求解一維熱傳導偏微分方程，結合自編的迭代優化程序，求出各時段的優化收斂的熱流密度。并求解出相應熱流密度下不同時刻的模面溫度，最后根據牛頓換熱定律可得該時刻所對應的表面換熱系數。如圖3所示。

圖3 Beck非線性估算法流程

2.1.3 FEM優化法

優化軟件Isight指導在FEM軟件Abaqus中建立好的仿真模型，通過給定IHTC的不斷試錯過程，使仿真溫度場與實驗值收斂吻合。分析中采用的優化方法為“Pointer”方法，優化進行時，“Pointer”方法能夠自行控制步長、迭代次數等優化參數實現最高的優化效率。

誤差函數為：

式(7)所示為零件溫度場仿真值與實測值最小方差，式(8)所示為模具已知測點位置溫度場仿真值與實測值的最小方差。

FEM優化法計算流程如圖4所示。

圖4 FEM優化流程

2.2 實驗結果和IHTC的求解

圖5所示為沖壓壓強1 MPa下采集的溫度曲線。由圖5可以看出：在極短時間內(約2.5 s)，樣件溫度由750 ℃降到約400 ℃(s點)，降溫速率突然減緩，此時發生馬氏體相變過程，釋放潛熱。模面溫度快速上升至約150 ℃后趨于穩定而后略微上升出現拐點，同樣由相變的原因來解釋。

1—試件溫度；2—測得模面溫度；3—測得模面下2 mm處溫度；4—測得模面下4 mm處溫度；5—測得模面下6 mm處溫度

3種方法所計算的界面換熱系數曲線如圖6所示。由圖6可見：熱平衡法與Beck估算法所得的IHTC曲線整體趨勢是一致的。在樣件和模具剛接觸時，樣件和模具逐漸進入理想的接觸狀態，然后隨其溫度下降，溫差縮小，熱傳導的驅動力不斷減小，IHTC減小。當溫度降到馬氏體轉變溫度s點(400 ℃附近)時相變發生，由此所導致的體積膨脹、相變塑性、相變表面浮凸[17]等因素綜合影響下會在一定程度上增大零件與模具的實際接觸面積，而相變過程中材料熱物性參數的改變、潛熱釋放等原因也會導致熱傳導驅動力增大，兩者綜合作用下最終導致了相變對IHTC正的增益效果。FEM優化法進行對沖壓壓強1 MPa下模具與板料間有效IHTC的估算結果為1.728 kW/(m2·K)。

1—界面換熱系數(Beck非線性估算法)；2—界面換熱系數(熱平衡法)；3—界面換熱系數(FEM優化法)

2.3 試驗結果分析

將熱平衡法和Beck非線性估算法所求換熱系數帶入Deform 3D軟件中，保證其他材料的力學參數與熱物性參數均為一致，帶入實驗的邊界條件進行瞬態熱傳導計算并模擬零件降溫過程。

3種方法下樣件淬火溫度仿真值與測量值之間的收斂情況如圖7(a)所示，利用Beck非線性估算法，反算求解的模具溫度與實測溫度吻合情況如圖7(b)所示。從整體上說，Beck非線性估算法能夠在一定程度上很好地反映模面溫度的變化，并且樣件淬火溫度變化的仿真值與實測值最為貼近。

(a) 3種方法下溫度仿真值與測量值之間的收斂情況；(b) Beck非線性估算法求解模溫度場吻合情況

圖8所示為誤差分析結果。由圖8可見：Beck非線性估算法能很好地反映出馬氏體相變對樣件和模具溫度的變化情況，吻合度較好，誤差控制在3.7%左右；而熱平衡法雖然也能近似的得出模具與樣件的實際溫度變化情況，但可以看出，用等效熱容來考慮相變的方式使得零件溫度下降到Ms點(圖8(a)中第4秒位置附近)以后，相對誤差略大，總體吻合誤差為7.5%左右；而FEM優化法，所求得的IHTC為1個等效值，仿真結果較前兩者，與測量值吻合度誤差超過了10%。

(a) 3種方法下零件淬火仿真值與真實值間的絕對誤差分析；(b) 多種誤差分析方法

2.4 沖壓壓強對界面換熱系數的影響分析

沖壓壓強是界面換熱系數的一個重要的因素。20世紀60年代，Cooper等[18?19]基于表面形貌的數學描述對接觸熱阻的計算和測試方法進行了研究。在不同表面進入名義上的接觸狀態后，其真實的接觸只發生在表面凸起點上，在凸起點的間隙存在大量空氣雜質，表面間的真實接觸面積也僅為名義接觸面積的1%~2%[20]，沖壓壓強越大，勢必會在一定程度上會增大接觸微凸體的變形，致使實際接觸面積增大，接觸熱阻減小，界面換熱系數也相應的增大。

試驗中通過調節液壓機沖壓壓力以獲取不同的壓強梯度，通過Beck非線性估算法得到了如圖9(a)所示不同壓強下樣件與45號模具鋼間IHTC隨溫度變化的實時曲線；為了更加清晰地探究IHTC與壓強間的關系，將每條IHTC曲線取均值處理，得到如圖9(b)所示的均值IHTC與壓強之間的關系。由圖9(b)可看出：均值IHTC整體上是隨沖壓壓強的增大而增大的。經過擬合，得出IHTC與壓強間存在近似的冪函數關系，按照均值IHTC與沖壓壓強的關系式，選用各種曲線擬合，得此關系擬合可靠度最高。

(a) 不同沖壓壓強下IHTC隨溫度變化曲線；(b) 不同壓強下均值IHTC變化趨勢

3 結論

1) Beck非線性估算法的求解精度最高，熱平衡法求解精度次之，而FEM優化法由于所求得的IHTC為一個有效值，根據仿真結果誤差最大。

2) 馬氏體相變對換熱會存在一定的影響，且有正向增益的效果。

3) 沖壓壓強是影響IHTC的重要因素，且與均值IHTC間存在近似冪函數關系。

4) 由于高強鋼板熱成形過程載荷較大，零件溫度變化劇烈，工業中應用廣泛的FEM優化法的精準度還有待進一步提高。

[1] Turetta A, Bruschi S, Ghiotti A. Investigation of 22MnB5 formability in hot stamping operations[J]. Journal of Materials Processing Technology, 2006, 177(1/2/3): 396?400.

[2] Ikeuchi K, Yanagimoto J. Valuation method for effects of hot stamping process parameters on product properties using hot forming simulator[J]. Journal of Materials Processing Technology, 2011, 211(8): 1441?1447.

[3] Merklein M, Lechler J, St?hr T. Investigation on the thermal behavior of ultra high strength born manganese steels within hot stamping[J]. International Journal of Material Forming, 2009, 2(S1): 259?262.

[4] Bosetti P. Interlaboratory comparison for heat transfer coefficient identification in hot stamping of high strength steels[J]. International Journal of Material Forming, 2010, 3(1): 817?820.

[5] Abdulhay B, Bourouga B, Dessain C, et al. Experimental study of heat transfer in hot stamping process[J]. International Journal of Material Forming, 2009, 2(1): 255?257.

[6] Abdulhay B, Bourouga B, Dessain C. Thermal contact resistance estimation: influence of the pressure contact and the coating layer during a hot forming process[J]. International Journal of Material Forming, 2012, 5(3): 183?197.

[7] Khandelwal M, Mench M M. Direct measurement of through-plane thermal conductivity and contact resistance in fuel cell materials[J]. Journal of Power Sources, 2006, 161(2): 1106?1115.

[8] Caron E, Daun K J, Wells M A. Experimental characterization of heat transfer coefficients during hot forming die quenching of boron steel[J]. Metallurgical and Materials Transactions B, 2013, 44(2): 332?343.

[9] Bai Q, Lin J, Zhan L, et al. An efficient closed-form method for determining interfacial heat transfer coefficient in metal forming[J]. International Journal of Machine Tools and Manufacture, 2012, 56: 102?110.

[10] 郭志鵬, 熊守美, 曺尚鉉, 等. 鋁合金壓鑄過程鑄件/鑄型界面換熱行為的研究: Ⅰ. 實驗研究和界面換熱系數求解[J]. 金屬學報, 2007, 43(11):1149?1154. GUO Zhipeng, XIONG Shoumei, CAO Shangxuan, et al. Study on heat transfer behavior at metal/die interface in aluminum Alloy die casing process: I. Experimental study and determination of the interfacial heat transfer coefficient[J]. Acta Metallurgica Sinica, 2007, 43(11): 1149?1154.

[11] 郭志鵬, 熊守美, 曹尚鉉, 等. 熱傳導反算模型的建立及其在求解界面熱流過程中的應用[J]. 金屬學報, 2007, 43(6): 607?611. GUO Zhipeng, XIONG Shoumei, CAO Shangxuan, et al. Development of an inverse heat transfer model and its application in the prediction of the interfacial heat flux[J]. Acta Metallurgica Sinica, 2007, 43(6): 607?611.

[12] Wendelstorf R, Spitzer K H, Wendelstorf J. Effect of oxide layers on spray water cooling heat transfer at high surface temperatures[J]. International Journal of Heat and Mass Transfer, 2008, 51(19/20): 4892?4901.

[13] Krishnan M, Sharma D G R. Determination of the interfacial heat transfer coefficient h in unidirectional heat flow by beck's non linear estimation procedure[J]. International Communications in Heat and Mass Transfer, 1996, 23(2): 203?214.

[14] 梁昆淼. 數學物理方程[M]. 北京: 高等教育出版社, 2002: 60?65 LIANG Kunmiao. Equations of mathematical physics[M]. Beijing: Higher Education Press, 2002: 60?65.

[15] Koistinen D P, Marburger R E. A general equation prescribing the extent of the austenite-martensite transformation in pure iron-carbone alloys and plain carbon steels[J]. Acta Metallurgica, 1959, 1(7): 59?60.

[16] 楊世銘, 陶文銓. 傳熱學[M].北京: 高等教育出版社, 2006: 69?166. YANG Shiming, TAO Wenquan, Heat transfer theory[M]. Beijing: Higher Education Press, 2006: 69?166.

[17] 劉宗昌, 計云萍, 任慧平. 近幾年馬氏體相變研究的進展[J]. 熱處理, 2011, 26(6): 1?8. LIU Zongchang, JI Yunping, REN Huiping. Research progress of martensite phase transformation in recent years[J]. Heat Treatment, 2011, 26(6): 1?8.

[18] Cooper M G, Mikic B B, Yovanovich M M. Thermal contact conductance[J]. International Journal of Heat and Mass Transfer, 1969, 12(3): 279?300.

[19] Yovanovich M M, Fenech H. Thermal contact conductance of nominally flat, rough surfaces in a vacuum environment[C]//Thermo Physics and Temperature Control of Spacecraft and Entry Vehicles. Monterey, California, 1966: 773?794.

[20] 趙蘭萍, 徐烈, 李兆慈, 等.固體界面間接觸導熱的機理和應用研究[J]. 低溫工程, 2000, 1(4): 29?34. ZHAO Lanping, XU Lie, LI Zhaoci, et al. Study on mechanism and application of thermal contact heat transfer between solid interfaces[J]. Cryogenics, 2000, 1(4): 29?34.

(編輯 楊幼平)

Comparison of methods for calculating interfacial heat transfer coefficient in hot forming process

CHANG Ying1, TANG Xinghui1, WANG Bin1, YING Liang1, WANG Cunyu2, ZHAO Kunmin1

1. State Key Laboratory of Structural Analysis for Industrial Equipment, School of Automotive Engineering, Dalian University of Technology, Dalian 116024, China; 2. Institute for Special Steels, Central Iron & Steel Research Institute, Beijing 100081, China

The real-time interfacial heat transfer coefficient (IHTC) of 22MnB5 steel in hot forming process was calculated by means of heat balance method, Beck’s non-linear estimation method and optimization-based finite element simulation. The results show an average error of 3.7% for Beck’s non-linear estimation method, 7.5% for heat balance method, and 10.3% for finite element simulation. The positive effect of martensitic transformation on interfacial heat transfer coefficient is obtained. Additionally, the relationship between the pressure and the heat transfer coefficient is approximated by a power function.

hot forming; 22MnB5 steel; interfacial heat transfer coefficient(IHTC); heat balance method; Beck’s non-linear estimation method; FEM optimization

10.11817/j.issn.1672-7207.2015.04.006

TG162.83

A

1672?7207(2015)04?1208?07

2014?04?13；

2014?06?20

國家自然科學基金資助項目(10932003，51101036)；遼寧省自然科學基金和教育部留學回國人員科研啟動基金聯合資助項目(優秀人才培育2014028001)(Projects (10932003，51101036) supported by the National Natural Science Foundation of China; Project (2014028001) supported by the Natural Science Foundation of Liaoning Province and China MOE Scientific Research Funds for Returned Overseas Scholar)

趙坤民，教授，博士生導師，從事金屬熱處理研究；E-mail：kmzhao@dlut.edu.cn