改進(jìn)的準(zhǔn)定常旋轉(zhuǎn)坐標(biāo)系阻尼動導(dǎo)數(shù)計(jì)算方法

吳耕宇，鮑君波，宋萬強(qiáng)，相倩

中國航空研究院，北京 100012

飛行器動導(dǎo)數(shù)是飛行器姿態(tài)控制系統(tǒng)設(shè)計(jì)、彈道設(shè)計(jì)及動態(tài)穩(wěn)定性分析中的關(guān)鍵參數(shù)[1-2]，在高超聲速飛行器[3]、可重復(fù)使用的跨大氣層軌道飛行器[4]等各類先進(jìn)飛行器的動態(tài)穩(wěn)定性分析中具有重要作用。計(jì)算流體力學(xué)（CFD）技術(shù)目前已廣泛應(yīng)用于動導(dǎo)數(shù)的計(jì)算。使用CFD 技術(shù)計(jì)算動導(dǎo)數(shù)主要有準(zhǔn)定常和非定常兩種方法，準(zhǔn)定常方法一般通過模擬定常拉升或勻速滾轉(zhuǎn)運(yùn)動計(jì)算俯仰、滾轉(zhuǎn)、偏航阻尼動導(dǎo)數(shù)，非定常方法一般通過模擬強(qiáng)迫振蕩等方法計(jì)算氣動力系數(shù)相對于俯仰角速度和迎角變化率的組合、氣動力系數(shù)相對于偏航角速度和側(cè)滑角變化率的組合等組合動導(dǎo)數(shù)。國內(nèi)外已有多位學(xué)者使用CFD 方法進(jìn)行動導(dǎo)數(shù)計(jì)算的研究。S.M.urman等[5]采用減縮頻率方法對基本尾翼導(dǎo)彈（BFM）、改進(jìn)的基本尾翼導(dǎo)彈（MBFM）和標(biāo)準(zhǔn)動力模型（SDM）三個標(biāo)模進(jìn)行動導(dǎo)數(shù)計(jì)算；Le Roy 等[6]對穩(wěn)定性和控制構(gòu)型（SACCON）翼身融合體布局進(jìn)行了靜導(dǎo)數(shù)和動導(dǎo)數(shù)的CFD計(jì)算和分析；Da Ronch等[7]使用Euler方程和雷諾時均納維-斯托克斯（RANS）方程分別對SDM 標(biāo)模和跨聲速巡航標(biāo)模（TCR）進(jìn)行了非定常動導(dǎo)數(shù)的計(jì)算，并進(jìn)行了較為詳細(xì)的分析；袁先旭等[8]使用時空二階精度的隱式迭代無波動、無自由參數(shù)的耗散差分（NND）算法對尖錐、鈍錐、彈道外形和飛船返回艙等典型再入飛行器進(jìn)行了俯仰靜、動導(dǎo)數(shù)的計(jì)算；范晶晶等[9]采用逐次超松弛時間推進(jìn)方法耦合求解非定常納維-斯托克斯（N-S）方程和強(qiáng)迫運(yùn)動方程，對美國國家航空咨詢委員會（NACA）開發(fā)的標(biāo)準(zhǔn)翼型（NACA0012）、彈道外形和有翼導(dǎo)彈等飛行器進(jìn)行了強(qiáng)迫俯仰振蕩的黏性動態(tài)流場求解；葉川等[10]分別采用基于強(qiáng)迫振蕩的非定常方法和基于定常拉升的準(zhǔn)定常旋轉(zhuǎn)坐標(biāo)系方法計(jì)算BFM 模型和水上飛機(jī)模型的組合動導(dǎo)數(shù)和阻尼動導(dǎo)數(shù)，該方法計(jì)算機(jī)體坐標(biāo)系下非零迎角滾轉(zhuǎn)阻尼動導(dǎo)數(shù)需要在計(jì)算穩(wěn)定坐標(biāo)系下滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)的同時計(jì)算穩(wěn)定坐標(biāo)系下偏航阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)；米百剛等[11]通過強(qiáng)迫振蕩方法和基于圓環(huán)網(wǎng)格的定常拉升方法辨識單獨(dú)動導(dǎo)數(shù)，并采用BFM 模型進(jìn)行驗(yàn)證；張一帆等[12]對歐洲動導(dǎo)數(shù)標(biāo)模（DLR-F12）標(biāo)準(zhǔn)模型進(jìn)行非定常強(qiáng)迫振蕩和準(zhǔn)定常旋轉(zhuǎn)坐標(biāo)系方法，獲得良好的組合動導(dǎo)數(shù)和阻尼動導(dǎo)數(shù)結(jié)果；相倩等[13]計(jì)算了BFM模型大迎角狀態(tài)下的非定常動導(dǎo)數(shù)，反映了流場非線性特征；朱海濤等[14]使用非定常強(qiáng)迫振蕩計(jì)算方法對TCR 標(biāo)模進(jìn)行了動導(dǎo)數(shù)計(jì)算的詳細(xì)研究；宋萬強(qiáng)等[15]采用在物面處添加法向擾動速度模擬強(qiáng)迫振蕩的方法計(jì)算DLR-F12 標(biāo)模的動導(dǎo)數(shù)，實(shí)現(xiàn)無網(wǎng)格變形的非定常動導(dǎo)數(shù)計(jì)算。非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)計(jì)算是動導(dǎo)數(shù)計(jì)算的一種特殊情況，國內(nèi)外對于該情況具有一定程度的研究，但相比俯仰阻尼動導(dǎo)數(shù)的計(jì)算研究較少。葉川等[16]利用渦格法和理論推導(dǎo)得到的動導(dǎo)數(shù)計(jì)算公式計(jì)算臨近空間長航時太陽能飛行器的非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)，并研究總體氣動參數(shù)對動導(dǎo)數(shù)的影響機(jī)理；岳杰順等[17]、王紀(jì)林等[18]采用多參考系方法計(jì)算非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)，將計(jì)算網(wǎng)格分為旋轉(zhuǎn)網(wǎng)格和固定網(wǎng)格兩部分，模型附近為旋轉(zhuǎn)網(wǎng)格，模型繞縱軸旋轉(zhuǎn)，這種方法主要應(yīng)用在彈箭等較為接近圓柱體構(gòu)型的滾轉(zhuǎn)動導(dǎo)數(shù)計(jì)算，對于飛機(jī)構(gòu)型來說，飛機(jī)旋轉(zhuǎn)90°時，迎角將變?yōu)閭?cè)滑角，該方法無法正確模擬；V.A.Bhagwandin[19]采用模型整體以恒定角速度旋轉(zhuǎn)的非定常方法計(jì)算了BFM和MBFM標(biāo)模在-5°～90°迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)等多個氣動導(dǎo)數(shù)，并與試驗(yàn)值進(jìn)行對比；劉偉等[20]采用基于簡諧振動的非定常方法計(jì)算非零迎角下的動導(dǎo)數(shù)，辨識得到滾轉(zhuǎn)組合動導(dǎo)數(shù)；B.Etkin[21]提出計(jì)算機(jī)體坐標(biāo)系下非零迎角滾轉(zhuǎn)阻尼動導(dǎo)數(shù)可通過計(jì)算穩(wěn)定坐標(biāo)系下滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)，同時計(jì)算穩(wěn)定坐標(biāo)系下偏航阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)求得，參考文獻(xiàn)[10]提及了該方法，計(jì)算了穩(wěn)定坐標(biāo)系下非零迎角滾轉(zhuǎn)阻尼動導(dǎo)數(shù)，并與機(jī)體坐標(biāo)系下的試驗(yàn)結(jié)果進(jìn)行對比。

