柔性后緣可變形機翼氣動彈性分析方法研究

趙仕偉，闞梓，李道春

北京航空航天大學，北京 100191

可變形機翼可以根據(jù)飛行環(huán)境和飛行任務的變化光滑、連續(xù)地改變自身形狀，使飛機在整個任務周期內(nèi)的總體性能達到最優(yōu)。然而其光滑連續(xù)變形的氣動外形特性導致非定常氣動力變化劇烈，需采用高精度氣動力分析方法。傳統(tǒng)計算流體力學（CFD）氣動分析方法[1]相比面元法[2]和降階方法計算精度高，但計算量大、耗時多，不便于開展氣動彈性分析。氣動力模型降階可建立一種精度高、適應性強的可變形機翼非定常氣動力模型，用于可變形機翼氣動彈性分析計算。

目前，常見的有兩類降階非定常空氣動力學建模方法：流場特征分析方法和系統(tǒng)辨識方法。流場特征分析方法針對整個流場進行模型降階分析，得到流場的特征狀態(tài)參數(shù)；包括正交分解（POD）方法[3-6]和諧波平衡（HB）方法[7-9]。基于流體計算模型的輸入響應獲得正交基，模型的狀態(tài)量可視為該組正交基的線性疊加，可以通過增加模型階數(shù)從而更準確地模擬實際流場運動特性。但該方法通常用于流場的分析，由于模型考慮整個流場狀態(tài)特性，模型較為復雜、階數(shù)較高，不便于開展氣動彈性分析。

系統(tǒng)辨識方法對流場的某些特定參數(shù)（如升力系數(shù)、力矩系數(shù)等）進行系統(tǒng)特征識別。非定常氣動力系統(tǒng)辨識常見的有Volterra級數(shù)方法[10-11]、線性狀態(tài)空間模型[12-13]、自回歸移動平均（ARMA）模型[14]、基于人工神經(jīng)網(wǎng)絡的替代模型[15-16]、徑向基函數(shù)（RBF）插值和Kriging模型[17-20]、支持向量機（SVM）方法[21-22]、狀態(tài)觀測器[23]、面向塊的Wiener 模型[24-25]等。

采用非定常氣動力降階模型，可以提高柔性后緣可變形機翼的氣動彈性響應計算效率，常見的模型降階方法有帶外輸入的自回歸模型（ARX）方法和Volterra 級數(shù)方法。ARX方法的非定常氣動力模型階數(shù)相比Volterra級數(shù)方法要低一個量級左右，計算效率比Volterra 級數(shù)高。因此，本文采用ARX方法建立三維非定常氣動力，與結(jié)構(gòu)狀態(tài)空間模型相耦合，得到柔性后緣可變形機翼氣動彈性系統(tǒng)的狀態(tài)空間模型，在此基礎上計算柔性后緣可變形機翼的顫振特性和陣風響應。

1 氣動彈性狀態(tài)空間模型

非定常氣動力ARX模型可以寫成如下表達式

式中，y(k)表示第k步的廣義氣動力系數(shù)，na、nb是系統(tǒng)廣義氣動力系數(shù)輸入和廣義位移輸入階數(shù)，可以將該系統(tǒng)改寫成離散狀態(tài)空間的形式

式中

狀態(tài)矢量如下

對于多自由度系統(tǒng)，通過歸一化可以轉(zhuǎn)換為模態(tài)坐標運動，結(jié)構(gòu)模態(tài)坐標下的動力學方程可表示如下

其中，M為系統(tǒng)質(zhì)量矩陣，歸一化處理后是單位矩陣；C為系統(tǒng)阻尼矩陣；K為該多自由度結(jié)構(gòu)系統(tǒng)剛度矩陣；x，x?，x?分別為該系統(tǒng)位移、速度、加速度。各階模態(tài)對應的廣義氣動力FA為動壓q、壓力系數(shù)Cp和結(jié)構(gòu)模態(tài)分析得到的各階振型矢量Si的乘積在各個氣動表面元上的積分

該系數(shù)可以Fluent 數(shù)值仿真獲取。通過編寫UDF 程序，在每一個時間步輸出對應的各階模態(tài)對應的廣義氣動力，從而得到各階模態(tài)輸入下的各階模態(tài)對應的廣義氣動力。

將上述模態(tài)坐標下的動力學方程寫成狀態(tài)空間表達形式，可以得到

式中

狀態(tài)矢量xs(t)為

將上述連續(xù)時間系統(tǒng)狀態(tài)空間結(jié)構(gòu)模型進行離散化處理，得到對應的離散狀態(tài)空間模型如下

其中：

式中，T為時間步長，上述過程為系統(tǒng)的結(jié)構(gòu)動力學模型的獲取流程，廣義氣動力為系統(tǒng)輸入，結(jié)構(gòu)彈性變形為系統(tǒng)輸出。

