范晨光,　楊翊仁

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

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

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

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

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

1　圓柱形扁殼的氣動彈性運動方程

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

(1)

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

E為楊氏模量;

ν為泊松比;

ρ為材料密度;

t為時間變量;

h為殼的厚度;

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

w0為法向初始位移；

f為Airy應力函數;

qγ為法向外力;

R為半徑;

x為軸向坐標;

θ為周向坐標；

考慮四邊簡支的圓柱形扁殼,其邊界條件為

(2)

式中:L為軸向長度;

θ0為周向角；

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

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

(3)

式中:Nx、Nθ、Nxθ、Nθx、Mx、Mθ、Mxθ、Mθx均為內力分量.

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

(4)

式中:q=1/2ρ0U2,ρ0為來流密度;U為來流速度;Ma為馬赫數.

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

(5)

得到無量綱氣動彈性運動方程為

(6)

相應無量綱邊界條件表示為

ε=0,ε=1時,v=w=w0=Nε=Mε=0,即

(7a)

φ=0,φ=1時,u=w=w0=Nφ=Mφ=0,即

(7b)

2　二維DQM離散系統(tǒng)運動方程

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

(8)

同樣地,在任一圓弧ε=εi上,i=1,2,…,N,函數η(ε,φ)在網格點(εi,φj)處,對φ的1～4階偏導可近似寫為

(9)

同理,對應力函數μ和初始撓度η0的各階偏導也可以寫做類似形式.其中,權系數A和B的確定參見文獻[16],在此網格點形式中取文獻[16]的非均勻網格點.

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

(10 a)

(10b)

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

ε=0,1時，

(11a)

φ=0,1時，

(11b)

寫成DQM的簡潔矩陣形式為

(12a)

(12b)

式中：

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

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

μ為不含邊界變量的無量綱Airy應力函數矩陣,

(13)

具有一般矩陣形式的氣動彈性方程可以寫成

(14)

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

(15)

式中：I為單位陣;

M為質量陣;

C為阻尼陣;

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

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

Rb為邊界變量的右剛度系數陣;

Lb為邊界變量的左剛度系數陣;

F0為轉換后的初始撓度系數陣.

式(14)的簡寫形式為

(16)

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

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

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

(17)

代入無量綱量

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

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

則式(17)為

(18)

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

3.1　初始撓度缺陷因子對顫振臨界動壓的影響

曲率是影響圓柱形扁殼顫振臨界動壓的一個重要因素.拋開曲率因素,單純地研究初始撓度缺陷因子的影響是片面的.在此,首先選取一個較小的曲率參數,來研究初始撓度缺陷因素k0、n0、m0對顫振臨界動壓的影響情況,然后,再與其它曲率參數的結果進行比較.取計算參數為[17]

(19)

(a) m0=1

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

半波數n0的變化.

由圖1(a)可知,當m0=1、n0=1時,隨著初始

撓度因子k0的增加,顫振臨界動壓呈現增大的趨勢,這種趨勢在m0=-1、n0=1時相反,說明當m0=1、n0=1時,k0的增加實際上使得扁殼曲率的增大,曲率的增大帶來了顫振臨界動壓的提高,而m0=-1、n0=1時,反之，圖2為兩者的對比情況.

由圖2可見,在這種特定參數下,顫振臨界動壓參數的變化與初始撓度缺陷因子成比例關系.當m0=1、n0>1時,k0的變化對顫振臨界動壓的影響變得復雜起來.n0=2、n0=4時,隨著k0的增加,顫振臨界動壓顯著下降,下降幅度高于n0=3時.當m0=1、n0>4時,k0增加,顫振臨界動壓變化不顯著.當m0=-1、n0>1時,n0=2、n0=4時,k0對顫振臨界動壓的影響與m0=1時相同,n0=3時,k0對顫振臨界動壓的影響也非常顯著.當m0=-1、n0>4時,k0增加,顫振臨界動壓變化不顯著.

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

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

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

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

3.2　不同曲率時初始撓度周向半波數對顫振臨界動壓的影響

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

由圖4可見,當n0=1時,增加初始撓度缺陷因子并沒有改變系統(tǒng)顫振臨界動壓,這是因為此處所取的曲率參數較大,初始撓度因子的增加并不能顯著改變扁殼的曲率,而影響系統(tǒng)顫振臨界動壓.當n0>1時,隨著初始缺陷因子的增大,系統(tǒng)顫振臨界動壓顯著減小,并且當n0>3時,這種影響依然明顯,與小曲率情況不同.這是由于隨著初始撓度因子的增大,產生耦合顫振的周向半波數比小曲率扁殼大,導致了較高的初始撓度周向半波數也會顯著影響顫振臨界動壓.

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

4　結　論

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

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

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