楊小亮，黃世林，楊 起，柴振霞，孫喜萬(wàn)，劉 偉

(1. 國(guó)防科技大學(xué) 空天科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410073； 2. 中國(guó)人民解放軍93108部隊(duì)， 黑龍江 齊齊哈爾 161007)

三角翼因結(jié)構(gòu)強(qiáng)度高、可提供非線性渦升力等優(yōu)點(diǎn)，是現(xiàn)代戰(zhàn)斗機(jī)主要采用的機(jī)翼構(gòu)型，其氣動(dòng)特性受到了學(xué)界的普遍關(guān)注，尤其是其動(dòng)態(tài)特性密切關(guān)系到現(xiàn)代戰(zhàn)斗機(jī)大攻角機(jī)動(dòng)飛行的安全[1]。大攻角條件下，三角翼會(huì)產(chǎn)生大幅度自維持的翼?yè)u滾現(xiàn)象。雖然早在20世紀(jì)中葉，人們就觀察到了翼?yè)u滾現(xiàn)象，但由于其復(fù)雜的非線性氣動(dòng)特性，直到20世紀(jì)80年代才陸續(xù)出現(xiàn)對(duì)三角翼?yè)u滾現(xiàn)象較為系統(tǒng)的試驗(yàn)研究[2-5]，而耦合求解N-S方程和Euler剛體動(dòng)力學(xué)方程[6-10]完成三角翼?yè)u滾現(xiàn)象的數(shù)值模擬則是二十多年來(lái)才開(kāi)始的事情。

準(zhǔn)確捕捉位于三角翼背風(fēng)面的非線性集中渦是數(shù)值模擬三角翼繞流時(shí)必須解決的問(wèn)題。近年來(lái)，大量文獻(xiàn)在模擬、分析及控制三角翼上渦流特性方面不遺余力[11-15]，進(jìn)展卓有成效[16]，而數(shù)值模擬渦流主導(dǎo)的三角翼?yè)u滾現(xiàn)象則主要集中在捕捉搖滾現(xiàn)象的定性研究。三角翼上集中渦形成和發(fā)展具有強(qiáng)非線性，故對(duì)三角翼構(gòu)型及影響因素極為敏感[17]，造成了數(shù)值模擬與實(shí)驗(yàn)之間、實(shí)驗(yàn)與實(shí)驗(yàn)之間往往存在明顯的分歧[5]。提升數(shù)值模擬和風(fēng)洞實(shí)驗(yàn)反映翼?yè)u滾現(xiàn)象物理本質(zhì)的水平，減小數(shù)值模擬和實(shí)驗(yàn)差距，必須研究翼?yè)u滾特性的影響因素。前緣構(gòu)型和滾轉(zhuǎn)軸安裝位置是實(shí)驗(yàn)研究必須陳述的關(guān)鍵細(xì)節(jié)，且文獻(xiàn)[4-5]表明不同的處理可導(dǎo)致?lián)u滾特性顯著不同。閻超等[18]研究了前緣形狀對(duì)三角翼上渦流特性的影響,但對(duì)搖滾特性影響的研究則相對(duì)較少。本文耦合求解N-S方程和剛體動(dòng)力學(xué)方程，針對(duì)80°后掠三角翼研究前緣構(gòu)型和轉(zhuǎn)軸安裝位置影響三角翼?yè)u滾特性的程度和規(guī)律。

1 模型及計(jì)算方法

1.1 三角翼模型及網(wǎng)格

研究采用80°后掠尖前緣平板三角翼模型，前緣分別采用下削尖、雙面削尖、上削尖三種方式，削尖角均為45°，尾緣不做削尖處理，垂直翼面，如圖1所示。為研究滾轉(zhuǎn)軸對(duì)細(xì)長(zhǎng)三角翼?yè)u滾特性的影響，基于下削尖前緣三角翼分別研究了滾轉(zhuǎn)軸安裝在四種不同位置時(shí)的搖滾特性，如圖2所示，滾轉(zhuǎn)軸位置分別安裝在下翼面下方5 mm處(記為Axel_DD)、下翼面上(記為Axel_D)、上翼面上(記為Axel_U)、上翼面上方5 mm處(記為Axel_UU)，四種安裝位置的滾轉(zhuǎn)軸均位于三角翼體對(duì)稱面上。網(wǎng)格生成采用OH型拓?fù)浣Y(jié)構(gòu)，周向O型、其余方向H型，以平衡描述精度和網(wǎng)格規(guī)模[19]。計(jì)算網(wǎng)格如圖3所示，網(wǎng)格上游邊界距離三角翼頂點(diǎn)2.5倍根弦長(zhǎng)，周向邊界距離體軸2.5倍根弦長(zhǎng)，下游邊界距離三角翼尾緣5倍根弦長(zhǎng)。網(wǎng)格數(shù)為112×226×50(流向×周向×法向)，翼尖、尾緣及翼前緣附近流動(dòng)較為復(fù)雜，局部適當(dāng)加密。

