超音速流-熱載作用下梯度多孔壁板的非線性振動特性分析

韓明君， 李金佩， 王 鵬

(蘭州理工大學 理學院, 蘭州 730050)

高超音速飛行器由于其具備高速度、高機動、突防能力強、被探測到概率低、作戰半徑和效能高等諸多優點，已經成為各大國下一代航天器、飛機、導彈等不可或缺的條件和主要的發展方向。憑借速度和高度的絕對優勢，高超音速飛行器可以在短時間內完成高難度的情報、監視、偵察和打擊任務，已經引起了國際社會對其廣泛的關注和討論。尤其是對于飛行器中的壁板結構，由于其自身的彈性變形會受到溫度荷載、氣動力、以及材料自身屬性的影響，在這些作用的影響下，壁板會產生復雜的非線性動力學響應，從而容易使壁板結構產生疲勞損傷，進而影響飛行器結構的疲勞壽命[1-2]。多孔金屬材料與傳統質密金屬材料相比，具有超輕、高比強、高比剛度、高強韌、高能量吸收等優良機械性能, 以及減震、散熱、吸聲、電磁屏蔽等特殊性質, 它兼具功能和結構雙重作用,是一種性能優異的多功能材料[3]。綜合來看，金屬多孔材料能大幅度提升高超音速飛行器的各項參數及性能，是一項值得關注的研究方向。

近年來，國內外許多學者對超聲速氣流中壁板非線性氣動彈性響應行為開展了許多研究。Cheng等[4]運用超高速氣流作用下壁板非線性顫振的有限元時域模態， 分析了混沌狀態下壁板的屈曲、極限環振蕩、周期以及時域圖和相位圖等。Ghoman等[5-6]基于有限元方法對壁板顫振問題進行了研究。Koo等[7]用有限元方法討論了阻尼對復合材料板顫振性能的影響。Ibrahim等[8-11]在超音速氣流中考慮了二維壁板的剪切變形，分析了壁板的氣動彈性模態。他們也考慮了平面內熱載荷對板殼結構橫向彎曲撓度的影響，研究了均勻溫度場中復合材料板殼的顫振和熱屈曲問題。李麗麗等[12]在非定常氣動力壁板模型中引進了熱荷載的影響，得到了熱壁板的顫振方程。楊智春等[13]基于Von Kármán非線性板理論和Galerkin方法建立超音速氣流中二維受熱壁板的氣動彈性模型。王廣勝等[14-15]對超音速流中飛行器的受熱平壁板和曲壁板非線性氣動彈性穩定性進行了深入研究，分析熱氣動彈性系統的顫振邊界特性以及不同的參數組合對系統顫振臨界動壓與穩定性的影響。周建等[16]用POD方法構造三維復合材料曲壁板顫振響應模態，通過數值積分方法計算三維復合材料曲壁板的顫振響應，計算結果與傳統的模態縮減法取得一致。梅冠華等[17]采用流-固耦合算法研究了曲壁板在跨/超聲速氣流下的氣動彈性特征，為高速飛行器的壁板設計和顫振抑制提供了依據。

綜上所述，近年來主要集中研究傳統質密金屬材料的氣動穩定性，對多孔金屬材料的氣動彈性鮮有研究。本文基于Von Kármán薄板大撓度理論和Kirchhoff假設，考慮孔隙沿矩形板厚度方向呈三種不同的分布方式，利用Galerkin法研究了溫差、孔隙率、氣動剛度系數對梯度多孔薄壁板氣動彈性的影響，并用Runge-Kutta法給出了系統的時間歷程圖、相位圖和龐加萊截面映射圖。

1 模型建立與方程處理

建立梯度多孔二維壁板在超音速氣流中的力學模型，如圖1所示。假設溫度在厚度方向均勻分布，溫度變化為ΔT的簡支無限展長的壁板模型，其長度為L,厚度為h，密度為ρ(z)，彈性模量為E(z)，α(z)為材料的熱膨脹系數，μ為沿矩形板厚度方向恒定的泊松比[18]。在其上表面有沿x方向的超音速氣流,流速為U∞，馬赫數為M∞，空氣密度為ρ∞。

圖1 運動模型圖Fig.1 Motion model

本文考慮了孔隙沿厚度方向呈均勻分布和兩種非均勻分布[19-21]

孔隙均勻分布

E(z)=E1(1-θγ)

(1)

(2)

α(z)=α1(1-amγ)

(3)

孔隙關于幾何中面均勻分布

(4)

(5)

(6)

