三角翼俯仰振蕩的非定常氣動力降階方法研究*

郭嘉瑞，關世璽，李立州，賈 凱，常 晶

(中北大學機電工程學院 太原 030051)

0 引言

大后掠三角翼能夠減小現代飛行器在高速飛行時的氣動阻力，并改善大迎角飛行的操縱性能。隨著高性能計算設備與高精度計算流體力學(CFD)方法的發展，國內外對三角翼在運動狀態下的非定常氣動力開展了大量數值模擬研究。Gordnier等[1]數值模擬了80°三角翼的滾轉運動，分析了動態流場中渦運動對非定常氣動特性的影響。Ekaterinar等[2-3]數值模擬了雙三角翼俯仰運動，計算得到的靜止狀態下物面壓力分布和俯仰運動的非定常時滯力矩均與實驗結果大致吻合。郭迪龍等[4]研究了俯仰滾轉耦合作用對三角翼非定常氣動力的影響。劉昕等[5]數值模擬了雙三角翼的俯仰運動，研究了減縮頻率、平均迎角和振幅對非定常氣動力的影響。

三角翼大迎角下非定常流場的高精度數值模擬需要花費較多的計算資源和計算周期。因此，需要研究一種高效、高精度的降階模型來代替CFD方法預測三角翼的非定常氣動特性。Suei Chin[6]等是用Fourier變換的數學模型來預測70°三角翼大振幅俯仰運動的非定常氣動力。Huang等[7]采用NIR-ISS方法，建立了大攻角下三角翼大振幅滾轉運動的非定常氣動力模型。史志偉等[8]采用非線性代數模型，建立了大振幅滾轉運動的數學模型并提取了動導數。孫海生等[9]開展了大迎角、大振幅運動的非定常氣動力建模方法的比較研究。上述建模研究成果多集中于大幅度運動時氣動力宏觀變化規律的預測，未涉及對小幅振蕩氣動力的線性與非線性分量的精確重構。

文中采用小波多分辨率分析方法對多輸入Volterra核進行壓縮，并預測三角翼俯仰振蕩的非定常氣動力。通過求解非定常Euler方程，得到了偽隨機信號激勵的三角翼氣動力與力矩響應，辨識了76°尖前緣三角翼在兩種典型迎角附近的俯仰振蕩氣動力與力矩。討論了一階、二階Volterra核能夠辨識氣動響應的線性與非線性分量的能力。

1 數值模擬方法

流動控制方程采用三維非定常Euler方程，在直角坐標下的守恒積分形式為：

(1)

式中：Q為單位體積內質量、動量和能量組成的守恒變量；F為對流矢通量；S為邊界外法向面積向量。利用有限體積法構造空間半離散格式，對流項采用二階精度的Roe格式，時間推進采用二階隱式LUSGS格式。動網格采用剛性旋轉法，嚴格遵守幾何守恒律[10]的要求。

2 非定常氣動力降階方法

2.1 單輸入Volterra級數

對于一個多輸入非線性時不變系統，對時域離散的Volterra級數可表示為如下形式：

(2)

式中：n為離散時間變量；h是系統的Volterra核；M1、M2為記憶長度。

傳統的非定常氣動力特點是一階Volterra級數[11]，其向量形式為：

y(n)=h0+uTh1

(3)

為增強模型辨識非線性系統的能力，將Volterra級數擴展至二階形式：

y(n)=h0+uTh1+uTh2u

(4)

2.2 Volterra核的小波壓縮

對Volterra一階、二階核采用離散小波變換，得到小波域下表達形式：

一階核：

(5)

二階核：

(6)

式中

分別為一維與二維小波變換矩陣；W1、W2分別為一維與二維小波重構矩陣，由Mallat算法得到[12-13]：

(7)

(8)

一階模型：

(9)

二階模型：

(10)

式中：

2.3 自適應QR-RLS算法

采用自適應QR-RLS算法[14]辨識小波域下的Volterra核。RLS算法用卡爾曼濾波的遞推公式推導出自適應濾波器權重適量更新方程，該算法具有比LMS算法更快的收斂速度，且對特征值擴展度不靈敏。其缺點是時間平均自相關矩陣的逆矩陣由于累計數值誤差而失去非負定性，算法將迅速發散。為了改善RLS算法的穩定性，引入QR分解算法直接對輸入矩陣進行三角分解，有效降低了自相關矩陣的條件數，提高了RLS算法的穩定性。QR-RLS算法具體步驟為：

1)k=0時，初始化矩陣U(0)：

(11)

式中δ為一個小量。

2) 對每一個k=1,2,…有：

(12)

矩陣Qθ為酉矩陣，通過連續Givens變換得到。矩陣Qθ將新的輸入信號向量旋轉到主對角線的上三角矩陣上。

3) 對期望輸出進行連續Givens變換，得到先驗誤差eq1:

(13)

4) 向后迭代求解濾波器系數:

U(k)η(k)=dq2(k)

(14)

5) 計算后驗誤差:

(15)

式中：γ為旋轉因子，λ為遺忘因子，是一個不大于1的數。

3 計算結果與分析

3.1 模型、網格與計算驗證

