蔡巧言 呂計男 郭力 王昕江

(1 中國運載火箭技術(shù)研究院,北京 100076；2 北京宇航系統(tǒng)工程研究所,北京 100076；3 中國航天空氣動力技術(shù)研究院,北京 100074）

0 引言

顫振作為復(fù)雜的氣動彈性現(xiàn)象，一旦發(fā)生往往在幾個周期內(nèi)就會使飛行器結(jié)構(gòu)發(fā)生破壞。國家軍用標(biāo)準(zhǔn)與民航適航條例中對飛行器顫振分析均有明確的要求。

升力式布局面對稱運載器外形復(fù)雜，往往由異型截面機身、大翼展機翼、多組氣動舵面組成，飛行剖面具有大空域、寬速域的特點，因此運載器面臨的氣動力環(huán)境非常復(fù)雜。特別是運載器返回減速需求使得此類飛行器返回飛行段多采用大攻角飛行形式，當(dāng)運載器大攻角飛行時運載器機身對舵面產(chǎn)生較大干擾，使得舵面前方來流不再表現(xiàn)為自由來流流動結(jié)構(gòu)，特別是運載器后部的V 尾幾乎完全處于前體尾流干擾區(qū)內(nèi)，進(jìn)一步增加了其繞流流場的復(fù)雜性。同時，運載器結(jié)構(gòu)廣泛使用先進(jìn)復(fù)合材料，V 尾彈性頻率低，模態(tài)振型復(fù)雜。在復(fù)雜的流動環(huán)境與結(jié)構(gòu)動力學(xué)特性共同作用下，V 尾結(jié)構(gòu)易與氣動力耦合發(fā)生顫振，因此需要在方案設(shè)計初期即對其顫振特性進(jìn)行準(zhǔn)確分析，減小飛行風(fēng)險。傳統(tǒng)顫振工程計算方法對小攻角條件下氣動控制面的亞聲速/超聲速顫振分析結(jié)果可滿足工程初步設(shè)計階段的精度要求。對于跨聲速域及復(fù)雜擾流環(huán)境，工程氣動力模型難以實現(xiàn)非定常氣動力的準(zhǔn)確描述。楊超[1]對高超聲速氣動彈性力學(xué)研究中的不同工程氣動力模型及CFD/CSD 時域計算方法進(jìn)行了綜述，給出了不同方法的適用范圍與計算優(yōu)劣。楊炳淵等[2]針對翼面大攻角顫振給出了一種當(dāng)?shù)鼗钊碚撚嬎惴椒?，使用算例證明了該方法在有限厚度翼面大攻角顫振中的可行性。金偉[3]對先進(jìn)戰(zhàn)斗機全動 V 尾抖振強度進(jìn)行了研究，使用RANS/LES 混合方法對大迎角飛行狀態(tài)下的V 尾擾流流場進(jìn)行計算，并使用CFD/CSD 方法對其抖振特性進(jìn)行了分析。對于升力式布局面對稱運載器超聲速大攻角返回飛行而言，受大攻角飛行下大范圍分離非定常流場復(fù)雜流動特性影響，傳統(tǒng)工程氣動力方法難以實現(xiàn)流場準(zhǔn)確計算。設(shè)計中需要采用基于計算流體力學(xué)/計算結(jié)構(gòu)動力學(xué)（CFD/CSD）耦合求解的數(shù)值方法，在精確模擬流場的基礎(chǔ)上研究運載器顫振特性。

本文分別采用基于CFD 和工程氣動力模型的顫振分析方法，對升力式運載器大攻角飛行狀態(tài)下V 尾的顫振特性進(jìn)行計算，對兩種方法計算結(jié)果進(jìn)行對比，分析了工程氣動力方法在此類飛行器大攻角顫振分析中的不足，提出類似問題需要在工程研制中引起關(guān)注。

1 顫振分析方法

常用的顫振仿真分析方法有工程氣動力[4,5]方法和CFD/CSD 方法[6-12]。

1.1 工程氣動力方法

時域顫振計算的控制方程可以表示為式（1）所示

式中，M為質(zhì)量矩陣，D為阻尼矩陣，K為剛度矩陣，u為結(jié)構(gòu)節(jié)點位移向量，Ma∞為來流馬赫數(shù)，k為減縮頻率，F(xiàn)a為作用在結(jié)構(gòu)上的氣動力。

根據(jù)工程氣動力模型，每個面元上的壓力值與面元下洗的關(guān)系可以表示成影響系數(shù)的形式，此時Fa可以表示為

