高超聲速飛行器熱環(huán)境與結(jié)構(gòu)傳熱的多場(chǎng)耦合數(shù)值研究

周印佳，孟松鶴＊，解維華，楊強(qiáng)

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

高超聲速飛行器熱環(huán)境與結(jié)構(gòu)傳熱的多場(chǎng)耦合數(shù)值研究

周印佳，孟松鶴＊，解維華，楊強(qiáng)

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

為了準(zhǔn)確預(yù)測(cè)高超聲速飛行器面臨的嚴(yán)峻氣動(dòng)熱／力環(huán)境以及結(jié)構(gòu)的熱力響應(yīng)，發(fā)展了高超聲速流動(dòng)與結(jié)構(gòu)傳熱耦合框架。采用分區(qū)求解方法，通過耦合界面的實(shí)時(shí)數(shù)據(jù)傳遞，實(shí)現(xiàn)了基于Navier－Stokes方程的高超聲速化學(xué)非平衡計(jì)算流體力學(xué)（CFD）求解器與結(jié)構(gòu)的熱力全耦合有限元法（FEM）求解器的多場(chǎng)耦合計(jì)算，建立了高超聲速飛行器的多場(chǎng)耦合數(shù)值分析方法。首先對(duì)經(jīng)典高超聲速圓柱繞流實(shí)驗(yàn)進(jìn)行了耦合計(jì)算，結(jié)果與實(shí)驗(yàn)值吻合良好。然后針對(duì)典型的超高溫陶瓷（UHTC）材料的耦合傳熱問題進(jìn)行了數(shù)值研究，考慮熱傳導(dǎo)效應(yīng)對(duì)氣動(dòng)熱環(huán)境和結(jié)構(gòu)熱響應(yīng)預(yù)測(cè)的影響，結(jié)果表明對(duì)于復(fù)雜外形且熱導(dǎo)率相對(duì)較高的UHTC材料，結(jié)構(gòu)內(nèi)部熱傳導(dǎo)對(duì)熱環(huán)境和表面溫度分布的影響不可忽略。最后針對(duì)UHTC材料熱物性（比熱和熱導(dǎo)率）非線性對(duì)高超聲速流動(dòng)傳熱過程的影響進(jìn)行了研究，結(jié)果表明當(dāng)比熱和熱導(dǎo)率處于合理的誤差范圍內(nèi)時(shí)，材料表面溫度響應(yīng)對(duì)其變化并不敏感。

高超聲速飛行器；多場(chǎng)耦合；氣動(dòng)熱；數(shù)值模擬；熱防護(hù)

高超聲速飛行器的快速發(fā)展給熱防護(hù)設(shè)計(jì)帶來了更為嚴(yán)峻的挑戰(zhàn)［1］。準(zhǔn)確地預(yù)測(cè)氣動(dòng)熱／力環(huán)境、結(jié)構(gòu)溫度和應(yīng)力狀態(tài)，能夠在提高飛行器安全性能的同時(shí)減小熱防護(hù)系統(tǒng)設(shè)計(jì)冗余，對(duì)提高飛行器的性能有著極為重要的意義。

傳統(tǒng)的高超聲速飛行器的熱／力載荷環(huán)境的預(yù)測(cè)與熱防護(hù)結(jié)構(gòu)性能的分析基本還處于分離狀態(tài)。這種方法實(shí)際上把多場(chǎng)耦合的事實(shí)人為地分割成多個(gè)獨(dú)立的物理場(chǎng)。在這種情況下，既無法得到精確的氣動(dòng)熱／力載荷環(huán)境，也無法正確評(píng)價(jià)熱防護(hù)材料及其結(jié)構(gòu)的服役特征。

高超聲速飛行熱防護(hù)設(shè)計(jì)是一個(gè)涉及到真實(shí)氣體效應(yīng)、耦合傳熱和結(jié)構(gòu)熱力響應(yīng)的復(fù)雜的多場(chǎng)耦合問題，必須采用多場(chǎng)耦合的方法求解。國外較早開展了多場(chǎng)耦合方法的研究，并在1987年首次發(fā)表了關(guān)于多場(chǎng)耦合分析研究方面的成果。Wieting等［2－3］做了高超聲速圓柱繞流試驗(yàn)，此后經(jīng)常被用于高超聲速流動(dòng)－傳熱耦合數(shù)值模擬的驗(yàn)證。Thornton等［4－5］采用流－熱－固一體化有限元法（FEM）研究了飛行器面板以及前緣的氣動(dòng)熱－結(jié)構(gòu)相互作用，得到了與實(shí)驗(yàn)相符的冷壁熱流。Loehner等［6］通過在每個(gè)時(shí)間步對(duì)耦合界面處的變量進(jìn)行插值，成功地將已有的計(jì)算流體力學(xué)（CFD）、計(jì)算熱力學(xué)和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)程序耦合在一起。近期的研究有，Miller等［7］發(fā)展了多物理場(chǎng)時(shí)間推進(jìn)的流動(dòng)－傳熱－結(jié)構(gòu)分區(qū)松耦合方法。該方法充分利用了流－熱－結(jié)構(gòu)的時(shí)間尺度差異，在每個(gè)傳熱分析時(shí)間步內(nèi)有多個(gè)結(jié)構(gòu)分析時(shí)間步，每個(gè)結(jié)構(gòu)時(shí)間步內(nèi)有多個(gè)流動(dòng)時(shí)間步，兼顧了計(jì)算精度和效率。

