基于CFD/CSD耦合的結(jié)構(gòu)幾何非線性靜氣動(dòng)彈性數(shù)值方法研究

聶雪媛， 黃程德， 楊國偉(中國科學(xué)院 力學(xué)研究所 流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京　100190)

?

基于CFD/CSD耦合的結(jié)構(gòu)幾何非線性靜氣動(dòng)彈性數(shù)值方法研究

聶雪媛， 黃程德， 楊國偉(中國科學(xué)院 力學(xué)研究所 流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京100190)

摘要：柔性飛行器在氣動(dòng)力作用下會(huì)發(fā)生大變形，產(chǎn)生結(jié)構(gòu)幾何非線性，線性小變形方法難以獲得準(zhǔn)確的氣動(dòng)彈性分析結(jié)果。基于RANs的三維N-S流場(chǎng)控制方程耦合非線性結(jié)構(gòu)靜力學(xué)方程時(shí)域分析方法，用于考慮結(jié)構(gòu)幾何非線性的靜氣動(dòng)彈性分析。該方法在結(jié)構(gòu)靜力學(xué)方程求解上采用非線性增量有限元方法進(jìn)行迭代求解，考慮結(jié)構(gòu)剛度矩陣隨結(jié)構(gòu)位形的變化，采用徑向基函數(shù)方法實(shí)現(xiàn)氣動(dòng)/結(jié)構(gòu)界面的數(shù)據(jù)交換和動(dòng)網(wǎng)格變形。在建立某型寬體客機(jī)復(fù)材機(jī)翼三維有限元模型的基礎(chǔ)上，對(duì)其靜氣動(dòng)彈性進(jìn)行了數(shù)值仿真，分析了線性結(jié)構(gòu)和考慮結(jié)構(gòu)幾何非線性的結(jié)構(gòu)在靜氣動(dòng)彈性作用下翼面扭轉(zhuǎn)、展向位移、垂向位移以及升力系數(shù)等物理量。算例結(jié)果表明，與線性結(jié)果相比，非線性結(jié)構(gòu)由于結(jié)構(gòu)幾何非線性的影響，在展向和垂向變形上兩者存在顯著差異。為準(zhǔn)確進(jìn)行柔性結(jié)構(gòu)的氣動(dòng)彈性分析，必須考慮結(jié)構(gòu)幾何非線性的影響。

關(guān)鍵詞：靜氣動(dòng)彈性；結(jié)構(gòu)幾何非線性；柔性飛行器

隨著對(duì)航空工業(yè)領(lǐng)域需求的發(fā)展，大型運(yùn)輸機(jī)，高空長航時(shí)無人機(jī)得以廣泛應(yīng)用。這類飛行器往往具有大的展弦比結(jié)構(gòu)，并大量使用復(fù)合材料，在實(shí)現(xiàn)更高的升阻比和更輕的結(jié)構(gòu)質(zhì)量的同時(shí)也增大了結(jié)構(gòu)的柔性，導(dǎo)致結(jié)構(gòu)在氣動(dòng)載荷作用下發(fā)生大變形。此時(shí)在對(duì)這類飛行器進(jìn)行氣動(dòng)彈性分析時(shí)，結(jié)構(gòu)的幾何非線性效應(yīng)不容忽視。

當(dāng)前對(duì)考慮結(jié)構(gòu)幾何非線性的氣動(dòng)彈性穩(wěn)定性分析通常采用氣動(dòng)-結(jié)構(gòu)相耦合的時(shí)域數(shù)值分析方法，根據(jù)氣動(dòng)力計(jì)算方法大致可分為兩大類，一類為線性氣動(dòng)力耦合梁模型的方法，由Patil等[1-3]提出，采用片條理論或者偶極子格網(wǎng)法計(jì)算氣動(dòng)力，研究幾何非線性對(duì)幾何精確梁模型的氣動(dòng)彈性穩(wěn)定性的影響。楊智春等[4]采用偶極子格網(wǎng)法和梁模型對(duì)考慮結(jié)構(gòu)幾何非線性效應(yīng)的顫振特性進(jìn)行研究，研究結(jié)果表明結(jié)構(gòu)幾何非線性對(duì)其動(dòng)力學(xué)特性影響顯著；另一類為非線性氣動(dòng)力耦合非線性梁模型，Smith等[5]采用Euler求解器耦合幾何精確梁模型對(duì)大展弦比機(jī)翼進(jìn)行了氣動(dòng)彈性分析，Garcia等[6]采用N-S方程計(jì)算氣動(dòng)力，對(duì)細(xì)長機(jī)翼進(jìn)行了非線性氣動(dòng)彈性分析。

