基于流-固耦合算法的跨/超聲速曲壁板氣動(dòng)彈性分析

梅冠華， 張家忠， 康 燦

(1. 江蘇大學(xué) 能源與動(dòng)力工程學(xué)院，江蘇 鎮(zhèn)江 212013； 2. 西安交通大學(xué) 能源與動(dòng)力工程學(xué)院，西安 710049)

?

基于流-固耦合算法的跨/超聲速曲壁板氣動(dòng)彈性分析

梅冠華1， 張家忠2， 康 燦1

(1. 江蘇大學(xué) 能源與動(dòng)力工程學(xué)院，江蘇 鎮(zhèn)江 212013； 2. 西安交通大學(xué) 能源與動(dòng)力工程學(xué)院，西安 710049)

采用流-固耦合算法研究了曲壁板在跨/超聲速氣流下的氣動(dòng)彈性特征。首先，給出了曲壁板的Von Kármán幾何大變形運(yùn)動(dòng)方程，并對(duì)其進(jìn)行了標(biāo)準(zhǔn)有限元離散。然后，簡(jiǎn)述了流動(dòng)的控制方程、數(shù)值解法、動(dòng)網(wǎng)格和流-固耦合方式。最后，對(duì)曲壁板的氣動(dòng)彈性響應(yīng)進(jìn)行了數(shù)值模擬和分析。結(jié)果表明彎曲造成壁板的初始?xì)鈩?dòng)載荷非零，使得其與平壁板的氣動(dòng)彈性特征大相徑庭：穩(wěn)定的曲壁板存在靜氣動(dòng)彈性變形；馬赫數(shù)為2時(shí)，曲壁板失穩(wěn)后，其顫振負(fù)向峰值遠(yuǎn)大于正向峰值，且彎曲高度的增大會(huì)誘發(fā)混沌型的振動(dòng)；馬赫數(shù)為0.8和0.9時(shí)，曲壁板僅會(huì)發(fā)生正向變形；馬赫數(shù)為1.2時(shí)，失穩(wěn)后曲壁板的顫振中心偏向了負(fù)方向。所得結(jié)果為高速飛行器的壁板設(shè)計(jì)和顫振抑制提供了依據(jù)，所提算法可推廣應(yīng)用于其它氣動(dòng)彈性問(wèn)題的數(shù)值分析。

壁板顫振；流-固耦合；氣動(dòng)彈性

當(dāng)前，隨著超聲速和高超聲速飛行器的研發(fā)熱潮，壁板氣動(dòng)彈性顫振引起了學(xué)者們的廣泛關(guān)注。壁板顫振是指在外部高速氣流的沖擊作用下飛行器表面的壁板結(jié)構(gòu)所誘發(fā)的自激振動(dòng)現(xiàn)象，這是典型的由慣性力、彈性力和氣動(dòng)力的共同作用而激發(fā)的氣動(dòng)彈性問(wèn)題。該類振動(dòng)往往具有很強(qiáng)的非線性動(dòng)力學(xué)特性，并將對(duì)飛行器的疲勞壽命、飛行性能、飛行安全和乘坐品質(zhì)帶來(lái)不利影響。該問(wèn)題的深入探討對(duì)于高速飛行器的壁板設(shè)計(jì)、顫振抑制和疲勞壽命估計(jì)都具有十分重要的意義[1-2]。

作為經(jīng)典的氣動(dòng)彈性問(wèn)題，壁板顫振已被眾多學(xué)者進(jìn)行了深入研究，精確高效的數(shù)值方法被不斷引入，如Galerkin方法[3-4]、有限元方法[5-6]、Rayleigh-Ritz方法[7]、時(shí)滯慣性流形方法[8-9]等。所研究的壁板材料也更加貼近實(shí)際，如復(fù)合材料層合壁板[10]、功能梯度材料壁板[11]等。雖然壁板顫振的研究已經(jīng)取得了豐碩的成果，然而在這些研究中，氣動(dòng)載荷大多采用簡(jiǎn)化氣動(dòng)力模型近似表達(dá)，如線性/非線性活塞理論[12]、線性化勢(shì)流理論等，雖然它們便于應(yīng)用，然而其適用范圍和精度都較為有限，尤其無(wú)法對(duì)非線性效應(yīng)特別強(qiáng)烈的跨聲速流動(dòng)進(jìn)行描述。此外，簡(jiǎn)化氣動(dòng)力理論給出的氣動(dòng)載荷是與固體位移、速度及其偏導(dǎo)數(shù)相關(guān)的函數(shù)，而并非是通過(guò)基于流動(dòng)控制方程的計(jì)算流體動(dòng)力學(xué)(CFD)所得出的，因此無(wú)法給出流場(chǎng)的細(xì)節(jié)信息，比如黏性邊界層、旋渦、流動(dòng)分離、激波等，進(jìn)而無(wú)法準(zhǔn)確分析流體與固體間的耦合機(jī)理。事實(shí)上，這些正是壁板顫振定量分析中的開放性問(wèn)題。