21世紀(jì)初，國內(nèi)學(xué)者首先對(duì)流場(chǎng)、熱、結(jié)構(gòu)一體化數(shù)值模擬方法展開研究。黃唐［8］和夏剛［9］等開展了二維流動(dòng)－傳熱耦合模擬，二者的研究中固體傳熱均采用FEM，流體則分別采用有限差分法和有限體積法。桂業(yè)偉［10］、耿湘人［11］、Zhao［12］和張兵［13］等相繼開展了流場(chǎng)、熱、結(jié)構(gòu)多場(chǎng)耦合計(jì)算研究和分析。國內(nèi)學(xué)者最近的研究有，Zhang等［14］開發(fā)了預(yù)測(cè)高超聲速氣動(dòng)熱環(huán)境的流動(dòng)－傳熱－結(jié)構(gòu)一體化分析方法，采用松耦合方法耦合了穩(wěn)態(tài)的CFD求解器與瞬態(tài)熱結(jié)構(gòu)力學(xué)求解器。董維中等［15］基于流場(chǎng)的非平衡Navier－Stokes方程、表面的能量守恒方程和內(nèi)部的熱傳導(dǎo)方程建立了表面溫度分布與氣動(dòng)熱的耦合計(jì)算方法，并考慮了流場(chǎng)的非平衡效應(yīng)、表面的熱輻射效應(yīng)、催化效應(yīng)和燒蝕效應(yīng)以及熱防護(hù)層內(nèi)部的熱傳導(dǎo)效應(yīng)。

由于多場(chǎng)耦合問題的復(fù)雜性，還需要進(jìn)一步開展分析方法研究，深刻把握防熱系統(tǒng)多場(chǎng)耦合規(guī)律及其效應(yīng)。需要指出的是，高超聲速飛行器面臨嚴(yán)峻的氣動(dòng)熱環(huán)境，材料在高溫下的熱物性（熱導(dǎo)率和比熱）非線性對(duì)耦合傳熱過程的影響目前尚未見到相關(guān)研究。

本文發(fā)展了高超聲速流動(dòng)和結(jié)構(gòu)傳熱的耦合框架，實(shí)現(xiàn)了高超聲速非平衡流動(dòng)求解器與結(jié)構(gòu)熱力全耦合的多場(chǎng)耦合計(jì)算，建立了高超聲速飛行器的多場(chǎng)耦合數(shù)值分析方法?？紤]熱傳導(dǎo)效應(yīng)對(duì)氣動(dòng)熱環(huán)境和結(jié)構(gòu)熱響應(yīng)預(yù)測(cè)的影響，并針對(duì)材料的比熱和熱導(dǎo)率非線性對(duì)高超聲速流動(dòng)傳熱過程的影響進(jìn)行了研究。

1 耦合分析模型與耦合策略

本文的耦合分析模型如圖1所示。在該模型中，流場(chǎng)采用同時(shí)求解連續(xù)、動(dòng)量、能量等守恒方程的耦合解格式得到熱流qw和壓力pw，結(jié)構(gòu)場(chǎng)采用熱力全耦合的方法進(jìn)行求解得到壁面溫度Tw和結(jié)構(gòu)位移us。這兩個(gè)獨(dú)立的模型通過在流－固耦合界面上實(shí)時(shí)交換參數(shù)來表征外部高超聲速流動(dòng)與內(nèi)部結(jié)構(gòu)響應(yīng)和熱響應(yīng)的耦合作用。熱流和壓力的聯(lián)合載荷影響結(jié)構(gòu)的熱力行為，相應(yīng)的壁面溫度和位移影響外部流場(chǎng)的氣動(dòng)熱力學(xué)行為。

圖1 多場(chǎng)耦合分析模型Fig．1 Multi－field coupling analysis model

采用分區(qū)求解方法完成對(duì)高超聲速流動(dòng)與結(jié)構(gòu)傳熱的耦合分析。流體和固體區(qū)域的求解器均為瞬態(tài)求解，每個(gè)求解器所需要的下一步計(jì)算的數(shù)據(jù)在耦合界面上反復(fù)交換。圖2中給出了耦合策略示意圖，詳細(xì)描述如下：

1）在等溫壁面邊界條件下（結(jié)構(gòu)溫度確定的定常壁溫）進(jìn)行穩(wěn)態(tài)流動(dòng)計(jì)算，結(jié)果作為瞬態(tài)耦合分析的初始條件。在耦合起始點(diǎn)，進(jìn)行瞬態(tài)流動(dòng)計(jì)算，得到初始時(shí)刻熱流qw和壁面壓力pw。

2）qw和pw在流體－固體耦合界面進(jìn)行插值，作為結(jié)構(gòu)求解器的邊界條件，在熱流和壓力聯(lián)合載荷下進(jìn)行熱－結(jié)構(gòu)耦合分析。一個(gè)時(shí)間步Δt后，結(jié)構(gòu)求解器得出壁面溫度Tw和結(jié)構(gòu)位移us。