已有對(duì)考慮結(jié)構(gòu)幾何非線性的氣動(dòng)彈性分析的研究中，結(jié)構(gòu)模型幾乎都是將三維機(jī)翼簡化為一直梁模型進(jìn)行分析，無法體現(xiàn)實(shí)際工程應(yīng)用中的三維模型的真實(shí)彈性變形[7]，也無法準(zhǔn)確地展現(xiàn)復(fù)合材料在結(jié)構(gòu)中的分布和對(duì)結(jié)構(gòu)變形的影響。此外，線性氣動(dòng)力模型無法應(yīng)用于存在強(qiáng)激波或流動(dòng)分離的流場(chǎng)，無法體現(xiàn)黏性效應(yīng)。為此，本文提出基于N-S方程耦合三維結(jié)構(gòu)有限元模型的方法用于分析結(jié)構(gòu)幾何非線性對(duì)氣動(dòng)彈性性能的影響。首先以文獻(xiàn)[2]中大展弦比機(jī)翼為算例，驗(yàn)證了該方法的準(zhǔn)確性，然后以某型號(hào)寬體客機(jī)的復(fù)材機(jī)翼作為研究對(duì)象，在跨聲速流場(chǎng)，分析比較其線性和非線性結(jié)構(gòu)達(dá)到靜平衡位置后的變形情況。為增強(qiáng)結(jié)構(gòu)的幾何非線性效應(yīng)，本文通過改變來流的攻角來突出線性和非線性結(jié)構(gòu)變形的差異。

1CFD/CSD耦合方法

本文采用氣動(dòng)結(jié)構(gòu)耦合的時(shí)域分析方法將流體控制方程和結(jié)構(gòu)動(dòng)力學(xué)方程耦合起來進(jìn)行求解。這主要是由于流場(chǎng)和結(jié)構(gòu)采用不同的坐標(biāo)體系描述，難以做到統(tǒng)一求解。非定常氣動(dòng)力的計(jì)算采用N-S方程實(shí)現(xiàn)，結(jié)構(gòu)動(dòng)力學(xué)方程采用非線性有限元增量法求解。此外，由于氣動(dòng)和結(jié)構(gòu)方程的數(shù)值求解都是在離散時(shí)間步上進(jìn)行的，需要相應(yīng)時(shí)刻的結(jié)構(gòu)位移和氣動(dòng)力，這就需要在氣動(dòng)和結(jié)構(gòu)之間建立數(shù)據(jù)交換，即結(jié)構(gòu)節(jié)點(diǎn)位移插值到氣動(dòng)物面網(wǎng)格點(diǎn)；氣動(dòng)物面氣動(dòng)力插值到結(jié)構(gòu)節(jié)點(diǎn)。本文采用徑向基函數(shù)實(shí)現(xiàn)界面數(shù)據(jù)插值和流場(chǎng)動(dòng)網(wǎng)格變形，考慮結(jié)構(gòu)幾何非線性的靜氣動(dòng)彈性分析流程如圖1所示。

圖1　結(jié)構(gòu)幾何非線性靜氣動(dòng)彈性分析流程Fig.1 Static aeroelastic analysis flow chart for structural geometric nonlinearity

1.1氣動(dòng)控制方程

氣動(dòng)力計(jì)算采用基于RANS的三維N-S控制方程， 守恒型的流動(dòng)方程可表達(dá)為，

