NACA0012翼型跨聲速強(qiáng)迫運(yùn)動(dòng)非定常氣動(dòng)力模型

張慶， 葉正寅

(1.西安航空學(xué)院 飛行器學(xué)院，陜西 西安 710077； 2.南洋理工大學(xué) 機(jī)械與航空工程學(xué)院，新加坡 639798； 3.西北工業(yè)大學(xué) 航空學(xué)院，陜西 西安 710072)

氣動(dòng)導(dǎo)數(shù)作為描述飛行器機(jī)動(dòng)飛行和受擾動(dòng)時(shí)氣動(dòng)特性的關(guān)鍵性氣動(dòng)參數(shù)，在飛行器氣動(dòng)性能、控制系統(tǒng)和總體設(shè)計(jì)中扮演著非常重要的作用[1-4]。在傳統(tǒng)的飛行動(dòng)力學(xué)相關(guān)問題的研究中，氣動(dòng)力的數(shù)據(jù)往往基于小擾動(dòng)線性疊加原理計(jì)算出來，在這種準(zhǔn)定常假設(shè)情況下，氣動(dòng)力僅僅表示為瞬時(shí)飛行狀態(tài)參數(shù)的函數(shù)，并且可以以一種簡單的解析函數(shù)關(guān)系式表示出來[2-5]。

但是，現(xiàn)代飛行器的飛行包線普遍向大迎角區(qū)域擴(kuò)展，在大迎角下飛機(jī)機(jī)動(dòng)飛行產(chǎn)生的三維非定常分離流和渦流使得空氣動(dòng)力呈現(xiàn)高度非線性特性，氣動(dòng)力和力矩不僅依賴于瞬時(shí)迎角、側(cè)滑角、姿態(tài)角等參數(shù), 而且與它們的時(shí)間歷程有關(guān), 因此原來使用的低階線性疊加模型將不再適用[5-6]。同時(shí)，由于機(jī)動(dòng)飛行狀態(tài)涵蓋了較大的迎角、側(cè)滑角、角速率的變化范圍，如果采用風(fēng)洞實(shí)驗(yàn)或是數(shù)值計(jì)算模擬，其時(shí)間成本和經(jīng)濟(jì)成本都難以接受[7-10]。因此有必要建立起較大飛行包線內(nèi)普適性較好的的非定常氣動(dòng)力模型[1,4]。

Etkin模型是目前動(dòng)導(dǎo)數(shù)求解時(shí)最常用的一種非定常氣動(dòng)力模型，Etkin模型物理意義明確，考慮了時(shí)間歷史效應(yīng)對(duì)氣動(dòng)導(dǎo)數(shù)的影響[3]。但是，在非定常氣動(dòng)力建模時(shí)，該模型中的各項(xiàng)氣動(dòng)導(dǎo)數(shù)對(duì)不同運(yùn)動(dòng)形式的非定常氣動(dòng)力的影響規(guī)律和適用程度尚不清楚。為此，本文結(jié)合Etkin氣動(dòng)力模型，研究了氣動(dòng)力關(guān)于迎角的一階和二階導(dǎo)數(shù)在氣動(dòng)力模型的作用，希望能精確地重構(gòu)出翼型單自由度或是耦合強(qiáng)迫運(yùn)動(dòng)過程中的非定常氣動(dòng)力，為未來發(fā)展高效的、可靠的氣動(dòng)力模型提供參考數(shù)據(jù)。

1 強(qiáng)迫運(yùn)動(dòng)非定常氣動(dòng)力模型

本文的計(jì)算采用課題組自己開發(fā)的柔性體動(dòng)力學(xué)問題求解軟件GMFlow[11-13]，其中流場(chǎng)求解部分采用基于SA模型的有限體積法[13]，強(qiáng)迫運(yùn)動(dòng)時(shí)的網(wǎng)格變形方法為彈簧網(wǎng)格變形方法[14-16]。為了驗(yàn)證求解方法的正確性，首先計(jì)算了標(biāo)準(zhǔn)算例NACA0012翼型強(qiáng)迫俯仰運(yùn)動(dòng)的非定常氣動(dòng)力變化情況，將計(jì)算結(jié)果與文獻(xiàn)中的計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果對(duì)比，對(duì)比結(jié)果見文獻(xiàn)[13]。

俯仰運(yùn)動(dòng)的運(yùn)動(dòng)規(guī)律可以描述為[15]：

α(t)=α0+Asin(ωt)=α0+Asin(2πft)

(1)