(a) 下削尖前緣(a) Down beveled leading edge

(b) 雙面削尖前緣(b) Double beveled leading edge

(c) 上削尖前緣(c) Up beveled leading edge圖1 三種不同前緣形狀的80°后掠三角翼模型Fig.1 Three 80° swept delta wing model with different leading edge configuration

(a) 滾轉(zhuǎn)軸位于下翼面下方(Axel_DD)(a) Roll-axis under the lower wing surface(Axel_DD)

(b) 滾轉(zhuǎn)軸位于下翼面(Axel_D)(b) Roll-axis on the lower wing surface(Axel_D)

(c) 滾轉(zhuǎn)軸位于上翼面(Axel_U)(c) Roll-axis on the upper wing surface(Axel_U)

(d)滾轉(zhuǎn)軸位于上翼面上方(Axel_UU)(d) Roll-axis above the upper wing surface(Axel_UU)圖2 四種不同滾轉(zhuǎn)軸位置示意圖Fig.2 Schema of four roll-axis installation of delta wing

圖3 80°后掠三角翼網(wǎng)格分布Fig.3 Grid distribution of the 80° swept delta wing

1.2 數(shù)值方法

采用非定常N-S方程組和滾轉(zhuǎn)動(dòng)力學(xué)方程描述三角翼單自由度搖滾運(yùn)動(dòng)[20]，貼體坐標(biāo)系(t,ξ,η,ζ)下的方程寫(xiě)為：

(1)

基于結(jié)構(gòu)網(wǎng)格的有限體積法，采用Roe格式和含雙時(shí)間步的LU-SGS方法離散非定常N-S方程組。滾轉(zhuǎn)動(dòng)力學(xué)方程采用時(shí)間二階精度的單邊差分離散，得到時(shí)間二階精度的滾轉(zhuǎn)差分方程：

(2)

非定常N-S方程組和滾轉(zhuǎn)動(dòng)力學(xué)方程的耦合求解采用松耦合模式，即由滾轉(zhuǎn)動(dòng)力學(xué)方程提供三角翼n+1時(shí)刻的位置，再由非定常N-S方程組計(jì)算n+1時(shí)刻的滾轉(zhuǎn)力矩系數(shù)。

遠(yuǎn)場(chǎng)采用無(wú)反射邊界條件，壁面邊界取絕熱壁，壓力條件為?p/?n|wall=-ρa(bǔ)wall·nwall，速度條件為無(wú)滑移條件，即V=Vwall。

1.3 驗(yàn)證與確認(rèn)

本文的計(jì)算基于國(guó)防科技大學(xué)劉偉教授團(tuán)隊(duì)開(kāi)發(fā)的飛行器動(dòng)態(tài)特性研究程序(Aircraft Dynamic Characteristics Research Program, ADCRP)。該程序針對(duì)多塊對(duì)接結(jié)構(gòu)網(wǎng)格，求解三維非定常RANS方程，經(jīng)過(guò)二十余年的發(fā)展，具備復(fù)雜飛行器流動(dòng)/運(yùn)動(dòng)耦合的多自由度非定過(guò)程的數(shù)值模擬能力。該程序采用系列標(biāo)模外形進(jìn)行了系統(tǒng)的驗(yàn)證與確認(rèn)計(jì)算，適用于三角翼等定、動(dòng)姿態(tài)非定常問(wèn)題的分析[20-21]，在多個(gè)型號(hào)項(xiàng)目中得到了應(yīng)用。該程序具備模擬細(xì)長(zhǎng)三角翼?yè)u滾現(xiàn)象的研究能力。

2 前緣構(gòu)型對(duì)搖滾特性的影響

