楊起帆 郭家寧 周兵 王際洲

摘?要：本文基于計算流體力學方法，通過強迫俯仰振動法及差分法對串列翼布局的小型低速飛行器進行了動導數的計算研究，又使用該方法對某串列翼布局的小型無人飛行器進行了縱向動導數的計算，驗證了CFD方法計算復雜外形飛行器動導數的可行性。本文選取了一例動導數標準模型Basic?Finner的算例進行計算，并與文獻數據進行了比較，計算結果表明，CFD方法估算動導數具有較高的可靠性與工程精度，進一步驗證了該方法的準確性與可靠性。串列翼布局的飛行器縱向動導數計算結果表明，CFD方法可以較好地對串列翼布局飛行器或其他非常規(guī)布局的飛行器進行動導數估算，為工程應用提供了可選解決方案。

關鍵詞：動導數；CFD；串列翼；非常規(guī)布局

文獻標識碼：A

飛行器動導數是指飛行器所受的氣動力與氣動力矩對飛行器運動參數的導數，如俯仰力矩系數隨俯仰角變化的動導數。飛行器動導數是求解小擾動方程特征根不可或缺的基本參數，而飛行器不同模態(tài)的阻尼比、自然頻率、周期以及倍幅時間都是由該模態(tài)的特征根決定的，因此飛行器的動導數直接決定了飛行器的飛行品質及響應。過去獲取動導數的方法包括飛行試驗、風洞試驗或工程估算。對于小型低速無人機而言，飛行試驗與風洞試驗周期太長，成本過高。而工程估算方法以一些經驗或半經驗公式為主導對飛行器的動導數進行估算，但是這些經驗或者半經驗公式都是基于常規(guī)布局發(fā)展而來，因此對于非常規(guī)布局的飛行器進行估算時會有較大的誤差，甚至產生錯誤的結果［1］。對于串列翼布局的飛行器而言，初步設計軟件如XFLR5就難以給出滿足需求的結果。

計算流體力學（CFD，Computational?Fluid?Dynamics）作為一門流體力學與數學的交叉學科，是通過數值方法來求解流體力學的控制方程。隨著近現代計算機的高速發(fā)展，使用計算流體力學作為手段來估算飛行器氣動特性已得到廣泛應用。CFD手段與傳統(tǒng)風洞試驗相比有成本低、速度快等優(yōu)點，在設計階段可以快速迭代優(yōu)化或者進行氣動外形的再設計。

來自NASA蘭利研究中心的Lawrence?L.Green等人［2］最早提出可以使用CFD手段計算飛行器的動導數。葉川、馬東立［3］使用強迫俯仰振動法與差分法對Finner導彈標準模型進行了計算，驗證了該方法的準確性。米百剛、詹浩與朱軍［4］在強迫俯仰振動法與差分法之外，進一步研究了準定常方法求解動導數。而孫濤、高正紅、黃江濤［5］研究了在使用強迫俯仰振動法與差分法求解動導數過程中減縮頻率對計算結果的影響。本文使用精度較高的強迫俯仰振動法與差分法對某串列翼布局的小型低速傾轉旋翼無人機進行了動導數計算。

1?控制方程與數值方法

1.1?控制方程

計算流體力學中對流動的控制方程為NS方程［6］。三維NS方程的守恒形式如下：

Ut+Fx+Gy+Hz=J（1）

其中：

U=ρ

ρu

ρv

ρw

ρ（e+V22）

F=ρu

ρu2+p－τxx

ρvu－τxy

ρvu－τxz

ρe+V22u+pu－kTx－uτxx－vτxy－wτxz

G=ρv

ρuv－τyx

ρv2+p－τyy

ρwv－τyz

ρe+V22v+pv－kTy－uτyx－vτyy－wτyz

H=ρw

ρuw－τzx

ρvw－τzy

ρw2+p－τzz

ρe+V22w+pw－kTz－uτzx－vτzy－wτzz

J=0

ρfx

ρfy

ρfz

ρufx+vfy+wfz+ρq·

其中列矢量U為流動的原始變量項，列矢量FGH為通量項，列矢量J代表源項。

1.2?數值方法

考慮到需要求解的飛行器設計巡航速度為0.1馬赫，因此求解過程中可以將流動近似為不可壓縮流動?？臻g離散采用有限體積法，數值計算采用針對不可壓流動的SIMPLE算法，其中動量方程與壓力方程都使用二階迎風格式。時間離散采用二階向后歐拉格式。湍流模型使用kωSST（剪切應力輸運）雙方程模型。

