肖中云，劉剛，江雄，王建濤

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

飛行器雙自由度耦合振蕩的數(shù)值模擬方法

肖中云*，劉剛，江雄，王建濤

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

針對模擬雙自由度(2DOF)運(yùn)動的“PQR”強(qiáng)迫振蕩試驗裝置，通過坐標(biāo)變換，將機(jī)構(gòu)運(yùn)動轉(zhuǎn)化為空間歐拉角進(jìn)行表示，實現(xiàn)對機(jī)構(gòu)俯仰/偏航、俯仰/滾轉(zhuǎn)和偏航/滾轉(zhuǎn)等3種耦合運(yùn)動的統(tǒng)一描述，采用預(yù)定運(yùn)動軌跡的動網(wǎng)格計算技術(shù)，構(gòu)造了基于非定常雷諾平均Navier-Stokes (URANS)方程的飛行器雙自由度振蕩數(shù)值模擬方法。圍繞某復(fù)雜構(gòu)型飛行器低速流動，在模型迎角10°、滾轉(zhuǎn)角±40°、偏航角±40°振蕩預(yù)設(shè)條件下分別對單自由度偏航、滾轉(zhuǎn)振蕩和雙自由度偏航/滾轉(zhuǎn)振蕩進(jìn)行了模擬，得到的氣動系數(shù)遲滯回線在形態(tài)和量值上與風(fēng)洞試驗結(jié)果吻合，證明當(dāng)前方法可以作為飛行器復(fù)雜耦合運(yùn)動的一種有效研究手段。針對偏航/滾轉(zhuǎn)耦合振蕩，從運(yùn)動形式和氣動特性上與單獨(dú)偏航、滾轉(zhuǎn)運(yùn)動進(jìn)行了對比分析，結(jié)果表明，耦合運(yùn)動在氣流角、繞體軸的旋轉(zhuǎn)角速度等方面與單自由度疊加效果不同，力矩系數(shù)遲滯回線存在曲線交叉的現(xiàn)象，表現(xiàn)出與單自由度振蕩不同的阻尼特性。在當(dāng)前模擬狀態(tài)下，偏航/滾轉(zhuǎn)振蕩流場以橫向分離渦運(yùn)動為主要特征，可以推測大迎角情況下縱橫向分離渦流動結(jié)構(gòu)更加復(fù)雜，耦合作用更強(qiáng)，需要進(jìn)一步在分離流模擬方法上開展研究。

雙自由度；大幅振蕩；耦合運(yùn)動；非定常；數(shù)值模擬

飛行器過失速機(jī)動是一種大迎角下的多自由度耦合運(yùn)動，氣動力呈現(xiàn)出嚴(yán)重的非對稱、非線性和非定常特性，由各種實驗數(shù)據(jù)線性疊加形成總的氣動系數(shù)的方法，理論上缺少依據(jù)，實際上會存在很大偏差。圍繞過失速機(jī)動問題，世界各航空發(fā)達(dá)國家都開發(fā)了相應(yīng)的風(fēng)洞試驗與理論分析方法[1-5]，目的是為了獲取飛行器在大迎角過失速區(qū)域的非線性、非定常氣動力特性，有助于更好理解動態(tài)空氣動力的復(fù)雜現(xiàn)象，以及幫助發(fā)展更好的飛行力學(xué)數(shù)學(xué)模型。

隨著對飛行器過失速機(jī)動能力的日益重視，國內(nèi)外發(fā)展了大量的風(fēng)洞動態(tài)試驗技術(shù)[6-8]，包括大迎角動導(dǎo)數(shù)、旋轉(zhuǎn)天平、旋轉(zhuǎn)振蕩、大幅振蕩試驗技術(shù)等?！癙QR”動態(tài)試驗系統(tǒng)是一種典型的大振幅雙自由度試驗裝置，可進(jìn)行大幅俯仰、偏航、滾轉(zhuǎn)振蕩及其雙自由度耦合振蕩、大幅俯仰拉升及俯仰和滾轉(zhuǎn)耦合拉升試驗[9]。大幅振蕩試驗是進(jìn)行飛行機(jī)動模擬，特別是過失速機(jī)動模擬的基本試驗數(shù)據(jù)。目前單自由度振蕩(俯仰、偏航和滾轉(zhuǎn))已經(jīng)開展了較多的研究，在獲取動導(dǎo)數(shù)、動態(tài)非定常時間遲滯效應(yīng)上十分有效，同時可以用于模擬簡單的機(jī)動飛行動作，如眼鏡蛇機(jī)動和機(jī)翼搖滾等。對于復(fù)雜機(jī)動飛行動作來說，運(yùn)動通常為多個自由度的耦合運(yùn)動，如Herbst機(jī)動，飛行器在大迎角下同時包含了偏航和滾轉(zhuǎn)的耦合運(yùn)動。這就要求不斷提高試驗?zāi)M的運(yùn)動相似性，從單自由度向雙自由度耦合、多自由度機(jī)動的方向發(fā)展。

