接觸熱阻對超聲速飛行器結構響應分析的影響

張 仡,成 竹,秦 強,陳 宏

(中國飛機強度研究所全尺寸飛機結構靜力/疲勞航空科技重點實驗室，西安 710065)

近年來，超聲速飛行器日益受到各國的廣泛關注[1]。與亞聲速飛行器不同，超聲速飛行器表面的氣動加熱現象更為顯著，對飛行器結構溫度與應力的影響更為突出。在氣動加熱過程中，熱量通過蒙皮傳遞到內部梁、肋等結構中，界面間的接觸熱阻會直接影響結構的溫度與應力。進行分析時，對接觸熱阻的影響估計過高會使結構過于笨重，而估計過低會帶來安全隱患。因此，必須準確判斷接觸熱阻對于結構響應的影響程度。

中外學者在接觸熱阻理論、仿真與試驗方面已經進行了大量研究。在理論方面，通過描述接觸面表面形貌與變形，并參照試驗數據，建立了一系列預測接觸熱阻的理論模型[2-4]。Fletcher[5]、Sridhar等[6]和Lambert等[7]分別綜述了部分理論模型。對接觸熱阻的理論研究拓展到微納尺度，提出了基于熱載子的理論模型[8]。在仿真與試驗方面，Mantelli等[9-10]、Song等[11]通過仿真與試驗研究了螺栓連接件中的接觸熱阻；Tirovic等[12]通過仿真與試驗研究了汽車剎車片界面的接觸熱阻；Milanez等[13]對低壓力下金屬間接觸熱阻進行了試驗研究并提出了相應的理論模型；Liu等[14-15]對C/C復合材料與GH600高溫合金之間的接觸熱阻進行了試驗研究，并通過數值仿真研究了接觸熱阻對疏導式熱防護結構防熱效果的影響。目前，對于接觸熱阻的研究主要集中在電子、衛星、核能等領域，而對于超聲速飛行中的接觸熱阻問題關注較少[13]。但考慮到熱問題在超聲速飛行領域中的重要地位，有必要開展接觸熱阻對超聲速飛行器結構響應分析影響的研究。

為此，以超聲速飛行器外表面沉頭螺栓局部模型為例，研究在數值仿真中，接觸熱阻對超聲速氣動載荷下結構響應的影響，為超聲速飛行器設計及安全評估提供參考。

1 接觸熱阻理論

在Cooper等[2]提出的(Copper-Mikic-Yovanovich)CMY接觸熱阻模型中，假設金屬表面粗糙峰的高度服從高斯分布，粗糙峰在金屬表面上隨機分布，粗糙峰變形均為塑性變形。基于以上假設，CMY模型中金屬界面間接觸熱導表示為

(1)

(2)

ψ(ε)=(1-ε)1.5, 0<ε<0.3

(3)

(4)

式中：hc為接觸熱導；n為接觸點密度；a為接觸點平均半徑；ks為等效熱導率；ψ為熱流約束參數；ε為接觸點相對尺寸；k1與k2分別為界面兩側材料的熱導率；Ar為接觸面實際接觸面積；Aa為接觸面名義接觸面積。

定義相對平均板間距λ=Y/σ，其中Y為實際平均板間距。對于塑性變形，Ar/Aa、n、a可用λ表示為

(5)

式(5)中：σ為等效均方根(RMS)表面粗糙度；m為等效粗糙峰斜率平均值。σ與m由界面兩側的表面形貌參數表示為

(6)

(7)

式中：σ1與σ2分別為界面兩側表面RMS粗糙度；m1與m2分別為界面兩側表面粗糙峰斜率平均值。

通過對金屬表面形貌的試驗研究，Lambert等[16]提出的m與σ的試驗關系式為

(8)

由式(1)～式(5)得到無量綱形式的接觸熱導Cc為

(9)

式(9)中,無量綱接觸熱導只與相對平均板間距λ有關。在CMY模型中，λ只能通過試驗獲得。為解決這一問題，Yovanovich[3]提出λ與P/Hc(其中，P為接觸面壓力，Hc為材料微硬度)的關系為

(10)

(11)

(12)

(13)

(14)

2 考慮接觸熱阻的局部建模

在熱力耦合接觸問題中，接觸熱阻使得結構應力場與溫度場之間存在復雜的耦合關系。根據式(9)、式(10)，面間壓力對接觸熱阻有直接影響，而接觸熱阻直接影響溫度場，從而使應力場間接對溫度場造成影響；而溫度場又直接影響應力場。為了避免間接耦合引起的計算誤差，采用直接耦合法進行熱力耦合分析。

