基于實(shí)測(cè)幾何缺陷的圓柱殼外壓承載能力分析

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

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

2武漢第二船舶設(shè)計(jì)研究所，湖北武漢430205

0 引 言

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

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

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

1 幾何缺陷的測(cè)量與表達(dá)

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

圖1 三維形貌測(cè)量系統(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

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

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

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

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

其矩陣形式寫為

其中，

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

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

式中，Λ為對(duì)角矩陣。

進(jìn)一步，式（2）可以變換為標(biāo)準(zhǔn)坐標(biāo)系下的式（4）

其中，

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

通過(guò)正交變換［10］，將二次曲面的一般方程轉(zhuǎn)換為標(biāo)準(zhǔn)方程。該變換消除了坐標(biāo)分量間的交叉項(xiàng)，通過(guò)校正附加的特征點(diǎn)，將測(cè)量坐標(biāo)系與標(biāo)準(zhǔn)坐標(biāo)系的方向?qū)R，同時(shí)保持曲面形狀不變。

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

圖4 圓柱殼坐標(biāo)系下初始幾何缺陷的散點(diǎn)圖Fig.4 Scattered points diagram of initial geometric defects in cylinder coordinate system

圓柱殼結(jié)構(gòu)的真實(shí)缺陷對(duì)結(jié)構(gòu)承載能力有著極大的影響，是多年來(lái)研究的熱點(diǎn)和難點(diǎn)問(wèn)題。國(guó)內(nèi)外學(xué)者通常將圓柱殼的真實(shí)缺陷（離面位移場(chǎng)w(x，y)）展開(kāi)為半波余弦或半波正弦的雙重傅里葉級(jí)數(shù)。

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

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

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

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

首先，本文將如圖4所示的圓柱殼坐標(biāo)系下的點(diǎn)云信息通過(guò)逆距離加權(quán)插值法進(jìn)行插值，獲得了如圖5所示的格柵化離面位移圖，從而將離面位移場(chǎng)表示為二維矩陣的形式。

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

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

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

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

圖7 含實(shí)測(cè)缺陷的有限元圓柱殼重構(gòu)模型Fig.7 FE model of reconfigured cylindrical shell with measured defects

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

本文所示完美模型是指不含實(shí)測(cè)缺陷且經(jīng)過(guò)完善的圓柱殼模型。該完美圓柱殼模型的名義半徑R是通過(guò)對(duì)實(shí)際結(jié)構(gòu)的實(shí)測(cè)形貌經(jīng)過(guò)特征提取處理獲得的擬合半徑。圖8所示為模型的邊界條件及載荷。

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

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

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

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

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

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

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

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

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

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

3）施加軸外壓載荷，并進(jìn)行計(jì)算分析。

分別引入如圖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)型缺陷時(shí)圓柱殼有限元模型的位移—載荷曲線Fig.11 Displacement-loading curves of cylindrical shell FE model considering mode defects

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

2.3 實(shí)測(cè)缺陷近似分析

由于實(shí)測(cè)缺陷真實(shí)地反映了圓柱殼結(jié)構(gòu)本身的初始幾何缺陷，所以基于結(jié)構(gòu)的實(shí)測(cè)缺陷開(kāi)展極限承載能力分析能夠非常精確地預(yù)測(cè)其極限承載能力。本文使用的方法為實(shí)測(cè)缺陷近似方法（Measured Imperfection Approach，MIA）。

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

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

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

2.4 兩種重構(gòu)方法的對(duì)比

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

圖13 網(wǎng)格收斂曲線和兩種方法的計(jì)算耗時(shí)對(duì)比Fig.13 Mesh convergence rate curves and time cost by two methods

3 結(jié) 論

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

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

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

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

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