翼面熱環境的并行迭代耦合方法及熱模態分析*

黃 杰， 姚衛星

(1.南京航空航天大學飛行器先進設計技術國防重點學科實驗室 南京，210016) (2.南京航空航天大學機械結構力學及控制國家重點實驗室 南京，210016)

引 言

隨著高超聲速飛行器的發展，飛行速度越來越快，飛行器氣動加熱問題日趨嚴重[1-3]，造成了飛行器結構溫度急劇升高。高超聲速飛行器翼面在熱環境下其剛度會發生變化，進而導致結構模態發生改變。若能準確計算翼面熱環境，并且分析熱環境下翼面的熱剛度，將會對高超聲速翼面熱結構設計產生重要意義。

飛行器高速飛行過程中，氣動熱會造成翼面結構溫度急劇升高，而翼面結構溫度升高后，邊界層內氣體與壁面之間的溫度梯度將減小，導致壁面熱流密度降低，即氣動加熱與結構傳熱之間存在強烈的耦合效應。早期傳統的翼面熱環境分析方法未考慮壁面溫度對熱流密度的影響[4]，直接將分析獲得的熱流密度作為邊界條件進行結構熱分析，計算結構溫度場，再根據溫度場評價結構的熱剛度。這種未考慮翼面溫度對氣動熱影響的方法必將造成翼面溫度場和熱模態分析的不準確。

現代計算流體力學(computational fluid dynamics, 簡稱CFD)和數值傳熱學(numerical heat transfer, 簡稱NHT)的發展為高超聲速翼面熱環境問題的精確耦合分析提供了可能。一些學者通過求解Navier-Stokes方程[5-6]分析翼面氣動熱，借助CFD-FASTRAN或ANSYS-FLUENT等CFD軟件或自編程序精確計算翼面熱流密度，并通過數據插值方法完成壁面熱流密度和壁面溫度的相互傳遞，以實現氣動熱與結構傳熱之間的耦合。

結構的熱剛度和熱模態研究一直是高速飛行器結構設計的重點。國內外一些學者對翼面熱模態已經開展了一些研究[7-9]，其主要研究成果為熱環境會導致結構固有頻率的降低。以往的研究中結構所處的熱環境大多為均勻溫度場，不能反映翼面結構的真實熱環境和熱剛度，且缺乏熱應力造成的附加幾何剛度和材料剛度的變化對結構熱剛度及熱模態貢獻的研究。

筆者探討了氣動熱和結構傳熱耦合模型，提出了一種高超聲速翼面熱環境分析的并行迭代耦合方法，并進行了圓管驗證試驗算例分析，最后分析了翼面的熱環境、熱剛度以及熱模態。

1 氣動熱數值算法

在不考慮體積力和內熱源的情況下，直角坐標系下的流體動力學N-S控制方程的積分形式為

(1)

其中：W為守恒向量；Fc為對流通量；Fv為黏性通量；dS為控制體邊界面；n為邊界面的單位外法線向量。

將式(1)按有限體積法進行空間離散可得

(2)

其中：Wi和Vi分別為控制體i的守恒向量和體積；NF為控制體邊界面的數目；ΔSN為第N個邊界面的面積。

由于對流通量Fc具有高度非線性特點，且集中體現了流場的對流特征，采用具有總變差衰減(total variation diminishing, 簡稱TVD)性質的無波動、無自由參數的耗散差分格式(non-oscillatory and non-free-parameter dissipation, 簡稱NND)[10]對其進行空間離散。半離散化的上風型NND格式為

(3)

(4)

(5)

為獲得單調解，采用完全迎風的二階MUSCL格式[11]離散分裂后的無黏通量,并采用minmod限制器使空間離散格式達到空間二階精度。流體控制方程中的黏性項采用中心格式進行空間離散，此外湍流模型采用兩方程Menter′s SST模型[12]。針對非定常問題的時間離散，在n+1時刻采用時間二階精度的隱式三點向后差分,得到二階精度的離散方程為