為在計算中引入接觸熱阻，需在接觸關系中，輸入一系列不同材料、壓力與溫度下的hc數值。計算過程中，通過線性插值獲得局部區域在給定材料、壓力與溫度下的hc，從而將接觸熱阻引入傳熱方程中，從而進行考慮接觸熱阻的熱力直接耦合計算。

問題的計算流程圖如圖1所示。首先進行常溫下應力分析，在接觸面間形成初始接觸狀態，并通過插值獲得常溫下接觸熱阻。將常溫下的應力場作為初始條件，開始直接熱力耦合計算，對問題進行迭代求解。每次計算中，需判斷計算前后的接觸狀態是否變化以及應力場與溫度場是否收斂。如不收斂，則需根據計算得到的接觸狀態、應力場與溫度場，重新計算材料的力熱參數及界面間的接觸熱阻，進行新的直接熱力耦合計算。如此反復迭代直至收斂，從而獲得考慮接觸熱阻的穩態結構響應。

2.1 幾何模型

選取超聲速飛行器前緣中，以螺栓為中心，邊長30 mm的方形局部模型進行研究。考慮到幾何與載荷的對稱性，取該方形模型的一半進行研究，如圖2所示。為降低計算量，不考慮沉頭螺栓與螺母間的接觸。

圖1 考慮接觸的熱力耦合計算流程圖Fig.1 Flow chart of thermo-mechanical coupling calculation considering contact

圖2 局部模型網格Fig.2 Grids of local model

圖2所示的模型中，蒙皮與肋的厚度為3 mm，梁的厚度為2 mm。沉頭螺栓公稱長度l=8 mm，螺紋規格為d=M5，其余尺寸參照《MJ螺紋十字槽100°沉頭螺栓》(GJB 3374/29)。螺母螺紋規格為D=M5，其余尺寸參照《抗剪型MJ螺紋高鎖螺母》(GJB 3377—1998)。

采用考慮傳熱的六面體減縮積分實體單元進行有限元仿真。由于螺栓孔邊、螺栓沉頭處與倒角處可能出現應力集中，對這些區域的網格進行了加密處理。

2.2 材料及其物理屬性

模型中，蒙皮、梁、肋的材料為TC4鈦合金，沉頭螺栓與螺母的材料為GH2132高溫合金。兩種材料的熱導率、線膨脹系數、彈性模量與泊松比如表1～表3所示。TC4鈦合金退火后的硬度為HBS 255～341，固溶時效硬度為HBS 293～361。GH 2132高溫合金硬度統計值為HBS 293～306。計算時取TC4與GH 2132硬度均為HBS 300。

表1 不同溫度下的材料熱導率

表2 不同溫度下的材料線膨脹系數

表3 不同溫度下的材料彈性模量及泊松比

2.3 模型中的接觸關系

模型中的接觸面包括：蒙皮與梁之間的接觸面；梁與肋之間的接觸面；蒙皮、梁、肋與螺栓及螺母的接觸面。根據接觸面兩側的材料，將接觸面分為TC4與TC4之間的接觸面、TC4與GH2132之間的接觸面。

給定所有接觸面間切向摩擦系數均為0.15。

為獲得界面間無量綱形式接觸熱導，需得知界面兩側材料的布氏硬度HB，界面兩側表面的RMS粗糙度σ1與σ2，以及界面平均壓力。首先，由較軟材料的布氏硬度HB，通過式(12)、式(13)算得維氏微硬度系數c1與維氏微硬度無量綱系數c2。之后，由σ1與σ2通過式(6)～式(8)分別算得等效RMS表面粗糙度σ與等效粗糙峰斜率平均值m。根據《GH2132 MJ螺栓螺紋和螺釘通用規范》(GJB 8618—2015)，參考機械設計中對于連接部位的相關要求，給定蒙皮、梁、肋及螺栓的表面粗糙度均為σ=3.2 μm。隨后，由σ、m、c1、c2通過式(11)、式(12)算得相對接觸壓力P/Hc與相對平均板間距λ。最后，將λ代入式(9)獲得無量綱形式的接觸熱導Cc。計算得到的Cc如表4所示。引入給定溫度下的等效熱導率ks，由式(9)最終求得對應溫度下的接觸熱導hc。

表4 不同壓力下的無量綱接觸熱導

2.4 邊界條件與載荷

