吳耕宇，鮑君波，宋萬強，相倩

中國航空研究院，北京 100012

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

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

1 計算方法

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

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

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

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

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

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

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

圖5 飛機運動速度矢量分解示意圖Fig.5 Aircraft velocity vector decomposition map

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

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

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

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

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

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

2.1 計算模型

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

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

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

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

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

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

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

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

以滾轉(zhuǎn)和偏航阻尼動導(dǎo)數(shù)為例，采用不同的迎角和不同的旋轉(zhuǎn)角速度計算滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)，最大旋轉(zhuǎn)角速度選用風(fēng)洞試驗時出現(xiàn)的最大旋轉(zhuǎn)角速度。以Euler 方程為例，不同滾轉(zhuǎn)角速度（p）下滾轉(zhuǎn)阻尼動導(dǎo)數(shù)（Clp）的計算結(jié)果以及不同偏航角速度（r）下偏航阻尼動導(dǎo)數(shù)（Cnr）的計算結(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)洞試驗時出現(xiàn)最大旋轉(zhuǎn)角速度時，滾轉(zhuǎn)阻尼動導(dǎo)數(shù)和偏航阻尼動導(dǎo)數(shù)隨旋轉(zhuǎn)角速度變化較小，且不存在明顯的增減趨勢，因此，可采用多個不同角速度的阻尼動導(dǎo)數(shù)計算結(jié)果取平均值得到最終的阻尼動導(dǎo)數(shù)。后續(xù)阻尼動導(dǎo)數(shù)的計算均采用不同旋轉(zhuǎn)角速度的阻尼動導(dǎo)數(shù)取平均值獲得。

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

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

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

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

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

2.3.3 計算結(jié)果對比

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

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

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

3.1 計算模型

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

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

3.2 定常計算結(jié)果

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

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

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

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

采用2.3.1和2.3.2節(jié)介紹的處理方法，對改進方法和參考文獻[10]方法的計算數(shù)據(jù)進行處理。將改進方法的計算結(jié)果、參考文獻[10]方法的計算結(jié)果及風(fēng)洞試驗結(jié)果進行對比。滾轉(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)洞試驗數(shù)據(jù)為組合動導(dǎo)數(shù)。

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

圖16 改進方法和文獻[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 改進方法、參考文獻[10]方法偏航阻尼動導(dǎo)數(shù)結(jié)果和風(fēng)洞試驗結(jié)果對比Fig.17 Yaw damping derivatives result comparison between our method,ref.[10]method and wind tunnel test

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

4 結(jié)論

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

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

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

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

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