含坑點腐蝕的殼體有限元方法

徐 強，萬正權

（中國船舶科學研究中心，江蘇 無錫 214082）

1 引 言

殼體在工程實際中得到大量應用，例如，航海和海洋工程中的艦船、潛艇和油氣管道等都廣泛地采用了各種形式的殼體結構。老齡殼體結構多數(shù)存在嚴重的腐蝕問題，坑點腐蝕是常見且危害極大的一種，坑點腐蝕屬于局部腐蝕（又稱非均勻腐蝕），是由于金屬表面物理和化學的不均勻性所引起的。嚴重的坑點腐蝕可能使金屬殼板穿透，引起流體泄漏、火災、爆炸等事故，因此，有必要開展坑點腐蝕對殼體結構強度的影響研究。近年來，國內(nèi)外學者對工程結構坑點腐蝕問題的重視程度在不斷上升，并開展了系列相關問題的研究。

TSCF[1]根據(jù)坑點腐蝕鋼板的彎曲試驗與理論推導，提出一個包括彎曲剛度、質(zhì)量、邊界條件、板尺寸規(guī)格等參數(shù)的腐蝕鋼板等效厚度計算公式。

Daidora[2]等人提出用腐蝕坑數(shù)據(jù)的平均值和最大值或者用腐蝕坑的個數(shù)和坑的最大深度來估計結構腐蝕后的剩余厚度，然后通過概率方法提出一種用于評估腐蝕厚度對板局部屈服和板屈曲影響的方法。

Flaks[3]提出一種數(shù)學方法，用于評估坑點腐蝕對受拉鋁板的極限強度的影響。通過對自然腐蝕鋁板的拉伸試驗，得到了一個考慮拉伸強度、屈服強度和硬度損失的系數(shù)。

Paik[4-5]等通過試驗和有限元模擬分析，研究了腐蝕和疲勞等因素對老齡船舶的極限強度的影響。提出以“極限強度破壞因子”來衡量坑點腐蝕引起的剖面面積減少對板構件強度的影響；并認為在評估結構強度時，坑點腐蝕所造成的截面積損失是比坑點腐蝕深度更值得關注的指標。

Paik[6]研究了受坑點腐蝕影響的鋼板構件在面內(nèi)剪切載荷作用下的極限強度特性。針對不同的坑點腐蝕強度和不同的板構件幾何特性，對受面內(nèi)剪切載荷的板構件進行了一系列的非線性有限元分析。通過對有限元計算結果進行回歸分析，得到了閉合形式的板構件的極限強度設計公式。

Dunbar[7]討論了坑點腐蝕對船舶結構中的典型板構件和加筋板構件的影響。用有限元方法對板的初始屈曲、極限破壞和后極限響應進行了研究分析。有限元模型考慮了幾何缺陷和殘余應力，并把有限元模擬結果與理論計算結果和有效的實驗測量進行對比。該文提出的受損船舶結構的完整性理論有助于發(fā)展更有效的船舶維修策略。

Nakai、Yamamoto[8]等開展了坑點腐蝕對散貨船貨艙肋骨局部強度的影響研究。對老齡散貨船和油船的腐蝕情況進行了調(diào)查統(tǒng)計，并在此基礎上分別對真實坑點腐蝕樣本和人工模擬坑點腐蝕樣本進行了拉伸、壓縮試驗，考察了不同的腐蝕分布情況和平均腐蝕厚度對結構強度的影響，并進行了有限元模擬分析。

Nakai[9]還研究了局部坑點腐蝕對工型材腹板斷裂的影響，研究發(fā)現(xiàn)，腐蝕的不同分布情況會顯著地影響腹板的斷裂行為。通過進一步的有限元分析發(fā)現(xiàn)，當坑點腐蝕腹板和均勻腐蝕腹板的平均厚度損失相同時，前者的極限強度要小于或等于后者的極限強度。

王燕舞、黃小平、崔維成[10]開展了船舶結構鋼海洋環(huán)境點蝕模型之最大點蝕深度時變模型研究。針對我國船舶結構常用炭鋼、低合金鋼的實海腐蝕試驗觀測數(shù)據(jù)建立了新型Weibull函數(shù)形式的點蝕最大深度模型。

王燕舞、吳曉源等[11]開展了船舶結構鋼海洋環(huán)境點蝕模型之實船蝕坑形態(tài)與徑深比時變模型研究?；邳c蝕生長原理與實船蝕坑形態(tài)檢測數(shù)據(jù)，對船舶受點蝕構件的兩類主要蝕坑的形態(tài)、形成機制及其隨船齡演化情況進行了描述。依據(jù)Yamamoto[8]實測散貨船貨艙區(qū)肋骨蝕坑數(shù)據(jù)并與點蝕深度模型[10]相結合，建立了相應的蝕坑徑深比時變模型。