而對于動導數標準模型finner導彈，由于其超聲速的特性，使用針對可壓縮流動的Roe通量分裂格式，空間上使用二階迎風格式。

2?動導數的計算方法

對動導數的求解主要使用強迫俯仰振動法與差分法相結合。

2.1?強迫俯仰振動法

基于小擾動簡化假設與氣動力線化解耦［1］，受擾動后的氣動力與力矩形式為（略去高階項）：

ΔCm=Cmuu+CmαΔα+CmθΔθ+Cmα·Δα·+

Cmqq+CmδeΔδe+CmδtΔδt+…（2）

對于受擾動的飛行狀態(tài)，式中u為無量綱化的來流速度變化量，俯仰角速度q為俯仰角變化量

u=ΔUU0??q=Δq

對于巡航狀態(tài)下的飛行器，來流速度u，當地水平面與飛機參考軸之間的夾角θ，升降舵偏轉角δe與發(fā)動機操縱參數δt的變化量均為0。同時對于以ω為頻率作小幅俯仰運動的飛行器而言，其俯仰角速度等同于其迎角變化率。因此將式（2）簡化為如下形式：

ΔCm=CmαΔα+Cmα·Δα·+CmqΔα·（3）

對于強迫俯仰運動，飛行器做周期性小幅俯仰運動，因此定義迎角變化如下：

Δα=αampsin（ωt）??Δα·=kαampcos（ωt）

其中αamp為振幅，k為減縮頻率：

k=ωlref2VSymboleB@

俯仰力矩系數隨時間變化的公式如下：

Cm=Cm0+ΔCm（4）

其中Cm0為無擾動情況下的俯仰力矩系數，ΔCm即為式（3）：

Cm=Cm0+Cmααampsin（ωt）+（Cmα·+Cmq）kαampcos（ωt）（5）

通過CFD方法可以求得俯仰力矩系數（Cm）隨時間變化的曲線，忽略其中受到初始條件影響的部分，對周期性變化的曲線部分進行最小二乘法回歸即可得到式（5）中的三個系數：無擾動下的俯仰力矩系數Cm0、俯仰力矩系數隨攻角變化的靜穩(wěn)定性導數Cmα以及俯仰組合導數（Cmα·+Cmq）。

2.2?差分法

使飛行器進行定常圓周運動，以定常拉升運動為例，改變飛行器的俯仰角但使其迎角保持不變。因此Δα與Δα·為0，簡化式（2）并代入式（4）可得：

Cm=Cm0+Cmqq-（6）

其中q-為無量綱化俯仰角速度：

q-=ql2V

對于圓周運動而言，線速度為角速度乘以半徑，V=qR。故：

q-=ql2qR=l2R

對半徑R選取不同的值進行多次CFD定常計算即可求得俯仰力矩系數對俯仰角速度的動導數Cmq。同時也可求得2.1節(jié)中俯仰組合導數中的俯仰力矩系數對迎角變化率的動導數Cmα·。

3?對某串列翼布局飛行器的計算

3.1?參數設置與網格劃分

該串列翼無人機幾何尺寸如圖1所示。網格劃分采用非結構化網格，無人機表面網格如圖2所示，總網格數為190萬。根據Thomas?H.Pulliam［7］的研究，對算例進行網格無關性與遠場無關性進行驗證。對于該低速低空無人機選取來流M∞=0.1，初始攻角α0=0°，減縮頻率k=0049，振幅αamp=2°，重心為俯仰平面內距離機頭前緣三倍弦長處，遠場距離飛行器40倍弦長處，采用滑移網格模擬飛行器的周期性振動。

3.2?俯仰組合導數求解

每時間步長設為0.01秒，最大迭代次數為50次，收斂條件限制連續(xù)性方程殘差小于1×10-6。求得圖3俯仰力矩系數隨時間變化的曲線。

從上圖可以看到俯仰力矩系數與瞬時迎角之間存在遲滯效應，這是非定常計算中出現的典型現象。圖4為該串列翼布局飛行器的遲滯環(huán)。可以看到部分計算結果受到了初始條件的影響，因此對于圖3中的俯仰力矩曲線截取不受初始條件影響，呈周期性變化的部分并進行前文2.1節(jié)所提及的系數回歸可得下文表1中的無擾動俯仰力矩系數Cm0，俯仰靜穩(wěn)定導數Cmα和俯仰組合導數（Cmα·+Cmq）。

