主動冷卻點(diǎn)陣夾層防熱結(jié)構(gòu)溫度響應(yīng)計算模型*

彭世彬，郭 瑞，馮上升，金 峰

（1.西安交通大學(xué) 機(jī)械結(jié)構(gòu)強(qiáng)度與振動國家重點(diǎn)實(shí)驗(yàn)室，西安 710049；2.西安交通大學(xué) 仿生工程與生物力學(xué)研究所，西安 710049；3.西安交通大學(xué) 生物醫(yī)學(xué)信息工程教育部重點(diǎn)實(shí)驗(yàn)室，西安 710049）

引言

點(diǎn)陣夾層結(jié)構(gòu)由于同時具備承受機(jī)械載荷和熱載荷的雙重功能而在航空航天領(lǐng)域具有廣闊的應(yīng)用前景[1-3]，如超高音速飛行器的熱防護(hù)部件[4]，火箭發(fā)動機(jī)的燃燒室[5]，以及航空母艦上的偏流板[6]等.現(xiàn)有關(guān)于點(diǎn)陣夾層結(jié)構(gòu)熱學(xué)特性的研究主要集中在穩(wěn)態(tài)傳熱，然而在上述應(yīng)用中，結(jié)構(gòu)表面的熱載荷通常是隨時間急劇變化的，因此，探究點(diǎn)陣夾層結(jié)構(gòu)在非穩(wěn)態(tài)熱載荷作用下的溫度響應(yīng)十分必要.

相比于孔結(jié)構(gòu)隨機(jī)排列的開孔金屬泡沫，具有周期性單胞的點(diǎn)陣夾層結(jié)構(gòu)在力學(xué)載荷作用下芯體變形以桿件拉伸為主導(dǎo)，因此在相同的孔隙率下具有更優(yōu)異的比強(qiáng)度/剛度[7].在傳熱方面，點(diǎn)陣結(jié)構(gòu)既增加了散熱面積，又可在芯體通道內(nèi)產(chǎn)生流動漩渦并促進(jìn)流體擾動，從而增大表面對流換熱系數(shù)；而且其內(nèi)部流動更加有序，比隨機(jī)開孔金屬泡沫具有更低的壓降[7].因此，點(diǎn)陣夾層結(jié)構(gòu)在承載-散熱多功能應(yīng)用場合更具綜合優(yōu)勢.Kim 等發(fā)現(xiàn)四面體夾芯板的散熱效率是空通道的6 倍左右[8].Lu 等利用PIV 粒子成像技術(shù)揭示了四面體點(diǎn)陣夾芯板中的流場分布，證明夾芯桿與面板之間形成的渦流是引起傳熱強(qiáng)化的主要原因[9].Yan 等通過沖壓金屬薄板制作了一種輕質(zhì)X 型點(diǎn)陣結(jié)構(gòu)并研究了其對流換熱性能，研究發(fā)現(xiàn)X 型點(diǎn)陣獨(dú)特的結(jié)構(gòu)形態(tài)導(dǎo)致大規(guī)模的螺旋主流，它與二次流相互作用提高了換熱能力[10-11].截至到目前，大部分研究集中在對點(diǎn)陣結(jié)構(gòu)內(nèi)部流動換熱機(jī)理的探究方面，而對其在非穩(wěn)態(tài)熱載荷作用下的溫度響應(yīng)鮮有研究.此外，在理論模擬方面，過去的研究主要針對均勻熱邊界條件，采用CFD 數(shù)值方法模擬一個或者一排單胞內(nèi)的流動傳熱[12]，然而對于工程中的非均勻熱邊界條件問題，需要的計算區(qū)域更大，計算非常耗時.現(xiàn)有理論建模一般是基于翅片理論，將點(diǎn)陣夾芯桿中的熱傳導(dǎo)簡化為表面對流的一維翅片導(dǎo)熱[8]，翅片根部的溫度認(rèn)為等于連接面板表面的平均溫度；實(shí)際上，由于芯體桿件根部的面積較小，熱量在透過面板的過程中逐漸向桿件根部聚集，導(dǎo)致桿件根部的溫度要低于連接面板表面的單元平均溫度，而現(xiàn)有的研究都忽略了這一重要影響因素.

綜上，本文以金字塔點(diǎn)陣結(jié)構(gòu)為研究對象，研究了主動冷卻下承受非穩(wěn)態(tài)熱載荷的點(diǎn)陣夾層結(jié)構(gòu)及內(nèi)部流體的溫度響應(yīng)問題，建立了夾層結(jié)構(gòu)導(dǎo)熱與冷卻劑對流耦合的理論模型，利用分離變量法推導(dǎo)了夾芯桿根部與連接面板之間收縮熱阻的近似解析解，并開展數(shù)值模擬驗(yàn)證了理論模型.