在超聲速飛行中，圖2所示的局部結構同時承受高溫空氣的氣動加熱以及飛行器其他部位傳遞來的剪力。因此，計算中對模型施加熱學邊界，以模擬氣動熱載荷；施加力學載荷與位移邊界，以模擬結構中的剪力。

模型中，載荷與邊界條件的位置如圖3所示。圖3中，A為模型外表面，在此處給定溫度邊界條件；B為蒙皮與肋的右側端面，在這兩個面上約束沿法向位移；C為肋內表面、肋右側端面與模型中間截面上的交點，在該點上給定固支邊界條件；D為模型內表面，在此處給定自然對流與輻射邊界條件，其中對流換熱系數為20 W/(m2·K)，環境溫度為 20 ℃，表面發射率為0.5；E為梁左側端面，在該面上定義沿法向指向面外的均布拉力；F為螺栓中垂直于軸向的截面，在該面上定義沿螺栓軸向大小為 1 932.25 N 的預緊力。

A為模型表面；B為蒙皮與肋的右側端面；C為肋內表面、肋右側端面與模型中間截面上的交點；D為模型內表面；E為梁左側端面；F為螺栓截面圖3 模型中邊界與加載位置Fig.3 Boundary and loading position in model

3 接觸熱阻影響分析

3.1 接觸熱阻對溫度場的影響

3.1.1 不考慮接觸熱阻時的溫度場

不考慮接觸熱阻時，應力對傳熱過程沒有顯著影響，結構溫度場只與熱學邊界條件有關。不同外表面溫度下的結構相對溫度分布如圖4所示。相對溫度為各點攝氏溫度與外表面攝氏溫度的比值。圖中，等溫線在所有交界面上連續，表明界面間不存在接觸熱阻。當外表面溫度升高時，結構內表面附近的相對溫度略有下降。這是由于溫升引起的結構熱阻下降不足以抵消內表面輻射換熱量上升對傳熱過程的影響。

圖4 不考慮接觸熱阻時的相對溫度分布Fig.4 Relative temperature distribution without contact thermal resistance

相同外表面溫度，拉力分別為5、50 MPa時的結構相對溫度分布如圖5所示。由圖5可知，即使拉力增大10倍，結構的相對溫度場分布仍未出現明顯差異。表明不考慮接觸熱阻時，結構應力場對溫度場沒有明顯影響。

3.1.2 考慮接觸熱阻時的溫度場

拉力為50 MPa時，不同外表面溫度下的結構相對溫度分布如圖6所示。由于接觸熱阻的影響，等溫線在交界面上不連續。由表5可知，在蒙皮與梁的交界面上，蒙皮內表面的最低溫度大于梁上表面的最高溫度,并且兩者的溫差隨著外表面溫度的升高不斷上升。梁與肋的交界面間也存在不低于 20 ℃ 的溫差。界面兩側的溫差表明接觸熱阻阻礙了蒙皮、梁、肋間的熱傳遞。同時，外表面溫度越高，界面兩側的溫差越大。即氣動加熱越強，接觸熱阻對結構傳熱的影響越顯著。

圖5 不考慮接觸熱阻時不同拉力下的相對溫度分布Fig.5 Relative temperature distribution under different tensile forces without contact thermal resistance

圖6 考慮接觸熱阻時的相對溫度分布Fig.6 Relative temperature distribution with contact thermal resistance

與不考慮接觸熱阻時的結構溫度場相比，考慮接觸熱阻時結構各部分溫度明顯下降，如表6所示，其中肋的溫度下降幅度最大。對比圖4、圖6發現越是接近內表面的位置，考慮接觸熱阻時溫度下降幅度越大。這是因為越接近內表面，傳熱路徑中的接觸面越多，接觸熱阻的對傳熱過程的影響更為突出。

表5 不同外表面溫度下交界面兩側溫度極值Table 5 Extreme temperature on both sides of interface at different outer surface temperatures

表6 外表面溫度500 ℃時結構各部分溫度極值

相同外表面溫度下，結構承受不同拉力時，肋中螺栓孔表面的溫度分布如圖7所示。由圖7可知，拉力增大10倍后，螺栓孔左側溫度略有上升，右側的溫度降低了約20℃。這是因為拉力的增大使得肋與螺栓在左側接觸壓力上升，右側接觸壓力下降，從而引起界面上接觸熱阻的變化，進而對溫度場產生影響。表明在考慮接觸熱阻時，應力場通過影響結構中的接觸關系，間接影響了溫度場分布。

