王偉，朱時洋，和衛(wèi)平

1海軍駐葫蘆島地區(qū)軍事代表室，遼寧葫蘆島125004

2武漢第二船舶設計研究所，湖北武漢430205

0 引 言

在工程結構中，圓柱殼得到廣泛應用，但其對缺陷具有明顯的敏感性，使屈曲載荷存在很大的折減和較大范圍的離散分布。在工程實際中，圓柱殼結構在生產制作、運輸?shù)冗^程中不可避免地會產生初始幾何缺陷［1］，且初始幾何缺陷對圓柱殼結構的承載能力有較大的影響，所以基于實測幾何缺陷分析圓柱殼的承載能力非常有必要。

國外很早就開展了針對此問題的相關研究。Southwell［2］通過對理論與實驗結果的對比，觀察到實測的屈曲載荷總是比對應模型的理論值小很多的物理規(guī)律。Von Karman 和 Tsien［3］通過求解圓柱殼的非線性大撓度方程，提出了更加普適的圓柱殼后屈曲求解方法。Koiter［4］首次提出了采用攝動理論研究分析線彈性結構初始后屈曲路徑的方法。Donnell和 Wan［5］將圓柱殼幾何缺陷引入到了后屈曲分析中，計算獲得的臨界屈曲載荷比完美/完善殼體的結果要低很多。在國內，有學者針對薄壁圓柱殼結構的屈曲問題，開展了大撓度方程求解方法［6］和薄壁圓柱殼彈塑性屈曲［7-8］的理論研究。上述研究表明，在殼體存在初始缺陷的情況下，圓柱殼的極限承載能力與實驗結果差異明顯，且實驗結果具有較大的隨機離散性。殼體的初始缺陷包括初始幾何缺陷、殼體厚度缺陷、加載載荷缺陷等，但是由于真實缺陷的種類有很多，已經超出了經典幾何缺陷的范疇，而非經典幾何缺陷則導致影響結構臨界屈曲載荷的因素更復雜多變，極大增加了殼體結構設計的難度。其他常見的缺陷還包括線性屈曲模態(tài)型缺陷（LBMI）、軸對稱缺陷（ASI）和實測缺陷（MI）。其中，線性屈曲模態(tài)型缺陷被廣泛應用于土木工程領域，用于對考慮了風載荷作用的油罐、軸壓作用的冷卻塔以及運載火箭等薄壁殼體結構進行屈曲分析。此外，線性屈曲模態(tài)型缺陷還被應用于減縮剛度分析。線性屈曲模態(tài)型缺陷被廣泛應用的原因之一在于借助通用的有限元軟件對結構進行線性屈曲分析，可以很方便地獲取該形式的缺陷。軸對稱缺陷現(xiàn)在更多地被用于開展缺陷敏感性對比分析。實測缺陷主要是通過接觸或者非接觸測量方法獲取殼體結構的表面形貌點云數(shù)據(jù)，并將該缺陷引入到有限元分析中。在實測缺陷研究領域，Arbocz和 Abramovich［9］做出了卓越貢獻，建立了初始缺陷數(shù)據(jù)庫。

國內對于在嚴酷運行環(huán)境下水下環(huán)肋圓柱殼結構的非線性屈曲理論研究和極限承載能力分析方面起步較晚，不僅分析手段單一，精度也不足，且試驗數(shù)據(jù)相對匱乏，在計及實測初始幾何缺陷的結構承載能力分析方面尤其如此。因此，為了測量某圓柱殼體的初始幾何缺陷，本文擬首先基于二次型變換方法提取殼體的初始幾何缺陷，然后利用雙重傅里葉級數(shù)方法重構有限元模型，分析初始幾何缺陷對圓柱殼外壓承載能力的影響，最后對比分析2種有限元模型重構方法在計算精度和效率方面存在的差異。

1 幾何缺陷的測量與表達

本文使用三維形貌測量系統(tǒng)獲取了試件表面的三維形貌結果，圖1所示為其點云形式。但是，在實際使用這些點云數(shù)據(jù)的過程中，點云數(shù)據(jù)所在的測量坐標系與試件仿真建模的坐標系并不一致，因此，還需要經過坐標變換處理以獲取期望的空間位置，如圖2所示。

圖1 三維形貌測量系統(tǒng)輸出的三維形貌數(shù)據(jù)Fig.1 3D topography data acquired from 3D measurement system

圖2 三維形貌數(shù)據(jù)的期望空間位置Fig.2 Expected spatial position of 3D topography data

