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

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

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

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

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

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

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

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

本文在進(jìn)行蜂窩夾芯腔體內(nèi)部傳熱機(jī)制分析的基礎(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)的計(jì)算方法.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 正六邊形蜂窩單元

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

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

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

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

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

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

2)第2類邊界條件.

由于蜂窩面板通常在熱防護(hù)系統(tǒng)最外層，而且溫度高，是對環(huán)境的輻射散熱表面，因此，要考慮表面對環(huán)境的輻射散熱;此外，蜂窩夾芯板兩側(cè)的蒙皮具有一定的質(zhì)量，雖然非常薄，如果忽略它的質(zhì)量可能會(huì)引入較大的誤差，因此，將蜂窩夾芯結(jié)構(gòu)考慮為一個(gè)具有一定質(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)過程的時(shí)間.

3)第三類邊界條件.

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

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

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

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

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

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

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

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

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

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

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

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

2.2 傳熱分析及實(shí)驗(yàn)驗(yàn)證

為了進(jìn)行實(shí)驗(yàn)驗(yàn)證分析，設(shè)計(jì)制造了高真空平板石墨加熱爐，可以實(shí)現(xiàn)對大尺寸(450 mm×450 mm×5 mm～80 mm)試件乃至熱防護(hù)系統(tǒng)整體樣件防、隔熱性能的測試，最高加熱溫度可達(dá)1 600℃、可控壓力范圍6.67×10－3～1.013 25×105Pa.爐體內(nèi)結(jié)構(gòu)示意圖如圖4所示.

圖4 加熱爐真空室示意圖

用于實(shí)驗(yàn)測試的試樣為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ù)的采集采用多路溫度測量儀，利用計(jì)算機(jī)窗口實(shí)時(shí)監(jiān)控和記錄實(shí)驗(yàn)數(shù)據(jù)，實(shí)驗(yàn)在氮?dú)猸h(huán)境下完成.

圖5 實(shí)驗(yàn)裝置簡圖

1)邊界條件.由于加熱表面直接與熱源接觸，邊界條件為施加的熱流密度即第2類邊界條件.而對于非加熱表面，為更接近實(shí)際條件，必須考慮蜂窩夾芯板蒙皮質(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］.蒙特卡羅模擬過程中每個(gè)表面單元發(fā)射能束數(shù)目取為105，能量方程迭代計(jì)算誤差取為10－5J.

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

圖6 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果比較

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

圖7 在50 s時(shí)沿厚度方向的溫度分布

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

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

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

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

3 結(jié)論

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

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

3)本文建立的輻射導(dǎo)熱耦合計(jì)算方法在揭示蜂窩夾芯板傳熱機(jī)制、計(jì)算精度和應(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ìn).一般殼體溫度場的有限元分析［J］.清華大學(xué)學(xué)報(bào)(自然科學(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é)報(bào)， 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ù)值模擬及試驗(yàn)研究［J］.理化檢驗(yàn)(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é)基金資助項(xiàng)目( 10772060，90916027);

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

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

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

(編輯 張 紅)