3）us通過界面插值進(jìn)行傳遞，更新流場(chǎng)網(wǎng)格，Tw更新流體求解的壁面溫度邊界條件，進(jìn)行下一步流場(chǎng)計(jì)算。經(jīng)過一個(gè)時(shí)間步Δt后得到新的qw和pw，提取出來用于下一時(shí)間步的結(jié)構(gòu)計(jì)算。

該過程不斷重復(fù)直到計(jì)算時(shí)間結(jié)束。

圖2 本文所用耦合策略Fig．2 Coupling strategy in this paper

2 非平衡流場(chǎng)計(jì)算

2．1 流動(dòng)方程

文中考慮高超聲速飛行器周圍的流場(chǎng)為化學(xué)非平衡、熱力學(xué)平衡的黏性可壓縮連續(xù)流動(dòng)，連續(xù)、動(dòng)量、能量守恒方程分別為對(duì)于包含組分混合或反應(yīng)的流動(dòng)，還需要組

分守恒方程，即

式（1）～式（4）中：p和ρ分別為混合氣體的壓力和密度；ν為速度向量；τ為應(yīng)力張量；E為比總能量；λ為氣體熱導(dǎo)率；h為焓；Ji為組分i的質(zhì)量擴(kuò)散通量；Yi為組分i的質(zhì)量分?jǐn)?shù)；Ri為組分i的凈生產(chǎn)率。

組分的分壓由組分密度和混合氣體的溫度求得，混合氣體的壓力由Dalton定律［16］給出。

2．2 化學(xué)動(dòng)力學(xué)模型

由于存在動(dòng)能耗散和激波，高超聲速飛行器周圍的空氣會(huì)達(dá)到極高的溫度，高溫使氣體離解甚至電離?；瘜W(xué)非平衡假設(shè)即特征化學(xué)反應(yīng)時(shí)間與流動(dòng)的特征時(shí)間相當(dāng)。采用Park和Lee［17］的5組分（O，N，NO，O2，N2）17反應(yīng)化學(xué)動(dòng)力學(xué)模型，并考慮第三體效應(yīng)。如表1所示，反應(yīng)機(jī)制包括3個(gè)分解反應(yīng)和2個(gè)交換反應(yīng)。Park動(dòng)力學(xué)模型中的常數(shù)Ar、βr和Er在表1中給出。

表1 Park 5組分動(dòng)力學(xué)模型及反應(yīng)速率參數(shù)Table 1 Reaction－rate parameters based on Park’s 5－species kinetic model

2．3 輸運(yùn)屬性

單一組分的黏性系數(shù)μi、熱導(dǎo)率λi和擴(kuò)散系數(shù)Dij由氣體分子動(dòng)力論［18］計(jì)算得到。

式中：R為普適氣體常數(shù)；Cpi為組分i的比熱；Mi為組分i的分子量；σij＝（σi＋σj）／2，σi為組分i的特征分子直徑；Ω為碰撞積分。

混合氣體的系數(shù)可以由半經(jīng)驗(yàn)Wilke公式［19］給出，即

式中：Xi為組分i的摩爾分?jǐn)?shù)。

混合氣體的多組分?jǐn)U散系數(shù)由式（11）的近似表達(dá)式計(jì)算。

2．4 熱力學(xué)模型

有限速率化學(xué)反應(yīng)流動(dòng)的計(jì)算需要每種單一組分的熱力學(xué)屬性（比熱Cpi和焓hi）。本文中只考慮流動(dòng)為化學(xué)非平衡、熱力學(xué)平衡狀態(tài)的反應(yīng)氣體，每種組分的熱力學(xué)屬性由當(dāng)?shù)仂o溫計(jì)算。高溫氣體的Cpi和hi可以通過溫度的多項(xiàng)式曲線擬合函數(shù)［20］計(jì)算：

式中：Ai（i＝1～6）為擬合系數(shù)。

混合氣體的屬性可由單一組分的屬性確定：

式中：NS為組分?jǐn)?shù)量。

2．5 數(shù)值方法

對(duì)于高超聲速化學(xué)非平衡流動(dòng)計(jì)算，連續(xù)、動(dòng)量、能量和組分方程采用有限體積法同時(shí)求解。已有研究［21］表明，不同數(shù)值格式對(duì)熱流計(jì)算精度有較大影響。通過大量數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，表明AUSM＋格式在熱流計(jì)算精確性方面具有優(yōu)勢(shì)。因此，本文計(jì)算采用AUSM＋格式。

3 熱－結(jié)構(gòu)控制方程

3．1 熱傳導(dǎo)

基于能量守恒定律和Fourier定律，可以得到結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)方程

式中：λs為固體熱導(dǎo)率；Cps為固體比熱容。

3．2 結(jié)構(gòu)響應(yīng)

當(dāng)固體結(jié)構(gòu)受到加熱時(shí)，它的體內(nèi)溫度會(huì)發(fā)生變化，這時(shí)在結(jié)構(gòu)內(nèi)部會(huì)有溫差存在。由于熱膨脹以及結(jié)構(gòu)的約束作用，固體結(jié)構(gòu)產(chǎn)生熱應(yīng)力，熱應(yīng)力又導(dǎo)致了熱變形的發(fā)生。對(duì)于二維結(jié)構(gòu)的響應(yīng)方程為