雖然基于三維形貌測量系統(tǒng)可以獲取圓柱殼表面的整體形貌點云信息，但卻無法直接獲取圓柱殼半徑、軸線方向等特征參數(shù)。因此，本文采用最小二乘法先將三維形貌測量系統(tǒng)輸出的圓柱殼表面整體形貌點云信息進行擬合，得到二次型曲面的一般方程，然后通過正交變換，將點云所在的測量坐標系轉換到有限元建模采用的標準參考坐標系中。此外，通過對二次型系數(shù)的分析，還可以直接獲取圓柱殼直徑的擬合值。

圖3所示為本文采用的坐標系。其中，O-XYZ為測量坐標系，O′-X′Y′Z′為標準坐標系（建模坐標系），O′-rθy（r，θ，y分別表示徑向、周向和軸坐標）為標準圓柱殼坐標系。圖中，A為自定義的標定點，本文取圓柱殼外表面第1象限與下端面的交點。

圖3 本文采用的坐標系定義Fig.3 Definition of the coordinate systems in this paper

在測量坐標系中，圓柱殼可以采用一般形式的三元二次曲面方程描述，即

其矩陣形式寫為

其中，

式中:A為二次型方程特征系數(shù)的對稱矩陣；X為自變量矩陣；b為一次項系數(shù)向量；a11～a33為系數(shù)；a44為常數(shù)項。

通過坐標值的交叉計算，可以將求解式（1）所示三元二次曲面方程轉化為多元線性擬合的問題。在擬合過程中，使用最小二乘法求解得到二次型方程的系數(shù)矩陣。由于A為實對稱矩陣，所以存在正交矩陣R0。

式中，Λ為對角矩陣。

進一步，式（2）可以變換為標準坐標系下的式（4）

其中，

式中，c為一次項系數(shù)向量。

通過正交變換［10］，將二次曲面的一般方程轉換為標準方程。該變換消除了坐標分量間的交叉項，通過校正附加的特征點，將測量坐標系與標準坐標系的方向對齊，同時保持曲面形狀不變。

將標準坐標系中的點云轉換到標準圓柱殼坐標系下，計算形貌點云中各散點與擬合曲面間的距離（離面位移w），如圖4所示。該離面位移場w(x，y)包含了殼體的初始幾何缺陷信息，也是后續(xù)研究的數(shù)據(jù)基礎。圖中：θ為周向角度；為歸一化的無量綱離面位移（(x，y)=w(x，y)/t；t為圓柱殼壁厚）；L為圓柱殼名義長度。

圖4 圓柱殼坐標系下初始幾何缺陷的散點圖Fig.4 Scattered points diagram of initial geometric defects in cylinder coordinate system

圓柱殼結構的真實缺陷對結構承載能力有著極大的影響，是多年來研究的熱點和難點問題。國內外學者通常將圓柱殼的真實缺陷（離面位移場w(x，y)）展開為半波余弦或半波正弦的雙重傅里葉級數(shù)。

在標準圓柱殼坐標系下，半波余弦和半波正弦的雙重傅里葉級數(shù)分別表示為：

式中：R為圓柱殼的名義半徑；Akl，Bkl為半波余弦傅里葉級數(shù)的系數(shù)；Ckl，Dkl為半波正弦傅里葉級數(shù)的系數(shù)；k，l分別表示軸向半波數(shù)和周向全波數(shù)；n1，n2為雙重傅里葉級數(shù)項數(shù)。

2 圓柱殼軸壓承載能力分析

為了獲取圓柱殼的實際屈曲載荷大小，需要將實測幾何缺陷引入到有限元模型中，以便將完美模型重構為含實測幾何缺陷的有限元模型。有限元模型重構方法包括了基于散點的有限元模型修正方法和基于傅里葉級數(shù)函數(shù)的有限元模型重構方法。后者在本文中稱為傅里葉級數(shù)法，即以傅里葉級數(shù)的形式引入實測缺陷。具體講，就是將實測幾何缺陷采用傅里葉級數(shù)形式表達后，再將該傅里葉級數(shù)引入到有限元模型中修改節(jié)點坐標，從而引入實測幾何缺陷。

首先，本文將如圖4所示的圓柱殼坐標系下的點云信息通過逆距離加權插值法進行插值，獲得了如圖5所示的格柵化離面位移圖，從而將離面位移場表示為二維矩陣的形式。

圖5 圓柱殼格柵化的離面位移圖Fig.5 Off-plane displacement diagram in grille form of cylindrical shell

然后，將該二維矩陣形式進行周期延拓，使用高斯濾波器對該二維矩陣進行濾波處理，以降低測量噪聲的干擾，得到如圖6所示高斯降噪后的離面位移結果。圖中，黑點為原始的散點數(shù)據(jù)。由圖6可見，經過格柵化及降噪處理后的離面位移場很好地反映了原始散點所表示的實測形貌。

