梅冠華,楊樹華,張家忠,孫旭,陳嘉輝

(1.西安交通大學(xué)能源與動(dòng)力工程學(xué)院, 710049, 西安; 2.沈陽鼓風(fēng)集團(tuán)股份有限公司, 110869, 沈陽)

用于跨/超聲速壁板顫振精確分析的流-固耦合有限元算法

梅冠華1,楊樹華2,張家忠1,孫旭1,陳嘉輝1

(1.西安交通大學(xué)能源與動(dòng)力工程學(xué)院, 710049, 西安; 2.沈陽鼓風(fēng)集團(tuán)股份有限公司, 110869, 沈陽)

為了精確和定量分析超聲速與跨聲速壁板的顫振特性,提出了一種基于有限元方法的流-固耦合算法,并用其研究了二維壁板顫振問題。首先,給出了壁板的von Kármán幾何大變形運(yùn)動(dòng)方程,以及高速氣流的歐拉控制方程。然后,采用標(biāo)準(zhǔn)有限元方法對(duì)壁板方程進(jìn)行空間離散,而對(duì)流動(dòng)控制方程的離散則運(yùn)用雙時(shí)間步長推進(jìn)的特征線分裂有限元方法,從而有效地消除了流場數(shù)值解的振蕩問題。隨后,采取松耦合算法實(shí)現(xiàn)了流體與固體間的數(shù)據(jù)傳遞。最后,運(yùn)用所提出的算法對(duì)超聲速和跨聲速氣流作用下壁板的氣動(dòng)彈性特性進(jìn)行了分析,考察了歸一化動(dòng)壓、預(yù)緊力和厚度比對(duì)系統(tǒng)特性的影響,并將該算法的分析結(jié)果與采用線性/非線性活塞理論和線性化勢流理論的經(jīng)典壁板顫振結(jié)果進(jìn)行了對(duì)比,證明該算法可以在較寬廣的馬赫數(shù)范圍內(nèi)給出氣動(dòng)力的精確描述,尤其適合于分析跨聲速氣流下的壁板氣動(dòng)彈性響應(yīng)。

壁板顫振;流-固耦合;特征線分裂算法;有限元方法;氣動(dòng)彈性

在特定參數(shù)下,飛行器表面暴露于高速氣流中的壁板結(jié)構(gòu)將產(chǎn)生自激振動(dòng)現(xiàn)象,常稱作壁板顫振。這類振動(dòng)是典型的氣動(dòng)彈性問題,往往是由慣性力、彈性力和氣動(dòng)力的共同作用而激發(fā)的,通常具有很強(qiáng)的非線性動(dòng)力學(xué)特性,并將嚴(yán)重影響飛行器的疲勞壽命、飛行性能、飛行安全和乘坐品質(zhì)。因此,深入研究壁板顫振問題對(duì)于高速飛行器的壁板設(shè)計(jì)、顫振抑制、疲勞壽命估計(jì)和顫振的合理利用都具有十分重要的意義[1-3]。

作為經(jīng)典問題,壁板顫振已由眾多學(xué)者進(jìn)行了深入細(xì)致的研究[4-10]。在大多數(shù)壁板顫振的研究中,壁板一般采用von Kármán幾何大變形理論描述,以便準(zhǔn)確捕捉系統(tǒng)的各類復(fù)雜響應(yīng),氣動(dòng)載荷大多采用簡化氣動(dòng)力模型近似表達(dá),如線性/非線性活塞理論[4,11]、線性化勢流理論[5-6]等,由此可推導(dǎo)出系統(tǒng)的偏微分控制方程。對(duì)于此類控制方程,常采用Galerkin方法、Rayleigh-Ritz方法、有限元方法(FEM)等將其離散為常微分控制方程,并運(yùn)用各種頻域或時(shí)域分析工具進(jìn)行研究。雖然這些簡化氣動(dòng)力理論使用起來較為方便,然而其適用范圍和精度都較為有限,一般而言,線性活塞理論適用于1.45時(shí),由于非線性效應(yīng)變得十分顯著,氣動(dòng)力常采用非線性活塞理論來描述。目前,尚沒有一個(gè)普遍適用于亞聲速、跨聲速、超聲速和高超聲速流動(dòng)的簡化氣動(dòng)力理論。此外,簡化氣動(dòng)力理論給出的氣動(dòng)載荷是與固體位移、速度及其偏導(dǎo)數(shù)相關(guān)的函數(shù),而并非是通過計(jì)算流體動(dòng)力學(xué)(CFD)所得出的,因此無法給出流場的細(xì)節(jié)信息,如黏性邊界層、旋渦、流動(dòng)分離、激波等,進(jìn)而無法準(zhǔn)確分析流體與固體間的耦合機(jī)理。事實(shí)上,這些正是壁板顫振定量分析中的開放性問題。

