氣動加熱環境下壁面熱響應的快速分析方法

洪 元，劉亮堂，楊立明，祝成民

(1 航天恒星科技有限公司，北京 100086；2 北京航空航天大學宇航學院，北京 102206)

0 引言

飛行器的熱防護設計需要計算殼體的溫度分布。計算的關鍵是得到準確的輸入熱流等邊界條件信息。對氣動加熱的分析方法主要包括經驗公式、數值分析和風洞實驗。風洞實驗得到的數據最準確，但實驗費用昂貴，且受條件限制，實驗的相似準則參數在很多情況下達不到實際飛行時的參數值，只能在少數情況下使用。經驗公式方法計算簡單，但是適用范圍窄。數值分析方法可以得到詳細的熱流分布數據，在通過實驗校正的情況下也有較高的精度，是目前使用比較廣泛的氣動熱分析方法。已有的研究結果表明，氣動熱計算會受到物體表面溫度的影響，而物體表面溫度是氣動加熱和物體內部傳熱與蓄熱互相平衡的結果。所以只有將氣動加熱與被加熱物體傳熱耦合計算才能保證計算結果與實際情況相符。但是，氣動熱數值分析計算量大、耗時長，造成耦合計算的代價太大，不適合工程設計。為了解決這個問題，研究人員提出了快速計算方法，但是這些方法還沒有在飛行器熱防護設計中實際應用，其有效與否仍需通過實踐檢驗。如何在保證分析精度的基礎上加快計算速度仍然是氣動加熱熱防護設計中急需解決的問題。

文中借鑒經驗公式方法的原理，提出了一種氣動熱與壁面傳熱的松耦合計算方法，該方法應用于氣動加熱熱防護工程設計取得了較好的效果。

1 理論描述

1.1 物理模型

假設殼體為等厚度平板結構，內部熱源為恒溫熱源，其外壁面承受均勻的氣動加熱，內壁面通過對流和輻射與內部熱源進行熱交換。由于殼體上下表面傳熱均勻，殼體內相同深度處的溫度在任何時刻都相同，所以只需考慮熱量沿深度方向的傳遞，殼體熱響應可簡化為一維非穩態問題，如圖1所示。

圖1 物理模型示意圖

1.2 氣動加熱計算模型

根據邊界層傳熱有關理論，通過邊界層對壁面的對流換熱可表示為：

=(-)

(1)

式中:氣體對壁面的熱流密度;,為邊界層外緣的氣體密度和速度;為斯坦頓數;為恢復焓;為壁面處氣體的焓。令=，則式(1)可簡化為：

=(-)

(2)

由式(2)可以看出，氣動加熱中的對流換熱計算可以歸結為如何確定恢復焓和換熱系數的問題。為了解決這一問題，需要使用數值分析方法求解氣動熱問題。首先設置氣動熱計算的壁面邊界條件為絕熱壁得到恢復焓；然后設置壁面為恒溫邊界得到參考溫度下的壁面熱流，將這兩個計算得到的和代入式(2)可求得到換熱系數值。在得到恢復焓和換熱系數后，再計算壁面溫度因內部傳熱而改變的最終情況。該熱流作為壁面熱傳導方程的邊界條件給出，從而實現氣動熱與物體內部傳熱的松耦合計算。

1.3 傳熱方程

根據上述物理模型假設，物體內部熱傳導控制方程為：

(3)

式中:為時間;為深度;(,)為物體內部深度為的點在時刻的溫度;=()為熱擴散系數;為物體的導熱系數;為物體密度;為物體的比熱容。

1.4 邊界條件

邊界條件包括外壁面邊界和內壁面邊界，具體如下：

1)外壁面邊界條件

殼體外壁面承受氣動加熱的作用，考慮到輻射散熱的影響，外壁面邊界為：

(4)

式中:=，為空氣定壓比熱容;為外表面輻射系數。

2)內壁面邊界條件

殼體內壁面邊界條件分為兩種情況：一種內壁面與恒溫熱源直接接觸，則有:

=，=

(5)

另一種內壁面與恒溫熱源通過對流和輻射進行換熱，則有:

(6)

式中:為殼體厚度;為恒溫熱源;為內壁面與恒溫熱源的對流換熱系數;為內表面輻射系數;為恒溫源輻射系數。

2 數值計算

2.1 計算方法