近年來(lái)，伴隨著計(jì)算機(jī)軟硬件技術(shù)的飛速發(fā)展，學(xué)者們開始嘗試將CFD和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(CSD)相結(jié)合，發(fā)展流-固耦合算法來(lái)分析壁板的氣動(dòng)彈性特征，如此不僅能更加精確地反映該問(wèn)題的本質(zhì)，而且為跨聲速氣流下該問(wèn)題的求解提供了有效手段。DAVIS等[13-14]基于Euler方程的氣動(dòng)力理論，采用流-固耦合算法分析了跨聲速氣流中的二維壁板顫振問(wèn)題，并與簡(jiǎn)化氣動(dòng)力理論所得結(jié)果進(jìn)行了對(duì)比，說(shuō)明了Euler氣動(dòng)理論在分析跨聲速壁板顫振問(wèn)題時(shí)的精確性。隨后，GORDINER等[15]采用Navier-Stokes方程求解氣動(dòng)載荷，對(duì)三維壁板顫振問(wèn)題進(jìn)行了研究，發(fā)現(xiàn)邊界層的存在推遲了顫振的發(fā)生。之后，HASHIMOTO等[16]同樣采用Navier-Stokes氣動(dòng)力理論分析了跨聲速三維壁板顫振問(wèn)題，在與相關(guān)實(shí)驗(yàn)結(jié)果[17]取得一致的同時(shí)，其研究還發(fā)現(xiàn)湍流邊界層不僅對(duì)于顫振可以起到穩(wěn)定作用，在特定條件下還能引起不穩(wěn)定效應(yīng)。國(guó)內(nèi)研究方面，竇怡彬等[18]采用流-固耦合算法對(duì)超聲速二維壁板顫振特性進(jìn)行了分析，發(fā)現(xiàn)流場(chǎng)內(nèi)激波和膨脹波的交替變化誘發(fā)了壁板的極限環(huán)顫振。梅冠華等[19-20]發(fā)展了一種基于有限元方法的流-固耦合算法，并用其詳細(xì)分析了超聲速和跨聲速二維壁板氣動(dòng)彈性問(wèn)題，初步分析了流-固耦合結(jié)果與簡(jiǎn)化氣動(dòng)理論所得結(jié)果間的差異及其形成原因，其結(jié)果與Davis的流-固耦合結(jié)果吻合較好。安效民等[21]采用流-固耦合算法對(duì)跨/超聲速二維壁板顫振進(jìn)行了詳細(xì)研究，并對(duì)超聲速三維壁板顫振進(jìn)行了探究，與他人結(jié)果取得了很好的一致。

目前，在壁板顫振的流-固耦合研究中關(guān)注的基本都是平面壁板，事實(shí)上，受幾何外形影響，飛行器表面的壁板大多呈現(xiàn)為一定程度的彎曲，而并非是理想的平面壁板。從理論上說(shuō)，該彎曲可以看作是對(duì)平面狀態(tài)的一種擾動(dòng)，而該擾動(dòng)的引入往往會(huì)給系統(tǒng)的特性帶來(lái)顯著影響。因此，采用流-固耦合算法研究壁板彎曲對(duì)其氣動(dòng)彈性特性的影響，這在理論上和實(shí)際上都具有重要意義。

為研究初始彎曲對(duì)壁板氣動(dòng)彈性特性的影響，也作為先前研究工作的進(jìn)一步深入，對(duì)文獻(xiàn)[19-20]中基于有限元方法的流-固耦合算法適當(dāng)改動(dòng)，在時(shí)域內(nèi)對(duì)曲壁板顫振問(wèn)題進(jìn)行了分析。首先，采用Von Kármán幾何大變形理論給出了曲壁板運(yùn)動(dòng)方程，并用標(biāo)準(zhǔn)有限元方法對(duì)該方程進(jìn)行離散。然后，簡(jiǎn)述了流動(dòng)控制方程、求解方法、動(dòng)網(wǎng)格和流-固耦合方式。最后，采用所發(fā)展的流-固耦合算法對(duì)超聲速和跨聲速氣流下曲壁板的氣動(dòng)彈性響應(yīng)行了分析，著重考察了無(wú)量綱動(dòng)壓、相對(duì)彎曲高度、馬赫數(shù)、厚度比對(duì)系統(tǒng)特性的影響。