本文通過將來流速度分解的處理手段將基于旋轉(zhuǎn)坐標(biāo)系的準(zhǔn)定常CFD 動導(dǎo)數(shù)計(jì)算方法[10]進(jìn)行改進(jìn)。通過將速度矢量分解為垂直于旋轉(zhuǎn)軸和平行于旋轉(zhuǎn)軸兩個分量，以垂直于旋轉(zhuǎn)軸的速度分量計(jì)算旋轉(zhuǎn)半徑并推算參考坐標(biāo)系的旋轉(zhuǎn)中心點(diǎn)，以平行于旋轉(zhuǎn)軸的速度分量作為來流速度矢量計(jì)算氣動力系數(shù)，通過計(jì)算兩個不同角速度的氣動力系數(shù)之差求得動導(dǎo)數(shù)。使用歐洲動導(dǎo)數(shù)標(biāo)準(zhǔn)模型（DLRF12）和跨聲速巡航標(biāo)模（TCR）進(jìn)行滾轉(zhuǎn)、俯仰和偏航動導(dǎo)數(shù)計(jì)算，并與風(fēng)洞試驗(yàn)數(shù)據(jù)、參考文獻(xiàn)[10]的準(zhǔn)定常計(jì)算方法和參考文獻(xiàn)[22]提供的其他CFD 計(jì)算結(jié)果進(jìn)行對比，以對改進(jìn)方法進(jìn)行驗(yàn)證。

1 計(jì)算方法

現(xiàn)有的基于旋轉(zhuǎn)參考坐標(biāo)系的準(zhǔn)定常阻尼動導(dǎo)數(shù)計(jì)算[10,12]一般通過模擬定常拉升/定常繞軸旋轉(zhuǎn)或者模擬旋轉(zhuǎn)軸與來流方向平行的旋轉(zhuǎn)運(yùn)動兩種途徑獲得。定常拉升如圖1所示，飛行器繞飛行器上方的旋轉(zhuǎn)軸以固定角速度q旋轉(zhuǎn)，旋轉(zhuǎn)半徑為R，來流速度為0，此時飛行器的運(yùn)動速度大小為