圖7 考慮接觸熱阻時肋中螺栓孔表面溫度分布Fig.7 Surface temperature distribution of bolt hole in rib with contact thermal resistance

3.2 接觸熱阻對應力場的影響

考慮接觸熱阻時，接觸熱阻通過影響溫度場從而影響結構熱應力與熱變形。根據表6，接觸熱阻對螺栓、梁、肋的溫度場影響較大。因此，接觸熱阻對這些部位的應力場也有較大影響，主要表現為應力分布與幅值的變化。

3.2.1 結構應力分布

結構中，螺栓孔表面的Mises應力較大。外表面溫度為500 ℃，結構承受5 MPa拉力時，螺栓孔表面應力分布如圖8所示。由圖8可知，不考慮接觸熱阻時，除沉頭位置的出現了應力集中外，螺栓孔表面其他部位的應力梯度較小。考慮接觸熱阻時，螺栓孔表面應力梯度增大，主要表現為蒙皮、梁、肋的交界部位的局部應力變化。

考慮接觸熱阻時，螺栓孔表面的應力分布變化是由溫度場變化引起的。外表面溫度為500 ℃，結構承受5 MPa拉力時，結構的變形情況如圖9所示。由圖4、圖6、表6可知，不考慮接觸熱阻時，蒙皮、梁、肋間溫度差低于20 ℃，沿x方向變形基本一致；考慮接觸熱阻時，蒙皮與肋之間最大溫差為129 ℃，肋沿X方向變形量不足蒙皮的2/3，導致螺栓轉動角度增加，使得螺栓與蒙皮、梁、肋間的接觸位置與壓力產生變化，導致螺栓孔表面的應力梯度上升。

因此，考慮接觸熱阻時，接觸熱阻通過影響結構溫度場，從而影響熱應力與熱變形，進而影響了結構應力分布。

圖8 承受5 MPa拉力時螺栓孔表面應力分布Fig.8 Stress distribution on bolt hole surface under 5 MPa tension

圖9 放大30倍的結構變形Fig.9 30-fold enlargement of structural deformation

3.2.2 結構應力水平

選取結構中各部分積分點上的最大Mises應力值代表其應力水平。不同的外表面溫度與拉力下，各部分的Mises應力最大值如圖10所示。圖10中，含(termal contact resistance)(TCR)的圖例代表考慮接觸熱阻時的情況。

圖10 結構各部分應力水平隨拉力變化情況Fig.10 Variation of stress level of structural parts with tensile force

由圖10可知，外表面溫度較低時，接觸熱阻對各部分應力水平的影響并不明顯。外表面溫度≥300 ℃時，接觸熱阻引起的結構溫度變化較為明顯，進而引起的應力變化也較為顯著。同時，隨著拉力的增大，接觸熱阻的影響也逐漸增大。其中，肋的應力水平變化最為顯著。當外表面溫度為500 ℃，拉力45 MPa時，由接觸熱阻引起的肋應力水平變化量達到89.6 MPa，相較于不考慮接觸熱阻時增長了23.0%。螺栓、蒙皮、梁中的應力水平變化率較小，均低于10%。不同溫度下各部件應力水平的最大變化量如表7所示。

表7 考慮接觸熱阻時各部分應力水平最大變化量Table 7 Maximum variation of stress levels in different parts with thermal contact resistance

4 結論

針對超聲速飛行中的結構響應問題，形成了考慮接觸問題的熱力直接耦合有限元計算方法，建立了考慮接觸熱阻的前緣局部有限元模型，研究了接觸熱阻對超聲速飛行器結構響應分析的影響，得到以下結論。

(1)在超聲速飛行器結構響應分析中，考慮接觸熱阻對溫度場有顯著影響，主要表現為改變了結構整體溫度分布，其中遠離高溫面的部位溫度變化最大。

(2) 在超聲速飛行器結構響應分析中，考慮接觸熱阻，使得應力場通過影響接觸關系從而間接影響溫度場，主要表現為螺栓孔及其附近溫度幅值的變化。

(3) 在超聲速飛行器結構響應分析中，當外表面溫度低于300 ℃時，接觸熱阻對結構應力場影響不大；當外表面溫度大于等于300 ℃時，接觸熱阻對結構應力場有顯著影響，主要表現為改變螺栓孔表面及附近的應力分布，并引起結構各部分應力水平的變化，其中肋的應力水平變化最大，其余部位的應力水平變化較小。