高超聲速全動(dòng)翼面全時(shí)域耦合分析方法及應(yīng)用

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

1. 航空工業(yè)沈陽飛機(jī)設(shè)計(jì)研究所，沈陽 110035

2. 南京航空航天大學(xué) 航空學(xué)院，南京 210016

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

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

氣動(dòng)加熱形成的熱環(huán)境與構(gòu)成氣動(dòng)彈性問題的氣動(dòng)力、慣性力和彈性力形成的耦合關(guān)系強(qiáng)弱程度不同[2]，大多數(shù)熱氣動(dòng)彈性問題的研究方法都是忽略了結(jié)構(gòu)振動(dòng)對溫度場的直接和間接影響，將高超聲速飛行器飛行過程分解為若干狀態(tài)點(diǎn)，在每個(gè)狀態(tài)點(diǎn)上固定來流參數(shù)計(jì)算結(jié)構(gòu)振動(dòng)響應(yīng)。

隨著技術(shù)的發(fā)展，先進(jìn)的高超聲速飛行器具有變軌飛行能力，能夠在沿軌道飛行的過程中控制翼面偏轉(zhuǎn)來改變姿態(tài)和飛行軌道。操縱面作為最可能出現(xiàn)氣動(dòng)彈性失穩(wěn)的部件之一，將會(huì)成為高超聲速飛行器熱氣動(dòng)彈性分析的重點(diǎn)對象[20]。應(yīng)當(dāng)建立能夠模擬高超聲速飛行器機(jī)動(dòng)過程中的操縱面氣動(dòng)力、氣動(dòng)熱、熱傳導(dǎo)和結(jié)構(gòu)振動(dòng)的相互作用的方法。

本文針對沿軌道運(yùn)動(dòng)的高超聲速全動(dòng)翼面的偏轉(zhuǎn)過程，在流場-結(jié)構(gòu)溫度場同步計(jì)算方法的基礎(chǔ)上，建立了氣動(dòng)力/氣動(dòng)熱/結(jié)構(gòu)溫度場/結(jié)構(gòu)振動(dòng)多物理場全時(shí)域耦合的計(jì)算分析方法。通過改變遠(yuǎn)場來流密度、壓強(qiáng)和溫度來模擬飛行軌道高度變化，采用移動(dòng)坐標(biāo)系和動(dòng)網(wǎng)格相結(jié)合的方式分別描述變速度飛行和翼面偏轉(zhuǎn)過程。通過坐標(biāo)系變換將翼面偏轉(zhuǎn)和振動(dòng)的網(wǎng)格變形量疊加，考慮翼面振動(dòng)和偏轉(zhuǎn)過程的耦合非定常效應(yīng)。通過結(jié)構(gòu)振動(dòng)時(shí)域響應(yīng)分析研究了結(jié)構(gòu)振動(dòng)響應(yīng)對熱流密度和溫度分布的影響。

1 多物理場全時(shí)域耦合計(jì)算流程

熱氣動(dòng)彈性問題涉及多物理場之間的耦合作用，難以采用統(tǒng)一方程求解，在工程上一般采用分區(qū)求解-邊界耦合的方法，包括氣動(dòng)力/氣動(dòng)熱計(jì)算、熱傳導(dǎo)計(jì)算以及熱結(jié)構(gòu)振動(dòng)響應(yīng)分析等。各物理場耦合關(guān)系緊密，如圖1所示。

圖1 熱氣動(dòng)彈性多物理場耦合關(guān)系Fig.1 Aerothermoelasticity multiphysics coupling relationship

作用在結(jié)構(gòu)表面的熱流會(huì)使結(jié)構(gòu)升溫，高溫會(huì)影響結(jié)構(gòu)的材料屬性，同時(shí)使結(jié)構(gòu)產(chǎn)生附加的熱應(yīng)力，改變結(jié)構(gòu)剛度特性，進(jìn)而影響結(jié)構(gòu)的振動(dòng)響應(yīng)特性。結(jié)構(gòu)表面的壓強(qiáng)和熱流分布對結(jié)構(gòu)振動(dòng)響應(yīng)比較敏感，結(jié)構(gòu)振動(dòng)特性的改變會(huì)通過改變熱流分布間接影響結(jié)構(gòu)溫度場。本文建立的全時(shí)域耦合計(jì)算流程如圖2所示。詳細(xì)流程如下：

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

