高超聲速全動翼面全時域耦合分析方法及應用

沈恩楠，郭同慶，吳江鵬，胡家亮，張桂江

1. 航空工業沈陽飛機設計研究所，沈陽 110035

2. 南京航空航天大學 航空學院，南京 210016

在高超聲速飛行器飛行過程中，激波和附面層的相互作用會使飛行器周圍的空氣溫度急劇升高，形成劇烈的氣動加熱環境。在高溫條件下，材料的性能會發生變化，結構內部會產生附加的熱應力，明顯影響飛行器的結構剛度。區別于一般的飛行器，高超聲速飛行器需要綜合考慮氣動熱及氣動彈性之間的耦合問題，即熱氣動彈性問題。熱氣動彈性問題屬于多學科綜合問題，難以采用統一方程求解[1]，大多數研究都是在合理假設下，采用簡化模型或者忽略其中次要的耦合部分，采用分區-邊界耦合的求解方法。

國內外有很多關于高超聲速熱氣動彈性問題的研究。楊超等[2]對高超聲速熱氣動彈性以及氣動伺服彈性的進展進行總結，提出了高超聲速飛行器在氣動彈性領域需要解決和關注的若干問題。Klock和Cesnik[3]采用活塞理論與Eckert參考焓理論分析了導彈的熱氣動彈性特性。以一定數量的固定溫度分布狀態作為參考狀態，分析得到這些狀態的模態振型，然后以這些模態振型為基礎來描述不同溫度分布狀態的模態振型，強調了溫度對模態的影響。Mcnamara等[4-5]將Eckert參考焓理論和計算流體力學(CFD)方法計算得到的熱流作用于結構有限元模型上，將熱結構振型插值到流場中，在時域內求解結構運動方程，模態保持不變。Lv等[6]采用松耦合方法研究了在固定來流參數條件下平板高超聲速機翼熱氣動彈性響應，計算得到不同時刻的結構溫度場與相應熱模態。Ye等[7]在給定的飛行狀態下，采用CFD與傳熱耦合計算程序得到結構溫度，然后通過模態分析獲得結構模態。采用松耦合方法確定高超聲速熱顫振臨界速度。郭同慶等[8]針對沿軌道運動的高超聲速翼面建立了同步計算方法，并對軌道狀態點的顫振特性進行分析，給出了狀態點的顫振邊界。張偉偉等[9]通過熱流與結構傳熱耦合計算獲得結構溫度分布與熱模態，并應用于熱氣動彈性仿真計算。吳志剛等[10]對高超聲速翼面結構、氣動力和氣動熱耦合動力學問題進行了研究。分析了結構在熱環境下的固有動力學特性。應用活塞理論計算高超聲速非定常氣動力，指出熱效應對小展弦比根部固支翼面的扭轉剛度影響很大，從而導致彎扭耦合形式的顫振臨界速度大幅下降。史曉鳴和楊炳淵[11]分析了瞬態加熱狀態下結構溫度場及固有振動特性。研究了加熱對結構固有振動特性的影響。葉獻輝和楊翊仁[12]根據Von Karman大變形應變位移關系與伽遼金法建立起熱環境下的壁板顫振計算方程。分析了不同溫度分布的壁板顫振特性。吳振強等[13]使用NASTRAN軟件建立分析模型，給出熱環境下壁板結構線性和非線性振動特性分析流程。分析了均勻和非均勻溫度分布狀態的壁板振動特性。王宏宏等[14]研究了變厚度導彈翼面模型在熱效應作用下的結構振動特性、結構熱應力以及材料彈性模量的變化影響了結構的固有頻率。認為材料參數的變化比熱應力的影響效果更明顯。楊享文等[15]采用CFD方法計算定常氣動熱，得到結構瞬態溫度場后分析模態。采用基于CFD技術的當地活塞理論作為氣動彈性分析的非定常氣動力模型。結果發現，結構在受到氣動加熱后，其固有頻率先快速減小，然后保持不變。固有頻率的改變對結構的顫振特性變化起決定性作用，顫振臨界速度變化規律與各階頻率之間的間距變化規律相同，均先快速減小后略有增大。劉成等[16]構造了不同轉捩位置的熱分布模型，研究了轉捩與舵面轉軸相對位置對高超聲速顫振特性的影響。采用Euler/Navier-Stokes(N-S)方程計算定常流場，積分得到當地活塞理論相關參數以及氣動力表達式，通過近似擬合得到結構溫度場。譚光輝等[17]提出了分步熱顫振工程分析方法，考慮翼型厚度的影響，研究了溫度對結構各階固有頻率的影響規律。李麗麗[18]建立并驗證了不同邊界條件下的壁板熱模態分析方法，得出溫度對固有頻率的影響趨勢。在活塞理論基礎上建立了翼面顫振方程的一般表達式。桂業偉等[19]研究了氣動力/熱與結構熱響應多場耦合問題的求解流程，提出了針對該耦合問題特有的時間-空間耦合概念。