考察前緣構(gòu)型對(duì)細(xì)長(zhǎng)三角翼?yè)u滾特性的影響。計(jì)算時(shí)，給定來(lái)流馬赫數(shù)Ma=0.2，基于根弦長(zhǎng)度的雷諾數(shù)Re=0.4×106，無(wú)量綱滾轉(zhuǎn)轉(zhuǎn)動(dòng)慣量Ixx=0.1，滾轉(zhuǎn)軸與x軸重合且位于上翼面，不計(jì)軸承阻尼。統(tǒng)一給定10°的初始滾轉(zhuǎn)角，研究不同前緣三角翼的滾轉(zhuǎn)響應(yīng)，分別對(duì)分岔攻角和搖滾振蕩的振幅進(jìn)行討論。

2.1 前緣形狀對(duì)分岔攻角的影響

文獻(xiàn)[22]指出，三角翼攻角α變化會(huì)導(dǎo)致其橫向穩(wěn)定性發(fā)生Hopf分岔，存在攻角的臨界值，攻角小于該值，三角翼橫向動(dòng)態(tài)穩(wěn)定，大于該值，三角翼發(fā)生橫向失穩(wěn)，形成大幅度自維持的極限環(huán)振蕩，該臨界值稱為分岔攻角。

對(duì)分岔攻角的數(shù)值模擬結(jié)果如圖4所示，顯示了下削尖、雙面削尖、上削尖三種前緣構(gòu)型的三角翼不同攻角下的滾轉(zhuǎn)角時(shí)間歷程。按照文獻(xiàn)描述，易于判斷下削尖前緣三角翼分岔攻角約為22.4°，雙面削尖前緣三角翼分岔攻角約為16.9°，上削尖前緣三角翼分岔攻角約為15.5°。其中：下削尖前緣三角翼能夠在較大的攻角范圍保持滾轉(zhuǎn)穩(wěn)定性，分岔攻角最大；上削尖前緣三角翼分岔攻角最小，即最容易發(fā)生橫向失穩(wěn)；雙面削尖前緣的三角翼橫向穩(wěn)定性介于二者之間。

(a) 下削尖前緣(a) Down beveled leading edge

(b) 雙面削尖前緣(b) Double beveled leading edge

(c) 上削尖前緣(c) Up beveled leading edge圖4 不同前緣形狀三角翼的滾轉(zhuǎn)角時(shí)間歷程Fig.4 Time history of roll angle of delta wing with different leading edge

為進(jìn)一步理解前緣構(gòu)型對(duì)三角翼動(dòng)態(tài)特性的影響，基于強(qiáng)迫簡(jiǎn)諧分析法計(jì)算三角翼滾轉(zhuǎn)穩(wěn)定性參數(shù)，如表1和表2所示。穩(wěn)定性參數(shù)為負(fù)則表示具有穩(wěn)定性，為正則表示不穩(wěn)定。由表1可知：不同前緣構(gòu)型的三角翼，在所研究的攻角范圍內(nèi)都是靜穩(wěn)定的，整體而言，靜穩(wěn)定性相當(dāng)。由表2可知：隨著攻角增大，滾轉(zhuǎn)阻尼導(dǎo)數(shù)由負(fù)變正，表明從滾轉(zhuǎn)動(dòng)穩(wěn)定變?yōu)闈L轉(zhuǎn)動(dòng)不穩(wěn)定，從動(dòng)態(tài)穩(wěn)定性的角度說(shuō)明了三角翼滾轉(zhuǎn)方向發(fā)生Hopf分岔的這一現(xiàn)象。其中，下削尖三角翼在22.4°、雙面削尖三角翼在16.9°、上削尖三角翼在15.5°時(shí)，滾轉(zhuǎn)阻尼導(dǎo)數(shù)的絕對(duì)值顯著小于其他角度的絕對(duì)值，這表明此時(shí)三角翼幾乎處于無(wú)阻尼狀態(tài)，正好對(duì)應(yīng)三角翼在跨過(guò)該攻角時(shí)的滾轉(zhuǎn)穩(wěn)定性性態(tài)的變化，分岔攻角在以上角度附近。

表1 三角翼滾轉(zhuǎn)靜穩(wěn)定性導(dǎo)數(shù)

表2 三角翼滾轉(zhuǎn)阻尼導(dǎo)數(shù)

2.2 前緣形狀對(duì)振幅特性的影響

上一節(jié)的研究表明，攻角大于分岔攻角以后，均會(huì)發(fā)生橫向失穩(wěn)，形成大幅度自維持極限環(huán)形式的搖滾現(xiàn)象。為進(jìn)一步考察不同前緣構(gòu)型的三角翼橫向失穩(wěn)以后的動(dòng)態(tài)特性，選取25°和30°攻角，給定0°初始滾轉(zhuǎn)角，研究前緣形狀對(duì)三角翼振幅特性的影響，結(jié)果如圖5所示。