① 讀入結(jié)構(gòu)初始外形(初始結(jié)構(gòu)變形def0一般為0)，以及初始溫度T0。流場與結(jié)構(gòu)溫度場同步推進(jìn)求解，得到下一時(shí)刻的流場與結(jié)構(gòu)溫度場。

② 將第①步計(jì)算得到的結(jié)構(gòu)表面壓強(qiáng)作為載荷作用在有限元模型表面，同時(shí)將計(jì)算得到的結(jié)構(gòu)溫度場分布T1賦值給有限元模型。

③ 考慮結(jié)構(gòu)瞬時(shí)溫度場的非線性影響，計(jì)算結(jié)構(gòu)瞬態(tài)響應(yīng)，得到t1時(shí)刻結(jié)構(gòu)響應(yīng)def1。

④ 將該時(shí)刻結(jié)構(gòu)振動(dòng)產(chǎn)生變形傳遞給流場物面邊界，流場快速生成動(dòng)態(tài)網(wǎng)格。

⑤ 重復(fù)①-④步，各物理場沿時(shí)間推進(jìn)計(jì)算。

2 各物理場計(jì)算與耦合方法

2.1 氣動(dòng)熱-結(jié)構(gòu)溫度場同步計(jì)算

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

對于結(jié)構(gòu)熱傳導(dǎo)，由于固體介質(zhì)中的熱傳導(dǎo)過程不存在對流現(xiàn)象，顯然直角坐標(biāo)系3個(gè)速度方向分量為0，即u=v=w=0，連續(xù)方程與動(dòng)量方程自然滿足。固體域能量方程為

(1)

式中：T為溫度；t為時(shí)間；k為第1層單元熱傳導(dǎo)系數(shù)；ρ為密度；C為比熱容；x、y、z為直角坐標(biāo)系3個(gè)方向分量；下標(biāo)s表示固體介質(zhì)。

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

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

(2)

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

2.2 氣動(dòng)力-結(jié)構(gòu)瞬態(tài)響應(yīng)耦合計(jì)算

得到某時(shí)刻的流場和結(jié)構(gòu)溫度場后，通過2個(gè)模型之間的映射關(guān)系實(shí)現(xiàn)第②步的結(jié)構(gòu)溫度場賦值和氣動(dòng)力插值，將有限體積模型的溫度場和物面邊界壓強(qiáng)傳遞給結(jié)構(gòu)有限元模型。

用于計(jì)算結(jié)構(gòu)溫度場的有限體積模型的網(wǎng)格密度與計(jì)算響應(yīng)的有限元模型的網(wǎng)格密度差別較大，溫度場和流場計(jì)算模型的網(wǎng)格要遠(yuǎn)密于有限元模型的單元。采用最鄰近搜索法建立有限體積模型到有限元模型的映射關(guān)系。在每一個(gè)物理時(shí)間步，對流場物面壓強(qiáng)和插值得到的結(jié)構(gòu)表面壓強(qiáng)進(jìn)行積分：

(3)

式中：pf為計(jì)算得到流場物面壓強(qiáng)；ps為插值得到的結(jié)構(gòu)表面壓強(qiáng)；Nf和Ns分別為物面氣動(dòng)網(wǎng)格總數(shù)和單元網(wǎng)格總數(shù)；Sf和Ss分別為流場物面網(wǎng)格單元面積和結(jié)構(gòu)單元面積；Γ為物面總面積。為了保證插值的守恒性，應(yīng)當(dāng)滿足如下關(guān)系式：

(4)

引入修正系數(shù)W，定義為

(5)

對作用在結(jié)構(gòu)單元表面的壓強(qiáng)進(jìn)行修正，有

(6)

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

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

步驟2提取剛度矩陣K。

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

步驟4如果收斂，則得到結(jié)構(gòu)的位移和應(yīng)力；如果不收斂，則重復(fù)步驟2～3，直至收斂。

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

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

(7)

f(xb)=db

(8)

(9)

其中：db為控制點(diǎn)處的函數(shù)值。令任意2點(diǎn)的范數(shù)表達(dá)式為

(10)

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

(11)

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

(12)

(13)

(14)

在計(jì)算過程中，雖然結(jié)構(gòu)表面點(diǎn)的位移在時(shí)刻發(fā)生變化，但是與氣動(dòng)網(wǎng)格物面點(diǎn)之間的相對位置關(guān)系是保持不變的，這樣在計(jì)算開始前就可以確定并儲(chǔ)存RBFs的插值函數(shù)。在后續(xù)的計(jì)算中，結(jié)構(gòu)發(fā)生任意變形都可以通過儲(chǔ)存的插值函數(shù)快速計(jì)算出流場物面網(wǎng)格的變形。