式中：σx、σy和τxy為固體結(jié)構(gòu)的應(yīng)力分量；而D為彈性矩陣；B和S分別為應(yīng)變和應(yīng)力矩陣；δe為位移矩陣。

在本文計(jì)算中，暫不考慮模量和泊松比等彈性參數(shù)隨溫度的變化。

固體區(qū)域內(nèi)，采用有限元法進(jìn)行熱－結(jié)構(gòu)分析。在本文研究中，應(yīng)力分析與溫度分布相關(guān)，而溫度分布對(duì)應(yīng)力為弱相關(guān)。但是在當(dāng)前工作中，仍然考慮應(yīng)力與溫度分布之間的雙向耦合關(guān)系，進(jìn)行了熱力全耦合分析。溫度用向后差分格式積分，非線性耦合系統(tǒng)采用牛頓法求解。牛頓迭代法采用非對(duì)稱雅克比矩陣，描述為如式（18）所示的矩陣方程組。

式中：ΔU和ΔT分別為修正后的位移和溫度增量；Kij（i，j＝U，T）為完全耦合雅克比矩陣的子陣；RU和RT分別為計(jì)算結(jié)構(gòu)力與熱的殘差。采用這種方式求解方程組，熱力與結(jié)構(gòu)問題是同時(shí)求解的。

4 耦合框架驗(yàn)證

為了驗(yàn)證耦合策略和數(shù)值方法，選擇了經(jīng)典的圓柱繞流實(shí)驗(yàn)作為算例。實(shí)驗(yàn)來流馬赫數(shù)、壓力和溫度分別為6．47、648Pa和241．5K，不銹鋼管長(zhǎng)度、直徑和厚度分別為0．609 6m、0．076 2m和0．012 7m。計(jì)算中采用標(biāo)準(zhǔn)鋼材的熱力學(xué)參數(shù)，其中：熱導(dǎo)率為16．27W／（m·K），比熱為502．48J／（kg·K），密度為8 030kg／m3，彈性模量為1．2×1013Pa，熱膨脹系數(shù)為1．68×10－51／K，泊松比為0．3。實(shí)驗(yàn)的更多細(xì)節(jié)由文獻(xiàn)［3］給出。

流體區(qū)域?yàn)榻Y(jié)構(gòu)網(wǎng)格，固體區(qū)域采用四邊形單元，二者在耦合界面處不匹配。邊界層網(wǎng)格進(jìn)行加密以保證計(jì)算結(jié)果的網(wǎng)格無關(guān)性，并使其具有足夠的分辨率。流場(chǎng)和圓柱的計(jì)算模型如圖3所示。計(jì)算中采用了SST k－ω 湍流模型。

圖4給出了不同耦合時(shí)間步長(zhǎng)時(shí)駐點(diǎn)處的溫度隨時(shí)間變化曲線，圖5為駐點(diǎn)熱流變化趨勢(shì)。可以看出，隨著加熱時(shí)間增長(zhǎng)，駐點(diǎn)溫度逐漸升高，而熱流變化趨勢(shì)則相反。另外，溫度和熱流都在初始的一小段時(shí)間內(nèi)變化劇烈，而后變化趨于緩和?；诓煌詈蠒r(shí)間步長(zhǎng)得到了不同的溫度分布結(jié)果，其中Δt＝0．000 1s時(shí)得到t＝2s時(shí)的

圖3 圓柱繞流計(jì)算模型Fig．3 Calculation model for flow over a cylinder

圖4 駐點(diǎn)溫度隨時(shí)間變化曲線Fig．4 Time histories of stagnation point temperature

圖5 駐點(diǎn)熱流隨時(shí)間變化曲線Fig．5 Time histories of stagnation point heat flux

結(jié)果（384．9K）與實(shí)驗(yàn)結(jié)果（約為388．72K）最為接近，而對(duì)于Δt＝0．001s和Δt＝0．01s，得到了略偏高的溫度結(jié)果?？梢缘贸鼋Y(jié)論，隨著時(shí)間步長(zhǎng)偏離某個(gè)合適的特定值逐漸增大，預(yù)測(cè)的結(jié)構(gòu)溫度相應(yīng)會(huì)逐漸偏高。以下取Δt＝0．000 1s的結(jié)果進(jìn)行分析。

圖6和圖7給出了不同時(shí)間點(diǎn)的溫度和熱流分布。在氣動(dòng)加熱開始后的前0．5s內(nèi)，溫度和熱流變化比較明顯。

圖8顯示，在熱力聯(lián)合載荷下，駐點(diǎn)處的位移和Mises應(yīng)力都隨時(shí)間逐漸增大，但是在當(dāng)前工況下位移量較小，不足以明顯影響流場(chǎng)。圖9給出了流體域中心線上的壓力和溫度分布，所預(yù)測(cè)的激波后溫度和壓力與實(shí)驗(yàn)值非常吻合。

圖8 結(jié)構(gòu)Mises應(yīng)力和位移Fig．8 Mises stress and displacement

圖9 駐點(diǎn)線溫度和壓力分布Fig．9 Fluid temperature and pressure distributions along stagnation line

5 材料熱物性對(duì)耦合傳熱的影響

