內壓碟形封頭塑性垮塌壓力預測方法

唐保威，徐 平

(浙江大學 應用力學研究所，杭州 310027)

0 引言

碟形封頭是壓力容器常見的部件，在石油、機械、化工、航空航天等領域有廣泛應用。因其使用環境復雜，為保障安全，碟形封頭的壁厚設計一般偏大，但這造成了材料浪費，增加了制造成本[1-3]。塑性垮塌失效是內壓碟形封頭設計需要考慮的重要失效模式，精準預測塑性垮塌壓力有助于設計碟形封頭時，在避免塑性垮塌的前提下減少用料，因此具有重要意義。

陸明萬等[4-9]介紹了基于材料硬化的真應力-真應變本構模型分析垮塌載荷的塑性分析方法，指出塑性分析方法的計算結果更為精確。ZHENG等[10-11]采用基于弧長法的非線性有限元技術，考慮實測的材料應力-應變曲線，預測一個與筒體焊接的冷沖壓標準碟形封頭塑性垮塌壓力，預測結果與試驗結果吻合好；用該方法探究了橢圓封頭各項參數對塑性垮塌壓力的影響，歸納出了橢圓封頭塑性垮塌壓力的預測公式，與試驗結果吻合良好。然而，基于有限元的數值模型要付出較大的計算工作量，過程更費時[4-7]。理論分析方面，UPDIKE等[12]將碟形封頭支撐環以上部分等效為球冠，基于von Mises各向同性硬化模型和大變形理論，采用Hill應變假設[13]推導公式計算了碟形封頭的塑性垮塌壓力，但其公式需要迭代計算，且支撐環半徑仍需采用有限元分析獲取，該方法無法獨立使用。

目前針對碟形封頭塑性垮塌壓力預測方法仍缺少深入系統的研究。本文從有限元分析、理論計算、試驗驗證三方面對該問題進行研究。首先，采用基于弧長法的非線性有限元技術，建立考慮應變硬化的碟形封頭塑性垮塌壓力預測數值模型；然后，在文獻[12]的研究基礎上，推導出可以獨立使用的碟形封頭塑性垮塌壓力計算公式，并將公式計算與有限元分析相互對比，以驗證預測方法的可靠性;最后，結合不同尺寸、材料和制造工藝的碟形封頭爆破試驗數據，對預測方法的準確性進行驗證分析。

1 有限元法

1.1 模型建立

利用ANSYS對內壓碟形封頭的塑性垮塌壓力進行有限元分析預測。工程上，碟形封頭通過法蘭與筒體連接，或直接與筒體焊接。以BLACHUT[14]試驗中的K3封頭(法蘭連接，見圖1)和MILLER等[15]試驗中的CBI1封頭(筒體焊接，見圖2)為例進行有限元分析，封頭的尺寸參數如表1所示。

圖1 K3封頭的對稱模型和網格劃分示意

圖2 CBI1封頭的對稱模型和網格劃分示意

建模時，由于碟形封頭在內壓下的幾何結構、載荷具有軸對稱特征，故建立軸對稱模型；采用Plane 182單元劃分網格；通過法蘭與筒體連接的封頭，在封頭頂部對稱軸線施加X方向的位移約束，封頭底部施加位移全約束；直接與筒體焊接的封頭，在封頭頂部對稱軸線施加X方向的位移約束，在筒體底部施加Y方向的位移約束。

表1 碟形封頭尺寸參數

材料模型采用真應力-真應變曲線(見圖3[14-15])。其中K3封頭模型采用BLACHUT[14]實測的真應力-真應變曲線。CBI1封頭則根據MILLER等[15]實測的材料彈性模量、屈服強度、抗拉強度的平均值，利用ASME BPVC Ⅷ.2—2019附錄3-D的本構關系確定其真應力-真應變曲線。

圖3 真應力-真應變曲線

1.2 計算結果

(a)壓力-應變曲線

采用基于弧長法的非線性有限元技術進行計算，弧長法可自動調整弧長增量，當曲線結構表現出非線性特征時，弧長會自動減小，迭代增多，故能很好地跟蹤壓力-應變曲線的突變點[16]。獲得K3和CBI1封頭模型在封頭頂點位置的壓力-von Mises應變曲線見圖4(a)，5(a)。

(a)壓力-應變曲線

由塑性垮塌的定義可知塑性垮塌壓力是結構能夠承載的最大壓力[4]，表現在“壓力-應變”圖像上即為曲線的最高點，故K3封頭和CBI1封頭的塑性垮塌壓力分別為50.3 MPa和2.0 MPa，兩者在塑性垮塌壓力下的等效塑性應變云圖如圖4(b)和圖5(b)所示。

2 理論方法

2.1 理論基礎