近年來,伴隨著計(jì)算機(jī)軟、硬件水平的飛速發(fā)展,流-固耦合算法被逐漸用于分析壁板顫振問題。其中,氣動(dòng)力由基于流體歐拉或Navier-Stokes控制方程的CFD方法獲得,同時(shí)將CFD與計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(CSD)耦合起來研究壁板的氣動(dòng)彈性響應(yīng)。流-固耦合算法不僅廣泛適用于亞聲速、跨聲速、超聲速和高超聲速流動(dòng),可以在幾乎所有馬赫數(shù)范圍內(nèi)給出氣動(dòng)力的精確表達(dá),而且還可以給出流場的細(xì)節(jié)描述,這些都是簡化氣動(dòng)力理論所不能實(shí)現(xiàn)的。使用流-固耦合算法,Davis求解了二維壁板顫振問題,發(fā)現(xiàn)了許多不同于簡化氣動(dòng)力理論的新現(xiàn)象[12-13];隨后,Gordiner研究了二維和三維壁板顫振問題,發(fā)現(xiàn)邊界層的存在推遲了顫振的發(fā)生[14];Hashimoto分析了跨聲速三維壁板顫振問題,在與他人實(shí)驗(yàn)結(jié)果取得一致的同時(shí),還發(fā)現(xiàn)湍流邊界層不僅對(duì)顫振可以起到穩(wěn)定作用,在特定條件下還能引起不穩(wěn)定效應(yīng)[15]。

在目前的壁板顫振的流-固耦合研究中,對(duì)流場的求解大多采用有限差分法(FDM)和有限體積法(FVM),而采用FEM的較少。雖然FEM的推導(dǎo)過程相對(duì)煩瑣,計(jì)算時(shí)間也稍長,然而FEM有其獨(dú)特的優(yōu)勢:①由于無須構(gòu)造交錯(cuò)網(wǎng)格,FEM的網(wǎng)格生成和數(shù)據(jù)結(jié)構(gòu)較為簡單;②在已經(jīng)生成的粗網(wǎng)格上,FEM可以通過選用高階單元來提高精度,而FDM和FVM則往往需要通過加密網(wǎng)格來提高精度;③由于FEM的普適性,可被廣泛用來處理多場耦合問題,尤其是流-固耦合問題。若能將流體和固體采用統(tǒng)一的FEM求解,就有可能實(shí)現(xiàn)流-固同步耦合求解,并在后續(xù)研究中引入系統(tǒng)自由度縮減方法[16],不僅能加快求解速度,更能提取出符合物理意義的氣動(dòng)彈性模態(tài)信息。

由于流體控制方程中的非線性對(duì)流項(xiàng)往往會(huì)引起數(shù)值解的振蕩,因此標(biāo)準(zhǔn)Galerkin FEM不能直接用于求解流體控制方程?；诖?眾多改進(jìn)的FEM被不斷提出,如流線迎風(fēng)Petrov-Galerkin(SUPG)方法、Taylor-Galerkin(TG)方法、Galerkin最小二乘(GLS)方法等,其中特征線分裂(CBS)方法由Zienkiewicz首先提出,并用于求解對(duì)流擴(kuò)散方程[17]。經(jīng)過眾多學(xué)者的不斷發(fā)展和完善,迄今為止,特征線分裂有限元方法(CBS-FEM)已被廣泛應(yīng)用于求解各類流動(dòng)問題,例如不可壓縮流動(dòng)、可壓縮流動(dòng)、含激波的高速流動(dòng)、湍流流動(dòng)等。近年來,將CBS-FEM發(fā)展應(yīng)用于求解動(dòng)邊界問題,如流-固耦合問題和自由液面流動(dòng)問題等,已成為新的研究熱點(diǎn)。

為了精確和定量分析超聲速與跨聲速壁板顫振特性,本文提出一種基于FEM的流-固耦合算法,并用來在時(shí)域內(nèi)對(duì)二維壁板顫振問題進(jìn)行分析。首先,分別采用von Kármán幾何大變形方程和歐拉方程來描述壁板變形和高速氣流。然后,對(duì)壁板控制方程采用標(biāo)準(zhǔn)FEM進(jìn)行空間離散,對(duì)流體控制方程則運(yùn)用雙時(shí)間步推進(jìn)的CBS-FEM進(jìn)行離散。隨后,采取松耦合算法實(shí)現(xiàn)流體與固體間的雙向耦合。最后,運(yùn)用所提出的算法對(duì)超聲速和跨聲速氣流下壁板的氣動(dòng)彈性響應(yīng)進(jìn)行分析,考察歸一化動(dòng)壓、預(yù)緊力和厚度比對(duì)系統(tǒng)特性的影響,并對(duì)壁板的穩(wěn)定性邊界進(jìn)行計(jì)算。

1 壁板顫振模型及其數(shù)值解法

若三維壁板的展向尺寸遠(yuǎn)大于其弦長,則可忽略展向效應(yīng),將其簡化為二維壁板,這給問題的分析帶來了很大便利。作為初步研究,先對(duì)二維壁板顫振問題進(jìn)行分析,其模型如圖1所示:二維平板置于剛性平面上,其上表面處于沿x方向的高速來流中,下表面對(duì)應(yīng)著空腔。流體的來流速度、馬赫數(shù)和密度分別為U∞、Ma∞和ρ∞,平板的長度、厚度、單位長度質(zhì)量、彈性模量和泊松比分別為a、h、ρm、E和μ,平板兩端初始拉伸面內(nèi)力為N0,上表面由流場所施加的氣動(dòng)載荷為p,下表面空腔壓力與來流壓力相同,皆為p∞。

