各向異性Kelvin泡沫胞元高徑比對其流動傳熱特性的影響

張 超，孔祥壯，杜雁霞，王 嫻，肖光明，*

(1.西安交通大學 航天航空學院，機械結構強度與振動國家重點實驗室，陜西省先進飛行器服役環境與控制重點實驗室，西安 710049；2.空天飛行空氣動力科學技術全國重點實驗室，綿陽 621000)

0 引 言

泡沫材料具有比表面積大、質量輕、高導熱系數等[1]優點，被廣泛應用于高超聲速飛行器熱管理系統中[2]，如用于高超聲速飛行器前緣主動冷卻結構的SiC泡沫[3]以及金屬泡沫[4]。然而，真實泡沫材料的微細結構極度復雜無序，且具有明顯的各向異性[5-6]，使得對泡沫材料內部流動傳熱特性的精確預測有較大困難。對于泡沫材料的微細結構表征，X射線斷層成像技術是較為精確的方法[7-9]。該方法能完全還原泡沫骨架微細結構及孔隙分布特征，有非常高的精度。然而，受限于加工過程中一些結構參數的不可控，如前文提到的孔隙結構，使X射線斷層成像技術無法高效地運用到結構優化設計的相關研究中。因此，為了進一步優化泡沫材料的傳熱性能，需要對真實泡沫材料微細結構進行關鍵特征提取，建立等效模型并對關鍵結構參數進行定量分析。

最初，許多學者不考慮泡沫材料的結構各向異性，采用各向同性的代表性單元表征泡沫材料微細結構。應用最為廣泛的是Kelvin泡沫[10]，該結構提取了泡沫材料中的韌帶、節點和孔3個關鍵參數，在泡沫材料的等效導熱系數理論預測問題上具有較高精度，最早由Boomsma和Poulikakos[11]將其運用到泡沫材料等效導熱系數的理論預測中，經過Dai等[12]對熱流傳遞方向的修正，形成一套完整的理論體系。在此基礎上，Yao等[13]通過改變節點形狀及分布進一步優化了該模型，并采用實驗進行驗證，最大相對誤差小于10%。上述研究中，泡沫結構各向異性對傳熱的影響被忽略；然而，自從Ranut等[14]通過ANSYS CFX計算真實泡沫結構沿3個正交方向的對流導熱系數和壓降，發現3個正交方向的導熱系數和壓降均存在較大差異，泡沫材料結構各向異性對傳熱各向異性的影響逐漸被一些學者關注。

泡沫結構各向異性產生的主要原因是在發泡過程中，由于重力和黏性力的影響，使泡沫骨架產生拉伸和變形[15]。而對于泡沫結構各向異性的描述，常用各個方向的孔徑長度比值進行定量表征[5,16]。在不考慮孔內介質的條件下，Iasiello等[8]采用形態學分析，統計出3個正交方向的胞元直徑，以高斯分布的中值作為各向異性的表征量。此外，采用COMSOL Multiphysics計算3個正交方向的等效導熱系數。結果表明，結構各向異性與導熱各向異性正相關。當考慮孔內流體流動時，Iasiello等[17]發現，當骨架自身為熱源對孔內流體進行加熱時，拉伸方向的傳熱系數最小，說明傳熱各向異性不僅與結構各向異性有關，也與邊界條件強相關。

上述關于各向異性的研究主要是針對真實結構的，然而，真實結構無法控制孔隙結構參數，無法有效對結構進行設計和優化。在此基礎上，本團隊[18]結合Dai等[12]的各向同性Kelvin泡沫理論模型和Iasiello等[8]的各向異性結構表征，通過理論推導，建立了各向異性Kelvin泡沫等效導熱系數理論預測公式，并且將理論預測結果和真實結構的數值計算結果進行對比，發現兩者吻合良好，進一步定量證明了導熱各向異性與結構各向異性正相關，并且指出Kelvin泡沫拉伸方向的導熱系數最大，且最大值數值可達到胞元最短方向等效導熱系數的1.47倍。