氣動加熱形成的熱環境與構成氣動彈性問題的氣動力、慣性力和彈性力形成的耦合關系強弱程度不同[2]，大多數熱氣動彈性問題的研究方法都是忽略了結構振動對溫度場的直接和間接影響，將高超聲速飛行器飛行過程分解為若干狀態點，在每個狀態點上固定來流參數計算結構振動響應。

隨著技術的發展，先進的高超聲速飛行器具有變軌飛行能力，能夠在沿軌道飛行的過程中控制翼面偏轉來改變姿態和飛行軌道。操縱面作為最可能出現氣動彈性失穩的部件之一，將會成為高超聲速飛行器熱氣動彈性分析的重點對象[20]。應當建立能夠模擬高超聲速飛行器機動過程中的操縱面氣動力、氣動熱、熱傳導和結構振動的相互作用的方法。

本文針對沿軌道運動的高超聲速全動翼面的偏轉過程，在流場-結構溫度場同步計算方法的基礎上，建立了氣動力/氣動熱/結構溫度場/結構振動多物理場全時域耦合的計算分析方法。通過改變遠場來流密度、壓強和溫度來模擬飛行軌道高度變化，采用移動坐標系和動網格相結合的方式分別描述變速度飛行和翼面偏轉過程。通過坐標系變換將翼面偏轉和振動的網格變形量疊加，考慮翼面振動和偏轉過程的耦合非定常效應。通過結構振動時域響應分析研究了結構振動響應對熱流密度和溫度分布的影響。

1 多物理場全時域耦合計算流程

熱氣動彈性問題涉及多物理場之間的耦合作用，難以采用統一方程求解，在工程上一般采用分區求解-邊界耦合的方法，包括氣動力/氣動熱計算、熱傳導計算以及熱結構振動響應分析等。各物理場耦合關系緊密，如圖1所示。

圖1 熱氣動彈性多物理場耦合關系Fig.1 Aerothermoelasticity multiphysics coupling relationship

作用在結構表面的熱流會使結構升溫，高溫會影響結構的材料屬性，同時使結構產生附加的熱應力，改變結構剛度特性，進而影響結構的振動響應特性。結構表面的壓強和熱流分布對結構振動響應比較敏感，結構振動特性的改變會通過改變熱流分布間接影響結構溫度場。本文建立的全時域耦合計算流程如圖2所示。詳細流程如下：

圖2 全時域多場耦合方法計算流程Fig.2 Process of transient multi-field coupling computation

① 讀入結構初始外形(初始結構變形def0一般為0)，以及初始溫度T0。流場與結構溫度場同步推進求解，得到下一時刻的流場與結構溫度場。

② 將第①步計算得到的結構表面壓強作為載荷作用在有限元模型表面，同時將計算得到的結構溫度場分布T1賦值給有限元模型。

③ 考慮結構瞬時溫度場的非線性影響，計算結構瞬態響應，得到t1時刻結構響應def1。

④ 將該時刻結構振動產生變形傳遞給流場物面邊界，流場快速生成動態網格。

⑤ 重復①-④步，各物理場沿時間推進計算。

2 各物理場計算與耦合方法

