董乃逸，賀小華

(南京工業(yè)大學 機械與動力工程學院，南京 211816)

0 引言

球殼封頭結(jié)構(gòu)簡單，受力均勻，廣泛應用于高壓容器、深海耐壓潛水器等承壓結(jié)構(gòu)[1]。因工藝需要，耐壓球殼需開設各類接管結(jié)構(gòu)，開孔接管破壞了球殼結(jié)構(gòu)的整體連續(xù)性，在外壓作用下，導致球殼承壓能力下降，易產(chǎn)生屈曲失穩(wěn)破壞。對球形封頭開孔接管結(jié)構(gòu)屈曲特性及臨界失穩(wěn)壓力進行研究意義重大。

對于外壓球殼的屈曲行為研究由來已久。早期針對理想球殼，ZOELLY[2]基于小變形理論，推導了外壓完整球殼的臨界失穩(wěn)壓力公式。由于Zoelly公式計算得到的結(jié)果與試驗值有較大差距，KARMAN[3]基于外壓球殼的軸對稱失穩(wěn)問題，提出了大變形理論公式以計算外壓球殼的臨界失穩(wěn)壓力，結(jié)果與試驗值相對較為接近。20世紀60年代，KRENZKE等[4]對200多個具有不同缺陷的耐壓球殼進行試驗，提出了臨界失穩(wěn)壓力的計算公式。目前，國內(nèi)外標準規(guī)范：GB/T 150—2011《壓力容器》、ASME Ⅷ-1和ASME Ⅷ-2，以及EN 13445-3:2014分別基于上述小變形理論及修正公式，提出了臨界失穩(wěn)壓力的計算公式。文獻[5]分析描述了球殼在均勻外壓下的前屈曲狀態(tài)，使用有限元計算球殼的臨界失穩(wěn)壓力和屈曲模式，顯示了球殼的后屈曲行為；文獻[6]研究了外壓半球封頭的彈塑性屈曲，提出了預測半球封頭臨界失穩(wěn)壓力的解析公式，并與數(shù)值結(jié)果進行對比。

由于球形封頭開孔接管結(jié)構(gòu)在其開孔處的應力狀態(tài)較為復雜，難以通過解析方法求解球形封頭開孔接管結(jié)構(gòu)臨界失穩(wěn)壓力。文獻[7]基于外壓球殼開孔區(qū)域的軸對稱性和等面積補強法，對大深度潛水器載人球殼開孔加強形式進行幾何變換，完善了球殼開孔加強理論計算方法，并對球殼進行外壓試驗，驗證了計算方法的有效性和適用性。文獻[8]運用有限元軟件ANSYS，分析了側(cè)壁加強、墊板加強、肘板加強這三種耐壓球殼開孔加強方式的極限強度并進行對比，結(jié)果表明側(cè)壁加強的極限強度最高；同時考慮到結(jié)構(gòu)性能和工藝性方面，提出了梯形加強這一新型的開孔加強形式。文獻[9]運用有限元軟件ANSYS對潛水器開孔耐壓殼半球封頭承載能力進行非線性分析，討論了開孔個數(shù)、孔徑大小、圍壁厚度3個參數(shù)對耐壓殼半球封頭承載力的影響，結(jié)果表明，開孔越多、孔徑越大、圍壁厚度越小都會使球形封頭承載能力下降。

本文系統(tǒng)討論開孔率di/Di、接管與封頭厚度比δet/δe及徑厚比Di/δe這3個無量綱結(jié)構(gòu)參數(shù)對球形封頭開孔接管結(jié)構(gòu)臨界失穩(wěn)壓力Pcr及屈曲特性的影響規(guī)律，為外壓球形封頭開孔接管結(jié)構(gòu)穩(wěn)定性研究及工程設計提供參考。

1 有限元模擬計算

1.1 材料性能參數(shù)

球形封頭材料采用Q345R，接管材料采用16Mn，常溫下的材料性能參數(shù)如表1所示。

表1 常溫下材料性能參數(shù)

1.2 有限元計算方案

