謝 萌程 攀

1.上海市浦東新區(qū)金科路5188號(hào)總體氣動(dòng)部；2.上海飛機(jī)設(shè)計(jì)研究院

基于兩種不同方法的民用客機(jī)靜氣動(dòng)彈性導(dǎo)數(shù)計(jì)算

謝 萌1程 攀2

1.上海市浦東新區(qū)金科路5188號(hào)總體氣動(dòng)部；2.上海飛機(jī)設(shè)計(jì)研究院

風(fēng)洞試驗(yàn)或是CFD計(jì)算過(guò)程中，一般使用剛性模型。這就存在一個(gè)問(wèn)題，就是如何根據(jù)剛性氣動(dòng)力修正到彈性氣動(dòng)力。本文使用兩種方法對(duì)大型民用客機(jī)的靜氣動(dòng)彈性導(dǎo)數(shù)修正系數(shù)做了計(jì)算分析，旨在為大型民用客機(jī)提供更加準(zhǔn)確的氣動(dòng)力修正系數(shù)。

隨著飛行器結(jié)構(gòu)重量的不斷減小、智能以及復(fù)合材料的大量運(yùn)用，飛行器氣動(dòng)彈性問(wèn)題在飛行器設(shè)計(jì)中也越來(lái)越突出。尤其是具有大展弦比機(jī)翼的民用飛機(jī)。

具有大展弦比機(jī)翼的民用大型客機(jī)在氣動(dòng)設(shè)計(jì)過(guò)程中給出的是在巡航飛行中的最佳外形，稱(chēng)為設(shè)計(jì)外形或者1G外形，而飛機(jī)的制造外形叫做型架外形。同時(shí)，這兩種外形又都存在剛性和彈性?xún)煞N情況。在實(shí)際情況中，飛機(jī)是一個(gè)彈性結(jié)構(gòu)，但是在風(fēng)洞試驗(yàn)過(guò)程中，由于模型的縮比，不是很方便根據(jù)縮比模擬真實(shí)的飛機(jī)結(jié)構(gòu)，當(dāng)然，純粹針對(duì)靜氣動(dòng)彈性的試驗(yàn)除外，但試驗(yàn)花費(fèi)較高，風(fēng)險(xiǎn)相對(duì)比較大，所以，在風(fēng)洞試驗(yàn)過(guò)程中，一般使用剛性模型。這就存在一個(gè)問(wèn)題，就是如何根據(jù)剛性氣動(dòng)力修正到彈性氣動(dòng)力，這種問(wèn)題針對(duì)CFD的氣動(dòng)力計(jì)算結(jié)果也存在。根據(jù)上面所描述的問(wèn)題來(lái)看，靜氣動(dòng)彈性導(dǎo)數(shù)計(jì)算非常重要。本文使用了基于線性和非線性氣動(dòng)力兩種方法做計(jì)算分析，目的是尋找更加高效且準(zhǔn)確的方法。

基于線性氣動(dòng)力分析的靜彈導(dǎo)數(shù)修正

飛行器彈性變形對(duì)定常氣動(dòng)載荷分布會(huì)產(chǎn)生影響，結(jié)構(gòu)彈性變形會(huì)引起氣動(dòng)載荷的重新分布，氣動(dòng)載荷的重新分布不僅會(huì)導(dǎo)致內(nèi)部結(jié)構(gòu)載荷和應(yīng)力的重新分布，也會(huì)改變空氣動(dòng)力穩(wěn)定性導(dǎo)數(shù)。靜氣動(dòng)彈性問(wèn)題分析的目的是獲得結(jié)構(gòu)和空氣動(dòng)力方面的數(shù)據(jù)，結(jié)構(gòu)方面的數(shù)據(jù)包括載荷、變形和應(yīng)力；空氣動(dòng)力方面的數(shù)據(jù)包括穩(wěn)定性和控制導(dǎo)數(shù)、配平條件、氣動(dòng)壓力。靜氣動(dòng)彈性導(dǎo)數(shù)修正系數(shù)計(jì)算分析是靜氣動(dòng)彈性專(zhuān)業(yè)中很重要的一部分工作，為操穩(wěn)專(zhuān)業(yè)和載荷專(zhuān)業(yè)提供重要數(shù)據(jù)。