1 曲壁板的運(yùn)動(dòng)方程及其數(shù)值解法

二維曲壁板顫振模型如圖1所示：曲壁板置于剛性平面上，其上表面處于沿x方向的高速氣流中，下表面對(duì)應(yīng)著空腔。流體的來(lái)流速度、馬赫數(shù)和密度分別為U∞、Ma∞和ρ∞。平板的長(zhǎng)度、厚度、單位長(zhǎng)度質(zhì)量、彈性模量和泊松比分別為a、h、ρm、E和μ。平板上表面由流場(chǎng)所施加的氣動(dòng)載荷為p，下表面空腔壓力與來(lái)流壓力相同，皆為p∞。

圖1 二維曲壁板顫振示意圖Fig.1 Schematic of two-dimensional curved panel flutter

假設(shè)曲壁板的曲率恒定，應(yīng)變-位移關(guān)系采用Von Kármán幾何大變形理論表示，則其控制方程為：

(1)

式中:板的彎曲剛度D=Eh3/[12(1-μ2)]，大變形產(chǎn)生的中面拉伸載荷Nx為：

(2)

式中:氣動(dòng)載荷Δp為平板的上下表面壓力差，即：Δp=p-p∞。式(1)和式(2)中的Rx為恒定的曲率半徑，對(duì)最大彎曲高度H?a的淺殼，可用拋物線方程近似描述其初始形狀，即：

(3)

則曲率半徑為Rx=a2/8H。壁板兩端為簡(jiǎn)支邊界條件，滿足w=0，?2w/?x2=0。

將壁板沿長(zhǎng)度方向均勻劃分為N個(gè)單元，則單元長(zhǎng)度l=a/N。對(duì)其中的任意單元，設(shè)其兩端整體坐標(biāo)分別為x1和x2，引入坐標(biāo)變換ξ=(x-x1)/l，將整體坐標(biāo)[x1,x2]轉(zhuǎn)化為局部坐標(biāo)[0,1]。在局部坐標(biāo)系下，采用Hermite單元，則單元上的位移w可以表達(dá)為w=Neqe，其中qe={w1θ1w2θ2}T為單元位移向量，Ne=[N1N2N3N4]為單元形函數(shù)向量，各形函數(shù)的具體表達(dá)式為:

N1=1-3ξ2+2ξ3，N2=l(ξ-2ξ2+ξ3)，

N3=3ξ2-2ξ3，N4=l(ξ3-ξ2)。

采用標(biāo)準(zhǔn)Galerkin有限元方法對(duì)控制式(1)進(jìn)行變分以及空間離散，可將其轉(zhuǎn)化為關(guān)于時(shí)間二階導(dǎo)數(shù)的常微分方程組：

(4)

式中:

(5)

式(4)和式(5)中的整體矩陣由單元矩陣組裝而成，即：

各單元矩陣的具體表達(dá)式為:

而q={w1,θ1,w2,θ2,…,wN+1,θN+1}T為壁板整體位移。

2 流動(dòng)控制方程及其數(shù)值解法

高速流動(dòng)采用Euler方程描述，其數(shù)值解法為雙時(shí)間步長(zhǎng)推進(jìn)的特征線分裂有限元方法。

壁板上方流場(chǎng)計(jì)算區(qū)域及網(wǎng)格剖分情況分別如圖2和圖3所示，所計(jì)算的流動(dòng)區(qū)域呈半圓形，其半徑為25倍的壁板弦長(zhǎng)。圓弧上為遠(yuǎn)場(chǎng)邊界條件，壁板上游和下游皆為剛性壁面，在無(wú)黏流動(dòng)中將其法向速度設(shè)為0即可，流-固交界面受壁板運(yùn)動(dòng)的影響，將其給定為(u-us)·n=0，其中u為流體速度，us為壁板運(yùn)動(dòng)速度，n為壁板節(jié)點(diǎn)的法向量。激波捕捉方法為基于壓力二階導(dǎo)數(shù)的方法。

