韓明君， 李金佩， 王 鵬

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

高超音速飛行器由于其具備高速度、高機(jī)動(dòng)、突防能力強(qiáng)、被探測(cè)到概率低、作戰(zhàn)半徑和效能高等諸多優(yōu)點(diǎn)，已經(jīng)成為各大國(guó)下一代航天器、飛機(jī)、導(dǎo)彈等不可或缺的條件和主要的發(fā)展方向。憑借速度和高度的絕對(duì)優(yōu)勢(shì)，高超音速飛行器可以在短時(shí)間內(nèi)完成高難度的情報(bào)、監(jiān)視、偵察和打擊任務(wù)，已經(jīng)引起了國(guó)際社會(huì)對(duì)其廣泛的關(guān)注和討論。尤其是對(duì)于飛行器中的壁板結(jié)構(gòu)，由于其自身的彈性變形會(huì)受到溫度荷載、氣動(dòng)力、以及材料自身屬性的影響，在這些作用的影響下，壁板會(huì)產(chǎn)生復(fù)雜的非線性動(dòng)力學(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]。綜合來(lái)看，金屬多孔材料能大幅度提升高超音速飛行器的各項(xiàng)參數(shù)及性能，是一項(xiàng)值得關(guān)注的研究方向。

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

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

1 模型建立與方程處理

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

圖1 運(yùn)動(dòng)模型圖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薄板大撓度理論，假定溫度場(chǎng)是準(zhǔn)定常的，并且只考慮板的橫向振動(dòng)，應(yīng)變位移關(guān)系為

(10)

(11)

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

(12)

(13)

對(duì)于無(wú)限展長(zhǎng)的二維壁板，可以忽略εy，由式(10)～式(13)可得

(14)

式中，NT為ΔT引起的溫度應(yīng)力?？紤]Kirchhoff假設(shè)，壁板的振動(dòng)微分方程為

(15)

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

(16)

(17)

(18)

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

(19)

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

(20)

(21)

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

(22)

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

(23)

(24)

(25)

對(duì)于孔隙均勻分布

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

(26)

對(duì)于孔隙幾何中面均勻分布

(27)

對(duì)于孔隙幾何中面非均勻分布

(28)

2.2 微分方程的無(wú)量綱化

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

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

(29)

(30)

(31)

考慮對(duì)邊簡(jiǎn)支壁板在x=0和x=1處的無(wú)量綱邊界條件，其形式如下

(32)

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

(33)

兩端簡(jiǎn)支的二維壁板的邊界條件為

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

(34)

ψ(4)+λ2ψ″=0

(35)

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

(36)

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

(37)

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

csinλ=0或λ=mπ

(38)

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

(39)

由式(39)可以得到

(40)

從而得到

(41)

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

在第一階臨界屈曲載荷Q1之前壁板都是穩(wěn)定的；隨著軸向力Q的增加，當(dāng)Q=Q1時(shí)，壁板開(kāi)始出現(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í)，壁板就會(huì)存在三個(gè)屈曲構(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不同的孔隙分布對(duì)壁板的整體剛度、溫度應(yīng)力和面內(nèi)力的影響不同而出現(xiàn)的情況。

4 Galerkin法求解振動(dòng)問(wèn)題

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

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

(42)

(43)

(44)

其中，

(45)

(46)

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

(47)

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

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

(48)

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

(49)

得到一個(gè)關(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可得無(wú)量綱臨界流速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í)無(wú)量綱臨界流速和無(wú)量綱臨界頻率等數(shù)據(jù)Tab.3 Critical velocity and critical frequency data at RT0=3π2/4

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

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

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

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

圖5 RT0=π2/4時(shí)的無(wú)量綱臨界流速和無(wú)量綱臨界頻率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)在無(wú)量綱動(dòng)壓到達(dá)P1時(shí)將發(fā)生顫振。當(dāng)無(wú)量綱動(dòng)壓小于P1，超音速流中受熱平壁板在平衡點(diǎn)是穩(wěn)定的，而當(dāng)無(wú)量綱動(dòng)壓大于P1時(shí)，超音速流中受熱平壁板在平衡點(diǎn)是不穩(wěn)定的。為了驗(yàn)證以上結(jié)論，計(jì)算了不同參數(shù)時(shí)對(duì)應(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)在無(wú)量綱動(dòng)壓到達(dá)P1時(shí)將發(fā)生顫振，當(dāng)無(wú)量綱動(dòng)壓小于P1時(shí)，平衡點(diǎn)的四個(gè)特征值全部具有負(fù)實(shí)部，超音速流中受熱平壁板在平衡點(diǎn)是穩(wěn)定的；而當(dāng)無(wú)量綱動(dòng)壓大于P1時(shí)，平衡點(diǎn)的四個(gè)特征值中，一對(duì)復(fù)特征值具有負(fù)實(shí)部；另一對(duì)復(fù)特征值具有正實(shí)部,超音速流中受熱平壁板在平衡點(diǎn)是不穩(wěn)定的，在其附近產(chǎn)生穩(wěn)定的極限環(huán)。隨著孔隙率和溫度應(yīng)力的增大都會(huì)讓無(wú)量綱動(dòng)壓P1減小，從而使壁板有穩(wěn)定變?yōu)椴环€(wěn)定的趨勢(shì)。

表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è)建立了多孔壁板在超音速流中振動(dòng)的控制微分方程并無(wú)量綱化，然后運(yùn)用Galerkin法對(duì)梯度多孔壁板進(jìn)行變換并對(duì)氣動(dòng)弦長(zhǎng)進(jìn)行積分得到非線性方程，再利用Hurwitz行列式將Hopf分岔的判定轉(zhuǎn)化為非線性方程的求根。最后分別分析了各參數(shù)值的變化對(duì)無(wú)量綱臨界頻率和無(wú)量綱臨界流速的影響，并驗(yàn)證了梯度多孔薄壁板氣動(dòng)彈性的穩(wěn)定性，得到主要結(jié)論如下：

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

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

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