目前，國內(nèi)外學者在對坑點腐蝕幾何形態(tài)的具體處理上還存在較大差異，如：將蝕坑處理成圓柱形、球冠形、圓錐形、矩形槽等[6-9，11-12]，這種差異主要是由腐蝕環(huán)境和金屬材料的不同所引起的。為了便于展開本文的研究工作，并考慮到蝕坑幾何形態(tài)的簡潔性，本文將坑點腐蝕取成圓柱形。

本文首先建立了坑點腐蝕的力學模型，將單邊坑點腐蝕殼板沿厚度方向分為腐蝕層和完好層，求出了腐蝕層的等效材料常數(shù)，開展了坑點腐蝕的應力集中分析。然后采用有限元理論、以超參數(shù)殼元為基礎推導了坑點腐蝕殼體單元，主要導出了單元的剛度矩陣和等效結點載荷向量的有限元表達格式，并在求解單元應力時考慮了坑點腐蝕應力集中的影響。最后開展了基于坑點腐蝕殼體單元的算例分析，并與ANSYS實體單元的計算結果進行了對比。

2 坑點腐蝕的力學模型

2.1 坑點腐蝕殼板的幾何模型

坑點腐蝕殼板的幾何模型如圖1所示，本文考慮的腐蝕形態(tài)為單邊腐蝕，殼板沿厚度方向分成兩層：腐蝕層（B）和完好層（C）。在B中，坑點腐蝕呈圓柱狀，沿殼板厚度方向穿透腐蝕層。腐蝕層（B）的厚度（即腐蝕坑的深度）取實測蝕坑深度的平均值。完好層（C）和腐蝕層（B）的彈性模量和泊松比分別為E，ν和Ec，νc，E，ν為材料常數(shù)，而后者為力學等效的彈性模量和泊松比。

腐蝕層的等效材料常數(shù)Ec，νc定義如下：

式中，P為拉伸載荷，S為腐蝕層沿拉伸方向的橫截面積，L、W分別為腐蝕層的長度和寬度，ΔL、ΔW分別為腐蝕層沿拉伸方向和橫向的變形。根據(jù)（1）式進一步得到無量綱化的等效材料常數(shù)：

本文通過有限元數(shù)值分析求kE、kν，腐蝕層的有限元模型如圖2所示，計算模型分為蝕坑有序分布和蝕坑隨機分布兩組。

經(jīng)過有限元數(shù)值計算并對結果進行回歸分析，得到腐蝕層的等效材料常數(shù)回歸公式如下：

式中DOP為腐蝕強度，它等于坑點腐蝕的面積與殼板面積的比值，即

其中，Api是第 i個蝕坑的面積，且n為蝕坑的個數(shù)，l、w分別為殼板的長度和寬度。

2.2 坑點腐蝕的應力集中分析

坑點腐蝕殼板的應力集中可分為薄膜應力集中和彎曲應力集中，如圖4。應力集中系數(shù)定義為腐蝕部位的應力峰值與基準應力的比值，基準應力取假設不存在腐蝕時腐蝕部位對應的應力值。

影響坑點腐蝕應力集中系數(shù)的因素包括蝕坑深度與殼板厚度之比（h/T h），雙向壓應力的比值（σ1/σ2），雙向彎矩的比值 （M1/M2）和坑點腐蝕強度（DO P），以及蝕坑的半徑/深度比（R/H）。

按照 0≤h/Th≤0.3，0.5≤σ2/σ1≤0.65，0≤DOP≤0.3，R/H=2 進行一系列有限元數(shù)值計算并對計算結果進行回歸分析，得到薄膜應力集中系數(shù)回歸公式如下：

按照 0≤h/Th≤0.3，0.24≤M1/M2≤0.3，0≤DOP≤0.3，R/H=2進行一系列有限元數(shù)值計算并對計算結果進行回歸分析，得到彎曲應力集中系數(shù)回歸公式如下：

坑點腐蝕應力集中系數(shù)回歸曲線如圖5所示。

3 坑點腐蝕殼體單元推導

3.1 單元剛度矩陣

坑點腐蝕殼體單元的剛度矩陣以超參數(shù)殼元的剛度矩陣為基礎，通過沿單元厚度方向腐蝕層和完好層的分段積分求得。若坐標系取自然坐標ε、η、ζ（如圖6所示），則坑點腐蝕殼體單元的剛度矩陣表達如下：

其中ζ0是腐蝕層和完好層分界點處的ζ坐標，h為腐蝕層厚度，t為單元厚度，D（E,）ν為彈性矩陣（腐蝕層對應的材料常數(shù)取Ec，νc），J是雅克比矩陣，B′是殼體局部坐標系下的應變矩陣（其推導過程及表達格式同超參數(shù)殼元[13]）。