圖1 二維平面壁板顫振示意圖

對(duì)于厚度遠(yuǎn)小于長度的壁板,即常見的薄板,可以采用Kirchhoff-Love假設(shè),僅考慮其橫向運(yùn)動(dòng)w(x,t)。應(yīng)變-位移關(guān)系采用von Kármán幾何大變形理論表示,則平板的控制方程為[4]

(1)

其中平板的彎曲剛度D=Eh3/[12(1-μ2)],大變形產(chǎn)生的中面拉伸載荷

(2)

氣動(dòng)載荷

Δp=p-p∞

(3)

壁板的邊界條件可以為簡支或固支。對(duì)于簡支端,其數(shù)學(xué)表達(dá)式為w=0,?2w/?x2=0;對(duì)于固支端,其數(shù)學(xué)表達(dá)式為w=0,θ=?w/?x=0。

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 FEM對(duì)控制方程(1)進(jìn)行變分以及空間離散,將其轉(zhuǎn)化為關(guān)于時(shí)間二階導(dǎo)數(shù)的常微分方程組

(4)

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

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

對(duì)于該動(dòng)力系統(tǒng),只要給定初始條件,采用自適應(yīng)Runge-Kutta積分方法數(shù)值求解即可。

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

在高速流動(dòng)問題中,黏性的影響相對(duì)較小,而可壓縮性的影響十分顯著,故為方便研究,可以采用可壓縮無黏Navier-Stokes方程,即歐拉方程來描述流場,其守恒形式為

連續(xù)方程

(5)

動(dòng)量方程

(6)

能量方程

(7)

補(bǔ)充方程

(8)

式中:ρ為流體密度;p為流體壓強(qiáng);e為單位質(zhì)量流體的內(nèi)能;Q為單位質(zhì)量流體的總能量;γ=1.4為空氣的比熱容比;xi、ui和Ui=ρui分別為坐標(biāo)、速度和守恒速度。在二維流動(dòng)中,下角標(biāo)i=1,2分別對(duì)應(yīng)著x和y方向。

流動(dòng)方程的求解采用基于特征線分裂算法的有限元方法,將流場變量沿著特征線運(yùn)動(dòng)方向進(jìn)行變換處理,從而可消除對(duì)流項(xiàng),并有效避免數(shù)值解的振蕩。當(dāng)前,在可壓縮流動(dòng)問題求解中廣泛采用的是單時(shí)間步長推進(jìn)的CBS-FEM,由于其大多采用完全顯式格式,使得時(shí)間步長的選取受限于收斂條件。此外,為方便計(jì)算,在求解過程中常進(jìn)行質(zhì)量集中處理,這種簡化對(duì)穩(wěn)態(tài)問題沒有任何影響,因?yàn)殡S著穩(wěn)態(tài)解的收斂,時(shí)間變化項(xiàng)也會(huì)隨之消失。然而,對(duì)于瞬態(tài)問題而言,該近似所帶來的誤差相當(dāng)嚴(yán)重,此時(shí)為獲得準(zhǔn)確解,可引入附加的迭代過程?；诖?本文采用雙時(shí)間步長CBS-FEM來求解歐拉方程。事實(shí)上這是隱式格式算法,在每個(gè)物理時(shí)間步的推進(jìn)過程中,通過引入一系列偽時(shí)間步的迭代過程,來獲得局部的偽定常收斂解[17-18]。這樣,真實(shí)時(shí)間步長的選取就不再受限于收斂條件,而且質(zhì)量集中帶來的誤差也可有效消除。雙時(shí)間步長CBS-FEM算法的具體流程如下。

第1步

(9)

第2步

(10)

第3步

(11)

第4步

(12)

修正速度、密度和能量

(ρQ)m+1=(ρQ)m+Δ(ρQ)

計(jì)算壓力

(13)

式(9)～式(13)中:上角標(biāo)n和n+1分別代表當(dāng)前時(shí)刻和下一時(shí)刻,m和m+1分別對(duì)應(yīng)前、后2個(gè)偽時(shí)刻;Δt為物理時(shí)間步長;Δτ為偽時(shí)間步長;θ1和θ2為松弛因子,一般選取θ1=1,θ2=0。

采用線性三角形單元剖分計(jì)算區(qū)域,則流體變量可表達(dá)為形函數(shù)和節(jié)點(diǎn)值相乘的形式

運(yùn)用標(biāo)準(zhǔn)Galerkin有限元方法,將式(9)～式(12)進(jìn)行變分和空間離散,并忽略高階項(xiàng),即可得到雙時(shí)間步長推進(jìn)CBS-FEM的離散格式。