25°攻角時(shí)，如圖5(a)所示，三種前緣構(gòu)型的三角翼均形成了極限環(huán)形式的等幅振蕩，下削尖前緣三角翼的振幅約為19.2°，雙面削尖前緣三角翼的振幅約為65.2°，上削尖前緣三角翼的振幅約為75.6°。其中，下削尖三角翼振幅最小，雙面削尖三角翼次之，上削尖三角翼的振幅最大，對(duì)比表2數(shù)據(jù)，三種前緣構(gòu)型的三角翼的滾轉(zhuǎn)阻尼導(dǎo)數(shù)均為正，從數(shù)值上而言，下削尖三角翼最小，雙面削尖三角翼次之，上削尖三角翼最大，三角翼自激滾轉(zhuǎn)振蕩的表現(xiàn)與動(dòng)態(tài)穩(wěn)定性參數(shù)預(yù)示的規(guī)律一致。

30°攻角時(shí)，如圖5(b)所示，三種機(jī)翼均最終形成了極限環(huán)形式的等幅振蕩，但雙面削尖和上削尖前緣的三角翼在發(fā)生翻轉(zhuǎn)后繞滾轉(zhuǎn)軸進(jìn)行滾轉(zhuǎn)振蕩。其中，雙面削尖前緣三角翼?yè)u滾振幅約為49.8°，上削尖前緣三角翼?yè)u滾振幅約為31.4°，下削尖前緣三角翼未發(fā)生翻轉(zhuǎn)，其振幅約為37.1°。

理論上講，雙面削尖前緣三角翼，上下對(duì)稱，翻轉(zhuǎn)后與自身同構(gòu)，上削尖前緣三角翼與下削尖前緣三角翼在結(jié)構(gòu)上相同，僅與來(lái)流在空間的相對(duì)位置不同。上削尖前緣三角翼翻轉(zhuǎn)后，空間位置與下削尖前緣三角翼相同。但由圖5(b)觀察到的結(jié)果，上削尖前緣三角翼翻轉(zhuǎn)后的振幅(31.4°)與下削尖前緣三角翼的振幅(37.1°)并不相同。經(jīng)分析，主要由于翻轉(zhuǎn)后滾轉(zhuǎn)軸位置不相同，對(duì)于雙面削尖的三角翼同樣存在該問(wèn)題，需要進(jìn)一步分析。

(a) α=25°

(b) α=30° 圖5 不同前緣構(gòu)型三角翼在大攻角下的搖滾特性Fig.5 Wing rock characteristics of delta wing with different leading edge in high angle of attack

3 滾轉(zhuǎn)軸位置對(duì)搖滾特性的影響

基于對(duì)三角翼?yè)u滾現(xiàn)象的理解，進(jìn)一步考察滾轉(zhuǎn)軸位置對(duì)細(xì)長(zhǎng)三角翼?yè)u滾特性的影響。計(jì)算時(shí)，給定來(lái)流馬赫數(shù)Ma=0.2，基于根弦長(zhǎng)度的雷諾數(shù)Re=0.4×106，無(wú)量綱滾轉(zhuǎn)轉(zhuǎn)動(dòng)慣量Ixx=0.1，統(tǒng)一指定10°的初始滾轉(zhuǎn)角，不計(jì)軸承阻尼。選取下削尖三角翼，指定不同滾轉(zhuǎn)軸位置，如圖2所示，研究滾轉(zhuǎn)軸位置對(duì)搖滾特性的影響，分別對(duì)分岔攻角和搖滾振幅進(jìn)行了討論。

3.1 滾轉(zhuǎn)軸位置對(duì)分岔攻角的影響

如圖6(a)所示，滾轉(zhuǎn)軸位于下翼面下方5 mm處時(shí)，三角翼從10°初始滾轉(zhuǎn)角釋放后，在25°攻角下?lián)u滾振幅呈衰減狀態(tài)，在30°攻角下?lián)u滾振幅呈增大狀態(tài)，而在27°攻角時(shí)搖滾振幅基本保持10°不變，說(shuō)明滾轉(zhuǎn)軸位于下翼面下方5 mm的三角翼的分岔攻角約為27°。由圖6(b)～(d)易知，滾轉(zhuǎn)軸位于下翼面的三角翼分岔攻角為26.5°，滾轉(zhuǎn)軸位于上翼面的三角翼分岔攻角為22.4°，滾轉(zhuǎn)軸位于上翼面上方5 mm處的三角翼的分岔攻角為17.4°。由此可見(jiàn)，隨著滾轉(zhuǎn)軸從下翼面下方5 mm處移動(dòng)到上翼面上方5 mm處，三角翼的分岔攻角逐漸減小?？梢缘贸鼋Y(jié)論，在所研究的滾轉(zhuǎn)軸位置范圍內(nèi)，滾轉(zhuǎn)軸位置越高，三角翼?yè)u滾的分岔攻角越小，越易于橫向失穩(wěn)，反之亦然。

