楊小亮，劉偉，趙云飛，劉 君

(國(guó)防科學(xué)技術(shù)大學(xué)航天與材料工程學(xué)院宇航科學(xué)與工程系，湖南 長(zhǎng)沙 410073)

0 引言

戰(zhàn)斗機(jī)大攻角條件下的動(dòng)態(tài)運(yùn)動(dòng)特性，一直是國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)之一，近幾十年來開展的大量研究主要采用了以下兩類思路:一類是不考慮滾轉(zhuǎn)運(yùn)動(dòng)的單自由度快速上仰或是強(qiáng)迫俯仰振動(dòng)，另一類是不考慮動(dòng)態(tài)俯仰角影響的單自由度自激搖滾。一般認(rèn)為，快速上仰運(yùn)動(dòng)能夠推遲戰(zhàn)斗機(jī)的失速攻角和提高戰(zhàn)斗機(jī)的最大升力系數(shù)［1］(如Su-27戰(zhàn)斗機(jī)的眼鏡蛇機(jī)動(dòng)Cobra Maneuver)，國(guó)外 Hoang［2］，Visbal［3］等以及國(guó)內(nèi)的白濤［4］、楊立芝［5］等都針對(duì)三角翼的俯仰運(yùn)動(dòng)進(jìn)行了研究。而在單自由度搖滾方面，Nguyen等1981年的一篇報(bào)告［6］，被公認(rèn)是對(duì)三角翼自激搖滾問題最早的實(shí)驗(yàn)研究，后來，Kandil領(lǐng)導(dǎo)的研究小組完成了三角翼自激搖滾的數(shù)值模擬［7］，在國(guó)內(nèi)，楊希明［8］、唐敏中［9］等人采用實(shí)驗(yàn)的方法研究了三角翼的自激搖滾，劉偉［10］研究了細(xì)長(zhǎng)機(jī)翼自激搖滾的數(shù)值模擬方法，張涵信［11］從非線性動(dòng)力學(xué)的觀點(diǎn)給出了細(xì)長(zhǎng)機(jī)翼產(chǎn)生自激搖滾的解析條件。然而真實(shí)的戰(zhàn)斗機(jī)上仰機(jī)動(dòng)飛行，往往導(dǎo)致橫向失穩(wěn)，發(fā)生非可控的自激滾轉(zhuǎn)，形成俯仰、滾轉(zhuǎn)耦合的運(yùn)動(dòng)，縱橫向運(yùn)動(dòng)交感耦合和運(yùn)動(dòng)學(xué)關(guān)聯(lián)效應(yīng)將會(huì)導(dǎo)致氣動(dòng)力/力矩的強(qiáng)非線性，使得戰(zhàn)斗機(jī)操穩(wěn)性能下降，危及飛行安全。目前，采用實(shí)驗(yàn)方法針對(duì)俯仰、滾轉(zhuǎn)耦合運(yùn)動(dòng)的研究［12］還不多，而對(duì)三角翼俯仰、滾轉(zhuǎn)雙自由度耦合運(yùn)動(dòng)的數(shù)值研究也僅僅局限于給定強(qiáng)迫俯仰、強(qiáng)迫滾轉(zhuǎn)運(yùn)動(dòng)條件下的模擬［13］，這樣的純強(qiáng)迫運(yùn)動(dòng)雖然能夠從一定程度上反映戰(zhàn)斗機(jī)氣動(dòng)力/力矩的遲滯特性，但對(duì)縱橫向運(yùn)動(dòng)交感耦合的模擬有其無法回避的局限性。

本文耦合求解Navier-Stokes方程和俯仰/滾轉(zhuǎn)運(yùn)動(dòng)方程研究80°三角翼在強(qiáng)迫俯仰、自由滾轉(zhuǎn)雙自由度條件下的運(yùn)動(dòng)特性，得到了與試驗(yàn)定性相符的結(jié)果，分析了俯仰/滾轉(zhuǎn)雙自由度耦合運(yùn)動(dòng)與單自由度滾轉(zhuǎn)及俯仰運(yùn)動(dòng)的差異，研究了多自由度耦合運(yùn)動(dòng)條件下自激滾轉(zhuǎn)運(yùn)動(dòng)產(chǎn)生的機(jī)理，給出了俯仰/滾轉(zhuǎn)耦合運(yùn)動(dòng)條件下三角翼的氣動(dòng)力/力矩特性。

