范晨光,　楊翊仁

(西南交通大學(xué)力學(xué)與工程學(xué)院, 四川 成都 610031)

薄殼作為一類常見(jiàn)的超聲速飛行器結(jié)構(gòu),在制造工藝流程中或服役階段會(huì)出現(xiàn)幾何或物理上的缺陷.初始幾何缺陷主要是指制造和裝配過(guò)程中產(chǎn)生的加工誤差、焊接殘余變形等初始撓度.對(duì)于近空間高速飛行薄殼而言,由于其惡劣的服役環(huán)境,初始缺陷可能會(huì)對(duì)結(jié)構(gòu)的安全性和耐久性產(chǎn)生嚴(yán)重的影響,導(dǎo)致其在遠(yuǎn)低于臨界動(dòng)壓的情況下顫振失穩(wěn),并由此帶來(lái)災(zāi)難性后果.

通常情況下,在材料尚未達(dá)到屈服強(qiáng)度前,薄殼會(huì)在靜力或動(dòng)力荷載作用下破壞,這種破壞呈現(xiàn)強(qiáng)烈的突然性,失穩(wěn)的形式主要包括屈曲和顫振.由于初始缺陷的存在,這類結(jié)構(gòu)可能會(huì)在低于其臨界失穩(wěn)速度的設(shè)計(jì)值時(shí)產(chǎn)生破壞.很早人們發(fā)現(xiàn)含有初始缺陷的薄殼其產(chǎn)生屈曲破壞的臨界荷載要比不含初始缺陷的薄殼低得多,并開(kāi)始重視薄殼穩(wěn)定性對(duì)初始缺陷的敏感度問(wèn)題.Koiter的研究表明[1-2],在薄殼存在初始撓度情況下,受軸壓的圓柱殼和與外壓下的球殼均有較大的缺陷敏感度.Arboz等[3]研究了圓柱殼地震載荷作用的缺陷敏感度問(wèn)題,得到了與Koiter相似的結(jié)論.Simitses[4]綜述了更多的圓柱殼穩(wěn)定性缺陷敏感度問(wèn)題,收集了129篇參考文獻(xiàn),詳細(xì)闡述了缺陷圓柱殼的屈曲和后屈曲問(wèn)題,涉及的缺陷主要有初始幾何缺陷、載荷偏心、材料或結(jié)構(gòu)缺陷.國(guó)內(nèi)在薄殼初始缺陷敏感度問(wèn)題上的研究也很多,沈惠申等[5]建立了圓柱薄殼外壓作用下的邊界層理論,同時(shí)考慮了前屈曲非線性,大撓度和初始幾何缺陷.杜啟端[6]對(duì)此做了較為全面的綜述.

與薄殼屈曲問(wèn)題一樣,在薄殼氣動(dòng)彈性力學(xué)問(wèn)題上,也是由于理論預(yù)測(cè)的線性系統(tǒng)顫振臨界動(dòng)壓與實(shí)驗(yàn)結(jié)果的顯著差異,促使研究者產(chǎn)生了對(duì)薄殼初始缺陷敏感度的研究動(dòng)機(jī).最早在薄殼氣動(dòng)彈性力學(xué)問(wèn)題上考慮初始缺陷影響的是Barr等[7]研究的模型為大撓度完全圓柱殼,缺陷類型為軸對(duì)稱撓度缺陷.他們的研究表明,初始撓度缺陷會(huì)嚴(yán)重影響圓柱殼的顫振特性,使其更加容易失穩(wěn).Amabili等[8]研究了軸對(duì)稱和非軸對(duì)稱兩種初始撓度缺陷類型對(duì)圓柱殼氣動(dòng)彈性問(wèn)題的影響,結(jié)果表明,圓柱殼顫振臨界動(dòng)壓對(duì)初始撓度缺陷的敏感性很大,會(huì)因存在缺陷而降低.近年來(lái),在與圓柱殼氣動(dòng)彈性系統(tǒng)相似的輸流圓柱殼動(dòng)力系統(tǒng)中,也開(kāi)始注重研究初始缺陷對(duì)其運(yùn)動(dòng)穩(wěn)定性的影響,Amabili[9-10]和Prado[11]等的研究一致表明,初始撓度缺陷對(duì)系統(tǒng)的穩(wěn)定性會(huì)產(chǎn)生不利影響,使其更加容易失穩(wěn).目前國(guó)內(nèi)還未開(kāi)展這方面的研究.

綜上所述,對(duì)于含有初始撓度缺陷的各向同性完全圓柱殼的流固耦合振動(dòng)及氣動(dòng)彈性問(wèn)題,目前研究較少,對(duì)于其它類型的薄殼結(jié)構(gòu),特別是含有初始撓度缺陷的圓柱形扁殼的氣動(dòng)彈性問(wèn)題,目前尚未有研究涉及.