一般而言,對(duì)于存在激波的流動(dòng)問題,如果直接進(jìn)行數(shù)值模擬,則解在激波附近會(huì)引起強(qiáng)烈的振蕩,影響收斂性。因此,必須添加合理的人工黏性項(xiàng),選用合適的激波捕捉方法來抑制振蕩并捕捉激波。本文采用基于壓力二階導(dǎo)數(shù)的激波捕捉方法[17,19]。

壁板上方流場計(jì)算區(qū)域及網(wǎng)格剖分如圖2和圖3所示,遠(yuǎn)場邊界位于25倍壁板弦長處。為了生成高質(zhì)量網(wǎng)格并獲得精確的氣動(dòng)載荷,將計(jì)算區(qū)域劃分為2個(gè)子區(qū)域:子域1為壁板上方0.5a高度的矩形區(qū)域,其余為子域2。對(duì)于子域1,流-固耦合界面和高度方向分別均布48和24個(gè)網(wǎng)格;對(duì)子域2,布置環(huán)繞子域1的48層C型網(wǎng)格。流場網(wǎng)格共計(jì)5 881個(gè)節(jié)點(diǎn)和11 520個(gè)三角形單元。

圖2 壁板上方流場計(jì)算區(qū)域(非等比例繪制)

(a)整體網(wǎng)格

(b)局部網(wǎng)格

遠(yuǎn)場邊界條件由Riemann無反射邊界條件給定[12]。壁板上游和下游皆為剛性壁面,在無黏流動(dòng)中將其法向速度設(shè)為0即可。流-固交界面雖然也為壁面邊界條件,然而由于壁板運(yùn)動(dòng)的影響,將其給定為(u-us)·n=0,其中u為流體速度,us為壁板運(yùn)動(dòng)速度,n為壁板節(jié)點(diǎn)的法向量。

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

隨著壁板的上下運(yùn)動(dòng),流-固交界面也在不斷運(yùn)動(dòng)。在該類動(dòng)邊界流動(dòng)問題的計(jì)算過程中,如果規(guī)定動(dòng)邊界上的節(jié)點(diǎn)隨邊界同步運(yùn)動(dòng),將使得動(dòng)邊界附近的網(wǎng)格單元發(fā)生扭曲變形甚至重疊,導(dǎo)致無法計(jì)算流場;如果規(guī)定流場網(wǎng)格固定不變,不隨邊界運(yùn)動(dòng)而移動(dòng),就會(huì)給邊界位置的捕捉帶來很大不便。常用的處理方法是采用動(dòng)網(wǎng)格技術(shù),讓動(dòng)邊界上的網(wǎng)格節(jié)點(diǎn)隨著邊界同步運(yùn)動(dòng),求解區(qū)域內(nèi)部的節(jié)點(diǎn)也隨之相應(yīng)調(diào)整,從而不僅可以準(zhǔn)確地描述動(dòng)邊界位置,還可以保證流體網(wǎng)格的質(zhì)量。

由于該流-固耦合問題的特殊性——壁板僅在y方向產(chǎn)生位移,而且在壁板上方的子域1內(nèi)均布了24層結(jié)構(gòu)化流體網(wǎng)格,故僅需將壁板節(jié)點(diǎn)的位移在子域1內(nèi)沿著y方向均分即可,這樣既能保證網(wǎng)格質(zhì)量,又不會(huì)在動(dòng)網(wǎng)格上耗費(fèi)大量的計(jì)算時(shí)間。作為驗(yàn)證,令壁板以w=0.05asin(2πx/a)形式彎曲,移動(dòng)后的網(wǎng)格如圖4所示。整體和子域1內(nèi)網(wǎng)格的最大單元偏斜度分別為0.568 2和0.485 7,遠(yuǎn)小于0.9,說明雖然壁板發(fā)生了較大的變形,然而網(wǎng)格的質(zhì)量仍得以較好地保持,這證實(shí)了動(dòng)網(wǎng)格方法的可靠性。

圖4 動(dòng)網(wǎng)格方法驗(yàn)證

(14)

式中:Δxi為舊網(wǎng)格節(jié)點(diǎn)至新網(wǎng)格節(jié)點(diǎn)的位移。

流-固耦合的方式為松耦合,在每個(gè)時(shí)間步流體與固體相互交換一次數(shù)據(jù),也稱為單步耦合,具體實(shí)施方法如下:

(1)將氣動(dòng)載荷Δp傳遞給壁板;

(2)將壁板響應(yīng)推進(jìn)至下一時(shí)刻;

(4)移動(dòng)流體網(wǎng)格;

(5)更新流場初值;

(6)將流場變量推進(jìn)到下一時(shí)刻;

(7)繼續(xù)下一個(gè)時(shí)間步直到計(jì)算結(jié)束。