(1)

式中,Q 為守恒向量，Gc和Gv分別為對(duì)流通量和黏性通量，S為控制體V的邊界面積，n 代表面的法向量，t代表物理時(shí)間。將式(1)按有限體積法進(jìn)行空間離散可得，

(2)

式中，VI代表第I個(gè)控制體單元,NF代表包圍第I個(gè)單元的所有面數(shù)，ΔSm代表第m個(gè)表面的法向面積。

N-S方程組空間離散分為無黏項(xiàng)和黏性項(xiàng)。黏性項(xiàng)的離散采用二階中心差分格式，無黏項(xiàng)是非線性對(duì)流的集中體現(xiàn)，因此無黏項(xiàng)離散是空間離散的研究重點(diǎn)。為捕捉激波并避免求解過程中的數(shù)值振蕩[8]，對(duì)無黏項(xiàng)的離散采用Roe格式，即

式中，A為Roe矩陣，上標(biāo)R，L分別代表相關(guān)變量來自單元界面的右邊和左邊。

對(duì)于等式(2)采用二階三點(diǎn)近似，進(jìn)行時(shí)間離散可得，

(4)

式中，RI稱為殘差，代表等式(2)的右邊項(xiàng)。

為提高時(shí)間推進(jìn)精度，在等式(4)中引入虛擬時(shí)間，最終可得到以下時(shí)間離散形式，

式中，上標(biāo)‘*’ 代表虛擬時(shí)刻對(duì)應(yīng)的物理量，m為虛擬時(shí)間步。本文在虛時(shí)間域上采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel) 隱格式對(duì)氣動(dòng)方程做隱式時(shí)間推進(jìn)，當(dāng)子迭代步m→∞時(shí)，Q*(m)→Q(n+1)即可求得流場(chǎng)在下一物理時(shí)間步的解。

1.2結(jié)構(gòu)靜力學(xué)求解

結(jié)構(gòu)線性分析認(rèn)為結(jié)構(gòu)的剛度不會(huì)隨結(jié)構(gòu)的位形改變，因此對(duì)于線性結(jié)構(gòu)靜力學(xué)方程的解可以通過結(jié)構(gòu)剛度矩陣直接求得。當(dāng)結(jié)構(gòu)柔度增大產(chǎn)生大變形，引起結(jié)構(gòu)的幾何非線性時(shí)，結(jié)構(gòu)的剛度矩陣成為隨結(jié)構(gòu)位形變化的函數(shù)，傳統(tǒng)的小變形假設(shè)理論不再適用。

在涉及幾何非線性問題的有限元方法中，通常采用增量法，將要求解的問題分為若干步，求解每一步的位移，并更新結(jié)構(gòu)剛度，每一步都是以上一步的解作為該步的初始解。由于結(jié)構(gòu)幾何非線性效應(yīng)，使得在建立結(jié)構(gòu)的平衡條件時(shí)必須考慮變形后的位形，這就涉及到所參考的平衡位形。在實(shí)際分析中，有兩種選擇[9]，一種是以上一時(shí)刻求得的位形作為參考位形，稱之為更新拉格朗日格式(Updated Lagrange Formulation, UL)；一種是以時(shí)間0的位形作為參考位形，稱之為完全拉格朗日格式(Total Lagrange Formulation, TL)。本文采用TL格式對(duì)幾何非線性問題進(jìn)行描述，現(xiàn)在假定在時(shí)間0到t的所有時(shí)間點(diǎn)的解已經(jīng)求得，當(dāng)前需要求解t+Δt時(shí)刻各物理量，利用虛功原理，結(jié)構(gòu)平衡方程可表示為，

(6)

(7)

(8)

(9)

其中，

將式(7)～式(9)代入式(6)，并對(duì)平衡方程進(jìn)行線性化處理，最終導(dǎo)出用于TL的靜力學(xué)求解方程為,

(10)

(11)