采用分塊結(jié)構(gòu)化網(wǎng)格剖分計(jì)算區(qū)域，子域1的上邊界為與彎曲壁板相對(duì)應(yīng)的彎曲邊界，這樣便于生成高質(zhì)量的網(wǎng)格，流場(chǎng)網(wǎng)格共計(jì)5 881個(gè)節(jié)點(diǎn)和11 520個(gè)三角形單元。

該部分的詳細(xì)敘述可參閱文獻(xiàn)[19-20]。

3 動(dòng)網(wǎng)格及流-固耦合方式

隨著流-固交界面的運(yùn)動(dòng)，流場(chǎng)內(nèi)部的計(jì)算網(wǎng)格也需進(jìn)行處理，由于壁板運(yùn)動(dòng)僅發(fā)生在y方向上，故僅需將壁板節(jié)點(diǎn)的位移在子域1內(nèi)沿著y方向均分即可，這樣既能保證網(wǎng)格質(zhì)量又不會(huì)在動(dòng)網(wǎng)格處理上耗費(fèi)大量的計(jì)算時(shí)間。

流-固耦合的方式為松耦合，在每個(gè)時(shí)間步流體傳遞氣動(dòng)載荷給壁板，壁板則傳遞位移和速度給流體，從而實(shí)現(xiàn)流動(dòng)和結(jié)構(gòu)間的雙向耦合。關(guān)于動(dòng)網(wǎng)格和流-固耦合方式的更詳細(xì)描述可參考文獻(xiàn)[19-20]。

圖2 流場(chǎng)計(jì)算區(qū)域(非等比例繪制)Fig.2 Computational domain for fluid field

圖3 流場(chǎng)網(wǎng)格Fig.3 Mesh for flow field

4 數(shù)值模擬結(jié)果和分析

采用FORTRAN語(yǔ)言實(shí)現(xiàn)上述算法，算法和代碼的驗(yàn)證工作已經(jīng)在文獻(xiàn)[19-20]中完成，故不再贅述。

4.1 超聲速氣流下的壁板響應(yīng)

若H/h=0，則曲壁板退化為平壁板。動(dòng)壓λ較小時(shí)，壁板將保持平面穩(wěn)定狀態(tài)，當(dāng)動(dòng)壓λ逐漸增大并越過(guò)臨界值λcr，系統(tǒng)將發(fā)生Hopf分岔，由平面穩(wěn)定狀態(tài)轉(zhuǎn)化為極限環(huán)振動(dòng)狀態(tài)，且極限環(huán)振幅隨λ的增大而增大，顫振近乎正負(fù)對(duì)稱。

與平壁板不同，受壁板彎曲的影響，H/h>0的曲壁板上下表面的初始?jí)毫Σ钪挡粸榱?，即其初始?xì)鈩?dòng)載荷非零，因此曲壁板將偏離其初始形狀，尋求新的靜態(tài)氣動(dòng)彈性平衡位置。以H/h=1和λ=200為例，系統(tǒng)的位移時(shí)間歷程如圖4所示，在初始?xì)鈩?dòng)載荷作用下，歷經(jīng)瞬態(tài)響應(yīng)的壁板達(dá)到了新的穩(wěn)態(tài)變形位置，并在該位置上和氣動(dòng)載荷達(dá)到了靜態(tài)平衡狀態(tài)。動(dòng)壓λ對(duì)壁板穩(wěn)態(tài)變形和氣動(dòng)載荷的影響如圖5所示：壁板變形以向下凹陷為主，僅后部略微向上拱起，隨著λ的增大最大彎曲位置逐漸后移且彎曲高度逐漸減??；在壁板前緣和尾緣附近存在兩個(gè)明顯的激波，前緣激波的強(qiáng)度隨著λ的增大而減弱，而后緣激波幾乎未受影響，壓力的整體分布則隨著λ的增大而趨向均勻。即，壁板變形和氣動(dòng)力之間的相互作用減弱了壁板的彎曲程度，也使得氣動(dòng)載荷的分布更加趨于均勻。

圖4 曲壁板的靜態(tài)氣動(dòng)彈性變形歷程(H/h=1，λ=200)Fig.4 Time history of static aeroelastic deformation for a curved panel (H/h=1, λ=200)

圖5 動(dòng)壓對(duì)曲壁板靜氣動(dòng)彈性特性的影響 (H/h=1)Fig.5 Effect of dynamic pressure on static aeroelastic behavior of curved panels (H/h=1)