與緊耦合和同步耦合算法相比,松耦合使用起來較為方便,僅需對(duì)CFD和CSD的程序進(jìn)行適當(dāng)修改即可,同時(shí)它的計(jì)算速度較快,是一種應(yīng)用較廣的雙向耦合算法。然而,由于在一個(gè)時(shí)間步長內(nèi)流體和固體僅相互交換一次數(shù)據(jù),且該數(shù)據(jù)并非是同一時(shí)刻的,因此兩者間存在著時(shí)間滯后,如果時(shí)間步長選取不當(dāng),很容易造成誤差累積并導(dǎo)致解的發(fā)散,故需要慎重選擇時(shí)間步長。經(jīng)過不斷嘗試,最終將時(shí)間步長選取為ΔtU∞/a=0.01,這樣既能保證解的準(zhǔn)確性,又不至于消耗過長的計(jì)算時(shí)間。對(duì)于該時(shí)間步長的獨(dú)立性驗(yàn)證,將在下一節(jié)進(jìn)行闡述。

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

4.1 網(wǎng)格無關(guān)性分析和時(shí)間步長獨(dú)立性驗(yàn)證

先分析流場網(wǎng)格的無關(guān)性,計(jì)算Ma∞=2時(shí)壁板在固定變形位置y=0.1asin(πx/a)上的穩(wěn)態(tài)擾流問題。計(jì)算3種不同疏密的網(wǎng)格,即第2節(jié)給出的標(biāo)準(zhǔn)網(wǎng)格、疏網(wǎng)格、密網(wǎng)格,在其壁板邊界分別均布48、24和72個(gè)單元,得到壁面上的相對(duì)壓力分布情況,如圖5所示,可見在壁板前緣和后緣存在著2個(gè)明顯的激波。標(biāo)準(zhǔn)網(wǎng)格與密網(wǎng)格所得出的結(jié)果并無明顯差異,說明采用標(biāo)準(zhǔn)網(wǎng)格可以得到足夠精度的流場。

圖5 流場網(wǎng)格無關(guān)性分析

圖6 固體網(wǎng)格無關(guān)性分析

圖7 流-固耦合算法中時(shí)間步長的獨(dú)立性分析

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

當(dāng)歸一化預(yù)緊力R=0時(shí),在較小的歸一化動(dòng)壓λ下,初始擾動(dòng)系統(tǒng)經(jīng)過瞬態(tài)振動(dòng)將恢復(fù)到初始位置,表現(xiàn)為平面穩(wěn)定狀態(tài)。如圖8所示λ=300時(shí)的壁板響應(yīng),經(jīng)過大約30至40個(gè)周期的瞬態(tài)振動(dòng),初始擾動(dòng)被完全衰減,系統(tǒng)恢復(fù)到了平面穩(wěn)定狀態(tài)。當(dāng)λ逐漸增大并越過臨界值λcr,系統(tǒng)將發(fā)生Hopf分岔,由平面穩(wěn)定狀態(tài)轉(zhuǎn)變?yōu)闃O限環(huán)顫振狀態(tài)。λ=500時(shí)的系統(tǒng)響應(yīng)如圖9a所示,穩(wěn)態(tài)極限環(huán)顫振的最大值和最小值分別為0.660和-0.679,后者的絕對(duì)值略大于前者,表明該振動(dòng)關(guān)于中性面是非對(duì)稱的,這意味著壁板的正向運(yùn)動(dòng)和負(fù)向運(yùn)動(dòng)對(duì)于流場的影響是不同的。

作為該振動(dòng)的詳細(xì)分析,圖9b和圖9c分別給出了在一個(gè)周期內(nèi)不同時(shí)刻的壁板變形和氣動(dòng)載荷。從圖9b中可以看出,雖然不同空間質(zhì)點(diǎn)的振幅不同,然而各質(zhì)點(diǎn)在時(shí)間上的運(yùn)動(dòng)規(guī)律是相同的,即它們同時(shí)達(dá)到最大或最小位移,這是駐波的典型特征,因此可認(rèn)為壁板振動(dòng)呈現(xiàn)出了駐波的形式,其中最大振幅位于x=0.75a附近,最小振幅位于x=0.3a附近。圖9c所示的壓力波動(dòng)也呈現(xiàn)出駐波的運(yùn)動(dòng)形式,壓力最大值約位于壁板末端x=a附近,且隨著壁板的上下運(yùn)動(dòng)在壁板尾緣附近表現(xiàn)為激波與膨脹波的交替變換。由于激波是流場中典型的非線性現(xiàn)象,在激波附近的流體被強(qiáng)烈壓縮,因此激波的存在對(duì)壁板的顫振具有重大影響。

(a)位移時(shí)間歷程

(b)一個(gè)周期內(nèi)的瞬時(shí)壁板位移

(c)一個(gè)周期內(nèi)的瞬時(shí)壁板表面壓力分布