對(duì)式(10)的求解采用增量法，在每一個(gè)增量步采用Newton-Raphson迭代算法，由于在每一次迭代中都需要進(jìn)行剛度矩陣的重新計(jì)算，當(dāng)?shù)諗亢螅M(jìn)入到下一個(gè)增量，直到增量載荷加載到t+ΔtQ。因此非線性結(jié)構(gòu)與線性結(jié)構(gòu)求解相比要花費(fèi)更多的計(jì)算時(shí)間。其具體求解過程可見圖1。

1.3數(shù)據(jù)交換和動(dòng)網(wǎng)格變形

為給氣動(dòng)控制方程提供運(yùn)動(dòng)邊界，同時(shí)給結(jié)構(gòu)方程提供節(jié)點(diǎn)載荷，以實(shí)現(xiàn)氣動(dòng)-結(jié)構(gòu)的耦合求解，需要在結(jié)構(gòu)和氣動(dòng)之間進(jìn)行數(shù)據(jù)轉(zhuǎn)換，即將結(jié)構(gòu)點(diǎn)位移插值到氣動(dòng)物面網(wǎng)格點(diǎn)，得到氣動(dòng)物面的網(wǎng)格變形，并將氣動(dòng)物面點(diǎn)的氣動(dòng)力變換為結(jié)構(gòu)節(jié)點(diǎn)載荷。此外，結(jié)構(gòu)變形也會(huì)帶動(dòng)整個(gè)空間流場(chǎng)的網(wǎng)格變形。由于該方法插值過程僅需要節(jié)點(diǎn)坐標(biāo)，而無需節(jié)點(diǎn)間連接信息，所以數(shù)據(jù)結(jié)構(gòu)簡單，是一種高效的數(shù)學(xué)插值方法[10]。在本文中采用基于徑向基函數(shù)[11](RBF)的插值方法用于上述氣動(dòng)-結(jié)構(gòu)之間的界面數(shù)據(jù)交換和動(dòng)網(wǎng)格變形，這是一種三維插值方法，其插值公式為，

(12)

式中,s(r)為在空間位置r處需要被估計(jì)的函數(shù)值，φ為徑向基函數(shù)，其函數(shù)類型表達(dá)可參見Wendland[12]，Ns是徑向基函數(shù)中心點(diǎn)個(gè)數(shù)，ri是中心點(diǎn)的空間位置，系數(shù)γi為插值系數(shù)。p(r)是一個(gè)可選的多項(xiàng)式，在界面數(shù)據(jù)交換時(shí)，其表達(dá)式為p(r)=a0+a1x+a2y+a3z以滿足力和力矩平衡的要求；在動(dòng)網(wǎng)格變形時(shí)，該表達(dá)式不出現(xiàn)在式(12)中。

RBF實(shí)現(xiàn)數(shù)據(jù)插值的過程可簡述如下。設(shè)

(13)

ax=(γ1…γNsa0a1a2a3)T

(14)

式中,Ns為用于RBF插值的結(jié)構(gòu)中心點(diǎn)個(gè)數(shù)，xi,yi,zi(i=1,…,Ns)為中心點(diǎn)坐標(biāo)。

由式(12)可知，結(jié)構(gòu)節(jié)點(diǎn)在x方向上的位移分量可表達(dá)為，

usx=Cssax

(15)

式中,Usx=(usx0000)T

根據(jù)對(duì)插值基函數(shù)的研究，本文選取的基函數(shù)如式(16)所述。

(16)

通過求解式(15)可得到系數(shù)矩陣式(14)，由此可計(jì)算出氣動(dòng)物面網(wǎng)格點(diǎn)的位移向量uax為，

uax=Casax

(17)

式中，

Na為氣動(dòng)物面點(diǎn)個(gè)數(shù)。

根據(jù)虛功原理，氣動(dòng)力Fa和結(jié)構(gòu)節(jié)點(diǎn)氣動(dòng)力Fs之間滿足以下關(guān)系

Fs·us=Fa·ua

(18)