通過ARX 方法建立離散形式的氣動力狀態(tài)空間模型如下

式中，u(k)為系統(tǒng)的廣義結(jié)構(gòu)位移；yA(k)為系統(tǒng)的廣義氣動力。

將基于ARX 方法建立的氣動力狀態(tài)空間模型和結(jié)構(gòu)狀態(tài)空間模型相耦合，可以得到整個氣動彈性系統(tǒng)的狀態(tài)空間模型，可以用于分析系統(tǒng)的氣動彈性特性，計算系統(tǒng)的氣動彈性響應。氣動彈性系統(tǒng)狀態(tài)空間模型表達式如下

上述狀態(tài)空間模型可以通過氣動彈性響應計算或者矩陣特征值來分析該氣動彈性系統(tǒng)的穩(wěn)定性。通過調(diào)整速度和調(diào)整動壓，從而獲得對應速度下的系統(tǒng)響應，當響應幅值隨著時間變化保持不變，該速度對應該馬赫數(shù)下的顫振速度。也可以求解狀態(tài)空間的特征值來分析該系統(tǒng)的穩(wěn)定性。對于離散系統(tǒng)，當所有特征值模長均小于等于1時，該系統(tǒng)穩(wěn)定；若至少有一個特征值模長大于1，該系統(tǒng)不穩(wěn)定?？梢詫㈦x散狀態(tài)空間模型轉(zhuǎn)化為連續(xù)狀態(tài)空間模型來分析系統(tǒng)穩(wěn)定性，當且僅當所有系統(tǒng)的特征值實部小于0，該系統(tǒng)穩(wěn)定。

通過ARX 方法同樣可以建立離散形式的陣風氣動力狀態(tài)空間模型，如下所示

式中，xg(k)為陣風輸入；u(k)是系統(tǒng)的廣義結(jié)構(gòu)位移；yg(k)是系統(tǒng)的陣風廣義氣動力。

將基于ARX 方法建立的氣動彈性狀態(tài)空間模型和陣風響應氣動力狀態(tài)空間模型相耦合，可以得到考慮陣風響應氣動彈性系統(tǒng)的狀態(tài)空間模型，用于分析系統(tǒng)的陣風響應。氣動彈性系統(tǒng)狀態(tài)空間表達式如下

2 氣動彈性特性分析方法驗證

為驗證基于ARX 方法建立的氣動力降階模型的有效性，本文采用Agard445.6 機翼模型。Agard445.6 幾何形狀如圖1 所示，機翼根部弦長0.559m，半展長0.762m，展弦比為1.65，梢根比為0.66，后掠角為45°，翼型為NACA 65A004。設置結(jié)構(gòu)網(wǎng)格最大尺寸為10mm，對機翼根部所在的平面使用固支約束，對該結(jié)構(gòu)模型進行模態(tài)分析。圖2為Agard445.6 機翼對應的前四階模態(tài)振型，前四階模態(tài)的振動頻率與實驗結(jié)果對比見表1，最大誤差為3.15%，有限元與試驗結(jié)果[26]吻合較好。

表1 前四階模態(tài)的振動頻率與試驗結(jié)果對比Table 1 Comparison between the vibration frequency of the first four modes and the experimental results

圖1 Agard445.6幾何形狀Fig.1 Geometry of Agard445.6

圖2 Agard445.6機翼前四階模態(tài)振型圖Fig.2 The first four modal shapes of Agard445.6 wing

圖3 所示為Agard445.6 機翼非結(jié)構(gòu)網(wǎng)格并導入Fluent中計算。遠場邊界位于距機翼20倍翼根弦長處，并設置為壓力遠場條件。機翼表面設置為無滑動靜態(tài)壁面條件，網(wǎng)格在翼型附近被細化。采用動態(tài)網(wǎng)格結(jié)合用戶自定義函數(shù)（UDF）接口編程，采用彈簧方法進行網(wǎng)格重構(gòu)實現(xiàn)機翼的模態(tài)位移連續(xù)變形。在數(shù)值計算中，采用雙精度求解器進行數(shù)值模擬，采用Spalart-Allmaras湍流模型。

圖3 Agard445.6機翼氣動網(wǎng)格Fig.3 Aerodynamic grid of Agard445.6 wing

圖4所示為顫振馬赫數(shù)和參考文獻[24]的對比，可以看出，二者的數(shù)值和趨勢均具有良好的一直性，最大誤差小于8%。

圖4 顫振馬赫數(shù)對比驗證Fig.4 Flutter Mach number verification

