二維功能梯度壁板熱顫振本征問題的精確解

代林桐，邢譽(yù)峰

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

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

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

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

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

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

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

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

1 顫振微分方程

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

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

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

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

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

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

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

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

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

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

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

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

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

應(yīng)力為

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

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

FGM薄板自由振動的動能為

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

由式(6）可得

軸力為

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

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

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

式中：q為動壓。

氣動力所做虛功δW的積分為

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

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

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

本文用到的簡支和固支邊界條件如表1所示。

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

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

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

2.1 精確本征解

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

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

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

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

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

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

或?qū)懗?/p>

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

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

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

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

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

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

2.2 顫振機(jī)理

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

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

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

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

選擇如表3中所示陶瓷和金屬材料組成的梯度材料。只考慮彈性模量E、熱膨脹系數(shù)α與溫度相關(guān)，熱傳導(dǎo)系數(shù)為Kc＝9.19 W／(m·K）和Km＝12.04 W／(m·K），其中c表示陶瓷材料，m表示金屬材料?？紤]非線性溫度場，上表面溫度350 K，下表面溫度為300 K，梯度指數(shù)n為5，弦長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）計算得到的各分層密度、彈性模量。

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

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

表6將所得功能梯度板的前兩階固有頻率與ABAQUS軟件計算所得結(jié)果進(jìn)行了對比。有限元網(wǎng)格的劃分為10×2 500，且采用平面應(yīng)力單元?？梢钥闯?，精確解與ABAQUS結(jié)果存在差異，當(dāng)功能梯度板的分層增加時，差異在減小。當(dāng)分層數(shù)為10層時，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，弦長a為0.3 m，跨厚比a／h＝200。不考慮溫度場作用，邊界條件為兩端簡支。Galerkin方法選用前三階簡支模態(tài)。表7給出了一、二階固有頻率ω1、ω2，以及顫振時對應(yīng)的顫振頻率ωf及臨界馬赫數(shù)Maf。可以看出，二者結(jié)果吻合較好，驗(yàn)證了精確顫振頻率和馬赫數(shù)的正確性。如果Galerkin方法選用更多階的模態(tài)，則二者結(jié)果將更加接近。

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

3.2 參數(shù)分析

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

3.2.1 梯度指數(shù)n對臨界動壓的影響

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

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

3.2.2 溫度對臨界顫振頻率的影響

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

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

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

3.2.3 邊界條件對臨界動壓和顫振頻率的影響

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

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

4 結(jié)論

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

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

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

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