同济大学机械与能源工程学院 陈 石 秦朝葵

边界元法介绍及其在阴极保护领域应用

同济大学机械与能源工程学院 陈 石 秦朝葵

数值模拟方法是对在实验研究方法遇到困难或难以展开时的补充。在对常见数值模拟方法的比较之后，得出边界元法在阴极保护领域的优势。边界元法作为一种新的数值模拟方法，其数学模型的建立需要一定合理的假设。数学模型的求解需要一定的边界条件，这些边界条件可分为三类，每一类别在相关环境中对应着一定的情况。在结果呈现中，由于参考电位的选取在不同行业有着不同的规定，因此在对相关标准进行对照比较时应搞清不同参比电极之间的电位关系并进行相应的转换。此外，还概括了边界元方法在阴极保护领域的使用现状。

数值模拟 阴极保护 边界元法 边界条件

0 前言

阴极保护技术是最常用、最有效的预防电化学腐蚀所导致危害的方法之一。通常被保护对象的现场环境情况比较复杂，如海洋中行驶的船舶、敷设在土壤中的钢质天然气管道等，一方面现场检测较为困难，另一方面由于船舶的形状不规则、天然气管道穿越不同性质的地形，影响因素多等因素，难以通过实验或现场测试的方法进行全面而详细的研究。作为实验研究的补充，数值模拟方法在阴极保护领域的应用研究就显得十分重要。边界元方法是一种上世纪60年代发展起来的新的数值模拟方法，已逐步应用在阴极保护的数值模拟领域。特别是天然气管道的阴极防护数值模拟仿真近年来受到很大的关注。

1 几种阴极保护数值模拟计算方法介绍与比较

常用的数值模拟方法包括有限差分法(FDM)、有限元法(FEM)和边界元法(BEM)。有限差分法是最原始的数值模拟方法，其用线性离散方法来计算函数的导数；有限元法是对有限差分法在三维复杂情况下的改进，但是在无限大或半无限大区域求解时，有限元法需对整个区域进行求解来得到目标区域的函数值，因此网格的数量和计算量较大。而边界元法仅需对目标区域的边界进行网格划分。图1和图2示出了边界元法与有限元法的对比。

图1 边界元法的网格划分示意

图2 有限元法的网格划分示意

前者对无底圆柱内部及周围无需进行网格划分，网格全部集中在圆柱的表面上；后者则需对整个圆柱体进行网格划分，图中显示了其中的三个截面，其内部还存在大量的网格。另外，若考虑周围的环境变化，还需要对其周围进行网格划分。

在阴极保护(Cathodic Protection，CP)模拟计算领域，FDM和FEM的共同特点是必须对全部区域进行网格划分，欲获得被保护体表面的电位分布，首先必须计算出电解质内部的电位分布，因此计算量和数据准备量巨大、计算精度不高是这两种方法处理该领域复杂三维空间问题的共同缺陷。BEM则因克服了FDM和FEM数据准备量大的缺陷，成为阴极保护辅助设计中最具广阔前景的数值计算方法。

2 边界元法的数学模型

2.1 边界元法的假设

数值模拟需要对实际事物及其所处环境进行合理的假设从而建立简化的数学模型进而求解，BEM也不例外，它需对研究对象电化学性质以及所处工作环境进行一定的简化，假设内容如下：

(1)电解质为均匀介质或分区域均匀介质：以天然气管道埋设的土壤环境为例，土壤的电导率受到不同环境的影响而有所变化。不同季节、不同温度、不同深度、不同含水量的土壤，其电导率均会有变动，但总体来说，对于具有相同性质同一分块区域的土壤，电导率变化是比较小的并且是各项同性的，他们对阴极保护电位的影响甚微，故对电解质或者同一分块区域的电解质，可假设电导率 k为常数。在BEASY软件中，同一Zone或同一Layer的电导率参数k为单一值，若要反映不同的土壤性质则要将土壤分为不同的Zone或Multi-layer，从而对不同区域进行电导率的定义。

(2)电位场为稳态场：阴极保护领域中，需考察的主要是长期极化后阴极表面的电位分布，其随时间变化极小，可以认为该电位场为稳态电位分布场。对于非稳态计算，需要根据时间输入相应的变化曲线，如极化曲线、钙沉积曲线等，但此类变化的时间单位为年，且经常将时间间隔定为 1年及以上，因此计算时可作为稳态问题来进行考虑。