3.2 單元等效結點載荷向量

（1）單元等效結點載荷向量計算公式

坑點腐蝕殼體單元的等效結點載荷向量在超參數(shù)殼元的基礎上引入了由腐蝕引起的偏心載荷，其表達公式如下：

式中N為超參數(shù)殼元的插值函數(shù)矩陣[13]，f為作用在單元上的體力，T為作用在單元上的面力。

（15）式中對應的A表達如下：

在求作用于η=±1或ζ=±1上的面力引起的等效結點載荷時，只需要將A中的積分變量作相應的輪換即可。

（15）式中

（16）式中

式中Pappi為單元第i個結點處因載荷偏心引起的附加彎矩，其中Ni（ε,η）為二維插值函數(shù)，ε,η取第i個結點處對應的坐標值，p為作用在ε=±1或η=±1面上的壓力，h為腐蝕層的厚度，w為單元沿ε或η方向的寬度。

（2）單元偏心載荷的推導

假設單元處于單向應力狀態(tài)，由于單元腐蝕層和完好層的剛度差異會導致載荷的不對稱，從而會引起附加彎矩Mapp，如圖7所示。若要保持單元不彎曲，則必須給單元施加一個與附加彎矩等大、反向的彎矩載荷M，此時，單元的腐蝕層和完好層的應變相等，根據(jù)材料的本構關系有：

式中，Ei為單元完好層的彈性模量，Ec為腐蝕層的等效彈性模量，ε為應變。

根據(jù)單元的力平衡，有

式中，σc、σi分別表示腐蝕層和完好層的應力。 聯(lián)立（21）～（22）式可得：

根據(jù)單元的彎矩平衡，有

聯(lián)立（23）～（24）式解得

在求單元等效結點載荷時，將附加彎矩均分到單元壓力面上的各個結點即可。

3.3 單元應力

若單元局部坐標系為x′y′z′（z′取單元的法線方向），則單元在局部坐標系中的應力表達如下：

其中D為彈性矩陣，如（12）式所示，ε′為單元在局部坐標系中的應變。

由于坑點腐蝕會引起應力集中，因此需要對應力作進一步修正。根據(jù)板殼理論，殼體結構中的主要應力為薄膜應力（中面應力）和板應力（彎曲應力），在此薄膜應力主要考慮拉、壓應力（忽略平錯應力），因此，在這里只修正 σx′、σy′。

令 kmx′、kmy′分別為薄膜應力在 x′、y′方向的應力集中系數(shù)，kbx′、kby′分別為彎曲應力在 x′、y′方向的應力集中系數(shù)。則應力修正結果如下：

式中，kmx′σx′mid和 kmy′σy′mid表示對薄膜應力的修正，和?表示對彎曲應力的修正，其中 σx′mid、σy′mid分別表示 x′、y′方向的中面應力，σx′top、σy′top分別表示 x′、y′方向的殼體頂部應力，t為殼體的厚度。

4 算 例

設帶坑點腐蝕的潛艇耐壓圓柱殼半徑R=250cm，殼板厚度t＝2.1cm，肋骨采用球扁鋼20a，肋骨間距 l=60cm，計算載荷 p＝4.26MPa，鋼材屈服極限 σs＝588MPa，腐蝕面積 A＝10×10cm2，坑點腐蝕半徑 r＝0.6cm、深度h＝0.3cm。根據(jù)腐蝕的分布位置將算例分成兩組：

（A）腐蝕區(qū)域位于跨中，腐蝕強度DOP＝15%，

（B）腐蝕區(qū)域位于肋骨，腐蝕強度DOP＝10%。

坑點腐蝕圓柱殼的有限元模型分別如圖8～10所示，考慮到模型的對稱性，計算模型沿周向取1/4，算例A中模型沿縱向取兩個肋距，算例B中模型沿縱向取一個肋距。

分別采用坑點腐蝕殼體單元（Pitting Corrosion Shell element,PCSE）和ANSYS實體單元（SOLID45）進行計算，得出潛艇耐壓圓柱殼的最大橫剖面正應力、最大縱剖面正應力和基于第四強度準則的最大等效應力的計算結果，如表1所示。

表1 坑點腐蝕殼元（PCSE）和ANSYS實體單元（SOLID45）的對比結果Tab.1 The comparison results of PCSE and SOLID45

由表1可知，除了算例A中坑點腐蝕單元和實體單元關于潛艇耐壓圓柱殼的最大等效應力的計算誤差偏大以外，其它特征應力的計算誤差均小于5%，滿足工程應用要求。算例A中的誤差較大的原因是坑點腐蝕單元將腐蝕區(qū)域的應力狀態(tài)處理成雙向應力狀態(tài)，且最大等效應力對應的主應力滿足；而坑點腐蝕區(qū)域的實際應力狀態(tài)為三向應力狀態(tài)，經(jīng)實體單元計算發(fā)現(xiàn)，最大等效應力對應的主應力滿足，σ2<<σ1，σ3<<σ1（對于算例 A：σ1＝－892MPa、σ2＝－38 MPa、σ3＝－4MPa），于是，根據(jù)第四強度準則有