穩(wěn)定性導(dǎo)數(shù)是飛行器設(shè)計(jì)過(guò)程中重要的設(shè)計(jì)約束，它是馬赫數(shù)和動(dòng)壓的函數(shù)。根據(jù)已知的配平自由度和配平變量通過(guò)求解平衡方程得到未知配平自由度和配平變量的過(guò)程稱(chēng)為配平分析，利用配平計(jì)算可得到穩(wěn)定性導(dǎo)數(shù)?？捎?jì)算出飛機(jī)各氣動(dòng)彈性導(dǎo)數(shù)的剛性及彈性值，彈性值與剛性值的比值即為各氣動(dòng)彈性導(dǎo)數(shù)的修正系數(shù)值。

基于VLM線性氣動(dòng)力計(jì)算方法

渦格法是一種比較實(shí)用的數(shù)值計(jì)算方法，基于求解拉普拉斯方程，其理論基礎(chǔ)和平板面元法相同。這種方法具有很強(qiáng)的工程運(yùn)用性，計(jì)算效率高，結(jié)果準(zhǔn)確。對(duì)于工程應(yīng)用來(lái)講，這種方法有極大的優(yōu)勢(shì)。

渦格法采用的計(jì)算模型是將機(jī)翼投影到其基本平面，翼型可以簡(jiǎn)化成由厚度和彎度來(lái)表示。

然后分別沿展向和弦向分布離散的馬蹄渦，使整個(gè)機(jī)翼用有限多個(gè)離散的馬蹄渦系來(lái)代替。

圖1 翼型簡(jiǎn)化示意圖

由于我們主要研究亞音速和跨音速情況，所以將每個(gè)馬蹄渦上面附著的渦線于該網(wǎng)格面元的1/4弦線重合，若是超音速情況，則位置有所不同。兩條自由渦線從1/4弦線的兩個(gè)端點(diǎn)沿x軸向下無(wú)限延伸。這些布了渦的網(wǎng)格就稱(chēng)為渦格，此時(shí)的氣動(dòng)模型稱(chēng)為渦格模型。在每個(gè)渦格3/4弦線的中點(diǎn)取為控制點(diǎn)，在這些點(diǎn)上計(jì)算全部離散馬蹄渦引起的誘導(dǎo)速度，需滿(mǎn)足翼面上無(wú)穿透速度的邊界條件。

假設(shè)共有n個(gè)網(wǎng)格，使用Biot－Savart定理求出第i個(gè)面元上單位強(qiáng)度的馬蹄渦對(duì)第j個(gè)控制點(diǎn)誘導(dǎo)速度，建立線性方程組第j個(gè)控制點(diǎn)處的總速度，在第j個(gè)控制點(diǎn)處需要滿(mǎn)足物面邊界條件，求解線性方程組，得到渦強(qiáng)Γ，由Kutta－Joukowski定理求出該面元所受的力F，然后求出總的氣動(dòng)力。

圖2 機(jī)翼渦格法模型示意圖

圖3 基礎(chǔ)坐標(biāo)系下邊界條件分析圖

圖4 機(jī)翼平板面元示意圖

圖5 基于VLM的機(jī)翼網(wǎng)格圖、控制點(diǎn)位置示意圖、面元?dú)鈩?dòng)力示意圖

機(jī)翼氣動(dòng)力模型

通過(guò)上面的原理，可以通過(guò)matlab程序代碼來(lái)實(shí)現(xiàn)對(duì)機(jī)翼氣動(dòng)力的求解。分為三個(gè)部分來(lái)完成。用一個(gè)簡(jiǎn)單的模型作為說(shuō)明。

首先是幾何輸入，輸入機(jī)翼在投影平面上的10個(gè)輪廓點(diǎn)，將點(diǎn)連接起來(lái)，構(gòu)成機(jī)翼平板。然后在機(jī)翼上劃分網(wǎng)格，得到所有組成機(jī)翼的面元。劃分出面元后，在每個(gè)面元上布馬蹄渦，并找出氣動(dòng)力作用點(diǎn)，以及下洗速度控制點(diǎn)。如下圖7示，紅色點(diǎn)表示的是控制點(diǎn)所在位置。通過(guò)Biotsavart定律以及Kutta－Joukowski定理求出每個(gè)面元的氣動(dòng)力，最后得到整個(gè)機(jī)翼的氣動(dòng)力。