1 模型及數(shù)值方法

1.1 三角翼模型

本文計(jì)算使用的80°后掠三角翼模型如圖1所示，模型的前緣及尾緣采用單面削尖處理。文獻(xiàn)［14］通過計(jì)算比較，認(rèn)為H型網(wǎng)格在模擬的準(zhǔn)確性上具有優(yōu)勢(shì)，O型網(wǎng)格能夠有效地節(jié)省周向網(wǎng)格，因此本文采用周向O型、其余方向H型的OH型網(wǎng)格進(jìn)行計(jì)算，翼面網(wǎng)格分布如圖1所示。

圖1 80°后掠三角翼模型Fig.1 80°delta wing model

1.2 控制方程

機(jī)翼的強(qiáng)迫俯仰運(yùn)動(dòng)給定正弦規(guī)律:

假定俯仰軸、滾轉(zhuǎn)軸均為慣量主軸，則三角翼滾轉(zhuǎn)運(yùn)動(dòng)方程為:

無量綱化(2)式并進(jìn)行簡(jiǎn)化得:

其中Q為守恒變量，Ek、Evk分別為對(duì)流通量和粘性通量，ξ、η、ζ為貼體坐標(biāo)系坐標(biāo)。

1.3 數(shù)值離散方法

采用空間二階NND格式和含雙時(shí)間步法的LUSGS方法有限體積離散NS方程。剛體滾轉(zhuǎn)運(yùn)動(dòng)方程采用時(shí)間二階精度的單邊差分離散、，得到時(shí)間二階精度的差分方程:

采用松耦合的方式聯(lián)立求解流動(dòng)控制方程和剛體動(dòng)力學(xué)方程，即由剛體動(dòng)力學(xué)方程提供三角翼“n+1”時(shí)刻的位置，再由流動(dòng)控制方程計(jì)算“n+1”時(shí)刻的滾轉(zhuǎn)力矩系數(shù)。

計(jì)算域遠(yuǎn)場(chǎng)采用無反射邊界條件，壁面邊界取絕熱壁假設(shè)，壓力條件為:，速度條件為考慮物面雙自由度運(yùn)動(dòng)的無滑移條件:

2 三角翼單自由度自激滾轉(zhuǎn)數(shù)值模擬

一些文獻(xiàn)完成了80°后掠三角翼單自由度自激滾轉(zhuǎn)的實(shí)驗(yàn)［8，9，15］和數(shù)值計(jì)算［10］，但由于模型的差異以及與實(shí)驗(yàn)條件不能完全一致，如前緣形狀、模型厚度、來流馬赫數(shù)、雷諾數(shù)等，特別是不同軸承摩擦阻尼的差異，導(dǎo)致研究得到的自激滾轉(zhuǎn)角振幅、頻率有一定差別，但定性結(jié)論基本一致。

本文首先計(jì)算了80°后掠三角翼的單自由度自激滾轉(zhuǎn)，根弦名義攻角30°，來流馬赫數(shù)Ma=0.35，根弦雷諾數(shù) Re=1.0 ×106，無量綱轉(zhuǎn)動(dòng)慣量 c1=0.65，不計(jì)機(jī)械阻尼，所得運(yùn)動(dòng)特性如圖2所示。振幅從0°逐漸增大，形成以0°為平均滾轉(zhuǎn)角、振幅38°左右極限環(huán)形式蕩，曲線振幅、頻率穩(wěn)定，與相關(guān)文獻(xiàn)的研究結(jié)果［10，15］定性相符。