式中，Ma∞為來流馬赫數(shù)，k為減縮頻率，ρ∞為來流密度，U∞為來流速度。A為氣動力影響系數(shù)矩陣，是減縮頻率與來流馬赫數(shù)的函數(shù)。w為面元下洗向量。實際分析中，一般通過模態(tài)坐標(biāo)投影，對模態(tài)坐標(biāo)下的氣動彈性控制方程進(jìn)行求解

式中，M,D,K分別為模態(tài)質(zhì)量矩陣、模態(tài)阻尼矩陣與模態(tài)剛度矩陣。q為廣義坐標(biāo)向量，Q=φTFa為廣義氣動力向量。φ為模態(tài)振型向量組成的振型矩陣。分析中采用基于統(tǒng)一升力面理論（ZONA7U）的工程氣動力模型。通過求解式(3)得到飛行器顫振特性。

1.2 CFD/CSD 耦合分析方法

CFD/CSD 耦合分析方法中CFD 計算模型對V 尾顫振特性的準(zhǔn)確預(yù)測至關(guān)重要。對于沒有氣動干擾的氣動控制面，只需考慮控制面運動產(chǎn)生的非定常氣動力，采用基于雷諾平均（RANS）模式精確模擬近壁面的流場結(jié)構(gòu)即可滿足分析要求。但對于大攻角飛行時處在飛行器前體干擾下的V 尾，僅模擬壁面附近的流場不能準(zhǔn)確反映前機身脫體渦對V 尾的影響，需采用更精確的流場模擬數(shù)值方法[3,7]。大渦模擬方法（LES）方法能精確求解大范圍的氣流分離和不穩(wěn)定流動，但由于計算效率低下的問題，不適用于實際工程問題分析。本文采用RANS/LES 的DES（Detached Eddy Simulation）混合算法，分析V 尾的顫振模式。DES 方法在氣動控制面壁面區(qū)域具有RANS方法的性能，在流動分離區(qū)表現(xiàn)為LES 方程的性能[3,10,12]。V 尾顫振具有模態(tài)頻率低、多模態(tài)耦合和響應(yīng)位移大的特點。本文在CFD/CSD 耦合求解方法引入預(yù)估矯正迭代方法進(jìn)一步提高分析精度，采用基于RBF 方法的網(wǎng)格變形技術(shù)保證求解魯棒性[8]。CFD/CSD 計算流程示意圖如圖1 所示。

圖1 CFD/CSD 計算流程示意圖Fig.1 CFD/CSD calculation process diagram

2 結(jié)構(gòu)動力學(xué)模型

V 尾顫振計算需要考慮所有和V 尾相關(guān)的彈性模態(tài)。使用Patran 軟件的Cquad4 和Cquad8 殼單元模擬蒙皮，使用殼單元與Cbar 梁單元模擬主要承力結(jié)構(gòu)，使用Conm2 集中質(zhì)量單元模擬主要質(zhì)量部件，建立V 尾有限元動力學(xué)模型，如圖2所示，設(shè)置自由-自由邊界條件。

圖2 V 尾有限元模型Fig.2 FEM model of V-tail

采用Lanczos 方法對有限元動力學(xué)模型進(jìn)行動力學(xué)分析，獲得所有和V 尾相關(guān)的彈性模態(tài)頻率和振型，其中前4 階為V 尾組合模態(tài)，第5 階為V 尾扭轉(zhuǎn)模態(tài)。CFD 顫振計算所用三維振型如表1 第3 行所示。使用薄板樣條TPS 方法（Tin Plate Spline Method）將所有和V 尾相關(guān)的彈性模態(tài)振型向面元網(wǎng)格插值，工程氣動力方法顫振計算所用振型如表1 第4 行所示。

表1 V 尾彈性模態(tài)振型Table 1 Modal shapes of V-tail

3 顫振分析

選取大攻角返回特征飛行狀態(tài)，馬赫數(shù)2.3攻角18°。分別采用工程氣動力/CSD 方法和CFD/CSD 方法開展顫振分析。

3.1 基于工程方法的顫振分析

采用ZONA7U 方法建立超聲速非定常氣動力模型，采用表1 中V 尾彈性模態(tài)，在頻域進(jìn)行顫振分析。分析得到各階模態(tài)“動壓-阻尼”曲線與“動壓-頻率”曲線如圖3 所示。

