金屬蜂窩夾芯板輻射導熱耦合問題

景 麗，王廣飛，唐紹鋒，梁 軍

(哈爾濱工業大學復合材料與結構研究所，哈爾濱 150001，jingli315@163.com)

金屬蜂窩夾芯板輻射導熱耦合問題

景 麗，王廣飛，唐紹鋒，梁 軍

(哈爾濱工業大學復合材料與結構研究所，哈爾濱 150001，jingli315@163.com)

針對金屬蜂窩夾芯板，研究了其在氣動加熱條件下的非穩態傳熱行為.基于高溫傳熱學原理，通過傳熱機制分析，建立了蜂窩夾芯板的導熱－輻射一維瞬態耦合傳熱數學物理模型.在此基礎上，采用控制容積法，結合蒙特卡羅法，發展形成了求解該類輻射導熱耦合傳熱問題的數值方法，并給出了蜂窩夾芯板當量熱導率和典型邊界條件下金屬蜂窩夾芯板瞬態溫度場、非加熱面熱響應的計算方法，進行了典型算例計算和實驗驗證.結果表明，所建立的數值計算模型在預報蜂窩結構熱響應方面是有效的，而且較之Swann-Pittman半經驗關系式，提高了計算精度.

蜂窩夾芯板;輻射導熱耦合;當量熱導率;數值預報

針對蜂窩夾芯板結構的傳熱研究大多是從無限大蜂窩夾芯板中選取周期分布的胞元，假定溫度沿板厚度方向線性分布［1］，或者二次函數分布［2］，或者某種多項式分布［3］，通過有限元計算結構的有效熱導率.近年來，在蜂窩夾芯結構熱分析理論和計算分析方面YI Long等［4］選取高階單元、采用高斯積分精確對單元表面變輻射熱流進行了研究.N.D.Kaushika等［5］推導了一種基于蜂窩芯層為灰體假設的物理模型，計算了蜂窩結構的輻射換熱量.Swann-Pittman半經驗關系式［6］已被做為一個標準模型用來計算蜂窩結構的熱傳導問題.

本文在進行蜂窩夾芯腔體內部傳熱機制分析的基礎上建立了蜂窩板導熱－輻射一維瞬態耦合傳熱數學物理模型.采用控制容積法，結合蒙特卡羅法［7］，發展形成了求解該類輻射導熱耦合傳熱問題的數值方法.給出了蜂窩夾芯板當量熱導率并揭示不同溫度下影響當量熱導率的主導因素，并給出了典型邊界條件下金屬蜂窩夾芯板瞬態溫度場、非加熱面熱響應的計算方法.

1 瞬態導熱－輻射耦合傳熱過程的理論模型

1.1 傳熱機制及當量熱導率的計算

蜂窩板由上、下蒙皮及中間蜂窩芯層構成，蜂窩板外表面在氣動熱流加熱下導致外表面溫度升高，外表面蒙皮將對外輻射能量，通過外表面蒙皮將大部分的熱量qR輻射回大氣環境，其余熱量qa以熱傳導的方式進入蜂窩板.當蜂窩尺寸在一定的范圍內時，通過蜂窩胞體腔內空氣對流傳遞的熱量可以忽略［8］，因此，蜂窩內部的熱量傳遞主要由傳導和輻射的耦合作用完成.如圖1所示，蜂窩板內存在3種熱量傳遞機制:1)蜂窩腔表面間的輻射換熱qr.2)固體壁面的熱傳導qsc.3)蜂窩腔內氣體的導熱qg.

圖1 蜂窩板內傳熱機理模型

計算蜂窩芯層當量熱導率，給出邊界條件為

式中:Tw1，Tw2分別為上、下邊界溫度，H為蜂窩夾芯板厚度，x=0面為非加熱面.

按上述邊界條件求解蜂窩傳熱控制方程，可以獲得穩態溫度場，進而計算得到蜂窩下底面上的導熱熱流密度大小為

式中:keff為蜂窩的有效導熱系數，根據蜂窩芯層固體和氣體兩相體積分數進行平均計算.文中普遍使用串聯模型，即

式中:ε為氣相所占有的體積分數，ks，kg分別為固相［9］和氣相［10－11］導熱系數.

下邊界面上的總熱流密度為導熱熱流密度與下邊界面上接收到的輻射熱流密度之和，即

式中:qrad為采用蒙特卡羅方法計算的下邊界面上的輻射熱流密度.

1.2 非穩態溫度場的計算