2.1 氣動熱-結構溫度場同步計算

采用基于CFD的同步計算方法[21]求解非定常N-S方程和結構熱傳導方程，實現圖2中的第①步。該方法計算的流場和結構溫度場模型是有限體積模型。高超聲速流場對流通量采用AUSM+(Advection Upstream Splitting Method+)格式進行離散，黏性通量項采用中心格式離散。

對于結構熱傳導，由于固體介質中的熱傳導過程不存在對流現象，顯然直角坐標系3個速度方向分量為0，即u=v=w=0，連續方程與動量方程自然滿足。固體域能量方程為

(1)

式中：T為溫度；t為時間；k為第1層單元熱傳導系數；ρ為密度；C為比熱容；x、y、z為直角坐標系3個方向分量；下標s表示固體介質。

其本質與黏性擴散方程相同，也采用中心格式離散。整場采用雙時間推進法，每個物理時間步內虛擬定常迭代過程采用隱式LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[21]。

在進行同步計算之前，要對流場和結構溫度場進行初始化。結構傳熱時間遠慢于流場建立的時間，通常在結構剛開始感受到流場傳熱之前，流場已經建立，可以將等溫壁邊界定常計算得到的收斂解作為初始時刻的流場初值，將等溫壁邊界的溫度值作為結構溫度場初值。在同步求解過程中，流-固交界面的溫度Tw可以采用熱流密度連續條件確定：

(2)

式中：df和ds分別為氣體和固體介質第1層單元中心到交界面之間的法向距離；Tf和Ts分別為氣體和固體介質第1層單元溫度；kf和ks分別為氣體和固體介質第1層單元熱傳導系數，均為溫度的函數。所用的流場網格與結構溫度場網格在交界面處一一對接，無需插值。

2.2 氣動力-結構瞬態響應耦合計算

得到某時刻的流場和結構溫度場后，通過2個模型之間的映射關系實現第②步的結構溫度場賦值和氣動力插值，將有限體積模型的溫度場和物面邊界壓強傳遞給結構有限元模型。

用于計算結構溫度場的有限體積模型的網格密度與計算響應的有限元模型的網格密度差別較大，溫度場和流場計算模型的網格要遠密于有限元模型的單元。采用最鄰近搜索法建立有限體積模型到有限元模型的映射關系。在每一個物理時間步，對流場物面壓強和插值得到的結構表面壓強進行積分：

(3)

式中：pf為計算得到流場物面壓強；ps為插值得到的結構表面壓強；Nf和Ns分別為物面氣動網格總數和單元網格總數；Sf和Ss分別為流場物面網格單元面積和結構單元面積；Γ為物面總面積。為了保證插值的守恒性，應當滿足如下關系式：

(4)

引入修正系數W，定義為

(5)

對作用在結構單元表面的壓強進行修正，有

(6)

高超聲速飛行器在飛行過程中，結構溫度分布與剛度矩陣是時刻變化的。采用考慮溫度分布與剛度變化的直接數值積分方法求解結構運動方程能夠得到結構的瞬態響應。對于線性問題，質量矩陣M與剛度矩陣K為常數，在已知載荷的情況下可以直接得到應力和位移。而非線性瞬態分析需要將載荷分解為多段載荷增量加載，具體步驟如下：

步驟1先施加總載荷的一部分載荷。

步驟2提取剛度矩陣K。

步驟3求解剛度矩陣和殘值向量，判斷是否收斂。

步驟4如果收斂，則得到結構的位移和應力；如果不收斂，則重復步驟2～3，直至收斂。

步驟5繼續增加載荷增量，重復步驟1～4，直至100%載荷加載。

采用徑向基函數方法(RBFs)實現結構表面變形的傳遞過程，將結構單元點的位移插值到流場物面網格點。該方法根據空間中一些函數值已知的點構造出一種插值函數，再利用構造出的函數插值得到其他函數值未知的點。RBFs構造的函數是一系列基函數的疊加：

(7)

f(xb)=db

(8)

(9)

其中：db為控制點處的函數值。令任意2點的范數表達式為

(10)