利用式(15)和式(17)即可得氣動(dòng)物面氣動(dòng)力到結(jié)構(gòu)節(jié)點(diǎn)氣動(dòng)了轉(zhuǎn)換關(guān)系為，

(19)

在獲得氣動(dòng)物面點(diǎn)位移后，通過貪婪算法[13]選擇某些氣動(dòng)物面點(diǎn)用于動(dòng)網(wǎng)格變形的插值計(jì)算，整個(gè)插值的過程與上述結(jié)構(gòu)節(jié)點(diǎn)到氣動(dòng)物面網(wǎng)格點(diǎn)位移插值相類似。需要指出的是，在動(dòng)網(wǎng)格變形時(shí)，式(13)～式(15)中的后四行數(shù)據(jù)需要去掉，這是因?yàn)榫W(wǎng)格變形和CFD/CSD 數(shù)據(jù)交換需要滿足的要求是不同的，CFD/CSD 數(shù)據(jù)變換矩陣需要保證準(zhǔn)確的力和力矩平衡準(zhǔn)則而包含線性項(xiàng)，而網(wǎng)格變換矩陣為了避免邊界網(wǎng)格的移動(dòng)和非物理力的傳遞而不需要線性項(xiàng)。另外在基函數(shù)的選取上，需要引入緊支半徑R，用于控制中心點(diǎn)的作用范圍，即基函數(shù)形式為

(20)

由于流場(chǎng)網(wǎng)格點(diǎn)數(shù)量巨大，將流場(chǎng)網(wǎng)格分塊后，采用并行計(jì)算實(shí)現(xiàn)動(dòng)網(wǎng)格變形。

2算例與分析

為驗(yàn)證本文中所采用的計(jì)算非線性結(jié)構(gòu)靜氣動(dòng)彈性算法的正確性，本文首先以諸多文獻(xiàn)在研究結(jié)構(gòu)非線性時(shí)通常所選取的結(jié)構(gòu)即半展弦比為16∶1的機(jī)翼作為研究對(duì)象，有關(guān)該結(jié)構(gòu)的具體參數(shù)見文獻(xiàn)[2]。計(jì)算在飛行高度20 km，飛行速度25 m/s，攻角為2°時(shí)的結(jié)構(gòu)靜平衡位置，并與已有文獻(xiàn)[5]結(jié)果進(jìn)行了比較，如圖2所示。其中面元法和歐拉法為文獻(xiàn)[5]計(jì)算時(shí)采用的方法，CFD為本文計(jì)算氣動(dòng)力的方法。

圖2　梁結(jié)構(gòu)靜氣動(dòng)彈性變形Fig.2 Static aeroelastic deformation of large span ratio beam

從圖2可以看出，采用本文所提出方法對(duì)非線性結(jié)構(gòu)進(jìn)行靜氣動(dòng)彈性分析，其線性結(jié)構(gòu)和非線性結(jié)構(gòu)的變形規(guī)律與文獻(xiàn)[5]的研究相似，即非線性結(jié)構(gòu)在垂向上的變形大于線性結(jié)構(gòu)。需要指出的是，由于結(jié)構(gòu)是按照梁結(jié)構(gòu)進(jìn)行建模，因此在計(jì)算線性結(jié)構(gòu)的變形時(shí)，軸向是沒有位移的，這也是將機(jī)翼簡單視為梁結(jié)構(gòu)的不足之處。