雙自由度振蕩是一種復(fù)合運(yùn)動，運(yùn)動參數(shù)包括平均角度、振幅、振蕩頻率、相位等，各參數(shù)的組合情況非常多，運(yùn)動規(guī)律及氣動特性十分復(fù)雜。文獻(xiàn)[10]較早開展了三角翼的俯仰/滾轉(zhuǎn)耦合運(yùn)動試驗研究，后續(xù)又開展了典型飛機(jī)的偏航/滾轉(zhuǎn)耦合氣動力試驗[11-12]，模型支撐迎角達(dá)到60°，將耦合運(yùn)動時的氣動特性與單自由度運(yùn)動疊加得到的氣動特性進(jìn)行了對比。文獻(xiàn)[13]介紹了生產(chǎn)型低速風(fēng)洞的雙自由度大幅振蕩試驗技術(shù)，在該技術(shù)支撐下開展了某復(fù)雜構(gòu)型飛機(jī)偏航/滾轉(zhuǎn)耦合運(yùn)動非定常氣動力特性試驗研究，獲得了重復(fù)性好、可靠性高的試驗數(shù)據(jù)[14]。以上這些試驗為研究飛行器大迎角氣動特性積累了寶貴經(jīng)驗，同時也為非定常氣動建模的發(fā)展提供了有力支撐[15-17]。

目前，用CFD方法模擬飛行器多自由耦合運(yùn)動的工作還比較少見。文獻(xiàn)[18]對二維平板機(jī)翼前沖、俯仰振蕩以及兩者的合成運(yùn)動進(jìn)行了模擬，研究大尺度渦結(jié)構(gòu)對氣動的影響。文獻(xiàn)[19]研究了低雷諾數(shù)下的機(jī)翼俯仰和下沉運(yùn)動，以此模擬微型飛行器的非定?？諝鈩恿π?yīng)。隨著計算機(jī)運(yùn)算速度的提高，非定常計算可采用的計算規(guī)模大幅提高，對分離渦流動、復(fù)雜湍流的模擬能力日益增強(qiáng)。發(fā)展多自由度大幅振蕩的數(shù)值模擬方法，有助于摸清復(fù)雜氣動特性背后的流動機(jī)理，對風(fēng)洞試驗起到相互驗證與校核的作用。本文基于描述姿態(tài)運(yùn)動的歐拉角方法和動網(wǎng)格數(shù)值計算，發(fā)展了模擬飛行器雙自由度大幅振蕩的數(shù)值模擬方法，探討了耦合運(yùn)動的特點(diǎn)及其與解耦運(yùn)動的關(guān)系，最后以偏航/滾轉(zhuǎn)振蕩為例，通過對某復(fù)雜外形飛行器的模擬驗證了本文方法。

1 “PQR”機(jī)構(gòu)的姿態(tài)運(yùn)動描述方法

圖1 俯仰/滾轉(zhuǎn)耦合運(yùn)動Fig.1 Pitch/roll coupled motion

圖2 偏航/滾轉(zhuǎn)耦合運(yùn)動Fig.2 Yaw/roll coupled motion

圖3 俯仰/偏航耦合運(yùn)動Fig.3 Pitch/yaw coupled motion

Aα=I3×3cosα+(1-cosα)e·eT-e×sinα

(1)

式中：I3×3為單位矩陣；e×定義為

將式(1)展開，可得

(2)

繞體軸旋轉(zhuǎn)的方向余弦矩陣為

(3)

式中：Cα代表cosα,Sα代表sinα,其余類似，兩自由度動態(tài)試驗包括了俯仰/滾轉(zhuǎn)、俯仰/偏航、偏航/滾轉(zhuǎn)3種類型，為了統(tǒng)一輸入輸出接口，將兩自由度運(yùn)動定義為“繞體軸旋轉(zhuǎn)+繞地面系下任意軸旋轉(zhuǎn)”兩次連續(xù)旋轉(zhuǎn)過程。用Xd和Xt表示地軸系和體軸系下的坐標(biāo)，Aφ,θ,ψ為繞體軸旋轉(zhuǎn)的方向余弦矩陣，得到坐標(biāo)變換關(guān)系為

等價于：

Xt=Aφ,θ,ψAαXd

(4)

式(4)說明了雙自由度運(yùn)動的坐標(biāo)變換包含了兩次矩陣左乘的運(yùn)算，注意此時運(yùn)算結(jié)果與乘法順序相關(guān)，矩陣的順序不能互換。最后得到的坐標(biāo)變換矩陣為

A=Aφ,θ,ψAα

(5)

下面考慮方向余弦矩陣向歐拉角的轉(zhuǎn)換，如果考慮俯仰角在-90°～+90°之間取值，滾轉(zhuǎn)角和偏航角在-180°～+180°之間取值，則方向余弦矩陣到歐拉角的轉(zhuǎn)換公式為

(6)

θ=arcsin(-A13)

(7)

(8)

(9)

ifu≠0 then