為考察系統(tǒng)對(duì)歸一化動(dòng)壓的敏感性,研究了λ對(duì)極限環(huán)顫振的影響,結(jié)果如圖10所示。由圖中可以看出,在λ>λcr時(shí)壁板發(fā)生顫振,振幅與壁板厚度量級(jí)相同,且隨λ的增大而增大。圖10中還展示了線性和非線性活塞理論所得的結(jié)果[21]:線性活塞理論所得出的振動(dòng)正負(fù)對(duì)稱,而非線性活塞理論所得振動(dòng)的負(fù)向峰值絕對(duì)值大于正向峰值,這與本文算法的結(jié)果是一致的。定量來說,本文算法所得出的臨界動(dòng)壓稍大,而振幅稍小。由于活塞理論沒有考察壁板不同時(shí)刻和不同空間質(zhì)點(diǎn)間的相互影響,是一種基于當(dāng)?shù)亓鲃?dòng)和當(dāng)前時(shí)刻的簡化氣動(dòng)力理論,而基于歐拉方程的流-固耦合算法則直接從流體控制方程著手,綜合考慮了上述因素,利用CFD方法可以精確計(jì)算氣動(dòng)載荷,因此這是兩者間存在差異的主要原因。

(a)位移時(shí)間歷程

(b)屈曲變形和壓力分布

在施加了預(yù)緊力的情況下,壁板更多形式的響應(yīng)將被激發(fā)出來。在較小的λ下,當(dāng)歸一化預(yù)緊力R逐漸增大并跨過臨界值時(shí),系統(tǒng)將發(fā)生靜態(tài)分岔,從平面穩(wěn)定狀態(tài)轉(zhuǎn)變?yōu)榍鸂顟B(tài)。例如:當(dāng)λ=50和R=2時(shí),壁板的位移時(shí)間歷程如圖11a所示,在動(dòng)壓和預(yù)緊力的聯(lián)合作用下,系統(tǒng)最終到達(dá)穩(wěn)定屈曲變形位置;壁板的最終變形和表面壓力分布如圖11b所示,最大變形位置大約位于x=0.6a附近,壁板的前緣和尾緣存在2個(gè)明顯的激波,其中尾緣處的激波比前緣處的激波更為劇烈,這可能是屈曲最大變形位置略向后偏移的結(jié)果。

(a)位移時(shí)間歷程

(b)相圖

(c)頻譜分析

在特定的λ和R下,系統(tǒng)可表現(xiàn)出非線性復(fù)雜振動(dòng)。例如:在λ=140和R=5時(shí),壁板的穩(wěn)態(tài)響應(yīng)如圖12a所示,雖然振動(dòng)為周期形式,然而振動(dòng)方式卻表現(xiàn)出了極強(qiáng)的非線性特性;在圖12b所示的相圖中,不難發(fā)現(xiàn)該振動(dòng)具有倍周期運(yùn)動(dòng)的特性;在圖12c所示的該振動(dòng)的頻譜中包含了6階頻率分量,分別約為15、30、45、60、75和90 Hz,其中前3階頻率所占比重較大,后3階頻率的比重較小,由于頻率間呈現(xiàn)出1至6的整數(shù)倍關(guān)系,因此系統(tǒng)表現(xiàn)出了倍周期運(yùn)動(dòng)的特征。

在R和λ所構(gòu)成的參數(shù)平面上,對(duì)系統(tǒng)的平面穩(wěn)定、屈曲、振動(dòng)響應(yīng)區(qū)域進(jìn)行了劃分,如圖13所示:隨著R從0開始逐漸增大,壁板的顫振臨界動(dòng)壓λcr逐漸降低,即預(yù)緊力的增加使得壁板更易發(fā)生顫振;隨著λ從0開始逐漸增大,壁板的屈曲臨界預(yù)緊力R逐漸增大,即氣動(dòng)載荷的增加使得壁板不易發(fā)生屈曲。此外,本文算法所得結(jié)果與線性活塞理論所得結(jié)果是一致的。

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

歸一化動(dòng)壓λ*對(duì)于顫振幅值的影響如圖15所示:當(dāng)λ*越過臨界值后,最大振幅隨著λ*的增加而增大,本文算法所得結(jié)果與文獻(xiàn)[5]采用線性化勢流理論所得結(jié)果吻合較好,然而隨著λ*的增大,兩者間的差距逐漸增大,這可能是由于線性化勢流理論是從擾動(dòng)速度和勢流理論的角度來推導(dǎo)氣動(dòng)力的,并不能準(zhǔn)確反映整體流場的狀況,與歐拉方程算得的氣動(dòng)力仍存在一定差距的緣故。

對(duì)于Ma∞<1(如Ma∞=0.8和Ma∞=0.9)的跨聲速氣流,系統(tǒng)失穩(wěn)后壁板將偏離初始位置,根據(jù)初始條件的不同到達(dá)正向或是負(fù)向屈曲變形位置。以位于海平面的鋁板為例,當(dāng)Ma∞=0.9和h/a=0.004時(shí),壁板的屈曲響應(yīng)如圖16所示,可以看到其正向屈曲變形略大于負(fù)向屈曲變形。

關(guān)于系統(tǒng)對(duì)厚度比的敏感性,圖17給出了在Ma∞=0.8和Ma∞=0.9時(shí)h/a對(duì)壁板正向和負(fù)向屈曲變形的影響:最大屈曲變形隨著h/a的增大而減小,這是因?yàn)檩^小的h/a意味著壁板相對(duì)剛度較小,更容易發(fā)生屈曲;Ma∞=0.8時(shí)正、負(fù)屈曲變形近乎對(duì)稱,而Ma∞=0.9時(shí)正向變形要大于負(fù)向變形,這說明Ma越接近1,系統(tǒng)變得越復(fù)雜。此外,本文算法所得結(jié)果與文獻(xiàn)[13]的結(jié)果相當(dāng)一致。