在對(duì)方法進(jìn)行驗(yàn)證后，本文以某帶短艙的寬體客機(jī)復(fù)材機(jī)翼作為研究幾何非線性靜氣動(dòng)彈性的模型。該機(jī)翼展弦比12，梢根比0.2，后掠角30°，翼根弦長約13 m。整個(gè)機(jī)翼由復(fù)合材料按照0°，±45°和90°進(jìn)行鋪層，以增大結(jié)構(gòu)變形時(shí)的彎扭耦合效應(yīng)。在彈性結(jié)構(gòu)有限元分析時(shí)認(rèn)為機(jī)身為剛性的，只考慮機(jī)翼的彈性變形，其有限元模型如圖3所示。流場(chǎng)計(jì)算網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，利用ICEM CFD，以全機(jī)半模外形為基礎(chǔ)生成O-H型多塊對(duì)接結(jié)構(gòu)網(wǎng)格，在機(jī)翼與機(jī)身連接處，采用H型網(wǎng)格，網(wǎng)格在機(jī)翼的前后緣方向以及翼根處做了適當(dāng)加密，以便更好地分辨流場(chǎng)細(xì)節(jié)。為更好的刻畫邊界層內(nèi)的流動(dòng)，在壁面附近布置了邊界層網(wǎng)格，其首層高度為0.05 mm，邊界層最大厚度為400 mm。計(jì)算網(wǎng)格規(guī)模約為1 700萬，流場(chǎng)計(jì)算網(wǎng)格如圖4所示。湍流模型采用兩方程渦黏性模型—Wilcoxk-ω模型。

圖3　結(jié)構(gòu)有限元模型Fig.3 Structural finite element model

圖4　氣動(dòng)表面網(wǎng)格Fig.4 Surface mesh and symmetry plane mesh

該算例的計(jì)算條件為：Ma=0.85，Re= 8.264×107，飛行高度11 km，大氣密度0.363 92 kg/m3。

通過計(jì)算不同來流攻角下結(jié)構(gòu)變形，分析結(jié)構(gòu)幾何非線性對(duì)靜氣動(dòng)彈性特性的影響。分別計(jì)算了1°,5°和8°攻角下，線性結(jié)構(gòu)和非線性結(jié)構(gòu)的靜平衡位置，圖5給出了沿展向在不同攻角下的變形。圖6給出了在8°攻角下線性結(jié)構(gòu)和非線性結(jié)構(gòu)在靜平衡位置的變形比較。

圖5　不同攻角下結(jié)構(gòu)沿展向變形比較Fig.5 Spanwise displacement for different attack angle at Mach=0.85

圖6　Ma=0.85 8°攻角下結(jié)構(gòu)變形Fig.6 Deformation comparison for linear and nonlinear structure at Mach 0.85 and AOA 8°

從上述兩圖可看出，與線性結(jié)構(gòu)相比，幾何非線性在展向和垂向上的變形都大于線性結(jié)構(gòu)，而且在本計(jì)算攻角范圍內(nèi)，隨著攻角的增大變形增大。這與Mian[14]中，線性結(jié)構(gòu)在展向沒有位移的結(jié)論有所不同，分析認(rèn)為Mian所建立的結(jié)構(gòu)實(shí)質(zhì)為單一的梁結(jié)構(gòu)，而在梁結(jié)構(gòu)的線性計(jì)算中，是不考慮展向變形的。

圖7給出了在8°攻角時(shí)，結(jié)構(gòu)在靜平衡位置處，沿展向各剖面的彎矩和扭矩。可以看出非線性結(jié)構(gòu)所受的彎矩大于線性結(jié)構(gòu)，導(dǎo)致非線性結(jié)構(gòu)垂向位移和展向位移大于線性結(jié)構(gòu)。

圖7　Ma=0.85 8°攻角下靜平衡位置彎扭矩比較Fig.7 Bending and torsion moments comparisons for linear and nonlinear structure at Mach 0.85 and AOA 8°

圖8給出了在0.85馬赫數(shù)下，線性和非線性結(jié)構(gòu)升力系數(shù)隨攻角的變化，可以看出在本算例中，線性結(jié)構(gòu)在平衡位置的升力小于非線性結(jié)構(gòu)。這可從翼面扭轉(zhuǎn)角來分析，以8°攻角為例，圖9給出了沿展向不同剖面的相對(duì)扭轉(zhuǎn)角。無論是線性還是非線性結(jié)構(gòu)，由于彈性變形，后掠機(jī)翼沿流向各個(gè)剖面都會(huì)產(chǎn)生負(fù)扭轉(zhuǎn)，且該扭轉(zhuǎn)角沿展向增加而增大。在本例中，線性結(jié)構(gòu)在靜氣動(dòng)彈性作用下各剖面產(chǎn)生的負(fù)扭轉(zhuǎn)角大于非線性結(jié)構(gòu)，導(dǎo)致線性結(jié)構(gòu)飛行攻角小于非線性結(jié)構(gòu)，因此其升力系數(shù)小于非線性結(jié)構(gòu)。