根據以上表述，如果將沿主流方向的尺寸作為直徑，加熱方向尺寸為高度，Kelvin泡沫的各向異性可用不同胞元高徑比來表征。胞元高徑比的改變影響各向同性Kelvin泡沫的傳熱特性，那么是否可利用胞元高徑比去設計和優化各向異性Kelvin的傳熱特性？對于該問題，Sun等[19]已經開展了部分工作，將Kelvin泡沫沿通道主流方向拉伸1.5倍和2.0倍，此時，胞元高徑比分別為2/3和1/2，最終結構的比表面積隨著胞元高徑比的減小而增大。在數值計算中，與Iasiello等[17]的邊界條件相同，骨架自身為熱源對孔內流體進行加熱，分析不同速度下胞元高徑比為1、2/3、1/2時Kelvin泡沫內的流動傳熱特性。結果表明，胞元高徑比為1/2時的導熱系數最低，傳熱系數與Kelvin泡沫比表面積成反比。該現象似乎與傳統認知相反，事實上，主要原因是Sun等的研究更偏向于工程應用，在計算中保持速度大小一致，在該條件下胞元高徑比越小，流道高度越低，導致雷諾數較低；而導熱系數隨著雷諾數的增大迅速增大，特別是在非達西區，Sun等[19]的計算又剛好位于非達西區。此外，Sun等未考慮骨架自身的導熱，將流-固邊界定義為等熱流邊界，只對孔隙中流體的流動傳熱進行數值模擬。事實上，在實際工程應用中，骨架的導熱增強是十分重要的。因此，在同時考慮泡沫骨架導熱和Kelvin泡沫內孔隙尺度強迫對流換熱的共軛傳熱時，胞元高徑比對傳熱性能的影響可在保持雷諾數一定的條件下作進一步分析。

基于上述問題和現狀，本研究在保證流態相同的條件下，以各向異性Kelvin泡沫為研究對象，開展不同胞元高徑比對各向異性Kelvin泡沫流動傳熱特性的影響研究，著重分析胞元高徑比變化對流動傳熱內在影響機制，為飛行器熱管理系統中泡沫結構設計提供技術支撐。

1 Kelvin泡沫幾何參數

本文選擇如圖1所示的5種不同泡沫進行研究。其中，Kelvin泡沫的長度x和寬度z保持一致，記為Kelvin泡沫的直徑，其尺寸均為D。高度H（y）與直徑D之比H/D則定義為胞元高徑比。當胞元高徑比H/D=1.0時，該結構為最為常用的各向同性Kelvin泡沫[20-22]。本文選取H/D= 0.5、0.75、1.5、2.0的各向異性Kelvin泡沫進行研究。

對于Kelvin泡沫，孔隙率是影響傳熱性能的主要參數之一，在本研究中孔隙率保持不變，因此通過調整Kelvin泡沫的韌帶半徑和節點尺寸，最終保證本研究的5種構型Kelvin泡沫的孔隙率ε均為0.912。此外，為了便于對比分析，保持流道方向長度D= 4 mm。經過上述尺寸設定，在孔隙率一定的條件下，5種構型Kelvin泡沫的比表面積αsf分別為1412.9、1189.2、1079.8、955.6、891.2 m2/m3，可以看出，隨著H/D不斷增大，比表面積不斷減小。表1給出了5種Kelvin結構的結構參數。

表1 Kelvin泡沫幾何參數Table 1 Geometrical parameter of Kelvin structure

2 計算模型與方法

2.1 計算模型