圖2 單自由度自由滾轉(zhuǎn)運(yùn)動(dòng)特性Fig.2 Characteristics of self-excited roll motion的 振

3 三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)雙自由度運(yùn)動(dòng)數(shù)值研究

對(duì)80°后掠三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)耦合運(yùn)動(dòng)的模擬采用圖1所示的網(wǎng)格，給定來流 Ma=0.35，Re=1.0 ×106，俯仰減縮頻率 k=0.005，研究俯仰角在 0°～40°之間以正弦規(guī)律運(yùn)動(dòng)時(shí)的自激滾轉(zhuǎn)特性。給定無量綱轉(zhuǎn)動(dòng)慣量c1=0.65，機(jī)構(gòu)阻尼c2=0，初始俯仰角θ0=20°，俯仰振幅 θm=20°，以及 0°初始滾轉(zhuǎn)角，模擬得到80°后掠三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)雙自由度耦合運(yùn)動(dòng)的姿態(tài)角時(shí)間歷程如圖3所示。在動(dòng)態(tài)俯仰角θ的影響下，滾轉(zhuǎn)角時(shí)間歷程曲線的振幅和頻率也隨之發(fā)生變化，大致呈現(xiàn)為:在俯仰運(yùn)動(dòng)的正行程(＞0)，由于攻角總體增大，自由滾轉(zhuǎn)的振幅增大、頻率加快;在俯仰運(yùn)動(dòng)的負(fù)行程(＜0)，由于攻角總體減小，自由滾轉(zhuǎn)的振幅減小、頻率變緩。與單自由度的自由滾轉(zhuǎn)(如圖2所示)比較發(fā)現(xiàn)，在雙自由度運(yùn)動(dòng)達(dá)到最大滾轉(zhuǎn)振幅時(shí)的頻率與單自由度的自激搖滾頻率相當(dāng)，而在俯仰振動(dòng)曲線的波谷區(qū)域，頻率明顯比單自由度自激滾轉(zhuǎn)的頻率低。在俯仰運(yùn)動(dòng)的正行程，機(jī)翼向上的速度分量導(dǎo)致等效攻角減小，延遲了渦破裂在翼面上發(fā)生;在俯仰運(yùn)動(dòng)的負(fù)行程，機(jī)翼向下的速度分量導(dǎo)致等效攻角增大，推遲了破裂渦流在翼面上的再附。

機(jī)翼俯仰運(yùn)動(dòng)速度對(duì)翼面上流動(dòng)結(jié)構(gòu)的形成具有延遲作用，導(dǎo)致最大振幅形成的位置滯后于最大俯仰角所在的位置(如圖3所示)。唐敏中等人［12］采用實(shí)驗(yàn)的方法，研究了80°三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)的運(yùn)動(dòng)特性，給定來流速度30m/s，俯仰頻率0.088Hz，折合無量綱俯仰減縮頻率約為0.007，得到的運(yùn)動(dòng)特性如圖4所示，其滾轉(zhuǎn)振幅滯后于俯仰角變化的效應(yīng)更為明顯。圖3(計(jì)算結(jié)果)顯示了與圖4(實(shí)驗(yàn)結(jié)果)定性一致的運(yùn)動(dòng)特性:隨著俯仰角的變化，滾轉(zhuǎn)角的振幅和頻率經(jīng)歷了增大、穩(wěn)定再減小的過程，滾轉(zhuǎn)振幅總體隨俯仰角的增大而增大，減小而減小。實(shí)驗(yàn)結(jié)果中值得注意的是，當(dāng)俯仰角回到零附近的區(qū)域時(shí)，機(jī)翼運(yùn)動(dòng)到約－17°滾轉(zhuǎn)角的位置突然停住，直至俯仰角再增大到5°左右時(shí)，又迅速起動(dòng)向0°滾轉(zhuǎn)角的位置運(yùn)動(dòng)，這可能和轉(zhuǎn)動(dòng)慣量、軸承摩擦阻尼以及翼面上的流態(tài)等因素有關(guān)，對(duì)該現(xiàn)象的進(jìn)一步研究，作者另文闡述。