圖8　升力系數(shù)隨攻角變化Fig.8 Lift coefficient curve with different angle of attack at Mach=0.85

圖9　沿展向剖面相對(duì)扭轉(zhuǎn)角Fig.9 Relative twist angle along spanwise

3結(jié)論

本文針對(duì)柔性飛行器在氣動(dòng)力作用下呈現(xiàn)出的結(jié)構(gòu)幾何非線性特點(diǎn)，采用非線性氣動(dòng)力耦合非線性結(jié)構(gòu)的時(shí)域分析方法，對(duì)考慮結(jié)構(gòu)幾何非線性效應(yīng)的靜氣動(dòng)彈性進(jìn)行了數(shù)值仿真。基于RANs的N-S流場(chǎng)求解器，在虛擬時(shí)間歩上采用LU-SGS 方法和 Newton-Raphson方法交替求解非線性氣動(dòng)力和非線性結(jié)構(gòu)變形的緊耦合方式，徑向基函數(shù)插值方法被用來進(jìn)行動(dòng)網(wǎng)格變形和耦合界面數(shù)據(jù)插值。運(yùn)用所發(fā)展的算法對(duì)已有文獻(xiàn)算例進(jìn)行了方法驗(yàn)證，然后對(duì)某寬體客機(jī)復(fù)材機(jī)翼的靜氣動(dòng)彈性開展研究，比較了線性和非線性結(jié)構(gòu)的靜氣動(dòng)彈性特性。針對(duì)本文算例，可以得到以下結(jié)論:

(1) 無論是線性還是非線性結(jié)構(gòu)，對(duì)于后掠角機(jī)翼，由于彈性變形的影響，機(jī)翼沿流向各個(gè)剖面會(huì)產(chǎn)生負(fù)的扭轉(zhuǎn)角，且隨著展向的增大，各剖面扭轉(zhuǎn)角增大，機(jī)翼飛行攻角減小；

(2) 非線性結(jié)構(gòu)因變形產(chǎn)生的各剖面負(fù)扭轉(zhuǎn)角小于線性結(jié)構(gòu)，因此線性結(jié)構(gòu)在靜平衡位置時(shí)的升力小于非線性結(jié)構(gòu)；

(3) 采用CFD與結(jié)構(gòu)耦合的方法進(jìn)行靜氣動(dòng)彈性分析，在結(jié)構(gòu)達(dá)到平衡位置時(shí)，線性結(jié)構(gòu)在展向上有變形；非線性結(jié)構(gòu)的展向位移大于線性結(jié)構(gòu)；在相同展向站位處，非線性結(jié)構(gòu)的垂向變形大于線性結(jié)構(gòu)。

參 考 文 獻(xiàn)

[ 1 ] Patil M J, Hodges D H. Limit cycle oscillations in high-aspect-ratio Wings[C]//AIAA-1099-1464,April12-15,1999,St.Louis,Mo,USA.

[ 2 ] Patil M J, Hodges D H. On the importance of aerodynamic and structural nonlinearities in aeroelastic behavior of high-aspect-ratio wings[J]. Journal of Fluids and Structures,2014,19(7):905-915.

[ 3 ] Patil M J, Hodges D H, Cesnik C E S. Characterizing the effects of geometrical nonlinearities on aeroelastic behavior of high-aspect-ratio wings[C]// International Forum on Aeroelasticity and Structural Dynamics, June 22-25,1999,Williamsburg, VA,USA.