本文采用如圖2所示的計算模型，參考Iasiello等[16]的研究，沿通道方向布置4個Kelvin泡沫，其尺寸為 4D×D×H，流體通道的整體尺寸為Lx×Ly×Lz=6D×D×H。其中，通道沿y方向的高度分別為0.5D、0.75D、1.0D、1.5D、2.0D。Kelvin泡沫設置為導熱體，其導熱系數為ks，W/(m2·K)；通道流體的導熱系數為kf，W/(m2·K)。此外，4個Kelvin泡沫的上壁面為加熱面，保持高溫Th，K；左壁面為進口，進口冷卻流體的溫度為Tc，K，速度為u0，m/s；右壁面為出口，設為壓力出口邊界；前后壁面為對稱面，其余壁面保持絕熱。

圖2 計算模型示意圖Fig.2 Schematic of numerical model

在對流動傳熱特性進行分析之前，需要控制一些變量以保證流動及傳熱特性相似。對于本文的強迫對流問題，最重要的參數是雷諾數。對于該結構，通常選取流道高度H為參考長度，其雷諾數的表達式為：

式中 υ為動力黏度，Pa·s。為了保證流動相似（即Re相同），對于不同胞元高徑比結構，需要設置不同的初速度。

2.2 計算方法與流程

本文計算采用混合熱格子Boltzmann法（hybridthermal lattice Boltzmann method, HTLBM），即流場采用多松弛格子Boltzmann法（multi-relaxation time lattice Boltzmann method, MRT-LBM），溫度場采用有限差分法（finite difference method, FDM）。

2.2.1 流場計算

MRT的演化方程為[23]：

其中：r為空間坐標系；δt為時間步長；f是速度空間的密度分布函數；M是N×N階轉化矩陣，用于將f從速度空間線性地轉化為矩空間的m；S是非負的松弛矩陣；m和meq是矩空間的分布函數和平衡態分布函數。m、f和M之間的關系滿足：

其中，m、f和meq的表達式為

本文采用D3Q19模型，其轉化矩陣M為：

轉換后，矩空間內m和meq分別為：

其中松弛系數的取值為：se=sε=sπ=sv=1/τ,sm=sq=8(2-sv)/(8-sv)。

對于流-固邊界，采用LBM反彈格式進行處理。本文采用均勻笛卡爾網格進行數值計算，流-固邊界與坐標點無法保證完全重合，直接反彈會產生一定誤差，因此通過提高網格規模來降低邊界格式帶來的誤差，最大計算網格數達1.16億。

2.2.2 溫度場計算

對于流體域，無內熱源、常物性穩態流體對流-擴散方程：

其中，αf為流體熱擴散系數，m2/s。采用以下參數進行無量綱處理：

得到：

分別采用二階迎風和中心差分格式離散式（11）中的對流項和擴散項：

其中：

式（13）中：

式（14）中：

時間項采用二階顯式Runge-Kutta方法進行離散：

式中，n為計算時刻點。

固體域的能量方程，即骨架導熱方程如下：

其中αs為固體域熱擴散系數，m2/s。時間項和擴散項離散方式與流體域相同，不再贅述。

對于流-固邊界，為了和流場計算邊界進行統一，不再對邊界進行單獨建模，采用虛擬節點來表示流-固界面，如圖3所示。圖中，（s,i）和（f,i）為同一節點，即s-f為虛擬點。通過上述定義，流-固界面的離散方程為：

圖3 流-固邊界示意圖Fig.3 Schematic of fluid-solid boundary

其中，kr=ks/kf，ks和kf分別為Kelvin泡沫骨架和流體的導熱系數。

2.2.3 計算流程

對于含有多孔介質的復雜流動、傳熱和結構導熱的多場耦合問題，采用如圖4所示流場和溫度場同步求解的方法進行數值計算，即流場和溫度場采用相同的時間步長。此外，復雜結構（圖1）的網格劃分及邊界識別方法可參考之前的研究[18]，這里不再贅述。

圖4 計算流程圖Fig.4 Flow chart of numerical simulation

為了充分發揮LBM天然并行性的優勢，采用4張Tesla A100 GPU卡并行加速，網格數最大為1.16億，可獲得2153.5 MLUPS的計算效率。