圖3分別用A、B、C標(biāo)識(shí)了80°后掠三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)耦合運(yùn)動(dòng)過程中俯仰運(yùn)動(dòng)第二周期中三個(gè)不同的滾轉(zhuǎn)周期，分別代表振幅增大、振幅極值和振幅衰減時(shí)運(yùn)動(dòng)特性，其滾轉(zhuǎn)力矩系數(shù)隨滾轉(zhuǎn)角的變化如圖5所示。在周期A的過程中(曲線A)，曲線呈順時(shí)針變化，說明機(jī)翼滾轉(zhuǎn)從流動(dòng)中汲取能量，因此振幅增大;在周期B的過程中(曲線B)，雙“8”字型曲線中部的內(nèi)環(huán)呈順時(shí)針變化，表征機(jī)翼從流動(dòng)中汲取能量，在曲線兩端的外環(huán)呈逆時(shí)針變化，表征機(jī)翼滾轉(zhuǎn)向流動(dòng)釋放能量，兩部分面積大致相當(dāng)，說明機(jī)翼在滾轉(zhuǎn)振幅極值附近機(jī)翼與流動(dòng)的能量交換大致平衡;在周期C的過程中(曲線C)，曲線呈逆時(shí)針變化，機(jī)翼滾轉(zhuǎn)向流動(dòng)釋放能量，振幅衰減。

圖6(a)、(b)分別是通過數(shù)值模擬和實(shí)驗(yàn)［12］得到的三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)雙自由度耦合運(yùn)動(dòng)中滾轉(zhuǎn)角的相平面圖。計(jì)算的結(jié)果和文獻(xiàn)［12］通過實(shí)驗(yàn)得到的結(jié)果定性一致，圖中縱橫交叉的相曲線顯示了三角翼雙自由度運(yùn)動(dòng)時(shí)，動(dòng)態(tài)俯仰角對(duì)三角翼自由滾轉(zhuǎn)運(yùn)動(dòng)特性的影響。

圖3 強(qiáng)迫俯仰、自由滾轉(zhuǎn)姿態(tài)角及滾轉(zhuǎn)力矩系數(shù)時(shí)間歷程Fig.3 Attitude angles and roll-moment coefficient vs.time in force-pitch and free-roll

圖7 俯仰/滾轉(zhuǎn)運(yùn)動(dòng)的氣動(dòng)力特性Fig.7 Aerodynamic characteristics of pitch/roll motion

三角翼自激滾轉(zhuǎn)運(yùn)動(dòng)的形成和維持主要受到滾轉(zhuǎn)力矩的影響，80°后掠三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)運(yùn)動(dòng)的滾轉(zhuǎn)力矩系數(shù)如圖3所示:起振的階段滾轉(zhuǎn)角在與其同號(hào)的滾轉(zhuǎn)力矩的作用下逐漸增大;在振幅相對(duì)穩(wěn)定的階段，滾轉(zhuǎn)力矩系數(shù)基本維持與機(jī)翼的自由滾轉(zhuǎn)同頻率的振蕩;當(dāng)俯仰角減小以后，由于俯仰角減小導(dǎo)致等效攻角顯著減小，流動(dòng)提供的滾轉(zhuǎn)力矩不足以維持機(jī)翼極限環(huán)形式的振蕩，振幅便開始減小;隨著俯仰角的進(jìn)一步減小，當(dāng)俯仰角趨近零時(shí)，由于等效攻角減小，滾轉(zhuǎn)力矩系數(shù)也趨于零，因此機(jī)翼在滾轉(zhuǎn)方向幾乎處于勻速滾轉(zhuǎn)狀態(tài)，滾轉(zhuǎn)角緩慢轉(zhuǎn)動(dòng)到較大的角度;當(dāng)俯仰角又從零開始逐漸增大，引起了等效攻角增大，使得起恢復(fù)作用的滾轉(zhuǎn)力矩迅速變大，將機(jī)翼很快拉回到零滾轉(zhuǎn)角附近，并進(jìn)一步增大，進(jìn)入了強(qiáng)迫俯仰運(yùn)動(dòng)下一周期的振蕩，如此循環(huán)形成了80°三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)的周期性運(yùn)動(dòng)特性。