有限元分析結(jié)構(gòu)為球形封頭徑向開孔接管結(jié)構(gòu)，如圖1所示。球形封頭厚度取定值δe=10 mm，取筒體直邊段長度L1=25 mm，接管外伸長度L2=200 mm，通過改變開孔率di/Di、接管與封頭厚度比δet/δe和封頭徑厚比Di/δe的值進行有限元模擬計算。di/Di分別取值0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8；δet/δe分別取值0.2,0.4,0.6,0.8,1.0,1.25,1.5,1.75,2.0；Di/δe分別取值50,100,150,200。有限元分析模型共計288組。

1.3 材料本構(gòu)關(guān)系

文獻[10]指出本構(gòu)關(guān)系對殼體結(jié)構(gòu)的臨界失穩(wěn)壓力影響甚少，因此,本文選取雙線性材料本構(gòu)模型進行有限元分析，塑性段斜率取E/100，其本構(gòu)關(guān)系如圖2所示。

圖1 球形封頭徑向開孔接管結(jié)構(gòu)

圖2 雙線性材料本構(gòu)模型

1.4 單元類型及邊界條件

分析結(jié)構(gòu)選擇Shell 181殼單元，該單元計算精度已在文獻[10-11]中證實。在滿足網(wǎng)格無關(guān)性要求下進行網(wǎng)格劃分，分析結(jié)構(gòu)的邊界條件為封頭直邊段施加固支約束，開孔接管端面施加環(huán)向約束以保證接管端面的圓度，封頭與接管均受外壓作用。球形封頭開孔接管有限元模型見圖3。

圖3 封頭開孔接管有限元模型

2 有限元模型驗證

2.1 完整球殼臨界失穩(wěn)壓力值比較

對球殼屈曲行為進行的研究歷程是從最初的不考慮缺陷，到后來發(fā)現(xiàn)初始缺陷對其屈曲耐壓特性有巨大的影響。實際工程中半球封頭是沖壓成型的，由于模具和工藝不可避免地會造成一定的缺陷，從而使半球封頭的實際形狀和完美的半球封頭存在一定的差距。本文選取文獻[12]中的4個試驗模型，模型參數(shù)分別為球殼直徑D=117.54,117.98,117.54,117.5 mm，計算厚度分別為0.422,0.423,0.406,0.415 mm，彈性模量E=1.93×105MPa，泊松比μ=0.28，屈服強度σs=205 MPa。分別從理論公式、試驗結(jié)果和有限元數(shù)值模擬三方面，對外壓球殼的臨界失穩(wěn)壓力進行計算對比。

小變形理論公式是基于完美幾何形狀且無缺陷的球殼，其臨界失穩(wěn)壓力的計算公式[2]：

Pcr=1.21E(δe/R)2

(1)

式中，E為材料彈性模量,MPa；δe為球殼計算厚度,mm；R為球殼半徑,mm。

大變形理論公式基于結(jié)構(gòu)的非線性因素影響，其臨界失穩(wěn)壓力的計算公式[3]：

Pcr=0.36516E(δe/R)2

(2)

本文有限元數(shù)值模擬時施加的缺陷方法有兩種，一種是采用“一致缺陷模態(tài)法”[13]設置初始幾何缺陷(一致缺陷模態(tài)法是基于最低階屈曲模態(tài)，從而模擬結(jié)構(gòu)的初始缺陷分布),通過ANSYS中的UPGEOM命令里的FACTOR值對模型統(tǒng)一施加10%的幾何缺陷，然后進行非線性求解；另一種是施加局部厚度缺陷。GB/T 150.4—2011《壓力容器 第4部分：制造、檢驗和驗收》中6.4.2規(guī)定球形封頭的最大形狀偏差外凸不超過1.25%Di，內(nèi)凹不超過0.625%Di。文獻[14-15]表明內(nèi)凹0.625%Di下臨界失穩(wěn)壓力低于外凸1.25%Di下臨界失穩(wěn)壓力，因此本文模型采用形狀偏差內(nèi)凹0.625%Di的球形封頭模型，并在封頭頂部最大壁厚減薄位置施加單個缺陷厚度減薄15%進行非線性分析。最終求解提取出位移最大點的載荷-位移圖，并根據(jù)零曲率準則得到臨界失穩(wěn)壓力值Pcr。

