耐壓圓柱殼體的臨界失穩(wěn)壓力預(yù)測

于 波，劉軍鵬，張博文，徐 蒙，邱中梁，段夢蘭

（1.中國石油大學(xué)安全與海洋工程學(xué)院，北京 102249；2.中國船舶科學(xué)研究中心，江蘇無錫 214082）

0 引 言

耐壓殼體結(jié)構(gòu)是深潛器、水下空間站、水下航行器等裝備的重要組成部分，主要有圓柱殼[1]、球殼[2]以及各種組合殼[3]等型式。圓柱殼由于幾何結(jié)構(gòu)簡單，制造加工相對比較容易，已逐漸發(fā)展成為最常見的水下耐壓殼結(jié)構(gòu)，其在外壓工程環(huán)境下失效形式主要有兩種：強度破壞和屈曲失穩(wěn)，強度破壞主要是材料達(dá)到了屈服強度而發(fā)生了塑性變形，屈曲失穩(wěn)是指耐壓結(jié)構(gòu)失去穩(wěn)定性而喪失承載外壓能力的過程。在實際工程中，圓柱殼屈曲失效的概率要遠(yuǎn)遠(yuǎn)大于強度破壞，且屈曲失穩(wěn)常在無明顯征兆的情況下突然發(fā)生[1]。因此，圍繞圓柱殼體結(jié)構(gòu)的屈曲問題，進(jìn)行極限承載力的預(yù)測具有重要的科學(xué)和工程意義。

屈曲問題本質(zhì)上屬于非線性理論范疇，目前解決屈曲問題的方法主要有靜力學(xué)法、能量法以及數(shù)值試驗法，三種方法的區(qū)別和聯(lián)系列于表1。

表1 三種求解方法的區(qū)別和聯(lián)系Tab.1 Differences and connections among three solutions

基于以上三類方法，國內(nèi)外學(xué)者分別從不同的角度對圓柱殼的屈曲穩(wěn)定性問題進(jìn)行了研究，并取得了重要的研究成果。von Mises 依據(jù)小撓度理論，給出了簡支邊界條件下薄壁圓柱殼在均勻側(cè)向外壓和軸向壓力作用下的圓柱殼失穩(wěn)壓力公式：

式中，E為彈性模量，μ為泊松比，t為圓柱殼壁厚，a為半徑，l為軸向長度，n為周向失穩(wěn)波數(shù)，可通過選取合適的周向失穩(wěn)波數(shù)n來得到Pcr的最小值[4]。Nash 采用能量法詳細(xì)計算了固支邊界條件下薄壁圓柱殼在靜水外壓作用下的臨界失穩(wěn)壓力值，結(jié)果與已知結(jié)論非常接近[5]。Pinna等[6]使用能量法研究了不同邊界條件下圓柱殼在靜水外壓作用下的屈曲問題，研究發(fā)現(xiàn)，對于中等長度的圓柱殼體，不同邊界條件下的失穩(wěn)壓力可通過對殼體約束端施加標(biāo)量乘子α來確定。

在數(shù)值試驗研究方面，Aghajari等[7]進(jìn)行了圓柱形鋼殼試樣在外壓作用下的壓潰試驗，發(fā)現(xiàn)試驗結(jié)果與有限元結(jié)果較為一致。江蘇科技大學(xué)的朱永梅團(tuán)隊對不同長徑比的304不銹鋼圓柱殼進(jìn)行了外壓下的屈曲試驗，同時進(jìn)行了不同長徑比下外壓圓柱殼的數(shù)值模擬，最終發(fā)現(xiàn)實驗結(jié)果與有限元結(jié)果的偏差為2%～9%[8]。

中國船級社CCS2013[9]中給出了圓柱殼的臨界失穩(wěn)壓力經(jīng)驗公式：

美國船級社ABS2012[10]也給出了計算圓柱殼臨界失穩(wěn)壓力的經(jīng)驗公式：

式中，Py是屈服載荷，Pm是von Mises失穩(wěn)壓力。

可以發(fā)現(xiàn)，前人采用的分析方法屬于拉格朗日解算系統(tǒng)，可通過應(yīng)力法或者位移法進(jìn)行求解，許多重要的基本結(jié)果可以從線性的小撓度理論來獲得，但在有關(guān)結(jié)構(gòu)穩(wěn)定性試驗中發(fā)現(xiàn)殼體的失穩(wěn)壓力與線性理論解相差較大[11]，主要是因為理論分析時，存在大量的假設(shè)條件，且并未考慮初始幾何缺陷值以及材料的非線性等因素。經(jīng)驗公式雖然能很好地指導(dǎo)實際工程項目，但由于基于工程經(jīng)驗而引入了大量的修正系數(shù)使得臨界失穩(wěn)壓力產(chǎn)生較大的誤差。因此本文擬在減少假設(shè)條件的基礎(chǔ)上，利用靜力學(xué)方法建立圓柱殼的臨界失穩(wěn)壓力預(yù)測模型，同時建立數(shù)值模型，該模型考慮初始幾何缺陷和材料的非線性因素，最后將理論模型與數(shù)值模型分別和試驗結(jié)果及前人模型進(jìn)行比較，證明本文模型的合理性和先進(jìn)性。