3 翼面偏轉(zhuǎn)與結(jié)構(gòu)振動(dòng)疊加的動(dòng)網(wǎng)格

采用移動(dòng)坐標(biāo)系描述飛行器的變速運(yùn)動(dòng)。流場方程建立在移動(dòng)坐標(biāo)系上，即氣流坐標(biāo)系(Oxayaza)，坐標(biāo)原點(diǎn)一般位于飛行器質(zhì)心，Oxa軸始終指向飛行器空速方向，Oya軸始終指向機(jī)翼展向。結(jié)構(gòu)運(yùn)動(dòng)方程建立在機(jī)體坐標(biāo)系(Oxbybzb)上，固聯(lián)于飛行器并隨飛行器一同運(yùn)動(dòng)，坐標(biāo)原點(diǎn)同樣位于飛行器質(zhì)心，Oxb軸位于機(jī)身軸線，或者平行于翼面的平均氣動(dòng)弦線，Oyb軸始終指向機(jī)翼展向。2個(gè)坐標(biāo)系的Oxbzb平面重合，下標(biāo)a表示該點(diǎn)在氣流坐標(biāo)系中的坐標(biāo)；b表示該點(diǎn)在機(jī)體坐標(biāo)系中的坐標(biāo)。

飛行器姿態(tài)的改變可以采用動(dòng)網(wǎng)格的方法實(shí)現(xiàn)，原點(diǎn)重合的2個(gè)三維坐標(biāo)系可以通過繞坐標(biāo)軸旋轉(zhuǎn)重合。繞x軸，y軸和z軸逆時(shí)針轉(zhuǎn)過α角的轉(zhuǎn)換矩陣可以分別表示為

(15)

(16)

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

(17)

結(jié)構(gòu)變形量在氣流坐標(biāo)系中可以表示為

Lydefb

(18)

式(18)為結(jié)構(gòu)單元點(diǎn)變形量從機(jī)體坐標(biāo)系到氣流坐標(biāo)系的變換方式。

將翼面偏轉(zhuǎn)得到的坐標(biāo)與坐標(biāo)變換后的變形量相加，可以得到當(dāng)迎角為α?xí)r，變形后的物面某網(wǎng)格點(diǎn)坐標(biāo)：

(19)

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

對計(jì)算流程進(jìn)行修正，如圖3所示。在第④步，將結(jié)構(gòu)振動(dòng)變形與增加翼面偏角引起的網(wǎng)格變形量疊加，傳遞給流場物面邊界，流場快速生成動(dòng)網(wǎng)格。

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

4 熱氣動(dòng)彈性穩(wěn)定性分析方法

采用全時(shí)域耦合分析方法計(jì)算飛行器全軌道運(yùn)動(dòng)過程的計(jì)算量十分龐大。一般情況下，軌道狀態(tài)點(diǎn)時(shí)間間隔比結(jié)構(gòu)振動(dòng)響應(yīng)特征時(shí)間大幾個(gè)數(shù)量級(jí)。結(jié)構(gòu)振動(dòng)若干周期對溫度場的影響較小?；谶@一假設(shè)，將同步計(jì)算方法與全時(shí)域耦合分析方法相結(jié)合，建立了沿軌道運(yùn)動(dòng)的高超聲速飛行器熱氣動(dòng)彈性穩(wěn)定性分析方法。后續(xù)計(jì)算結(jié)果也對這一假設(shè)進(jìn)行了驗(yàn)證。計(jì)算流程如圖4所示。詳細(xì)流程如下：

圖4 沿軌道運(yùn)動(dòng)的高超聲速飛行器多場全時(shí)域耦合熱氣動(dòng)彈性響應(yīng)計(jì)算及穩(wěn)定性分析流程Fig.4 Aerothermoelasticity response computation and stability analysis of an orbital moving hypersonic vehicle

① 首先采用能夠考慮飛行高度、姿態(tài)、馬赫數(shù)變化的流場-結(jié)構(gòu)溫度場同步計(jì)算方法得到軌道各個(gè)狀態(tài)點(diǎn)的流場與結(jié)構(gòu)溫度場。