(3)模型遵从欧姆定律：由电化学腐蚀的原因和基本特点可以得知此假设是合理的。电解质中两点的电位差与两点之间的电流密度值存在线性关系。

(4)模型区域内电解质服从电中性原理：对于电解质中特定的闭合曲面，流进的电流大小等于流出的电流大小，其本身不储存或者释放电荷，电解质中任何一处的净电荷为零。即：电介质中的电位场为无源场。

2.2 数学模型

对均匀的各项同性介质，由于电解质(土壤)中呈电中性，电流的流动遵循拉普拉斯方程：

要获得式(1)在特定情况下的解还需要两个限定条件。其中一个是普适性的，即电流和电位之间的关系服从电工学上的欧姆定律：

式(1)、(2)中：

u——电位，mV；

j——电流密度，mA/m2；

k——电导率，S/m；下文相同。

另外一个即为边界条件，可以分为三种类型，笔者分别将其归纳为第一类、第二类和第三类边界条件。

2.2.1 第一类边界条件

给定边界的电位值：

对于无外加电流的阴极保护系统，阳极与阴极管道的电位即为保持其本身自然电位不变的恒定电位。

对于无穷远处电位假定为零，即u=0(单位：mV)，属于给定边界电位值的情况。

另外，在恒电位控制的外加电流阴极保护系统中，被控制点电位的边界条件也属于此类。

2.2.2 第二类边界条件

给定边界的电流密度值：

对于一般的外加电流阴极保护系统中的辅助阳极，输出的阳极电流是定值i=Const(单位：mA)。

对于特定的系统及特定形状的辅助阳极，相应的电流密度值为定值，即j=i/A=Const (单位：mA/m2，其中 A为辅助阳极表面面积，单位：m2)。在相关软件中也可直接将阳极的边界条件设置为一定的电流密度值，即j=Const (单位：mA/m2)。

在土壤表面，空气的电阻率无限大，其电导率相对于土壤可视为无穷小，因此地表面可作绝缘处理，其电流密度值为零，即j=-k▽u=0(单位：mA/m2)。一般可将地表面处理为几何学上的对称面(Symmetry plane)，在计算时表示该面的法向通量(电流密度值)为零。另外对于土壤电解质中一般的绝缘面也符合该类条件。

另外，无穷远处法向电流密度亦假定为零，即j=-k▽u=0(单位：mA/m2)，属于给定电流密度值的情况。

2.2.3 第三类边界条件

给定边界的电位与电流密度之间的关系：

所谓金属材质的“极化曲线”即属于该种情况，在BEASY软件中默认采用的该类边界条件即为软件自带，一般是反应材料极化电位与法向电流密度之间关系的实测曲线。

另外在以往的相关研究中，还有一些对于该类边界条件的处理方法[3][7]包括：假定极化电位不受电流密度的影响；假定极化的电位与法向电流密度存在线性函数关系，即j=au+b(其中a，b为常数，j，u为电流密度与电位，下同)；基于相关理论确定极化电位与电流密度之间的函数关系式j=f(u)(包括一些经验公式)。对于这几种情况，均可通过人为地在有关软件中编辑材料的极化曲线进而实现。

边界元方法的求解是用格林公式对微分方程积分后再离散处理，最终通过求解关于电位和电流密度的线性方程组系统。可表示为

式(3)中：H、G分别为通过边界元方法定义的系数矩阵；E、I分别为节点电位矩阵和电流密度矢量，其矩阵中每一项分别为待求的电位和电流密度。

3 边界元方法及相关软件的使用

对于边界元法使用在阴极保护领域，国内外也做了相关的探索，其中包括指导性的介绍、针对电化学腐蚀本身的使用以及针对特定领域或特定技术(包括油气管道的阴极保护)的使用。

侯志强、邢少华为了准确评价被保护结构物的保护效果，概括了使用BEASY软件所需要的外加给定条件。在阴极保护系统中这些外加条件一般包括：牺牲阳极以及被保护结构材料的极化特性曲线、涂层状态、电解质的电导率或电阻率以及被保护结构的网格划分形式和方法等。