圖7是80°后掠三角翼單自由度強(qiáng)迫俯仰氣動(dòng)力特性與俯仰/滾轉(zhuǎn)耦合雙自由度運(yùn)動(dòng)氣動(dòng)力特性的比較(曲線A是單自由度俯仰結(jié)果)。三角翼俯仰、滾轉(zhuǎn)雙自由度運(yùn)動(dòng)時(shí)，攻角的變化主要受俯仰角影響，如圖7(a)所示，總體變化趨勢(shì)是隨俯仰角的增大而增大，減小而減小，而自由滾轉(zhuǎn)運(yùn)動(dòng)的影響使得攻角減小，也就是說，攻角的變化是俯仰、滾轉(zhuǎn)運(yùn)動(dòng)共同作用的結(jié)果。升力主要受攻角的影響，升力的變化總體和攻角的變化一致，因此俯仰運(yùn)動(dòng)的過程中伴隨的滾轉(zhuǎn)運(yùn)動(dòng)將會(huì)導(dǎo)致飛行器升力顯著減小，如圖7(b)所示，考慮自由滾轉(zhuǎn)影響下的機(jī)翼俯仰運(yùn)動(dòng)，與單自由度俯仰運(yùn)動(dòng)相比，上仰引起平均失速攻角增大的幅度要小，能夠達(dá)到的最大升力系數(shù)也要小，這對(duì)于機(jī)動(dòng)飛行的安全具有重要影響，必須引起足夠的重視。圖7(c)、(d)分別是側(cè)力和滾轉(zhuǎn)力矩的變化，與單自由度的俯仰運(yùn)動(dòng)相比，俯仰滾轉(zhuǎn)耦合運(yùn)動(dòng)產(chǎn)生了顯著不等于零的側(cè)力和滾轉(zhuǎn)力矩，說明在滾轉(zhuǎn)運(yùn)動(dòng)的影響下，機(jī)翼的橫向穩(wěn)定性變差，側(cè)力的存在還會(huì)導(dǎo)致機(jī)翼的側(cè)滑運(yùn)動(dòng)，應(yīng)進(jìn)一步研究。

4 結(jié)論

通過耦合求解NS方程和滾轉(zhuǎn)運(yùn)動(dòng)方程研究了80°三角翼強(qiáng)迫俯仰、自由滾轉(zhuǎn)雙自由度運(yùn)動(dòng)特性，得到了與試驗(yàn)定性一致的結(jié)果，分析了雙自由度耦合運(yùn)動(dòng)與單自由度滾轉(zhuǎn)運(yùn)動(dòng)的差異以及一些運(yùn)動(dòng)現(xiàn)象產(chǎn)生的機(jī)理，并研究了雙自由度耦合運(yùn)動(dòng)的氣動(dòng)力特性，得到了以下結(jié)論:

(1)機(jī)翼的俯仰運(yùn)動(dòng)速度對(duì)機(jī)翼自由滾轉(zhuǎn)運(yùn)動(dòng)的振幅具有遲滯作用，使得滾轉(zhuǎn)振幅的變化落后于俯仰角的變化;

(2)在動(dòng)態(tài)俯仰角的影響下，三角翼的自由滾轉(zhuǎn)不能形成極限環(huán)形式的自由滾轉(zhuǎn)振蕩，滾轉(zhuǎn)運(yùn)動(dòng)的振幅大致隨俯仰角的增大而增大，減小而減小;