{α=arctan(w/u)

β=arctan(-(vcosα)/u)

elseifw≠0 then

{α=(π/2)sign(w)

β=arctan(-v/w)

elseifw=0 then

{α=π/2

β=(π/2)sign(-v)

2 數(shù)值計算方法

流動控制方程采用任意拉格朗日-歐拉(Arbitrary Lagrangian Eulerian，ALE)形式下的可壓縮非定常雷諾平均Navier-Stokes(URANS)方程，允許網(wǎng)格的任意運(yùn)動和變形，在絕對坐標(biāo)系下的積分形式為

(10)

式中：Q為守恒變量;Ω表示控制體的體積；S為控制體邊界面的面積;n為單位法向矢量；VΩ為邊界面上的網(wǎng)格運(yùn)動速度；HI和HV分別為對流通量和黏性通量。非定常計算采用雙時間步方法，其中內(nèi)層時間離散采用了點(diǎn)隱式LU-SGS方法，真實時間采用了三層二階精度格式離散，在每個真實時間步內(nèi)通過內(nèi)層時間迭代使氣動力積分值收斂。黏性離散采用全Navier-Stokes方程方法，湍流采用的是兩方程k-ω剪切應(yīng)力輸運(yùn)(SST)模型進(jìn)行模擬。由于是計算低速大迎角問題，為降低低速流動的剛性，流動控制方程采用了低速預(yù)處理方法，同時為加速流場收斂還采用了3重V循環(huán)的多重網(wǎng)格方法，每一重網(wǎng)格的迭代步數(shù)為5步，每個真實時間步內(nèi)多重網(wǎng)格循環(huán)的最大次數(shù)為200步。在當(dāng)前參數(shù)設(shè)置下，基本上都能保證子迭代過程的氣動力收斂，如果在迭代過程中密度殘差下降到兩個量級以上，則被認(rèn)為是子迭代收斂，提前退出循環(huán)。

計算模型如圖4所示，平均氣動弦長(參考長度)為0.426 9 m，來流風(fēng)速為25 m/s，雷諾數(shù)Re=1.71×106，第一層網(wǎng)格距離物面5×10-6m，滿足湍流附面層模擬y+=Ο(1)的要求，網(wǎng)格采用多塊對接網(wǎng)格(包括232塊)，單元總數(shù)約為1 228萬。

圖4 計算模型和滾轉(zhuǎn)振蕩瞬時流場Fig.4 Calculation model and instantaneous flow field of roll oscillation

在開展動態(tài)計算之前，首先通過靜態(tài)流場對計算方法進(jìn)行校核，模擬迎角α=0°～80°。圖5給出的是計算與試驗的升力系數(shù)CL比較，可以看到，在所有迎角范圍內(nèi)計算與試驗值吻合良好，最大升力系數(shù)約為1.48，對應(yīng)的迎角約為35°。圖6給出的是阻力系數(shù)CD曲線比較，可以看到，在迎角α=40°以前，計算與試驗值吻合較好，在α>40°以后，計算值較試驗值偏高，最大偏差量約為10%。圖7給出了俯仰力矩系數(shù)Cm的比較，也可以看到兩者之間吻合較好，并且可以看到該飛行器模型在中小迎角范圍內(nèi)俯仰方向是靜不穩(wěn)定的，俯仰力矩隨迎角增大而增大，而在大迎角下(α>30°)是俯仰靜穩(wěn)定的。總的來說，當(dāng)前計算在規(guī)律和趨勢上和試驗值都比較一致，反映了低速大迎角流場的主要特征。由于目前公認(rèn)的URANS方法模擬大分離流動存在不足，所以認(rèn)為在失速迎角以后計算與試驗值存在一定偏差是合理的，在當(dāng)前計算中主要表現(xiàn)在阻力值的預(yù)測上。

圖5 升力系數(shù)比較Fig.5 Comparison of lift coefficients

圖6 阻力系數(shù)比較Fig.6 Comparison of drag coefficients

圖7 俯仰力矩系數(shù)比較Fig.7 Comparison of pitching moment coefficients

3 計算結(jié)果與分析

3.1 偏航、滾轉(zhuǎn)單自由度大幅振蕩

在計算偏航/滾轉(zhuǎn)雙自由度振蕩之前，首先對偏航振蕩和滾轉(zhuǎn)振蕩分別進(jìn)行計算，單自由度滾轉(zhuǎn)振蕩的運(yùn)動規(guī)律用正弦函數(shù)表示，φ=40°·sin(2πft),振蕩頻率f=0.8 Hz。圖4給出的是計算模型和滾轉(zhuǎn)振蕩瞬時流場，對應(yīng)瞬時滾轉(zhuǎn)角為34.6°，分別給出了橫向截面的馬赫數(shù)Ma云圖、空間流線與渦量等值面，可以看到在當(dāng)前位置下，飛行器左右側(cè)旋渦的位置和強(qiáng)度都發(fā)生改變，模擬得到的物面邊界層與空間旋渦的形態(tài)分布合理。

圖8給出了滾轉(zhuǎn)和偏航力矩系數(shù)(Cl,Cn)的計算值與試驗值比較，可以看到計算與試驗值吻合很好，進(jìn)一步驗證了當(dāng)前計算方法的合理性。盡管計算流體力學(xué)存在計算坐標(biāo)系、飛機(jī)體軸系之分，這里為敘述方便將角度定義的方向和氣動力矩的方向都統(tǒng)一到機(jī)體坐標(biāo)系下。在轉(zhuǎn)動方向與力矩方向定義一致的前提下，可以根據(jù)力矩曲線的斜率或方向來判斷飛機(jī)的力矩特性[10]。具體來說，靜態(tài)情況下如果力矩系數(shù)的斜率為負(fù)，說明在該方向上的力矩是靜穩(wěn)定的，否則為靜不穩(wěn)定；在動態(tài)情況下，如果遲滯回線沿逆時針方向旋轉(zhuǎn)，則說明在該方向上的運(yùn)動對力矩起阻尼作用，如果遲滯回線沿順時針方向旋轉(zhuǎn)，則說明在該方向上的運(yùn)動對力矩起發(fā)散作用。

圖8 滾轉(zhuǎn)振蕩的滾轉(zhuǎn)與偏航力矩系數(shù)(α=10°)Fig.8 Roll and yaw moment coefficients in roll oscillation (α=10°)

如圖8所示，對于滾轉(zhuǎn)力矩來說，從最小滾轉(zhuǎn)角對應(yīng)的力矩點(diǎn)和最大滾轉(zhuǎn)角對應(yīng)的力矩點(diǎn)之間連線，可以得出滾轉(zhuǎn)力矩呈左上右下的走勢，即正的滾轉(zhuǎn)角對應(yīng)的滾轉(zhuǎn)力矩為負(fù)，說明飛機(jī)在滾轉(zhuǎn)方向上是靜穩(wěn)定的。但是對當(dāng)前滾轉(zhuǎn)大幅振蕩來說，滾轉(zhuǎn)力矩形成了很大的遲滯圈，即滾轉(zhuǎn)力矩產(chǎn)生了很強(qiáng)的動態(tài)非定常效應(yīng)。圖中標(biāo)注了遲滯回線的旋轉(zhuǎn)方向，可以看到滾轉(zhuǎn)力矩曲線呈逆時針走向，并且滾轉(zhuǎn)力矩環(huán)呈O型，這說明在整個滾轉(zhuǎn)角歷程中滾轉(zhuǎn)運(yùn)動對滾轉(zhuǎn)力矩的變化起阻尼作用。圖中同時還給出了偏航力矩的遲滯回線，可以看到當(dāng)前偏航力矩并非小量，說明飛行器滾轉(zhuǎn)和偏航之間存在較強(qiáng)的耦合作用，同時偏航力矩遲滯回線的旋轉(zhuǎn)為逆時針方向，意味著滾轉(zhuǎn)引起的偏航交叉動導(dǎo)數(shù)起阻尼作用。在振蕩過程中，當(dāng)飛行器處于平衡位置時，圖8中滾轉(zhuǎn)力矩遲滯圈包圍面積很大，對比上行和下行狀態(tài)，當(dāng)飛行器處于平衡位置時，兩者的姿態(tài)角和氣流迎角、側(cè)滑角都一致，角速度不同、時間歷程不同，造成了滾轉(zhuǎn)力矩系數(shù)的較大差異。

圖9 偏航振蕩的滾轉(zhuǎn)與偏航力矩系數(shù) (α=10°)Fig.9 Roll and yaw moment coefficients in yaw oscillation (α=10°)

單自由度偏航振蕩的運(yùn)動規(guī)律用正弦函數(shù)表示，ψ=40°sin(2πft),振蕩頻率f=0.8 Hz。圖9給出的是偏航振蕩的計算與試驗值比較，包含滾轉(zhuǎn)和偏航力矩系數(shù)，可以看到計算與試驗在趨勢與量值上都比較一致，可以說計算正確反映了偏航運(yùn)動下的橫航向力矩特性變化，不足之處是在較大偏航角下，計算的偏航力矩與滾轉(zhuǎn)力矩峰值較試驗值偏小，其余地方兩者都吻合較好。對于偏航力矩來說，在偏航角ψ=-20°～20°區(qū)間，偏航力矩呈左高右低的走向，說明此時偏航力矩與運(yùn)動方向相反，有利于飛行器的穩(wěn)定；但是在ψ<-20°和ψ>20°的區(qū)間，偏航力矩的走向與運(yùn)動方向相同，不利于偏航方向的穩(wěn)定。由于偏航力矩主要是飛行器垂尾產(chǎn)生，在較大的偏航角下，垂尾效率下降甚至失速，會導(dǎo)致偏航穩(wěn)定的效率下降。從總的曲線變化趨勢上看，偏航力矩遲滯環(huán)呈逆時針變化，動導(dǎo)數(shù)起阻尼作用。另外，由偏航引起的滾轉(zhuǎn)力矩遲滯環(huán)呈順時針變化，偏航引起的滾轉(zhuǎn)交叉動導(dǎo)數(shù)起發(fā)散作用。

3.2 偏航/滾轉(zhuǎn)耦合大幅振蕩

偏航/滾轉(zhuǎn)耦合的運(yùn)動機(jī)構(gòu)見圖2，機(jī)構(gòu)偏角定義為

(11)

圖10 偏航/滾轉(zhuǎn)運(yùn)動中的歐拉角、迎角和側(cè)滑角變化 Fig.10 Euler angle，angle of attack and sideslip angle variations in yaw/roll coupled motion

設(shè)定滾轉(zhuǎn)和偏航的振蕩幅值φA=ψA=40°,振蕩頻率均為f=0.8 Hz,初始相位相同，即飛行器從平衡位置開始向右偏航、向右滾轉(zhuǎn)。圖10給出的是偏航/滾轉(zhuǎn)耦合運(yùn)動下的歐拉角、迎角和側(cè)滑角變化曲線，可以看到，在運(yùn)動過程中俯仰角θ始終為0°，偏航角ψ和滾轉(zhuǎn)角φ在±40°之間變化，偏航與滾轉(zhuǎn)方向的歐拉角與機(jī)構(gòu)運(yùn)動偏角相等。圖中同時給出了迎角、側(cè)滑角的變化，其中平衡位置的來流迎角α=10°，可以看到飛行器在姿態(tài)角發(fā)生變化的同時，氣流角也在發(fā)生變化，其中迎角在10°～36°之間變化，側(cè)滑角則在±22°之間變化。

圖11給出的是體軸系的3個角速度變化曲線，可以看到，除了偏航和滾轉(zhuǎn)方向的角速度外，俯仰方向上也產(chǎn)生了角速度，并且俯仰角速度的周期是其他兩個方向角速度周期的一半。由此可以看到，偏航/滾轉(zhuǎn)耦合運(yùn)動并不是偏航運(yùn)動和滾轉(zhuǎn)運(yùn)動的簡單疊加，其中包括了多個動態(tài)參數(shù)的變化，根據(jù)線性理論假設(shè)，將氣動力/力矩系數(shù)通常表示為飛行狀態(tài)參數(shù)的函數(shù)[21]，即

Cj0(α,β,V∞,h,δi)+Cjd(α,β,V∞,h,δi,p)+

Cjd(α,β,V∞,h,δi,q)+…

(12)

圖12給出的是偏航/滾轉(zhuǎn)振蕩的升力系數(shù)CL、阻力系數(shù)CD和側(cè)向力系數(shù)Cy比較，從圖中可以看到，升力系數(shù)和阻力系數(shù)呈“∞”形變化，側(cè)向力系數(shù)呈雙“∞”形變化。計算值與試驗值總體上吻合較好，一定程度說明本文算法的正確性，不足之處是最大振蕩角度下計算的升力系數(shù)值略小于試驗值，同時也導(dǎo)致了該角度下的阻力系數(shù)較試驗值小。

圖11 偏航/滾轉(zhuǎn)運(yùn)動中的體軸系角速度變化Fig.11 Angular velocity variations of body axial system in yaw/roll coupled motion

滾轉(zhuǎn)力矩和偏航力矩變化是偏航/滾轉(zhuǎn)耦合振蕩重點(diǎn)關(guān)注的對象。圖13給出了滾轉(zhuǎn)力矩系數(shù)隨滾轉(zhuǎn)角變化曲線。可以看到當(dāng)前力矩曲線呈兩個“∞”形分布，在滾轉(zhuǎn)角-18°<φ<18° 區(qū)域，曲線為逆時針旋轉(zhuǎn)，說明滾轉(zhuǎn)動導(dǎo)數(shù)起阻尼作用，而在φ<-18°和φ>18°的兩個區(qū)域，曲線為順時針旋轉(zhuǎn)，說明此時滾轉(zhuǎn)動導(dǎo)數(shù)起發(fā)散作用。

圖14給出了耦合振蕩下的偏航力矩系數(shù)變化曲線，可以看到當(dāng)前曲線的旋轉(zhuǎn)方向為逆時針，說明耦合振蕩引起的偏航動導(dǎo)數(shù)起阻尼作用。圖13和圖14分別給出了滾轉(zhuǎn)和偏航力矩系數(shù)的風(fēng)洞試驗值，從計算值與試驗值的比較來看，兩者在曲線形態(tài)和量值上都吻合較好，說明當(dāng)前采用的計算方法是合理的。

圖12 偏航/滾轉(zhuǎn)振蕩的升力、阻力和側(cè)向力系數(shù)比較Fig.12 Comparison of lift, drag and side force coefficients in yaw/roll oscillation

圖13 偏航/滾轉(zhuǎn)振蕩的滾轉(zhuǎn)力矩系數(shù)Fig.13 Roll moment coefficients in yaw/roll oscillation

圖14 偏航/滾轉(zhuǎn)振蕩的偏航力矩系數(shù)Fig.14 Yaw moment coefficients in yaw/roll oscillation

圖15 耦合運(yùn)動與單獨(dú)運(yùn)動疊加的偏航力矩系數(shù)比較Fig.15 Comparison of yaw moment coefficients between couple motion and single freedom motion addition

下面考慮耦合運(yùn)動與單獨(dú)運(yùn)動疊加結(jié)果的比較，圖15給出的是偏航力矩系數(shù)曲線，可以看到耦合結(jié)果與疊加結(jié)果具有較高的相似性，并且與單獨(dú)偏航運(yùn)動的結(jié)果相似，可見該飛行器的偏航力矩主要由偏航運(yùn)動產(chǎn)生，相比之下滾轉(zhuǎn)振蕩產(chǎn)生的偏航力矩較小。圖16給出的是滾轉(zhuǎn)力矩系數(shù)曲線，可以看到，當(dāng)前耦合運(yùn)動的力矩系數(shù)在大的趨勢上與圖9偏航振蕩比較相近，呈左下右上的分布形態(tài)。經(jīng)過簡單疊加得到的力矩曲線在較大滾轉(zhuǎn)角范圍內(nèi)(-32°<φ<32°)都是沿逆時針方向旋轉(zhuǎn)，在φ>32°和φ<-32°之后曲線出現(xiàn)了交叉，但包圍的面積很小。出現(xiàn)交叉的原因是滾轉(zhuǎn)運(yùn)動和偏航運(yùn)動的力矩曲線旋轉(zhuǎn)方向相反，具體來說，滾轉(zhuǎn)運(yùn)動產(chǎn)生的滾轉(zhuǎn)力矩系數(shù)曲線沿逆時針方向旋轉(zhuǎn)，偏航運(yùn)動產(chǎn)生的滾轉(zhuǎn)力矩系數(shù)沿順時針方向旋轉(zhuǎn)。與單獨(dú)運(yùn)動疊加對比，耦合運(yùn)動在小滾轉(zhuǎn)角下表現(xiàn)出了類似的形態(tài)，力矩曲線沿逆時針方向旋轉(zhuǎn)，不同的是當(dāng)滾轉(zhuǎn)角較大時(φ>18°和φ<-18°)，力矩曲線出現(xiàn)了交叉，旋轉(zhuǎn)方向變?yōu)轫槙r針，并且每個順時針旋轉(zhuǎn)曲線包圍的面積較大，與逆時針旋轉(zhuǎn)所包圍的面積相當(dāng)。

另外，由于耦合運(yùn)動的滾轉(zhuǎn)力矩曲線形態(tài)發(fā)生改變，當(dāng)飛行器從平衡位置向右側(cè)滾轉(zhuǎn)時，隨滾轉(zhuǎn)角增大時，滾轉(zhuǎn)力矩加速上揚(yáng)，在φ=34°附近達(dá)到最大值，該值較單獨(dú)運(yùn)動和單獨(dú)運(yùn)動疊加對應(yīng)的滾轉(zhuǎn)力矩要高出很多，力矩的作用方向是進(jìn)一步加速飛行器的滾轉(zhuǎn)運(yùn)動。

圖16 耦合運(yùn)動與單獨(dú)運(yùn)動疊加的滾轉(zhuǎn)力矩系數(shù)比較Fig.16 Comparison of roll moment coefficients between couple motion and single freedom motion addition

圖17 耦合運(yùn)動與單獨(dú)運(yùn)動疊加的側(cè)滑角比較Fig.17 Comparison of side slip angles between couple motion and single freedom motion addition

為了弄清力矩曲線變化的原因，首先對耦合運(yùn)動的姿態(tài)角與氣流角變化進(jìn)行了分析，從圖10可以看到，耦合運(yùn)動除了偏航和滾轉(zhuǎn)方向的姿態(tài)角變化以外，可能對滾轉(zhuǎn)力矩產(chǎn)生較大影響的是側(cè)滑角的變化。為進(jìn)一步分析，圖17對單獨(dú)運(yùn)動和耦合運(yùn)動的側(cè)滑角進(jìn)行了比較，可以看見，單獨(dú)偏航運(yùn)動的側(cè)滑角在[-40°,40°] 區(qū)間變化，單獨(dú)滾轉(zhuǎn)運(yùn)動的側(cè)滑角在[-6.4°,6.4°] 區(qū)間變化，但方向相反，耦合運(yùn)動的側(cè)滑角約在[-22°,22°] 區(qū)間變化。從圖9知道，偏航運(yùn)動當(dāng)偏航角在25°之前產(chǎn)生的滾轉(zhuǎn)力矩單調(diào)上升，25°以后開始下降，這其中起主導(dǎo)作用的就是側(cè)滑角的變化。而對于耦合運(yùn)動來說，側(cè)滑角始終在±25°內(nèi)變化，該角度范圍不足以使垂尾發(fā)生失速，使?jié)L轉(zhuǎn)力矩得以持續(xù)上升，這也就解釋了圖16耦合振蕩的滾轉(zhuǎn)力矩系數(shù)峰值增大的原因。從本文分析可以看到，非定常氣動系數(shù)的影響因素很多，包括模型姿態(tài)角、氣流角、對應(yīng)角速度等，非定常運(yùn)動的時間歷程效應(yīng)也是影響因素之一，這些因素的疊加導(dǎo)致了耦合運(yùn)動氣動特性的復(fù)雜性，同時給氣動系數(shù)建模帶來了難度。

4 結(jié) 論

采用URANS方法模擬預(yù)定運(yùn)動軌跡,建立了專門針對飛行器雙自由度大振幅振蕩的數(shù)值模擬方法，完成了某復(fù)雜外形飛行器靜態(tài)、偏航振蕩、滾轉(zhuǎn)振蕩和偏航/滾轉(zhuǎn)耦合振蕩的模擬，結(jié)合試驗數(shù)據(jù)對力矩特性的分布規(guī)律進(jìn)行了分析。

1) “PQR”試驗裝置在進(jìn)行雙自由度試驗時，模型方位與旋轉(zhuǎn)順序無關(guān)，在用CFD進(jìn)行預(yù)定軌跡計算時，將用戶輸入用一次繞任意軸旋轉(zhuǎn)疊加一次繞體軸旋轉(zhuǎn)表示，輸出用坐標(biāo)變換的歐拉角及繞體軸的角速度表示，可以實現(xiàn)對任意俯仰/偏航、俯仰/滾轉(zhuǎn)、偏航/滾轉(zhuǎn)等耦合運(yùn)動的描述。