由有限元計算結果可以看出，內壓作用下碟形封頭球冠區向外變形，過渡段向內變形，因此碟形封頭在變形前后存在環向應變為零(橫向位移為零)的位置，見圖6。UPDIKE等[12]曾將環向應變為零位置及以上部分的封頭等效為具有支撐環的球冠(即將圖6中箭頭所指向的記號以上的封頭等效為圖7的模型)，并基于Hill應變假設[13]，計算內壓碟形封頭塑性垮塌壓力。

(a)K3封頭

公式推導在下列假設下進行：(1)工程上碟形封頭徑厚比大于10，認為封頭處于膜應力狀態；(2)等效球冠在變形中各處曲率始終保持一致，即等效球冠始終為球狀；(3)Hill應變假設[13](等效球冠頂點附近的微元的瞬時路徑與等效球冠輪廓截面垂直)；(4)碟形封頭所用材料大多為金屬，適用冪指數本構模型，且由于金屬材料的彈性段對塑性垮塌壓力幾乎沒影響，采用Swift本構關系；(5)材料不可壓縮。

圖7 等效球冠模型示意

根據圖7可以得到幾何關系：

h=Rφ(1-cosβ)

(1)

a=Rφsinβ

(2)

(3)

式中，h為等效球冠的高度，mm；Rφ為等效球冠的曲率半徑，mm；β為等效球冠球心和支撐環連線與對稱軸的夾角，(°)；a為等效球冠半徑(支撐環半徑)，mm。

Mises等效應力和等效應變：

σe=σθ=σφ

(4)

εe/2=εθ=εφ=-εz/2

(5)

式中,σe為Mises等效應力，MPa；σθ為環向應力，MPa；σφ為經向應力，MPa；εe為Mises等效應變；εθ為環向應變；εφ為經向應變；εz為法向應變。

t=t0exp(εZ)=t0exp(-εe)

(6)

根據Hill應變假設[13]:

(7)

式(5)(7)積分得到：

(8)

Swift本構模型[12]：

σe=C(b+εe)n

(9)

式中,C,b,n為Swift本構關系的材料參數。

薄膜理論：

(10)

2.2 本構關系

Swfit本構關系中的材料參數C,b和n可以通過下式[12]計算：

n=b+ln(1+eu)

(11)

su(1+eu)=C[b+ln(1+eu)]n=Cnn

(12)

sy(1+ey)=C[n-ln(1+eu)+ln(1+ey)]n

(13)

聯立式(11)～(13)可得到n的非線性方程：

(14)

式中,e為工程應變；ε為真應變;s為工程應力，MPa；σ為真應力，MPa；下標y為屈服點；下標u為抗拉強度點。

由式(14)首先求得n，進而得到b和C。其中，屈服強度對應的工程應變ey=sy/E。對于抗拉強度的工程應變eu，可以先確定屈強比K(K=sy/su)，采用ASME BPVC Ⅷ.2—2019附錄3-D的本構關系確定真應變εu(見式(15))，再按工程應變與真應變關系ε=ln(1+e)計算eu。

(15)

由此，獲得K3封頭所用的材料參數C,b和n分別為750，0.002 5，0.172 9；CBI1封頭所用的材料參數C,b和n分別為945，0.010 3，0.217 3，對應的Swift真應力-真應變曲線見圖8。

圖8 K3封頭和CBI1封頭的Swift曲線

2.3 塑性垮塌壓力

塑性垮塌壓力是結構能夠承載的最大壓力，表現在“壓力-應變”圖像上即為曲線的最高點，故封頭塑性垮塌條件[17]：

(16)

為求得封頭塑性垮塌壓力，避免迭代運算，需進一步獲得封頭所受壓力與等效應變的關系式。聯立式(3)(5)(7)，得等效球冠的高度：

(17)

聯立式(3)(6)(10)(17)，即可得封頭所受壓力P：

P=4t0σe[exp(-c-2.5εe)

-exp(-2c-3εe)a2]1/2

(18)

結合式(9)和式(18)，得封頭所受壓力P與等效應變εe的關系式：

P=4t0C(b+εe)n[exp(-c-2.5εe)

-exp(-2c-3εe)a2]1/2

(19)

利用式(19)可將封頭所受壓力P與等效應變εe繪制出關系曲線(見圖10)，可知曲線上最高點的壓力為塑性垮塌壓力。

(20)

對式(9)求導并與式(20)比較得：

(21)

(22)

式(22)可以直接求得碟形封頭的塑性垮塌壓力。對于一個確定的封頭，t0,C,b,n,c,a均為常量，其中a只是理論上的概念，如何獲得a是下一步計算的關鍵。