圖5所示為動壓對應的馬赫數(shù)為0.29時Agard445.6機翼各階模態(tài)時域響應，從圖5中可以看出，各階模態(tài)位移隨著時間變化位移逐漸減少呈收斂趨勢，說明在該動壓下該機翼尚未達到顫振。通過分析離散狀態(tài)空間模型可以看出，最大特征值模長為0.994，小于1，該系統(tǒng)穩(wěn)定。通過分析連續(xù)狀態(tài)空間模型的特征值也可以看出，在該狀況下所有特征值的實部均小于0，該系統(tǒng)穩(wěn)定。圖6所示為動壓對應的馬赫數(shù)為0.3時Agard445.6機翼各階模態(tài)時域響應，從圖6 中可以看出，各階模態(tài)位移隨著時間變化位移基本保持不變，說明動壓對應的馬赫數(shù)為0.3 時非常接近顫振速度。圖7所示為動壓對應的馬赫數(shù)為0.31時Agard445.6機翼各階模態(tài)時域響應，從圖7中可以看出，各階模態(tài)位移隨著時間變化位移逐漸增大，說明該馬赫數(shù)已經(jīng)超過顫振速度了。通過分析離散狀態(tài)空間模型可以看出，最大特征值模長為1.0003，該系統(tǒng)發(fā)散。通過分析連續(xù)狀態(tài)空間模型的特征值也可以看出，在該狀況下所有特征值存在特征值0.2917±91.8737i實部大于0，該系統(tǒng)發(fā)散。

圖5 Ma=0.29時Agard445.6機翼各階模態(tài)時域響應Fig.5 Time domain response of various modes of Agard445.6 wing at Ma=0.29

圖6 Ma=0.3時Agard445.6機翼各階模態(tài)時域響應Fig.6 Time domain response of various modes of Agard445.6 wing at Ma=0.3

圖7 Ma=0.31時Agard445.6機翼各階模態(tài)時域響應Fig.7 Time domain response of various modes of Agard445.6 wing at Ma=0.31

3 柔性后緣可變形機翼氣動彈性分析計算

3.1 顫振分析計算

本節(jié)考慮柔性后緣機翼在不同偏角下的顫振速度計算，僅考慮氣動外形變化對氣動彈性特性的影響，柔性后緣可變形機翼的流體網(wǎng)格如圖8（a）所示。弦向變化范圍為翼型60%到后緣部分，這一段翼型中弧線變形采用作拋物線軌跡變彎，定義60%翼型處中弧線所在的位置和變形后的后緣點的連線與初始翼型中弧線的夾角為后緣偏角β，向下偏為正偏轉(zhuǎn)方向。

為節(jié)約網(wǎng)格繪制時間，對不同柔性后緣偏角情況采用同一套流體網(wǎng)格進行CFD 仿真計算。在計算某一特定偏角時，可以基于UDF程序預先將該部分后緣偏轉(zhuǎn)至給定的角度，圖8（b）所示為機翼柔性后緣偏轉(zhuǎn)6°時流體網(wǎng)格的變化圖。等到氣動力響應基本穩(wěn)定后，輸入對應的各階位移激勵，從而得到對應偏角下的各階廣義氣動力系數(shù)響應。

圖8 機翼柔性后緣偏轉(zhuǎn)前后網(wǎng)格變化Fig.8 Grid changes before and after flexible trailing edge deflection

在訓練馬赫數(shù)為0.9 時，后緣偏角為0°、2°、4°、6°時的顫振馬赫數(shù)分別為0.202、0.203、0.211、0.224?？梢钥闯觯S著后緣偏角的增大，顫振馬赫數(shù)相比無偏角情況有一定的提高。

3.2 陣風響應分析計算

取馬赫數(shù)為0.7 工況下計算1-cos 陣風響應，圖9 所示為陣風尺寸分別為60m、90m、140m 時柔性后緣偏角為6°下的機翼各階模態(tài)位移響應?？梢钥闯?，在該工況下，飛機遭遇陣風響應隨著陣風尺寸的增大，模態(tài)位移響應呈下降趨勢。

圖9 模態(tài)位移時域響應Fig.9 Time domain response of modal displacement

4 結(jié)論

本文提出一種基于氣動降階模型建立柔性后緣可變形機翼氣動彈性分析方法，應用于柔性后緣可變形機翼的顫振特性和陣風響應分析計算。通過結(jié)構(gòu)有限元方法得到各階模態(tài)，將各階模態(tài)變形導入CFD 計算，得到對應的氣動力，基于ARX 方法得到氣動狀態(tài)空間模型，耦合結(jié)構(gòu)狀態(tài)空間模型從而建立可變形機翼的氣動彈性狀態(tài)空間模型，用于可變形機翼氣動彈性顫振特性分析和陣風響應計算。研究結(jié)果表明，隨著后緣偏角的增大，顫振馬赫數(shù)相比無偏角情況有一定的提高。在馬赫數(shù)為0.7的工況下，隨著陣風尺寸的增大，模態(tài)位移響應呈下降趨勢。本文尚未考慮柔性后緣結(jié)構(gòu)變形對氣動彈性特性的影響，可在后續(xù)工作中綜合考慮結(jié)構(gòu)氣動影響的氣動彈性特性。