[ 4 ] 楊智春, 張惠，谷迎松,等.考慮幾何非線性效應(yīng)的大展弦比機(jī)翼氣動(dòng)彈性分析[J].振動(dòng)與沖擊，2014,33(16):72-75.

YANG Zhi-chun, ZHANG Hui, GU Ying-song, et al. Aeroelastic analysis of the high aspect ratio wing considering the geometric nonlinearity[J]. Journal of Vibration and Shock, 2014,33(16):72-75.

[ 5 ] Simth M J, Patil M J,Hodges D H. CFD-based analysis of nonlinear aeroelastic behavior of high-aspect ratio wing [C]// AIAA 2001-1582,April 16-19,2001,Seattle, WA USA.

[ 6 ] Garcia J A,Guruswamy G P. Aeroelastic analysis of transonic wings using Navier-Stokes equations and a nonlinear beam finite element model[C]//AIAA 1999-1215,April 12-15,1999, St. Louis,MO,USA.

[ 7 ] Palacios R,Cesnik C E S. Static nonlinear aeroelasticity of flexible slender wings in compressible flow[C]//AIAA 2005-1945,April 18-21,2005,Austin,TX,USA.

[ 8 ] Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics,1981,43:357-372.

[ 9 ] 王勖成.有限單元法[M].北京：清華大學(xué)出版社，2002.

[10] 張偉偉，高傳強(qiáng)，葉正寅. 氣動(dòng)彈性計(jì)算中網(wǎng)格變形方法研究進(jìn)展[J]. 航空學(xué)報(bào),2014(35):303-319.

ZHANG Wei-wei,GAO Chuan-qiang,YE Zheng-yin. Research progress on mesh deformation method in computational aeroelasticity[J]. Acta Aeronautica et Astronautica Sinica,2014(35):303-319.

[11] Michler A K. Aircraft control surface deflection using RBF-based mesh deformation[J]. International Journal for Numerical Methods in Engineering,2001,88:996-1007.

[12] Beckert A, Wendland H. Multivariate interpolation for fluid-structure interaction problems using radial basis functions[J]. Aerospace Science and Technology,2001,5(2):125-134.

[13] Rendall T C S,Allen C B. Reduced surface point selection options for efficient mesh deformation using radial basis functions[J]. Journal of Computational Physics,2010,229: 2810-2820.

[14] Mian H H, Wang Gang, Ye Zheng-yin. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach[J]. Journal of Fluids and Structures, 2014,49:186-201.

Numerical analysis for aeroelastic with structural geometrical nonlinearity using a CFD/CSD-coupled method

NIEXue-yuan,HUANGCheng-de，YANGGuo-wei(Key Laboratory for Mechanics in Fluid Solid Coupling Systems of Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China)

Abstract:Flexible aircrafts tend to undergo large deformations under aerodynamic forces.As a result, structural geometric nonlinearity occurs.The linear small-deformation theory cannot provide accurate results for analyzing static aeroelasticity.Fluid structure-coupling method based on three-dimensional RANs, Navier-Stokes equations and nonlinear static equation is used for static aeroelastic analysis with structural geometric nonlinearity.The mentioned approach adopts the nonlinear incremental finite-element method to solve nonlinear static equations with assembled structure stiffness matrixes.Moreover, RBF method is used for data interpolation and mesh deformation.Based on the multi-material wing finite-element model, the numerical simulations were made to analyze the static aeroelastic behavior.Comparisons of twist angles, vertical displacements, spanwise displacements and life coefficients between linear and nonlinear structures were made.The results show that geometric nonlinearity cannot be neglected for predicting accurate static aeroelastic behavior for large, flexible airplanes.

Key words:static aeroelastic analysis; geometric nonlinearity; flexible aircrafts

中圖分類號(hào):V211

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.08.008

通信作者楊國偉 男，博士，研究員，1967年生

收稿日期：2015-07-28修改稿收到日期：2015-09-24

基金項(xiàng)目：國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2014AA110501)

第一作者 聶雪媛 女，博士，助理研究員，1978年生

E-mail:gwyang@imech.ac.cn