2) 當(dāng)前計算采用了低速預(yù)處理、多重網(wǎng)格等加速收斂技術(shù)，保證非定常計算每個子迭代過程氣動力的充分收斂，得到的計算結(jié)果合理，在模型平衡位置迎角10°，滾轉(zhuǎn)±40°偏角、偏航±40°偏角下的結(jié)果與風(fēng)洞試驗值吻合一致。

3) 飛行器雙自由度運(yùn)動不是兩個單自由度運(yùn)動的簡單疊加，氣動特性具有單自由度運(yùn)動的部分特點(diǎn)，但遲滯回線的形態(tài)更加復(fù)雜，甚至表現(xiàn)出不同的運(yùn)動阻尼特性。

4) 本文偏航/滾轉(zhuǎn)耦合振蕩的平衡位置迎角為10°，流場以橫向分離渦運(yùn)動為主要特征，可以推測大迎角情況下縱橫向分離渦流動結(jié)構(gòu)更加復(fù)雜，耦合作用更強(qiáng)，需要進(jìn)一步在大分離流模擬方法上開展研究。

本文模擬外形為1∶16風(fēng)洞縮尺模型，來流風(fēng)速為25 m/s，基于平均氣動弦長的雷諾數(shù)為百萬量級，當(dāng)前條件下本文URANS方法取得了與試驗值比較吻合的結(jié)果，真實飛行條件下飛行器速度更快、雷諾數(shù)更高，對計算方法構(gòu)成挑戰(zhàn)。下一步有必要開展真實飛行條件下的計算研究，依托高雷諾數(shù)風(fēng)洞試驗來驗證CFD方法；同時分離渦流動模擬是當(dāng)前的熱門研究方向，開展這類方法研究將有助于解決大迎角大分離流動的動態(tài)模擬問題。