蜂窩結構具有周期性，可以將其中一個蜂窩單元作為研究對象，如圖2所示.以壁面的中間面作為絕熱邊界面.金屬蜂窩結構尺寸都不大，計算的Ra數一般很小，可以忽略自然對流換熱作用.沿蜂窩板厚度方向的一維非穩態傳熱控制方程為

式中:T為熱力學溫度，t為時間，ρ，c，k分別為蜂窩芯的有效密度、有效比熱與有效導熱系數.S為蜂窩腔體內表面間輻射換熱源項，采用蒙特卡羅方法計算.其中，蜂窩芯的有效密度根據兩相材料密度和體積分數進行平均為

式中:ρs為蜂窩芯金屬密度，ρg為蜂窩內氣體密度，As，Ag分別為蜂窩單元橫截面內固相和氣相所占有的面積.

圖2 正六邊形蜂窩單元

蜂窩芯的有效比熱根據兩相比熱和體積分數進行平均，即

式中:cs，cg分別為固相和氣相的比熱.

蜂窩內某一控制體單元ΔVi的輻射換熱源項為

式中:σ為黑體輻射常數，RD*ji為歸一化輻射傳遞系數，對控制體ΔVi上參與輻射換熱的面元ΔAj有RD*ji= εj·ΔAj·RDji，其中，RDji為輻射傳遞系數，表示單位時間內蜂窩腔內面積單元ΔAj發射的熱輻射能最終被單元ΔAi吸收的部分，Tj，Ti分別為面積單元ΔAj、ΔAi的溫度.

在3類邊界條件下求解蜂窩夾芯結構非穩態傳熱控制方程，得到溫度場，實現蜂窩夾芯結構傳熱分析.為方便起見，初始條件均為

式中:Tw1，Tw2分別為上、下邊界溫度，H為蜂窩夾芯板厚度，x=0面為非加熱面.

2)第2類邊界條件.

由于蜂窩面板通常在熱防護系統最外層，而且溫度高，是對環境的輻射散熱表面，因此，要考慮表面對環境的輻射散熱;此外，蜂窩夾芯板兩側的蒙皮具有一定的質量，雖然非常薄，如果忽略它的質量可能會引入較大的誤差，因此，將蜂窩夾芯結構考慮為一個具有一定質量，并且對環境輻射散熱的邊界，它的質量等于蜂窩夾芯板的蒙皮質量，則給定環境熱流密度的第2類邊界條件為

式中:qBD為外界環境對邊界的熱流密度，εBD為邊界對外界環境的發射率，qradint為邊界面向蜂窩腔內部的輻射熱流密度，qcondint為邊界面蜂窩腔內部的導熱熱流密度，ρ，cp，L分別為蒙皮的密度、比熱、厚度，τ為非穩態過程的時間.

3)第三類邊界條件.

由于蜂窩面板的非加熱面是對外界流體的對流換熱表面，基于第2類邊界條件，考慮蜂窩夾芯板蒙皮質量和對流換熱的第3類邊界條件可以寫為

式中:αBD為外界流體對邊界的對流換熱系數，Tf為外界流體溫度.

1.3 輻射導熱耦合換熱控制方程的數值求解

蜂窩板內輻射導熱耦合傳熱數值計算分為溫度場迭代求解與輻射換熱場的計算求解.

在溫度場迭代求解過程中，將能量方程中的輻射換熱處理為源項，使輻射－導熱耦合傳熱能量方程具有導熱微分方程的形式，從而將輻射－導熱耦合傳熱計算分解為輻射換熱源項計算、含源項的導熱微分方程求解兩個子過程，通過這兩個子過程迭代，實現耦合傳熱計算.采用控制容積法［12］對計算域進行均勻離散，邊界條件的處理采用補充邊界節點離散代數方程［13］方法，由于輻射源項的非線性，在開發熱分析計算模塊過程中采用了欠松弛迭代［14］，以保證數值計算的穩定.

在輻射換熱場的計算過程中，采用蒙特卡羅法將熱輻射的傳輸過程分解為發射、反射、吸收、散射等一系列獨立的子過程，并建立每個過程的概率模型.對每個單元(面元和體元)進行一定數目的能束抽樣，跟蹤、統計每束能束的歸宿(被介質和界面吸收或從系統中逸出)，從而獲得輻射傳遞因子.通過引入各離散單元之間的輻射傳遞系數，將輻射傳遞過程求解與溫度場求解分離，從而大大降低輻射換熱的求解復雜性，提高了輻射導熱耦合換熱數值計算的計算效率.

2 數值結果與驗證

2.1 當量熱導率的數值求解

蜂窩芯結構參數如表1所示，上、下面板與蜂窩芯材料均為 Haynes 214，氣體壓力 1.013×105Pa，計算中光子抽樣個數為105，溫度場迭代誤差為10－5K.

