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

(哈爾濱工業(yè)大學(xué)復(fù)合材料與結(jié)構(gòu)研究所，哈爾濱 150001，jingli315@163.com)

金屬蜂窩夾芯板輻射導(dǎo)熱耦合問題

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

(哈爾濱工業(yè)大學(xué)復(fù)合材料與結(jié)構(gòu)研究所，哈爾濱 150001，jingli315@163.com)

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

蜂窩夾芯板;輻射導(dǎo)熱耦合;當(dāng)量熱導(dǎo)率;數(shù)值預(yù)報

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

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

1 瞬態(tài)導(dǎo)熱－輻射耦合傳熱過程的理論模型

1.1 傳熱機制及當(dāng)量熱導(dǎo)率的計算

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

圖1 蜂窩板內(nèi)傳熱機理模型

計算蜂窩芯層當(dāng)量熱導(dǎo)率，給出邊界條件為

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

按上述邊界條件求解蜂窩傳熱控制方程，可以獲得穩(wěn)態(tài)溫度場，進而計算得到蜂窩下底面上的導(dǎo)熱熱流密度大小為

式中:keff為蜂窩的有效導(dǎo)熱系數(shù)，根據(jù)蜂窩芯層固體和氣體兩相體積分?jǐn)?shù)進行平均計算.文中普遍使用串聯(lián)模型，即

式中:ε為氣相所占有的體積分?jǐn)?shù)，ks，kg分別為固相［9］和氣相［10－11］導(dǎo)熱系數(shù).

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

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

1.2 非穩(wěn)態(tài)溫度場的計算

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

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

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

圖2 正六邊形蜂窩單元

蜂窩芯的有效比熱根據(jù)兩相比熱和體積分?jǐn)?shù)進行平均，即

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

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

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

在3類邊界條件下求解蜂窩夾芯結(jié)構(gòu)非穩(wěn)態(tài)傳熱控制方程，得到溫度場，實現(xiàn)蜂窩夾芯結(jié)構(gòu)傳熱分析.為方便起見，初始條件均為

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

2)第2類邊界條件.

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

式中:qBD為外界環(huán)境對邊界的熱流密度，εBD為邊界對外界環(huán)境的發(fā)射率，qradint為邊界面向蜂窩腔內(nèi)部的輻射熱流密度，qcondint為邊界面蜂窩腔內(nèi)部的導(dǎo)熱熱流密度，ρ，cp，L分別為蒙皮的密度、比熱、厚度，τ為非穩(wěn)態(tài)過程的時間.

3)第三類邊界條件.

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

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

1.3 輻射導(dǎo)熱耦合換熱控制方程的數(shù)值求解

蜂窩板內(nèi)輻射導(dǎo)熱耦合傳熱數(shù)值計算分為溫度場迭代求解與輻射換熱場的計算求解.

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

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

2 數(shù)值結(jié)果與驗證

2.1 當(dāng)量熱導(dǎo)率的數(shù)值求解

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

表1 蜂窩幾何參數(shù)

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

圖3 不同發(fā)射率下蜂窩當(dāng)量熱導(dǎo)率

2.2 傳熱分析及實驗驗證

為了進行實驗驗證分析，設(shè)計制造了高真空平板石墨加熱爐，可以實現(xiàn)對大尺寸(450 mm×450 mm×5 mm～80 mm)試件乃至熱防護系統(tǒng)整體樣件防、隔熱性能的測試，最高加熱溫度可達1 600℃、可控壓力范圍6.67×10－3～1.013 25×105Pa.爐體內(nèi)結(jié)構(gòu)示意圖如圖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的區(qū)域作為數(shù)據(jù)采集區(qū)，在試樣和水冷板之間鋪設(shè)了9只K型熱電偶及3只熱流傳感器，如圖5所示.在試樣的上表面放置了3只K型熱電偶，位置與圖5中熱流傳感器的位置相同.溫度和熱流數(shù)據(jù)的采集采用多路溫度測量儀，利用計算機窗口實時監(jiān)控和記錄實驗數(shù)據(jù)，實驗在氮氣環(huán)境下完成.

圖5 實驗裝置簡圖

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

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

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

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

圖6 計算結(jié)果與實驗結(jié)果比較

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

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

對于厚度為H的試樣，一維熱傳導(dǎo)達到穩(wěn)態(tài)后，若熱邊、冷邊溫度及通過試樣的熱流分別為:Theat、Tcool、q，則試樣的當(dāng)量熱導(dǎo)率為

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

圖8 蜂窩夾芯板的當(dāng)量熱導(dǎo)率的比較

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

3 結(jié)論

1)通過傳熱機制分析，建立了蜂窩夾芯板的導(dǎo)熱－輻射瞬態(tài)耦合傳熱數(shù)學(xué)物理模型及數(shù)值求解方法，經(jīng)過典型算例計算和實驗驗證表明該方法是正確的，適合金屬熱防護系統(tǒng)蜂窩夾芯結(jié)構(gòu)的傳熱分析.

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

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

［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］.清華大學(xué)學(xué)報(自然科學(xué)版)， 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］.沙洲職業(yè)工學(xué)院學(xué)報， 2004，7(1):12－21.

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

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

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

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

［17］陳勇，高德平.金屬蜂窩平板加熱過程的數(shù)值模擬及試驗研究［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.

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

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

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

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

(編輯 張 紅)