系數ai，b0，b1，b2，b3可以通過求解下列方程組得到：

(11)

其中:Q為基函數φij組成的方陣；a、b為附加條件系數向量；P為與附加系數向量相乘的坐標矩陣，具體表達式為

(12)

(13)

(14)

在計算過程中，雖然結構表面點的位移在時刻發生變化，但是與氣動網格物面點之間的相對位置關系是保持不變的，這樣在計算開始前就可以確定并儲存RBFs的插值函數。在后續的計算中，結構發生任意變形都可以通過儲存的插值函數快速計算出流場物面網格的變形。

3 翼面偏轉與結構振動疊加的動網格

采用移動坐標系描述飛行器的變速運動。流場方程建立在移動坐標系上，即氣流坐標系(Oxayaza)，坐標原點一般位于飛行器質心，Oxa軸始終指向飛行器空速方向，Oya軸始終指向機翼展向。結構運動方程建立在機體坐標系(Oxbybzb)上，固聯于飛行器并隨飛行器一同運動，坐標原點同樣位于飛行器質心，Oxb軸位于機身軸線，或者平行于翼面的平均氣動弦線，Oyb軸始終指向機翼展向。2個坐標系的Oxbzb平面重合，下標a表示該點在氣流坐標系中的坐標；b表示該點在機體坐標系中的坐標。

飛行器姿態的改變可以采用動網格的方法實現，原點重合的2個三維坐標系可以通過繞坐標軸旋轉重合。繞x軸，y軸和z軸逆時針轉過α角的轉換矩陣可以分別表示為

(15)

(16)

通過結構運動方程計算得到的結構變形是在機體坐標系中的變形，需要轉換至氣流坐標系。結構有限元模型中的某一單元點變形前后在機體坐標系中的坐標表示為(xb,yb,zb)和(x′b,y′b,z′b)。其變形量表示成向量的形式為defb=[x′b-xb,y′b-yb,z′b-zb]T。當空氣來流迎角為α時，不考慮側滑角，機體坐標系繞Oyb逆時針旋轉角α即可與氣流坐標系重合。因此，翼面結構單元變形應該順時針旋轉α轉換至氣流坐標系。當Oyb與Oya不重合時，應當先在Oxbzb平面內平移，使Oyb與Oya重合。以(xa,ya,za)，(x′a,y′a,z′a)分別表示單元點變形前后，在氣流坐標系中的坐標，則有：

(17)

結構變形量在氣流坐標系中可以表示為

Lydefb

(18)

式(18)為結構單元點變形量從機體坐標系到氣流坐標系的變換方式。

將翼面偏轉得到的坐標與坐標變換后的變形量相加，可以得到當迎角為α時，變形后的物面某網格點坐標：

(19)

每一時刻，在氣流坐標系中，變形后的物面網格坐標都可以通過式(19)表示。得到物面坐標后，整個流場的網格坐標均可求出。相鄰時刻的網格坐標相減即可得到在氣流坐標系中物面網格的變形量。

對計算流程進行修正，如圖3所示。在第④步，將結構振動變形與增加翼面偏角引起的網格變形量疊加，傳遞給流場物面邊界，流場快速生成動網格。

圖3 修正后的多場耦合方法計算流程Fig.3 Revised process of multi-field coupling computation

4 熱氣動彈性穩定性分析方法

采用全時域耦合分析方法計算飛行器全軌道運動過程的計算量十分龐大。一般情況下，軌道狀態點時間間隔比結構振動響應特征時間大幾個數量級。結構振動若干周期對溫度場的影響較小。基于這一假設，將同步計算方法與全時域耦合分析方法相結合，建立了沿軌道運動的高超聲速飛行器熱氣動彈性穩定性分析方法。后續計算結果也對這一假設進行了驗證。計算流程如圖4所示。詳細流程如下：

圖4 沿軌道運動的高超聲速飛行器多場全時域耦合熱氣動彈性響應計算及穩定性分析流程Fig.4 Aerothermoelasticity response computation and stability analysis of an orbital moving hypersonic vehicle