式中：α0是初始位置處的迎角；A是簡諧振動(dòng)的振幅；ω是簡諧振動(dòng)的圓頻率；f是簡諧振動(dòng)的頻率。

本文定義減縮頻率為：

(2)

式中C是翼型的弦長。在本文中，強(qiáng)迫運(yùn)動(dòng)時(shí)自由來流的馬赫數(shù)為0.755，翼型弦長為1.0 m，強(qiáng)迫運(yùn)動(dòng)的減縮頻率為0.081 4。俯仰運(yùn)動(dòng)的初始迎角為0.016°，俯仰振幅為2.51°。

圖1(a)是強(qiáng)迫俯仰運(yùn)動(dòng)時(shí)的升力系數(shù)和關(guān)于1/4弦點(diǎn)的俯仰力矩系數(shù)隨時(shí)間的變化曲線，圖中計(jì)算了3個(gè)周期的氣動(dòng)力，由圖可知，在第1個(gè)計(jì)算周期的初始階段，計(jì)算的結(jié)果收斂性較差，這主要是由于定常計(jì)算的步數(shù)不足。從第2個(gè)周期開始，力和力矩系數(shù)已經(jīng)達(dá)到了較好的諧振性，可以認(rèn)為計(jì)算結(jié)果已經(jīng)收斂。因此，為了減小計(jì)算量，本文的所有強(qiáng)迫運(yùn)動(dòng)過程都只計(jì)算了3個(gè)運(yùn)動(dòng)周期。

圖1(b)是翼型強(qiáng)迫沉浮運(yùn)動(dòng)時(shí)的力和力矩系數(shù)變化情況，其運(yùn)動(dòng)規(guī)律為：

z(t)=z0+zmsin(ωt)

(3)

式中：z0=0是初始位置處的縱向位移；zm=0.1 m是沉浮運(yùn)動(dòng)的振幅。

考慮洗流影響，在沉浮運(yùn)動(dòng)的任一時(shí)刻，瞬時(shí)迎角為：

α(t)=α0-ωzmcos(ωt)/V∞

(4)

圖1(c)是翼型強(qiáng)迫俯仰/沉浮耦合運(yùn)動(dòng)時(shí)的力和力矩系數(shù)變化情況，其運(yùn)動(dòng)規(guī)律為式(1)和式(3)疊加。對(duì)比圖1可知，雖然耦合運(yùn)動(dòng)形式是俯仰和沉浮運(yùn)動(dòng)的疊加，但是耦合運(yùn)動(dòng)的氣動(dòng)力和力矩并不等于俯仰運(yùn)動(dòng)和沉浮運(yùn)動(dòng)的簡單疊加，這也說明了翼型強(qiáng)迫運(yùn)動(dòng)時(shí)氣動(dòng)力的非線性遲滯特性比較復(fù)雜，并不是簡單的線性疊加關(guān)系。

圖1 不同運(yùn)動(dòng)過程升力和力矩系數(shù)隨時(shí)間變化Fig.1 History of lift/moment coefficients in different motions

1.1 一階氣動(dòng)模型

根據(jù)Etkin氣動(dòng)力模型[2-3]，強(qiáng)迫運(yùn)動(dòng)過程中的非定常氣動(dòng)力可以表示為：

ΔCj=Cj-Cj0=CjαΔα+

(5)

(6)

俯仰運(yùn)動(dòng)的非定常氣動(dòng)力可以表示為：

(7)

此處需要注意，由于ΔCj是相對(duì)于初始位置的變化量，因此右側(cè)是(cos(ωt)-1)而不是cos(ωt)。所以：

(8)

根據(jù)強(qiáng)迫沉浮運(yùn)動(dòng)時(shí)運(yùn)動(dòng)規(guī)律可知：

(9)

沉浮運(yùn)動(dòng)的非定常氣動(dòng)力可以表示為：

(10)

所以：

(11)

1.2 二階氣動(dòng)模型

根據(jù)Etkin氣動(dòng)力模型[2-3]，非定常氣動(dòng)力可以表示為：

(12)

(13)

由沉浮運(yùn)動(dòng)的氣動(dòng)力變化規(guī)律可以得到：

(14)

2 氣動(dòng)力建模結(jié)果比較

為了定量考察這些氣動(dòng)力模型對(duì)強(qiáng)迫運(yùn)動(dòng)過程非定常遲滯效應(yīng)模擬的適用程度，本節(jié)對(duì)比了這些氣動(dòng)力模型的計(jì)算結(jié)果與直接采用CFD進(jìn)行計(jì)算得到的結(jié)果。