將文獻[12]試驗結(jié)果、公式(1)解、公式(2)解及本文有限元特征值、有限元非線性解1(一致缺陷模態(tài)法)以及有限元非線性解2(施加局部厚度缺陷)列于表2，并將其數(shù)據(jù)作成圖4，以便于直觀比較。

表2 球形封頭臨界失穩(wěn)壓力值Pcr比較

圖4 球形封頭臨界失穩(wěn)壓力值Pcr比較

由表2和圖4可以看出，有限元特征值與小變形公式(1)解結(jié)果較為接近，由于這兩種方法均未考慮幾何非線性以及缺陷的存在，導致其結(jié)果遠大于其他解；大變形公式(2)解略大于試驗值，計算結(jié)果偏于保守；有限元非線性解1及有限元非線性解2兩者均考慮了結(jié)構(gòu)的幾何大變形和初始幾何缺陷，因此其值均與試驗結(jié)果較為接近。考慮到一致缺陷模態(tài)法結(jié)果的有效性及操作的方便性，后續(xù)有限元模型均采用此法設置初始幾何缺陷。

2.2 開孔球殼臨界失穩(wěn)壓力結(jié)果比較

選用文獻[16]中按1∶10比例自制的不銹鋼縮比試驗模型，模型結(jié)構(gòu)如圖5(a)所示。該模型參數(shù)為：球殼外半徑R0=100 mm，球殼厚度t=6.3 mm，接管內(nèi)徑d1=50 mm，接管外徑d=77 mm，接管厚度t1=13.5 mm，內(nèi)伸高度h1=32.2 mm，外伸高度h2=32.2 mm，彈性模量E=1.95×105MPa，泊松比μ=0.3，屈服強度σs=300 MPa。試驗測得的縮比模型臨界失穩(wěn)壓力為20.5 MPa，采用一致缺陷模態(tài)法，施加10%的初始缺陷進行有限元非線性求解，有限元模型如圖5(b)所示，得到臨界失穩(wěn)壓力為19.223 MPa，與試驗值比較，誤差為6.64%。兩者誤差基本在工程允許范圍內(nèi)，由此也說明本文非線性屈曲分析方法的有效性。

(a)開孔球殼圍壁結(jié)構(gòu)示意

(b)開孔球殼圍壁有限元模型

3 球形封頭開孔接管結(jié)構(gòu)有限元非線性屈曲分析

3.1 臨界失穩(wěn)壓力削弱系數(shù)K

為了分析討論開孔接管結(jié)構(gòu)對球形封頭臨界失穩(wěn)壓力的影響，通過引入削弱系數(shù)K來反映開孔接管的作用，K值為球形封頭開孔接管結(jié)構(gòu)的臨界失穩(wěn)壓力Pcr與球形封頭未開孔結(jié)構(gòu)的臨界失穩(wěn)壓力Pcrs(有限元非線性模擬值)之比，即K=Pcr/Pcrs。

3.1.1 接管與封頭厚度比δet/δe與K值關(guān)系

圖6示出了不同徑厚比Di/δe及不同開孔率di/Di下，臨界失穩(wěn)壓力削弱系數(shù)K與δet/δe的分布圖?？梢钥闯?，不同徑厚比及不同開孔率下，K值與δet/δe變化規(guī)律基本相同；同一δet/δe下，K值隨di/Di的增大而減??；當δet/δe較小時，K值明顯偏低，隨著δet/δe增大，K值顯著增大，當δet/δe≥1.0后，K值隨δet/δe增大變化趨于平緩。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

3.1.2 開孔率di/Di與K值關(guān)系