1　圓柱形扁殼的氣動(dòng)彈性運(yùn)動(dòng)方程

假設(shè)圓柱形扁殼的初始徑向位移(初撓度)為w0,該初始撓度并不引起應(yīng)力的變化.那么,不考慮非線性項(xiàng)時(shí),帶有初撓度的圓柱形扁殼的線性小撓度運(yùn)動(dòng)方程為[8-14]

(1)

式中:D=Eh3/12(1-ν2)為彎曲剛度;

E為楊氏模量;

ν為泊松比;

ρ為材料密度;

t為時(shí)間變量;

h為殼的厚度;

w為法向位移,以向外為正;

w0為法向初始位移；

f為Airy應(yīng)力函數(shù);

qγ為法向外力;

R為半徑;

x為軸向坐標(biāo);

θ為周向坐標(biāo)；

考慮四邊簡(jiǎn)支的圓柱形扁殼,其邊界條件為

(2)

式中:L為軸向長(zhǎng)度;

θ0為周向角；

v、u分別為周向和軸向位移.

在此將u=0和v=0寫為

(3)

式中:Nx、Nθ、Nxθ、Nθx、Mx、Mθ、Mxθ、Mθx均為內(nèi)力分量.

超聲速軸向流作用下的圓柱殼氣動(dòng)力計(jì)算采用帶有曲率修正的活塞理論公式[15](也適用于圓柱形扁殼),考慮初撓度后可以寫為[8]

(4)

式中:q=1/2ρ0U2,ρ0為來(lái)流密度;U為來(lái)流速度;Ma為馬赫數(shù).

將式(4)代入式(1)中,并引入無(wú)量綱量

(5)

得到無(wú)量綱氣動(dòng)彈性運(yùn)動(dòng)方程為

(6)

相應(yīng)無(wú)量綱邊界條件表示為

ε=0,ε=1時(shí),v=w=w0=Nε=Mε=0,即

(7a)

φ=0,φ=1時(shí),u=w=w0=Nφ=Mφ=0,即

(7b)

2　二維DQM離散系統(tǒng)運(yùn)動(dòng)方程

依據(jù)微分求積法(以下簡(jiǎn)稱DQM)二維問(wèn)題的離散化方法[16],在0≤ε≤1,0≤φ≤1的區(qū)域內(nèi),分別在ε和φ方向上置入N和M個(gè)網(wǎng)格點(diǎn),于是在平行于ε的任一直線φ=φj上,j=1,2,…,M,函數(shù)η(ε,φ)在網(wǎng)格點(diǎn)(εi,φi)處,對(duì)ε的1～4階偏導(dǎo)可近似寫為

(8)

同樣地,在任一圓弧ε=εi上,i=1,2,…,N,函數(shù)η(ε,φ)在網(wǎng)格點(diǎn)(εi,φj)處,對(duì)φ的1～4階偏導(dǎo)可近似寫為

(9)

同理,對(duì)應(yīng)力函數(shù)μ和初始撓度η0的各階偏導(dǎo)也可以寫做類似形式.其中,權(quán)系數(shù)A和B的確定參見(jiàn)文獻(xiàn)[16],在此網(wǎng)格點(diǎn)形式中取文獻(xiàn)[16]的非均勻網(wǎng)格點(diǎn).

根據(jù)式(8)、(9),可以將無(wú)量綱形式的圓柱形扁殼氣動(dòng)彈性方程式(6)寫為DQM離散后的形式，即

(10 a)

(10b)

邊界條件式(7)的DQM形式為

ε=0,1時(shí)，

(11a)

φ=0,1時(shí)，

(11b)

寫成DQM的簡(jiǎn)潔矩陣形式為

(12a)

(12b)

式中：

η0為不含邊界變量的無(wú)量綱初始撓度矩陣;

η為不含邊界變量的無(wú)量綱位移矩陣，

μ為不含邊界變量的無(wú)量綱Airy應(yīng)力函數(shù)矩陣,

(13)

具有一般矩陣形式的氣動(dòng)彈性方程可以寫成

(14)

式(14)中,因不考慮初始撓度引起的應(yīng)力,顯然含有μ0的系數(shù)項(xiàng)皆為0.引入邊界條件式(11),將含有初始撓度的矩陣η0先行代入,求得相應(yīng)的系數(shù)陣,連同代入邊界條件后含有所有位移矩陣和應(yīng)力矩陣的方程(在此不再列出)一起進(jìn)行“Kronecker”矩陣轉(zhuǎn)換.經(jīng)過(guò)“Kronecker”轉(zhuǎn)換,原矩陣形式的變量矩陣可以展開(kāi)為一維的列向量.同時(shí),將帶有初始撓度的系數(shù)陣也展開(kāi)為一維列陣,然后將其轉(zhuǎn)化為對(duì)角陣,變點(diǎn)乘為叉乘,這樣,對(duì)于“Kronecker”轉(zhuǎn)換后的一維列向量,重新構(gòu)造得到新的系數(shù)陣,新的系數(shù)陣可以寫為