圖2分別是采用一階和二階Etkin氣動(dòng)力模型計(jì)算得到的強(qiáng)迫俯仰運(yùn)動(dòng)、強(qiáng)迫沉浮運(yùn)動(dòng)以及耦合運(yùn)動(dòng)的氣動(dòng)力與采用CFD方法得到的氣動(dòng)力遲滯曲線的對(duì)比圖。由圖2(a)可知，對(duì)于強(qiáng)迫俯仰運(yùn)動(dòng)，采用二階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)與CFD計(jì)算值完全重合，而采用一階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)誤差隨著迎角的增加而增大，在最大迎角位置比CFD計(jì)算值大50%。對(duì)于俯仰力矩系數(shù)，一階模型的誤差較大，而二階模型的結(jié)果與CFD計(jì)算值雖然不像升力系數(shù)那樣完全重合，但是吻合程度也較好。

圖2 不同運(yùn)動(dòng)過程升力和力矩系數(shù)遲滯曲線Fig.2 Comparison of lift/moment coefficients indifferent motions

表1 俯仰運(yùn)動(dòng)不同位置非定常氣動(dòng)力分布情況

圖2(b)是采用一階和二階Etkin氣動(dòng)力模型計(jì)算得到的強(qiáng)迫沉浮運(yùn)動(dòng)的氣動(dòng)力與采用CFD方法得到的氣動(dòng)力遲滯曲線的對(duì)比圖。由圖2(b)可知，對(duì)于強(qiáng)迫沉浮運(yùn)動(dòng)，采用二階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)與CFD計(jì)算值幾乎重合，而采用一階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)誤差較大，在最大縱向位移位置比CFD計(jì)算值大90%。對(duì)于俯仰力矩系數(shù)，一階模型的誤差較大，而二階模型的結(jié)果與CFD計(jì)算值雖然不像升力系數(shù)那樣完全重合，但是吻合程度也較好。

圖2(c)是采用一階和二階Etkin氣動(dòng)力模型計(jì)算得到的強(qiáng)迫俯仰/沉浮耦合運(yùn)動(dòng)的氣動(dòng)力與采用CFD方法得到的氣動(dòng)力遲滯曲線的對(duì)比圖。由圖2(c)可知，對(duì)于強(qiáng)迫耦合運(yùn)動(dòng)，采用二階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)與CFD計(jì)算值完全重合，而采用一階氣動(dòng)導(dǎo)數(shù)得到的升力系數(shù)誤差較大，在最大縱向位移位置比CFD計(jì)算值大145%。對(duì)于俯仰力矩系數(shù)，一階模型的誤差較大，而二階模型的結(jié)果與CFD計(jì)算值雖然不像升力系數(shù)那樣完全重合，但是吻合程度也較好。

表2 沉浮運(yùn)動(dòng)不同位置非定常氣動(dòng)力分布情況

3 結(jié)論

1)不論是強(qiáng)迫俯仰運(yùn)動(dòng)、沉浮運(yùn)動(dòng)，還是俯仰/沉浮耦合運(yùn)動(dòng)，將氣動(dòng)導(dǎo)數(shù)拓展至迎角和俯仰角的二階導(dǎo)數(shù)，都可以十分精確地重現(xiàn)出強(qiáng)迫運(yùn)動(dòng)過程中的非定常升力變化情況。

2)由于俯仰力矩的遲滯曲線并不是簡單的橢圓形，二階模型計(jì)算出的強(qiáng)迫運(yùn)動(dòng)過程的俯仰力矩與CFD計(jì)算值的吻合程度不像升力那么好，說明俯仰力矩的模型要比升力更加復(fù)雜。

3)俯仰/沉浮耦合運(yùn)動(dòng)的非定常氣動(dòng)力并不是俯仰運(yùn)動(dòng)和沉浮運(yùn)動(dòng)的簡單疊加，說明精確的氣動(dòng)力建模還需要深入考慮其他變量的影響。

本文的研究結(jié)果表明，Etkin氣動(dòng)力模型對(duì)于非線性較強(qiáng)的氣動(dòng)力建模仍然具有較好的適用性，但是，對(duì)于三維流動(dòng)以及接近失速迎角情況下的非定常氣動(dòng)力的建模，需要更加深入地討論馬赫數(shù)、減縮頻率、更高階導(dǎo)數(shù)以及交叉導(dǎo)數(shù)在非定常氣動(dòng)力模型中的作用。