② 讀入同步計(jì)算方法的結(jié)果，以同步計(jì)算方法得到的流場和結(jié)構(gòu)溫度場作為相應(yīng)狀態(tài)點(diǎn)的初場，采用全時(shí)域耦合方法計(jì)算軌道相鄰狀態(tài)點(diǎn)之間的結(jié)構(gòu)振動(dòng)響應(yīng)。在推進(jìn)至下一狀態(tài)點(diǎn)前，通過結(jié)構(gòu)振動(dòng)瞬態(tài)響應(yīng)曲線可以判斷結(jié)構(gòu)穩(wěn)定與否。

5 算例與分析

5.1 計(jì)算模型

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

表1 軌道狀態(tài)點(diǎn)參數(shù)

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

5.2 多場耦合計(jì)算

采用同步計(jì)算方法計(jì)算翼面沿軌道運(yùn)動(dòng)的結(jié)構(gòu)溫度。流場計(jì)算網(wǎng)格如圖6所示，計(jì)算模型如圖7和圖8所示。流場網(wǎng)格與同步計(jì)算傳熱的結(jié)構(gòu)有限體積網(wǎng)格邊界一一對接。

圖6 翼面空間多塊結(jié)構(gòu)網(wǎng)格Fig.6 Multi-block structured grid of wing

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

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

以同步計(jì)算方法得到的第3、第4和第5個(gè)狀態(tài)點(diǎn)流場和結(jié)構(gòu)溫度場為初場[8]，采用瞬態(tài)多場耦合方法計(jì)算翼面振動(dòng)響應(yīng)。物理時(shí)間步為1×10-4s。圖9所示的4個(gè)角點(diǎn)為監(jiān)測點(diǎn)。

圖9 監(jiān)測點(diǎn)位置Fig.9 Monitoring points location

圖10～圖12為狀態(tài)點(diǎn)3到狀態(tài)點(diǎn)6之間各個(gè)監(jiān)測點(diǎn)的振動(dòng)響應(yīng)曲線。在狀態(tài)點(diǎn)3到狀態(tài)點(diǎn)4之間，所有監(jiān)測點(diǎn)均保持在一定幅值振動(dòng)；在狀態(tài)點(diǎn)4到狀態(tài)點(diǎn)5之間，所有監(jiān)測點(diǎn)位移曲線在初期均保持在一定幅值振動(dòng)，在11.5 s左右，振動(dòng)幅度不斷增加，逐步出現(xiàn)發(fā)散趨勢。在狀態(tài)點(diǎn)5與狀態(tài)點(diǎn)6之間，各個(gè)監(jiān)測點(diǎn)位移曲線均呈發(fā)散狀態(tài)，結(jié)構(gòu)振動(dòng)幅度不斷增加，結(jié)構(gòu)發(fā)生顫振。最終在12.5 s左右由于振動(dòng)幅值過大，網(wǎng)格發(fā)生破壞，終止計(jì)算。1號(hào)與8號(hào)監(jiān)測點(diǎn)位于翼根處，振動(dòng)幅值較小；3號(hào)與9號(hào)監(jiān)測點(diǎn)位于翼梢處，振動(dòng)幅值較大。

圖10 狀態(tài)點(diǎn)3-4各監(jiān)測點(diǎn)響應(yīng)Fig.10 Transient response of monitoring points between states 3 and 4

圖11 狀態(tài)點(diǎn)4-5各監(jiān)測點(diǎn)響應(yīng)Fig.11 Transient response of monitoring points between states 4 and 5

圖12 狀態(tài)點(diǎn)5-6各監(jiān)測點(diǎn)響應(yīng)Fig.12 Transient response of monitoring points between states 5 and 6

為研究結(jié)構(gòu)振動(dòng)對熱流密度和結(jié)構(gòu)溫度場的影響，采用同步計(jì)算方法計(jì)算剛體模型沿相應(yīng)軌道運(yùn)動(dòng)過程中的熱流密度和溫度，并與本文的計(jì)算結(jié)果對比。圖13與圖14為1號(hào)和3號(hào)監(jiān)測點(diǎn)的熱流密度變化曲線。從整體的趨勢來看，在狀態(tài)點(diǎn)3到狀態(tài)點(diǎn)5之間，熱流密度波動(dòng)幅值并不明顯。在接近狀態(tài)點(diǎn)5時(shí)，隨著結(jié)構(gòu)振動(dòng)幅值逐漸增加，熱流密度波動(dòng)的幅值也逐漸增加。在狀態(tài)點(diǎn)5與狀態(tài)點(diǎn)6之間，1號(hào)和3號(hào)監(jiān)測點(diǎn)的響應(yīng)幅值呈發(fā)散趨勢，熱流密度的波動(dòng)幅度更大。3號(hào)監(jiān)測點(diǎn)的響應(yīng)波動(dòng)的幅值約為1號(hào)監(jiān)測點(diǎn)幅值的4倍，相同時(shí)刻的熱流密度波動(dòng)的幅值也約為1號(hào)點(diǎn)的4倍。

