基于CFD方法對串列翼飛行器動導(dǎo)數(shù)的計算

楊起帆 郭家寧 周兵 王際洲

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

關(guān)鍵詞：動導(dǎo)數(shù)；CFD；串列翼；非常規(guī)布局

文獻標識碼：A

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

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

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

1?控制方程與數(shù)值方法

1.1?控制方程

計算流體力學(xué)中對流動的控制方程為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?數(shù)值方法

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

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

2?動導(dǎo)數(shù)的計算方法

對動導(dǎo)數(shù)的求解主要使用強迫俯仰振動法與差分法相結(jié)合。

2.1?強迫俯仰振動法

基于小擾動簡化假設(shè)與氣動力線化解耦［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，當(dāng)?shù)厮矫媾c飛機參考軸之間的夾角θ，升降舵偏轉(zhuǎn)角δe與發(fā)動機操縱參數(shù)δt的變化量均為0。同時對于以ω為頻率作小幅俯仰運動的飛行器而言，其俯仰角速度等同于其迎角變化率。因此將式（2）簡化為如下形式：

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

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

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

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

k=ωlref2VSymboleB@

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

Cm=Cm0+ΔCm（4）

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

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

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

2.2?差分法

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

Cm=Cm0+Cmqq-（6）

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

q-=ql2V

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

q-=ql2qR=l2R

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

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

3.1?參數(shù)設(shè)置與網(wǎng)格劃分

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

3.2?俯仰組合導(dǎo)數(shù)求解

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

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

3.3?俯仰阻尼導(dǎo)數(shù)求解

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

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

4?有翼導(dǎo)彈動導(dǎo)數(shù)計算

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

4.1?計算參數(shù)與算例選取

Finner導(dǎo)彈模型及表面網(wǎng)格如圖7所示：

選取算例中以下參數(shù)均與原文［8］參數(shù)保持一致，馬赫數(shù)M∞=2.1，減縮頻率k=0.00316，xcg=5d，其中d為彈體直徑。網(wǎng)格采用非結(jié)構(gòu)化網(wǎng)格，總網(wǎng)格量為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)的每時間步長設(shè)為0.01秒，最大迭代次數(shù)為50次，收斂條件限制連續(xù)性方程殘差小于1×10-6。

4.2?俯仰組合動導(dǎo)數(shù)計算結(jié)果

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

計算結(jié)果如表2所示：

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

結(jié)語

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

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

參考文獻：

［1］班度·N.帕瑪?shù)?飛機的性能、穩(wěn)定性、動力學(xué)與控制［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技術(shù)計算飛行器動導(dǎo)數(shù)［J］.北京航空航天大學(xué)學(xué)報，2013（2）：196200.

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

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

［6］小約翰·D.安德森.計算流體力學(xué)：基礎(chǔ)與應(yīng)用［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—?），男，漢族，浙江嘉興人，研究方向：計算流體力學(xué)。