2.3 流動傳熱特性參數定義

對于該類問題，衡量流動特性的參數有很多種，應用最為廣泛的是流道進出口壓降、滲透率和摩擦系數。

進出口壓降的計算公式為：

式中，pin和pout分別為進口、出口壓力，Pa。

由于壓降的方向是沿流動方向，所以滲透率同樣沿流動方向進行計算。根據達西定律，滲透率的表達式為：

其中，Kf為滲透率， μm2；υ 為動力黏度，Pa·s； 〈ux〉為沿流道流過泡沫的平均速度，m/s； ?p為壓力梯度，Pa。

對于本研究，采用LBM進行流場求解，因此，以x方向為例，最終的滲透率計算公式為[24]：

其中，Lf為通道長度，m；Nf為Lf長度對應的格子點數。

摩擦系數的定義如下：

其中， ?p為 壓差，Pa；Dh為水力直徑，m，在本文中，流道截面為正方形，水力直徑等于流道截面的邊長Lf；其中， 〈ux〉的求解見公式（22）。

為了定量化地表征泡沫材料在不同方向的傳熱性能，定義換熱系數h和體導熱系數hv如下：

其中，q為熱流密度，W/m2；Th為加熱面溫度，K；Tf為泡沫材料內部平均溫度，K。

此外，泡沫結構需要保持高導熱系數和低壓降。然而，當hv很高時，泡沫結構總是伴隨著高壓降。因此，需要對這些參數之間進行權衡，既滿足傳熱需求，又在合理的壓降范圍。一些學者采用綜合評價指標j/f1/3來衡量泡沫結構綜合性能指標。

其中，j為Colburn因子，Nuv為努塞爾數：

病史摘要:患者17歲,因“反復下腹脹痛1月,加重3天”入院?；颊?2歲初潮,月經基本規則,量中,伴痛經。近1月反復出現下腹脹痛,可忍受,伴腰痛,無異常陰道流血流液,自行購買止痛藥物后緩解。3天前患者自感上述癥狀加重,伴發熱,遂來我院門診就診。未婚有性生活1年,未孕。

2.4 程序驗證

為了驗證本文數值計算方法，采用如圖2所示的計算模型進行數值驗證。其中，Kelvin泡沫胞元高徑比選為1。為了與Xia等[25]的實驗結果進行對比，Kelvin泡沫骨架為銅，流體為空氣。本文計算結果和文獻對比結果如圖5所示。結果表明，整體上看，本研究的計算結果與文獻中的結果整體上吻合良好，但在局部存在一定差異，與Xia等[26]的實驗最大相對誤差為7.63%。主要原因是本文采用的是Kelvin結構，而Xia等[25]采用的是真實泡沫銅，兩者微細結構存在一定差異。

圖5 體傳熱系數hv計算值與實驗值的對比Fig.5 Comparison of calculated and experimental values of body heat transfer coefficient hv

2.5 網格無關性驗證

由于本文中的模型尺寸不相同，在計算中分別對H/D= 0.5、0.75、1.0、1.5、2.0五組模型在Re= 1000的工況下進行網格無關性考核。結果如表2所示，表中采用不同條件下的高度H的網格數為參考。從表中可知，對于H/D= 0.5、0.75、1.0、1.5、2.0的五組模型，當高度方向的網格數分別達到154、140、156、206、182時，hv不再受網格變化的影響，因此，H/D= 0.5、0.75、1.0、1.5、2.0時， 高度方向的網格數分別取154、140、156、206、182。

表2 不同 H/D 和 高度方向網格數下的hvTable 2 The value of hv in different grids andH/D

3 結果與討論

本文在Sun等[20]的研究基礎上，重點分析Kelvin泡沫胞元高徑比對其流動傳熱特性的影響。與Sun等的研究不同的是，本文不僅保持了相同的雷諾數，而且考慮了多孔介質的復雜流動、傳熱和結構導熱的多場耦合，同時分析了結構參數Kelvin泡沫胞元高徑比的變化對流動傳熱的影響。相比較于Sun等的研究，本文的研究更全面地考慮了多種因素對流動傳熱的影響。