I.A. Metwally等人对阴极保护干扰进行分类，并用边界元软件研究了不同参数变化对每一干扰类别的影响并对其进行分析。

邢少华、李相波等运用基于BEM算法的软件计算了以铁作为牺牲阳极对不同材质、型号的铜管材的有效保护距离大小，并用实际测量结果验证了软件的计算结果，证明铁阳极在350 mm范围内对铜质管路起到良好的牺牲阳极保护作用。

J.X.Jia等人用对数和线性边界条件的边界元程序来研究电腐蚀问题，其中线性条件运用线性或分段线性方法实现，对数条件运用分段线性方法实现。

姜润翔、胡英娣、龚沈光利用应用边界元法的软件计算出船舶电场信号作为仿真数据对基于点电源的静电场深度换算方法的有效性进行了检验。

刘福国、马桂君、李响对小型导管架式海洋平台外加电流阴极保护系统进行数值模拟分析研究，并进行物理模型实验。

姚华介绍采用数值模拟软件系统得到地下管线及储罐罐底的阴极保护电位和电流密度分布及其干扰因素影响的计算结果。

I.A. Metwally等人用现场测试以及边界元法模拟评估脉冲阴极保护系统的有效性和适用性，从而寻找找出了杂散电流的源头，并提出了降低腐蚀电流的实用性措施。

以BEASY为例，边界元法的最终计算结果可以视觉呈现出电极表面垂直方向的电流密度、被保护装置的保护电位以及各点的电位值。其中电流密度值即为绝对值的大小，保护电位及电位值为相对值，其大小离不开参考电位的选取，保护电位即为阴极电位与其邻近电解质电位之差，而电解质的电位为其各点与参比电极的电位差。参比电极的使用已扩展至海洋、淡水、污水、熔盐、高温油气、土壤等环境。其中Ag/AgCl电极被认为是高温高压电化学试验中最耐用的参比电极，也被广泛应用在船舶海洋电化学腐蚀领域。

对于在埋地燃气管道电化学腐蚀的检测与研究领域，比较常用的是 Cu/CuSO4参比电极，有关标准的规定值也是针对 Cu/CuSO4参比电极设定。然而在 BEASY软件中，参考电位选择的是Ag/AgCl(海水)电极的电位。因此了解不同参比电极之间的电位关系显得非常必要，图3示出了常用的参比电极电位相对于标准氢电极的大小。

图3 常用参比电极相对于标准氢电极的电位值(25℃)

利用该软件计算出的电位结果(如图4)对有关标准进行参照时要进行一定的换算，即电位的大小要减去Cu/CuSO4参比电极与Ag/AgCl参比电极电位之差0.05 V(50mV)。

目前国内曾将边界元法用于舰艇阴极保护系统或船舶、海洋平台方面的仿真优化设计，而在输油气管道上的研究仍十分有限，其大部分研究仅限于实验室模拟。

根据相关软件的帮助文件，该类软件在油气管道阴极保护领域可实现的功能包含两个方面：预测管道阴极保护系统性能，主要为预测保护等级、评估各设计变量的重要性、预计牺牲阳极寿命、评估系统在不同环境和设计条件下的性能；预测不同系统间的相互影响，主要为不同来源杂散电流腐蚀、已有牺牲阳极对现有管道设备的干扰等。

图4 BEASY电位计算结果示意

4 结语

边界元法有着完整的数学模型和计算理论体系，凭借其在无限大或半无限大介质分析中计算量方面的优越性，在阴极保护领域可以得到很好地应用。

边界元法及其相关软件对于现在受轨道交通杂散电流影响天然气管道的数值模拟能够很大程度上简化预测和监管的工作量，带来一定的经济性。

有关数值模拟软件针对管道的建模提供了简化的工具和方法，减小了使用人员的工作量，但在使用时应注意参考电位的选取和大小关系的数量转换。

Introduction of Boundary Element Method and Its Application in Cathode Protection

School of Mechanical Engineering Tongji University Chen Shi Qin Chaokui

Numerical simulation is the supplement of experiments in meeting difficulties, and advantages of the boundary element method (BEM) in cathode protection (CP) have been found in comparing with other common numerical simulation methods. The mathematical model has been introduced, and the present situation of the application of BEM software in CP has been discussed as well.

numerical simulation, cathode protection, boundary element method, boundary condition