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

趙仕偉，闞梓，李道春

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

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

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

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

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

1 氣動彈性狀態空間模型

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

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

式中

狀態矢量如下

對于多自由度系統，通過歸一化可以轉換為模態坐標運動，結構模態坐標下的動力學方程可表示如下

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

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

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

式中

狀態矢量xs(t)為

將上述連續時間系統狀態空間結構模型進行離散化處理，得到對應的離散狀態空間模型如下

其中：

式中，T為時間步長，上述過程為系統的結構動力學模型的獲取流程，廣義氣動力為系統輸入，結構彈性變形為系統輸出。

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

式中，u(k)為系統的廣義結構位移；yA(k)為系統的廣義氣動力。

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

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

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

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

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

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

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

表1 前四階模態的振動頻率與試驗結果對比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機翼前四階模態振型圖Fig.2 The first four modal shapes of Agard445.6 wing

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

圖3 Agard445.6機翼氣動網格Fig.3 Aerodynamic grid of Agard445.6 wing

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

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

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

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

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

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

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

3.1 顫振分析計算

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

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

圖8 機翼柔性后緣偏轉前后網格變化Fig.8 Grid changes before and after flexible trailing edge deflection

在訓練馬赫數為0.9 時，后緣偏角為0°、2°、4°、6°時的顫振馬赫數分別為0.202、0.203、0.211、0.224。可以看出，隨著后緣偏角的增大，顫振馬赫數相比無偏角情況有一定的提高。

3.2 陣風響應分析計算

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

圖9 模態位移時域響應Fig.9 Time domain response of modal displacement

4 結論

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