圖7示出了不同徑厚比Di/δe及不同接管與封頭厚度比δet/δe下，臨界失穩(wěn)壓力削弱系數(shù)K與開孔率di/Di的分布圖。可以看出，Di/δe與δet/δe一定時，K值隨著di/Di增大而減??；當δet/δe較小時，K值隨di/Di增大顯著降低；而當δet/δe較大時，K值隨di/Di增大、降低幅度較小；當δet/δe≥1.0時，同一di/Di下，δet/δe對臨界失穩(wěn)壓力削弱系數(shù)K的影響較小。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

3.2 計算模型失穩(wěn)模態(tài)分析

由圖6可以看出，臨界失穩(wěn)壓力削弱系數(shù)K與接管與封頭厚度比δet/δe變化有個轉(zhuǎn)折點，當δet/δe≥1.0時，K值隨δet/δe增大變化不顯著。同樣由圖7可以看出，當δet/δe≥1.0時，同一開孔率下，不同計算模型臨界失穩(wěn)壓力削弱系數(shù)K值較為接近。為進一步討論球形封頭開孔接管結(jié)構(gòu)屈曲特性，圖8示出徑厚比Di/δe=150時部分計算模型的失穩(wěn)模態(tài)圖。

圖8 球形封頭中心開孔接管結(jié)構(gòu)失穩(wěn)模態(tài)圖(Di/δe=150)

由圖8(a)可以看出，較小開孔率(di/Di=0.4)，當接管較薄(δet/δe=0.2)時，計算模型失穩(wěn)位置位于接管；當δet/δe增大到0.6時，由圖8(b)可知，計算模型失穩(wěn)位置由接管轉(zhuǎn)移到封頭-接管連接處；當δet/δe增至2.0時，由圖8(c)可知，計算模型失穩(wěn)位置位于封頭處。

由圖8(d)～(f)可看出，當接管較薄(δet/δe=0.2)時，計算模型失穩(wěn)位置位于接管，隨著接管壁厚增大，計算模型失穩(wěn)位置逐漸由接管轉(zhuǎn)移到封頭。由此，解釋了圖6中臨界失穩(wěn)壓力削弱系數(shù)K與δet/δe分布圖中轉(zhuǎn)折點的由來。當接管與封頭厚度比δet/δe低于圖6中轉(zhuǎn)折點時，接管較薄，承壓能力較弱，計算模型失穩(wěn)位置位于接管，此時分析結(jié)構(gòu)臨界失穩(wěn)壓力較小；當δet/δe增大到一定值時，接管承壓能力超過封頭，計算模型失穩(wěn)位置轉(zhuǎn)到封頭，此時分析結(jié)構(gòu)臨界失穩(wěn)壓力較大，進一步由圖6可以看出，隨著開孔率di/Di增大，失穩(wěn)位置由接管轉(zhuǎn)到封頭所需接管與封頭厚度比越大。本文計算參數(shù)范圍內(nèi)，當δet/δe≥1.0時，計算模型失穩(wěn)位置均位于封頭上。

3.3 外壓球形封頭開孔接管結(jié)構(gòu)臨界接管與封頭厚度比(δet/δe)c

由圖6可以看出，臨界失穩(wěn)壓力削弱系數(shù)K與δet/δe變化存在轉(zhuǎn)折點。當δet/δe低于轉(zhuǎn)折點對應值時，計算模型失穩(wěn)位置位于接管；而當δet/δe高于轉(zhuǎn)折點對應值時，計算模型失穩(wěn)位置位于封頭。本文將失穩(wěn)位置由接管轉(zhuǎn)向封頭的最小接管與封頭厚度比稱為臨界接管與封頭厚度比 (δet/δe)c。表3列出本文計算參數(shù)范圍下的臨界接管與封頭厚度比(δet/δe)c。

表3 臨界接管與封頭厚度比(δet/δe)c

由表3可以看出，臨界接管與封頭厚度比(δet/δe)c是3組定值(0.6,0.8,1.0)，其與徑厚比及開孔率之間的關(guān)聯(lián)式如下：

(3)