當(dāng)動(dòng)壓λ進(jìn)一步增大并跨過(guò)臨界值λcr=248，系統(tǒng)發(fā)生了Hopf分岔，表現(xiàn)為極限環(huán)顫振。圖6給出了在λ=300下系統(tǒng)的響應(yīng)，經(jīng)過(guò)大約10個(gè)周期的瞬態(tài)振蕩壁板進(jìn)入到了穩(wěn)定的顫振狀態(tài)，該振動(dòng)并不是關(guān)于中性面正負(fù)對(duì)稱的，負(fù)向峰值的絕對(duì)值大約為正向峰值的兩倍，而平壁板的振動(dòng)則近乎是上下對(duì)稱的。這是因?yàn)榍诎迨窃谙掳嫉撵o態(tài)氣動(dòng)彈性變形位置上失穩(wěn)發(fā)生顫振的，振動(dòng)的中心應(yīng)在該穩(wěn)態(tài)變形位置附近，因此振動(dòng)呈現(xiàn)出了負(fù)向峰值較大的特性。動(dòng)壓λ對(duì)顫振峰值的影響如圖7所示，顫振負(fù)向峰值的絕對(duì)值較大，并隨著λ的增大而增大，顫振的正向峰值較小，且?guī)缀醣３衷谝粋€(gè)穩(wěn)定值上，圖7中還用虛線標(biāo)出了顫振的中點(diǎn)，可以清晰地發(fā)現(xiàn)其位于壁板失穩(wěn)前靜態(tài)氣動(dòng)彈性變形隨λ變化曲線的延長(zhǎng)線上，即恰是由于此平衡位置的失穩(wěn)才導(dǎo)致了曲壁板的非對(duì)稱型顫振。

圖6 曲壁板顫振的位移時(shí)間歷程 (H/h=1，λ=300)Fig.6 Time history of displacement for flutter of a curved panel (H/h=1, λ=300)

圖7 動(dòng)壓對(duì)曲壁板顫振正負(fù)峰值的影響 (H/h=1)Fig.7 Effect of dynamic pressure on positive and negative peak amplitudes of curved panel flutter (H/h=1)

進(jìn)一步，增大壁板的彎曲高度至H/h=2，其穩(wěn)態(tài)變形特性與H/h=1的曲壁板相同，不再重復(fù)。隨著動(dòng)壓的增大并越過(guò)臨界值λcr=196，壁板突然進(jìn)入到了混沌運(yùn)動(dòng)狀態(tài)，圖8(a)給出了λ=200時(shí)壁板的穩(wěn)態(tài)響應(yīng)，其表現(xiàn)為看似雜亂無(wú)章實(shí)則又有一定規(guī)律的運(yùn)動(dòng)狀態(tài)，同樣地，壁板的負(fù)向運(yùn)動(dòng)要比正向運(yùn)動(dòng)更為劇烈。將無(wú)量綱時(shí)間步長(zhǎng)變換為真實(shí)物理時(shí)間步長(zhǎng)，并對(duì)壁板的穩(wěn)態(tài)響應(yīng)進(jìn)行頻譜分析，所得結(jié)果如圖8(b)所示，可見該振動(dòng)包含了眾多的頻率組分。作為定量判定混沌的指標(biāo)，Lyapunov指數(shù)描述了系統(tǒng)臨近軌道間的距離在長(zhǎng)時(shí)間演變過(guò)程中指數(shù)發(fā)散程度的強(qiáng)弱。采用C-C算法計(jì)算了最佳延遲時(shí)間和重構(gòu)維數(shù)，對(duì)壁板穩(wěn)態(tài)響應(yīng)時(shí)間序列進(jìn)行了相空間重構(gòu)，并在此相空間上使用Wolf所提出的算法[22]得到了系統(tǒng)的最大Lyapunov指數(shù)，為0.088大于0，故可認(rèn)為系統(tǒng)處于混沌運(yùn)動(dòng)狀態(tài)。隨著λ的繼續(xù)增大，系統(tǒng)又出現(xiàn)了穩(wěn)態(tài)的周期振蕩，如圖9所示的λ=400時(shí)壁板的穩(wěn)態(tài)響應(yīng)。可見初始彎曲作為非線性擾動(dòng)對(duì)系統(tǒng)的特性可造成顯著影響，隨彎曲高度的增大壁板更加復(fù)雜的運(yùn)動(dòng)形式被不斷激發(fā)出來(lái)。

圖8 混沌運(yùn)動(dòng)狀態(tài)的曲壁板 (H/h=2，λ=200)Fig.8 Chaotic motion of a curved panel (H/h=2, λ=200)