1 理論模型

1.1 模型介紹

圖1 為本文研究的金字塔型點(diǎn)陣金屬夾芯結(jié)構(gòu)的對流換熱示意圖.考慮該金字塔型點(diǎn)陣夾芯結(jié)構(gòu)在主動冷卻作用下承受非穩(wěn)態(tài)熱流密度時的溫度響應(yīng)問題.外界熱量由上面板輸入，下面板底部處于絕熱狀態(tài)，冷卻液以初始溫度Tcool= 300 K 流入夾芯板帶走熱量.上、下面板厚度為tb= 0.8 mm，夾芯高度為H= 5 mm，夾芯桿直徑為d= 1 mm，夾芯桿與下面板的夾角為θ = 45°.

1.2 非穩(wěn)態(tài)傳熱控制方程

1.3 控制方程離散化

1.3.1 結(jié)構(gòu)體控制方程離散化

針對點(diǎn)陣結(jié)構(gòu)周期排列的特點(diǎn)，以流動方向上每個單胞的上、下面板為面板離散的最小單元，以單胞內(nèi)流體區(qū)域?yàn)榱黧w域的最小離散單元，每個單胞內(nèi)的夾芯桿沿著長度方向進(jìn)行區(qū)域離散，點(diǎn)陣夾層結(jié)構(gòu)及內(nèi)部流體的離散示意圖如圖2所示.

圖2 點(diǎn)陣夾層結(jié)構(gòu)及內(nèi)部流體的離散示意圖Fig.2 Discretization of the lattice sandwich structure with internal fluid

1.3.2 流體控制方程離散化

在MATLAB 中對結(jié)構(gòu)和流體的非穩(wěn)態(tài)傳熱代數(shù)方程進(jìn)行迭代求解，具體的過程如下：

1)選取合適的時間步長Δt和空間步長Δη；

2)為結(jié)構(gòu)域和流體域賦予初始溫度(300 K)；

3)計算下一時刻上、下面板和夾芯桿及流體域離散方程中的附加熱源項(xiàng)；

4)計算下一時刻結(jié)構(gòu)域和流體域中離散單元內(nèi)的溫度；

5)將下一時刻結(jié)構(gòu)域和流體域的溫度作為初始溫度，重復(fù)步驟3)～5)，直至整個非穩(wěn)態(tài)過程結(jié)束.

計算在滿足收斂條件的情況下，取Δt= 0.004 s，Δη = 0.35 mm 足夠滿足計算精度.

1.4 確定對流換熱系數(shù)

上述離散后的代數(shù)方程中出現(xiàn)的兩個對流換熱系數(shù)hb和hfin是未知的，因此理論模型并不完整.為此，采取基于單胞的周期性流動傳熱模型進(jìn)行仿真求解，具體做法如下.

7）如果 len（Lborder[]）＞len（Rborder[]），就把 w-1加入到右邊界列表中（防止沒有檢測到最后一個字符的右邊界）。進(jìn)入步驟8）。

首先根據(jù)結(jié)構(gòu)對稱性建立只有半個單胞的周期性流動傳熱模型.將流體進(jìn)出口設(shè)置為周期性邊界條件，流體進(jìn)口的平均速度在u0= 0.2～1.0 m/s 范圍內(nèi)變化，固體域及流體域左右側(cè)均設(shè)置為對稱邊界，流固交界面設(shè)為無滑移耦合傳熱邊界條件，單胞上面板給定恒定的熱流密度128500 W/m2，其他壁面設(shè)置為絕熱邊界.在Fluent Meshing 中對當(dāng)前模型進(jìn)行網(wǎng)格劃分，網(wǎng)格類型采用多面體型，對壁面邊界層處進(jìn)行加密處理確保y+<1，湍流模型選用k-ω SST，該模型在計算復(fù)雜分離流時具有優(yōu)越性，能夠很好地預(yù)測剪切流動中的流動分離現(xiàn)象[14]，因此非常適用于當(dāng)前流動模型中漩渦流的模擬.采用SIMPLE 算法耦合壓力和速度進(jìn)行數(shù)值分析，對流項(xiàng)的離散采用二階迎風(fēng)格式，動量方程的迭代收斂準(zhǔn)則為10?3，能量方程的迭代收斂準(zhǔn)則為10?7.圖3中展示了模型的計算域、邊界條件和網(wǎng)格.