表1 蜂窩幾何參數

研究輻射換熱對蜂窩夾芯結構當量熱導率的影響，當蜂窩腔發射率分別取值為0.86、0.50和0. 10，上、下邊界溫度相差10 K時當量熱導率隨溫度變化的計算結果對比如圖3所示.對比顯示隨著溫度的升高，輻射換熱的影響會逐步增強，如不考慮它的影響，就會給計算結果帶來很大誤差.

圖3 不同發射率下蜂窩當量熱導率

2.2 傳熱分析及實驗驗證

為了進行實驗驗證分析，設計制造了高真空平板石墨加熱爐，可以實現對大尺寸(450 mm×450 mm×5 mm～80 mm)試件乃至熱防護系統整體樣件防、隔熱性能的測試，最高加熱溫度可達1 600℃、可控壓力范圍6.67×10－3～1.013 25×105Pa.爐體內結構示意圖如圖4所示.

圖4 加熱爐真空室示意圖

用于實驗測試的試樣為INCONEL 617蜂窩夾芯板，尺寸為200 mm×200 mm×5 mm，蒙皮厚度為0.12 mm，單元胞壁厚度為0.076 mm，蜂窩直徑5.6 mm，輸出控制熱流為40 kW/m2.為了避免邊界的影響，將試樣中心125 mm×12 5 mm的區域作為數據采集區，在試樣和水冷板之間鋪設了9只K型熱電偶及3只熱流傳感器，如圖5所示.在試樣的上表面放置了3只K型熱電偶，位置與圖5中熱流傳感器的位置相同.溫度和熱流數據的采集采用多路溫度測量儀，利用計算機窗口實時監控和記錄實驗數據，實驗在氮氣環境下完成.

圖5 實驗裝置簡圖

1)邊界條件.由于加熱表面直接與熱源接觸，邊界條件為施加的熱流密度即第2類邊界條件.而對于非加熱表面，為更接近實際條件，必須考慮蜂窩夾芯板蒙皮質量和與對流換熱的影響，因此采用了第3類邊界條件，根據水平放置平板自然對流換熱關聯式［15］，與下冷水板的對流換熱系數取4.6 W/(m2·K).

2)初始條件.室內空氣溫度取293 K，初始溫度為293 K.

蜂窩板外表面發射率取0. 86，內表面發射率按照被氧化了的鋁表面發射率取值為0.3［16］.蒙特卡羅模擬過程中每個表面單元發射能束數目取為105，能量方程迭代計算誤差取為10－5J.

蜂窩夾芯板熱分析計算模塊給出的蜂窩板上、下蒙皮表面溫度響應與實驗結果對比如圖6所示，對比顯示計算結果與實驗結果吻合很好.

圖6 計算結果與實驗結果比較

圖7給出了數值計算所得50 s時沿厚度x方向的溫度分布，其中，x=0 mm為非加熱面，x=5 mm為加熱面.從圖7可以看出，溫度分布呈現很規則的線性規律，其趨勢同文獻［17］中吻合.

圖7 在50 s時沿厚度方向的溫度分布

對于厚度為H的試樣，一維熱傳導達到穩態后，若熱邊、冷邊溫度及通過試樣的熱流分別為:Theat、Tcool、q，則試樣的當量熱導率為

圖8給出了蜂窩板當量熱導率計算模塊給出的不同溫度下采用相同熱邊界條件時，當量熱導率計算結果、實驗結果與Swann-Pittman半經驗公式計算結果的比較.

圖8 蜂窩夾芯板的當量熱導率的比較

結果表明，實驗結果預報結果在溫度較低時存在一定誤差，相對誤差最大為24%，本文認為引起誤差的來源主要包括:熱流值、溫度值以及厚度值的測量誤差.超過573 K后預報值與實驗值相對誤差小于6%，兩者吻合較好.較之根據數值結果擬合的Swann-Pittman半經驗公式，輻射導熱的耦合計算采用熱網絡法和蒙特卡羅法，在計算精度、應用范圍上都有提高.

3 結論

1)通過傳熱機制分析，建立了蜂窩夾芯板的導熱－輻射瞬態耦合傳熱數學物理模型及數值求解方法，經過典型算例計算和實驗驗證表明該方法是正確的，適合金屬熱防護系統蜂窩夾芯結構的傳熱分析.

2)根據計算結果分析，隨著溫度的升高，夾芯結構內部輻射換熱的影響會逐步增強，與實驗結果吻合較好;當溫度超過573 K后，本方法預報值比Swann-Pittman半經驗公式計算結果更接近實驗值，相對誤差小于6%.