圖9 曲壁板周期振動(dòng)的位移時(shí)間歷程 (H/h=2，λ=400)Fig.9 Time history of displacementfor a periodic oscillation of curved panel (H/h=2, λ=400)

4.2 跨聲速氣流下的壁板響應(yīng)

對(duì)于H/h=0的平壁板，在Ma∞=0.8或0.9時(shí)，隨h/a的減小壁板由平面穩(wěn)定狀態(tài)轉(zhuǎn)化為屈曲狀態(tài)，依據(jù)不同的初始條件到達(dá)正向或是負(fù)向的穩(wěn)定變形位置。變形量隨著h/a的減小而增大，這是因?yàn)楸诎逶奖。鋭偠仍叫?，越容易發(fā)生變形。Ma∞=0.8時(shí)壁板正負(fù)屈曲近乎對(duì)稱，而Ma∞=0.9時(shí)正向變形要略大于負(fù)向變形。

對(duì)于H/h>0的曲壁板，以H/h=1和h/a=0.002為例，在Ma∞=0.8時(shí)壁板的穩(wěn)態(tài)變形和壓力系數(shù)分布如圖10所示，由于壁板上下表面的壓力差值為兩端高中間低的對(duì)稱分布形式，且以負(fù)值為主，故壁板表現(xiàn)為左右對(duì)稱的上凸變形。接下來(lái)，考察了相對(duì)厚度h/a和相對(duì)彎曲高度H/h對(duì)壁板穩(wěn)態(tài)變形的影響，并與平壁板的正向屈曲結(jié)果進(jìn)行了對(duì)比，將其匯總在圖11中。與平壁板相比，曲壁板僅會(huì)出現(xiàn)上凸的穩(wěn)態(tài)變形，而并不會(huì)發(fā)生下凹的穩(wěn)態(tài)變形，這是因?yàn)槌跏細(xì)鈩?dòng)載荷以負(fù)值為主，迫使壁板運(yùn)動(dòng)到了向上的穩(wěn)態(tài)變形位置。另外，曲壁板的變形量要小于平壁板，且穩(wěn)態(tài)變形隨著彎曲高度的增大而減小，這說(shuō)明曲壁板抵抗穩(wěn)態(tài)變形的能力要強(qiáng)于平壁板。由于氣動(dòng)載荷隨著Ma∞的增加而增強(qiáng)，故Ma∞=0.9下的變形要大于Ma∞=0.8下的變形。

圖10 曲壁板穩(wěn)態(tài)變形和氣動(dòng)載荷分布 (海平面鋁板，Ma∞=0.8，h/a=0.002，H/h=1)Fig.10Stable deformation and aerodynamic load of a curved panel (aluminum panel at sea level, Ma∞=0.8, h/a=0.002, H/h=1)

圖11 h/a和H/h對(duì)曲壁板穩(wěn)態(tài)變形的影響 (海平面鋁板)Fig.11 Effect of h/a and H/h on stable deformation of curved panels (aluminum panel at sea level)

圖12 曲壁板的靜氣動(dòng)彈性變形和壓力系數(shù) (Ma∞=1.2，H/h=1，λ*=40)Fig.12 Static aeroelastic deformation and pressure coefficient of a curved panel (Ma∞=1.2, H/h=1, λ*=40)

圖13 一個(gè)周期內(nèi)的壁板瞬態(tài)變形 (Ma∞=1.2，H/h=1，λ*=100)Fig.13 Transient deformation of panel in one period (Ma∞=1.2, H/h=1, λ*=100)

圖14 極限環(huán)振幅隨動(dòng)壓變化 (Ma∞=1.2，H/h=1)Fig.14 Effect of dynamic pressure on limit cycle oscillation amplitude (Ma∞=1.2, H/h=1)

5 結(jié) 論

采用基于特征線分裂有限元方法的流-固耦合算法，分析了超聲速和跨聲速氣流作用下二維曲壁板的氣動(dòng)彈性特性。著重考察了馬赫數(shù)、相對(duì)彎曲高度和動(dòng)壓對(duì)系統(tǒng)特性的影響，通過(guò)對(duì)數(shù)值模擬結(jié)果進(jìn)行分析，所得主要結(jié)論如下：

(1)與平壁板不同，由于壁板彎曲造成的初始?xì)鈩?dòng)載荷不為零，曲壁板將偏離其初始狀態(tài)，達(dá)到靜態(tài)氣動(dòng)彈性變形位置。