致 謝

本研究得到了汪清研究員、劉志濤副研究員的幫助和指導(dǎo)，在此表示感謝。

[1] NELSON R C, PELLETIER A. The unsteady aerodynamics of slender wings and aircraft undergoing large amplitude maneuvers[J]. Progress in Aerospace Sciences, 2003, 39(2): 185-248.

[2] KRAMER B R. Experimental evaluation of superposition techniques applied to dynamic aerodynamics: AIAA-2002-0700[R]. Reston, VA: AIAA, 2002.

[3] LAN C E, LI J, YAU W, et al. Longitudinal and lateral-directional coupling effects on nonlinear unsteady aeodynamic modeling from flight data: AIAA-2002-4804[R]. Reston, VA： AIAA, 2002.

[4] WANG Z, LAN C E, BRANDON J M. Fuzzy logic modeling of lateral-directional unsteady aerodynamics：AIAA-1999-4012[R]. Reston, VA： AIAA, 1999.

[5] JENKINS J E. Nonlinear aerodynamic characteristics of a 65 degree delta wing in rolling motion: Implications for testing and flight mechanics analysis：AIAA-1997-0742[R].Reston, VA： AIAA, 1997.

[6] OWENS D B, BRANDON J M, CROOM M A, et al. Overview of dynamic test techniques for flight dynamics research at NASA LaRC: AIAA-2006-3146[R]. Reston, VA： AIAA, 2006.

