含坑点腐蚀的壳体有限元方法

徐 强，万正权

（中国船舶科学研究中心，江苏 无锡 214082）

1 引 言

壳体在工程实际中得到大量应用，例如，航海和海洋工程中的舰船、潜艇和油气管道等都广泛地采用了各种形式的壳体结构。老龄壳体结构多数存在严重的腐蚀问题，坑点腐蚀是常见且危害极大的一种，坑点腐蚀属于局部腐蚀（又称非均匀腐蚀），是由于金属表面物理和化学的不均匀性所引起的。严重的坑点腐蚀可能使金属壳板穿透，引起流体泄漏、火灾、爆炸等事故，因此，有必要开展坑点腐蚀对壳体结构强度的影响研究。近年来，国内外学者对工程结构坑点腐蚀问题的重视程度在不断上升，并开展了系列相关问题的研究。

TSCF[1]根据坑点腐蚀钢板的弯曲试验与理论推导，提出一个包括弯曲刚度、质量、边界条件、板尺寸规格等参数的腐蚀钢板等效厚度计算公式。

Daidora[2]等人提出用腐蚀坑数据的平均值和最大值或者用腐蚀坑的个数和坑的最大深度来估计结构腐蚀后的剩余厚度，然后通过概率方法提出一种用于评估腐蚀厚度对板局部屈服和板屈曲影响的方法。

Flaks[3]提出一种数学方法，用于评估坑点腐蚀对受拉铝板的极限强度的影响。通过对自然腐蚀铝板的拉伸试验，得到了一个考虑拉伸强度、屈服强度和硬度损失的系数。

Paik[4-5]等通过试验和有限元模拟分析，研究了腐蚀和疲劳等因素对老龄船舶的极限强度的影响。提出以“极限强度破坏因子”来衡量坑点腐蚀引起的剖面面积减少对板构件强度的影响；并认为在评估结构强度时，坑点腐蚀所造成的截面积损失是比坑点腐蚀深度更值得关注的指标。

Paik[6]研究了受坑点腐蚀影响的钢板构件在面内剪切载荷作用下的极限强度特性。针对不同的坑点腐蚀强度和不同的板构件几何特性，对受面内剪切载荷的板构件进行了一系列的非线性有限元分析。通过对有限元计算结果进行回归分析，得到了闭合形式的板构件的极限强度设计公式。

Dunbar[7]讨论了坑点腐蚀对船舶结构中的典型板构件和加筋板构件的影响。用有限元方法对板的初始屈曲、极限破坏和后极限响应进行了研究分析。有限元模型考虑了几何缺陷和残余应力，并把有限元模拟结果与理论计算结果和有效的实验测量进行对比。该文提出的受损船舶结构的完整性理论有助于发展更有效的船舶维修策略。

Nakai、Yamamoto[8]等开展了坑点腐蚀对散货船货舱肋骨局部强度的影响研究。对老龄散货船和油船的腐蚀情况进行了调查统计，并在此基础上分别对真实坑点腐蚀样本和人工模拟坑点腐蚀样本进行了拉伸、压缩试验，考察了不同的腐蚀分布情况和平均腐蚀厚度对结构强度的影响，并进行了有限元模拟分析。

Nakai[9]还研究了局部坑点腐蚀对工型材腹板断裂的影响，研究发现，腐蚀的不同分布情况会显著地影响腹板的断裂行为。通过进一步的有限元分析发现，当坑点腐蚀腹板和均匀腐蚀腹板的平均厚度损失相同时，前者的极限强度要小于或等于后者的极限强度。

王燕舞、黄小平、崔维成[10]开展了船舶结构钢海洋环境点蚀模型之最大点蚀深度时变模型研究。针对我国船舶结构常用炭钢、低合金钢的实海腐蚀试验观测数据建立了新型Weibull函数形式的点蚀最大深度模型。

王燕舞、吴晓源等[11]开展了船舶结构钢海洋环境点蚀模型之实船蚀坑形态与径深比时变模型研究。基于点蚀生长原理与实船蚀坑形态检测数据，对船舶受点蚀构件的两类主要蚀坑的形态、形成机制及其随船龄演化情况进行了描述。依据Yamamoto[8]实测散货船货舱区肋骨蚀坑数据并与点蚀深度模型[10]相结合，建立了相应的蚀坑径深比时变模型。

目前，国内外学者在对坑点腐蚀几何形态的具体处理上还存在较大差异，如：将蚀坑处理成圆柱形、球冠形、圆锥形、矩形槽等[6-9，11-12]，这种差异主要是由腐蚀环境和金属材料的不同所引起的。为了便于展开本文的研究工作，并考虑到蚀坑几何形态的简洁性，本文将坑点腐蚀取成圆柱形。