本文中某大型客機(jī)的機(jī)翼氣動(dòng)面元模型如圖8所示。

單梁結(jié)構(gòu)有限元分析

針對(duì)大展弦比飛機(jī)，將結(jié)構(gòu)簡(jiǎn)化為單梁模型，建立結(jié)構(gòu)的單梁結(jié)構(gòu)有限元模型，結(jié)構(gòu)模型中包括了機(jī)翼、機(jī)身、垂尾、平尾、短艙及吊掛等主要部件，結(jié)構(gòu)有限元用變截面梁的單元形式來(lái)實(shí)現(xiàn)，將經(jīng)過(guò)集中處理過(guò)的質(zhì)量節(jié)點(diǎn)通過(guò)RBE3卡與相鄰的2個(gè)剛軸節(jié)點(diǎn)連接。機(jī)翼單梁的結(jié)構(gòu)點(diǎn)變換為上反平面內(nèi)原始結(jié)構(gòu)點(diǎn)在新剛軸上的投影點(diǎn)。

機(jī)身部分，前機(jī)身和后機(jī)身結(jié)構(gòu)為梁模型，其中后機(jī)身為兩段式拐折梁，機(jī)身和機(jī)翼連接處以及機(jī)身和尾翼連接處，為模擬真實(shí)的結(jié)構(gòu)連接關(guān)系，采用了柔度陣模型進(jìn)行模擬。

機(jī)翼機(jī)身使用剛性連接方式，同時(shí)平尾和機(jī)身之間還有輔助支撐的彈簧連接單元。另外由于質(zhì)量節(jié)點(diǎn)與結(jié)構(gòu)節(jié)點(diǎn)不重合，考慮加權(quán)平均的方法將質(zhì)量及結(jié)構(gòu)點(diǎn)關(guān)聯(lián)起來(lái)。

配平分析及靜彈導(dǎo)數(shù)系數(shù)修正計(jì)算

穩(wěn)定性導(dǎo)數(shù)是馬赫數(shù)和動(dòng)壓的函數(shù)，對(duì)于無(wú)約束氣動(dòng)彈性方程（無(wú)約束穩(wěn)定性導(dǎo)數(shù)），簡(jiǎn)寫(xiě)的靜氣動(dòng)彈性配平方程，是配平參數(shù)向量，即配平變量，是支持點(diǎn)的剛體運(yùn)動(dòng)加速度向量，即為配平變量，、是各相關(guān)系數(shù)的組合，通常配平自由度和配平變量可以分為兩組：第一組由已知的配平自由度和配平變量組成，在靜彈分析中經(jīng)常使用的詞語(yǔ)“給定機(jī)動(dòng)條件”就是由那些給定的配平變量和配平自由度定義的；第二組由余下的配平自由度和配平變量組成，是未知的，需要通過(guò)上面的配平方程進(jìn)行求解。

圖6 全機(jī)單梁結(jié)構(gòu)有限元圖

圖7 不同高度下無(wú)尾隨馬赫數(shù)Ma的變化

基于非線性氣動(dòng)力分析的靜彈導(dǎo)數(shù)分析

基于線性氣動(dòng)力的計(jì)算分析具有高效率的絕對(duì)優(yōu)勢(shì)，因?yàn)獒槍?duì)不同馬赫數(shù)和高度，在配平分析以及導(dǎo)數(shù)計(jì)算中會(huì)有上百個(gè)狀態(tài)，但是，考慮到氣動(dòng)力的非線性性，針對(duì)具有特征性的狀態(tài)，例如針對(duì)本文中的大型客機(jī)，選取了0.785和0.82兩個(gè)馬赫數(shù)，采用了基于非線性氣動(dòng)力的CFD計(jì)算和耦合方法來(lái)實(shí)現(xiàn)對(duì)靜氣彈導(dǎo)數(shù)的計(jì)算。

若使用高效的靜氣動(dòng)彈性分析軟件，例如NASTRAN，運(yùn)用簡(jiǎn)化的線性氣動(dòng)力計(jì)算方法，可以非常迅速的得到彈性模型和剛性模型的穩(wěn)定性導(dǎo)數(shù)比值，但是對(duì)于科研工作來(lái)說(shuō)，還需要探索更多的方法，這就是本文中的第二種方法，使用基于非定常氣動(dòng)力計(jì)算，耦合結(jié)構(gòu)模型，來(lái)得到想要的結(jié)果。

