環(huán)量控制翼型非定常氣動力建模

雷玉昌，張登成，張艷華

(1.空軍工程大學(xué) 研究生院，西安 710051；2.中國人民解放軍95034部隊，百色 533600；3.空軍工程大學(xué) 航空工程學(xué)院，西安 710038）

在航空技術(shù)應(yīng)用中，使用環(huán)量控制技術(shù)作為主動流動控制的解決方案越來越引起人們的注意［1-2］。環(huán)量控制技術(shù)是指在翼型后緣表面開縫以形成沿著物面切向的射流，用以增加沿著翼型表面的環(huán)量，進而增加升力［3-4］。風(fēng)洞試驗和數(shù)值仿真結(jié)果表明，環(huán)量控制技術(shù)能夠大幅度提高翼型升力，并在高升力條件下改善升阻比［5-6］。環(huán)量控制技術(shù)最初就是因為其卓越的高升力性能而引起關(guān)注，經(jīng)過八十多年的發(fā)展，該技術(shù)逐漸在無尾飛行器的氣動力控制、風(fēng)力渦輪機的載荷控制、降低噪聲、降低飛行器結(jié)構(gòu)質(zhì)量及提高隱身性能等方面表現(xiàn)出優(yōu)越的性能［7-10］。特別是在無尾飛行器的控制方面，逐漸表現(xiàn)出了替代傳統(tǒng)控制舵面的發(fā)展趨勢。英國BAE系統(tǒng)公司先后試飛了“惡魔(Demon）”、“Magma”等技術(shù)驗證無人機，利用環(huán)量控制技術(shù)進行了滾轉(zhuǎn)操作，驗證了環(huán)量控制射流替代傳統(tǒng)舵面的可行性［11-12］。雖然已經(jīng)證明了環(huán)量控制技術(shù)的相關(guān)優(yōu)越性，但要實現(xiàn)可用的環(huán)量控制射流，還有很多關(guān)鍵理論和技術(shù)問題需要研究。其中考慮迎角改變時環(huán)量控制射流對于非定常氣動力的瞬態(tài)反應(yīng)是必不可少的。因為一方面環(huán)量控制射流并不能像傳統(tǒng)舵面一樣，在迎角改變時與翼型保持同步偏轉(zhuǎn)，而是存在時間遲滯效應(yīng)；另一方面高動量系數(shù)情況下會導(dǎo)致機翼提前失速，導(dǎo)致翼型在中小迎角下就會出現(xiàn)失速分離情況。因此，傳統(tǒng)飛行器的非定常氣動力模型將無法適用于環(huán)量控制飛行器。而飛機過失速機動時會產(chǎn)生流動分離和漩渦破碎，使得氣動力和氣動力矩呈現(xiàn)高度的非線性和非定常特性［13-14］。建立環(huán)量控制翼型在機動過程中的非定常氣動力模型，對于準確理解射流作用機理及進行飛行控制系統(tǒng)設(shè)計具有重要意義。

目前，關(guān)于環(huán)量控制技術(shù)的研究中并未建立起高精度的非定常氣動力模型。Loth和Englar等［15-16］給出了基于定常實驗數(shù)據(jù)的曲線擬合公式，其中，文獻［15］給出的經(jīng)驗公式在較低迎角、單射流情況下與可用實驗數(shù)據(jù)有較好的一致性。Hoholis［17］針對二維環(huán)量控制翼型的氣動參數(shù)進行了曲線擬合，并添加了穩(wěn)定性導(dǎo)數(shù)以凸顯其非定常效應(yīng)，結(jié)果表明，該方法對升力和力矩系數(shù)預(yù)測結(jié)果較好，對阻力系數(shù)預(yù)測較差，而且并不能準確預(yù)測非定常效應(yīng)的氣動力遲滯效應(yīng)。Krukow和Dinkler［18］在進行環(huán)量控制機翼的氣動彈性研究中利用高階多項式擬合得到了定常氣動力模型。Semaan等［19］使用LASSO方法辨識了周期性射流和脈沖型射流疊加后對機翼升力系數(shù)的增量影響。綜上所述，目前對于環(huán)量控制飛行器氣動力模型的建立多集中在定常模式下，對非定常氣動力模型的研究較為匱乏。而傳統(tǒng)飛行器非定常氣動力建模的發(fā)展已經(jīng)歷經(jīng)了幾十年，逐漸發(fā)展形成了氣動導(dǎo)數(shù)模型［20］、狀態(tài)空間模型［21］和微分方程模型［22］等數(shù)學(xué)模型和模糊邏輯模型［23］、神經(jīng)網(wǎng)絡(luò)模型［24］等人工智能模型。這些模型都能在一定程度上描述傳統(tǒng)飛行器在大迎角下的非定常氣動力變化趨勢。

因此，本文借助高精度的數(shù)值仿真技術(shù)，分析環(huán)量控制翼型非定常流動狀態(tài)?；谖⒎址匠棠Ｐ?，添加環(huán)量控制射流參數(shù)源項。通過兩步線性回歸參數(shù)辨識方法建立適用于環(huán)量控制翼型的非定常氣動力數(shù)學(xué)模型，并對高動量系數(shù)下的非線性影響進行修正。通過對任意條件下的非定常流動充分建模，用于環(huán)量控制技術(shù)的非定常狀態(tài)研究。

1 數(shù)值計算方法與驗證

1.1 計算模型和網(wǎng)格

超臨界翼型有著鈍前緣、大厚弦比的特點，內(nèi)部能夠容納環(huán)量控制技術(shù)所需要的供氣機構(gòu)或管道等，因而成為環(huán)量控制技術(shù)良好的研究翼型。本文采用修形過的超臨界翼型進行相關(guān)研究。翼型相對厚度為17%，弦長c＝240 mm。采用文獻［25］中的修形方法，將距后緣25%弦長的部分修形處理，按照后緣半徑r／c＝2%繪制柯恩達后緣曲面，射流口高度h＝0.001c。圖1為修形過的超臨界翼型示意圖。計算區(qū)域選取翼型弦長的30倍，網(wǎng)格剖分時采用結(jié)構(gòu)網(wǎng)格，生成O型網(wǎng)格拓撲結(jié)構(gòu)，邊界層內(nèi)第1層的網(wǎng)格高度約為1×10－5m，射流口附近第1層網(wǎng)格高度為2×10－6m，保證第1層高度的y＋均小于1，以滿足黏性底層的計算要求。網(wǎng)格總數(shù)約為58萬，進行非定常計算時采用FLUENT網(wǎng)格動態(tài)重構(gòu)技術(shù)，圖2為翼型網(wǎng)格劃分及非定常計算過程迎角α從－5°變化到5°翼型后緣網(wǎng)格變形情況?？梢钥闯觯W(wǎng)格適應(yīng)情況較好。

圖1 修形后的超臨界翼型Fig.1 Supercritical airfoil after modification

圖2 翼型計算網(wǎng)格Fig.2 Airfoil’s computational grid

1.2 計算方法和驗證

本文的數(shù)值模擬方法采用二維雷諾平均Navier-Stokes方程，其積分形式為

式中：Q為流動變量；Ω為控制體；FC和FV分別為無黏性通量和黏性通量；n為控制體面的外法線向量。采用k-εSST湍流模型，該模型對于有較大逆壓梯度的邊界層流動、分離預(yù)測性能較好。采用有限體積法離散控制方程，黏性通量采用二階迎風(fēng)格式離散。遠場邊界為壓力遠場，壁面邊界采用無滑移壁面條件。

本文采用的數(shù)值仿真軟件為商用軟件FLUENT。采用的求解器為基于Simple方法的Pressure-Based Solver。對求解不可壓縮流動及微可壓流動效果較好。該方法對于非定常問題的流動計算，主要借助隱式時間積分方案，在每個計算時間步內(nèi)進行迭代，在計算步之內(nèi)的內(nèi)迭代計算與普通穩(wěn)態(tài)問題的內(nèi)迭代計算是相類似的，在得到本時間步內(nèi)的收斂解之后，轉(zhuǎn)入下一個時間步進行繼續(xù)計算。

為了驗證本文所采用數(shù)值計算方法的準確性，分別對標模NACA0012翼型強迫俯仰振動下的非定常氣動力數(shù)據(jù)及環(huán)量控制翼型在射流吹氣下的定常氣動力數(shù)據(jù)進行對比驗證。

對于標模NACA0012翼型而言，計算過程中強迫翼型繞重心做俯仰振蕩，其迎角變化規(guī)律為

式中：減縮頻率k＝ωc／(2V）∞，c為弦長；V∞為來流速度為無量綱時間；α0為基準迎角；αm為振幅。減縮頻率用來描述翼型的振蕩頻率。計算條件：k＝0.081 4，α0＝0.016°，αm＝2.51°。本文計算時采用的物理時間步長t＝0.01T，內(nèi)迭代步數(shù)為50，其中T表示翼型非定常運動周期。為加快收斂，初始解從中心迎角處的定常收斂解開始計算。將計算結(jié)果與文獻［26］中Landon的實驗結(jié)果和Uygun的數(shù)值仿真結(jié)果進行對比，結(jié)果如圖3所示。圖中：CL為升力系數(shù)，α為迎角。在部分迎角下，仿真結(jié)果與實驗數(shù)據(jù)存在差異，但是在總體上與實驗數(shù)據(jù)是一致的，能夠反映非定常流動時的氣動力變化趨勢。

圖3 NACA0012翼型非定常升力系數(shù)仿真與實驗對比Fig.3 Comparison of simulation and experimental results of unsteady lift coefficient of NACA0012 airfoil

對于環(huán)量控制翼型而言，動量系數(shù)是影響環(huán)量控制中后緣射流控制效果的重要參數(shù)，升阻特性受動量系數(shù)影響較大，其定義為

式中：m˙為射流出口處的質(zhì)量流量；Vjet為射流速度；ρ∞為來流密度；V∞為來流速度；S為二維翼型單位展長的面積，在數(shù)值上等于弦長c。假定氣流等熵膨脹至射流出口，遠前方來流靜壓為射流出口靜壓。通過調(diào)節(jié)射流出口速度，調(diào)節(jié)射流出口處的質(zhì)量流量和動量系數(shù)。本文將使用上述參數(shù)作為后緣射流的邊界條件，調(diào)整參數(shù)以進行相關(guān)研究。對0°迎角下不同吹氣動量系數(shù)下的環(huán)量控制翼型進行數(shù)值模擬。計算條件為：來流速度為30 m／s，來流溫度為293.15 K，雷諾數(shù)為5×105，動量系數(shù)在0～0.06之間變化。將仿真結(jié)果與文獻［25］中的實驗數(shù)據(jù)進行對比，仿真采用穩(wěn)態(tài)方法，結(jié)果如圖4所示，仿真結(jié)果與實驗數(shù)據(jù)基本吻合。綜上所述，本文采用的數(shù)值計算方法應(yīng)當(dāng)能夠較好地模擬環(huán)量控制翼型的非定常氣動力數(shù)據(jù)。

圖4 環(huán)量控制翼型升力系數(shù)仿真與實驗對比Fig.4 Comparison of simulation and experimental results of circulation control airfoil’s lift coefficient

2 定常氣動力插值

借助數(shù)值仿真或風(fēng)洞試驗手段得到的翼型以迎角作為自變量的氣動參數(shù)往往具有極大的離散性和不完整性，因此，在進行氣動力和力矩建模時需要對離散的氣動力數(shù)據(jù)進行擬合或插值，使之成為光滑的連續(xù)函數(shù)。本節(jié)基于Kriging模型對定常狀態(tài)下的環(huán)量控制翼型氣動力進行插值處理，并與傳統(tǒng)氣動導(dǎo)數(shù)模型的預(yù)測精度進行對比。

2.1 K riging模型

Kriging模型最初是由南非采礦工程師Krige在1951年提出的一種空間估計技術(shù)，最初應(yīng)用于礦床儲量計算和誤差估計，經(jīng)過幾十年的發(fā)展，逐漸在地質(zhì)、氣象、航空航天、汽車等領(lǐng)域得到發(fā)展和應(yīng)用［27］。Kriging模型由于其卓越的非線性函數(shù)插值預(yù)測能力和誤差估計功能，正在受到越來越多的關(guān)注。

設(shè)m1，m2，…，mn為數(shù)量集上的一系列變量值。z(m1），z(m2），…，z(mn）為相應(yīng)的參數(shù)值。則Kriging模型定義變量 m0處的插值結(jié)果z*(m0）為已知樣本值的線性加權(quán)，即

因此，通過求解加權(quán)系數(shù)λi就可以得到數(shù)量集中任意位置的參數(shù)估計值。關(guān)于Kriging模型的詳細介紹和推導(dǎo)參見文獻［27］。

影響環(huán)量控制翼型定常氣動力的主要因素包括迎角α和動量系數(shù)Cμ。仿真計算了α為－5°～15°，Cμ為0～0.06之間的氣動力參數(shù)。

環(huán)量控制翼型的氣動力參數(shù)呈現(xiàn)高度的非線性特征，迎角α和動量系數(shù)Cμ對氣動力的影響高度耦合，翼型失速與動量系數(shù)大小密切相關(guān)，動量系數(shù)越大，翼型失速迎角越小。以其氣動力參數(shù)可作為兩者的空間分布變量。選取樣本數(shù)量集α為－5°～15°，Cμ為0～0.06之間的離散氣動力參數(shù)。對樣本數(shù)量集進行歸一化后利用Kriging模型分不同方向進行函數(shù)計算，其預(yù)測結(jié)果和均方根誤差結(jié)果如圖5和圖6所示。

圖5 Kriging模型預(yù)測結(jié)果Fig.5 Prediction result of Kriging model

圖6 Kriging模型預(yù)測誤差Fig.6 Prediction error of Kriging model

從整體上來看，Kriging模型對于環(huán)量控制翼型的定常氣動力插值結(jié)果精度較高，特別是對于采樣點范圍內(nèi)的內(nèi)插結(jié)果誤差值低于1×104。對于邊界值誤差較高，這反映了模型對氣動力參數(shù)的外插能力需要進一步提高。Kriging模型能夠?qū)崿F(xiàn)環(huán)量控制翼型氣動力參數(shù)設(shè)計的光滑連續(xù)插值，且具有較高的預(yù)測精度。

2.2 傳統(tǒng)氣動導(dǎo)數(shù)模型

傳統(tǒng)的氣動導(dǎo)數(shù)模型是基于線性疊加原理，假設(shè)氣動力可由飛行器迎角、舵面偏角、角速度等狀態(tài)變量線性疊加，其本質(zhì)是氣動力的泰勒展開式分解后只保留其一階表達式，保留的各狀態(tài)變量系數(shù)分別被稱為氣動靜、動導(dǎo)數(shù)。對于定常氣動力而言，僅涉及氣動靜導(dǎo)數(shù)。為了提高氣動導(dǎo)數(shù)模型在大迎角非線性情況下的預(yù)測精度，一般采用高階展開式進行擬合，以升力系數(shù)為例可以表示為

將射流的影響分為兩部分：一部分設(shè)置為初始狀態(tài)量，包含在系數(shù)f0中，另一部分用于修正迎角帶來的氣動力影響，包含在系數(shù)f1和f2中。將系數(shù)fi設(shè)置為射流動量系數(shù)的三次多項式函數(shù)：

分別在不同的吹氣動量系數(shù)下驗證Kriging模型與傳統(tǒng)氣動導(dǎo)數(shù)模型的準確性，將2種模型的預(yù)測結(jié)果與CFD數(shù)值計算值進行比較，圖7給出了不同動量系數(shù)下2種模型的預(yù)測情況。

圖7 Kriging模型與傳統(tǒng)氣動導(dǎo)數(shù)模型對比結(jié)果Fig.7 Comparison between Kriging model and aerodynamic derivative model

從圖7中可以看出，在采樣范圍的動量系數(shù)(Cμ為0～0.06）作用下，Kriging模型較傳統(tǒng)的氣動導(dǎo)數(shù)模型預(yù)測精度有明顯提升，預(yù)測結(jié)果更加接近CFD計算值，并能準確預(yù)測失速迎角范圍；對于采樣范圍之外的參數(shù)，Kriging模型仍然能保持一定的預(yù)測精度，預(yù)測效果較傳統(tǒng)氣動導(dǎo)數(shù)模型有較大提升。

3 非定常氣動力建模

環(huán)量控制射流相較于傳統(tǒng)舵面而言可以大幅度提高翼型升力、改善升阻比，這使得環(huán)量控制飛行器往往在小迎角下就能夠獲得傳統(tǒng)飛行器在大迎角下的氣動力特性。因此，本文主要進行中低迎角(－5°～10°）下的環(huán)量控制翼型非定常氣動力建模。

3.1 環(huán)量控制翼型非定常氣動力模型

在機翼小振幅強迫振蕩風(fēng)洞試驗中，若還未出現(xiàn)深動態(tài)失速現(xiàn)象，則根據(jù)小擾動假設(shè)，氣動力參數(shù)可以由高階線性狀態(tài)方程表示，包括定常和非定常兩部分［28］。其中定常氣動力反映了機翼振蕩過程中的氣動力快速響應(yīng)部分，非定常氣動力部分則反映了機翼運動導(dǎo)致的流場內(nèi)渦系移動產(chǎn)生的氣動力遲滯效應(yīng)。已有的研究表明［29］，環(huán)量控制射流通過改變翼型表面壓力分布、控制分離點移動而改變氣動力系數(shù)，這與迎角對氣動力系數(shù)的作用機理一致。對于環(huán)量控制翼型而言，動量系數(shù)與迎角共同影響翼型的氣動力系數(shù)，并且呈現(xiàn)高度耦合。以縱向運動為例，氣動力參數(shù)可以表示為

式中：Ci為氣動力力矩系數(shù)；Catt(α，Cμ）為氣流在附著流動的狀態(tài)下，即無分離假設(shè)下的氣動力系數(shù)；Cα˙att(α，Cμ）為機翼定常俯仰下的氣動力導(dǎo)數(shù)，為迎角變化率；Cdyn為非定常氣動力增量。

上述模型適用于尚未出現(xiàn)深動態(tài)失速時的氣動力系數(shù)變化情況，其中，Cdyn采用以下線性微分方程式表示：

式中：τ為非定常運動中的特征時間常數(shù)；Cst(α，Cμ）為實際定常狀態(tài)下氣動力系數(shù)，即Kriging模型預(yù)測值；ΔC(α，Cμ）為由于氣流分離導(dǎo)致的定常氣動力增量。

對于環(huán)量控制翼型而言，假設(shè)在迎角俯仰過程中射流動量系數(shù)不變，這種假設(shè)是合理的，目前對于傳統(tǒng)舵面飛行器在大迎角下非定常氣動力模型的建立也是基于舵面偏角不變的情況。

3.2 非定常氣動力模型參數(shù)辨識

借助經(jīng)典的強迫振動法計算諧振運動過程中的非定常氣動力，強迫翼型繞重心做小幅度的俯仰振蕩，其迎角變化規(guī)律為

式中：Δα為迎角改變幅度。

非定常氣動力系數(shù)可以看做時間ˉt的連續(xù)函數(shù)，因此，可以將翼型小幅度俯仰振蕩過程中的氣動力系數(shù)變化按照減縮頻率進行傅里葉分解，以升力系數(shù)為例，保留一階形式可以表示為

式(10）表示了氣動力系數(shù)可類似于“小擾動”假設(shè)，由基準系數(shù)加上擾動系數(shù)構(gòu)成。其中，C0(α0，Cμ0）表示中心迎角對應(yīng)下的基準氣動力系數(shù)；Cα(α0，Cμ0，k）稱為同相導(dǎo)數(shù)；Cα˙(α0，Cμ0，k）稱為異相導(dǎo)數(shù)。

基于小振幅情況下的線性假設(shè)：

該假設(shè)表示在小振幅情況下實際靜態(tài)氣動力系數(shù)與氣流不分離假設(shè)下靜態(tài)氣動力系數(shù)之量與迎角變化近似呈現(xiàn)線性關(guān)系。ΔCα(α0，Cμ0）表示該差值與迎角變化之間的導(dǎo)數(shù)。

借助歐拉方程求解微分方程式(7），即求解下述方程組：

使用微分算子法容易解得方程組(12）的解，將式(8）、式(9）的迎角變化規(guī)律代入式(6）中，與式(10）的相關(guān)系數(shù)作對比后整理得到以下關(guān)系：

由上述推導(dǎo)可知，同相導(dǎo)數(shù)并不等同于靜導(dǎo)數(shù)，異相導(dǎo)數(shù)也并不等同于動導(dǎo)數(shù)。

但是當(dāng)τ非常小時，此時同相導(dǎo)數(shù)和異相導(dǎo)數(shù)分別退化為傳統(tǒng)線性氣動導(dǎo)數(shù)模型的靜導(dǎo)數(shù)Cαst(α0，Cμ0）和動導(dǎo)數(shù)Cα˙att(α0，Cμ0），即如下關(guān)系式：

式中：Cα(α0，Cμ0）＝Cαst(α0，Cμ0）。

此時，Cdyn(0）＝0，C0(α0，Cμ0）＝Cst(α0，Cμ0）。

τ定量地表征了翼型非定常氣動力效應(yīng)的顯著程度。解式(13）可得

由于在同一中心迎角、不同減縮頻率對應(yīng)下的同相／異相導(dǎo)數(shù)數(shù)據(jù)點落在同一條直線上，而直線的斜率則表示特征時間常數(shù)的負值：－τ。

對于某一特定中心迎角下的小振幅強迫振動數(shù)據(jù)而言，采用最小二乘法可得到對應(yīng)中心迎角下的同相／異相導(dǎo)數(shù)數(shù)據(jù)點，部分數(shù)據(jù)點如圖8所示。

圖8 不同迎角、減縮頻率下的同相／異相導(dǎo)數(shù)Fig.8 In-phase and out-phase derivatives at different angles of attack and reduced frequencies

由上述分析可知，對于特定的小振幅流動狀態(tài)(特定的迎角和動量系數(shù)）下，τ是固定的，計算得到的不同中心迎角、動量系數(shù)下的特征時間常數(shù)如圖9所示。

圖9 時間常數(shù)隨動量系數(shù)的變化曲線Fig.9 Variation curve of time constantwith momentum coefficient

特征時間常數(shù)值受迎角和動量系數(shù)影響較大?？傮w來看，τ隨Cμ增大呈現(xiàn)先降低后增加的趨勢，這說明低動量系數(shù)下，翼型的非定常氣動效應(yīng)減弱，但是當(dāng)動量系數(shù)增加到一定值后，翼型的非定常氣動力效應(yīng)增強，并且隨著迎角的增大，該臨界動量系數(shù)逐漸降低。當(dāng)α0＝0°時，在全動量系數(shù)范圍內(nèi)時間常數(shù)趨于0，氣動力模型可近似用傳統(tǒng)線性氣動導(dǎo)數(shù)模型代替。當(dāng)α0＝10°時，在全動量系數(shù)范圍內(nèi)時間常數(shù)都較大，翼型處于臨近動態(tài)失速或失速后狀態(tài)，此時翼型的非定常氣動力特性主要由后緣渦系的遲滯運動導(dǎo)致，傳統(tǒng)線性氣動導(dǎo)數(shù)模型不再適用。

對于τ較小的情況，基于傳統(tǒng)線性氣動導(dǎo)數(shù)模型，使用不同迎角和動量系數(shù)下非定常運動的若干離散點，借助最小二乘法辨識靜導(dǎo)數(shù)Cαst(α0，Cμ0）和動導(dǎo)數(shù)Cα˙att(α0，Cμ0）。

對于τ較大的情況，由于無法使用傳統(tǒng)的線性氣動導(dǎo)數(shù)模型，在獲得τ(α，Cμ）后，仍需要辨識同相／異相導(dǎo)數(shù)中的Cαatt(α0，Cμ0）、ΔCα(α0，Cμ0）和Cα˙att(α0，Cμ0）。

為了便于描述，將式(13）用以下方程描述：

式中：y1i、y2i、b0、a0、c0分別為式(13）中的相關(guān)參數(shù)。為了辨識相關(guān)參數(shù)，引入以下誤差函數(shù)：

式中：n為同一中心迎角和動量系數(shù)、不同頻率下的仿真實驗次數(shù)；y1i和y2i分別為特定小振幅流動狀態(tài)下的同相／異相導(dǎo)數(shù)；b0、a0和c0為參數(shù)辨識后的相應(yīng)值。該方法的核心思想是兩步線性回歸法，即先根據(jù)非定常運動的離散點辨識得到y(tǒng)1i和y2i，再根據(jù)得到的y1i和y2i辨識b0、a0和c0。辨識過程中，因為存在中心迎角、動量系數(shù)、減縮頻率3個自變量，直接進行參數(shù)辨識需要耗費極大的仿真計算量，所以本節(jié)針對射流動量系數(shù)與迎角的相互作用特點，對b0、a0和c0做合理優(yōu)化，提高計算效率。

對于b0而言，Jones［25］曾指出動量系數(shù)對氣動力的影響呈現(xiàn)2個比較明顯的線性階段，其流動狀態(tài)滿足無分離假設(shè)，因此在低動量系數(shù)下b0可以用以下關(guān)系式代替：

式中：CαCμst(α0）為Cαst(α0）對Cμ的導(dǎo)數(shù)，可直接由Kriging模型插值得到，待辨識參量變?yōu)橛?。采用這樣的簡化方法在低迎角線性階段可直接退化為Cαst(α0，Cμ0），與無分離情況下的氣動導(dǎo)數(shù)一致。

對于c0而言，當(dāng)自由來流不發(fā)生變化時，迎角俯仰導(dǎo)致翼型表面產(chǎn)生繞重心的旋轉(zhuǎn)速度，與外流相對流動產(chǎn)生下洗效應(yīng)，而后緣射流并不會直接隨著翼型俯仰，而是由于后緣壓力場的變化時刻處于動態(tài)平衡之中。由定常旋轉(zhuǎn)產(chǎn)生的準定常氣動力增量主要在小迎角時起作用，與迎角變化率近似呈線性關(guān)系。在較大迎角范圍內(nèi)，其引起的準定常氣動力增量遠小于非定常氣動力增量。將動量系數(shù)作為修正參數(shù)，可以將c0表示為

對于a0而言，a0＝Cst(α，Cμ）－Catt(α，Cμ），因此只需辨識b0即可。在辨識得到τ和ΔC(α，Cμ）后，根據(jù)非定常氣動力增值周期性變化的特性，即Cdyn(0）＝Cdyn(2π／ω），即可求得初值Cdyn(0）。

以誤差函數(shù)E最小為優(yōu)化目標函數(shù)，以升力系數(shù)為例，得到的Cαatt(α0，Cμ0）曲線如圖10所示。

圖10 Cαatt(α0，Cμ0）的辨識結(jié)果Fig.10 Identification results of Cαatt(α0，Cμ0）

3.3 非定常氣動力模型與仿真數(shù)據(jù)對比驗證

為了驗證上述非定常氣動力模型的準確性，分別對環(huán)量控制翼型在不同迎角、動量系數(shù)、振幅下的預(yù)測結(jié)果進行對比分析。預(yù)測結(jié)果取升力系數(shù)CL和俯仰力矩系數(shù)Cm進行對比分析。計算條件為：來流速度為30 m／s，來流溫度為293.15 K。中心迎角、動量系數(shù)、振幅、減縮頻率驗證如圖11、圖12所示。

如圖11所示，為了得到非定常氣動力模型對大振幅情況下的預(yù)測結(jié)果，將大振幅的流動情況分段處理，每一段的氣動力變化情況都可以用小振幅情況下迎角對應(yīng)的同相／異相導(dǎo)數(shù)表示，分段處理后可以將大振幅情況下的預(yù)測結(jié)果近似表示為

圖11 小振幅流動狀態(tài)下結(jié)果對比Fig.11 Comparison of results under small-amplitude flow

將分段處理的αm取為1°，得到非定常氣動力模型對于大振幅情況下的預(yù)測結(jié)果，如圖12所示。

圖12 大振幅流動狀態(tài)下結(jié)果對比Fig.12 Comparison of results under large-amp litude flow

可以看出，對于升力系數(shù)而言，線性非定常氣動力模型對小振幅情況下流動狀態(tài)有較好的預(yù)測精度，特別是在0°迎角下，此時升力系數(shù)線性特征明顯，在較大振幅情況下(見圖11(a）、圖11(b））也能夠取得較好的預(yù)測精度。但是在較大迎角、較高動量系數(shù)的流動狀態(tài)下，如圖12(c）所示，此時升力系數(shù)非線性特征明顯，線性非定常氣動力模型難以取得較高的預(yù)測精度。對于俯仰力矩系數(shù)而言，線性非定常氣動力模型在各個動量系數(shù)下都能較好地預(yù)測力矩系數(shù)變化趨勢。這是因為俯仰力矩系數(shù)主要受翼型上表面分離渦影響，當(dāng)翼型上表面分離渦移動時，對應(yīng)的負壓區(qū)逐漸擴大或者縮小，導(dǎo)致俯仰力矩系數(shù)變化。高動量系數(shù)下翼型上表面不再分離，導(dǎo)致遲滯效應(yīng)逐漸減弱。值得一提的是，在出現(xiàn)動態(tài)失速分離前，高動量系數(shù)具有減弱氣動力遲滯效應(yīng)的優(yōu)勢，如圖11(f）、圖12(f）所示，此時俯仰力矩的變化幾乎不存在遲滯效應(yīng)，與定常狀態(tài)下的變化幾乎一致。在出現(xiàn)動態(tài)失速后，同樣能在線性小迎角階段減弱遲滯效應(yīng)，如圖12(b）、圖12(c）所示。

3.4 非定常氣動力模型非線性修正

當(dāng)振蕩幅度增大并出現(xiàn)深動態(tài)失速的情況時，上述線性非定常氣動力模型將會在失速段失去預(yù)測精度。這是因為在大振幅情況下進行分段處理完全忽略了洗流時差的累計效應(yīng)，在非線性較為明顯的階段將完全喪失其預(yù)測精度。為了改善非定常氣動力模型在動態(tài)失速情況下的預(yù)測精度，采用文獻［30］的方法將非定常氣動力增量式(7）改寫為以下非線性形式：

式中：ki(α，Cμ）＝1／τi(α，Cμ），該方法采用多次逼近的思想，避免了大振幅過程中僅用單一時間常數(shù)τ線性表示增量的誤差，同時充分考慮了累計效應(yīng)。利用上述線性非定常氣動力模型辨識得到的Catt(α，Cμ）和Cα˙att(α，Cμ），可以獲得線性假設(shè)下CL的變化，與實際CL的差值即非定常氣動力增量，與線性假設(shè)下得到的Cdyn相疊加即實際情況下的Cdyn，保持ΔC(α，Cμ）和k1(α，Cμ）不變，利用最小二乘法即可辨識得到ki(α，Cμ），采用這樣的假設(shè)能夠有效降低參數(shù)辨識的難度。

上述線性非定常氣動力模型已經(jīng)能夠滿足無射流情況下中低迎角的大振幅預(yù)測精度，因此，非線性項主要用來提高高動量系數(shù)情況下的預(yù)測精度。為了提高辨識效率，本文僅對包含Kriging模型中靜態(tài)失速點附近的流動狀態(tài)進行參數(shù)辨識。采用三次非線性形式描述非定常氣動力增量，辨識得到的k2(α，Cμ）和k3(α，Cμ）如圖13所示。

圖13 ki(α，Cμ）辨識結(jié)果Fig.13 Identification result of ki(α，Cμ）

非線性項修正后的非定常氣動力模型進行大振幅流動情況下的升力系數(shù)預(yù)測，結(jié)果如圖14所示，其結(jié)果較線性模型有明顯提升。

圖14 非線性非定常氣動力模型結(jié)果對比(Cμ ＝0.05）Fig.14 Comparison of results of nonlinear unsteady aerodynamic models(Cμ ＝0.05）

4 結(jié)論

本文主要進行了環(huán)量控制翼型非定常氣動力的建模，通過上述分析，可以得到以下結(jié)論：

1）基于Kriging模型完成了環(huán)量控制翼型的定常氣動力插值，能夠精確的反映其定常氣動力變化，插值精度要優(yōu)于傳統(tǒng)氣動導(dǎo)數(shù)模型。

2）基于微分方程模型，通過添加射流參數(shù)源項，完成了適用于環(huán)量控制翼型的線性微分方程建模。對Cαatt(α0，Cμ0）和Cα˙att(α0，Cμ0）等參數(shù)簡化處理，提高辨識效率。

3）對高動量系數(shù)情況下的非定常氣動力模型進行非線性修正，能夠較為準確的反映不同動量系數(shù)下的非定常流動特性。

4）對于升力系數(shù)而言，線性非定常氣動力模型能夠滿足小振幅和低動量系數(shù)下大振幅流動狀態(tài)；對于俯仰力矩系數(shù)而言，線性非定常氣動力模型能夠滿足大部分流動狀態(tài)。

環(huán)量控制翼型的非定常氣動力耦合特性復(fù)雜，仍需要在以下方面進一步研究：①高動量系數(shù)下缺乏足夠精度的數(shù)學(xué)模型來描述后緣渦的分離變化，需要對流場特性進行更深入的研究，進而優(yōu)化模型結(jié)構(gòu)。提高在高動量系數(shù)下非線性非定常氣動力模型預(yù)測效果。②關(guān)于環(huán)量控制翼型的非定常氣動力建模研究仍處于起步階段，如何有效的結(jié)合已經(jīng)較為成熟的大迎角非定常氣動力數(shù)學(xué)模型是未來研究的途徑之一。