孔隙關于幾何中面非均勻分布

(7)

(8)

(9)

2 基本方程的建立及處理

2.1 建立控制方程

基于Von Kármán薄板大撓度理論，假定溫度場是準定常的，并且只考慮板的橫向振動，應變位移關系為

(10)

(11)

式中，u,v和w分別為x,y和z軸的位移。根據Hooke定律，應變表示為

(12)

(13)

對于無限展長的二維壁板，可以忽略εy，由式(10)～式(13)可得

(14)

式中，NT為ΔT引起的溫度應力?？紤]Kirchhoff假設，壁板的振動微分方程為

(15)

式中：Mx為彎矩；Nx為面內力；mx為板的單位長度的質量；qa為氣動力。其中

(16)

(17)

(18)

將E(z),ρ(z)和α(z)的表達式和式(10)、式(14)和式(16)～ 式(18)代入振動微分方程式(15)中得到

(19)

式中：D0為矩形板的最大抗彎剛度；ek1為對參變量的影響，k=1,2,3為不同孔隙分布，由于多孔壁板的邊界為二端簡支，則u|x=0.L=0,邊界處的彎矩為零

(20)

(21)

又由于面內力為常數，對其取x方向的平均值，由式(14)和式(21)，振動微分方程式(19)中的面內力可以表示為

(22)

將式(22)代入振動微分方程式(19)中，可以得到多孔壁板的振動微分方程

(23)

(24)

(25)

對于孔隙均勻分布

e11=1-θγ,e12=1-θγ,

(26)

對于孔隙幾何中面均勻分布

(27)

對于孔隙幾何中面非均勻分布

(28)

2.2 微分方程的無量綱化

引入如下無量綱參數

將無量綱參數代入多孔壁板的振動微分方程式(22)中，得到無量綱微分方程

(29)

(30)

(31)

考慮對邊簡支壁板在x=0和x=1處的無量綱邊界條件，其形式如下

(32)

3 Galerkin法求解屈曲和靜態分岔問題

(33)

兩端簡支的二維壁板的邊界條件為

ψ|x=0,1=0，ψ″|x=0,1=0

(34)

ψ(4)+λ2ψ″=0

(35)

無量綱微分方程式(35)的通解為

(36)

將式(36)代入二維壁板的邊界條件式(34)中可得

(37)

式(36)中ψ(x)應具有非零解，因此線性齊次方程的系數行列式為零化簡可得特征方程，由式(37)的行列式可得特征方程的解為

csinλ=0或λ=mπ

(38)

λ2=Q-Γ=Q-3(ek2/ek1)c2λ2

(39)

由式(39)可以得到

(40)

從而得到

(41)

為了研究熱應力對二維壁板的影響，本文通過算例展開討論。取壁板的材料特性參數E1=78.55 GPa，μ=0.3，L=0.5 m，h=0.002 m。選取壁板1/4長度處進行計算，得到其前三階彎曲構型的靜態分岔圖如圖2所示。

在第一階臨界屈曲載荷Q1之前壁板都是穩定的；隨著軸向力Q的增加，當Q=Q1時，壁板開始出現分岔，發生屈曲,當Q>Q1時，原先在x=0處的平衡點變為不平衡點; 隨著軸向力Q的繼續增加，當Q大于第二階臨界屈曲載荷Q2時,壁板存在三種平衡：靜平衡構型、屈曲構型和新的屈曲構型；當軸向載荷Q繼續增加，大于第三階臨界屈曲載荷Q3時，壁板就會存在三個屈曲構型。

由圖2和表1計算可知：當孔隙率θ=0時，三種不同孔隙分布壁板的臨界屈曲載荷數值相同的，并且與王廣勝的研究結果一致；隨著孔隙率的增大，三種不同孔隙分布壁板的前三階彎曲構型的臨界屈曲載荷也隨之增加，這是壁板中的孔隙削弱了壁板的整體剛度，同時也減小了壁板中的溫度應力和面內力，所以前三階彎曲構型的臨界屈曲載荷隨著孔隙率θ的增大而增加。

圖2 不同孔隙分布的前三階彎曲構型的靜態分岔圖Fig.2 Static bifurcation diagrams of the first three order bending configurations with different pore distributions

表1 不同孔隙率前三階彎曲構型的靜態分岔臨界屈曲載荷Tab.1 The static bifurcation critical buckling load of the first three order bending configuration with different porosity