本文首先建立了坑点腐蚀的力学模型，将单边坑点腐蚀壳板沿厚度方向分为腐蚀层和完好层，求出了腐蚀层的等效材料常数，开展了坑点腐蚀的应力集中分析。然后采用有限元理论、以超参数壳元为基础推导了坑点腐蚀壳体单元，主要导出了单元的刚度矩阵和等效结点载荷向量的有限元表达格式，并在求解单元应力时考虑了坑点腐蚀应力集中的影响。最后开展了基于坑点腐蚀壳体单元的算例分析，并与ANSYS实体单元的计算结果进行了对比。

2 坑点腐蚀的力学模型

2.1 坑点腐蚀壳板的几何模型

坑点腐蚀壳板的几何模型如图1所示，本文考虑的腐蚀形态为单边腐蚀，壳板沿厚度方向分成两层：腐蚀层（B）和完好层（C）。在B中，坑点腐蚀呈圆柱状，沿壳板厚度方向穿透腐蚀层。腐蚀层（B）的厚度（即腐蚀坑的深度）取实测蚀坑深度的平均值。完好层（C）和腐蚀层（B）的弹性模量和泊松比分别为E，ν和Ec，νc，E，ν为材料常数，而后者为力学等效的弹性模量和泊松比。

腐蚀层的等效材料常数Ec，νc定义如下：

式中，P为拉伸载荷，S为腐蚀层沿拉伸方向的横截面积，L、W分别为腐蚀层的长度和宽度，ΔL、ΔW分别为腐蚀层沿拉伸方向和横向的变形。根据（1）式进一步得到无量纲化的等效材料常数：

本文通过有限元数值分析求kE、kν，腐蚀层的有限元模型如图2所示，计算模型分为蚀坑有序分布和蚀坑随机分布两组。

经过有限元数值计算并对结果进行回归分析，得到腐蚀层的等效材料常数回归公式如下：

式中DOP为腐蚀强度，它等于坑点腐蚀的面积与壳板面积的比值，即

其中，Api是第 i个蚀坑的面积，且n为蚀坑的个数，l、w分别为壳板的长度和宽度。

2.2 坑点腐蚀的应力集中分析

坑点腐蚀壳板的应力集中可分为薄膜应力集中和弯曲应力集中，如图4。应力集中系数定义为腐蚀部位的应力峰值与基准应力的比值，基准应力取假设不存在腐蚀时腐蚀部位对应的应力值。

影响坑点腐蚀应力集中系数的因素包括蚀坑深度与壳板厚度之比（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 进行一系列有限元数值计算并对计算结果进行回归分析，得到薄膜应力集中系数回归公式如下：

按照 0≤h/Th≤0.3，0.24≤M1/M2≤0.3，0≤DOP≤0.3，R/H=2进行一系列有限元数值计算并对计算结果进行回归分析，得到弯曲应力集中系数回归公式如下：

坑点腐蚀应力集中系数回归曲线如图5所示。

3 坑点腐蚀壳体单元推导

3.1 单元刚度矩阵

坑点腐蚀壳体单元的刚度矩阵以超参数壳元的刚度矩阵为基础，通过沿单元厚度方向腐蚀层和完好层的分段积分求得。若坐标系取自然坐标ε、η、ζ（如图6所示），则坑点腐蚀壳体单元的刚度矩阵表达如下：

其中ζ0是腐蚀层和完好层分界点处的ζ坐标，h为腐蚀层厚度，t为单元厚度，D（E,）ν为弹性矩阵（腐蚀层对应的材料常数取Ec，νc），J是雅克比矩阵，B′是壳体局部坐标系下的应变矩阵（其推导过程及表达格式同超参数壳元[13]）。

3.2 单元等效结点载荷向量

（1）单元等效结点载荷向量计算公式

坑点腐蚀壳体单元的等效结点载荷向量在超参数壳元的基础上引入了由腐蚀引起的偏心载荷，其表达公式如下：

式中N为超参数壳元的插值函数矩阵[13]，f为作用在单元上的体力，T为作用在单元上的面力。

（15）式中对应的A表达如下：

在求作用于η=±1或ζ=±1上的面力引起的等效结点载荷时，只需要将A中的积分变量作相应的轮换即可。

（15）式中

（16）式中

式中Pappi为单元第i个结点处因载荷偏心引起的附加弯矩，其中Ni（ε,η）为二维插值函数，ε,η取第i个结点处对应的坐标值，p为作用在ε=±1或η=±1面上的压力，h为腐蚀层的厚度，w为单元沿ε或η方向的宽度。