主要考慮D,Ri,r,sy,su對a的影響，記Ri/D=m,r/D=l，引入強度常數Q=100 MPa對sy和su進行無量綱處理，記sy/Q=syQ,su/Q=suQ。希望得到a=f(D,m,l,syQ,suQ)的表達式，為此利用有限元對不同條件下a的值進行仿真分析。尺寸方面考慮歐盟規范EN 13445-3UnfiredPressureVessel、ASME規范ASME BPVC.Ⅷ.2—2019Rulesforconstructionofpressurevessels,division2,alternativerules和我國GB/T 150—2011《壓力容器》中關于m,l的規定，在D=200～5 000 mm,m=0.7～1.0，l=0.05～0.35范圍內進行有限元仿真；材料的sy,su則參考GB 713—2014《鍋爐和壓力容器用鋼板》中常見牌號的sy和su，并為了研究規律適當擴大范圍，取sy=230～650 MPa，su/sy=1.1～2.6。仿真結果顯示a和D,m,l,syQ,suQ存在著極強的規律性。將所有仿真結果用SPSS(Statistical Product and Service Solutions)軟件進行非線性回歸，選擇Levenberg-Marquardt方法對回歸參數進行最小二乘估計，得到：

-0.011syQ+0.008suQ

(23)

該式可以在D=200～5000 mm，m=0.7～1.0，l=0.05～0.35，sy=230～650 MPa，su/sy=1.1～2.6的范圍內使用。將仿真得到的a與公式計算得到的a對比，如圖9所示，兩者平均誤差僅為-0.14%，標準差為1.94%。式(23)可以獨立于有限元獲得支撐環半徑a的值，使本節得到的理論公式具備了使用價值。

圖9 仿真得到的a與公式計算得到的a的比較

K3封頭模型和CBI1封頭模型等效球冠半徑a(以內徑為基準)、材料參數C,b和n見表2。圖10示出了K3封頭模型和CBI1封頭模型的壓力-等效應變曲線，從而得到K3封頭和CBI1封頭模型的塑性垮塌壓力分別為53.44和2.15 MPa。

表2 等效球冠參數

圖10 K3和CBI1封頭的塑性垮塌壓力

3 試驗驗證

為驗證碟形封頭塑性垮塌壓力預測的有限元法和公式計算法，從文獻中查閱總結16個碟形封頭的爆破試驗數據，見表3。對表中的碟形封頭分別按第1章的預測模型和第2章的公式法計算塑性垮塌壓力的有限元解和公式解，并與封頭爆破壓力試驗值進行比較。

表3 碟形封頭試驗數據匯總及誤差對比

從表3可以看出，碟形封頭塑性垮塌壓力的有限元解與公式解的相對誤差范圍為-9.6%～8.7%，平均誤差僅為2.47%，表明碟形封頭塑性垮塌壓力的有限元解與公式解的吻合程度很高，驗證了本文預測方法的可靠性。且公式法只需知道碟形封頭的D,Ri,r,t0,sy,su，便可完成計算，效率遠高于有限元，當面臨大量碟形封頭需要預測其塑性垮塌壓力時，公式法無疑會是一種便捷、高效的方法。

Blachut系列封頭(編號為K)，有限元解比試驗值高20%左右，公式解比試驗值高20%以上。該系列封頭采用鋼錠機加工制造，而鋼錠中心存在縮孔殘余，即封頭頂部存在縮孔缺陷，導致試驗封頭承載性能大大降低[14,19]。另3個Blachut系列封頭(編號為1,2,3)的縮孔缺陷相對較小，有限元解、公式解均與試驗值吻合較好。CBI系列封頭分瓣加工，封頭存在較多焊接接頭，且CBI1封頭在過渡段焊縫處破裂、CBI2封頭破口平行于球冠拼焊縫,在焊接接頭熱影響區破裂，導致分瓣加工的封頭爆破壓力試驗值低于有限元解和公式解。對于冷沖壓和冷旋壓的封頭，封頭質量較高，塑性垮塌壓力的有限元解、公式解與試驗結果吻合良好?？傮w來說，碟形封頭塑性垮塌壓力的有限元解、公式解與試驗值吻合良好，驗證了碟形封頭塑性垮塌壓力有限元計算和計算公式在理想條件下的準確性。

4 結論

(1)采用基于弧長法的非線性有限元技術，建立了考慮應變硬化的碟形封頭塑性垮塌壓力有限元預測數值模型。推導出碟形封頭塑性垮塌壓力的公式計算方法，并使之能夠獨立使用。碟形封頭塑性垮塌壓力的有限元模型和公式計算的預測結果吻合良好，表明預測方法可靠。

(2)匯總了16個不同碟形封頭爆破試驗數據。公式法、有限元法計算結果與試驗結果吻合良好，驗證了碟形封頭塑性垮塌壓力的有限元計算和公式計算在理想條件下的準確性。為建立考慮材料應變硬化的、防止碟形封頭塑性垮塌的設計方法提供技術支撐。