本算例主要分析超高溫陶瓷（UHTC）雙錐的耦合傳熱問題，計(jì)算對(duì)象為德國L2K高溫風(fēng)洞中的軸對(duì)稱UHTC雙錐。雙錐幾何尺寸如圖10

圖10 雙錐幾何示意圖Fig．10 Schematic of double cone geometry

所示，后部有銅質(zhì)夾具支撐。來流馬赫數(shù)、壓力和溫度分別為4．57，272Pa和764K，混合氣體中N2、O2和O組分的質(zhì)量分?jǐn)?shù)分別為0．739、0．113和0．092。材料表面為完全催化壁面，其輻射率為0．85。更多試驗(yàn)細(xì)節(jié)和計(jì)算結(jié)果見文獻(xiàn)［22］和文獻(xiàn)［23］。

為了能夠精確地預(yù)測(cè)結(jié)構(gòu)的熱流和溫度分布，采用了溫度相關(guān)的材料熱物性［24］，如圖11所示。本次計(jì)算采用ZrB2－30vol．%SiC的材料物性，其余參數(shù)用于后面對(duì)比分析熱物性對(duì)傳熱影響的研究。計(jì)算過程中的流動(dòng)為層流狀態(tài)。

由前面的算例可知，耦合時(shí)間步長(zhǎng)Δt＝0．001s可以得到比較準(zhǔn)確的結(jié)果并具有較高的效率，因此計(jì)算采用耦合時(shí)間步長(zhǎng)Δt＝0．001s，計(jì)算時(shí)間總長(zhǎng)為90s。壁面設(shè)為完全催化壁面，即所有原子均在壁面發(fā)生復(fù)合，則壁面處的組分濃度與來流組分相同，即（Yi）w＝Y(jié)i，∞。

圖12給出的是本文方法計(jì)算的距離前緣4mm和35mm處的壁面溫度，與文獻(xiàn)［12］中的數(shù)值結(jié)果及實(shí)驗(yàn)測(cè)得結(jié)果的對(duì)比。在x＝4mm處，本文計(jì)算所得壁面溫度結(jié)果低于實(shí)驗(yàn)值。在t＜40s時(shí)計(jì)算結(jié)果與實(shí)驗(yàn)值比較接近，但是在t＞40s時(shí)，與實(shí)驗(yàn)值偏差逐漸加大，直到t＝90s時(shí)最大相差約91K。而在x＝35mm處，本文結(jié)果與實(shí)驗(yàn)值和文獻(xiàn)值［22］都吻合得較好。造成頭部最終溫度偏低的原因可能與采用熱力平衡模型有關(guān)。

圖11 超高溫陶瓷（UHTC）材料隨溫度變化的熱物性Fig．11 Temperature dependent thermophysical properties of ultra high temperature ceramics（UHTC）

圖12 壁面溫度隨時(shí)間變化曲線Fig．12 Time histories of wall temperature

傳統(tǒng)的熱防護(hù)系統(tǒng)一般是良好的隔熱體，熱導(dǎo)率較小，并能夠有效地輻射熱量。因此計(jì)算熱載荷與表面溫度時(shí)通常忽略固體內(nèi)部熱傳導(dǎo)，對(duì)流傳熱與向外輻射熱之間形成局部熱平衡。此時(shí)計(jì)算得到的局部輻射平衡壁面溫度能夠滿足工程需求。而對(duì)于如UHTC這類熱導(dǎo)率較高的新型防熱材料，表面對(duì)流加熱一部分向固體內(nèi)部傳導(dǎo)，另一部分輻射回大氣。當(dāng)系統(tǒng)達(dá)到穩(wěn)態(tài)時(shí)，形成全局輻射熱平衡。因此，對(duì)于超高溫陶瓷類前緣結(jié)構(gòu)，局部熱輻射平衡被打破，需要考慮結(jié)構(gòu)內(nèi)部熱傳導(dǎo)。

圖13給出了是否考慮熱傳導(dǎo)情況下表面溫度分布和熱流對(duì)比?？梢钥吹?，考慮熱傳導(dǎo)以后，雙錐的頭部和肩部的大量熱量傳導(dǎo)到溫度較低的中間部分和尾部，從而使頭部和肩部的溫度顯著降低，其中駐點(diǎn)處溫度比輻射平衡狀態(tài)的溫度降低了500K左右。由于頭部和肩部溫度變化幅度較大，因此導(dǎo)致其熱流增大明顯；而熱流相對(duì)較小的中部和尾部，表面熱流由于溫度升高而略有降低。因此，對(duì)于外部環(huán)境變化劇烈或者熱導(dǎo)率較高且外形復(fù)雜的結(jié)構(gòu)，結(jié)構(gòu)內(nèi)部熱傳導(dǎo)對(duì)熱環(huán)境和表面溫度分布的影響不可忽略。

圖13 有無熱傳導(dǎo)情況下雙錐表面溫度和熱流分布對(duì)比Fig．13 Comparison of surface temperature and surface heat flux distribution with and without heat conduction

在高超聲速流動(dòng)引起的嚴(yán)峻氣動(dòng)熱環(huán)境中，材料熱物性非線性會(huì)對(duì)傳熱過程產(chǎn)生顯著影響。目前對(duì)于該問題研究的相關(guān)報(bào)道還不多見。下面主要針對(duì)材料熱導(dǎo)率、比熱容等熱物性的非線性對(duì)高超聲速流動(dòng)傳熱過程的影響進(jìn)行了研究。