[7] GILLARD W J. AFRL F-22 dynamic wind tunnel test results: AIAA-1999-04015[R]. Reston, VA： AIAA, 1999.

[8] 李偉. 飛機(jī)大振幅動態(tài)實驗及其非定常氣動特性研究[D]. 南京: 南京航空航天大學(xué), 2008.

LI W. Unsteady aerodynamic experiments and characteristics for the aircraft oscillation in large amplitude[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008 (in Chinese).

[9] 李周復(fù). 風(fēng)洞特種試驗技術(shù)[M]. 北京: 航空工業(yè)出版社, 2010.

LI Z F. Special test technology in wind tunnel[M]. Beijing: Aviation Industry Press, 2010 (in Chinese).

[10] 黃達(dá), 吳根興. 三角翼俯仰滾轉(zhuǎn)耦合運(yùn)動氣動特性研究[J]. 航空學(xué)報, 1999, 20(6): 485-488.

HUANG D, WU G X. Investigation of unsteady aerodynamic characteristics for a delta wing oscillating in large amplitude pitching roll motion[J]. Acta Aeronautica et Astronautica Sinica, 1999, 20(6): 485-488 (in Chinese).

[11] 黃達(dá), 吳根興. 飛機(jī)偏航-滾轉(zhuǎn)耦合運(yùn)動非定?？諝鈩恿嶒瀃J]. 南京航空航天大學(xué)學(xué)報, 2005, 37(4): 408-411.