表3中列出的臨界接管與封頭厚度比(δet/δe)c基本對應于圖6中曲線的轉(zhuǎn)折點。由圖6可以看出，當δet/δe<(δet/δe)c時，臨界失穩(wěn)壓力削弱系數(shù)K明顯較小，結(jié)構(gòu)承載能力較弱；當δet/δe≥(δet/δe)c時，K值相對較大，且K值隨δet/δe變化不顯著。實際工程設計中，建議參考表3中給出的臨界接管與封頭厚度比，取δet/δe=(δet/δe)c，可以在滿足結(jié)構(gòu)輕量化設計基礎上，獲得較優(yōu)的結(jié)構(gòu)承載能力。

3.4 臨界失穩(wěn)壓力削弱系數(shù)K值與結(jié)構(gòu)參數(shù)的經(jīng)驗關(guān)系式

基于臨界失穩(wěn)壓力削弱系數(shù)K隨Di/δe，di/Di，δet/δe的變化曲線，進一步擬合出K值與3個結(jié)構(gòu)參數(shù)的經(jīng)驗關(guān)系式。由于所研究參數(shù)δet/δe變化范圍(0.2～2.0)較大，且K值隨δet/δe的變化可分為兩個階段。為了保證擬合精度，采用分段擬合方法。當δet/δe位于[0.2,0.6]區(qū)間內(nèi)，采用以下關(guān)系式進行擬合：

(4)

式中，a1～a7為待定系數(shù)。

基于有限元模擬所得數(shù)據(jù)點，最終擬合結(jié)果為：

(5)

當δet/δe位于(0.6,2.0]區(qū)間內(nèi)，基于有限元模擬所得數(shù)據(jù)點，最終擬合結(jié)果為：

(6)

為了驗證所得關(guān)系式的可靠性，進一步給出當開孔率di/Di=0.2,0.5,0.8時，利用以上經(jīng)驗關(guān)系式所得擬合K值與實際有限元K值數(shù)據(jù)點對比圖，如圖9所示。圖9中補充計算了δet/δe=0.7時相應有限元模型的臨界壓力削弱系數(shù)K，并與第二區(qū)間段公式擬合值進行比較。由圖9可以看出，不同區(qū)間內(nèi)擬合所得K值與實際有限元K值兩者較為接近，且擬合曲線能夠反映出K值隨不同結(jié)構(gòu)參數(shù)的變化趨勢。

(a)Di/δe=50

(b)Di/δe=100

(c)Di/δe=150

(d)Di/δe=200

計算表明,除了δet/δe=0.2時個別點誤差較大，其他點公式擬合值和有限元模擬值誤差較小，所有分析計算點平均誤差小于10%。

4 結(jié)論

分析討論3個結(jié)構(gòu)參數(shù)：徑厚比Di/δe、開孔率di/Di、接管與封頭厚度比δet/δe對外壓球形封頭開孔接管結(jié)構(gòu)臨界失穩(wěn)壓力Pcr的影響。研究結(jié)果如下。

(1)球形封頭臨界失穩(wěn)壓力削弱系數(shù)K值隨著開孔率di/Di增大而減小。

(2)接管與封頭厚度比較小時，臨界失穩(wěn)壓力削弱系數(shù)K值越小，計算模型失穩(wěn)位置位于接管，當接管與封頭厚度比增至一定值時，K值顯著增大，計算模型失穩(wěn)位置位于封頭。

(3)開孔率越大，失穩(wěn)位置由接管轉(zhuǎn)到封頭所需接管與封頭厚度比越大。本文計算參數(shù)范圍內(nèi)，當δet/δe≥1.0時，計算模型失穩(wěn)位置均位于封頭上。

(4)提出了臨界接管與封頭厚度比(δet/δe)c，實際工程設計中，當接管與封頭厚度比取臨界接管與封頭厚度比時，可以在滿足結(jié)構(gòu)輕量化設計基礎上，獲得較優(yōu)的結(jié)構(gòu)承載能力。

(5)建立了臨界失穩(wěn)壓力削弱系數(shù)K值與不同結(jié)構(gòu)參數(shù)的經(jīng)驗關(guān)系式，對比發(fā)現(xiàn)所得關(guān)系式預測精度較高。