(15)

式中：I為單位陣;

M為質(zhì)量陣;

C為阻尼陣;

Rd為不含邊界變量的右剛度系數(shù)陣;

Ld為不含邊界變量的左剛度系數(shù)陣;

Rb為邊界變量的右剛度系數(shù)陣;

Lb為邊界變量的左剛度系數(shù)陣;

F0為轉(zhuǎn)換后的初始撓度系數(shù)陣.

式(14)的簡(jiǎn)寫形式為

(16)

3　系統(tǒng)顫振臨界動(dòng)壓分析

一般情況下,圓柱形扁殼的初始撓度可以通過(guò)事先精確測(cè)量獲得,形成離散的測(cè)量數(shù)據(jù),直接代入式(14)～(16)中進(jìn)行變換求解得到初始撓度系數(shù)矩陣,然后進(jìn)行方程求解.根據(jù)圓柱形扁殼的幾何特征和受力特性,在此,考慮幾種特殊情況下的初始撓度,研究在這些初始撓度情況下,圓柱形扁殼線性系統(tǒng)的顫振臨界動(dòng)壓.

假設(shè)圓柱形扁殼的初始撓度形式為[12-13]

(17)

代入無(wú)量綱量

ε=x/L,η0=w0/L,φ=θ/θ0,

k0=κA0/L,κ=L/h,

則式(17)為

(18)

令k0為初始撓度缺陷因子,n0、m0分別為周向半波數(shù)和軸向半波數(shù).在給定k0、n0、m0的條件下,η0為ε、φ的確定函數(shù),因此可以通過(guò)公式求導(dǎo)得到系數(shù)矩陣.一般的情況下,η0為離散的數(shù)據(jù),可利用式(8)、(9)進(jìn)行數(shù)值求導(dǎo)得到系數(shù)矩陣.

3.1　初始撓度缺陷因子對(duì)顫振臨界動(dòng)壓的影響

曲率是影響圓柱形扁殼顫振臨界動(dòng)壓的一個(gè)重要因素.拋開(kāi)曲率因素,單純地研究初始撓度缺陷因子的影響是片面的.在此,首先選取一個(gè)較小的曲率參數(shù),來(lái)研究初始撓度缺陷因素k0、n0、m0對(duì)顫振臨界動(dòng)壓的影響情況,然后,再與其它曲率參數(shù)的結(jié)果進(jìn)行比較.取計(jì)算參數(shù)為[17]

(19)

(a) m0=1

(b) m0=-1圖1　小曲率α=20、θ0=0.05時(shí),初始撓度缺陷因子對(duì)顫振臨界動(dòng)壓的影響Fig.1　Effect of initial deflection imperfection factor on the flutter critical aerodynamic parameter at small curvature wherein α=20, θ0=0.05

半波數(shù)n0的變化.

由圖1(a)可知,當(dāng)m0=1、n0=1時(shí),隨著初始

撓度因子k0的增加,顫振臨界動(dòng)壓呈現(xiàn)增大的趨勢(shì),這種趨勢(shì)在m0=-1、n0=1時(shí)相反,說(shuō)明當(dāng)m0=1、n0=1時(shí),k0的增加實(shí)際上使得扁殼曲率的增大,曲率的增大帶來(lái)了顫振臨界動(dòng)壓的提高,而m0=-1、n0=1時(shí),反之，圖2為兩者的對(duì)比情況.

由圖2可見(jiàn),在這種特定參數(shù)下,顫振臨界動(dòng)壓參數(shù)的變化與初始撓度缺陷因子成比例關(guān)系.當(dāng)m0=1、n0>1時(shí),k0的變化對(duì)顫振臨界動(dòng)壓的影響變得復(fù)雜起來(lái).n0=2、n0=4時(shí),隨著k0的增加,顫振臨界動(dòng)壓顯著下降,下降幅度高于n0=3時(shí).當(dāng)m0=1、n0>4時(shí),k0增加,顫振臨界動(dòng)壓變化不顯著.當(dāng)m0=-1、n0>1時(shí),n0=2、n0=4時(shí),k0對(duì)顫振臨界動(dòng)壓的影響與m0=1時(shí)相同,n0=3時(shí),k0對(duì)顫振臨界動(dòng)壓的影響也非常顯著.當(dāng)m0=-1、n0>4時(shí),k0增加,顫振臨界動(dòng)壓變化不顯著.

圖2　n0=1時(shí),顫振臨界動(dòng)壓參數(shù)隨初始缺陷因子的變化Fig.2　Flutter critical aerodynamic parameter vs initial deflection imperfection factor at n0=1

