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

韓明君， 李金佩， 王 鵬

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

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

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

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

1 模型建立與方程處理

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

圖1 運(yùn)動模型圖Fig.1 Motion model

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

孔隙均勻分布

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

(1)

(2)

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

(3)

孔隙關(guān)于幾何中面均勻分布

(4)

(5)

(6)

孔隙關(guān)于幾何中面非均勻分布

(7)

(8)

(9)

2 基本方程的建立及處理

2.1 建立控制方程

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

(10)

(11)

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

(12)

(13)

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

(14)

式中，NT為ΔT引起的溫度應(yīng)力。考慮Kirchhoff假設(shè)，壁板的振動微分方程為

(15)

式中：Mx為彎矩；Nx為面內(nèi)力；mx為板的單位長度的質(zhì)量；qa為氣動力。其中

(16)

(17)

(18)

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

(19)

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

(20)

(21)

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

(22)

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

(23)

(24)

(25)

對于孔隙均勻分布

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

(26)

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

(27)

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

(28)

2.2 微分方程的無量綱化

引入如下無量綱參數(shù)

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

(29)

(30)

(31)

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

(32)

3 Galerkin法求解屈曲和靜態(tài)分岔問題

(33)

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

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

(34)

ψ(4)+λ2ψ″=0

(35)

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

(36)

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

(37)

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

csinλ=0或λ=mπ

(38)

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

(39)

由式(39)可以得到

(40)

從而得到

(41)

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

在第一階臨界屈曲載荷Q1之前壁板都是穩(wěn)定的；隨著軸向力Q的增加，當(dāng)Q=Q1時(shí)，壁板開始出現(xiàn)分岔，發(fā)生屈曲,當(dāng)Q>Q1時(shí)，原先在x=0處的平衡點(diǎn)變?yōu)椴黄胶恻c(diǎn); 隨著軸向力Q的繼續(xù)增加，當(dāng)Q大于第二階臨界屈曲載荷Q2時(shí),壁板存在三種平衡：靜平衡構(gòu)型、屈曲構(gòu)型和新的屈曲構(gòu)型；當(dāng)軸向載荷Q繼續(xù)增加，大于第三階臨界屈曲載荷Q3時(shí)，壁板就會存在三個屈曲構(gòu)型。

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

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

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

當(dāng)0<θ<0.3時(shí)，均勻分布比幾何中面非均勻分布先到達(dá)臨界點(diǎn)，出現(xiàn)分岔；當(dāng)0.3≤θ≤0.5時(shí)，幾何中面非均勻分布比均勻分布先到達(dá)臨界點(diǎn)，出現(xiàn)分岔；當(dāng)0<θ≤0.5，幾何中面均勻分布比其他二種分布都先到達(dá)臨點(diǎn)，出現(xiàn)分岔，這是因?yàn)槿N不同的孔隙分布對壁板的整體剛度、溫度應(yīng)力和面內(nèi)力的影響不同而出現(xiàn)的情況。

4 Galerkin法求解振動問題

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

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

(42)

(43)

(44)

其中，

(45)

(46)

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

(47)

4.2 系統(tǒng)發(fā)生Hopf分岔的邊界曲線

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

(48)

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

(49)

得到一個關(guān)于λ的四次方程

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

(50)

其中，

a0=1，

(51)

計(jì)算相應(yīng)的各階Hurwitz行列式

Δ1=a1

(52)

(53)

(54)

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

(55)

(56)

(57)

設(shè)向量X=(x1,x2,x3,x4和Y=(y1,y2,y3,y4)T分別是矩陣A的屬于特征值iω的歸一化的左、右特征向量。根據(jù)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)

此時(shí)系統(tǒng)發(fā)生Hopf分岔，即在U=Ucr時(shí)超音速流中的壁板發(fā)生顫振。

4.3 數(shù)值算例

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

表2 某鋁合金機(jī)械性能參數(shù)Tab.2 Mechanical property parameters of an aluminum alloy

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

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

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

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

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

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

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

4.4 穩(wěn)定性驗(yàn)證及分析

超音速流中受熱壁板系統(tǒng)在平衡點(diǎn)X0(0,0,0,0)處發(fā)生Hopf分岔。即超音速流中受熱壁板系統(tǒng)在無量綱動壓到達(dá)P1時(shí)將發(fā)生顫振。當(dāng)無量綱動壓小于P1，超音速流中受熱平壁板在平衡點(diǎn)是穩(wěn)定的，而當(dāng)無量綱動壓大于P1時(shí)，超音速流中受熱平壁板在平衡點(diǎn)是不穩(wěn)定的。為了驗(yàn)證以上結(jié)論，計(jì)算了不同參數(shù)時(shí)對應(yīng)的平衡點(diǎn)處Jacobi矩陣的特征值。

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

由表4、表5、表6和圖6可以得出并驗(yàn)證：超音速流中受熱壁板系統(tǒng)在無量綱動壓到達(dá)P1時(shí)將發(fā)生顫振，當(dāng)無量綱動壓小于P1時(shí)，平衡點(diǎn)的四個特征值全部具有負(fù)實(shí)部，超音速流中受熱平壁板在平衡點(diǎn)是穩(wěn)定的；而當(dāng)無量綱動壓大于P1時(shí)，平衡點(diǎn)的四個特征值中，一對復(fù)特征值具有負(fù)實(shí)部；另一對復(fù)特征值具有正實(shí)部,超音速流中受熱平壁板在平衡點(diǎn)是不穩(wěn)定的，在其附近產(chǎn)生穩(wěn)定的極限環(huán)。隨著孔隙率和溫度應(yīng)力的增大都會讓無量綱動壓P1減小，從而使壁板有穩(wěn)定變?yōu)椴环€(wěn)定的趨勢。

表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 三種分布的時(shí)間歷程圖、相位圖和龐加萊映射圖Fig.6 Time history diagram, phase diagram and Poincare map of three distributions

5 結(jié) 論

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

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

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

(3) 超音速流中受熱壁板系統(tǒng)在無量綱動壓到達(dá)P1時(shí)將發(fā)生顫振；當(dāng)無量綱動壓小于P1時(shí)，超音速流中受熱平壁板在平衡點(diǎn)是穩(wěn)定的；而當(dāng)無量綱動壓大于P1時(shí)，超音速流中受熱平壁板在平衡點(diǎn)是不穩(wěn)定的，在其附近產(chǎn)生穩(wěn)定的極限環(huán)。隨著孔隙率和溫度應(yīng)力的增大都會讓無量綱動壓P1減小，從而使壁板有穩(wěn)定變?yōu)椴环€(wěn)定的趨勢。