圖6 圓柱殼高斯降噪后的離面位移圖Fig.6 Off-plane displacement diagram of cylindrical shell after Gaussian noise reduction

最后，將求得的半波余弦傅里葉級數(shù)的系數(shù)Akl，Bkl以文件形式保存，編寫Python腳本讀取該系數(shù)，并在ABAQUS有限元模型中使用該系數(shù)表示的傅里葉級數(shù)引入實測缺陷，實現(xiàn)含實測缺陷的有限元模型重構。圖7所示為重構形成的有限元模型，考慮到顯示效果，實測缺陷的幅值被放大了30倍。

圖7 含實測缺陷的有限元圓柱殼重構模型Fig.7 FE model of reconfigured cylindrical shell with measured defects

2.1 完美模型的承載能力數(shù)值分析

本文所示完美模型是指不含實測缺陷且經過完善的圓柱殼模型。該完美圓柱殼模型的名義半徑R是通過對實際結構的實測形貌經過特征提取處理獲得的擬合半徑。圖8所示為模型的邊界條件及載荷。

圖8中模型底端為固支，即約束模型底端框上節(jié)點的3個線位移、3個轉角自由度。上端僅放松軸向位移自由度，約束其他5個自由度。模型外表面施加均布壓力載荷，端面施加相當?shù)妮S向載荷，并考慮結構的非線性變形。通過網格收斂性分析，合理設置網格尺寸。材料屬性采用雙線性彈塑性模型。

圖8 有限元圓柱殼模型邊界條件Fig.8 Boundary conditions of cylindrical shell FE model

計算結果顯示，該完美圓柱殼有限元模型的極限承載能力Pperfect=0.835 7（基于理論極限承載能力歸一化的無量綱壓力載荷，下同），其位移—載荷曲線如圖9所示。由圖9可見，OA線段為圓柱殼的線性前屈曲階段，其承載能力隨著位移載荷的增加而呈線性增大，在A點（臨界失穩(wěn)狀態(tài)）處達到最大值后發(fā)生局部殼板失穩(wěn)；隨著施加位移載荷的繼續(xù)增加，承載能力逐漸下降并在B點（壓潰后某一大變形狀態(tài)）逐漸穩(wěn)定。

2.2 模態(tài)型缺陷的近似分析

圖9 完美圓柱殼有限元模型的位移—載荷曲線Fig.9 Displacement-loading curve of perfect cylindrical shell FE model

由于模態(tài)或形如模態(tài)的初始幾何缺陷對結構的軸壓臨界屈曲載荷有著巨大影響，因此模態(tài)型缺陷近似方法也是圓柱殼結構承載能力預測分析時所采用的重要方法之一。該方法主要針對初始缺陷形式并不十分清楚的結構，并通常使用結構的一階模態(tài)來模擬結構的初始幾何缺陷。由于該方法的計算結果不一定與實驗數(shù)據(jù)吻合得很好，故一般用于計算結構的下限承載能力。

這里，使用2.1節(jié)的圓柱殼有限元模型對模型施加歸一化位移載荷，并建立線性屈曲分析步。為了引入屈曲模態(tài)缺陷，修改線性屈曲分析的INP文件，以增加節(jié)點位移向量的輸出。完成線性屈曲分析后，新建Riks分析步，并在對應的INP文件中增加引入線性屈曲分析結果的命令，然后開始極限承載能力的計算。相關流程及命令如下：

1）進行線性屈曲分析前，在INP文件中添加以下命令。計算完成后，工作目錄下會生成一個擴展名為fil的文件，該文件將會在下一步被用到。

2）使用以下命令引入上一步的模態(tài)型缺陷文件：

其中，F(xiàn)ileName為上一步fil文件的文件名（不含擴展名）。StepNumber必須與上一步模型中的相同。ScalingFactor一般為1。

3）施加軸外壓載荷，并進行計算分析。

分別引入如圖10所示的一階局部模態(tài)型缺陷和一階總體模態(tài)型缺陷，缺陷幅值分別為0.2t和0.002 5R，得到如圖11所示圓柱殼位移—載荷曲線。

圖10 圓柱殼的一階模態(tài)型缺陷Fig.10 The first-order mode defects of cylindrical shell

圖11 考慮模態(tài)型缺陷時圓柱殼有限元模型的位移—載荷曲線Fig.11 Displacement-loading curves of cylindrical shell FE model considering mode defects