1 理論模型

1.1 基本假設(shè)

前人在利用靜力法求解圓柱殼的應(yīng)變—位移幾何關(guān)系時，基于不同的假設(shè)條件得到了不同的應(yīng)變—位移關(guān)系表達(dá)式，因此，有必要對相關(guān)假設(shè)條件進(jìn)行歸納整合：

（1）中面的法線在變形過程中保持不變，滿足直法線假設(shè)；

（2）對于殼體的運動學(xué)關(guān)系，假設(shè)從中面到任一點的距離z不受殼體變形的影響，并忽略z方向的應(yīng)力σz；

（3）所有的位移都很小，滿足小變形假設(shè)；

（4）前屈曲狀態(tài)為薄膜無矩狀態(tài)，不考慮屈曲前的彎曲等因素的影響。

依據(jù)整合后的假設(shè)條件推導(dǎo)圓柱殼任意一點的應(yīng)變—位移幾何關(guān)系。

1.2 推導(dǎo)過程

根據(jù)文獻(xiàn)[12]的內(nèi)容，可得圓柱殼中面的應(yīng)變—位移關(guān)系式為

針對圓柱殼內(nèi)任一點的應(yīng)變—位移關(guān)系，在圓柱殼橫剖截面上采用微元法（圖1）進(jìn)行位移關(guān)系的分析。A為圓柱殼體半徑為r的中線上一點，其環(huán)向和徑向位移的變化分別為v和w；B為殼體內(nèi)距離中線為z的任意一點，且z＜t/2，t為圓柱殼厚度，其環(huán)向和徑向位移變化分別為vB和wB。

圖1 圓柱殼橫剖截面微元體Fig.1 Microelement body of cross section of cylindrical shell

由微元體圖可以看出，弧長vx與A點的環(huán)向位移v有相同的圓心角，計算可得

圓柱殼體變形過程中中線的偏轉(zhuǎn)角為β，同理可得

同理，在縱剖截面上（圖2）對位移關(guān)系進(jìn)行分析：

圖2 圓柱殼縱剖截面微元體Fig.2 Microelement body of longitudinal section of cylindrical shell

圖中，A點軸向和徑向位移的變化分別為u和w；B點環(huán)向和徑向位移變化分別為uB和wB，根據(jù)微元體幾何關(guān)系，有

對于B點的徑向位移wB，綜合考慮圓柱殼體縱剖面和橫剖面的相關(guān)幾何關(guān)系，則有

為了得到殼體上任一點B的應(yīng)變—位移關(guān)系，將B點的位移uB、vB和wB代替中面的位移u、v、w，同時將B點的半徑r+z代替中面的半徑r，得到圓柱殼任一點的應(yīng)變—位移關(guān)系表達(dá)式

根據(jù)胡克定律，可得

殼體的內(nèi)力和內(nèi)力矩分別為各項應(yīng)力在殼體橫截面上沿厚度方向的積分，則有

引入圓柱殼失穩(wěn)時的軸向、環(huán)向以及徑向的位移函數(shù)：

式中，U、V和W為任意常數(shù),λ=mπ/L，L為圓柱殼的長度，m為圓柱殼屈曲軸向半波數(shù)，n為環(huán)向波數(shù)。

分析位移函數(shù)可以發(fā)現(xiàn)，當(dāng)x=0 和x=L時，v=w=0,u≠0。因此可得圓柱殼的邊界條件為：兩端切向與徑向位移被限制，軸向方向沒有施加約束。

最終可得齊次線性方程組：

由此，采用靜力學(xué)方法將圓柱殼的屈曲失穩(wěn)問題轉(zhuǎn)化為本征值問題，求解圓柱殼的臨界失穩(wěn)壓力，要保證方程組有非零解，即其系數(shù)行列式為零：

式中，