當0<θ<0.3時，均勻分布比幾何中面非均勻分布先到達臨界點，出現分岔；當0.3≤θ≤0.5時，幾何中面非均勻分布比均勻分布先到達臨界點，出現分岔；當0<θ≤0.5，幾何中面均勻分布比其他二種分布都先到達臨點，出現分岔，這是因為三種不同的孔隙分布對壁板的整體剛度、溫度應力和面內力的影響不同而出現的情況。

4 Galerkin法求解振動問題

4.1 非線性常微分方程組的處理

利用Galerkin法，將簡支邊界的位移函數表示為正弦函數的線性組合

(42)

(43)

(44)

其中，

(45)

(46)

非線性常微分方程式(46)寫成矩陣形式

(47)

4.2 系統發生Hopf分岔的邊界曲線

參考文獻[22-24]利用Hurwitz行列式，將Hopf分岔的判定轉化為非線性方程求根的問題，解決Hopf分岔的代數判據和計算。假設系統的一個平衡點為X0(0,0,0,0)，則該平衡點處的Jacobi矩陣為

(48)

對應的特征方程為det(JA-λI)=0，即

(49)

得到一個關于λ的四次方程

a0λ4+a1λ3+a2λ2+a3λ+a4=0

(50)

其中，

a0=1，

(51)

計算相應的各階Hurwitz行列式

Δ1=a1

(52)

(53)

(54)

由Δ3=0可得無量綱臨界流速Ucr為

(55)

(56)

(57)