(6)

引入虛擬時間項進行內迭代求解，并采用一階前差處理得到

(7)

其中：Δτ和Δt分別為虛擬時間步長和物理時間步長，稱為雙時間步長法[13]；p和n分別為虛擬時間迭代步和物理時間迭代步。

虛擬時間步上的內迭代可采用LU-SGS格式[14]求解，當p→∞時虛擬時間項趨近于零，式(7)的定常解即為二階精度的非定常解。

采用隱式殘差光順技術[15]加速收斂，并且針對高超聲速氣動熱問題，流體導熱系數和黏性系數通常采用Sutherland公式或分子動力學計算，其對物面熱流密度的計算精度有重要影響。

2 熱分析數值算法

在無體積熱源的假設下，結構瞬態熱傳導的控制方程為

(8)

其中：ρ0為結構材料密度；c為材料比熱容；kx，ky和kz分別為材料3個方向的導熱系數，比熱容和導熱系數一般為溫度的函數。

針對本研究熱防護系統的熱分析問題，其外表面邊界條件為壁面熱流密度Qaero和壁面熱輻射量Qrad，其表達式分別為

其中：?T/?n為壁面法向溫度梯度；r為壁面熱輻射率；玻爾茲曼常數σ=5.67×10-8W/(m2·K4)；Twall為壁面溫度；Tat為大氣環境溫度。

對式(8)進行有限元離散，得到總體合成矩陣求解方程為

(11)

針對n～n+1時間步，用Gaierkin格式離散得到

(12)

求解式(12)即可得到結構各個時刻的溫度場。

3 熱剛度與熱模態

(13)

其中：B為單元幾何矩陣；D為熱環境下單元彈性矩陣，其與結構彈性模量E及泊松比μ有關。

當結構溫度改變時，彈性模量E的改變會引起D的變化。

此外，氣動加熱下結構內部溫度分布通常是非均勻的，即結構存在溫度梯度，這將會使結構內部產生熱應力。熱應力會導致結構產生附加的初始應力剛度矩陣(幾何剛度)，從而改變結構整體剛度。

結構在溫度載荷作用下產生初始應變ε0，在熱應力σ作用下產生彈性應變D-1σ，結構的總應變為兩者之和

ε=ε0+D-1σ

(14)

熱應力σ可由式(14)得到

σ=D(ε-ε0)

(15)

在熱應力作用下結構單位體積內的應變能密度U可表示為

(16)

由式(16)可得到有限單元內熱應力產生的應變能Ue為

(17)

將式(15)代入式(17)，并考慮彈性矩陣D的對稱性得

(18)

應變ε和幾何矩陣B的關系為

ε=Bδe

(19)

其中：δe為單元節點位移向量。

將式(19)代入式(18)右端的第1項，得到

(22)

其中：Nx，Ny和Nxy分別為在各方向由熱載荷引起的薄膜應力。

由以上分析可知，熱環境下結構熱模態可通過求解以下廣義特征值問題獲得。

[(KT+KS)-ω2M]φ=0

(23)

其中：M為結構質量矩陣；KT和KS分別為結構總體材料剛度矩陣和總體附加幾何剛度矩陣。

4 熱環境分析的并行迭代耦合方法

高超聲速翼面的氣動加熱效應會造成結構溫度急劇升高，而結構溫度升高后，邊界層內氣體與壁面之間的溫度梯度將會減小，造成壁面熱流密度的降低，即翼面氣動熱與結構傳熱之間存在強烈的耦合效應，如圖1所示。

筆者采用有限體積法離散流場并求解翼面氣動熱，而結構熱傳導采用有限元法離散和求解，并采用如圖2所示的并行迭代耦合方法進行翼面熱環境分析。其特點為：a.在任意迭代步內FVM(或FEM)求解過程中壁面溫度(或壁面熱流密度)不變，即準靜態假設；b.流場和結構傳熱均在每個迭代分析步結束后進行數據交換，以保證耦合分析的協調性和時間精度;c.流場采用子循環分析，且子循環迭代次數為n。