圖1 飛行器定常拉升Fig.1 Aircraft steady state lifting

在相對飛行器固定的機(jī)體坐標(biāo)系下，定常拉升運(yùn)動時飛行器的運(yùn)動速度大小和方向、旋轉(zhuǎn)角速度大小和方向、迎角和側(cè)滑角均保持不變，因此可用準(zhǔn)定常方法進(jìn)行動導(dǎo)數(shù)計(jì)算。取不同的角速度q1和q2分別計(jì)算相應(yīng)的氣動力或力矩系數(shù)Ci1和Ci2，其中Ci為升力系數(shù)CL、俯仰力矩系數(shù)Cm或其他氣動力或力矩系數(shù)，則相應(yīng)的動導(dǎo)數(shù)Ciq可由式（2）計(jì)算[12]

式中：Cref根據(jù)所計(jì)算動導(dǎo)數(shù)的不同，為縱向參考長度或參考展長，V∞為來流速度，A一般取2，部分計(jì)算模型為與風(fēng)洞試驗(yàn)處理一致取1[12]。定常繞軸旋轉(zhuǎn)選取飛機(jī)左側(cè)或右側(cè)的旋轉(zhuǎn)軸，飛機(jī)繞旋轉(zhuǎn)軸定常運(yùn)動，來流速度為0，動導(dǎo)數(shù)的計(jì)算方法與定常拉升類似，通過計(jì)算兩個不同旋轉(zhuǎn)角速度的氣動力系數(shù)，通過式（2）求得動導(dǎo)數(shù)，通常用于計(jì)算偏航阻尼動導(dǎo)數(shù)。旋轉(zhuǎn)軸與來流方向平行的旋轉(zhuǎn)運(yùn)動，在機(jī)體坐標(biāo)系下，飛行器運(yùn)動速度大小和方向、旋轉(zhuǎn)角速度大小和方向、迎角和側(cè)滑角均保持不變，也可使用基于旋轉(zhuǎn)坐標(biāo)系的準(zhǔn)定常方法計(jì)算，通過計(jì)算兩個不同旋轉(zhuǎn)角速度的氣動力系數(shù)，通過式（2）求得動導(dǎo)數(shù)。

模擬定常拉升/定常繞軸旋轉(zhuǎn)或者模擬旋轉(zhuǎn)軸與來流方向平行的旋轉(zhuǎn)運(yùn)動兩種處理途徑只能計(jì)算部分情況的阻尼動導(dǎo)數(shù)。一種典型的無法使用上述兩種處理途徑直接計(jì)算的情況是計(jì)算非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)。參考文獻(xiàn)[10]中計(jì)算非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)需要通過模擬旋轉(zhuǎn)軸與來流方向平行的旋轉(zhuǎn)運(yùn)動途徑計(jì)算穩(wěn)定坐標(biāo)系下滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)，并通過模擬定常拉升/定常繞軸旋轉(zhuǎn)運(yùn)動途徑計(jì)算穩(wěn)定坐標(biāo)系下偏航阻尼動導(dǎo)數(shù)和交叉動導(dǎo)數(shù)。穩(wěn)定坐標(biāo)系下的滾轉(zhuǎn)運(yùn)動、穩(wěn)定坐標(biāo)系下的偏航運(yùn)動和機(jī)體坐標(biāo)系下的滾轉(zhuǎn)運(yùn)動分別如圖2～圖4 所示。計(jì)算完成后，采用式（3）[10,21]進(jìn)行機(jī)體坐標(biāo)系下動導(dǎo)數(shù)的求解。

圖2 在穩(wěn)定坐標(biāo)系中勻速滾轉(zhuǎn)Fig.2 Rolling in stability coordinate frame

圖3 在穩(wěn)定坐標(biāo)系中繞軸偏航Fig.3 Yawing in stability coordinate frame

圖4 在機(jī)體坐標(biāo)系中勻速滾轉(zhuǎn)Fig.4 Rolling in body coordinate frame

1.1 非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)計(jì)算

假設(shè)迎角不為0，側(cè)滑角為0，將飛機(jī)相對于地面運(yùn)動速度矢量V分解為平行于機(jī)體坐標(biāo)系x軸方向Vx和垂直于機(jī)體坐標(biāo)系x軸方向Vz兩個分量。假定初始狀態(tài)風(fēng)軸系和地面坐標(biāo)系的坐標(biāo)軸平行，則俯仰角和迎角相同，且飛機(jī)沿地面坐標(biāo)系x軸運(yùn)動，由于迎角為飛機(jī)體軸系和風(fēng)軸系的夾角，如圖5所示，有