HUANG D, WU G X. Experiment on fighter oscillating in large amplitude yaw-roll motion[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2005, 37(4): 408-411 (in chinese) .

[12] HUANG D, WU G X. Unsteady rolling moment characteristics for fighter oscillating with yawing-rolling coupled motion[J]. Journal of Aircraft, 2006, 43(5): 1570-1573.

[13] 劉春明, 趙志軍, 卜忱. 低速風(fēng)洞雙自由度大幅振蕩試驗技術(shù)[J]. 航空學(xué)報, 2016, 37(8): 2417-2425.

LIU C M, ZHAO Z J, BU C， et al. Double degree-of-freedom large amplitude oscillation test technology in low speed wind tunnel[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2417-2425 (in Chinese).

[14] 楊文, 卜忱, 眭建軍. 某復(fù)雜構(gòu)型飛機(jī)偏航-滾轉(zhuǎn)耦合運(yùn)動非定常氣動力特性實驗研究[J]. 實驗流體力學(xué), 2016, 30(3): 61-65.

YANG W, BU C, GUI J J. Investigation of the unsteady aerodynamic characteristics of a fighter with complex configuration undergoing yaw-roll coupling oscillation motion[J]. Journal of Experiments in Fluid Mechanics, 2016, 30(3): 61-65 (in Chinese).

[15] 孫海生, 張海酉, 劉志濤. 大迎角非定常氣動力建模方法研究[J]. 空氣動力學(xué)學(xué)報, 2011, 29(6): 733-737.

SUN H S, ZHANG H Y, LIU Z T. Comparative evaluation of unsteady aerodynamics modeling approaches at high angle of attack[J]. Acta Aerodynamica Sinica, 2011, 29(6): 733-737 (in Chinese).

[16] 汪清, 錢煒祺, 丁娣. 飛機(jī)大迎角非定常氣動力建模研究進(jìn)展[J]. 航空學(xué)報, 2016, 37(8): 2332-2347.

WANG Q, QIAN W Q, DING D. A review of unsteady aeeodynamic modeling of aircrafts at high angles of attack[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2332-2347 (in Chinese).

[17] 史志偉, 黃達(dá), 吳根興, 等. 耦合運(yùn)動非定常氣動模型對飛機(jī)飛行特性仿真的影響[J]. 航空學(xué)報, 2008, 29(6): 1424-1428.

SHI Z W, HUANG D, WU G X, et al. Effects of coupled motion unsteady aerodynamic model on flight characteristics simulation of aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(6): 1424-1428 (in Chinese).

[18] JANTZEN R, TAIRA K. Aerodynamic force modeling for unsteady wing maneuvers: AIAA-2014-0223[R]. Reston, VA: AIAA, 2014.

[19] MICHAEL V L. The high frequency, high-amplitude pitch problem: Airfoils, plates and wings: AIAA-2009-3686[R]. Reston, VA: AIAA, 2009.

[20] 張帆, 曹喜濱, 鄒經(jīng)湘. 一種新的全角度四元數(shù)與歐拉角的轉(zhuǎn)換算法[J]. 南京理工大學(xué)學(xué)報, 2002, 26(4): 376-380.

ZHANG F, CAO X B, ZOU J X. A new large-scale transformation algorithm of quaternion to Euler angle[J]. Journal of Nanjing University of Science and Technology, 2002, 26(4): 376-380 (in Chinese).

[21] WANG F J, CHEN L. Numerical prediction of stability derivatives for complex configurations[J]. Procedia Engineering, 2015, 99(1): 1561-1575.

Numericalsimulationmethodsfor2DOFcoupledoscillationsofaircraft

XIAOZhongyun*,LIUGang,JIANGXiong,WANGJiantao

ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Focusingonthe“PQR”forcedoscillationtestapparatusformodelingTwo-Degree-of-Freedom(2DOF)motions,coordinatetransformationisperformedtoobtainunifiedEulerangledescriptionsofthepitch/yaw,pitch/rollandyaw/rollcoupledmotions.BasedonmovinggridsforprescribedmotionandUnsteadyReynoldsAveragedNavier-Stokes(URANS)equations,anumericalalgorithmisdevelopedtosimulatecoupledoscillationofaircrafts.Surroundinglowspeedflowsofacomplexaircraft,single-degreeoscillationsofyawandrollandyaw/rollcoupledoscillationarerespectivelycomputedatangleofattackof10°andamplitudeofoscillationof40°.Thenumericalresultsagreewellwiththeexperimentaldata,showingthatcurrentalgorithmhaspromisingpotentialsinsimulationofcomplicatedcoupledmotionsofaircrafts.Asformovingprocessandaerodynamics,theyaw/rollcoupledmotioniscomparedwiththesuperpositionofsingleyawandrolloscillations,whichshowsthatdifferencesexistinbodyangularvelocitiesandflowangles,andthemomenthysteresiscurveoftheyaw/rollcoupledmotionshowscrossings,indicatingtotallydifferentdampingcharacteristicsfromsingleDOFoscillation.Incurrentyawandrolloscillations,flowsarecharacterizedbylaterallyseparatedvortexes.Itisassumedthatathighangleofattack,whenflowsaremorecomplicatedwithlateralandlongitudinalvortexesandstrongcoupledeffects,morestudiesshouldbedoneonnumericalmethodsofhighlyseparatedflows.

2DOF;largeamplitudeoscillation;coupledmotion;unsteady;numericalsimulation

2017-04-21；Revised2017-05-16；Accepted2017-06-07；Publishedonline2017-06-161118

URL：http://hkxb.buaa.edu.cn/CN/html/20171111.html

s：NationalBasicResearchProgramofChina(613221)；NationalNaturalScienceFoundationofChina(11572341);NationalKeyResearchandDevelopmentProgram(2016YFB0200704)

.E-mailscxiaozy@sina.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121345

V211.3

A

1000-6893(2017)11-121345-10

2017-04-21；退修日期2017-05-16；錄用日期2017-06-07；< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2017-06-161118

http://hkxb.buaa.edu.cn/CN/html/20171111.html

國家“973”計劃(613221)； 國家自然科學(xué)基金(11572341)； 國家重點(diǎn)研發(fā)計劃(2016YFB0200704)

.E-mailscxiaozy@sina.cn

肖中云，劉剛，江雄，等．飛行器雙自由度耦合振蕩的數(shù)值模擬方法J． 航空學(xué)報,2017,38(11)：121345.XIAOZY,LIUG,JIANGX,etal.Numericalsimulationmethodsfor2DOFcoupledoscillationsofaircraftJ.ActaAeronauticaetAstronauticaSinica,2017,38(11):121345.

(責(zé)任編輯：李明敏)