3.3?俯仰阻尼導數求解

為了求解俯仰阻尼導數，需要對作勻速圓周運動的飛行器進行計算，采用多重參考系法，計算結果如圖5所示：

根據2.2節(jié)的推導，無量綱俯仰角速度q-只與旋轉半徑R有關。但是為了保證飛行器的飛行速度V不發(fā)生改變，因此在改變半徑進行計算時也需要改變俯仰角速度q與之匹配。計算結果圖6所示，通過求解圖6曲線的斜率可求得俯仰力矩系數對俯仰角速率的動導數Cmq，同時也可從3.2節(jié)中所求得的俯仰組合導數中解耦求得俯仰力矩系數對迎角變化率的動導數Cmα·。

4?有翼導彈動導數計算

為驗證CFD方法的可靠性，對國際動導數標準模型finner導彈選取算例進行計算，計算模型、計算參數及計算結果的對比參照Erdal?Oktay［8］等人的報告。

4.1?計算參數與算例選取

Finner導彈模型及表面網格如圖7所示：

選取算例中以下參數均與原文［8］參數保持一致，馬赫數M∞=2.1，減縮頻率k=0.00316，xcg=5d，其中d為彈體直徑。網格采用非結構化網格，總網格量為203萬，計算遠場位于45d，計算模型采用sstkw湍流模型，隱式時間推進，振幅選取αamp=3°。計算非穩(wěn)態(tài)解的方法與原文［8］保持一致，即計算0°攻角下的穩(wěn)態(tài)解，再以穩(wěn)態(tài)解為初始條件計算攻角隨時間變化的非穩(wěn)態(tài)解。非穩(wěn)態(tài)的每時間步長設為0.01秒，最大迭代次數為50次，收斂條件限制連續(xù)性方程殘差小于1×10-6。

4.2?俯仰組合動導數計算結果

俯仰力矩系數隨時間變化曲線如圖8所示：

計算結果如表2所示：

從選取的算例可以看出本文選取的攻角變化速率、攻角變化幅度與原文［8］所使用的攻角變化速率、攻角變化幅度有所不同，依然能較好地預測俯仰靜穩(wěn)定導數與俯仰組合動導數。

結語

本文使用CFD方法對串列翼布局的小型低速無人機進行了縱向的靜穩(wěn)定性導數和一些動導數的計算。同時為了驗證CFD計算結果的可靠性，本文對國際動導數標準模型basic?finner導彈進行了計算，計算結果表明CFD計算方法對獲取飛行器的動導數具有較好的可靠性和較高的工程精度。

因此，對于具有復雜氣動外形或者非常規(guī)布局的飛行器，若不適合使用工程估算、半經驗方法或風洞實驗，則可以在飛行器基本氣動外形確定后使用CFD方法計算這些導數，用于判斷飛行器的飛行品質與響應，以及設計飛機的飛行控制與導航系統(tǒng)。

參考文獻：

［1］班度·N.帕瑪迪.飛機的性能、穩(wěn)定性、動力學與控制［M］.第2版.北京：航空工業(yè)出版社，2013.

［2］Lawrence?L.Green，Angela?M.Spence，Patrick?C.Murphy.Computational?Methods?for?Dynamic?Stability?and?Control?Derivatives［R］.AIAA，20040015?January?2004.

［3］葉川，馬東立.利用CFD技術計算飛行器動導數［J］.北京航空航天大學學報，2013（2）：196200.

［4］米百剛，詹浩，朱軍.基于CFD數值仿真技術的飛行器動導數計算［J］.空氣動力學學報，2014（6）：834839.

［5］孫濤，高正紅，黃江濤.基于CFD的動導數計算及減縮頻率影響分析［J］.飛行力學，2011，29（4）：1518.

［6］小約翰·D.安德森.計算流體力學：基礎與應用［M］.北京：航空工業(yè)出版社，2018.

［7］Thomas?H.Pulliam.Solution?methods?in?Computational?Fluid?Dynamics［R］.Moffett?Field：NASA?Ames?research?center，1986.

［8］Erdal?Oktay.CFD?predictions?of?dynamic?derivatives?for?missiles［R］.AIAA，20020276，2002.

作者簡介：楊起帆（1994—?），男，漢族，浙江嘉興人，研究方向：計算流體力學。