由圖11可見，OA段處于線性前屈曲階段，圓柱殼的承載能力與位移載荷基本呈線性關系；隨著位移載荷的增加，圓柱殼進入非線性后屈曲階段（AB段），此階段圓柱殼的承載能力與位移載荷也基本呈線性關系，不過結構的靜力學剛度明顯下降；隨著廣義位移的增加，圓柱殼的最大承載能力 Plocal，Pglobal分別降為0.812 5和0.685 3，隨后殼板發(fā)生失穩(wěn)，且這2種模態(tài)型缺陷導致的殼板變形存在一定的差異。

2.3 實測缺陷近似分析

由于實測缺陷真實地反映了圓柱殼結構本身的初始幾何缺陷，所以基于結構的實測缺陷開展極限承載能力分析能夠非常精確地預測其極限承載能力。本文使用的方法為實測缺陷近似方法（Measured Imperfection Approach，MIA）。

這里，同樣使用與2.1節(jié)相同的有限元模型。采用面向大型圓柱殼結構的形貌測量方法，對圓柱殼進行了形貌測量、分析及表達，并通過基于散點或傅里葉級數(shù)函數(shù)的有限元模型修正方法引入實測缺陷，獲得了前文圖7所示的含實測缺陷的圓柱殼有限元模型。

使用散點法和傅里葉級數(shù)法分別實現(xiàn)圓柱殼有限元模型極限承載能力分析后的結果顯示，在考慮實測缺陷后，得到如圖12所示采用2種方法得出的圓柱殼極限承載能力（PPoint和PFourier）分別為0.824 2和0.823 8。由圖12可見，散點法和傅里葉級數(shù)法計算的結果幾乎完全吻合。在OA段，圓柱殼處于線性前屈曲階段，其承載能力隨著位移載荷的增加呈線性增加；當最大承載能力達到0.824 2后結構發(fā)生局部失穩(wěn)；在非線性后屈曲階段（AB段），結構逐漸從局部失穩(wěn)轉化為肋間殼板環(huán)向失穩(wěn)，承載能力降低至0.4左右。

圖12 考慮實測缺陷時圓柱殼有限元模型的位移—載荷曲線Fig.12 Displacement-loading curve of cylindrical shell FE model considering measured defects

2.4 兩種重構方法的對比

使用散點法和傅里葉級數(shù)法這兩種有限元模型重構方法引入相同的實測缺陷，對圓柱殼有限元模型進行重構后進行極限承載分析，得到兩種重構方法的性能對比如圖13所示。圖中，實線表示通過兩種方法確定的不同網格尺寸下的圓柱殼屈曲載荷（歸一化無量綱結果），虛線表示兩種方法導入缺陷過程的耗時。由圖13可見，兩種方法確定的不同網格尺寸下的圓柱殼屈曲載荷幾乎沒有差別，使用散點法引入缺陷的過程耗時遠遠大于傅里葉級數(shù)法。隨著網格數(shù)的增加，引入缺陷過程的耗時呈指數(shù)增加，兩者的差異更明顯。其中，傅里葉級數(shù)法耗時最高為散點法的48%，平均耗時為散點法的32%?？梢?，采用傅里葉級數(shù)法的總體計算效率比散點法的效率高。

圖13 網格收斂曲線和兩種方法的計算耗時對比Fig.13 Mesh convergence rate curves and time cost by two methods

3 結 論

本文針對含初始幾何缺陷的圓柱殼外壓承載能力進行了數(shù)值分析，得到如下主要結論：

1）基于二次型的正交變換及傅里葉級數(shù)的初始幾何缺陷表達方法，可以解決測量坐標系到建模坐標系不統(tǒng)一的問題，形成基于實測幾何缺陷的有限元模型重構方法。

2）傅里葉級數(shù)法的計算效率明顯高于散點法，在缺陷引入過程最少可節(jié)省52%的時長，且隨著計算規(guī)模的增加，效率優(yōu)勢會更加明顯。

3）基于模態(tài)型缺陷的圓柱殼極限承載能力與基于實測缺陷的圓柱殼極限承載能力相比更加保守，可以作為極限承載能力的下限值。但是，實測缺陷為多種模態(tài)型缺陷的組合，基于實測缺陷的圓柱殼極限承載能力更接近實際值，能夠為圓柱殼結構極限承載能力的精確分析及結構優(yōu)化設計提供指導。

4）與總體模態(tài)型缺陷相比，圓柱殼極限承載力對局部模態(tài)型缺陷更加敏感，在加工建造過程中，需嚴格控制由焊接變形、加工制造等導致的肋間殼板形位偏差，保證圓柱殼的承載能力。