楊起帆 郭家寧 周兵 王際洲

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

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

文獻(xiàn)標(biāo)識碼：A

飛行器動導(dǎo)數(shù)是指飛行器所受的氣動力與氣動力矩對飛行器運(yùn)動參數(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ù)的方法包括飛行試驗(yàn)、風(fēng)洞試驗(yàn)或工程估算。對于小型低速無人機(jī)而言，飛行試驗(yàn)與風(fēng)洞試驗(yàn)周期太長，成本過高。而工程估算方法以一些經(jīng)驗(yàn)或半經(jīng)驗(yàn)公式為主導(dǎo)對飛行器的動導(dǎo)數(shù)進(jìn)行估算，但是這些經(jīng)驗(yàn)或者半經(jīng)驗(yàn)公式都是基于常規(guī)布局發(fā)展而來，因此對于非常規(guī)布局的飛行器進(jìn)行估算時會有較大的誤差，甚至產(chǎn)生錯誤的結(jié)果［1］。對于串列翼布局的飛行器而言，初步設(shè)計(jì)軟件如XFLR5就難以給出滿足需求的結(jié)果。

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

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

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

1.1?控制方程

計(jì)算流體力學(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為流動的原始變量項(xiàng)，列矢量FGH為通量項(xiàng)，列矢量J代表源項(xiàng)。

1.2?數(shù)值方法

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

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

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

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

2.1?強(qiáng)迫俯仰振動法

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

Δ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飛機(jī)參考軸之間的夾角θ，升降舵偏轉(zhuǎn)角δe與發(fā)動機(jī)操縱參數(shù)δt的變化量均為0。同時對于以ω為頻率作小幅俯仰運(yùn)動的飛行器而言，其俯仰角速度等同于其迎角變化率。因此將式（2）簡化為如下形式：

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

對于強(qiáng)迫俯仰運(yùn)動，飛行器做周期性小幅俯仰運(yùn)動，因此定義迎角變化如下：

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

2.2?差分法

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

Cm=Cm0+Cmqq-（6）

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

q-=ql2V

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

q-=ql2qR=l2R

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

3?對某串列翼布局飛行器的計(jì)算

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

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

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

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

從上圖可以看到俯仰力矩系數(shù)與瞬時迎角之間存在遲滯效應(yīng)，這是非定常計(jì)算中出現(xiàn)的典型現(xiàn)象。圖4為該串列翼布局飛行器的遲滯環(huán)?？梢钥吹讲糠钟?jì)算結(jié)果受到了初始條件的影響，因此對于圖3中的俯仰力矩曲線截取不受初始條件影響，呈周期性變化的部分并進(jìn)行前文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ù)，需要對作勻速圓周運(yùn)動的飛行器進(jìn)行計(jì)算，采用多重參考系法，計(jì)算結(jié)果如圖5所示：

根據(jù)2.2節(jié)的推導(dǎo)，無量綱俯仰角速度q-只與旋轉(zhuǎn)半徑R有關(guān)。但是為了保證飛行器的飛行速度V不發(fā)生改變，因此在改變半徑進(jìn)行計(jì)算時也需要改變俯仰角速度q與之匹配。計(jì)算結(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ù)計(jì)算

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

4.1?計(jì)算參數(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萬，計(jì)算遠(yuǎn)場位于45d，計(jì)算模型采用sstkw湍流模型，隱式時間推進(jìn)，振幅選取αamp=3°。計(jì)算非穩(wěn)態(tài)解的方法與原文［8］保持一致，即計(jì)算0°攻角下的穩(wěn)態(tài)解，再以穩(wěn)態(tài)解為初始條件計(jì)算攻角隨時間變化的非穩(wěn)態(tài)解。非穩(wěn)態(tài)的每時間步長設(shè)為0.01秒，最大迭代次數(shù)為50次，收斂條件限制連續(xù)性方程殘差小于1×10-6。

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

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

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

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

結(jié)語

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

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

參考文獻(xiàn)：

［1］班度·N.帕瑪?shù)?飛機(jī)的性能、穩(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ù)計(jì)算飛行器動導(dǎo)數(shù)［J］.北京航空航天大學(xué)學(xué)報(bào)，2013（2）：196200.

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

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

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