共進(jìn) 行 4 組計(jì) 算，分 別 為：①ZrB2－SiC，Cp（T）和k（T）均采用ZrB2－30vol．%SiC的物性，作為基礎(chǔ)物性和參照對(duì)象；②ZrB2，Cp（T）和k（T）均采用ZrB2的物性，考察熱導(dǎo)率和比熱物性同時(shí)作用的影響；③ZrB2－SiC＿Cp，k（T）取基礎(chǔ)物性，Cp（T）采用 ZrB2－30vol．%SiC的質(zhì)量分?jǐn)?shù)平均計(jì)算得到的比熱值，考察比熱在合理的物性偏差內(nèi)的影響；④ZrB2－SiC＿k，Cp（T）取基礎(chǔ)物性，k（T）采用ZrB2的物性，考察熱導(dǎo)率對(duì)傳熱的影響。計(jì)算所用的各物性參數(shù)如圖11所示。

圖14給出的是以上不同物性參數(shù)情況下計(jì)算的表面溫度對(duì)比。可以看出在初始的一段時(shí)間內(nèi)物性參數(shù)的變化對(duì)熱響應(yīng)的影響甚微，基本可以忽略不計(jì)，差別在一定時(shí)間后即結(jié)構(gòu)達(dá)到一定溫度（大于800K）時(shí)開始顯現(xiàn)。整體上看，在90s時(shí)間內(nèi)熱導(dǎo)率和比熱對(duì)材料表面最終溫度的影響有限。

通過曲線ZrB2－SiC與曲線 ZrB2－SiC＿k的對(duì)比可以看出，較高的熱導(dǎo)率使頭部熱量得以向低溫區(qū)轉(zhuǎn)移，使其溫度降低了約40K，變化量為3．6%；而在溫度較低的肩部，各組計(jì)算結(jié)果差異很小，熱導(dǎo)率造成的差異僅為8K，降低了1%。

通過曲 線 ZrB2－SiC 與曲線 ZrB2－SiC＿Cp的對(duì)比可以看出，使用質(zhì)量分?jǐn)?shù)平均比熱值計(jì)算的表面溫度低于對(duì)照值。這是因?yàn)橛少|(zhì)量分?jǐn)?shù)平均預(yù)測(cè)的比熱值高于作為基礎(chǔ)物性的ZrB2－30vol．%SiC實(shí)際測(cè)量值，因此每升高一度所需要的能量更多，因此在同樣來流條件下，高比熱容能夠降低材料表面的溫度。但是從結(jié)果可以看出，材料表面溫度對(duì)合理范圍內(nèi)的比熱值變化并不敏感，在超出比熱測(cè)量誤差±3%范圍、比熱值平均相差50J／（kg·K）的情況下，造成的溫度差異在頭部最大為11K，偏差為0．98%，肩部?jī)H為11．3 K，偏差為1．4%。

圖14 不同熱物性參數(shù)作用時(shí)的表面溫度分布對(duì)比Fig．14 Comparison of surface temperature distribution for different thermophysical properties

最后，曲線ZrB2－SiC與曲線ZrB2的對(duì)比顯示，ZrB2的熱導(dǎo)率高于ZrB2－SiC，使其頭部溫度有低于ZrB2－SiC的趨勢(shì)，而ZrB2的比熱低于ZrB2－SiC，有使其溫度高于ZrB2－SiC的趨勢(shì)。二者共同作用，其表面溫度介于比熱和熱導(dǎo)率單獨(dú)作用的曲線之間。

6 結(jié) 論

1）傳統(tǒng)的先流場(chǎng)再結(jié)構(gòu)的耦合方法無法得到精確的氣動(dòng)熱力載荷環(huán)境，也無法正確評(píng)價(jià)熱防護(hù)材料及其結(jié)構(gòu)的服役特征。采用高超聲速流動(dòng)和結(jié)構(gòu)傳熱耦合的方法，對(duì)高超聲速飛行器的氣動(dòng)熱力環(huán)境和結(jié)構(gòu)熱力響應(yīng)的預(yù)測(cè)更符合物理實(shí)際，保證計(jì)算精度。

2）對(duì)于如超高溫陶瓷（UHTC）這類熱導(dǎo)率較高的新型防熱材料，必須考慮熱傳導(dǎo)效應(yīng)對(duì)氣動(dòng)熱環(huán)境和結(jié)構(gòu)熱響應(yīng)預(yù)測(cè)的影響?？紤]熱傳導(dǎo)效應(yīng)后表面溫度和熱流相差分別可達(dá)500K和330kW／m2。

3）對(duì)UHTC類材料熱物性非線性對(duì)高超聲速流動(dòng)傳熱過程的影響進(jìn)行了研究，結(jié)果表明比熱和熱導(dǎo)率在合理誤差范圍內(nèi)的變化對(duì)材料表面溫度響應(yīng)的影響有限。在比熱變化±3%范圍、比熱值平均相差50J／（kg·K）的情況下，造成的溫度差異最大僅為0．98%

