代林桐，邢譽(yù)峰

(北京航空航天大學(xué) 固體力學(xué)所，北京 100083）

顫振是指彈性結(jié)構(gòu)在氣流中由于受到氣動(dòng)力、彈性力和慣性力的耦合作用而發(fā)生的振幅不衰減的自激振動(dòng)，對(duì)結(jié)構(gòu)的使用壽命影響很大，甚至危及結(jié)構(gòu)的安全［1］。

隨著飛行馬赫數(shù)的不斷提高，飛行器面臨的氣動(dòng)環(huán)境越來(lái)越嚴(yán)峻。熱氣動(dòng)彈性是近年來(lái)最熱門(mén)的方向之一，研究成果不斷涌現(xiàn)。作為有著優(yōu)良的力學(xué)性能和耐熱特性的功能梯度材料(FGM）板，研究其氣彈及熱氣動(dòng)彈性問(wèn)題具有重要的理論和應(yīng)用價(jià)值。這也給氣動(dòng)彈性問(wèn)題的研究帶來(lái)了新的挑戰(zhàn)，那就是在分析中必須考慮高馬赫數(shù)下氣動(dòng)加熱所引起的熱效應(yīng)對(duì)飛行器結(jié)構(gòu)氣動(dòng)彈性的影響［2］。

Praveen和Reddy［3］對(duì)功能梯度壁板進(jìn)行了非線性瞬態(tài)熱彈性分析，并討論了溫度場(chǎng)對(duì)壁板響應(yīng)的影響。苑凱華和邱志平［4］利用有限元法建立了壁板在熱效應(yīng)下的運(yùn)動(dòng)微分方程，基于非線性動(dòng)力學(xué)模型研究了對(duì)壁板施加控制后對(duì)顫振極限環(huán)幅值的影響。Hosseini等［5］對(duì)超聲速氣流中功能梯度平板的熱氣彈問(wèn)題進(jìn)行了深入分析，建立了基于馮·卡門(mén)薄板大變形理論的結(jié)構(gòu)方程并考慮了溫度變化對(duì)壁板材料性能的影響，通過(guò)一階活塞理論模擬氣動(dòng)力，并采用Galerkin方法對(duì)壁板的控制方程進(jìn)行求解。李麗麗和趙永輝［6］研究了熱環(huán)境下四邊固支壁板結(jié)構(gòu)頻率特性的變化，進(jìn)而利用p-k法進(jìn)行了顫振分析，研究表明，熱效應(yīng)對(duì)結(jié)構(gòu)的頻率特性有很大的影響，并影響顫振邊界。Shahverdi和Khalafi［7］利用廣義微分求積的方法研究了功能梯度曲板的熱氣動(dòng)彈性行為。黃小林等［8］基于復(fù)合材料薄板理論和有氣流偏角的氣動(dòng)壓力的一階活塞模型，用Galerkin方法分析了氣流偏角、熱環(huán)境等因素對(duì)FGM板固有頻率和顫振臨界速壓的影響。

壁板顫振分析最終歸結(jié)為利用不同的結(jié)構(gòu)和氣動(dòng)力模型求解顫振微分方程的問(wèn)題，其解法包括數(shù)值解法和解析解法。

已有文獻(xiàn)主要采用Galerkin方法［9］及Rayleigh-Ritz方法［10］分析壁板的非線性顫振特性。雖然這2種方法精度比較高，但僅能從數(shù)值解的角度得到顫振特性。而解析解法，如直接求解法［11］、半逆法［12］，這幾種方法都能直接從求解本征方程入手來(lái)研究顫振特性。

基于活塞理論，Li和Song［11］采用直接求解法求得了Kirchhoff和Mindlin板在不同邊界下的顫振邊界，但沒(méi)有考慮熱效應(yīng)的影響。Sun和Xing［13］根據(jù)經(jīng)典板理論和一階活塞理論建立了二維層合板在簡(jiǎn)支、固支和自由邊界條件下的氣動(dòng)彈性模型，得到了壁板顫振本征解的統(tǒng)一顯式形式，并與Galerkin方法進(jìn)行了比較，但也未考慮超聲速下氣動(dòng)熱的影響。

目前，鮮見(jiàn)到公開(kāi)發(fā)表的關(guān)于功能梯度壁板熱顫振本征問(wèn)題精確解的工作。針對(duì)該問(wèn)題，本文基于經(jīng)典薄板理論和一階活塞理論求得了不同邊界下壁板熱顫振本征問(wèn)題的精確解并分析了其顫振特性。

1 顫振微分方程

考慮如圖1所示的二維FGM壁板模型，其弦長(zhǎng)為a，展長(zhǎng)為無(wú)限長(zhǎng)，厚為h，NTx相當(dāng)于熱應(yīng)力的作用，T為溫度。坐標(biāo)系建立在板的中面，原點(diǎn)在板的中心處。壁板的上表面為陶瓷層，下表面為金屬層。

圖1 兩端簡(jiǎn)支壁板幾何模型Fig.1 Geometric model of simply supported panel at both ends

FGM板的彈性模量E、熱膨脹系數(shù)α、熱傳導(dǎo)系數(shù)K及密度ρ等物理屬性量皆按照式(1）進(jìn)行估計(jì)［14］：

式中：P(z）為板內(nèi)任意一點(diǎn)的材料參數(shù)(彈性模量、泊松比、剪切模量等）；Pt和Pb分別為上表面(用t表示）和下表面(用b表示）對(duì)應(yīng)的材料參數(shù)；n為材料梯度指數(shù)。

功能梯度材料通常工作在一些高溫環(huán)境，其不同屬性對(duì)溫度的敏感程度是不一樣的。設(shè)彈性模量E和熱膨脹系數(shù)α隨溫度的變化規(guī)律為［15］

式中：Pi(i＝－1，0，1，2，3）為材料某一物理屬性的溫度相關(guān)系數(shù)，其值是唯一的。

由式(2）可知，要獲得有效物理屬性，除了知道其對(duì)應(yīng)的溫度系數(shù)以外，還需要知道結(jié)構(gòu)溫度場(chǎng)分布。當(dāng)材料熱傳導(dǎo)系數(shù)為常數(shù)時(shí)，由一維熱傳導(dǎo)穩(wěn)態(tài)方程所得到的是一個(gè)線性的溫度場(chǎng)。對(duì)于功能梯度材料組成的板而言，當(dāng)熱傳導(dǎo)系數(shù)k(i）是厚度坐標(biāo)的函數(shù)而不再是一個(gè)常數(shù)時(shí)，這時(shí)得到一個(gè)非線性的溫度場(chǎng)。該非線性溫度場(chǎng)滿足如下一維熱傳導(dǎo)方程：

相關(guān)文獻(xiàn)給出了其近似理論解［16］。

假定壁板上表面受氣動(dòng)力作用，另一側(cè)為靜止空氣。氣體流速為V，密度為ρa(bǔ)，馬赫數(shù)為Ma。對(duì)于本文后續(xù)選擇的Si3N4和SUS304作為組成成分的FGM板，泊松比υ沿著厚度方向的變化對(duì)結(jié)果幾乎沒(méi)有影響，因此本文在后續(xù)計(jì)算中取泊松比為2種材料的泊松比的均值υ＝0.28［17］，稱為等效泊松比。

下面基于經(jīng)典薄板理論和Ham ilton變分原理，推導(dǎo)二維壁板顫振控制微分方程。在薄板理論中，位移分量可以表示為

式中：u0為FGM 矩形薄板幾何中面x方向的位移；w為FGM薄板的撓度。

在二維壁板顫振問(wèn)題分析中，展長(zhǎng)假設(shè)為無(wú)限長(zhǎng)，因此不考慮y方向位移的作用。在下面公式推導(dǎo)中，采用的是單位寬度。

FGM薄板的應(yīng)變可以表示為

應(yīng)力為

式中：εx(T）為溫度應(yīng)變，即

熱環(huán)境中FGM矩形薄板的總應(yīng)變能為

FGM薄板自由振動(dòng)的動(dòng)能為

對(duì)于功能梯度板來(lái)說(shuō)，由于材料不均勻，材料中性面與幾何中面并不重合，設(shè)二者之間距離為z0。在中性面上位移u＝0，從式(4）可得

由式(6）可得

軸力為

本文等效泊松比為常數(shù)，因此利用軸力等于0得到

實(shí)際上，即使泊松比不是常數(shù)，中性面也是存在的，這不同于三維壁板情況［18］。

根據(jù)一階活塞理論［19］，氣動(dòng)力的形式為

式中：q為動(dòng)壓。

氣動(dòng)力所做虛功δW的積分為

根據(jù)Ham ilton原理的廣義形式為

可以得到二維功能梯度壁板在超聲速流下的運(yùn)動(dòng)微分方程，如下：

其中：Deq為等效剛度；為熱應(yīng)力。

本文用到的簡(jiǎn)支和固支邊界條件如表1所示。

表1 簡(jiǎn)支和固支邊界條件Table 1 Boundary conditions for sim p le support and clam p

下面將根據(jù)邊界條件求方程(17）的本征精確解，并分析梯度指數(shù)、溫度場(chǎng)對(duì)顫振邊界的影響。

2 顫振精確本征解和顫振機(jī)理

2.1 精確本征解

設(shè)撓度w有如下分離變量的形式：

式中：顫振的穩(wěn)定性取決于本征值Ω，其實(shí)部β代表幅值，而虛部ω代表振動(dòng)頻率。若β＝0，壁板將發(fā)生顫振。

將式(19）代入控制微分方程(17）可得

式(20）即是二維壁板顫振的本征微分方程。根據(jù)邊界條件求解方程(20）可以得到顫振問(wèn)題的模態(tài)函數(shù)φ和頻率方程。引入無(wú)量綱量ξ＝x／a。為了求得本征微分方程的通解，φ可寫(xiě)為

式中：A為待定系數(shù)；λ為x方向的空間本征根。

將式(21）代入方程(20）可得

或?qū)懗?/p>

這即是二維壁板顫振本征代數(shù)方程，其中R、p和k都是無(wú)量綱系數(shù)，具體形式為

求解方程(23）可以得到二維壁板的空間本征根，為一對(duì)實(shí)數(shù)根和一對(duì)復(fù)數(shù)根：

則本征函數(shù)φ的通解可表示為

將本征根(25）代入本征代數(shù)方程(23），得到如下重要關(guān)系式：

當(dāng)給定壁板和氣流參數(shù)時(shí)，便可確定R與p，因此α1、β1和?中只有一個(gè)獨(dú)立變量，這里選?為獨(dú)立變量。同時(shí)k也取決于?，根據(jù)式(24）可知Ω取決于k。將式(19）與式(26）代入表1所示的邊界條件，可求得頻率方程與本征函數(shù)的系數(shù)(見(jiàn)表2），進(jìn)而求出時(shí)間本征值Ω，其反應(yīng)了壁板的振動(dòng)特性。

表2 二維壁板顫振頻率方程和本征函數(shù)Table 2 Eigensolutions of two-dim ensional panel flutters

2.2 顫振機(jī)理

顫振的發(fā)生必然是因?yàn)楸诎逭駝?dòng)過(guò)程中有氣動(dòng)力的作用。下面通過(guò)分析顫振方程和本征根Ω的性質(zhì)來(lái)判斷壁板為何會(huì)發(fā)生顫振。

3 算例與影響因素參數(shù)分析

首先，針對(duì)固有頻率和顫振參數(shù)，把精確解與有限元分析結(jié)果和Galerkin方法結(jié)果進(jìn)行比較，以驗(yàn)證精確解的正確性，再對(duì)影響顫振特性的參數(shù)進(jìn)行分析。

3.1 典型算例的計(jì)算結(jié)果比較

選擇如表3中所示陶瓷和金屬材料組成的梯度材料。只考慮彈性模量E、熱膨脹系數(shù)α與溫度相關(guān)，熱傳導(dǎo)系數(shù)為Kc＝9.19 W／(m·K）和Km＝12.04 W／(m·K），其中c表示陶瓷材料，m表示金屬材料。考慮非線性溫度場(chǎng)，上表面溫度350 K，下表面溫度為300 K，梯度指數(shù)n為5，弦長(zhǎng)a為0.5 m，跨厚比a／h＝250。

表3 功能梯度材料彈性常數(shù)Table 3 Elastic constants of FGM

1）與有限元結(jié)果的比較

現(xiàn)有有限元商業(yè)軟件如ANSYS、ABAQUS等都未提供成熟的功能梯度材料分析模塊。功能梯度板最主要的結(jié)構(gòu)特點(diǎn)在于其材料特性在板厚度方向不斷變化，因而在有限元模型中，將功能梯度板認(rèn)為成層合板。層合板的每層均為各向同性材料，每層的材料參數(shù)由該層所處的位置確定。考慮了總層數(shù)分別為5層和10層2種情況。表4及表5給出了根據(jù)式(1）計(jì)算得到的各分層密度、彈性模量。

表4 分層為5層時(shí)密度和彈性模量Table 4 Density and m odulus of elasticity for 5 layers

表5 分層為10層時(shí)密度和彈性模量Table 5 Density and m odulus of elasticity for 10 layers

表6將所得功能梯度板的前兩階固有頻率與ABAQUS軟件計(jì)算所得結(jié)果進(jìn)行了對(duì)比。有限元網(wǎng)格的劃分為10×2 500，且采用平面應(yīng)力單元。可以看出，精確解與ABAQUS結(jié)果存在差異，當(dāng)功能梯度板的分層增加時(shí)，差異在減小。當(dāng)分層數(shù)為10層時(shí)，SC邊界情況的第2階頻率之間的差異已經(jīng)小于0.413%。由此可以推斷，建立的FGM有限元模型是可行的，所得精確固有頻率是正確的。

表6 不同邊界條件下FGM 板頻率Table 6 Frequency of FGM p late under different boundary conditions

2）與Galerkin方法結(jié)果的比較

在顫振理論分析領(lǐng)域，Galerkin方法的應(yīng)用是非常廣泛的。下面把精確解與Galerkin結(jié)果進(jìn)行比較，以驗(yàn)證所得顫振結(jié)果的正確性?？紤]一均勻壁板，其剛度為D＝148.73 N·m2，密度為ρ＝1 600 kg／m3，弦長(zhǎng)a為0.3 m，跨厚比a／h＝200。不考慮溫度場(chǎng)作用，邊界條件為兩端簡(jiǎn)支。Galerkin方法選用前三階簡(jiǎn)支模態(tài)。表7給出了一、二階固有頻率ω1、ω2，以及顫振時(shí)對(duì)應(yīng)的顫振頻率ωf及臨界馬赫數(shù)Maf?？梢钥闯?，二者結(jié)果吻合較好，驗(yàn)證了精確顫振頻率和馬赫數(shù)的正確性。如果Galerkin方法選用更多階的模態(tài)，則二者結(jié)果將更加接近。

表7 本文解與Galerkin方法結(jié)果的對(duì)比Table 7 Com parison of result between present method and Galerkin method

3.2 參數(shù)分析

下面討論梯度指數(shù)、溫度場(chǎng)、邊界條件等對(duì)FGM板的臨界動(dòng)壓和頻率的影響。

3.2.1 梯度指數(shù)n對(duì)臨界動(dòng)壓的影響

功能梯度材料屬于非均勻材料，其組分是隨坐標(biāo)變化的。由式(1）可知，梯度指數(shù)n決定了2種材料成分各占多少。對(duì)于上表面是陶瓷，下表面是金屬的壁板，當(dāng)n為0時(shí)，壁板為各向同性的純陶瓷壁板。當(dāng)n為無(wú)窮大時(shí)表示純金屬壁板。采用表3所示材料，考慮兩邊簡(jiǎn)支邊界條件，a／h＝250，在不考慮溫度場(chǎng)的情況下研究無(wú)量綱臨界動(dòng)壓λ*隨梯度指數(shù)n的變化關(guān)系，λ*＝ρa(bǔ)V2a3／(MaD0），D0＝Emh3／［12(1－υ2）］，Em表示金屬材料在常溫T＝300 K下的彈性模量。如圖2所示，臨界動(dòng)壓值隨梯度指數(shù)的增加出現(xiàn)先下降較快，而后下降變緩的現(xiàn)象，當(dāng)指數(shù)較大時(shí)，臨界動(dòng)壓值趨于穩(wěn)定。

圖2 臨界動(dòng)壓值隨梯度指數(shù)變化曲線Fig.2 Critical dynamic pressure versus gradient index

3.2.2 溫度對(duì)臨界顫振頻率的影響

圖3 臨界顫振頻率隨溫度的變化曲線Fig.3 Critical flutter frequency versus temperature

當(dāng)ΔT＝100 K、n＝5、a／h＝250時(shí)，在兩端邊界為簡(jiǎn)支的條件下，圖4給出FGM板在均勻溫度場(chǎng)及非線性溫度場(chǎng)下的一階頻率與二階頻率隨著馬赫數(shù)的變化情況。可以得出，2種溫度場(chǎng)下均發(fā)生頻率耦合型顫振，在均勻溫度場(chǎng)下，壁板的臨界顫振頻率為408.24 Hz，臨界馬赫數(shù)為2.90。而在非線性溫度場(chǎng)下，壁板的臨界顫振頻率為423.83 Hz，較均勻溫度場(chǎng)時(shí)提高了3.82%，臨界馬赫數(shù)為3.02，較均勻溫度場(chǎng)時(shí)提高了4.14%。

圖4 兩種熱環(huán)境下頻率隨馬赫數(shù)的變化曲線Fig.4 Frequency versus Mach number in two thermal environments

3.2.3 邊界條件對(duì)臨界動(dòng)壓和顫振頻率的影響

表8為簡(jiǎn)支、固支及其組合邊界條件下的功能梯度板的臨界顫振頻率和臨界動(dòng)壓的比較。此處分析了均勻溫度場(chǎng)及非線性溫度場(chǎng)下，3種邊界條件下的無(wú)量綱臨界顫振頻率及臨界動(dòng)壓的大小，在相同條件下，三者的λ*與ω*的關(guān)系為：兩端固支大于一端簡(jiǎn)支一端固支，后者又大于兩端簡(jiǎn)支。

表8 不同邊界下FGM 板的臨界顫振頻率和臨界動(dòng)壓的比較Table 8 Com parison of flutter frequency and critical dynam icpressure of FGM plate under different boundary conditions

4 結(jié)論

基于經(jīng)典薄板理論和一階活塞理論建立了氣動(dòng)力作用下二維功能梯度壁板的熱顫振模型，采用分離變量法得到了顫振問(wèn)題的精確解。結(jié)論如下：

1）系統(tǒng)剛度由結(jié)構(gòu)彈性剛度和氣動(dòng)剛度組成。從數(shù)學(xué)上而言，顫振現(xiàn)象的發(fā)生是由于氣動(dòng)剛度致使系統(tǒng)本征根變成了復(fù)數(shù)，其實(shí)部決定系統(tǒng)振動(dòng)是衰減、等幅還是發(fā)散。

2）功能梯度材料能夠有效提高熱環(huán)境下壁板的顫振邊界。

3）在相同條件下，3種邊界條件下的臨界顫振頻率及臨界動(dòng)壓大小關(guān)系為：兩端固支大于一端簡(jiǎn)支一端固支，后者又大于兩端簡(jiǎn)支。