（2）单元偏心载荷的推导

假设单元处于单向应力状态，由于单元腐蚀层和完好层的刚度差异会导致载荷的不对称，从而会引起附加弯矩Mapp，如图7所示。若要保持单元不弯曲，则必须给单元施加一个与附加弯矩等大、反向的弯矩载荷M，此时，单元的腐蚀层和完好层的应变相等，根据材料的本构关系有：

式中，Ei为单元完好层的弹性模量，Ec为腐蚀层的等效弹性模量，ε为应变。

根据单元的力平衡，有

式中，σc、σi分别表示腐蚀层和完好层的应力。 联立（21）～（22）式可得：

根据单元的弯矩平衡，有

联立（23）～（24）式解得

在求单元等效结点载荷时，将附加弯矩均分到单元压力面上的各个结点即可。

3.3 单元应力

若单元局部坐标系为x′y′z′（z′取单元的法线方向），则单元在局部坐标系中的应力表达如下：

其中D为弹性矩阵，如（12）式所示，ε′为单元在局部坐标系中的应变。

由于坑点腐蚀会引起应力集中，因此需要对应力作进一步修正。根据板壳理论，壳体结构中的主要应力为薄膜应力（中面应力）和板应力（弯曲应力），在此薄膜应力主要考虑拉、压应力（忽略平错应力），因此，在这里只修正 σx′、σy′。

令 kmx′、kmy′分别为薄膜应力在 x′、y′方向的应力集中系数，kbx′、kby′分别为弯曲应力在 x′、y′方向的应力集中系数。则应力修正结果如下：

式中，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。根据腐蚀的分布位置将算例分成两组：

（A）腐蚀区域位于跨中，腐蚀强度DOP＝15%，

（B）腐蚀区域位于肋骨，腐蚀强度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中的误差较大的原因是坑点腐蚀单元将腐蚀区域的应力状态处理成双向应力状态，且最大等效应力对应的主应力满足；而坑点腐蚀区域的实际应力状态为三向应力状态，经实体单元计算发现，最大等效应力对应的主应力满足，σ2<<σ1，σ3<<σ1（对于算例 A：σ1＝－892MPa、σ2＝－38 MPa、σ3＝－4MPa），于是，根据第四强度准则有

综上可知，坑点腐蚀壳体单元与ANSYS实体单元的计算结果基本吻合，表明坑点腐蚀壳体单元的计算结果是可信的。

5 结 语

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

（1）建立了坑点腐蚀的力学模型。将单边坑点腐蚀壳板沿厚度方向分为腐蚀层和完好层两层，对腐蚀层进行力学等效，求出了腐蚀层的等效弹性模量和等效泊松比；开展了坑点腐蚀的应力集中分析，坑点腐蚀壳板的应力集中分为薄膜应力集中和弯曲应力集中，分别给出了薄膜应力集中系数和弯曲应力集中系数的计算公式。

（2）以超参数壳元为基础推导了坑点腐蚀壳体单元。导出了坑点腐蚀壳体单元的刚度矩阵和等效结点载荷向量的有限元表达格式，单元刚度矩阵通过沿壳体厚度方向的分段积分求得，在积分时腐蚀层的材料常数取等效的材料常数；在计算单元等效结点载荷向量时考虑了由腐蚀引起的偏心载荷；在求解单元应力时考虑了应力集中的影响，并对应力进行了相应修正。

（3）以坑点腐蚀壳体单元为基础开展了算例分析，并采用ANSYS实体单元的计算结果作为依据进行了对比，对比结果表明坑点腐蚀壳体单元的计算结果是可信的。

本文的坑点腐蚀壳体单元是基于线弹性理论推导的，该单元只适用于壳体结构的线弹性分析，实际结构问题多数属于非线性（材料非线性和几何非线性）问题，此时线弹性理论不再适用，因此在后续研究工作中应当考虑非线性因素。另外，由于壳体的实际蚀坑形态并非规则的圆柱形，因此本文的研究成果要应用到工程实际还需要更进一步、深入的研究工作。

[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]王燕舞,黄小平,崔维成.船舶结构钢海洋环境点蚀模型研究之一:最大点蚀深度时变模型[J].船舶力学,2007,11(4):577-586.

[11]王燕舞,吴晓源,张雨华,黄小平,崔维成.船舶结构钢海洋环境点蚀模型研究之二:实船蚀坑形态与径深比时变模型[J].船舶力学,2007,11(5):735-743.

[12]王燕舞,崔维成.考虑腐蚀影响的船舶结构可靠性研究现状与展望[J].船舶力学,2007,11(2):307-320.

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