① 首先采用能夠考慮飛行高度、姿態、馬赫數變化的流場-結構溫度場同步計算方法得到軌道各個狀態點的流場與結構溫度場。

② 讀入同步計算方法的結果，以同步計算方法得到的流場和結構溫度場作為相應狀態點的初場，采用全時域耦合方法計算軌道相鄰狀態點之間的結構振動響應。在推進至下一狀態點前，通過結構振動瞬態響應曲線可以判斷結構穩定與否。

5 算例與分析

5.1 計算模型

模型外形如圖5所示[8]。根部半弦長位置采用鉸鏈連接，扭轉剛度為1×107N·m/(°)。弦向中部厚度為弦長C的0.025倍。前緣鈍頭半徑為2.5×10-3m。結構溫度場初始溫度為293.15 K。 物面第1層網格高度為4×10-7m。材料密度為7 900 kg/m3，比熱容為500 J/(kg·K)，熱傳導系數為16.3 W/(m·K)。軌道狀態點參數如表1所示。

表1 軌道狀態點參數

圖5 翼面外形Fig.5 Configuration of the wing

5.2 多場耦合計算

采用同步計算方法計算翼面沿軌道運動的結構溫度。流場計算網格如圖6所示，計算模型如圖7和圖8所示。流場網格與同步計算傳熱的結構有限體積網格邊界一一對接。

圖6 翼面空間多塊結構網格Fig.6 Multi-block structured grid of wing

圖7 用于同步方法計算傳熱的有限體積模型Fig.7 Finite volume model for synchronization algorithm

圖8 用于響應分析的有限元模型Fig.8 Finite element model for response analysis

以同步計算方法得到的第3、第4和第5個狀態點流場和結構溫度場為初場[8]，采用瞬態多場耦合方法計算翼面振動響應。物理時間步為1×10-4s。圖9所示的4個角點為監測點。

圖9 監測點位置Fig.9 Monitoring points location

圖10～圖12為狀態點3到狀態點6之間各個監測點的振動響應曲線。在狀態點3到狀態點4之間，所有監測點均保持在一定幅值振動；在狀態點4到狀態點5之間，所有監測點位移曲線在初期均保持在一定幅值振動，在11.5 s左右，振動幅度不斷增加，逐步出現發散趨勢。在狀態點5與狀態點6之間，各個監測點位移曲線均呈發散狀態，結構振動幅度不斷增加，結構發生顫振。最終在12.5 s左右由于振動幅值過大，網格發生破壞，終止計算。1號與8號監測點位于翼根處，振動幅值較小；3號與9號監測點位于翼梢處，振動幅值較大。

圖10 狀態點3-4各監測點響應Fig.10 Transient response of monitoring points between states 3 and 4

圖11 狀態點4-5各監測點響應Fig.11 Transient response of monitoring points between states 4 and 5

圖12 狀態點5-6各監測點響應Fig.12 Transient response of monitoring points between states 5 and 6

為研究結構振動對熱流密度和結構溫度場的影響，采用同步計算方法計算剛體模型沿相應軌道運動過程中的熱流密度和溫度，并與本文的計算結果對比。圖13與圖14為1號和3號監測點的熱流密度變化曲線。從整體的趨勢來看，在狀態點3到狀態點5之間，熱流密度波動幅值并不明顯。在接近狀態點5時，隨著結構振動幅值逐漸增加，熱流密度波動的幅值也逐漸增加。在狀態點5與狀態點6之間，1號和3號監測點的響應幅值呈發散趨勢，熱流密度的波動幅度更大。3號監測點的響應波動的幅值約為1號監測點幅值的4倍，相同時刻的熱流密度波動的幅值也約為1號點的4倍。

圖13 1號監測點熱流隨時間變化曲線Fig.13 Time history of heat flux comparison on monitoring point 1

圖15與圖16所示為1號和3號監測點溫度隨時間變化曲線。結構溫度場對熱流的變化并不敏感，在接近狀態點5時，1號和3號監測點熱流的波動發散幅度已經比較明顯，而溫度曲線波動的幅值仍然很小。

