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

代林桐，邢譽峰

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

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

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

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

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

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

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

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

1 顫振微分方程

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

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

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

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

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

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

由式(2）可知，要獲得有效物理屬性，除了知道其對應的溫度系數以外，還需要知道結構溫度場分布。當材料熱傳導系數為常數時，由一維熱傳導穩態方程所得到的是一個線性的溫度場。對于功能梯度材料組成的板而言，當熱傳導系數k(i）是厚度坐標的函數而不再是一個常數時，這時得到一個非線性的溫度場。該非線性溫度場滿足如下一維熱傳導方程：

相關文獻給出了其近似理論解［16］。

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

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

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

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

FGM薄板的應變可以表示為

應力為

式中：εx(T）為溫度應變，即

熱環境中FGM矩形薄板的總應變能為

FGM薄板自由振動的動能為

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

由式(6）可得

軸力為

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

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

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

式中：q為動壓。

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

根據Ham ilton原理的廣義形式為

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

其中：Deq為等效剛度；為熱應力。

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

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

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

2 顫振精確本征解和顫振機理

2.1 精確本征解

設撓度w有如下分離變量的形式：

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

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

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

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

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

或寫成

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

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

則本征函數φ的通解可表示為

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

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

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

2.2 顫振機理

顫振的發生必然是因為壁板振動過程中有氣動力的作用。下面通過分析顫振方程和本征根Ω的性質來判斷壁板為何會發生顫振。

3 算例與影響因素參數分析

首先，針對固有頻率和顫振參數，把精確解與有限元分析結果和Galerkin方法結果進行比較，以驗證精確解的正確性，再對影響顫振特性的參數進行分析。

3.1 典型算例的計算結果比較

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

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

1）與有限元結果的比較

現有有限元商業軟件如ANSYS、ABAQUS等都未提供成熟的功能梯度材料分析模塊。功能梯度板最主要的結構特點在于其材料特性在板厚度方向不斷變化，因而在有限元模型中，將功能梯度板認為成層合板。層合板的每層均為各向同性材料，每層的材料參數由該層所處的位置確定。考慮了總層數分別為5層和10層2種情況。表4及表5給出了根據式(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軟件計算所得結果進行了對比。有限元網格的劃分為10×2 500，且采用平面應力單元。可以看出，精確解與ABAQUS結果存在差異，當功能梯度板的分層增加時，差異在減小。當分層數為10層時，SC邊界情況的第2階頻率之間的差異已經小于0.413%。由此可以推斷，建立的FGM有限元模型是可行的，所得精確固有頻率是正確的。

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

2）與Galerkin方法結果的比較

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

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

3.2 參數分析

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

3.2.1 梯度指數n對臨界動壓的影響

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

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

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

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

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

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

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

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

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

4 結論

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

1）系統剛度由結構彈性剛度和氣動剛度組成。從數學上而言，顫振現象的發生是由于氣動剛度致使系統本征根變成了復數，其實部決定系統振動是衰減、等幅還是發散。

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

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