熱傳導控制方程，即式(3)為二階拋物型偏微分方程，可以采用經典的顯式中心差分格式求解。

如圖2所示，將由空間坐標和時間坐標構成的計算域劃分為均勻網格。空間方向的網格間距為，時間方向上的網格間距為。圖中圓點代表上一時間步的已知量，叉代表待求量。

圖2 差分網格示意圖

經典顯式中心差分格式的計算公式為：

,+1=(-1,-2,++1,)+,

(7)

式中,=，稱為網格比。根據偏微分方程數值穩定性分析，當≤12時，數值格式是穩定的。

2.2 邊界條件處理

1)外壁面邊界

對于外壁面邊界條件，設：

則根據式(4)有(0,+1)=0。采用牛頓迭代法求解該方程：

(8)

式中′(·)為(·)的導數。迭代過程中選擇1,+1為初始條件。

2)內壁面邊界

對于內壁面邊界條件，設邊界節點所對應=。如果內壁面直接與恒溫熱源接觸，則有：

,+1=

(9)

如果內表面通過對流和輻射與恒溫熱源進行熱交換，設：

則根據式(6)有(,+1)=0。采用牛頓迭代法求解該方程：

(10)

式中′(·)為(·)的導數。

3 計算結果及有效性分析

為驗證文中分析方法的有效性，以型號1、型號2兩種導彈的天線窗為例，選取其經過飛行驗證的樣本數據與文中計算結果進行對比分析。兩導彈的天線窗均采用石英聚酰亞胺材料，其物性參數如表1所示。

表1 天線窗材料物性參數

3.1 案例1

型號1導彈天線窗厚度為5 mm，包含37 s、140 s、240 s三種飛行工況：

1)飛行工況1

飛行時間37 s，計算結果如圖3所示。從圖中可以看出，整個氣動加熱過程存在兩個峰值，本模型計算結果與樣本數據基本一致，二者最大偏差為24.5 ℃，其相對誤差約為4.9%。

圖3 型號1天線窗壁溫計算結果對比圖(飛行工況1)

2)飛行工況2

飛行時間140 s，計算結果如圖4所示。從圖中可以看出，整個氣動加熱過程存在3個峰值，本模型計算結果與樣本數據基本一致，二者最大偏差為5 ℃，其相對誤差約為2%。

圖4 型號1天線窗壁溫計算結果對比圖(飛行工況2)

3)飛行工況3

飛行時間240 s，計算結果如圖5所示。從圖中可以看出，整個氣動加熱過程存在兩個峰值，本模型計算結果與樣本數據基本一致，二者最大偏差為7 ℃，其相對誤差約為3.3%。

圖5 型號1天線窗壁溫計算結果對比圖(飛行工況3)

3.2 案例2

導彈天線窗厚度為7 mm，其飛行時間為418 s，計算結果如圖6所示。

從圖中可以看出，整個氣動加熱過程存在3個峰值，本模型計算結果與樣本數據基本一致，二者最大偏差為24 ℃，其相對誤差約為5%。

圖6 型號2天線窗壁溫計算結果對比圖

3.3 結果分析

通過對比可以看出，文中模型的計算結果與樣本數據整體符合性較好，計算誤差可控制在5%以內，滿足工程設計的需求，證明了文中模型的正確性及工程應用的可行性。

此外，文中模型采用氣動熱與物體熱傳導的松耦合計算方式，其計算資源和計算時間主要集中在恢復焓和恒壁溫熱流的計算；傳熱計算基于一維傳熱模型完成，資源占比較小，計算時間一般不超過1 min。在氣動熱防護設計時，由于氣動外形變化很小，工程上可假定恢復焓和恒壁溫熱流等氣動熱環境數據一定。在這種情況下，每次迭代設計采用文中模型可在幾分鐘內完成計算。相比于氣動熱與熱傳導的緊耦合計算模型，文中模型可將計算效率提高一個量級，非常適合氣動熱防護設計的工程應用。

4 結論

文中建立的松耦合計算方法可以正確評估氣動加熱對飛行器殼體等相關部件的熱響應問題；使用恢復焓和傳熱系數解耦氣動熱與物體熱傳導計算的方法具有極高的計算效率，相比于傳統緊耦合計算方法效率提高一個量級，滿足氣動熱防護工程設計需求，并為高超聲速飛行器的氣動熱防護優化設計提供技術支撐。