［1］ 彭治雨，石義雷，龔紅明，等．高超聲速氣動(dòng)熱預(yù)測(cè)技術(shù)及發(fā)展趨勢(shì)［J］．航空學(xué)報(bào)，2015，36（1）：325－345．PENG Z Y，SHI Y L，GONG H M，et al．Hypersonic aeroheating prediction technique and its trend of development［J］．Acta Aeronautica et Astronautica Sinica，2015，36（1）：325－345（in Chinese）．

［2］ WIETING A R．Experimental study of shock wave interference heating on a cylindrical leading edge：NASA－TM－100484［R］．Washington，D．C．：NASA，1987．

［3］ WIETING A R，HOLDEN M S．Experimental study of shock wave interference heating on a cylindrical leading edge at Mach 6and 8：AIAA－1987－1511［R］．Reston：AIAA，1987．

［4］ THORNTON E A，DECHAUMPHAI P．Finite element prediction of aerothermal－structural interaction of aerodynamically heated panels：AIAA－1987－1610［R］．Reston：AIAA，1987．

［5］ DECHAUMPHAI P，THORNTON E A，WIETING A R．Flow－thermal－structural study of aerodynamically heated leading edges［J］．Journal of Spacecraft and Rockets，1989，26（4）：201－209．

［6］ LOEHNER R，YANG C，CEBRAL J，et al．Fluid－structure－thermal interaction using a loose coupling algorithm and adaptive unstructured grids：AIAA－1998－2419［R］．Reston：AIAA，1998．

［7］ MILLER B A，CROWELL A R，MCNAMARA J J．Loosely coupled time－marching of fluid－thermal－structural interactions： AIAA－2013－1666［R］． Reston： AIAA，1998．

［8］ 黃唐，毛國良，姜貴慶，等．二維流場(chǎng)、熱、結(jié)構(gòu)一體化數(shù)值模擬［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2000，18（1）：115－119．HUANG T，MAO G L，JIANG G Q，et al．Two dimensional coupled flow－thermal structural numerical simulation［J］．Acta Aerodynamica Sinica，2000，18（1）：115－119（in Chinese）．

［9］ 夏剛，劉新建，程文科，等．鈍體高超聲速氣動(dòng)加熱與結(jié)構(gòu)熱傳遞耦合的數(shù)值計(jì)算［J］．國防科技大學(xué)學(xué)報(bào)，2003，25（1）：35－39．XIA G，LIU X J，CHENG W K，et al．Numerical simulation of coupled aeroheating and solid heat penetration for a hypersonic blunt body［J］．Journal of National University of Defense Technology，2003，25（1）：35－39（in Chinese）．

［10］ 桂業(yè)偉，袁湘江．類前緣防熱層流場(chǎng)與熱響應(yīng)耦合計(jì)算研究［J］．工程熱物理學(xué)報(bào)，2002，23（6）：733－735．GUI Y W，YUAN X J．Numerical simulation on the coupling phenomena of aerodynamic heating with thermal response in the region of the leading edge［J］．Journal of Engineering Thermophysics，2002，23（6）：733－735（in Chinese）．

［11］ 耿湘人，張涵信，沈清，等．高速飛行器流場(chǎng)和固體結(jié)構(gòu)溫度場(chǎng)一體化計(jì)算新方法的初步研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20（4）：422－427．GENG X R，ZHANG H X，SHEN Q，et al．Study on an integrated algorithm for the flowfields of high speed vehicles and the heat transfer in solid structures［J］．Acta Aerodynamica Sinica，2002，20（4）：422－427 （in Chinese）．

［12］ ZHAO X L，SUN Z X，TANG L S，et al．Coupled flowthermal－structural analysis of hypersonic aerodynamically heated cylindrical leading edge［J］．Engineering Applications of Computational Fluid Mechanics，2011，5（2）：170－179．

［13］ 張兵，韓景龍．多場(chǎng)耦合計(jì)算平臺(tái)與高超聲速熱防護(hù)結(jié)構(gòu)傳熱問題研究［J］．航空學(xué)報(bào)，2011，32（3）：400－409．ZHANG B，HAN J L．Multi－field coupled computing platform and thermal transfer of hypersonic thermal protection structures［J］．Acta Aeronautica et Astronautica Sinica，2011，32（3）：400－409（in Chinese）．

［14］ ZHANG S T，CHEN F，LIU H．Integrated fluid－thermal－structural analysis for predicting aerothermal environment of hypersonic vehicles：AIAA－2014－1394［R］．Reston：AIAA，2014．

［15］ 董維中，高鐵鎖，丁明松，等．高超聲速飛行器表面溫度分布與氣動(dòng)熱耦合數(shù)值研究［J］．航空學(xué)報(bào)，2015，36（1）：311－324．DONG W Z，GAO T S，DING M S，et al．Numerical study of coupled surface temperature distribution and aerodynamic heat for hypersonic vehicles［J］．Acta Aeronautica et Astronautica Sinica，2015，36（1）：311－324（in Chinese）．

［16］ SILBERBERG M．Chemistry：The molecular nature of matter and change［M］．New York：McGraw－Hill，2009．

［17］ PARK C，LEE S H．Validation of multi－temperature nozzle flow code NOZNT：AIAA－1993－2862［R］．Reston：AIAA，1993．

［18］ FRANCESE A．Numerical and experimental study of UHTC materials for atmospheric re－entry［D］．Napoli：University of Napoli，2007：94－95．