圖5 飛機(jī)運(yùn)動速度矢量分解示意圖Fig.5 Aircraft velocity vector decomposition map

式中：α為迎角，ib，kb分別代表機(jī)體坐標(biāo)系的單位矢量?？紤]一個較小的單位時間Δt，在Δt時間內(nèi)，飛機(jī)運(yùn)動距離為|VΔt|i，i代表地面坐標(biāo)系x方向的單位矢量。記ω為飛機(jī)旋轉(zhuǎn)角速度，則飛機(jī)旋轉(zhuǎn)角度為ωΔt，下一個Δt時間，飛機(jī)運(yùn)動距離則為

式中：i，j，k分別代表地面坐標(biāo)系的單位矢量。當(dāng)單位時間Δt無限小時，飛機(jī)實(shí)際進(jìn)行一個螺旋運(yùn)動，此時，飛機(jī)的運(yùn)動方向和機(jī)體方向的夾角固定，飛機(jī)運(yùn)動速度大小固定為|V|，且飛機(jī)在單位時間繞機(jī)體坐標(biāo)系x軸的旋轉(zhuǎn)角度固定，即迎角和滾轉(zhuǎn)角速度保持不變，飛機(jī)側(cè)滑角則固定為0。由于飛機(jī)運(yùn)動速度大小保持不變，速度方向在機(jī)體坐標(biāo)系下固定，迎角、側(cè)滑角和滾轉(zhuǎn)角速度均保持不變，因此在機(jī)體坐標(biāo)系下，這個過程的流動是穩(wěn)態(tài)的，流場狀態(tài)不隨時間發(fā)生變化。該螺旋運(yùn)動過程中，飛機(jī)整體的前進(jìn)速度為速度平行于機(jī)體坐標(biāo)系x軸方向的分量，滾轉(zhuǎn)角速度矢量為ω，繞軸線速度為|V|sinα，計(jì)算可得半徑為

因此該過程可用基于旋轉(zhuǎn)坐標(biāo)系的定常方法計(jì)算，來流速度為|V|cosα，方向與旋轉(zhuǎn)軸同向，旋轉(zhuǎn)半徑為，旋轉(zhuǎn)中心點(diǎn)為表示力矩中心點(diǎn)。此時可通過準(zhǔn)定常計(jì)算方法計(jì)算非零迎角下的滾轉(zhuǎn)運(yùn)動氣動力。計(jì)算兩個不同滾轉(zhuǎn)角速度的氣動力，通過式（2）即可得到非零迎角下的滾轉(zhuǎn)動導(dǎo)數(shù)。

1.2 任意方向阻尼動導(dǎo)數(shù)計(jì)算方法

對于任意迎角、側(cè)滑角和任意俯仰角速度、偏航角速度和滾轉(zhuǎn)角速度及其線性組合，也可采取類似方法計(jì)算動導(dǎo)數(shù)。首先建立一個坐標(biāo)系，x軸為機(jī)體旋轉(zhuǎn)軸，按右手螺旋定則與飛機(jī)旋轉(zhuǎn)方向重合，在該坐標(biāo)系下將自由來流方向矢量V分解為平行于該坐標(biāo)系x軸Vx和垂直于該坐標(biāo)系x軸Vn兩個分量，計(jì)算旋轉(zhuǎn)半徑為，其中ω為俯仰角速度、偏航角速度和滾轉(zhuǎn)角速度的線性組合，旋轉(zhuǎn)中心點(diǎn)P=其中PM為力矩中心點(diǎn)，即可計(jì)算任意旋轉(zhuǎn)軸和來流方向下的氣動力。通過計(jì)算兩個不同角速度下的氣動力和式（2）求出動導(dǎo)數(shù)。

2 歐洲動導(dǎo)數(shù)標(biāo)模（DLR-F12）動導(dǎo)數(shù)計(jì)算

2.1 計(jì)算模型

歐洲動導(dǎo)數(shù)標(biāo)模（DLR-F12）是由德國航空航天中心（DLR）研制的標(biāo)準(zhǔn)空氣動力學(xué)模型，機(jī)翼參考面積S=0.4441m2，平均氣動弦長c=0.252625m，展長b=2.03852m，力矩參考點(diǎn)x坐標(biāo)xM=1.04882m，z坐標(biāo)zM=-0.03029m。該模型在德國-荷蘭風(fēng)洞（DNW）進(jìn)行了一系列風(fēng)洞試驗(yàn)[22]，模擬風(fēng)洞試驗(yàn)來流速度V∞=70m/s，雷諾數(shù)Re=1.2×106情形。對DLR-F12 標(biāo)準(zhǔn)模型進(jìn)行非結(jié)構(gòu)網(wǎng)格的劃分，其中Euler方程計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共416598個三角形，計(jì)算域共5871138個四面體，RANS方程計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共232108個三角形，計(jì)算域共7263616個三棱柱和2746623個四面體。計(jì)算條件與試驗(yàn)條件相同，即馬赫數(shù)Ma=0.20597，雷諾數(shù)Re=1.2×106。采用Euler方程和RANS方程進(jìn)行氣動力和滾轉(zhuǎn)、俯仰和偏航阻尼動導(dǎo)數(shù)CFD 計(jì)算，其中RANS 方程采用剪切應(yīng)力輸運(yùn)（SST）湍流模型。DLRF12標(biāo)準(zhǔn)模型的Euler表面網(wǎng)格劃分如圖6所示。