圖13 1號(hào)監(jiān)測點(diǎn)熱流隨時(shí)間變化曲線Fig.13 Time history of heat flux comparison on monitoring point 1

圖15與圖16所示為1號(hào)和3號(hào)監(jiān)測點(diǎn)溫度隨時(shí)間變化曲線。結(jié)構(gòu)溫度場對熱流的變化并不敏感，在接近狀態(tài)點(diǎn)5時(shí)，1號(hào)和3號(hào)監(jiān)測點(diǎn)熱流的波動(dòng)發(fā)散幅度已經(jīng)比較明顯，而溫度曲線波動(dòng)的幅值仍然很小。

圖15 1號(hào)監(jiān)測點(diǎn)溫度隨時(shí)間變化Fig.15 Time history of temperature comparison on monitoring point 1

圖16 3號(hào)監(jiān)測點(diǎn)溫度隨時(shí)間變化對比Fig.16 Time history of temperature comparison on monitoring point 3

在狀態(tài)點(diǎn)5-6，監(jiān)測點(diǎn)的熱流密度發(fā)生明顯變化，3號(hào)監(jiān)測點(diǎn)熱流密度波動(dòng)幅值約占熱流峰值的10%左右，而溫度變化并不明顯，相比于剛體模型的監(jiān)測點(diǎn)溫度只下降了0.3%左右。結(jié)構(gòu)振動(dòng)過程對前緣駐點(diǎn)的局部熱流密度的影響比較明顯，而對溫度分布影響較小。

1號(hào)和3號(hào)監(jiān)測點(diǎn)的熱流密度和溫度曲線說明，結(jié)構(gòu)發(fā)生顫振或者較大幅度的擺動(dòng)之前，結(jié)構(gòu)的振動(dòng)對結(jié)構(gòu)溫度場影響較小。在結(jié)構(gòu)振動(dòng)對溫度場產(chǎn)生明顯影響的時(shí)候，結(jié)構(gòu)振動(dòng)已經(jīng)呈發(fā)散趨勢。采用同步計(jì)算方法計(jì)算得到的結(jié)構(gòu)溫度場已經(jīng)能夠滿足精度需求。由此可以認(rèn)為，結(jié)構(gòu)振動(dòng)若干周期對溫度場的影響較小的假設(shè)是合理的。本文建立的先利用流場-結(jié)構(gòu)溫度場同步計(jì)算方法獲得軌道全部狀態(tài)點(diǎn)的溫度，再采用多場全時(shí)域耦合方法計(jì)算狀態(tài)點(diǎn)之間的結(jié)構(gòu)振動(dòng)響應(yīng)，可以用于沿軌道運(yùn)動(dòng)飛行器的熱氣動(dòng)彈性響應(yīng)計(jì)算和穩(wěn)定性分析。

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

圖17 1號(hào)和8號(hào)監(jiān)測點(diǎn)響應(yīng)對比Fig.17 Transient response comparison between monitoring point 1 and 8

圖18 3號(hào)和9號(hào)監(jiān)測點(diǎn)響應(yīng)對比Fig.18 Transient response comparison between monitoring point 3 and 9

“凍結(jié)”模態(tài)的熱顫振方法計(jì)算得到的顫振臨界點(diǎn)在軌道狀態(tài)點(diǎn)4-5之間[8]。多場耦合方法得到的顫振臨界點(diǎn)在同樣狀態(tài)點(diǎn)4-5之間，在11.5 s左右監(jiān)測點(diǎn)響應(yīng)開始發(fā)散，在穩(wěn)定性方面2種方法計(jì)算結(jié)果基本一致。

6 結(jié) 論

本文在流場-結(jié)構(gòu)溫度場同步計(jì)算方法的基礎(chǔ)上，建立了氣動(dòng)力/氣動(dòng)熱/結(jié)構(gòu)溫度場/結(jié)構(gòu)振動(dòng)多物理場全時(shí)域耦合的熱氣動(dòng)彈性分析方法。以全動(dòng)翼面為研究對象，通過分析對比得到以下結(jié)論：

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

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