(a) 滾轉(zhuǎn)軸位于下翼面以下(Axel_DD)(a) Roll-axis under the lower wing surface(Axel_DD)

(b) 滾轉(zhuǎn)軸位于下翼面(Axel_D)(b) Roll-axis on the lower wing surface(Axel_D)

(c) 滾轉(zhuǎn)軸位于上翼面(Axel_U)(c) Roll-axis on the upper wing surface(Axel_U)

(d) 滾轉(zhuǎn)軸位于上翼面以上(Axel_UU)(d) Roll-axis above the upper wing surface(Axel_UU)圖6 不同滾轉(zhuǎn)軸位置的三角翼滾轉(zhuǎn)角時(shí)間歷程曲線Fig.6 Roll angle time historical solution of delta wing installed at different roll-axes

3.2 滾轉(zhuǎn)軸位置對(duì)振幅特性的影響

進(jìn)一步研究滾轉(zhuǎn)軸位置對(duì)三角翼振幅特性的影響，給定30°、40°、50°和60°攻角，模擬四種滾轉(zhuǎn)軸位置的三角翼?yè)u滾振幅特性，計(jì)算結(jié)果如圖7所示。

從圖7中可以看出，四種滾轉(zhuǎn)軸位置下?lián)u滾振幅隨攻角變化的規(guī)律基本一致，30°攻角時(shí)，四種滾轉(zhuǎn)軸位置的三角翼均形成了極限環(huán)自激搖滾，滾轉(zhuǎn)軸位置越靠上振幅越大；40°攻角左右建立最大振幅，隨著攻角進(jìn)一步增大，由于受到渦破裂的影響，搖滾振幅隨攻角增大而減??；當(dāng)攻角增加到60°以后，四種滾轉(zhuǎn)軸位置的三角翼均難以形成大幅度的搖滾振幅。

圖7 搖滾振幅隨攻角變化曲線Fig.7 Amplitude of roll varies with attack angle

4 結(jié)論

耦合求解三維非定常N-S方程和歐拉剛體動(dòng)力學(xué)方程組，研究了前緣構(gòu)型與滾轉(zhuǎn)軸位置對(duì)80°后掠平板三角翼?yè)u滾特性的影響，得到以下結(jié)論：

1)前緣削尖方式不改變?nèi)且頇M向動(dòng)態(tài)穩(wěn)定性的性態(tài)，但影響三角翼?yè)u滾的分岔攻角和搖滾振幅；不同前緣構(gòu)型三角翼的自激搖滾特性與動(dòng)態(tài)穩(wěn)定性參數(shù)表征的穩(wěn)定性一致。

2)對(duì)于轉(zhuǎn)軸位于上表面的三角翼，攻角25°時(shí),前緣下削尖時(shí)搖滾振幅最小，上削尖時(shí)搖滾振幅最大；攻角30°時(shí)，雙面削尖前緣和上削尖前緣的三角翼自激搖滾時(shí)發(fā)生了翻轉(zhuǎn)，翻轉(zhuǎn)后繞滾轉(zhuǎn)軸等幅振蕩，雙面削尖前緣與上削尖前緣三角翼先翻轉(zhuǎn)后再形成極限環(huán)形式的等幅振蕩。

3)滾轉(zhuǎn)軸的安裝位置影響三角翼的搖滾特性，滾轉(zhuǎn)軸位置越高，三角翼?yè)u滾的分岔攻角越小，極限環(huán)滾轉(zhuǎn)振蕩的振幅越大，反之亦然。

三角翼?yè)u滾現(xiàn)象研究多年依然是熱點(diǎn)問(wèn)題，對(duì)其產(chǎn)生和維持機(jī)理已有較多深刻的認(rèn)識(shí)，但對(duì)敏感因素的分析還不足，旋渦等流動(dòng)結(jié)構(gòu)與飛行器之間作用的過(guò)程和機(jī)理等問(wèn)題仍值得進(jìn)一步的分析。