解析法及數值法在求解壓力容器缺陷應力強度因子中的應用 ①

□□ 張巨銀，王文達，王學平，牟龍龍，趙志琦

(1.甘肅省特種設備檢驗檢測研究院，甘肅 蘭州 730050；2.寧夏特種設備檢驗檢測院，寧夏 銀川 750000)

引言

GB/T 19624—2019《在用含缺陷壓力容器安全評定》[1]于2020年1月1日正式頒布實施。該標準以裂紋張開位移和應力強度因子為主要參量，以彈塑性雙判據法為基礎，對壓力容器缺陷進行安全評定。

應力強度因子為表征裂紋尖端應力場強度的物理量，是含缺陷壓力容器安全評定的重要內容，其求解方式主要有解析法和數值法。GB/T 19624—2019中詳細地區(qū)分了壓力容器的缺陷類型，并規(guī)定了各種缺陷的安全評估方法，同時也加入了更多缺陷評估中間參量的計算模型和方法。在平面缺陷的評估中，明確地指出了多種情況下的應力強度因子計算模型，評定人員可以根據相應的模型來便捷地計算相應的工程問題。

盡管目前標準中的計算模型可以覆蓋部分常見工程結構，但對于一些較為特殊的問題還是無從下手，而利用XFEM技術計算裂紋尖端應力強度因子則具有更好的適應性。XFEM技術可以快速準確地建立含缺陷結構的有限元計算模型。

本文以平板裂紋模型為例，分別采用GB/T 19624—2019中的解析法和XFEM數值法對相同工況下裂紋尖端的應力強度因子進行求解，并對兩種方法所得計算結果的差異性進行對比分析。

1 XFEM技術的介紹

擴展有限元法(XFEM)自Belytschko T等[2]首次提出以來，在斷裂力學領域的研究中得到了一定程度的普及與發(fā)展。ZHUANG Z等[3]基于XFEM計算出的數值結果與實驗數據進行比較，證明了XFEM計算的準確性及效率；郭歷倫等[4]對制約擴展有限元發(fā)展的兩個問題及其修正方法進行了討論；陳華等[5]基于XFEM計算了裂尖應力變化及裂紋擴展期間的能量變化，體現了擴展有限元法的準確性以及其在裂紋擴展領域得天獨厚的優(yōu)勢，為工程實際應用提供參考。

擴展有限元法的基本思路是利用不同的單元位移模型來對裂紋路徑上單元的位移不連續(xù)性進行描述。根據與裂紋擴展路徑的位置關系，有限元模型中的單元類型被劃分為三種：常規(guī)單元、裂紋貫穿單元以及裂紋處單元，如圖1所示。

圖1 數值模型單元示意圖

三種單元的位移模式見式(1)～(3)。

ui(x)=∑Niui

(1)

uj(x)=∑(Niui+NjH(x)uj)

(2)

uk(x)=∑(Niui+NkφSuk)

(3)

式中，i、j、k分別對應常規(guī)單元、裂紋貫穿單元以及裂紋處單元，Ni、Nj和Nk分別為三種單元的位移形函數，ui、uj以及uk為相應單元的節(jié)點位移。H(x)與φS為Heaviside函數與裂尖增強函數，分別用于描述裂紋面與裂紋尖端的不連續(xù)性，具體見式(4)和(5)。

式中，r與θ為以裂紋尖端的原點的極坐標系坐標值。

通常，XFEM中裂紋尖端的應力強度因子會采用相互作用積分法來計算，其求解的相應過程在文獻[6-7]中有較為詳細的說明，這里就不再贅述。本文的數值計算過程在ALOF軟件中進行，它是基于XFEM技術開發(fā)的國產有限元軟件，可以快速計算復雜三維模型上裂紋的應力強度因子，在一些工程案例中得到了良好的應用。

2 解析解計算模型

GB/T19624—2019《在用含缺陷壓力容器安全評定》提供了常見含缺陷結構缺陷尖端部位的應力強度因子K的計算公式，在此主要介紹標準中平板(板寬2W，板長2L，板厚B)表面裂紋和埋藏裂紋的應力強度因子求解的計算模型。