氣動(dòng)力計(jì)算

這種基于CFD計(jì)算直接耦合的方法要結(jié)合型架外形的設(shè)計(jì)中的工作，通過(guò)給定指定馬赫數(shù)下不同的攻角，得到相應(yīng)的氣動(dòng)力系數(shù)，由這些攻角以及氣動(dòng)力系數(shù)，就可以得到關(guān)于攻角的導(dǎo)數(shù)，同時(shí)計(jì)算出剛性情況下的導(dǎo)數(shù)，這樣就知道了彈性模型和剛性模型的比值。這種方法，原理上很直接，優(yōu)點(diǎn)明顯，缺點(diǎn)也很明顯，就是耗時(shí)。

本文使用商業(yè)軟件CFX求解雷諾平均N－S方程，采用K－ω SST湍流模型。CFX是耦合隱式求解器，通過(guò)合理設(shè)置邊界條件和求解參數(shù)，來(lái)保證計(jì)算精度和收斂性，得到符合實(shí)際問(wèn)題的解。

這種基于貼體的結(jié)構(gòu)網(wǎng)格的雷諾平均N－S方程，作為主控方程，其在直角坐標(biāo)系Oxyz中表示為：

公式的右邊表示了流場(chǎng)中的粘性項(xiàng)。

分析說(shuō)明

圖8 0.785馬赫數(shù)2度攻角時(shí)機(jī)翼變形過(guò)程圖

圖9 0.785馬赫數(shù)各攻角升力系數(shù)值

下面以為例來(lái)說(shuō)明?；诜蔷€性氣動(dòng)力，直接耦合計(jì)算的方法，以型架外形為基礎(chǔ)，通過(guò)耦合迭代，當(dāng)機(jī)翼的變形趨于穩(wěn)定時(shí)，取出一系列的的值，其斜率即為該狀態(tài)下的彈性導(dǎo)數(shù)值。這種方法，通過(guò)精確計(jì)算每一步的氣動(dòng)力，在耦合結(jié)構(gòu)，產(chǎn)生結(jié)構(gòu)位移，機(jī)翼外形發(fā)生變化，再次進(jìn)行精確氣動(dòng)力計(jì)算，耦合結(jié)構(gòu)，依次這種迭代，直至變形穩(wěn)定，從理論上講，這種方法已經(jīng)盡力做到了精確，但是對(duì)于一個(gè)點(diǎn)的導(dǎo)數(shù)，需要進(jìn)行大量的計(jì)算，以及人工干預(yù)，非常的耗時(shí)耗力。下面的結(jié)果是馬赫數(shù)0.785，約11000m高度下的結(jié)果。

圖9中各剛性值是型架外形的各狀態(tài)結(jié)果，由圖9和圖8可以明顯看出，一個(gè)導(dǎo)數(shù)修正系數(shù)所需要的工作量。

本文選擇0.785馬赫數(shù)，約11000m高度這個(gè)比較典型的狀態(tài)來(lái)說(shuō)明，此時(shí)，可以知道此時(shí)升力攻角導(dǎo)數(shù)修正系數(shù)的值為0.9004。

計(jì)算分析

在飛機(jī)設(shè)計(jì)過(guò)程中，需要許多狀態(tài)下的導(dǎo)數(shù)修正系數(shù)的結(jié)果，由于上面的結(jié)果中沒(méi)有所需高度時(shí)的結(jié)果，以其它各狀態(tài)結(jié)果值為基礎(chǔ)，插值得到0.785馬赫數(shù)，所需高度下的導(dǎo)數(shù)修正系數(shù)的結(jié)果為0.910。這樣來(lái)看，基于線性氣動(dòng)力計(jì)算所得到的導(dǎo)數(shù)結(jié)果是可用的。

氣動(dòng)力的靜彈修正是飛機(jī)設(shè)計(jì)中非常重要的一部分工作，如何使用更加合適的方法獲得更加準(zhǔn)確的氣動(dòng)力是今后工作中需要繼續(xù)研究的問(wèn)題。