(2)在超聲速氣流下，曲壁板失穩(wěn)后出現(xiàn)了極限環(huán)顫振，且與平壁板近乎對(duì)稱的振動(dòng)不同，由于曲壁板是從靜態(tài)氣動(dòng)彈性平衡位置失穩(wěn)進(jìn)入顫振的，故其顫振的負(fù)向峰值的絕對(duì)值遠(yuǎn)大于正向峰值。隨著壁板彎曲高度的增加，其非線性特性增強(qiáng)，混沌形式的顫振被激發(fā)了出來(lái)。

(3)在跨聲速氣流下，當(dāng)馬赫數(shù)為0.8和0.9時(shí)，曲壁板僅可發(fā)生正向的靜態(tài)氣動(dòng)彈性變形，而不會(huì)像平壁板那樣出現(xiàn)正負(fù)屈曲變形，且變形量隨彎曲高度的增大而減小，隨馬赫數(shù)的增大而增大。當(dāng)馬赫數(shù)為1.2時(shí)，曲壁板失穩(wěn)后表現(xiàn)為極限環(huán)顫振，且與超聲速情況類似，曲壁板的負(fù)向振動(dòng)更為劇烈。

作為初步研究，流動(dòng)采用Euler方程描述，并未考慮氣動(dòng)加熱效應(yīng)，且著重考察了壁板響應(yīng)特性。未來(lái)將采用Navier-Stokes方程，引入湍流模型和壁板導(dǎo)熱方程，發(fā)展流-固-熱耦合求解算法，對(duì)該耦合系統(tǒng)內(nèi)的物質(zhì)輸運(yùn)、動(dòng)量遷移和能量傳遞特征進(jìn)行深入剖析。

[1] MCNAMARA J J, FRIEDMANN P P, POWELL K G, et al. Aeroelastic and aerothermoelastic behavior in hypersonic flow [J]. AIAA Journal, 2008, 46(10): 2591-2610.

[2] LAMORTE N, FRIEDMANN P P. Hypersonic aeroelastic and aerothermoelastic studies using computational fluid dynamics [J]. AIAA Journal, 2014, 52(9): 2062-2078.

[3] DOWELL E H. Nonlinear oscillations of a fluttering plate [J]. AIAA Journal, 1966, 4(7): 1267-1275.

[4] DOWELL E H. Nonlinear flutter of curved plates [J]. AIAA Journal, 1969, 7(3): 424-431.

[5] OLSON M D. Finite element approach to panel flutter [J]. AIAA Journal, 1967, 5(12): 226-227.

[6] OLSON M D. Some flutter solutions using finite element [J]. AIAA Journal, 1970, 8(4): 747-752.

[7] YE W L, DOWELL E H. Limit cycle oscillation of a fluttering cantilever plate [J]. AIAA Journal,1991,29(11): 1929-1936.

[8] 梅冠華, 張家忠, 席光. 基于時(shí)滯慣性流形的二維平面壁板非線性氣動(dòng)彈性分析 [J]. 振動(dòng)與沖擊,2012,31(10): 141-146. MEI Guanhua, ZHANG Jiazhong, XI Guang. Nonlinear aeroelastic analysis of a two-dimensional panel based on inertial manifolds with delay [J]. Journal of Vibration and Shock, 2012, 31(10): 141-146.

[9] MEI Guanhua, ZHANG Jiazhong, WANG Zhuopu. Numerical analysis of panel flutter on inertial manifolds with delay [J]. Journal of Computational and Nonlinear Dynamics-Transactions of ASME, 2013, 8(2): 021009.

[10] GRAY C E, MEI C. Large amplitude finite element flutter analysis of composite panels in hypersonic flow [J]. AIAA Journal, 1993, 31(6): 1090-1099.

[11] 李凱倫, 張家忠. 功能梯度材料薄板的熱氣動(dòng)彈性數(shù)值分析方法及特性研究 [J]. 宇航學(xué)報(bào), 2013, 34(9): 1177-1186. LI Kailun, ZHANG Jiazhong. Numerical analysis method and aerothermoelastic behaviors of temperature dependent functional graded panels [J]. Journal of Astronautics, 2013, 34(9): 1177-1186.

[12] ASHLEY H, ZARTARIAN G. Piston theory—a new aerodynamic tool for the aeroelastician [J]. Journal of the Aeronautical Science, 1956, 23(12): 1109-1118.

[13] DAVIS G A, BENDIKSEN O O. Unsteady transonic two-dimensional Euler solutions using finite elements [J]. AIAA Journal, 1993, 31(6): 1051-1059.