綜上可知，坑點腐蝕殼體單元與ANSYS實體單元的計算結果基本吻合，表明坑點腐蝕殼體單元的計算結果是可信的。

5 結 語

通過本文的研究工作，取得如下研究成果：

（1）建立了坑點腐蝕的力學模型。將單邊坑點腐蝕殼板沿厚度方向分為腐蝕層和完好層兩層，對腐蝕層進行力學等效，求出了腐蝕層的等效彈性模量和等效泊松比；開展了坑點腐蝕的應力集中分析，坑點腐蝕殼板的應力集中分為薄膜應力集中和彎曲應力集中，分別給出了薄膜應力集中系數(shù)和彎曲應力集中系數(shù)的計算公式。

（2）以超參數(shù)殼元為基礎推導了坑點腐蝕殼體單元。導出了坑點腐蝕殼體單元的剛度矩陣和等效結點載荷向量的有限元表達格式，單元剛度矩陣通過沿殼體厚度方向的分段積分求得，在積分時腐蝕層的材料常數(shù)取等效的材料常數(shù)；在計算單元等效結點載荷向量時考慮了由腐蝕引起的偏心載荷；在求解單元應力時考慮了應力集中的影響，并對應力進行了相應修正。

（3）以坑點腐蝕殼體單元為基礎開展了算例分析，并采用ANSYS實體單元的計算結果作為依據(jù)進行了對比，對比結果表明坑點腐蝕殼體單元的計算結果是可信的。

本文的坑點腐蝕殼體單元是基于線彈性理論推導的，該單元只適用于殼體結構的線彈性分析，實際結構問題多數(shù)屬于非線性（材料非線性和幾何非線性）問題，此時線彈性理論不再適用，因此在后續(xù)研究工作中應當考慮非線性因素。另外，由于殼體的實際蝕坑形態(tài)并非規(guī)則的圓柱形，因此本文的研究成果要應用到工程實際還需要更進一步、深入的研究工作。

[1]TSCF.Experimental and theoretical investigation of the strength of corroded hull elements[R].Project 300,Report No.84-3438,Tanker Structure Co-operative Forum,1984.

[2]Daidora JC,Parente J,Orisamolu IR,Ma KT.Residual strength assessment of pitted plate panels[R].SSC-394,Ship Structure Committee,1997.

[3]Flaks V Y.Correlation of pitting corrosion of aluminum plates and reduction of load-beating capacity under tension[J].Fiziko-Khimicheskaya Mekhanika Materialov,1987,14(1):89-93.

[4]Paik J K.Ultimate strength of ships time-variant risk assessment of aging ship staking account to general/pit corrosion,fatigue cracking and local dent damage[R].Technical Report RD 2002-11,American Bureau of Shipping,November 2002.

[5]Paik J K,Wang G,Thayamballi A K,Lee J M.Time-dependent risk assessment of aging ships accounting for general/pit corrosion,fatigue cracking and local denting[C].SNAME Trans,2003,111.

[6]Paik J K,Lee J M,Ko M J.Ultimate shear strength of plate elements with pit corrosion wastage[J].Thin-Walled Structures,2004,42:1161-1176.

[7]Dunbar T E,Pegg N,Taheri F,Jiang L.A computational investigation of the effects of localized corrosion on plates and stiffened panels[J].Marine Structures,2004,17:385-402.

[8]Nakai T,Matsushita H,Yamamoto N,Arai H.Effect of pitting corrosion on local strength of hold frames of bulk carriers(1st report)[J].Marine Structure,2004,17:403-432.

[9]Nakai T,Matsushita H,Yamamoto N,Arai H.Effect of pitting corrosion on local strength of hold frames of bulk carriers(2nd report)-lateral distortional buckling and local face buckling[J].Marine Structure,2004,17:612-641.

[10]王燕舞,黃小平,崔維成.船舶結構鋼海洋環(huán)境點蝕模型研究之一:最大點蝕深度時變模型[J].船舶力學,2007,11(4):577-586.

[11]王燕舞,吳曉源,張雨華,黃小平,崔維成.船舶結構鋼海洋環(huán)境點蝕模型研究之二:實船蝕坑形態(tài)與徑深比時變模型[J].船舶力學,2007,11(5):735-743.

[12]王燕舞,崔維成.考慮腐蝕影響的船舶結構可靠性研究現(xiàn)狀與展望[J].船舶力學,2007,11(2):307-320.

[13]王勖成.有限單元法[M].北京:清華大學出版社,2003:406-413.