3.1 流動特性分析

圖6為Re=100 時、不同H/D下的流線圖。從圖中可以看出，隨著H/D的不斷增大，流線分布逐漸均勻，即骨架柱體對流線的擾動越來越小。特別是在通道的尾端，當H/D=0.5時，尾端流線受骨架的擾動發生明顯偏轉，但當H/D增大到2.0時，這種偏轉現象變得平緩，流線逐漸變得平直。原因主要是骨架柱體對流體的擾動程度主要受壓力和雷諾數的影響，隨著流通面積的逐漸增大，相同流態下（Re相等），壓差變得越來越小。為了進一步分析，沿流道中心（y=0.5H，z=0.5D），取三種條件下的壓力沿x方向的分布如圖7所示。為了便于比較，采用通道長度6D對x進行無量綱化。

圖6 不同H/D 時的流線圖Fig.6 The streamlines at different H/D

圖7 H/D = 0.5、1.0、2.0時流道中心（ y = 0.5H、z = 0.5D ）上的壓力隨x的變化Fig.7 Variation of pressure on the center of the flow channel(y = 0.5H and z = 0.5D) with x at H/D = 0.5, 1.0, 2.0

從圖7中可以看出，在x方向的相同位置，以進口壓力為參考，隨著H/D的不斷減小，壓差逐漸增大。此外，在骨架位置壓差的分布呈周期性下降，取一個單胞分析，H/D= 0.5、1.0和2.0時的壓差分別為0.00066p0、0.00135p0和0.00374p0。H/D從2.0減小到0.5，壓差增大55.67倍。因此，該結果也進一步解釋了圖6中，流線由不均勻分布轉為均勻分布的原因。此外，H/D=2.0時，壓力隨著x增大是連續下降的，骨架喉道區域的前端壓力幾乎不變，尾部是產生壓降的主要部位。隨著H/D減 小，特別是H/D=0.5時，壓力隨著x增大在骨架的喉道區域出現明顯的脈動，盡管產生壓降的主要部位仍位于喉道的尾部，但最大壓力同樣為喉道的尾部區，說明喉道尾部發生明顯阻塞。當H/D=2.0時，喉道前中部位置的壓力幾乎不隨x的增大發生明顯變化，當H/D減小到0.5時，喉道前中部的壓力隨x的增大出現先減后增的趨勢。產生該現象的主要原因是尾部的阻塞使得前部壓力降低，當流體流過喉道區域，流道的擴張使阻塞效應變小，使得流體壓力降低。而使壓力產生波動的主要原因是隨著H/D的減小，喉道內Kelvin泡沫頂部骨架和底部骨架的距離減小，骨架間流體的相互擾動也隨即增大，使得壓力沿流道方向的分布不再均勻。

圖8為不同H/D時z=0.5D截面上的無量綱速度分布云圖。在前文中提到，為了保證流態相似，在計算過程中保持雷諾數不變，因此無量綱速度的范圍保持一致?？梢钥闯?，z=0.5D截面上的速度分布隨著H/D的減小發生明顯的變化。圖中紅色區域為流速較高區域，不同H/D時，其分布形狀均近似為三角形。隨著H/D的減小，喉道進口位置的高度越來越低，三角形逐漸被壓扁，由近似為等邊三角形變化為明顯等腰的銳角三角形，且流道方向的銳度也越來越大，原因主要是喉道進口處出現的阻塞效應隨著喉道高度的減小而增大。此時，流體在喉道進口處被擠壓的強度最高，因此，相同喉道長度下，壓差隨著H/D的減小而增大，導致高流速區域面積占比增大。此外，定義等腰三角形的底角為 θu，如圖8所示。

表3 不同H/D時的θuTable 3 The value of θu at different H/D