圖1 氣動熱與結構傳熱耦合模型Fig.1 The coupled model for aerodynamic heating and structural heat transfer

圖2 并行迭代耦合方法Fig.2 Parallel iterative coupled method

氣動熱和結構傳熱的并行迭代耦合分析流程如圖3所示，其主要步驟為：

1) 對流場分析和結構傳熱分析進行內存分配，建立相應的數值分析模型，定義來流條件和初始溫度T0，進行定常流場的計算，并將初始熱流密度傳遞給結構傳熱模型；

2) 進行N=i～N=i+1步的求解，流場和結構傳熱分析同時進行，分析結束后輸出熱流密度和結構溫度場；

3) 判斷熱流密度或結構溫度場是否收斂，若收斂，結束耦合分析；若未收斂，則進行熱流密度和翼面表面溫度的數據傳遞；

4) 返回步驟2，進行下一個迭代步的求解，直到熱流密度或結構溫度場收斂，結束耦合分析。

圖3 并行迭代耦合分析流程Fig.3 Parallel iterative coupled analysis process

5 基于虛擬空間的數據插值算法

并行迭代耦合分析方法需要進行熱流密度和壁面溫度的插值計算，這是由于熱分析模型網格尺寸通常遠大于流體網格尺寸，即耦合界面節點不一致，如圖4所示，故需要在耦合面進行插值以完成數據傳遞。

圖4 FVM節點與FEM節點的對應關系Fig.4 Nodes relationship between FVM and FEM

本研究分析涉及到熱流密度和結構壁面溫度的數據插值，其中的核心技術是要保證耦合面上熱流密度通量的守恒性。

(24)

其中：Q為FVM網格面的熱流密度；q為與FVM網格節點對應的FEM網格面的熱流密度。

筆者采用虛擬空間插值方法進行熱流密度的數據傳遞，如圖5所示，其主要步驟為：

1) 將各學科耦合面上節點從物理空間(x,y,z)通過坐標變換映射到二維虛擬空間(u,v)，x=x(u,v)，y=y(u,v)和z=z(u,v)，即將三維曲面上的空間節點轉換到二維虛擬空間平面上；

2) 在物理空間中搜索任意FEM網格節點ζi(x,y,z)附近的FVM網格節點ηi(x,y,z)，并將其轉換到虛擬空間得到ζi(u,v)和ηi(u,v)；

3) 將FVM網格節點坐標ηi(u,v)和相應的熱流Qi(u,v)帶入到如下的二次插值函數Q(u,v)，采用最小二乘法求解插值函數的系數ai(i=1,2,…,10)；

Q(u,v)=a1u3+a2v3+a3u2v+a4uv2+a5u2+

a6v2+a7uv+a8u+a9v+a10

(25)

4) 將所有FEM網格節點ζi(u,v)帶入到已知系數ai的插值函數Q(u,v)中，即可求得FEM網格節點插值熱流密度qi(u,v)，并通過式(24)進行熱流密度通量的守恒性檢驗。

圖5 數據插值分析流程Fig.5 The analysis process for data interpolation

6 圓管驗證算例

筆者采用NASA的圓管風洞試驗模型[16]進行以上并行迭代耦合方法的驗證。其中：不銹鋼圓管內徑R1=25.4 mm;外徑R2=38.1mm;密度ρ=8 030 kg/m3;導熱系數k=16.72 W/(m·K);比熱容c=502.48 J/(kg·K);圓管初始溫度T0=294.4K;來流馬赫數Ma=6.47;來流溫度T=241.5 K;來流壓強P=648.1 Pa;攻角α=0°。建立了二維分析模型，如圖6所示，CFD模型壁面第1層網格高度Δh=1×10-5m，圓管結構采用四節點平面單元模擬。采用虛擬空間插值方法進行熱流密度和壁面溫度的數據的傳遞，分析類型為瞬態分析，耦合分析時間步長Δt=1×10-4s，分析總時間ttotal=2 s。