圖18顯示了Ma∞=0.8～1.4時(shí)跨聲速氣流作用下壁板的穩(wěn)定性邊界。在Ma∞=1.0的聲速氣流附近,臨界動(dòng)壓明顯降低,穩(wěn)定性曲線在此處形成了一個(gè)“凹坑”,通常稱為“跨聲速凹坑”,因此對(duì)壁板顫振的研究設(shè)計(jì)者來說,聲速附近的臨界動(dòng)壓是需要仔細(xì)考察的。本文算法所獲得的穩(wěn)定性邊界與文獻(xiàn)[13]的結(jié)果亦相當(dāng)一致。

(a)正向屈曲

(b)負(fù)向屈曲Ma∞=0.9, h/a=0.004, 海平面鋁板

(a)Ma∞=0.8

(b)Ma∞=0.9

圖18 壁板的穩(wěn)定性邊界

5 結(jié)論與展望

為了精確定量分析超聲速和跨聲速壁板的顫振特性,本文提出了一種基于有限元方法的流-固耦合算法,并用其研究了二維壁板顫振問題。該算法成功消除了流體方程中對(duì)流項(xiàng)所引發(fā)的數(shù)值振蕩,并實(shí)現(xiàn)了流體和固體間的雙向數(shù)據(jù)交換。通過對(duì)數(shù)值結(jié)果進(jìn)行分析,所得主要結(jié)論如下:

(1)在馬赫數(shù)為2的超聲速氣流下,當(dāng)歸一化動(dòng)壓越過臨界值后,系統(tǒng)由平面穩(wěn)定狀態(tài)轉(zhuǎn)變?yōu)闃O限環(huán)顫振狀態(tài),振幅隨動(dòng)壓增大而增大,且負(fù)向峰值的絕對(duì)值略大于正向峰值。對(duì)典型極限環(huán)顫振的分析發(fā)現(xiàn):一個(gè)周期內(nèi)壁板變形和氣動(dòng)載荷皆呈現(xiàn)駐波的運(yùn)動(dòng)形式;在歸一化預(yù)緊力與歸一化動(dòng)壓的共同作用下,系統(tǒng)可表現(xiàn)為屈曲或非簡諧復(fù)雜振動(dòng);激波對(duì)壁板顫振特性的影響較為顯著。

(2)在馬赫數(shù)為1.2的跨聲速氣流下,壁板響應(yīng)與馬赫數(shù)為2的系統(tǒng)類似,只是此時(shí)系統(tǒng)對(duì)外界擾動(dòng)反應(yīng)較快,達(dá)到穩(wěn)態(tài)響應(yīng)所需的瞬態(tài)過程較短。

(3)在馬赫數(shù)為0.8和0.9的亞聲速氣流下,壁板失穩(wěn)后表現(xiàn)為正向或負(fù)向的屈曲變形,正、負(fù)位置由初始條件所決定,且屈曲變形隨著壁板厚度比的減小而增大。

(4)對(duì)馬赫數(shù)在0.8至1.4之間的系統(tǒng)的穩(wěn)定性邊界進(jìn)行了計(jì)算,精確捕捉到了“跨聲速凹坑”現(xiàn)象。

(5)通過與采用線性/非線性活塞理論和線性化勢流理論的經(jīng)典壁板顫振結(jié)果進(jìn)行對(duì)比,表明本文算法可在更為寬廣的馬赫數(shù)范圍內(nèi)給出氣動(dòng)載荷的精確描述,尤其適合于分析跨聲速氣流下的壁板氣動(dòng)彈性特性。

作為初步探索,本文對(duì)流動(dòng)采用了歐拉方程描述,下一步將采用完全Navier-Stokes方程,并考察湍流的影響,以更加準(zhǔn)確地描述流場和氣動(dòng)載荷。此外,流-固耦合的計(jì)算過程耗費(fèi)了大量的計(jì)算時(shí)間,未來將引入有效的系統(tǒng)自由度縮減方法,在保證計(jì)算精度的同時(shí)實(shí)現(xiàn)節(jié)省計(jì)算時(shí)間的目的。本文的工作僅對(duì)固體的氣動(dòng)彈性特性進(jìn)行了分析,今后將對(duì)流動(dòng)的穩(wěn)定性與失穩(wěn)機(jī)理展開深入研究。

[1] 梅冠華, 張家忠. 時(shí)滯慣性流形在三維壁板顫振數(shù)值分析中的應(yīng)用 [J]. 西安交通大學(xué)學(xué)報(bào), 2011, 45(9): 40-46. MEI Guanhua, ZHANG Jiazhong. Numerical analysis of 3-D panel flutter by inertial manifolds with delay [J]. Journal of Xi’an Jiaotong University, 2011, 45(9): 40-46.