圖9 流道中心（y = 0.5H、z = 0.5D）上的速度沿x分布Fig.9 Distribution of velocity along x at H/D = 0.5, 1.0, 2.0

為了進一步分析流場的分布特性，圖10給出了不同H/D時z=0.5D截面上的無量綱渦量分布。從圖中可知，除了流道進口位置外，對稱渦主要存在于各個喉道的進出口位置，其余僅在H/D= 0.5時，在結構的每個骨架的迎風面形成明顯的對稱渦。

圖10 H/D = 0.5、0.75、1.0、1.5、2.0時，z = 0.5D截面上，無量綱渦量分布云圖Fig.10 Dimensionless vorticity distribution on z = 0.5D surface at H/D = 0.5, 0.75, 1.0, 1.5, 2.0

與流場特性相似，在喉道位置渦中心線與豎直方向夾角發生明顯變化，定義該夾角為 θω，如圖10所示。不同H/D條 件下的 θω如表4所示，與速度場中的 θu相似，隨著H/D的 減小， θω也逐漸增大，增大的趨勢也較為接近。

表4 不同H/D條件下的θωTable 4 The value of θω at different H/D

對于該類Kelvin泡沫，無量綱參數D2/Kf是評價其流動特性的重要參數，其表達式中運用了滲透率的倒數，D2/Kf越小，Kelvin泡沫的滲透性越好，因此認為D2/Kf可表征流經Kelvin泡沫產生的阻力特性。為了更好地對比不同H/D條件下的滲透特性，分別計算Re=10 、 100、1000條件下的D2/Kf，如圖11所示。不同Re條件下，D2/Kf均隨著H/D的增大而減小，表明隨著H/D的增大，滲透性越好，主要原因是隨著H/D的增大，Kelvin泡沫沿著高度方向被拉伸，流動截面逐漸增大，壓差逐漸減小，流動阻力越來越小，滲透性也就越來越好；相同H/D條件下，D2/Kf隨著Re的增大而增大，表明流速越大，相同結構下流動阻力越來越大，骨架的滲透性越差；隨著H/D的增大，不同Re間D2/Kf的差異減小，表明H/D和Re對D2/Kf的影響均是非線性的。進一步說明，對于各向異性Kelvin泡沫，結構參數H/D與流態表征Re在實際工程設計中均不可忽視。

圖11 不同Re下D2/Kf隨H/D的變化Fig.11 The change of D2/Kf with H/D at Re = 10, 100 and 1 000

3.2 傳熱特性分析

3.1節中重點分析了Kelvin泡沫H/D對其流動特性的影響，本節將重點分析H/D對傳熱特性的影響。

圖12為Re=100 條件下，H/D= 0.5、1.0、2.0時，Kelvin骨架表面溫度分布以及3個方向中截面溫度分布。由于3個樣品的高度不同，選擇同一位置的單個細胞進行相互比較具有一定的局限性[20]，因此在進行骨架對比分析時選擇4個胞元。3種樣品的骨架溫度分布相似，溫度較低的部位主要位于流速較高的喉部。此外，當H/D=2.0時，由于頂部骨架背風處存在明顯的低速區，該區域的空氣溫度高于胞元其他區域。當H/D=0.5時，背風區域減小，低溫區域不僅位于喉部位置，骨架通道位置也形成明顯的低溫區。

圖12 Re = 100時Kelvin骨架表面以及x 、y 、z 方向中截面溫度分布：（a）H/D = 0.5；（b）H/D = 1.0；（c）H/D = 2.0Fig.12 The temperature distribution on the surface of the Kelvin skeleton，and middle sections in x, y and z directions at(a) H/D = 0.5，(b) H/D = 1.0 and (c) H/D = 2.0 when Re = 100

為了進一步分析上述現象，給出不同H/D條件下，z=0.5D截面上的溫度分布如圖13（a）所示。隨著H/D的增大，骨架底部區域的溫度場差異較小，但頂部區域的差異較為明顯。頂部區域的差異主要體現在以下3個方面：