圖6 流場和結構模型網格Fig.6 The meshes for fluid and structural models

圖7為2 s時刻流場和圓管結構的溫度云圖。圖8和圖9分別為2 s時刻圓管外壁面的相對熱流密度和相對溫度(量綱為1)的分析情況。從圖中可以觀察到數值分析結果與試驗結果吻合良好，駐點處的熱流密度分析值Qstag為657 kW/m2，試驗值為670 kW/m2， 兩者的相對誤差為1.94%。此外, 駐點處溫度計算值Tstag為441 K，試驗值為465 K，兩者相對誤差為5.16%。通過本算例驗證了熱環境分析的并行迭代耦合分析方法的正確性和數據插值方法的精度。

圖7 2 s時刻流體和結構溫度場Fig.7 The temperature fields for fluid and structural models at two second

圖8 2 s時刻壁面熱流密度分布Fig.8 The distribution for wall heat flux at two second

圖9 2 s時刻壁面溫度分布情況Fig.9 The distribution for wall temperature at two second

圖10 小展弦比翼面的平面和剖面Fig.10 Platform and cross-sectional views of the low aspect ratio wing

7 翼面熱環境和熱模態分析

7.1 分析模型

筆者選取小展弦比翼面為分析模型，如圖10所示。其中，來流馬赫數Ma=6，飛行高度H=60 km，攻角α=1.5°，翼面結構初始溫度為T0=300 K，劃分了六面體CFD網格，網格總量約100萬。為了獲得網格無關性的壁面熱流密度，壁面第1層網格高度小于1×10-5m。圖11為流體CFD網格和翼面FEM網格。從圖中可觀察到這兩套網格翼面節點并非一一對應，故采用虛擬空間插值方法進行熱流密度和壁面溫度的數據傳遞。

圖11 流體和結構模型網格Fig.11 The meshes for fluid and structural models

翼面前緣為碳/碳復合材料結構，其導熱系數k1=42 W/(m·K)，密度ρ1=2 000 kg/m3，彈性模量E1=95 GPa，熱膨脹系數a1=4×10-6K；其余部位為鈦合金結構，其密度為ρ2=4.4×103kg/m3，而導熱系數k2,彈性模量E2和熱膨脹系數a2隨溫度變化情況如圖12和圖13所示。翼面外表面熱輻射率ε=0.8，大氣環境溫度Tat=247.02 K，并行迭代耦合分析中流場分析的子循環步n設置為30。

圖12 鈦合金導熱系數隨溫度變化情況Fig.12 The thermal conductivity of titanium alloy varying with temperature

圖13 鈦合金彈性模量及熱膨脹系數隨溫度變化情況Fig.13 The elasticity modulus and thermal expansion coefficient of titanium alloy varying with temperature

7.2 翼面熱環境

通過氣動熱與結構傳熱的并行迭代耦合分析方法獲得了翼面結構的穩態溫度場，如圖14(a)所示。從圖中可以觀察到翼面前緣溫度最高，從翼面前緣往后緣溫度逐漸降低。傳統的非耦合方法分析得到的翼面溫度場如圖14(b)所示。

非耦合分析方法獲得的翼面最高和最低溫度分別為1 711.5 K和928.1 K，耦合分析方法獲得的翼面最高和最低溫度分別為1 369.3 K和785.0 K，如表1所示。可見，非耦合方法分析結果偏高，無法準確評價翼面的熱環境，這是因為非耦合方法未考慮翼面溫度升高引起熱流密度降低的影響。

圖14 翼面穩態溫度場Fig.14 The steady temperature field of wing表1 耦合和非耦合翼面最高和最低溫度對比

Tab.1 The maximum and minimum temperature of wingby the coupled and uncoupled methodsK

翼面溫度耦合非耦合Tmax1 369.31 711.5Tmin785.0928.1

圖15 翼面前緣點和后緣點溫度收斂歷程Fig.15 History of convergence for leading and trailing edges of wing