圖6 DLR-F12模型Euler表面網(wǎng)格劃分及其局部放大示意圖Fig.6 Euler surface mesh division and local zoom-in of DLR-F12 model

本文分別對改進(jìn)方法和參考文獻(xiàn)[10]的方法進(jìn)行CFD求解，獲得定常氣動力以及滾轉(zhuǎn)、俯仰和偏航阻尼動導(dǎo)數(shù)。

2.2 定常氣動力計(jì)算結(jié)果

采用以上模型進(jìn)行Euler 方程和RANS 方程的定常氣動力計(jì)算，并與風(fēng)洞試驗(yàn)數(shù)據(jù)及其他CFD計(jì)算結(jié)果進(jìn)行對比。計(jì)算迎角為-5°～8°（間隔為1°）。計(jì)算結(jié)果及其與風(fēng)洞試驗(yàn)和其他CFD 計(jì)算結(jié)果的升力系數(shù)（CL）和俯仰力矩系數(shù)（Cm）對比如圖7所示。圖中風(fēng)洞試驗(yàn)數(shù)據(jù)和其他CFD計(jì)算結(jié)果由參考文獻(xiàn)[22]提供。

圖7 定常氣動力計(jì)算結(jié)果及其與風(fēng)洞試驗(yàn)和其他CFD計(jì)算結(jié)果對比Fig.7 Steady aerodynamic coefficient computation vs wind tunnel test vs other CFD results

從計(jì)算結(jié)果可以看出，升力系數(shù)與風(fēng)洞試驗(yàn)和其他CFD計(jì)算結(jié)果基本吻合，Euler方程俯仰力矩系數(shù)與風(fēng)洞試驗(yàn)結(jié)果存在一個截距的差異，RANS 方程俯仰力矩系數(shù)在迎角較小時與風(fēng)洞試驗(yàn)結(jié)果存在一定差異，在迎角較大時與風(fēng)洞試驗(yàn)結(jié)果較為接近，與其他CFD 計(jì)算結(jié)果基本吻合。定常氣動力計(jì)算結(jié)果驗(yàn)證了本文使用的CFD 計(jì)算程序的有效性。

2.3 動導(dǎo)數(shù)計(jì)算結(jié)果

2.3.1 不同旋轉(zhuǎn)角速度的動導(dǎo)數(shù)計(jì)算結(jié)果

以滾轉(zhuǎn)和偏航阻尼動導(dǎo)數(shù)為例，采用不同的迎角和不同的旋轉(zhuǎn)角速度計(jì)算滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)，最大旋轉(zhuǎn)角速度選用風(fēng)洞試驗(yàn)時出現(xiàn)的最大旋轉(zhuǎn)角速度。以Euler 方程為例，不同滾轉(zhuǎn)角速度（p）下滾轉(zhuǎn)阻尼動導(dǎo)數(shù)（Clp）的計(jì)算結(jié)果以及不同偏航角速度（r）下偏航阻尼動導(dǎo)數(shù)（Cnr）的計(jì)算結(jié)果如圖8所示。

圖8 不同角速度下滾轉(zhuǎn)阻尼和偏航阻尼動導(dǎo)數(shù)結(jié)果Fig.8 Roll damping derivatives and yaw damping derivatives with different angular velocity

可以看出，在旋轉(zhuǎn)角速度不大于風(fēng)洞試驗(yàn)時出現(xiàn)最大旋轉(zhuǎn)角速度時，滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)隨旋轉(zhuǎn)角速度變化較小，且不存在明顯的增減趨勢，因此，可采用多個不同角速度的阻尼動導(dǎo)數(shù)計(jì)算結(jié)果取平均值得到最終的阻尼動導(dǎo)數(shù)。后續(xù)阻尼動導(dǎo)數(shù)的計(jì)算均采用不同旋轉(zhuǎn)角速度的阻尼動導(dǎo)數(shù)取平均值獲得。

2.3.2 參考文獻(xiàn)[10]計(jì)算方法的結(jié)果及其處理