圖3 周期性流動傳熱數(shù)值模擬：(a)計算域及邊界條件；(b)網(wǎng)格Fig.3 Numerical simulation of periodic flow and heat transfer:(a)the computational domain and boundary conditions;(b)the mesh

定義Red，hb和hfin為[9]

圖4 不同Red 下，夾芯桿表面的對流換熱系數(shù)Fig.4 Heat transfer coefficients on the surface of lattice struts under different Red values

1.5 求解收縮熱阻

面板-點(diǎn)陣之間的收縮熱阻的求解可以等效為以下問題：如圖5所示，一塊矩形板頂部承受均勻熱流，四周絕熱，底部分別受到不同程度的對流冷卻作用，即在面板內(nèi)表面的中間與翅片連接界面上的等效對流換熱系數(shù)為heff，而其他直接接觸冷卻流體區(qū)域的表面對流換熱系數(shù)為hb.由于模型的對稱性，可取圖5 中模型的1/4 進(jìn)行求解，因此該問題的控制方程和邊界條件可表示為

圖5 求解收縮熱阻的等效模型Fig.5 The equivalent model for solving the constriction thermal resistance

對于一般的情況，由于方程的復(fù)雜性，可在MATLAB 中進(jìn)行迭代求解，具體過程如下：

1)選擇所需的擴(kuò)展項(xiàng)數(shù)量(N×N)；

2)令cmn=dmn= 0(下表面均勻換熱)作為初始值；

3)將dmn代入式(45)和(47)計算新的A，B和cmn；

4)將cmn代入式(46)得到新的dmn；

經(jīng)過驗(yàn)證，N= 30 足夠滿足求解精度，并設(shè)置兩次迭代之間dmn的變化小于0.001.

2 數(shù)值模型

為驗(yàn)證當(dāng)前理論模型的準(zhǔn)確性，對所示的點(diǎn)陣結(jié)構(gòu)進(jìn)行數(shù)值模擬，考慮到計算資源的限制，沿流動方向的單胞數(shù)量取n= 10，但對理論計算模型n可以遠(yuǎn)大于10.結(jié)構(gòu)材料設(shè)為鋼，其密度ρs= 8 030 kg/m3，熱導(dǎo)率ks=16.27 W/(m·K)，定壓比熱容cps= 502.48 J/(kg·K).冷卻劑為航空煤油JP-7[18]，密度ρf= 800 kg/m3，定壓比熱容cpf= 2575 J/(kg·K)，熱導(dǎo)率kf= 0.11 W/(m·K)，黏度μf= 1.984 × 10?4kg/(m·s).

圖6 展示了模型的計算域及邊界條件：本文中首先按照1.4 小節(jié)中的方法對單胞進(jìn)行周期性流動傳熱模擬，為避免入口效應(yīng)，入口邊界上的速度直接讀取周期性流動傳熱模擬得到的周期邊界上的速度分布，出口采用壓力邊界，上面板外表面采用非穩(wěn)態(tài)熱流邊界，時間步長設(shè)置為1 s，兩側(cè)為對稱邊界，流固交界面采用無滑移耦合傳熱邊界條件，其余面均為絕熱邊界.網(wǎng)格類型、劃分方式及湍流模型均與單胞周期性模擬一致.為確保計算結(jié)果與網(wǎng)格數(shù)量無關(guān)，取加熱表面平均溫度為監(jiān)測值進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證，測試網(wǎng)格數(shù)從120 萬增加到400 萬，發(fā)現(xiàn)在網(wǎng)格數(shù)量超過340 萬時計算結(jié)果已無明顯變化，如圖7所示，因此最終采取的網(wǎng)格數(shù)為340 萬.

圖6 含多個單胞點(diǎn)陣夾層結(jié)構(gòu)數(shù)值模擬計算域及邊界條件Fig.6 The numerical simulation domain and boundary conditions for the multi-cell lattice sandwich structure

圖7 加熱面平均溫度隨網(wǎng)格數(shù)量的變化曲線Fig.7 The variation of the average temperature of heating surface with the number of meshes

3 結(jié)果與討論

3.1 理論與仿真對比

本文算例中采用可重復(fù)使用運(yùn)載器(RLV)進(jìn)入大氣時表面產(chǎn)生的瞬時熱通量作為外加熱源，以確定理論模型的適用性.RLV 表面的氣動加熱曲線如圖8所示，瞬時熱通量從0 s 變化到2200 s，其中qmax= 128500 W/m2是熱通量的最大強(qiáng)度，2200 s 是再入過程的最后時間.RLV 表面的熱通量分布假定為分段多項(xiàng)式，各階段具體的表達(dá)式見文獻(xiàn)[19]中的表1.