賦予圓柱殼相關(guān)幾何參數(shù)，假設(shè)不同的失穩(wěn)波數(shù)m和n值，就可以得到對應(yīng)的圓柱殼的失穩(wěn)壓力，將失穩(wěn)壓力值用對數(shù)形式表示，結(jié)果如圖3 所示。在由（m，n）的不同組合求出的失穩(wěn)壓力中，選取最小的那個失穩(wěn)壓力，即為圓柱殼的臨界失穩(wěn)壓力Pcr，對應(yīng)的（m，n）值就是臨界失穩(wěn)波數(shù)，從圖3中可以看出，此圓柱殼的臨界失穩(wěn)波數(shù)是m=1，n=3，對數(shù)坐標(biāo)下結(jié)果值為-0.039，實際的臨界失穩(wěn)壓力為0.914 MPa。

圖3 不同失穩(wěn)波數(shù)下的圓柱殼失穩(wěn)壓力Fig.3 Instability pressure of cylindrical shell under different instability wave numbers

2 數(shù)值模型

分別建立不同長徑比的圓柱殼有限元模型，尺寸參數(shù)列于表2。

表2 不同長徑比的圓柱殼尺寸Tab.2 Cylindrical shell sizes with different aspect ratios

圓柱殼采用304 不銹鋼材料，E=176.05 GPa，μ=0.25，σs=323.18 MPa。分析中，設(shè)置外壓為P0=1 MPa，作用在圓柱殼整個外表面，邊界條件為：兩端周向和徑向位移被限制，軸向方向沒有約束，即

U1=U2=UR3=0，U3≠0

Riks弧長分析法可以追蹤殼體結(jié)構(gòu)的非線性平衡路徑，獲得殼體結(jié)構(gòu)的載荷—位移曲線，并且可以評估殼體的真實臨界失穩(wěn)壓力和后屈曲行為[13]。本文采用特征值屈曲—Buckle 算法以及非線性屈曲—弧長分析法進(jìn)行圓柱殼的屈曲研究。

特征值屈曲分析中，進(jìn)行了不同長徑比下的圓柱殼模態(tài)分析，將其與理論推導(dǎo)出的最小失穩(wěn)波數(shù)進(jìn)行了對比，如圖4所示?；￠L分析法中，考慮材料的非線性影響，引入由304不銹鋼真實應(yīng)力和真實應(yīng)變確定的材料塑性參數(shù)；對于初始幾何缺陷，常將線性屈曲的第1階本征模態(tài)作為最壞幾何缺陷引入。對于薄殼結(jié)構(gòu)，缺陷幅值一般取為1%～10%的殼體厚度，為了探究不同初始缺陷幅值對臨界失穩(wěn)壓力的影響規(guī)律，建立不同缺陷值下的圓柱殼模型，進(jìn)行了有限元驗證，結(jié)果如圖5 所示。其中，橫坐標(biāo)代表不同壁厚百分比下的初始缺陷幅值，縱坐標(biāo)代表圓柱殼的臨界屈曲值?？梢园l(fā)現(xiàn)，當(dāng)初始缺陷幅值增大到殼體厚度的10%后，臨界屈曲值變化趨于平緩，故缺陷幅值定為10%的殼體厚度，最終可以得到圓柱殼的平衡路徑，如圖6 所示，僅提供長徑比為1.5 情況下的結(jié)果。

圖4 不同長徑比下的失穩(wěn)波數(shù)對比Fig.4 Comparison of instability wave numbers under different aspect ratios

圖5 初始缺陷幅值對臨界失穩(wěn)壓力的影響Fig.5 Effect of initial defect amplitude on critical buckling load

圖6 圓柱殼的非線性屈曲平衡路徑Fig.6 Nonlinear buckling equilibrium path of cylindrical shell

可以看出，隨著弧長值的增加，荷載幾乎呈線性增加，直至達(dá)到臨界失穩(wěn)壓力3.401 MPa，超過該峰值，荷載急劇下降，隨后緩慢下降。同理，可以得到其他長徑比下的結(jié)果，見圖7。

圖7 理論結(jié)果與有限元結(jié)果對比圖Fig.7 Comparison between theoretical results and finite element results

由圖可知，在圓柱殼長徑比較小時，理論結(jié)果與有限元模擬結(jié)果存在一定的誤差，隨著長徑比的不斷增大，理論值與有限元結(jié)果逐漸趨于吻合。

3 案例分析與討論

為了驗證理論結(jié)果與有限元模擬結(jié)果，本文取江蘇科技大學(xué)朱永梅團(tuán)隊[8]的304不銹鋼圓柱殼的耐壓試驗作為對比驗證，E=176 050 MPa，μ=0.25，σs=323.15 MPa，試驗相關(guān)數(shù)據(jù)列于表3。表中最后一列是換算之后的圓柱殼中面半徑r中。

表3 試驗相關(guān)數(shù)據(jù)Tab.3 Relevant data of the experiment

