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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

1 理論模型

1.1 基本假設

前人在利用靜力法求解圓柱殼的應變—位移幾何關系時，基于不同的假設條件得到了不同的應變—位移關系表達式，因此，有必要對相關假設條件進行歸納整合：

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

（2）對于殼體的運動學關系，假設從中面到任一點的距離z不受殼體變形的影響，并忽略z方向的應力σz；

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

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

依據(jù)整合后的假設條件推導圓柱殼任意一點的應變—位移幾何關系。

1.2 推導過程

根據(jù)文獻[12]的內容，可得圓柱殼中面的應變—位移關系式為

針對圓柱殼內任一點的應變—位移關系，在圓柱殼橫剖截面上采用微元法（圖1）進行位移關系的分析。A為圓柱殼體半徑為r的中線上一點，其環(huán)向和徑向位移的變化分別為v和w；B為殼體內距離中線為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有相同的圓心角，計算可得

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

同理，在縱剖截面上（圖2）對位移關系進行分析：

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

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

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

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

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

殼體的內力和內力矩分別為各項應力在殼體橫截面上沿厚度方向的積分，則有

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

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

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

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

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

式中，

賦予圓柱殼相關幾何參數(shù)，假設不同的失穩(wěn)波數(shù)m和n值，就可以得到對應的圓柱殼的失穩(wěn)壓力，將失穩(wěn)壓力值用對數(shù)形式表示，結果如圖3 所示。在由（m，n）的不同組合求出的失穩(wěn)壓力中，選取最小的那個失穩(wěn)壓力，即為圓柱殼的臨界失穩(wěn)壓力Pcr，對應的（m，n）值就是臨界失穩(wěn)波數(shù)，從圖3中可以看出，此圓柱殼的臨界失穩(wěn)波數(shù)是m=1，n=3，對數(shù)坐標下結果值為-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。分析中，設置外壓為P0=1 MPa，作用在圓柱殼整個外表面，邊界條件為：兩端周向和徑向位移被限制，軸向方向沒有約束，即

U1=U2=UR3=0，U3≠0

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

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

圖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

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

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

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

3 案例分析與討論

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

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

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

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

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

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

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

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

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

4 結 語

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

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

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