為驗(yàn)證本文提出的方法的正確性，對參考文獻(xiàn)[10]提出的方法進(jìn)行CFD計(jì)算，將計(jì)算結(jié)果進(jìn)行處理并和改進(jìn)方法對比。以滾轉(zhuǎn)阻尼動導(dǎo)數(shù)為例，對于非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)的計(jì)算，參考文獻(xiàn)[10]采用的方法是計(jì)算繞風(fēng)軸系x軸滾轉(zhuǎn)的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)及交叉動導(dǎo)數(shù)，同時計(jì)算繞風(fēng)軸系z軸偏航的偏航阻尼動導(dǎo)數(shù)及交叉動導(dǎo)數(shù)，使用式（3）進(jìn)行計(jì)算結(jié)果的處理，求得體軸系下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)。以滾轉(zhuǎn)角速度p=91.61(°)/s 為例，參考文獻(xiàn)[10]計(jì)算方法的繞風(fēng)軸滾轉(zhuǎn)力矩系數(shù)（Cl）和偏航力矩系數(shù)（Cn）計(jì)算結(jié)果及其繞體軸動導(dǎo)數(shù)處理結(jié)果分別見表1和表2。

表1 參考文獻(xiàn)［10］方法動導(dǎo)數(shù)計(jì)算結(jié)果Table 1 Dynamic derivatives computation results using ref.［10］method

表2 參考文獻(xiàn)［10］方法動導(dǎo)數(shù)處理結(jié)果Table 2 Dynamic derivatives process results using ref.［10］method

類似地，其他狀態(tài)下的阻尼動導(dǎo)數(shù)也可使用以上方法處理得到。

2.3.3 計(jì)算結(jié)果對比

將本文計(jì)算方法的計(jì)算結(jié)果、風(fēng)洞試驗(yàn)結(jié)果、參考文獻(xiàn)[10]方法的計(jì)算結(jié)果和其他CFD計(jì)算結(jié)果進(jìn)行對比。滾轉(zhuǎn)阻尼動導(dǎo)數(shù)、升力系數(shù)對俯仰角速度q動導(dǎo)數(shù)CLq、俯仰阻尼動導(dǎo)數(shù)Cmq和偏航阻尼動導(dǎo)數(shù)對比結(jié)果分別如圖9～圖12所示，圖中風(fēng)洞試驗(yàn)數(shù)據(jù)和其他CFD 計(jì)算結(jié)果由參考文獻(xiàn)[22]提供，其中風(fēng)洞試驗(yàn)數(shù)據(jù)為組合動導(dǎo)數(shù)。參考文獻(xiàn)[22]中未提供改變側(cè)滑角的CLq和Cmq的風(fēng)洞試驗(yàn)數(shù)據(jù)和其他CFD計(jì)算結(jié)果。由于改變迎角的俯仰阻尼動導(dǎo)數(shù)計(jì)算使用參考文獻(xiàn)[10]的方法也可通過兩次計(jì)算得到一個狀態(tài)的動導(dǎo)數(shù)，本文并未進(jìn)行改進(jìn)，因此這一情況只將改進(jìn)方法與風(fēng)洞試驗(yàn)數(shù)據(jù)和其他CFD計(jì)算結(jié)果進(jìn)行對比。

圖9 改進(jìn)方法、參考文獻(xiàn)[10]方法、風(fēng)洞試驗(yàn)和其他CFD計(jì)算的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)結(jié)果對比Fig.9 Roll damping derivatives result and comparison of our method,ref.[10]method,wind tunnel test and other CFD results

圖10 不同迎角下改進(jìn)方法、風(fēng)洞試驗(yàn)和其他CFD計(jì)算的俯仰運(yùn)動主要動導(dǎo)數(shù)結(jié)果對比Fig.10 Main pitch dynamic derivatives result and comparison between our method,wind tunnel test and other CFD results with different angle of attack

圖11 不同側(cè)滑角下改進(jìn)方法和參考文獻(xiàn)[10]方法計(jì)算的俯仰運(yùn)動主要動導(dǎo)數(shù)結(jié)果對比Fig.11 Main pitch dynamic derivatives result comparison of our method and ref.[10]method with different angle of sideslip

圖12 改進(jìn)方法、參考文獻(xiàn)[10]方法、風(fēng)洞試驗(yàn)和其他CFD計(jì)算的偏航阻尼動導(dǎo)數(shù)結(jié)果對比Fig.12 Yaw damping derivatives result and comparison of our method,ref.[10]method,wind tunnel test and other CFD results