GB/T19624規(guī)定在計算應力強度因子時，所取用的應力是缺陷部位的主應力，計算該主應力時采用線彈性計算方法，并假設結構中不存在缺陷。沿板厚方向利用應力線性化規(guī)則計算薄膜應力σm與彎曲應力σB，計算方法見式(6)和式(7)。

式中，σ1、σ2可分別取模型厚度方向兩個表面的第一主應力。

裂紋尖端的應力強度因子表達式見式(8)。

式中，a為裂紋半長，fm、fb分別為裂紋構形因子。

對于半橢圓表面裂紋，其構形因子計算見式(9)～(12)。

(9)

式中，c為裂紋深度，B為平板厚度，上標A與B分別為裂紋深度處K的裂紋構形因子與裂紋長度方向兩端處K的裂紋構形因子。

對于橢圓埋藏裂紋，其構形因子計算見式(13)～(16)。

(14)

(15)

(16)

式中，c為裂紋高度的一半，B為平板厚度，上標A與B分別為裂紋深度處K的裂紋構形因子與裂紋長度方向兩端處K的裂紋構形因子。

3 解析解模型與結果

依據GB/T19624的計算模型來計算平板上表面裂紋的應力強度因子。在通用有限元軟件中建立如圖2所示的正方體平板模型，其中，L=160mm，B=30mm，左端為固定邊界，右端施加140MPa的拉力。材料彈性模量取210 000MPa，泊松比取0.3。

圖2 數值模型

在計算完成后，沿板厚方向取上下兩個面的第一主應力σ1、σ2，其值均為143MPa，進一步根據應力線性化規(guī)則和公式(6)、(7)計算得到σm與σB，分別為143MPa與0。

平板上的裂紋的尺寸見表1。

表1 裂紋尺寸 mm

將以上數據代入標準第二部分的計算公式中即可計算得到裂紋的K1值，計算結果見表2。

表2 解析解計算結果 N/mm3/2

4 數值解模型與結果

采用XFEM技術實現平板裂紋的應力強度因子數值法求解，主要步驟如下：

(1)建立正方形平板數值模型，其尺寸和受力狀況均與上文解析法求解部分所建立的模型保持一致。

(2)建立一個半橢圓平面型裂紋，裂紋的短半軸長度為5mm，長半軸長度為8mm。采用三角形網格對其進行離散。

(3)建立一個橢圓埋藏型裂紋，裂紋的短半軸長度為2.5mm，長半軸長度為8mm。采用三角形網格對其進行離散。

(4)分別將平板模型和兩種裂紋模型進行耦合，形成含缺陷模型。圖3為耦合得到的含半橢圓裂紋缺陷平板的數值模型；圖4為表面裂紋和埋藏裂紋分別在平板模型壁厚方向所處狀態(tài)的局部視圖。

圖3 含半橢圓裂紋缺陷的平板模型

圖4 裂紋在平板壁厚方向所處位置的示意圖

(5)完成計算，直接輸出應力強度因子值。表面裂紋最大應力強度因子計算結果為：KI=451.16N/mm3/2。埋藏裂紋最大應力強度因子計算結果為：KI=326.54N/mm3/2。

表3為解析計算結果與XFEM數值計算結果的對比，可以看出兩種方法計算結果的最大相對誤差在10%左右，數值法結果普遍小于解析法結果。同時，數值法結果更加求準，解析法結果更加保守，更有利于保證含缺陷承壓設備的安全使用。當然XFEM數值法在計算含缺陷壓力容器的應力強度因子中也具有較好的適用性。

表3 規(guī)范計算結果與數值解最大值對比

5 結論

5.1GB/T19624—2019中解析解計算流程明確，算法直接，但其所用公式復雜，并且需要借助有限元軟件對模型前期的受力狀態(tài)進行求解。

5.2XFEM技術對求解應力強度因子有良好適應性，其建模簡單，計算結果輸出方便，需要使用者掌握一定的有限元基礎知識。

5.3 兩種方法計算平板表面裂紋模型應力強度因子結果的相對誤差在10%左右，可以達到工程計算精度要求，應視具體情況選用合適的方法進行計算。