圖16 耦合和非耦合分析翼面熱流密度分布情況Fig.16 The heat flux distributions of wing for coupled and uncoupled analysis

圖17 翼面最高溫度隨馬赫數的變化情況Fig.17 The maximum temperature of wing varying with Mach number

圖15為并行迭代耦合分析過程中翼面平均氣動弦長(Cmac)位置前緣點與后緣點的溫度收斂情況。從圖中可觀察到前緣點溫度收斂慢于后緣點，但迭代15步后兩點均基本收斂。圖16為翼面平均氣動弦長位置耦合分析和非耦合分析的熱流密度分布情況。可知耦合分析時翼面溫度的升高導致了熱流密度的降低，且翼面迎風面熱流密度高于背風面。翼面最高溫度隨馬赫數變化情況如圖17所示。馬赫數從5增加到9，翼面最高溫度從1 116.3 K上升到2 061.7 K，且數據點的線性度很好，可根據此規律估算其他馬赫數下的翼面最高溫度。

7.3 翼面熱模態

筆者進行了300 K和高溫環境下的模態分析，其中翼根部分為固支邊界條件。圖18為300K溫度下的翼面前4階固有振型，分別為一階彎曲、一階扭轉、二階彎曲和二階扭轉，并且前4階固有頻率分別為13.44，33.58，48.54和70.58 Hz。高溫對翼面模態的影響體現在兩方面：a.高溫會使翼面產生熱應力(預應力)，預應力會造成翼面結構產生附加幾何剛度KS，存在幾何非線性；b.溫度的升高會導致結構材料彈性模量降低(圖8)，即翼面材料剛度KT會降低。

為了研究KT和KS對翼面結構剛度和模態的影響，建立了相應的有限元模型，分析結果如表2所示。從表中結果可知,熱環境下翼面各階固有頻率相對300 K時均有所降低，且在熱環境下前4階固有頻率分別降低了15.3%,13.3%,12.0%和11.7%。其中,KS對固有頻率的影響很小，而KT卻對固有頻率影響很大，即熱環境下翼面剛度的主要影響因素為高溫引起的材料剛度的降低。

圖18 300 K情況下翼面固有振型Fig.18 The natural modes of vibration under 300 K

表2 熱模態影響因素分析

Tab.2 The analysis of influence factor for thermal modeHz

項目f1f2f3f4300 K13.4433.5848.5470.58KS13.3333.2649.6571.06KT11.4829.3441.7761.98KT +KS11.3829.1042.7362.29

筆者還分析了馬赫數對翼面結構固有頻率的影響，其中定義頻率比r(量綱為1)為熱環境下的頻率fi與300K溫度下基準頻率fi,300K的比值，分析結果如圖19所示。可見， 隨著馬赫數的增加， 前4階頻率比均降低，并且1階頻率下降最快，即低階固有頻率比高階固有頻率下降得更快。

圖19 翼面頻率比隨馬赫數變化情況Fig.19 The frequency ratios of wing varying with Mach number

(26)

8 結 論

1) 提出了一種分析高超聲速翼面熱環境的并行迭代耦合分析方法，氣動熱采用CFD分析，結構傳熱采用FEM分析，兼顧計算效率和計算精度，并利用獲得的翼面溫度場進行了結構熱模態分析。

2) 進行了圓管試驗驗證算例分析，2 s時刻的駐點溫度計算值與試驗值相對誤差為5.16%，從而驗證了本研究熱環境分析并行迭代耦合方法。

3) 相對傳統的非耦合方法，迭代耦合分析方法考慮了翼面溫度升高引起熱流密度降低的因素，能準確評價翼面的熱環境，而翼面最高溫度與馬赫數具有良好的線性關系。

4) 熱環境導致了翼面結構剛度及固有頻率的降低，且材料彈性模量的降低是其主要影響因素。此外隨著馬赫數的增加，低階固有頻率比高階固有頻率下降得更快。