可以看出，在滾轉(zhuǎn)、俯仰、偏航三個方向上，無論是Euler方程還是RANS方程，改進(jìn)方法的計(jì)算結(jié)果均與參考文獻(xiàn)[10]方法的計(jì)算結(jié)果高度吻合。此外，改進(jìn)方法計(jì)算結(jié)果均與其他CFD計(jì)算結(jié)果相差不大，滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)基本吻合，但偏航阻尼動導(dǎo)數(shù)的RANS方程計(jì)算結(jié)果隨迎角變化的趨勢與風(fēng)洞試驗(yàn)結(jié)果存在一定差異。俯仰運(yùn)動的升力系數(shù)對俯仰角速度動導(dǎo)數(shù)和俯仰阻尼動導(dǎo)數(shù)計(jì)算結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)存在系統(tǒng)性偏差，但計(jì)算結(jié)果隨迎角的變化趨勢與風(fēng)洞試驗(yàn)相同。本文主要針對參考文獻(xiàn)[10]的阻尼動導(dǎo)數(shù)計(jì)算方法進(jìn)行改進(jìn)，風(fēng)洞試驗(yàn)和CFD的系統(tǒng)性偏差是一個專門的研究領(lǐng)域，不在本文研究范圍內(nèi)。

3 跨聲速巡航標(biāo)模（TCR）動導(dǎo)數(shù)計(jì)算

3.1 計(jì)算模型

跨聲速巡航標(biāo)模（TCR）是由瑞典薩伯（SAAB）公司研制的標(biāo)準(zhǔn)空氣動力學(xué)模型，機(jī)翼參考面積S=0.3056m2，平均氣動弦長c=0.2943m，展長b=1.1165m，力矩參考點(diǎn)x坐標(biāo)xM=0.87475m。該模型在俄羅斯中央空氣流體力學(xué)研究院（TsAGI）進(jìn)行了一系列風(fēng)洞試驗(yàn)[23]，風(fēng)洞試驗(yàn)來流速度V∞=40m/s，雷諾數(shù)Re=7.9×105。對TCR標(biāo)準(zhǔn)模型進(jìn)行非結(jié)構(gòu)網(wǎng)格的劃分，其中Euler方程定常計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共172550個三角形，計(jì)算域共2280959個四面體，動導(dǎo)數(shù)計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共159740 個三角形，計(jì)算域共2046552個四面體；RANS方程定常計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共172382 個三角形，計(jì)算域共5516224 個三棱柱和1979151 個四面體，動導(dǎo)數(shù)計(jì)算網(wǎng)格機(jī)身和機(jī)翼表面共138814個三角形，計(jì)算域共4442048個三棱柱和3074231個四面體。計(jì)算條件采用與試驗(yàn)條件相同的計(jì)算條件，即馬赫數(shù)Ma=0.1179，雷諾數(shù)Re=7.9×105。采用Euler 方程和RANS 方程進(jìn)行滾轉(zhuǎn)阻尼動導(dǎo)數(shù)CFD 計(jì)算，其中RANS 方程采用SST湍流模型。TCR標(biāo)準(zhǔn)模型的表面網(wǎng)格劃分及其局部放大如圖13所示。

圖13 TCR模型表面網(wǎng)格劃分及其局部放大示意圖Fig.13 Euler surface mesh division and local zoom-in of TCR model

參考文獻(xiàn)[24]中已采用TCR 模型進(jìn)行Euler 方程和RANS方程的定常氣動力計(jì)算，本文采用與參考文獻(xiàn)[24]相同的計(jì)算軟件和計(jì)算網(wǎng)格，分別對改進(jìn)方法和參考文獻(xiàn)[10]的方法進(jìn)行CFD 求解，獲得滾轉(zhuǎn)、俯仰和偏航阻尼動導(dǎo)數(shù)。

3.2 定常計(jì)算結(jié)果

參考文獻(xiàn)[24]的定常計(jì)算結(jié)果及其與風(fēng)洞試驗(yàn)的法向力系數(shù)（CN）和俯仰力矩系數(shù)對比如圖14所示。

圖14 氣動力系數(shù)的定常計(jì)算結(jié)果及其與風(fēng)洞試驗(yàn)對比Fig.14 Steady computation of aerodynamic coefficient and comparison with wind tunnel test

結(jié)果顯示，法向力系數(shù)計(jì)算結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)高度吻合，俯仰力矩系數(shù)計(jì)算結(jié)果在-6°～2°迎角工況下與風(fēng)洞試驗(yàn)數(shù)據(jù)基本吻合，在迎角≥4°情況下存在系統(tǒng)性偏差。由于本文主要將改進(jìn)的阻尼動導(dǎo)數(shù)計(jì)算方法和參考文獻(xiàn)[10]的阻尼動導(dǎo)數(shù)計(jì)算方法進(jìn)行對比，因此計(jì)算結(jié)果滿足本文研究要求。

3.3 動導(dǎo)數(shù)計(jì)算結(jié)果