圖3 各階模態(tài)“動壓-阻尼/頻率”曲線Fig.3 Dynamic pressure-damping/frequency curve

分析“動壓-阻尼”曲線與“動壓-頻率”曲線可見，工程方法計算得到的顫振動壓為11.3kPa。其中，第4 階模態(tài)動壓-阻尼曲線過零。通過模態(tài)剔除分析，確定第4 階、第5 階V 尾模態(tài)為顫振的主要參與模態(tài)。

3.2 基于CFD 方法的顫振分析

相同來流狀態(tài)馬赫數(shù)2.3、攻角18°，采用CFD 與V 尾彈性模態(tài)耦合的分析方法開展V 尾顫振分析。CFD 計算采用DES 模型，結(jié)構(gòu)模態(tài)阻尼取0，采用變動壓分析方法開展CFD/CSD 計算。通過各階廣義模態(tài)位移的時域響應(yīng)，判斷顫振動壓與發(fā)散模態(tài)。來流動壓11.3kPa 時，CFD/CSD 計算結(jié)果顯示（圖4）V 尾第4 階、第5 階彈性模態(tài)均未發(fā)生顫振發(fā)散。不斷提高來流動壓進(jìn)行顫振計算，直到來流動壓達(dá)到22.2kPa 時，V尾第4 階、第5 階彈性模態(tài)的廣義坐標(biāo)才出現(xiàn)小幅近似等幅振動。來流動壓22.2kPa 時，V 尾5階廣義模態(tài)位移時域響應(yīng)曲線如圖5 所示。其中，第1～3 階廣義模態(tài)位移時域響應(yīng)均為收斂狀態(tài)，第4 階模態(tài)呈現(xiàn)緩慢發(fā)散趨勢，第5 階模態(tài)呈現(xiàn)極限環(huán)等幅震蕩趨勢。

圖4 廣義模態(tài)位移時域響應(yīng)（11.3kPa）Fig.4 Generalized displacements in time domain(11.3kPa)

圖5 廣義模態(tài)位移時域響應(yīng)（22.2kPa）Fig.5 Generalized displacements in time domain (22.2kPa)

不同動壓第4 階廣義位移對比如圖6 所示，可見，11.3kPa 來流動壓時，第4 階廣義位移呈現(xiàn)周期性緩慢衰減狀態(tài)，阻尼比約0.2%；22.2kPa來流動壓時第4 階廣義位移呈線性小幅增長狀態(tài)，阻尼比約-0.4%。在顫振點附近，阻尼呈現(xiàn)非線性特征。

圖6 不同動壓下V 尾廣義位移Fig.6 Generalized displacements in time domain of V-tail

3.3 顫振特性差異分析

工程氣動力方法與CFD 方法預(yù)示顫振動壓分別為11.3kPa 和22.2kPa，兩者相差將近一倍，V 尾顫振特性差別較大。對流場進(jìn)行精確分析，V 尾三個不同站位（Y=0.68/1.29/1.69m）截面流場當(dāng)?shù)貏訅喝鐖D7 所示，其中，Y為飛行器全局坐標(biāo)系法向。自由來流經(jīng)運載器頭部脫體激波后，流場壓力上升，速度下降。經(jīng)過機身產(chǎn)生的膨脹波后流動加速，壓力下降。V 尾附近流動壓縮產(chǎn)生激波，使得V尾附近流場處于局部亞聲速狀態(tài)。由于大攻角逆壓梯度存在、激波附面層干擾，在V 尾前部產(chǎn)生大范圍分離，進(jìn)而影響V 尾附近流場分布，流動結(jié)構(gòu)非常復(fù)雜。在左側(cè) V 尾Y=0.68/1.29/1.69m（圖7）的展向站位上，分別取靠近前緣上表面、后緣上表面、前緣下表面與后緣下表面處的4 個特征點，對特征點處的來流動壓進(jìn)行測量。隨著展向站位Y 值增加，前緣上表面特征點處的當(dāng)?shù)貏訅悍謩e為1.13kPa、17.75kPa與19.32kPa；后緣上表面特征點處的當(dāng)?shù)貏訅悍謩e為6.61kPa、8.86kPa 與7.29kPa；前緣下表面特征點處的當(dāng)?shù)貏訅悍謩e為22.67kPa、18.12kPa 與20.50kPa；后緣下表面當(dāng)?shù)貏訅悍謩e為12.95kPa、14.09kPa 與16.88kPa。三個站位截面當(dāng)?shù)貏訅号c自由來流動壓22.2kPa 存在較大差異；此外，隨著展向站位變化，各特征點處當(dāng)?shù)貏訅褐档拇笮∫舶l(fā)生變化。

