改進的準定常旋轉坐標系阻尼動導數計算方法

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

中國航空研究院，北京 100012

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

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

1 計算方法

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

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

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

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

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

圖2 在穩定坐標系中勻速滾轉Fig.2 Rolling in stability coordinate frame

圖3 在穩定坐標系中繞軸偏航Fig.3 Yawing in stability coordinate frame

圖4 在機體坐標系中勻速滾轉Fig.4 Rolling in body coordinate frame

1.1 非零迎角下的滾轉阻尼動導數計算

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

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

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

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

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

1.2 任意方向阻尼動導數計算方法

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

2 歐洲動導數標模（DLR-F12）動導數計算

2.1 計算模型

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

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

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

2.2 定常氣動力計算結果

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

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

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

2.3 動導數計算結果

2.3.1 不同旋轉角速度的動導數計算結果

以滾轉和偏航阻尼動導數為例，采用不同的迎角和不同的旋轉角速度計算滾轉阻尼動導數和偏航阻尼動導數，最大旋轉角速度選用風洞試驗時出現的最大旋轉角速度。以Euler 方程為例，不同滾轉角速度（p）下滾轉阻尼動導數（Clp）的計算結果以及不同偏航角速度（r）下偏航阻尼動導數（Cnr）的計算結果如圖8所示。

圖8 不同角速度下滾轉阻尼和偏航阻尼動導數結果Fig.8 Roll damping derivatives and yaw damping derivatives with different angular velocity

可以看出，在旋轉角速度不大于風洞試驗時出現最大旋轉角速度時，滾轉阻尼動導數和偏航阻尼動導數隨旋轉角速度變化較小，且不存在明顯的增減趨勢，因此，可采用多個不同角速度的阻尼動導數計算結果取平均值得到最終的阻尼動導數。后續阻尼動導數的計算均采用不同旋轉角速度的阻尼動導數取平均值獲得。

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

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

表1 參考文獻［10］方法動導數計算結果Table 1 Dynamic derivatives computation results using ref.［10］method

表2 參考文獻［10］方法動導數處理結果Table 2 Dynamic derivatives process results using ref.［10］method

類似地，其他狀態下的阻尼動導數也可使用以上方法處理得到。

2.3.3 計算結果對比

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

圖9 改進方法、參考文獻[10]方法、風洞試驗和其他CFD計算的滾轉阻尼動導數結果對比Fig.9 Roll damping derivatives result and comparison of our method,ref.[10]method,wind tunnel test and other CFD results

圖10 不同迎角下改進方法、風洞試驗和其他CFD計算的俯仰運動主要動導數結果對比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 不同側滑角下改進方法和參考文獻[10]方法計算的俯仰運動主要動導數結果對比Fig.11 Main pitch dynamic derivatives result comparison of our method and ref.[10]method with different angle of sideslip

圖12 改進方法、參考文獻[10]方法、風洞試驗和其他CFD計算的偏航阻尼動導數結果對比Fig.12 Yaw damping derivatives result and comparison of our method,ref.[10]method,wind tunnel test and other CFD results

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

3 跨聲速巡航標模（TCR）動導數計算

3.1 計算模型

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

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

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

3.2 定常計算結果

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

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

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

3.3 動導數計算結果

采用2.3.1和2.3.2節介紹的處理方法，對改進方法和參考文獻[10]方法的計算數據進行處理。將改進方法的計算結果、參考文獻[10]方法的計算結果及風洞試驗結果進行對比。滾轉阻尼動導數、法向力系數動導數CNq、俯仰阻尼動導數和偏航阻尼動導數對比結果分別如圖15～圖17 所示，其中滾轉阻尼動導數、偏航阻尼動導數的風洞試驗數據為組合動導數。

圖15 改進方法、參考文獻[10]方法滾轉阻尼動導數結果和風洞試驗結果對比Fig.15 Roll damping derivatives result comparison between our method,ref.[10]method and wind tunnel test

圖16 改進方法和文獻[10]方法法向力系數相對俯仰角速度動導數和俯仰阻尼動導數結果對比Fig.16 Pitch rate-normal force coefficient derivatives and pitch damping derivatives result comparison between our method and ref.[10]method

圖17 改進方法、參考文獻[10]方法偏航阻尼動導數結果和風洞試驗結果對比Fig.17 Yaw damping derivatives result comparison between our method,ref.[10]method and wind tunnel test

可以看出，在計算TCR 標準模型時，改進方法的計算結果在滾轉、俯仰、偏航三個方向上均與參考文獻[10]方法高度吻合。本文的計算方法計算得到的動導數和風洞試驗數據存在系統性偏差，現象為風洞試驗數據值圍繞計算結果跳躍。造成差異可能的原因有風洞試驗結果包括了側滑角變化率的影響，以及風洞試驗和計算的精度原因。

4 結論

本文對基于旋轉參考坐標系的阻尼動導數CFD 準定常計算方法進行了改進，并使用歐洲動導數標模（DLRF12）和跨聲速巡航標模（TCR）進行驗證。主要結論如下：

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

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

（3）該方法的計算結果均與其他CFD計算結果基本一致。DLR-F12 標模滾轉阻尼動導數和偏航阻尼動導數與風洞試驗基本吻合，俯仰阻尼動導數與風洞試驗存在系統性偏差；由于風洞試驗結果包括了側滑角變化率的影響，以及風洞試驗和計算的精度等可能的原因，TCR 標模動導數與風洞試驗存在偏差。

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