設向量X=(x1,x2,x3,x4和Y=(y1,y2,y3,y4)T分別是矩陣A的屬于特征值iω的歸一化的左、右特征向量。根據XA=iωA,YA=iωA,XA=1，可得

(58)

(59)

A2=33π4ek1ek4+π4ek1-9π2ek3ek4RT0,

(60)

(61)

ζ1,2(P1)=α1(P1)±iω(P1)

(62)

式中，α1(P1)=0,ω(P1)>0。于是X(P1)A(P1)·Y(P1)=α1(P1)±iω(P1)，得到

(63)

a1>0,a2>0,…,an>0
Δn-1=0,Δi>0(i=n-3,n-5,…)

(64)

此時系統發生Hopf分岔，即在U=Ucr時超音速流中的壁板發生顫振。

4.3 數值算例

選取楊智春教授等研究的某鋁合金壁板材料模型的計算參數，其機械性能參數如表2所示。

表2 某鋁合金機械性能參數Tab.2 Mechanical property parameters of an aluminum alloy

假設飛行器的飛行高度為1.1 km，此時空氣密度0.364 kg/m3，音速295.065 m/s，馬赫數為5。當孔隙率為0時，代入數據可以得到各參數的數值，表3中本文退化數據與曹麗娜研究中的數據基本一致。

表3 RT0=3π2/4時無量綱臨界流速和無量綱臨界頻率等數據Tab.3 Critical velocity and critical frequency data at RT0=3π2/4

取相同材料特性參數,得到隨著孔隙率增長不同分布和不同溫度應力的壁板的無量綱臨界流速和無量綱臨界頻率。

由圖3、圖4、圖5可知，當孔隙率一定時，隨著溫度應力的減小，三種不同分布的壁板發生顫振的無量綱臨界流速在不斷的增大；當溫度應力一定時，隨著孔隙率的增長三種不同分布壁板的無量綱臨界流速都在下降，其中均勻分布中壁板無量綱臨界流速隨著孔隙率的增加而下降的程度最大，幾何中面非均勻分布壁板無量綱臨界流速隨著孔隙率的增加而下降的程度次之，幾何中面均勻分布壁板無量綱臨界流速隨著孔隙率的增加而的程度最小。壁板無量綱臨界頻率隨著不同孔隙率和溫度應力的變化比較復雜。當孔隙率一定時，隨著溫度應力的減小，三種不同分布的壁板無量綱臨界頻率都在在不斷增大；隨著孔隙率的增長三種不同分布壁板的無量綱臨界頻率也與它初始溫度應力有關系，當溫度應力小于3π2/4時，幾何中面均勻分布壁板無量綱臨界頻率隨著孔隙率的增加而增大，均勻分布壁板和幾何中面非均勻分布壁板的無量綱臨界頻率隨著孔隙率的增加而減?。划敎囟葢Υ笥?π2/4時，幾何中面均勻分布壁板和幾何中面非均勻分布壁板的無量綱臨界頻率隨著孔隙率的增加而增大，均勻分布壁板的無量綱臨界頻率隨著孔隙率的增加先增加后減小，這是由于對三種不同孔隙分布壁板的整體剛度、溫度應力、等效質量和面內力減弱程度不同，三種不同孔隙分布壁板的無量綱臨界流速和無量綱臨界頻率變化趨勢也不一樣。

圖3 RT0=5π2/4時的無量綱臨界流速和無量綱臨界頻率Fig.3 Dimensionless critical velocity and dimensionless critical frequency at RT0=5π2/4

圖4 RT0=3π2/4時的無量綱臨界流速和無量綱臨界頻率Fig.4 Dimensionless critical velocity and dimensionless critical frequency at RT0=3π2/4

圖5 RT0=π2/4時的無量綱臨界流速和無量綱臨界頻率Fig.5 Dimensionless critical velocity and dimensionless critical frequency at RT0=π2/4

4.4 穩定性驗證及分析

超音速流中受熱壁板系統在平衡點X0(0,0,0,0)處發生Hopf分岔。即超音速流中受熱壁板系統在無量綱動壓到達P1時將發生顫振。當無量綱動壓小于P1，超音速流中受熱平壁板在平衡點是穩定的，而當無量綱動壓大于P1時，超音速流中受熱平壁板在平衡點是不穩定的。為了驗證以上結論，計算了不同參數時對應的平衡點處Jacobi矩陣的特征值。

當θ=0.05,P1=190,RT0=3π2/4時，初值取(0.005, 0.005,0.01,0.01)，給出了系統的時間歷程圖、相位圖和龐加萊截面映射圖。

由表4、表5、表6和圖6可以得出并驗證：超音速流中受熱壁板系統在無量綱動壓到達P1時將發生顫振，當無量綱動壓小于P1時，平衡點的四個特征值全部具有負實部，超音速流中受熱平壁板在平衡點是穩定的；而當無量綱動壓大于P1時，平衡點的四個特征值中，一對復特征值具有負實部；另一對復特征值具有正實部,超音速流中受熱平壁板在平衡點是不穩定的，在其附近產生穩定的極限環。隨著孔隙率和溫度應力的增大都會讓無量綱動壓P1減小，從而使壁板有穩定變為不穩定的趨勢。

表4 θ=0.5, RT0=π2/4的Jacobi矩陣的特征值Tab.4 Eigenvalues of Jacobi matrices of θ=0.5, RT0=π2/4

表5 θ=0.5, RT0=3π2/4的Jacobi矩陣的特征值Tab.5 Eigenvalues of Jacobi matrices of θ=0.5, RT0=3π2/4

表6 θ=0.3,RT0=3π2/4的Jacobi矩陣的特征值Tab.6 Eigenvalues of Jacobi matrices of θ=0.3,RT0=3π2/4

圖6 三種分布的時間歷程圖、相位圖和龐加萊映射圖Fig.6 Time history diagram, phase diagram and Poincare map of three distributions

5 結 論

本文基于von Kármán薄板大撓度理論和Kirchhoff假設建立了多孔壁板在超音速流中振動的控制微分方程并無量綱化，然后運用Galerkin法對梯度多孔壁板進行變換并對氣動弦長進行積分得到非線性方程，再利用Hurwitz行列式將Hopf分岔的判定轉化為非線性方程的求根。最后分別分析了各參數值的變化對無量綱臨界頻率和無量綱臨界流速的影響，并驗證了梯度多孔薄壁板氣動彈性的穩定性，得到主要結論如下：

(1) 隨著孔隙率的增大，三種不同孔隙分布壁板的前三階彎曲構型的臨界屈曲載荷也隨之增加，這是由于壁板中的孔隙削弱了壁板的整體剛度，同時也減小了壁板中的溫度應力和面內力，前三階彎曲構型的臨界屈曲載荷隨著孔隙率的增大而增加。

(2) 當孔隙率一定時，隨著溫度應力的減小，三種不同分布的壁板發生顫振的無量綱臨界流速在不斷的增大；當溫度應力一定時，隨著孔隙率的增長三種不同分布壁板的無量綱臨界流速都在下降；當孔隙率一定時，隨著溫度應力的減小，三種不同分布的壁板無量綱臨界頻率都在在不斷增大。

(3) 超音速流中受熱壁板系統在無量綱動壓到達P1時將發生顫振；當無量綱動壓小于P1時，超音速流中受熱平壁板在平衡點是穩定的；而當無量綱動壓大于P1時，超音速流中受熱平壁板在平衡點是不穩定的，在其附近產生穩定的極限環。隨著孔隙率和溫度應力的增大都會讓無量綱動壓P1減小，從而使壁板有穩定變為不穩定的趨勢。