趙仕偉，闞梓，李道春

北京航空航天大學(xué)，北京 100191

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

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

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

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

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

非定常氣動(dòng)力ARX模型可以寫成如下表達(dá)式

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

式中

狀態(tài)矢量如下

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

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

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

將上述模態(tài)坐標(biāo)下的動(dòng)力學(xué)方程寫成狀態(tài)空間表達(dá)形式，可以得到

式中

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

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

其中：

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

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

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

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

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

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

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

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

2 氣動(dòng)彈性特性分析方法驗(yàn)證

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

表1 前四階模態(tài)的振動(dòng)頻率與試驗(yàn)結(jié)果對(duì)比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機(jī)翼前四階模態(tài)振型圖Fig.2 The first four modal shapes of Agard445.6 wing

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

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

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

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

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

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

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

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

3 柔性后緣可變形機(jī)翼氣動(dòng)彈性分析計(jì)算

3.1 顫振分析計(jì)算

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

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

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

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

3.2 陣風(fēng)響應(yīng)分析計(jì)算

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

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

4 結(jié)論

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