模型為一尖前緣、76°后掠三角翼，這與文獻[15]中描述的風洞模型一致，根弦長作歸一化處理。網格拓補結構為O-H型，適當加密了前緣及物面法向。沿流向、周向、法向的網格點分布為81×129×44，且第一層網格距物面1×10-4倍弦長。圖1給出了網格空間結構。

通過三角翼大振幅俯仰運動風洞實驗驗證了數值方法的非定常計算能力[15]。來流Ma數為0.3，迎角變化為α=0.024t，旋轉中心為2/3根弦長處。以零迎角時的定常流場為初場，時間步長取(1×10-3) s，模擬迎角從0°～80°的俯仰運動。圖2(a)、圖2(b)比較了非定常升力、阻力系數的計算值和實驗值。可以看出，計算程序對三角翼俯仰運動具有模擬的能力。當迎角大于40°時，計算值與實驗值之間具有一定差異。這是由于迎角過大時三角翼前緣渦發生破裂，而Euler方程對渦破裂的模擬有所不足，但不失研究氣動力降階方法的可行性。

3.2 非定常氣動力降階

用上述小波壓縮的Volterra核，在25°和35°迎角附近，辨識出三角翼的俯仰振蕩非定常氣動力。最優多層偽隨機信號[16]作為輸入的迎角變化信號。為提高辨識效率，對輸入信號進行濾波，并使其帶寬覆蓋俯仰運動的減縮頻率。

為驗證辨識結果的準確性，通過降階模型分別計算3°與5°幅值的俯仰振蕩運動，并與CFD結果進行對比，比對的運動狀態為：

狀態1:α=25°+3°×sin(5t)

(16)

狀態2:α=25°+5°×sin(5t)

(17)

1)25°迎角情形

從圖2(a)可以看出，在25°迎角附近，三角翼俯仰運動的非定常升力系數保持較好的線性特性，因此用一階Volterra級數辨識非定常氣動力。Volterra核記憶長度為32，小波壓縮后的實際辨識個數為24。圖3給出了輸入信號激勵下的升力系數與力矩系數的CFD計算結果。圖4給出了一階Volterra核辨識結果，可以看出，Volterra核衰減得較快，而振蕩幅值較弱。

圖5、圖6分別比較了CFD和降階模型得到的非定常升力系數和俯仰力矩系數的時域和頻域結果。可以看出，兩種運動狀態的非線性氣動特性相對較弱。當振幅為3°時，降階模型和CFD的計算結果在時域和頻域中基本一致。當振幅增加到5°時，氣動響應的非線性主頻的能量略有增加。一階Volterra級數準確地捕捉了非定常氣動響應的線性主頻。在峰值處兩者的時域結果略有不同，是因為降階模型無法表征非定常氣動響應的非線性分量。

2)35°迎角情形

從圖2(a)所示的升力系數曲線可以看出，在35°迎角附近，非定常升阻力系數開始出現明顯的非線性特性。非線性氣動響應分別用一階、二階Voltterra級數辨識。一階核記憶長度為128，小波壓縮后的實際辨識個數為96。二階核記憶長度為32，對應的辨識參數個數為1 024，小波壓縮后的實際辨識個數為576。圖7給出了在輸入信號激勵下的升力系數與力矩系數的CFD計算結果。與25°迎角相比，升力和俯仰力矩響應的平均值更高且振幅相對較小。這表明隨著迎角的增大，前緣渦吸力增加，但所產生的非線性特性抑制了氣動響應的振幅。圖8給出了一階Volterra核辨識結果，可見與25°迎角相比，一階核的收斂比較緩慢，且振蕩幅值比較明顯。圖9,圖10分別給出了升力系數和俯仰力矩系數的二階Volterra核的辨識結果。

圖11、圖12比較了兩種運動狀態下，一階、二階Volterra級數與CFD計算得到的氣動力時域和頻域響應。可以看出，與25°迎角相比，氣動響應的非線性主頻能量均大幅增加。由于一階Volterra核無法辨識出非線性分量的頻率和幅值，因此一階模型所得結果與CFD計算結果有較大差異。而二階Volterra核能較為準確地捕捉非線性分量，辨識效果明顯優于一階模型。

4 結論

Euler方程用于數值模擬76°尖前緣三角翼的俯仰運動。利用小波壓縮Volterra級數，辨識三角翼正弦振蕩的非定常氣動力，得出如下結論：

1) 隨著平均迎角的增大，三角翼前緣渦能量和渦吸力的非線性效應顯著增加。渦能量的增加增大了氣動響應的均值，而渦吸力的非線性效應抑制了氣動響應的振幅。

2) 一階Volterra核具有精確的能力來辨識氣動響應線性分量的頻率和幅值。但無法辨識出非線性分量，對小迎角小振幅下的氣動響應辨識效果較好。二階Volterra核能較為準確的捕捉到非定常氣動力和力矩的二次非線性分量，明顯提高了在較大迎角和振幅下，對三角翼俯仰振蕩非定常氣動響應的預測能力。

3) 小波變換有助于減少一階和二階Volterra核的實際識別參數的個數，提高識別效率。