圖7 V 尾典型截面流場云圖Fig.7 Dynamic pressure contour of typical section

對比工程氣動力方法與CFD 方法得到的壓力系數(shù)（圖8）。根據(jù)CFD 計算結(jié)果，左側(cè)V 尾背風(fēng)面的壓力系數(shù)與其迎風(fēng)面相比可以近似視為均勻分布。此時，截取左側(cè)V 尾不同展向站位截面，并提取截面前緣、中部與后緣特征點處的表面壓力系數(shù)，與背風(fēng)面對應(yīng)位置壓力系數(shù)相減，可以得到不同特征點處CFD 計算得到的表面壓差系數(shù)。從根部特征截面到梢部特征截面，CFD計算左側(cè)V尾前緣特征點處的表面壓差系數(shù)分別為0.96、0.14 和0.03；中部特征點處的表面壓差系數(shù)分別為0.61、0.23 和0.07；后緣特征點處的表面壓差系數(shù)分別為0.51、0.28 和0.11。

不同特征點處CFD 計算的壓力系數(shù)顯示，左側(cè)V尾前緣高壓力系數(shù)區(qū)域隨著展向站位的增加不斷向后緣擴展，這與工程氣動力方法計算得到的表面壓差系數(shù)分布規(guī)律存在較大差異。

圖8 為兩種方法得到的V 尾表面壓力云圖。經(jīng)分析，在計算狀態(tài)下，運載器頭部激波、機身壓縮作用，在V 尾前方已出現(xiàn)激波（自身產(chǎn)生激波），使其表面局部出現(xiàn)亞聲速區(qū)，并出現(xiàn)分離區(qū)，工程氣動力方法已不能反映該流動結(jié)構(gòu)。CFD方法比較準(zhǔn)確的預(yù)示了當(dāng)?shù)氐牧鲌鼋Y(jié)構(gòu)，對展向和弦向流動引起的壓力變化預(yù)示準(zhǔn)確，較好的描述了前緣高壓區(qū)和尾緣及翼尖的低壓區(qū)的位置。

不同馬赫數(shù)、飛行攻角前體對V 尾的干擾特性不同，V 尾的當(dāng)?shù)貏訅?壓力分布變化規(guī)律也不同。采用固定馬赫數(shù)（Ma2.3）和固定來流動壓（22.2kPa）變換攻角（8°/18°/28°）的方式研究攻角的影響（圖9）。整體規(guī)律上看，在所計算攻角范圍內(nèi)，V 尾都會發(fā)生小阻尼發(fā)散。不同的是，隨著攻角增大，V 尾更早進(jìn)入等幅震蕩狀態(tài)，這主要是不同迎角下，前體對V 尾的干擾不同，V尾局部亞聲速區(qū)域范圍不同導(dǎo)致。

圖9 攻角對顫振的影響Fig.9 Influence of angle of attack

總之，通過對流場進(jìn)行精確分析，發(fā)現(xiàn)工程氣動力方法與CFD 方法預(yù)示顫振動壓產(chǎn)生較大偏差的直接原因是兩種分析方法得到的壓力分布和當(dāng)?shù)貏訅翰灰恢聦?dǎo)致，工程方法沒有反映出機身波系及大攻角分離對舵面前方來流流場的干擾影響。

4 結(jié)論

本文通過分析升力式布局面對稱飛行器V尾大攻角飛行狀態(tài)下的顫振特性，發(fā)現(xiàn)工程氣動力頻域方法和CFD/CSD 時域方法在預(yù)示中存在較大差異，工程氣動力方法分析結(jié)果過于保守，易使總體設(shè)計發(fā)生誤判，采取額外的剛度補償，增加結(jié)構(gòu)設(shè)計代價。通過對升力式布局面對稱運載器大攻角擾流流場開展分析，認(rèn)為當(dāng)?shù)貏訅汉蛪毫Ψ植嫉牟町愂菍?dǎo)致工程氣動力和CFD 方法存在差異的主要原因。工程設(shè)計中需要辨識飛行剖面內(nèi)氣動干擾嚴(yán)重的區(qū)域，在此區(qū)域需采用基于CFD 的數(shù)值方法開展精細(xì)的流動特性分析才能對顫振特性開展準(zhǔn)確的預(yù)示與評估。