將試驗結(jié)果分別與理論結(jié)果以及有限元結(jié)果進(jìn)行了對比，如圖8 所示。由圖8 可以看出，試驗結(jié)果在L/R=3 時存在異常，不作分析討論。三條結(jié)果曲線在L/R較小時差別較大，隨著L/R的增大，三條曲線趨于吻合，并對此進(jìn)行了誤差的定量分析。

圖8 不同長徑比下的理論、試驗與有限元結(jié)果對比Fig.8 Comparison among theoretical,experimental and finite element results under different aspect ratios

由圖9 可以看出，試驗結(jié)果與理論結(jié)果在L/R較小時，誤差值較大，且在L/R=1 時，誤差值達(dá)到30%，隨著L/R的不斷增大，誤差值逐漸減小，在L/R=4 時，誤差值降為0.6%。試驗結(jié)果與有限元模擬結(jié)果的誤差曲線較為平緩，誤差值基本都在10%以內(nèi)，試驗結(jié)果能很好地驗證有限元的模擬結(jié)果。同時可以發(fā)現(xiàn)，當(dāng)L/R值較小時，數(shù)值模擬結(jié)果更加貼近試驗結(jié)果，而當(dāng)L/R值較大時，理論模型則更具優(yōu)勢。

圖9 不同長徑比下試驗與理論、有限元結(jié)果的誤差Fig.9 Errors of experimental and theoretical and finite element results under different aspect ratios

運用靜力學(xué)方法求解圓柱殼的臨界失穩(wěn)壓力時，會涉及到求解應(yīng)變—位移幾何關(guān)系，F(xiàn)lügge[12]、Donnel[14]和Kraus[15]分別基于一定的理論及假設(shè)得到了相應(yīng)的幾何關(guān)系表達(dá)式。為了更進(jìn)一步證明理論結(jié)果的準(zhǔn)確性，本文基于不同的幾何關(guān)系表達(dá)式，得到了不同的臨界失穩(wěn)壓力，將其與理論結(jié)果進(jìn)行了比較，如圖10 所示。

由圖10 可以看出，相比于其他理論結(jié)果，本文理論模型得到的臨界失穩(wěn)壓力更為貼近試驗結(jié)果，雖然在長徑比較小時仍存在一定的誤差，但誤差值低于30%，且隨著長徑比的增大，理論解能很好地預(yù)測試驗結(jié)果。

圖10 理論結(jié)果與其他模型結(jié)果的比較Fig.10 Comparison of theoretical results with results from other models

4 結(jié) 語

本文針對圓柱殼的屈曲失穩(wěn)問題，總結(jié)了分析屈曲問題的常用方法：靜力學(xué)方法、能量法以及數(shù)值試驗法，并指出了其在分析過程中存在的局限性。本文基于整合后的假設(shè)條件，運用微元法，推導(dǎo)了圓柱殼任意一點的應(yīng)變—位移關(guān)系式，繼而應(yīng)用靜力學(xué)方法，將圓柱殼的屈曲失穩(wěn)問題化為求解方程組的特征值問題，根據(jù)穩(wěn)定性條件求解系數(shù)行列式為零的解，最終得到了圓柱殼的最小失穩(wěn)波數(shù)以及對應(yīng)的臨界失穩(wěn)壓力。

為了驗證理論推導(dǎo)出的圓柱殼臨界失穩(wěn)壓力，本文利用有限元軟件ABAQUS 對長徑比為1～7 的不同圓柱殼進(jìn)行了屈曲模擬，進(jìn)行特征值屈曲分析時，在相同的邊界條件下，得到了圓柱殼的一階屈曲模態(tài)值與理論計算的最小失穩(wěn)波數(shù)一致的結(jié)論；進(jìn)行弧長分析時，驗證了初始幾何缺陷幅值對圓柱殼臨界失穩(wěn)壓力的影響規(guī)律，根據(jù)規(guī)律將缺陷幅值取為殼體壁厚的10%，得到了圓柱殼的臨界失穩(wěn)壓力值。同時，通過對比得到臨界失穩(wěn)壓力的理論結(jié)果與有限元結(jié)果的偏差為4%～32%，且隨圓柱殼長徑比的增大而減小，有限元結(jié)果較好地驗證了圓柱殼理論屈曲解的正確性。

通過將理論計算結(jié)果和有限元模擬結(jié)果與試驗結(jié)果對比發(fā)現(xiàn)，在長徑比較小情況下數(shù)值模型結(jié)果更加貼近試驗結(jié)果，但在長徑比較大情況下，本文給出的理論模型結(jié)果則與試驗結(jié)果更加吻合。在長徑比為1的情況下，數(shù)值模型和理論模型的結(jié)果均與試驗結(jié)果存在一定的誤差，這還有待進(jìn)一步研究。