(3)三角翼的俯仰、滾轉(zhuǎn)雙自由度耦合運(yùn)動(dòng)，在俯仰運(yùn)動(dòng)的波谷區(qū)域，起恢復(fù)作用的氣動(dòng)力矩較小，在不計(jì)軸承機(jī)械阻尼的情況下，機(jī)翼勻速滾轉(zhuǎn)至較大滾轉(zhuǎn)角的位置，與試驗(yàn)結(jié)果有明顯差異，應(yīng)進(jìn)一步考慮轉(zhuǎn)動(dòng)慣量、軸承摩擦以及湍流等多種因素的影響;

(4)俯仰耦合滾轉(zhuǎn)運(yùn)動(dòng)將導(dǎo)致升力顯著降低，對(duì)機(jī)動(dòng)飛行安全不利，同時(shí)俯仰、滾轉(zhuǎn)運(yùn)動(dòng)的耦合導(dǎo)致機(jī)翼的橫向穩(wěn)定性變差，真實(shí)的機(jī)動(dòng)飛行可能還耦合了包括側(cè)滑在內(nèi)的其他形式的運(yùn)動(dòng)，應(yīng)進(jìn)一步研究。

致 謝:張涵信院士對(duì)本文的工作提出了許多寶貴的建議，在此表示衷心的感謝。

［1］ERICSSON L E.Delta wing vortex dynamics in high rate oscillations in pitch［R］.AIAA paper 97－0327，1997.

［2］HOANG N T，REDINIOTIS O K，TELIONIS D P.3-D LDV measurements over a delta wing in pitch-up motion.［R］.AIAA paper 93－0185，1993.

［3］VISBAL M R.Computational study of vortex breakdown on a pitching delta wing［R］.AIAA paper 93 －2974，1993.

［4］白濤，呂志詠，楊曉峰.動(dòng)態(tài)三角翼的非定常氣動(dòng)特性和壓力分布相關(guān)研究［J］.流體力學(xué)實(shí)驗(yàn)與測(cè)量，2004，18(2):38－42.

［5］楊立芝，高正紅.繞三角翼縱向俯仰大迎角氣動(dòng)特性計(jì)算研究［J］.航空學(xué)報(bào)，2003，24(5):413－416.

［6］NGUYEN L E，YIN L P，CHAMBERS J R.Self-induced wing rock of slender wings［R］.AIAA paper 81 － 1883，1981.

［7］KANDIL O A，MENZIES M A.Effective control of computationally simulated wing rock in subsonic flow［R］.AIAA paper 97－0831，1997.

［8］楊希明，程厚梅，張祖庚.大攻角80°前緣后掠平板三角翼?yè)u滾實(shí)驗(yàn)與數(shù)值模擬［A］.空氣動(dòng)力學(xué)研究文集(第四卷)［C］.北京，1994，142 －148.

［9］唐敏中，王鐵，張偉，徐忠利.低速三角翼滾擺試驗(yàn)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1997，15(4):436 －443.

［10］劉偉.細(xì)長(zhǎng)機(jī)翼?yè)u滾機(jī)理的非線性動(dòng)力學(xué)分析及數(shù)值模擬方法研究［D］.［博士學(xué)位論文］，國(guó)防科學(xué)技術(shù)大學(xué)研究生院，2004.

［11］張涵信，劉偉，謝昱飛，葉友達(dá).后掠三角翼的搖滾及其動(dòng)態(tài)演化問題［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2006，24(1):5 －9.

［12］唐敏中，張偉，何宏麗.俯仰、滾擺耦合復(fù)雜流場(chǎng)試驗(yàn)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19(1):47 －55.

［13］KANDIL O A，MENZIES M A.Coupled rolling and pitching oscillation effects on transonic shock-induced vortex-breakdown flow of a delta wing［R］.AIAA paper 96 －0828，1996.

［14］劉剛，周鑄，黃勇，陳作斌.三角翼大迎角繞流數(shù)值模擬中網(wǎng)格的影響研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2004，22(4):481－485.

［15］ALAIN PELLETIER.A study of the nonlinear aerodyanmic characteristics of a slender double-delta wing in roll［D］.［Ph.D.thesis］.Notre Dame University，April 1998.