［19］ ANDERSON J D．Hypersonic and high temperature gas dynamics［M］．New York：McGraw－Hill，2006：697－699．

［20］ GUPTA R N，YOS M，THOMPSON R A．A review of reaction rates and thermodynamic and transport properties for the 11－species air model for chemical and thermal nonequilibrium calculations to 30000K：NASA TM－101528［R］．Washington，D．C．：NASA，1989：13－14．

［21］ 閻超，禹建軍，李君哲．熱流CFD計(jì)算中格式和網(wǎng)格效應(yīng)若干問題研究［J］．空氣動(dòng)力學(xué)報(bào)，2006，24（1）：125－130．YAN C，YU J J，LI J Z．Scheme effect and grid dependency in CFD computations of heat transfer［J］．Acta Aerodynamica Sinica，2006，24（1）：125－130（in Chinese）．

［22］ FERRERO P，D’AMBROSIO D．A numerical method for conjugate heat transfer problems in hypersonic flows：AIAA－2008－4247［R］．Reston：AIAA，2008．

［23］ MURTY M C，MANNA P，CHAKRABORTY D．Conjugate heat transfer analysis in high speed flows［J］．Proceedings of the Institution of Mechanical Engineers，Part G：Journal of Aerospace Engineering，2013，227（10）：1672－1681．

［24］ ZIMMERMANN J W．Thermophysical properties of ZrB2and ZrB2－SiC ceramics［J］．Journal of the American Ceramic Society，2008，91（5）：1551－2916．

Multi－field coupling numerical analysis of aerothermal environment and structural heat transfer of hypersonic vehicles

ZHOU Yinjia，MENG Songhe＊，XIE Weihua，YANG Qiang
Center for Composite Materials，Harbin Institute of Technology，Harbin 150080，China

A coupling framework of hypersonic flow，heat transfer and structural response is presented in this paper in order to accurately predict aerodynamic environment，extreme aerothermal environment，as well as thermal and structural response for hypersonic flight．Multi－field coupling analysis is implemented in conjunction with the hypersonic chemical nonequilibrium computational fluid dynamics（CFD）solver and the fully coupled thermo－structural finite element method（FEM）solver by using partition algorithm，with a real time data exchange between non－matched meshes．This coupling method is validated by comparison with the experiment of a turbulent flow over a circular cylinder and good agreements with experiment are achieved．Coupling analysis of ultra high temperature ceramics（UHTC）is also conducted，and the effects of thermal conductivity on the prediction of aerothermal environment and structural thermal response have been considered．The results show that the effects of structural internal heat conduction on aerothermal environment and surface temperature distribution cannot be neglected for a structure with high thermal conductivity and complex geometry．At last，the effects of non－linear thermophysical properties（specific heat and thermal conductivity）of UHTC on hypersonic flow and heat transfer process have been studied．The results show that the surface temperature of structure is not sensitive to thermal conductivity and specific heat when thermal conductivity and specific heat are in a limit of allowable error．

hypersonic vehicles；multi－field coupling；aeroheating；numerical simulation；thermal protection

2015－09－08；Revised：2015－11－22；Accepted：2016－02－15；Published online：2016－02－23 09：00

URL：www．cnki．net／kcms／detail／11．1929．V．20160223．0900．008．html

s：National Natural Science Foundation of China（11272107）；National Basic Research Program of China（2015CB655200）

V434＋．11

A

1000－6893（2016）09－2739－10

10．7527／S1000－6893．2016．0040

2015－09－08；退修日期：2015－11－22；錄用日期：2016－02－15；網(wǎng)絡(luò)出版時(shí)間：2016－02－23 09：00

www．cnki．net／kcms／detail／11．1929．V．20160223．0900．008．html

國家自然科學(xué)基金（11272107）；國家“973”計(jì)劃（2015CB655200）

＊通訊作者．Tel．：0451－86417560 E－mail：mengsh＠hit．edu．cn

周印佳，孟松鶴，解維華，等．高超聲速飛行器熱環(huán)境與結(jié)構(gòu)傳熱的多場(chǎng)耦合數(shù)值研究［J］．航空學(xué)報(bào)，2016，37（9）：2739－2748．ZHOU Y J，MENG S H，XIE W H，et al．Multif－ield coupling numerical analysis of aerothermal environment and structural heat transfer of hypersonic vehicles［J］．Acta Aeronautica et Astronautica Sinica，2016，37（9）：27392－748．

周印佳 男，博士研究生。主要研究方向：高超聲速氣動(dòng)熱與熱防護(hù)、高超聲速多物理場(chǎng)耦合。Tel．：0451－86402432 E－mail：12B918024＠hit．edu．cn孟松鶴 男，博士，教授，博士生導(dǎo)師。主要研究方向：極端環(huán)境與材料響應(yīng)的耦合問題、新型熱防護(hù)系統(tǒng)及材料的設(shè)計(jì)與優(yōu)化、先進(jìn)復(fù)合材料結(jié)構(gòu)分析與評(píng)價(jià)。Tel．：0451－86417560 E－mail：mengsh＠hit．edu．cn

＊Corresponding author．Tel．：0451－86417560 E－mail：mengsh＠hit．edu．cn