圖15 1號監測點溫度隨時間變化Fig.15 Time history of temperature comparison on monitoring point 1

圖16 3號監測點溫度隨時間變化對比Fig.16 Time history of temperature comparison on monitoring point 3

在狀態點5-6，監測點的熱流密度發生明顯變化，3號監測點熱流密度波動幅值約占熱流峰值的10%左右，而溫度變化并不明顯，相比于剛體模型的監測點溫度只下降了0.3%左右。結構振動過程對前緣駐點的局部熱流密度的影響比較明顯，而對溫度分布影響較小。

1號和3號監測點的熱流密度和溫度曲線說明，結構發生顫振或者較大幅度的擺動之前，結構的振動對結構溫度場影響較小。在結構振動對溫度場產生明顯影響的時候，結構振動已經呈發散趨勢。采用同步計算方法計算得到的結構溫度場已經能夠滿足精度需求。由此可以認為，結構振動若干周期對溫度場的影響較小的假設是合理的。本文建立的先利用流場-結構溫度場同步計算方法獲得軌道全部狀態點的溫度，再采用多場全時域耦合方法計算狀態點之間的結構振動響應，可以用于沿軌道運動飛行器的熱氣動彈性響應計算和穩定性分析。

1號和8號監測點的響應曲線對比如圖17所示。可以看出曲線的波動頻率與幅值基本一致，相位角始終差180°，可以認為這2個點的振動規律由鉸鏈扭轉模態起主導作用；3號和9號監測點的響應曲線(圖18)波動的頻率與幅值也基本一致，存在一定的相位差，這兩點的運動規律由一階彎曲模態起主導作用。與“凍結”模態熱顫振方法[8]得出的鉸鏈扭轉和一階彎曲模態耦合顫振的結論是一致的。

圖17 1號和8號監測點響應對比Fig.17 Transient response comparison between monitoring point 1 and 8

圖18 3號和9號監測點響應對比Fig.18 Transient response comparison between monitoring point 3 and 9

“凍結”模態的熱顫振方法計算得到的顫振臨界點在軌道狀態點4-5之間[8]。多場耦合方法得到的顫振臨界點在同樣狀態點4-5之間，在11.5 s左右監測點響應開始發散，在穩定性方面2種方法計算結果基本一致。

6 結 論

本文在流場-結構溫度場同步計算方法的基礎上，建立了氣動力/氣動熱/結構溫度場/結構振動多物理場全時域耦合的熱氣動彈性分析方法。以全動翼面為研究對象，通過分析對比得到以下結論：

1) 與同步計算方法相比，全時域耦合方法能夠模擬結構振動對流場和結構溫度場的影響。監測點的熱流密度發生明顯變化，3號監測點波動幅值占熱流峰值的10%左右。而溫度變化并不明顯，相比于剛體模型，監測點溫度只下降了0.3%左右。高超聲速飛行器結構發生顫振或者較大幅度的擺動之前，結構的振動對結構溫度場影響較小。在對溫度場能夠產生明顯影響時，結構振動響應已經呈發散趨勢。采用同步計算方法計算得到的結構溫度場已經能夠滿足精度需求。所建立的方法可以用于熱氣動彈性穩定性分析。

2) 計算得到的顫振臨界點在4-5號狀態點之間，顫振形式為鉸鏈扭轉模態與一階彎曲模態的耦合顫振，與“凍結”模態的熱顫振方法結果一致。“凍結”模態的熱顫振計算只能離散地分析每個軌道狀態點的顫振特性，顫振臨界點的判斷精度依賴于軌道狀態點的選取。全時域耦合方法在時域內推進求解流場、結構溫度場和結構振動響應，溫度場計算能夠考慮結構振動的影響，能獲得沿軌道運動的飛行器比較明確的顫振臨界點。在應用方面，當飛行器處于流場參數和結構溫度場變化并不劇烈的飛行狀態時，比如巡航階段，可以采用 “凍結”模態的熱顫振計算方法；當飛行器處于流場參數和結構溫度場變化劇烈的階段，比如機動過程，需要開展響應計算時，可以采用全時域耦合分析方法。