采用2.3.1和2.3.2節(jié)介紹的處理方法，對改進(jìn)方法和參考文獻(xiàn)[10]方法的計(jì)算數(shù)據(jù)進(jìn)行處理。將改進(jìn)方法的計(jì)算結(jié)果、參考文獻(xiàn)[10]方法的計(jì)算結(jié)果及風(fēng)洞試驗(yàn)結(jié)果進(jìn)行對比。滾轉(zhuǎn)阻尼動導(dǎo)數(shù)、法向力系數(shù)動導(dǎo)數(shù)CNq、俯仰阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)對比結(jié)果分別如圖15～圖17 所示，其中滾轉(zhuǎn)阻尼動導(dǎo)數(shù)、偏航阻尼動導(dǎo)數(shù)的風(fēng)洞試驗(yàn)數(shù)據(jù)為組合動導(dǎo)數(shù)。

圖15 改進(jìn)方法、參考文獻(xiàn)[10]方法滾轉(zhuǎn)阻尼動導(dǎo)數(shù)結(jié)果和風(fēng)洞試驗(yàn)結(jié)果對比Fig.15 Roll damping derivatives result comparison between our method,ref.[10]method and wind tunnel test

圖16 改進(jìn)方法和文獻(xiàn)[10]方法法向力系數(shù)相對俯仰角速度動導(dǎo)數(shù)和俯仰阻尼動導(dǎo)數(shù)結(jié)果對比Fig.16 Pitch rate-normal force coefficient derivatives and pitch damping derivatives result comparison between our method and ref.[10]method

圖17 改進(jìn)方法、參考文獻(xiàn)[10]方法偏航阻尼動導(dǎo)數(shù)結(jié)果和風(fēng)洞試驗(yàn)結(jié)果對比Fig.17 Yaw damping derivatives result comparison between our method,ref.[10]method and wind tunnel test

可以看出，在計(jì)算TCR 標(biāo)準(zhǔn)模型時，改進(jìn)方法的計(jì)算結(jié)果在滾轉(zhuǎn)、俯仰、偏航三個方向上均與參考文獻(xiàn)[10]方法高度吻合。本文的計(jì)算方法計(jì)算得到的動導(dǎo)數(shù)和風(fēng)洞試驗(yàn)數(shù)據(jù)存在系統(tǒng)性偏差，現(xiàn)象為風(fēng)洞試驗(yàn)數(shù)據(jù)值圍繞計(jì)算結(jié)果跳躍。造成差異可能的原因有風(fēng)洞試驗(yàn)結(jié)果包括了側(cè)滑角變化率的影響，以及風(fēng)洞試驗(yàn)和計(jì)算的精度原因。

4 結(jié)論

本文對基于旋轉(zhuǎn)參考坐標(biāo)系的阻尼動導(dǎo)數(shù)CFD 準(zhǔn)定常計(jì)算方法進(jìn)行了改進(jìn)，并使用歐洲動導(dǎo)數(shù)標(biāo)模（DLRF12）和跨聲速巡航標(biāo)模（TCR）進(jìn)行驗(yàn)證。主要結(jié)論如下：

（1）該方法直接計(jì)算滾轉(zhuǎn)阻尼動導(dǎo)數(shù)時考慮了速度相對旋轉(zhuǎn)軸存在垂直分量的影響，計(jì)算任意旋轉(zhuǎn)方向和任意來流方向的一個計(jì)算狀態(tài)的動導(dǎo)數(shù)均只需兩次CFD計(jì)算，相比參考文獻(xiàn)[10]的旋轉(zhuǎn)坐標(biāo)系CFD計(jì)算方法可以減少一半的計(jì)算量。

（2）該方法計(jì)算非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)等特殊情況下，均與參考文獻(xiàn)[10]的計(jì)算結(jié)果高度吻合，可以替代參考文獻(xiàn)[10]方法計(jì)算非零迎角下的滾轉(zhuǎn)阻尼動導(dǎo)數(shù)等特殊情況。

（3）該方法的計(jì)算結(jié)果均與其他CFD計(jì)算結(jié)果基本一致。DLR-F12 標(biāo)模滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)與風(fēng)洞試驗(yàn)基本吻合，俯仰阻尼動導(dǎo)數(shù)與風(fēng)洞試驗(yàn)存在系統(tǒng)性偏差；由于風(fēng)洞試驗(yàn)結(jié)果包括了側(cè)滑角變化率的影響，以及風(fēng)洞試驗(yàn)和計(jì)算的精度等可能的原因，TCR 標(biāo)模動導(dǎo)數(shù)與風(fēng)洞試驗(yàn)存在偏差。

（4）該方法可用于Euler方程和RANS方程CFD計(jì)算，既適用于飛機(jī)概念設(shè)計(jì)階段快速迭代的計(jì)算需求，也滿足飛機(jī)初步設(shè)計(jì)和詳細(xì)設(shè)計(jì)階段高精度計(jì)算需求。