蔡巧言 呂計(jì)男 郭力 王昕江

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

0 引言

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

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

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

1 顫振分析方法

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

1.1 工程氣動(dòng)力方法

時(shí)域顫振計(jì)算的控制方程可以表示為式（1）所示

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

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

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

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

1.2 CFD/CSD 耦合分析方法

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

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

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

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

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

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

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

3 顫振分析

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

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

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

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

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

3.2 基于CFD 方法的顫振分析

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

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

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

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

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

3.3 顫振特性差異分析

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

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

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

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

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

不同馬赫數(shù)、飛行攻角前體對V 尾的干擾特性不同，V 尾的當(dāng)?shù)貏?dòng)壓/壓力分布變化規(guī)律也不同。采用固定馬赫數(shù)（Ma2.3）和固定來流動(dòng)壓（22.2kPa）變換攻角（8°/18°/28°）的方式研究攻角的影響（圖9）。整體規(guī)律上看，在所計(jì)算攻角范圍內(nèi)，V 尾都會(huì)發(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)工程氣動(dòng)力方法與CFD 方法預(yù)示顫振動(dòng)壓產(chǎn)生較大偏差的直接原因是兩種分析方法得到的壓力分布和當(dāng)?shù)貏?dòng)壓不一致導(dǎo)致，工程方法沒有反映出機(jī)身波系及大攻角分離對舵面前方來流流場的干擾影響。

4 結(jié)論

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