[2] 楊智春, 夏巍, 孫浩. 高速飛行器壁板顫振的分析模型和分析方法 [J]. 應(yīng)用力學(xué)學(xué)報(bào), 2006, 23(4): 537-542. YANG Zhichun, XIA Wei, SUN Hao. Analysis of panel flutter in high speed flight vehicles [J]. Chinese Journal of Applied Mechanics, 2006, 23(4): 537-542.

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

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

[5] DOWELL E H. Nonlinear oscillations of a fluttering plate: II [J]. AIAA Journal, 1967, 5(10): 1856-1862.

[6] DOWELL E H. A review of the aeroelastic stability of plates and shells [J]. AIAA Journal, 1970, 8(1): 385-399.

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

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

[9] MEI Chuh, ABDEL-MOTAGLY K, CHEN R. Review of nonlinear panel flutter at supersonic and hypersonic speeds [J]. ASME Applied Mechanics Reviews, 1999, 52(10): 321-332.

[10]CHENG Guangfeng, MEI Chuh. Finite element modal formulation for hypersonic panel flutter analysis with thermal effects [J]. AIAA Journal, 2004, 42(4): 687-695.

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

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

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

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

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

[16]ZHANG Jiazhong, REN Sheng, MEI Guanhua. Model reduction on inertial manifolds for N-S equations approached by multilevel finite element method [J]. Communications in Nonlinear Science and Numerical Simulation, 2011, 16(1): 195-205.

[17]ZIENKIEWICZ O C, TAYLOR R L, NITHIARASU P. The finite element method for fluid dynamics [M]. 6th ed. Singapore: Elsevier Pte Ltd., 2009: 195-211.

[18]NITHIARASU P. An efficient artificial compressibility (AC) scheme based on the characteristic based split (CBS) method for incompressible flows [J]. International Journal for Numerical Methods in Engineering, 2003, 56(13): 1815-1845.

[19]NITHIARASU P, ZIENKIEWICZ O C, SATYASAI B V K, et al. Shock capturing viscosities for the general fluid mechanics algorithm [J]. International Journal for Numerical Methods in Fluids, 1998, 28(9): 1325-1353.

[20]孫旭, 張家忠. 具有運(yùn)動(dòng)邊界的不可壓縮黏性流動(dòng)的CBS有限元解法 [J]. 西安交通大學(xué)學(xué)報(bào), 2011, 45(1): 99-104. SUN Xu, ZHANG Jiazhong. A characteristic-based split-FEM scheme for incompressible viscous flow with moving boundaries [J]. Journal of Xi’an Jiaotong University, 2011, 45(1): 99-104.

[21]EASTEP F E, MCINTOSH S C. Analysis of nonlinear panel flutter and response under random excitation or nonlinear aerodynamic loading [J]. AIAA Journal, 1971, 9(3): 411-418.

(編輯 葛趙青)

AFluid-StructureCouplingAlgorithmBasedonFiniteElementMethodforPreciseAnalysisofTransonicandSupersonicPanelFlutter

MEI Guanhua1,YANG Shuhua2,ZHANG Jiazhong1,SUN Xu1,CHEN Jiahui1

(1. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2. Shenyang Blower Works Group Corporation, Shenyang 110869, China)

To analyze the supersonic and transonic panel flutter behavior quantitatively and accurately, a fluid-structure coupling algorithm based on the finite element method (FEM) is proposed for the two-dimensional panel flutter problem. First, the von Kármán’s large deformation theory is adopted to model the panel, and the high speed air flow is approached by the Euler equations. Then, the equation of panel is discretized spatially by the standard FEM, and the equations of fluid are discretized by the characteristic-based split finite element method (CBS-FEM) with dual time stepping, thus the numerical oscillation often encountered in numerical simulation of fluid flow can be eliminated. Furthermore, a loose coupling algorithm is applied to the data exchange between the fluid and the structure. Finally, the proposed algorithm is used to investigate the aeroelastic behavior of the panel in supersonic and transonic air flows and the influences of the non-dimensional dynamic pressure, pre-tightening force and thickness ratio on the system. The results are compared with those of the classical panel flutter analyses using linear/nonlinear piston theory and linearized potential flow theory. It shows that the proposed algorithm enables to obtain accurate aerodynamic pressure in a wide range of Mach numbers, especially for the analysis of panel aeroelasticity in transonic air flows.

panel flutter; fluid-structure interaction; characteristic-based split method; finite element method; aeroelasticity

10.7652/xjtuxb201401013

2013-03-29。 作者簡介: 梅冠華(1984—),男,博士生;張家忠(通信作者),男,教授,博士生導(dǎo)師。 基金項(xiàng)目: 國家“973計(jì)劃”資助項(xiàng)目(2012CB026002);國家“863計(jì)劃”資助項(xiàng)目(2012AA052303)。

時(shí)間: 2013-10-22 網(wǎng)絡(luò)出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131022.1113.003.html

O323

:A

:0253-987X(2014)01-0073-11