[14] DAVIS G A. Transonic aeroelasticity solutions using finite elements in an arbitrary Larangian-Eulerian formulation [D]. Los Angeles:University of California, 1994.

[15] GORDINER R E, FITHEN R. Coupling of a nonlinear finite element structural method with a Navier-Stokes solver [J]. Computers & Structures, 2003, 81(2): 75-89.

[16] HASHIMOTO A, AOYAMA T. Effects of turbulent boundary layer on panel flutter [J]. AIAA Journal, 2009, 47(12): 2785-2791.

[17] MUBLSTEIN L, GASPERS P A, RIDDLE D W. An experimental study of the influence of the turbulent boundary layer on panel flutter[R]. NASA Paper No. TN D-4486, 1968.

[18] 竇怡彬, 徐敏, 蔡天星, 等. 基于CFD/CSD耦合的二維壁板顫振特性研究 [J]. 工程力學(xué), 2011, 28(6): 176-188. DOU Yibin, XU Min, CAI Tianxing, et al. Investigation of a two-dimensional panel flutter based on CFD/CSD coupling method [J]. Engineering Mechanics, 2011, 28(6): 176-181.[19] 梅冠華, 楊樹華, 張家忠, 等. 用于跨/超聲速壁板顫振精確分析的流-固耦合有限元算法 [J]. 西安交通大學(xué)學(xué)報(bào), 2014, 48(1): 73-83. MEI Guanhua, YANG Shuhua, ZHANG Jiazhong,et al. A fluid-structure coupling algorithm based on finite element method for precise analysis of transonic and supersonic panel flutter [J]. Journal of Xi’an Jiaotong University, 2014, 48(1): 73-83.

[20] MEI Guanhua, ZHANG Jiazhong, XI Guang, et al. Analysis of supersonic and transonic panel flutter using a fluid-structure coupling algorithm[J]. Journal of Vibration and Acoustics-Transactions of ASME, 2014, 136(3): 031013.

[21] 安效民, 胥偉, 徐敏. 非線性壁板顫振分析 [J]. 航空學(xué)報(bào), 2015, 36(4): 1119-1127. AN Xiaomin, XU Wei, XU Min. Analysis of nonlinear panel flutter [J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(4): 1119-1127.

[22] WOLF A, SWIFT J B, SWINNEY H L, et al. Determining Lyapunov exponents from a time series [J]. Physica D: Nonlinear Phenomena, 1985, 16(3): 285-317.

Aeroelastic analysis of curved panels in transonic/supersonic airflowbased on a fluid-structure coupling algorithm

MEI Guanhua1, ZHANG Jiazhong2, KANG Can1

(1. School of Energy and Power Engineering, Jiangsu University, Zhenjiang 212013, China; 2. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

A fluid-structure coupling algorithm was used to analyze the aeroelastic behaviors of curved panels in transonic/supersonic airflow. According to the Von Karman’s large deformation theory, the governing equation of curved panels was presented, and it was discretized by the standard finite element method. The governing equations of fluid, numerical method, moving mesh technique and fluid-structure coupling way were introduced briefly. Numerical simulations and analyses were carried out to study the aeroelastic behaviors of curved panels. The results demonstrate that the curvature causes nonzero initial aerodynamic load on panels, which brings about greatly different aeroelastic features for curved panels compared with flat panels. Static aeroelastic deformations exist on curved panels in steady state. At Mach number of 2, the asymmetric flutter will be born as curved panels lose their stability. As the curvature height increases, the chaotic flutter can be induced. At Mach number of 0.8 and 0.9, curved panels bear only positive aeroelastic deformation. At Mach number of 1.2, with the stability lost, curved panels flutter more violently in the negative direction. The results obtained could guide the panel design and flutter suppression for high performance flight vehicles. The presented algorithm could be extended to numerically analyze other aeroelastic problems.

panel flutter; fluid-structure coupling; aeroelasticity

國(guó)家973計(jì)劃(2012CB026002)；國(guó)家科技支撐計(jì)劃(2013BAF01B02)；江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目；江蘇大學(xué)高級(jí)人才科研啟動(dòng)基金(15JDG155)

2015-09-16 修改稿收到日期：2015-11-09

梅冠華 男，博士，1984年11月生

張家忠 男，博士，講師，教授，博士生導(dǎo)師，1968年1月生

O323

A

10.13465/j.cnki.jvs.2016.22.009