圖3給出3種不同情況下無(wú)量綱頻率Ω特征值虛部隨無(wú)量綱動(dòng)壓λ變化曲線.

(a) k0=0(b) n0=2, m0=-1, k0=0.1(c) n0=2, m0=-1, k0=0.5圖3 n0=2時(shí),無(wú)量綱特征值虛部隨動(dòng)壓變化Fig.3 Imaginary parts of non-dimensional eigenvalues vs aerodynamic pressure at n0=2

由圖3可見(jiàn),隨著初始撓度缺陷系數(shù)的增加,系統(tǒng)產(chǎn)生顫振的頻率階數(shù)發(fā)生了變化,由先前的1～2階模態(tài)耦合顫振,變成了2～3階模態(tài)耦合顫振.這是由于隨著初始撓度因子的增大,耦合顫振的周向半波數(shù)發(fā)生了變化,并由此導(dǎo)致了顫振臨界動(dòng)壓的顯著變化.對(duì)于不同的n0,這種變化具有差異性.這說(shuō)明初始撓度的形式會(huì)影響產(chǎn)生顫振的周向半波數(shù).

3.2　不同曲率時(shí)初始撓度周向半波數(shù)對(duì)顫振臨界動(dòng)壓的影響

對(duì)于不同曲率的圓柱形扁殼,初始撓度缺陷對(duì)顫振臨界動(dòng)壓的影響是不同的,這尤其表現(xiàn)在不同的初始周向半波數(shù)n0上.圖4給出了α=5,θ0=0.2 圓柱形扁殼的研究結(jié)果，其它計(jì)算參數(shù)同式(19).

由圖4可見(jiàn),當(dāng)n0=1時(shí),增加初始撓度缺陷因子并沒(méi)有改變系統(tǒng)顫振臨界動(dòng)壓,這是因?yàn)榇颂幩〉那蕝?shù)較大,初始撓度因子的增加并不能顯著改變扁殼的曲率,而影響系統(tǒng)顫振臨界動(dòng)壓.當(dāng)n0>1時(shí),隨著初始缺陷因子的增大,系統(tǒng)顫振臨界動(dòng)壓顯著減小,并且當(dāng)n0>3時(shí),這種影響依然明顯,與小曲率情況不同.這是由于隨著初始撓度因子的增大,產(chǎn)生耦合顫振的周向半波數(shù)比小曲率扁殼大,導(dǎo)致了較高的初始撓度周向半波數(shù)也會(huì)顯著影響顫振臨界動(dòng)壓.

(a) m0=1(b) m0=-1圖4 大曲率α=5、θ0=0.2時(shí),初始撓度缺陷因子對(duì)顫振臨界動(dòng)壓的影響Fig.4 Effect of initial deflection imperfection factor on the critical flutter aerodynamic parameter at large curvature wherein α=5, θ0=0.2

4　結(jié)　論

(1) 當(dāng)圓柱形扁殼的曲率較小時(shí),單一的沿外法線方向的初始撓度會(huì)顯著提高系統(tǒng)的顫振臨界動(dòng)壓,而單一的沿內(nèi)法線方向的初始撓度則會(huì)降低系統(tǒng)的顫振臨界動(dòng)壓,這種影響與初始撓度因子成線性關(guān)系.而對(duì)于曲率較大的圓柱形扁殼,這種單一的初始撓度類型對(duì)系統(tǒng)顫振臨界動(dòng)壓的影響并不顯著.

(2) 當(dāng)初始撓度包含的周向半波數(shù)大于1時(shí),不同曲率的圓柱形扁殼顫振臨界動(dòng)壓受初始撓度缺陷因子的影響不同.當(dāng)曲率較小時(shí),周向半波數(shù)n0=2，3，4時(shí),隨著k0的增加,顫振臨界動(dòng)壓顯著下降,n0>4時(shí),參數(shù)k0的增加帶來(lái)的顫振臨界動(dòng)壓變化并不顯著.當(dāng)曲率較大時(shí),隨著初始缺陷因子的增大,系統(tǒng)顫振臨界動(dòng)壓顯著減少,并且當(dāng)n0>4時(shí),這種影響依然明顯.

(3) 對(duì)于特定的初始撓度缺陷類型,隨著初始撓度缺陷系數(shù)的增加,其含有的某些周向半波數(shù)會(huì)使系統(tǒng)產(chǎn)生耦合顫振的頻率階數(shù)發(fā)生變化,并由此導(dǎo)致了顫振臨界動(dòng)壓的顯著變化.這說(shuō)明初始撓度的形式會(huì)影響產(chǎn)生顫振的周向半波數(shù).這對(duì)飛行器顫振設(shè)計(jì)而言是不利的.