圖8 RLV 再入過程中表面的入射熱通量隨時間的變化曲線Fig.8 The incident heat flux profile vs.the re-entry time on the RLV surface

將整個點(diǎn)陣結(jié)構(gòu)的最大溫度Ts,max(沿流動方向最后一個單胞的上面板外表面平均溫度)以及流體出口的溫度Tf,out作為考核目標(biāo)并進(jìn)行理論與仿真結(jié)果對比.圖9(a)為不同流體平均進(jìn)口速度(u0= 0.2～1.0 m/s)下點(diǎn)陣結(jié)構(gòu)的最大溫度隨時間的變化曲線，從圖中可以看出，理論結(jié)果與CFD 仿真結(jié)果吻合良好，最大誤差低于1%.圖9(b)為流體出口溫度隨時間的變化曲線，在當(dāng)前的驗(yàn)證模型中，由于總尺寸較小，同時冷卻劑流速較大，導(dǎo)致與入口溫度相比，出口溫度并無明顯上升，但在整個非穩(wěn)態(tài)過程中理論結(jié)果與CFD 結(jié)果高度吻合.

圖9 結(jié)構(gòu)最大溫度與流體出口溫度的變化曲線：(a)結(jié)構(gòu)最大溫度；(b)流體出口溫度Fig.9 Variations of the maximum temperature of the sandwich structure and the outlet fluid temperature with time:(a)the maximum temperature of the sandwich structure;(b)the outlet fluid temperature

3.2 收縮熱阻對模型精度的影響

圖10 給出了u0= 0.2 m/s 時，不同的恒定熱流密度下點(diǎn)陣結(jié)構(gòu)所能達(dá)到的最大溫度，從圖中可以看出，忽略收縮熱阻使得計算結(jié)果偏低，并且隨著熱流密度不斷增大，忽略收縮熱阻使得計算結(jié)果造成的誤差呈放大趨勢，因此在外界熱載荷較高時不能忽略收縮熱阻的影響.

圖10 不同熱流密度下結(jié)構(gòu)的最大溫度Fig.10 The maximum temperature of the structure under different heat flux densities

3.3 導(dǎo)熱率和面板厚度對收縮熱阻的影響

圖11 和12 分別給出了在不同流體進(jìn)口速度下收縮熱阻隨導(dǎo)熱率和面板厚度的變化曲線.結(jié)果表明，在不同的流速下，收縮熱阻隨著導(dǎo)熱率以及面板厚度的減小而增大，并且隨著流速增大而減小，因此在冷卻劑流速比較小時，忽略收縮熱阻會造成較大的誤差.

圖11 不同流體進(jìn)口速度下收縮熱阻隨導(dǎo)熱率的變化曲線Fig.11 Variations of the constriction thermal resistance with the thermal conductivity under different inlet fluid velocities

圖12 不同流體進(jìn)口速度下收縮熱阻隨面板厚度的變化曲線Fig.12 Variations of the constriction thermal resistance with the facesheet thickness under different inlet fluid velocities

4 結(jié)論

本文建立了可用于非均勻、非穩(wěn)態(tài)復(fù)雜熱載荷下快速預(yù)測點(diǎn)陣夾層結(jié)構(gòu)溫度響應(yīng)的非穩(wěn)態(tài)傳熱理論模型，首次考慮了夾芯桿根部與連接面板之間收縮熱阻對點(diǎn)陣夾層結(jié)構(gòu)傳熱的影響，并用數(shù)值方法進(jìn)行了驗(yàn)證.得到的結(jié)論如下：

1)理論模型能夠準(zhǔn)確預(yù)測金字塔點(diǎn)陣結(jié)構(gòu)以及內(nèi)部流體的溫度場變化，理論結(jié)果與模擬結(jié)果吻合良好，最大誤差不超過1%，在考核算例下計算耗時僅為CFD 模擬的4%.

2)隨著外加熱流密度不斷增大，忽略收縮熱阻使得計算結(jié)果造成的誤差不斷增大.

3)收縮熱阻隨著導(dǎo)熱率以及面板厚度的減小而增大，隨著冷卻流體進(jìn)口流速增大而減小.因此，在高熱流密度、薄面板厚度以及低冷卻劑流速條件下預(yù)測點(diǎn)陣夾層結(jié)構(gòu)傳熱特性，必須考慮面板與點(diǎn)陣芯體之間的收縮熱阻.