圖13 Re = 100時，（a）不同H/D下，z = 0.5D截面的溫度分布；（b）H/D = 0.5、1.0、2.0時流道中心（y = 0.5H、z =0.5D）溫度沿x的分布Fig.13 Re = 100, (a) The temperature distribution on z = 0.5D surface under different H/D; (b) Temperature distribution of the center of the flow channel (y = 0.5H and z = 0.5D) along x at H/D= 0.5, 1.0 and 2.0

首先，喉道前部位置，當H/D= 0.5時，骨架迎風面幾乎無高溫區域，隨著H/D的增大，骨架沿高度方向的長度增大，迎風區域骨架邊緣的溫度明顯降低；與此同時，骨架背風區域的低溫區域占比變大，原因是背風區域形成了明顯的低速區，導致傳熱性能降低。

其次，在喉道的中間區域，當H/D為0.5時，在第一個骨架區域，等溫線與主流方向幾乎平行，但隨著H/D的增大，等溫線形成了明顯的拱形且凹面向下分布（圖13（a）中紅色虛線）。該現象說明胞元高徑比H/D影 響了傳熱主導方式。當H/D=2.0時，由于背風面存在流動低速區，導致強迫對流換熱強度較小。由于骨架和流體的導熱系數比為1000，骨架導熱起主導作用，熱量優先沿骨架傳遞，因此骨架周邊的溫度明顯高于四周，等溫線在上部區域形成明顯的拱形；但當H/D減小到0.5時，背風面的低速區域明顯減小，此時，強迫對流占主導，拱形區域逐漸消失；但在單胞尾部喉道位置，H/D=2.0時的骨架周邊溫度卻高于H/D=0.5時的溫度。為了進一步定量分析這種現象，引入場協同理論，對于對流傳熱問題，場協同角 β的表達式為：

根據場協同理論， β是流動方向與溫度梯度方向的夾角， β越小傳熱性能越優。H/D=0.5和2.0時，z= 0.5D截面上Kelvin泡沫區域的場協同角分布如圖14所示。重點關注通道上部區域，單喉道前部位置（圖14中黑色虛線標記），H/D=0.5 時 的 β明顯小于H/D=2.0 的 β。與之相反的是，單喉道尾部位置（圖14中紅色虛線標記），H/D=0.5 時 的 β 明顯大于H/D=2.0的 β，這就可以解釋為何Kelvin泡沫前部和尾部喉道位置，H/D=0.5和2.0時的溫度分布呈現明顯差異。

圖14 Re = 100時，z = 0.5D截面上Kelvin泡沫區域的場協同角分布：（a）H/D = 0.5；（b）H/D = 2.0Fig.14 The distribution of field-synergy angles on the z = 0.5D surface of the Kelvin structure for (a) H/D = 0.5 and (b) H/D = 2.0 with Re = 100

最后，高溫區域與主流方向的夾角呈現明顯差異。為了便于對比，取通道尾部溫度為0.3向通道進口連線（圖13（a）中黑色虛線），對通道的高溫區域和低溫區域進行近似劃分。取分界線與主流方向夾角為 θT?，如圖13（a）所示。從圖中可以看出， θT?隨著H/D的增大而增大；與之相反的是，高溫區域面積占比隨著H/D的增大而明顯減小。

為了進一步定量分析，沿流道中心（y=0.5H，z=0.5D），取三種條件下（H/D= 0.5、1.0和2.0）的無量綱溫度沿x方向的分布如圖13（b）所示。從圖中可以看出，三種條件下，溫度沿x方向的變化趨勢相類似，呈明顯的波動特性。在每個單喉道內，三種條件下溫度曲線發生波動的位置保持一致，分別位于喉道進口、中間位置和喉道尾部位置。此外，除了進口位置，H/D=0.5 時 溫度最高，隨著H/D的增大，溫度越來越低，這與圖13（a）中高溫面積增大的原因相同，主要是喉道高度變低，流動由骨架導熱主導變為對流傳熱占主導，使得對流傳熱性能提高，中心位置的溫度升高。對流傳熱性能的提高也會導致整體喉道前部溫度升高，這也是 θT?降低的原因，因此隨著H/D的增大，Kelvin泡沫的傳熱性能逐漸降低。為了進一步驗證，根據公式（27），計算不同H/D條 件下的Nuv如圖15所示。

