劉明峰, 徐 典, 倪卓凡, 李逸豪, 李 銳
(大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析優(yōu)化與CAE軟件全國(guó)重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024)
圓柱殼以其輕質(zhì)屬性和優(yōu)異的力學(xué)性能,被廣泛應(yīng)用于航空航天[1-2]、船舶工程[3]、土建工程[4]等領(lǐng)域,其結(jié)構(gòu)失穩(wěn)可能導(dǎo)致嚴(yán)重的安全問(wèn)題,這使得屈曲分析成為圓柱殼結(jié)構(gòu)設(shè)計(jì)中的必要環(huán)節(jié).軸向載荷作用下的屈曲是圓柱殼最重要的失效模式之一,許多學(xué)者采用各類(lèi)數(shù)值方法對(duì)其屈曲問(wèn)題開(kāi)展了深入研究,例如有限單元法[5]、有限差分法[6]、Galerkin方法[7-9]、Rayleigh-Ritz法[10]、微分求積法[11]、無(wú)網(wǎng)格法[12]等.
在數(shù)值方法興盛的同時(shí),對(duì)解析解的研究仍然具有重要意義:解析解不僅可以作為檢驗(yàn)各種數(shù)值方法的對(duì)比基準(zhǔn),還能夠快速揭示不同參量對(duì)結(jié)構(gòu)力學(xué)行為的影響,從而指導(dǎo)結(jié)構(gòu)的高效設(shè)計(jì).然而,對(duì)于圓柱殼屈曲問(wèn)題的解析解研究并不多見(jiàn),例如:Zou和Foster[13]基于Flügge的各向同性圓柱殼線(xiàn)性理論,推導(dǎo)了軸壓和外壓共同作用下圓柱殼屈曲的一般解,并應(yīng)用于正交各向異性圓柱殼的求解;張俊霖等[14]使用辛方法開(kāi)展了吸濕老化影響下天然纖維增強(qiáng)圓柱殼的屈曲分析;桂夷斐和馬建敏[15]根據(jù)Hamilton變分原理求解了彈性介質(zhì)中受軸向沖擊載荷作用的圓柱殼的屈曲;Sun等[16]發(fā)展了一種等效硬殼理論,用以預(yù)測(cè)準(zhǔn)各向同性?shī)A層圓柱殼的力學(xué)行為;Chen等[17]利用狀態(tài)空間法分析了受壓功能梯度彈性空心圓柱的分叉屈曲問(wèn)題.已有研究大多聚焦封閉圓柱殼的屈曲分析,對(duì)于工程中亦常見(jiàn)的開(kāi)口圓柱殼的屈曲問(wèn)題則少有討論,而板殼力學(xué)中的經(jīng)典解析解法——Navier解法[18]和Lévy解法[19]又僅適用于四邊簡(jiǎn)支和對(duì)邊簡(jiǎn)支情況(即Lévy型開(kāi)口圓柱殼),而對(duì)于非Lévy型情況則難以處理.
本文基于筆者近年來(lái)針對(duì)板殼力學(xué)問(wèn)題提出的一種新解析方法:辛疊加方法[20-22],獲得了經(jīng)典解法難以直接求解的典型非Lévy型正交各向異性開(kāi)口圓柱殼屈曲問(wèn)題的解析解.該方法直接從正交各向異性圓柱殼的基本方程出發(fā),將問(wèn)題導(dǎo)入到Hamilton體系中,然后將非Lévy型邊界下的原問(wèn)題拆分為兩個(gè)子問(wèn)題,通過(guò)分離變量和辛本征展開(kāi)求出子問(wèn)題的解,進(jìn)一步通過(guò)疊加求出原問(wèn)題的解.在基于辛疊加方法的求解過(guò)程中,不需要預(yù)先選取試函數(shù)(如位移函數(shù)),而是通過(guò)逐步理性推導(dǎo),嚴(yán)格得到問(wèn)題的解,因而突破了半逆法的限制,可以求解傳統(tǒng)框架下不易處理的殼體問(wèn)題.通過(guò)數(shù)值算例驗(yàn)證了本文求解的正確性,并進(jìn)一步開(kāi)展了定量的參數(shù)分析,研究結(jié)果可以為開(kāi)口圓柱殼的屈曲設(shè)計(jì)提供理論參考.
如圖1所示,開(kāi)口圓柱殼承受軸向力N1的作用,柱殼的厚度為δ,曲率半徑為R,軸向(α方向)的彈性模量為E1、Poisson比為ν21,環(huán)向(β方向)的彈性模量為E2,Poisson比為ν12,剪切模量為G12.柱殼沿軸向、環(huán)向和徑向的位移函數(shù)分別記為u,v和w.
根據(jù)Kirchhoff-Love假設(shè),殼體內(nèi)部任意一點(diǎn)的應(yīng)變分量可以表示為
e={eα,eβ,eαβ}T=ε+γχ,
(1)
其中ε={εα,εβ,εαβ}T={?u/?α,?v/?β+w/R,?u/?β+?v/?α}T,εα和εβ為中面內(nèi)各點(diǎn)沿α和β方向的線(xiàn)應(yīng)變,εαβ為面內(nèi)切應(yīng)變;γ為殼體中面法線(xiàn)方向的坐標(biāo);χ={χα,χβ,χαβ}T={-?2w/?α2,-?2w/?β2,-?2w/?α?β}T,χα和χβ為中面內(nèi)各點(diǎn)主曲率的改變,χαβ為扭率的改變.
圖1 開(kāi)口圓柱殼示意圖Fig. 1 Schematic diagram of a cylindrical shell
對(duì)于正交各向異性圓柱殼,物理方程為
(2)
其中
b11=E1/(1-ν21ν12),b12=E1ν21/(1-ν21ν12),b22=E2/(1-ν21ν12),b66=G12.
進(jìn)一步,基于Donnell薄殼理論[23],殼體所受內(nèi)力和彎矩可以表示為
(3)
其中FT1和FT2分別為垂直于α和β軸的截面內(nèi)單位寬度上的拉壓力,FT12和FT21分別為相應(yīng)截面內(nèi)的平錯(cuò)力,M1和M2分別為彎矩,M12和M21分別為扭矩.
等效剪力V1,V2可以表示為
(4)
其中FS1和FS2分別為垂直于α和β軸的截面內(nèi)的橫向剪力.
殼體的平衡方程可以表示為
(5)
利用1.1小節(jié)所示基本方程,可將開(kāi)口圓柱殼屈曲問(wèn)題導(dǎo)入Hamilton體系.定義
(6)
由式(1)—(6)可以得到
(7)
(8)
(9)
其中
(10)
采用類(lèi)似的方法,定義?w/?α=θ1,對(duì)于子問(wèn)題2,Hamilton對(duì)偶方程為
(11)
(12)
其中
(13)
這里
在辛幾何求解框架下,可以采用分離變量[25]使得
Z=A(α)B(β),
(14)
(15)
其中μ為本征值,A(α)為本征向量.
式(15)中第二式的通解可以寫(xiě)為
(16)
其中λi是微分方程的特征值,待定系數(shù)Cji和Dji(j=1,2,…,8;i=1,2,3,4)并不相互獨(dú)立.將式(16)代回到式(15)的第二式,可以得到系數(shù)之間的關(guān)系:
(17)
其中
考慮如下邊界條件的開(kāi)口圓柱殼:
(18)
如圖2(a)所示的原問(wèn)題邊界條件為:α=0和α=a邊固支,用“C”表示,β=0和β=b邊簡(jiǎn)支,但并不滿(mǎn)足Lévy型邊界條件[23],用“S′”表示.將原問(wèn)題拆分成兩個(gè)子問(wèn)題:子問(wèn)題1是四邊簡(jiǎn)支的開(kāi)口圓柱殼施加了非零的軸向位移u|β=0,b≠0,子問(wèn)題2是四邊簡(jiǎn)支的開(kāi)口圓柱殼施加了非零的環(huán)向位移v|α=0,a≠0和彎矩M1|α=0,a≠0.兩個(gè)子問(wèn)題中的Lévy型簡(jiǎn)支邊用“S”表示.
將式(16)代入α方向的對(duì)邊簡(jiǎn)支邊界條件,即
FT1|α=0,a=0,v|α=0,a=0,w|α=0,a=0,M1|α=0,a=0,
(19)
可以得到
sinh(aλ1)sinh(aλ2)sinh(aλ3)sinh(aλ4)=0.
(20)
求解μI-H的行列式展開(kāi)可以得到重根零本征值μ00=0和本征值μim(i=1,2,…,8;m=1,2,3,…),對(duì)應(yīng)的本征向量為
A00=[1,0,0,0,0,0,0,0]T,
(21)
A01=[0,0,0,0,b66δ,0,0,0]T,
(22)
以及
(23)
其中, 式(21)和式(22)分別由AA00=0和HA01=A00得到,am=mπ/a(m=1,2,3,…),uim,vim,wim,Kim分別為
(24)
齊次方程(8)的解可以寫(xiě)為
Z=[A00,A01,A11(α),A21(α),…,A81(α),…,A1m(α),A2m(α),…,A8m(α),…]B(β),
(25)
其中B(β)=[B00(β),B01(β),B11(β),B21(β),…,B81(β),…,B1m(β),B2m(β),…,B8m(β),…]T,由式(15)中的第一式可以得到
B01(β)=c01,B00(β)=c01β+c00,Bim(β)=cimeμimβ.
(26)
因此,子問(wèn)題1的位移解可以寫(xiě)為
(27)
其中c00,c01和cim為待定系數(shù).
(a) 原問(wèn)題 (b) 子問(wèn)題1 (c) 子問(wèn)題2 (a) The original problem (b) Subproblem 1 (c) Subproblem 2
在β=0和β=b處施加級(jí)數(shù)形式的位移,相關(guān)的邊界條件可以表示為
(28)
將式(27)代入式(28),待定系數(shù)cim得以由u0m和ubm表示.
對(duì)于子問(wèn)題2,代入相應(yīng)的邊界條件,同理可以得出位移解:
(29)
在α=0和α=a處施加級(jí)數(shù)形式的位移,相應(yīng)的邊界條件表示為
(30)
將式(27)和式(29)中兩個(gè)子問(wèn)題的解進(jìn)行疊加,為了與式(18)所對(duì)應(yīng)的原問(wèn)題的解等價(jià),還需要滿(mǎn)足的邊界條件為
(31)
將式(27)和式(29)代入式(31),利用三角級(jí)數(shù)展開(kāi),具體地,由FT12+(M12/R)|α=0=0可以得到
(32)
(33)
其中n=1,2,3,….
由FT12+(M12/R)|α=a=0可以得到式(32)和
(34)
其中n=1,2,3,….
由?w/?α|α=0,a=0可以得到
(35)
和
(36)
其中n=1,2,3,….
由FT21|β=0,b=0可以得到
(37)
(38)
以及
(39)
其中m=1,2,3,….
下面給出開(kāi)口圓柱殼屈曲問(wèn)題的數(shù)值算例.選取由石墨和環(huán)氧樹(shù)脂組成的正交各向異性復(fù)合材料[26],其材料屬性為E1=128 GPa,E2=11 GPa,ν12=0.25,G12=4.48 GPa.
表1為屈曲載荷解的收斂性研究,所研究的開(kāi)口圓柱殼縱向尺寸與環(huán)向尺寸以及半徑相同(a=b=R=1 m),δ分別為0.01 m和0.005 m.顯然,在兩種厚度下,所有計(jì)算結(jié)果在級(jí)數(shù)取20項(xiàng)時(shí)均收斂到5位有效數(shù)字(表中粗體數(shù)據(jù)),表明本文方法得到的解析解收斂性非常好.由于文獻(xiàn)中缺乏可供對(duì)比的解析解,因此將本文的計(jì)算結(jié)果與有限元軟件ABAQUS得到的結(jié)果進(jìn)行對(duì)比,其中有限元的參數(shù)按照如下設(shè)置:采用S4R單元,網(wǎng)格大小為殼最小面內(nèi)尺寸的1/100.表2和表3分別展示了不同尺寸下開(kāi)口圓柱殼的屈曲載荷.當(dāng)δ=0.01 m時(shí),本文結(jié)果與有限元的差別均在5%以?xún)?nèi);當(dāng)δ=0.005 m時(shí),與有限元的差別均在1%以?xún)?nèi),驗(yàn)證了本文解析結(jié)果的正確性.除得到屈曲載荷之外,本文方法還得到了開(kāi)口圓柱殼的屈曲模態(tài),包括高階模態(tài).以a=b=R=1 m,δ=0.005 m的圓柱殼為例,圖3給出了其前十階屈曲模態(tài)與有限元結(jié)果的對(duì)比,結(jié)果高度吻合.
表1 a=b=R=1 m時(shí),開(kāi)口圓柱殼前十階屈曲載荷的收斂性研究(單位: kN/m)
圖3 開(kāi)口圓柱殼的前十階屈曲模態(tài)Fig. 3 The 1st 10 buckling modes of the cylindrical shell
表2 b=1 m, δ=0.01 m時(shí),開(kāi)口圓柱殼的前十階屈曲載荷(單位: kN/m)
表3 b=1 m, δ=0.005 m時(shí),開(kāi)口圓柱殼的前十階屈曲載荷(單位: kN/m)
圖4 不同長(zhǎng)度下,開(kāi)口圓柱殼臨界屈曲載荷隨厚度的變化曲線(xiàn)
利用所得解析解進(jìn)一步考查主要參數(shù)對(duì)開(kāi)口圓柱殼屈曲結(jié)果的影響.圖4給出了在不同長(zhǎng)度下厚度對(duì)臨界屈曲載荷的影響(取b=R=1 m).可以看到,隨著厚度增加臨界載荷均顯著增大,而且在殼的長(zhǎng)度較短時(shí),這種效應(yīng)更加明顯.圖5給出了兩種厚度下曲率1/R對(duì)臨界載荷的影響(取a=b=1 m),在曲率較小時(shí),殼的形態(tài)接近于平板,隨著曲率的增大,對(duì)軸向剛度的貢獻(xiàn)加強(qiáng),因而提高了殼的承載能力.上述結(jié)果可以為開(kāi)口圓柱殼的抗屈曲設(shè)計(jì)提供參考.
由于高階偏微分控制方程的求解困難,開(kāi)口圓柱殼的解析求解本就是一項(xiàng)具有挑戰(zhàn)性的課題,而對(duì)于非Lévy邊界的正交各向異性開(kāi)口圓柱殼的屈曲問(wèn)題,求解難度無(wú)疑會(huì)進(jìn)一步增加.本文從Donnell殼體理論的基本方程出發(fā),導(dǎo)出了正交各向異性開(kāi)口圓柱殼屈曲的Hamilton對(duì)偶方程,利用辛疊加方法,將原問(wèn)題拆分,采用分離變量、辛本征展開(kāi)求解得到子問(wèn)題的解,最后根據(jù)疊加后的解與原問(wèn)題解的等價(jià)性要求得到關(guān)于待定常數(shù)的齊次線(xiàn)性方程組,根據(jù)方程組有非零解的條件得到屈曲載荷解,進(jìn)而得到相應(yīng)的屈曲模態(tài).由于文獻(xiàn)中缺乏對(duì)應(yīng)的解析解,因此采用有限元數(shù)值解對(duì)本文結(jié)果進(jìn)行了對(duì)比驗(yàn)證,同時(shí)還利用本文的解研究了殼的尺寸對(duì)屈曲結(jié)果的影響.辛疊加方法在求解過(guò)程中無(wú)須假定試函數(shù),兼?zhèn)淞诵练椒ǖ睦硇院童B加法的規(guī)律性,有望推廣到更多殼體問(wèn)題的求解當(dāng)中.