3)本文建立的輻射導熱耦合計算方法在揭示蜂窩夾芯板傳熱機制、計算精度和應用范圍上都有提高.

［1］TENEK L，ARGYRIS J，OEBERG F.A multilayer composite triangular element for steady-state conduction/convection/radiation heat transfer in complex shells［J］.Computer Methods in Applied Mechanics and Engineering， 1994，120(3/4):271－301.

［2］王勖成，唐永進.一般殼體溫度場的有限元分析［J］.清華大學學報(自然科學版)， 1989，29(5):103－111.

［3］SURANA K S，ORTH N J.Completely hierarchical P－version anisymmetric shell element for nonlinear heat conduction in laminated composites［J］.Computers ＆Structures， 1992，46(5):777－789.

［4］YI Long，PENG Yun，SUN Qin.Research of the higher-order finite element arithme-tic for radiation exchange［J］.Chinese Journal of Aeronautics， 2006，9(3):197－202.

［5］KAUSHIKA N D，ARULANANTHAM M.Radiative heat transfer across transparent honey-comb insulation materials［J］.Heat Mass Transfer， 1995，22(5):751 －760.

［6］DARYABEIGI K.Heat transfer in adhesively bonded honeycomb core panels［J］.Journal of Thermophysics and Heat Transfer， 2002，16(2):217－221.

［7］HOWELL J R.Monte Carlo method in radiative heat transfer［J］.Journal of Heat Transfer， 1998，120(3)，547－560.

［8］KAUSHIKA N D，PADMAPRIYA R，SINGH T P.Convection theory of honeycomb and slat-devices for solar energy applications［J］.Energy Res Technol， 1992，35(2):126－146.

［9］KAMRAN D.Effective Thermal Conductivity of High Temperature Insulations for Reusable Launch Vehicles［R］.Washington，DC:NASA Technical memorandum，1999:208972.

［10］KAMRAN D.Thermal Analysis and Design of Multilayer Insulation for Re-entry Aerodynamic Heating［R］.Easton，PA:AIAA，2001:2001－2834.

［11］POTEET C，KHAJEEL H A，HSU S Y.Preliminary Thermal-Mechanical Sizing of Metallic TPS［R］.Easton，PA:AIAA，2002:2002－0505.

［12］魯懷敏.用控制容積法求解連鑄方坯溫度場［J］.沙洲職業工學院學報， 2004，7(1):12－21.

［13］楊世銘.傳熱學［M］.北京:高等教育出版社，2006:167－169.

［14］葛德彪.平面介質波導折射率剖面重建的欠松弛迭代方法［J］.光子學， 1996，25(5):439－445.

［15］楊世銘.傳熱學［M］.北京;高等教育出版社，2006:212.

［16］楊世銘.傳熱學［M］.北京:高等教育出版社，2006:369.

［17］陳勇，高德平.金屬蜂窩平板加熱過程的數值模擬及試驗研究［J］.理化檢驗(2)物理分冊， 2003，39(5):234－236.

Radiation and conduction coupling problems of honeycomb sandwich panel

JING Li，WANG Guang-fei，TANG Shao-feng，LIANG Jun
(Centre for Composite Materials and Structures，Harbin Institute of Technology，Harbin 150001，China，jingli315@163.com)

In this paper，the non-steady-state heats transfer behavior of metal honeycomb sandwich panel under aerodynamic heating was studied.Based on high temperature heat transfer theory and mechanism analysis，one-dimensional transient heat transfer mathematical physics model of honeycomb sandwich panel taking the combination of heat conduction and radiation into account was given.The control volume method and Monte Carlo method were adopted and then a numerical computation method was developed to calculate the coupling problem of heat conduction and radiation.Computation methods for effective thermal conductivities of honeycomb sandwich panel，the transient temperature field and heat response of non-heated surface under typical boundary conditions were given as well.An experiment was designed to validate this computation method.The result shows that the numerical computation model given in this paper is valid to predict the thermal performance of honeycomb sandwich structures and the values calculated are more precise than those calculated by Swann-Pittman semi-empirical relationship.

honeycomb sandwich panel;combination of heat conduction and radiation;effective thermal conductivities;numerical prediction

V435

A

0367－6234(2010)05－0827－05

2009－06－22.

國家自然科學基金資助項目( 10772060，90916027);

黑龍江省杰出青年基金資助項目(JC2006-13).

景 麗(1984—)，女，碩士研究生;

梁 軍(1969—)，男，教授，博士生導師.

(編輯 張 紅)