圖15 Re = 10、100、1 000時Nuv和αsf隨H/D的變化Fig.15 The change of Nuv and αsf with H/D at Re = 10, 100, 1 000

從圖15中可以看出，在三種Re條件下，Kelvin泡沫的Nuv均隨著H/D的 增大明顯降低，并且，Nuv隨H/D的 變化曲線是非線性的。與此同時，隨著H/D的增大，不同Re之間Nuv的差異越來越大。H/D從0.5增大到2.0，Re=10 、 100和1000時，Nuv分別降低77%、80%和84%，說明Re越大，H/D對Nuv的影響越大。此外，圖15給出了不同H/D條件下的比表面積。從圖中可知，隨著H/D的增大，比表面積減小，這與Nuv隨H/D的變化規律保持一致。對于多孔介質材料，比表面積越大，流-固界面換熱效率越高，這也可進一步解釋Nuv隨著H/D的增大而降低的原因。

通過前述分析可知，對于不同H/D的Kelvin泡沫而言，降低H/D能很大程度上提高傳熱效率，但與此同時會使流動阻力增大。因此，需要采用綜合評價指標對該問題作進一步分析。 本文采用常用的j/f1/3（見公式（26））來評價流動傳熱性能，不同H/D條件下的j/f1/3如圖16所示。從圖中可以看出，相同H/D條件下，隨著Re的增大，j/f1/3減小，這與Sun等[19,26]以及Yang等[9]的結果一致。此外，在不同Re條件下，隨著H/D的 增大，j/f1/3也減小，與圖15中Nuv的變化趨勢相同。結果表明，盡管H/D的增大使得流動阻力減小，但Nuv的降低占主導地位，使得綜合評價指標j/f1/3減小，因此整體上看，降低Kelvin泡沫胞元高徑比有利于結構的綜合傳熱性能。

圖16 Re = 10、100、1 000時，j/f 1/3隨H/D的變化Fig.16 The change of j/f 1/3 with H/D at Re = 10, 100 , 1 000

4 結 論

本文以各向異性Kelvin泡沫為研究對象，分析了其胞元高徑比對流動傳熱特性的影響，主要結論如下：

1）不同Re下，各向異性Kelvin泡沫流動阻力（D2/Kf）及努塞爾數（Nuv）均隨著胞元高徑比（H/D）的減小而增大。當Re= 10、100和1000時，D2/Kf分別增大203%、208%和305%，Nuv分別增大77%、80%、84%。盡管不同Re下流動阻力增大的百分比均大于Nuv增大的百分比，但綜合傳熱因子j/f1/3仍隨著H/D的減小而增大，即小胞元高徑比的Kelvin泡沫具有更好的綜合傳熱性能。

2）不同H/D下 ，綜合傳熱因子j/f1/3隨著Re的增大而減小，表明Kelvin泡沫在低雷諾數下的綜合傳熱性能更優。

3）本研究在考慮泡沫骨架導熱且保證雷諾數一定的條件下獲得的規律與Sun等[19]不考慮骨架導熱且保持速度一定時得到的結果完全相反，表明骨架導熱和Re是該類問題中不能忽略的重要影響因素。

本研究假設飛行器在無重力環境下服役，未考慮重力的影響。實際上，在一些情況下，如在飛行器起落和在微重力環境下飛行時，重力的方向和大小對流動傳熱的影響不可忽略，這也是作者下一步的研究工作。