基于边界元法的站场阴极保护数值模拟软件设计与开发

廖柯熹,马晓磊,张中放,周卫军,张 瑶,任驰远(. 西南石油大学 石油与天然气工程学院，成都 60500； . 中国石油天然气集团公司 塔里木油田分公司油气运销部，库尔勒 84000； . 西南石油大学 理学院，成都 60500)

基于边界元法的站场阴极保护数值模拟软件设计与开发

廖柯熹1,马晓磊1,张中放2,周卫军2,张 瑶2,任驰远3
(1. 西南石油大学 石油与天然气工程学院，成都 610500； 2. 中国石油天然气集团公司 塔里木油田分公司油气运销部，库尔勒 841000； 3. 西南石油大学 理学院，成都 610500)

边界元法在阴极保护电位分布计算方面有着明显的优势，采用边界元法对数值模型边界进行离散，计算节点少，计算速率快。站场阴极保护数值模拟软件基于边界元法对站场阴极保护电位分布进行计算，得到可视化结果，通过该结果可以发现屏蔽和干扰位置，对辅助阳极位置进行优化，从而解决站场阴极保护系统腐蚀区域不能预测的困难，使站场阴极保护系统工程更加合理有效。

边界元法；站场；阴极保护；电位分布模拟

20世纪60年代国外开始进行区域性阴极保护的研究和应用，我国从20世纪70年代末、80年代初开始在油田和部分输油站尝试采用区域性阴极保护技术，到90年代中期该技术已相对成熟并开始推广应用。2001年11月，鄯乌输气管道，鄯善首站成为第一个实施区域性阴极保护的压气站[1]。

目前，现场阴极保护施工过程中普遍采用试凑法，保护效果不太理想，往往需要进行后期整改。随着数值模拟技术在阴极保护领域的深入应用，使得复杂环境中阴极保护效果的预测和低成本分析研究成为可能。与传统方法相比，数值模拟技术在管道和站场阴极保护各种影响因素的干扰趋势和规律研究方面具有无法比拟的优势，同时可用于阳极位置优化、阴极保护效果的预测和全面评价[2]。

目前，可以进行区域阴极保护系统电位分布计算的软件主要是基于有限元法和边界元法进行计算，基于有限元法进行计算的软件有FLUENT、COMSOL等，网格划分复杂，计算节点多，不适用于现场情况。而BEASY是基于边界元法进行计算的软件，在进行网格划分时，只需对边界进行划分，计算量小，可以节省计算时间。

由于阴极保护计算通常只关心材料外表面和地表的参数分布情况，相较于有限元法，边界元法具有只对边界进行离散化的特点，更适用于计算阴极保护等与腐蚀相关的问题。为此，本工作设计开发了基于边界元法的阴极保护电位分布计算软件——站场阴极保护数值模拟软件(Numerical Simulation Software Of Cathodic Protection，简称NSSCP)，以期可靠分析与评价站场阴极保护系统的有效性。

1 阴极保护电位数值计算模型

埋地管道阴极保护电位控制方程和储罐底板阴极保护电位控制方程的推导过程类似，在此仅考虑储罐底板阴极保护的电位控制方程，可以用Poisson方程表示为

2φ

式中：δ(X-Xe)表示以阳极点X(x,y,z)为中心的狄拉克函数，其性质见式(2)。

式中：φ是施加阴极保护后的极化电位值；x，y和z分别是计算点与储罐底板中心的距离；Ie为阳极点Xe(x,y,z)输出的电流(即点源强度)。当所研究的恒定电场中存在恒定的电源，电位满足poisson方程，在给定边界条件时，可通过求解该方程得到电位的分布。

公式(1)仅用来描述储罐底板阴极保护的电位分布，但是如果想要求得具体解，就需要附加边界条件。数值模型的边界为Ω(研究区域)和总边界Γ(Γ表面包围3 Ω),总边界Γ包括无限大地面表面边界Γd、储罐底板外表面边界Γc、虚设半无限域土壤球冠边界Γ∞。

边界条件简化：地表面边界Γd与大气相接触，外加阴极保护电流只能在土壤介质中流动，不能流入大气，所以把地表面当作绝缘面；假设无穷远处的电位为零，忽略阳极对半无限大土壤表面的影响，所以在无穷远处表面的电流密度也为零；阴极的边界条件采用描述电流密度J和电位φ之间函数关系的极化函数来定义

J=f(φ

式中：φeq是储罐底板在土壤环境中的自腐蚀电位，f(φ-φeq)为阴极极化函数。

综上所述得到数学模型及边值问题如下

2 边界积分方程推导与处理

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

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

2.1 边界积分方程推导

以P′为圆心，以微小ε为半径，在边界上建立微小半球模型，设Q′点位于边界上。通过积分方程式把域内的点用P、Q表示，边界上的点用P′、Q′的位势值φ以及通量值q联系起来，对区域Ω内任意点均成立。在阴极保护电位分布模型中，边界积分方程式为

C(P′)φ(P′)

所研究的模型属于半无限域，很难在虚设半无限域土壤球冠边界Γ∞和无限大地面Γd表面边界上划分网格。因此根据模型已知的边界条件，对模型的边值问题进行简化。在边界Γ∞处，满足函数值和函数的法向导数都为零，即：

代入后可以发现，不需要在边界Γ∞上划分单元和计算。在边界Γd上：

代入后不能完全消除该边界的积分项，因此需要对边界Γd进行单元划分。由于Γd是边界无限大地面表面，不能确定应该取多大的边界来划分单元，同时划分单元后也会增大计算量。通过镜像法可在Γd边界，寻找一个特定的基本解φ*，使边界Γd无需划分。

根据Poisson方程的基本解得到的模型为

φ*(b)=φ*(A,B)+φ*(A′,B′)=

根据模型基本解，边界积分方程可简化为

选用模型的基本解时，无需对边界Γd划分单元，地表面的边界条件自动满足，所以只需对阴极的边界进行划分。

2.2 边界离散

边界元法将边界上的积分方程式离散化，并归结为代数方程组求得近似解。通过插值函数得到各个边界单元上的函数值和函数的法向导数值。

按单元求和的离散形式改写为按节点求和的边界积分方程离散格式，即

由于j点也是边界上的点，所以必然与某I节点为同一点，所以合并同类项后，公式简化为

式(14)中:

对于m个节点，可以得到联立的一次方程组为

式中:H与G表示系数矩阵，Φ表示边界单元节点处的函数值的列向量，Q表示边界单元节点处函数法向导数值的列向量，按照已给定的边界条件，由式(16)可以计算出边界上节点的未知函数值与其法向导数值。

2.3 辅助阳极处理

根据δ函数的性质可知，域内存在一个集中源点X1时

式中:φ*(P′,A1)是对应A1的基本解，得

当域内存在多个集中源点An(n=1LN)时，得

φ*(P′,An)b(Q)dΩ(Q)=

积分项Bj不会产生未知量。

2.4 非线性边界条件处理

在数值模拟过程中，边界条件的确定及其处理方式直接影响计算结果的准确性，软件总共有4种边界条件，包括极化曲线边界、电位边界、电流边界、阻抗边界，其中最主要的边界条件是非线性极化曲线边界。

由于阴极和阳极的电位和电流密度的关系并非简单的线性关系，而是非线性关系，在模型处理过程中必须对其进行处理。储罐底板表面受到多种不同因素的影响，极化函数不能只用单独的环境或者材料的影响函数来表达，为得到极化函数定义的阴极边界条件J=f(φ-φeq)，通过试验测得极化曲线作为边界条件是非常有必要的。

通过对测得的极化曲线进行拟合，可以确定电流密度与极化电位之间的定量关系。主流的拟合方法有Newton-Raphson迭代法和分段线性拟合法。

由于被保护金属表面受到多种不同因素的影响，所以极化函数不能只用单独的环境或者材料的影响函数来表达，为得到极化函数定义的阴极边界条件，综合考虑现场实际状况，通过试验测得被保护金属在土壤中的电位分布，为了取得良好的收敛性和计算精度，该软件采用分段拟线性化法处理非线性边界条件。拟合得到电流密度J与电位φ的之间的线性关系

J=aφ

推导得

将式(22)代入下列矩阵方程

移项整理得

矩阵最终化为AX=F的形式，辅以迭代的方法，可以求解得到电位值大小。

3 软件开发

基于边界元算法开发了“站场阴极保护数值模拟软件(NSSCP)”，该软件可以实现三维站场阴极保护系统的模拟，获得可视化电位分布规律，从而解决现有站场和新改扩建站场内埋地金属构筑物密集所造成的屏蔽和干扰问题，并防止站场内阴极保护电流对未纳入阴极保护管道的干扰影响。相较于有限元法计算,该软件的优点就是简单快捷。

3.1 软件开发流程和计算流程

在软件的设计开发过程中，使用Access 2007数据库作为后台数据库，并基于Visual Studio C# 2010设计了软件用户界面。软件共划分出阴极保护系统参数存储、边界条件处理、电位计算及数据导出四个主模块。图1和图2分别是软件开发的流程图和软件计算流程图。

图1 软件开发流程Fig. 1 Software development process

图2 软件计算流程Fig. 2 Software calculation process

3.2 阴极保护系统参数存储和建模

用户可以在“工程参数”界面输入工程名称、站场名称、土壤电阻、站场简介、模拟目的、精度选择(高精度、一般精度、较低精度)、精度参数(最大迭代次数和收敛误差阈值)等基础数据。土壤电阻率会随地域、季节和天气发生变化，在设置土壤电阻率时，应考虑这些因素的影响。在“实体对象”界面中可以进行腐蚀实体对象建模，储罐、阳极、埋地管道需要设置保护层系数、直径、埋深、位置坐标、边界条件等参数。

3.3 电位分布模拟计算

采用对于求解区域进行自动剖分网格的方法来划分单元，将形成的边界单元的坐标保存在相应的输入文件中。计算系数矩阵，调用边界条件，求解线性方程组，完成迭代。

软件在使用时最重要的步骤就是进行参数收集与导入，在软件的主界面应先进行阴极极化曲线边界参数的导入；其次是土壤电阻率的导入，需要考虑不同季节导致的土壤电阻率的变化，根据实测结果进行导入；最后导入的是实体腐蚀对象的参数(包括地理位置、几何尺寸等)，通过一系列的计算便可以得到电位分布。

干扰分析菜单栏可以实现对计算结果查看和分析，得到整个站场节点电位分布数据列表,见图3。

干扰问题是站场阴极保护系统实施的一个难题。

图3 数值计算结果查看界面Fig. 3 View interface of numerical results

对油气输送站场进行阴极保护，由于是将整个站区内的所有地下金属结构全部纳入保护系统，因此一般不会产生内部干扰，但是由于站场阴极保护电流通常远大于干线保护电流，故常常对管道干线等外部结构及其阴极保护系统造成干扰。一般系统电流输出越大，干扰也会越严重。通过该软件干扰分析产生的干扰数据可以有效判断出是否发生发生干扰，提高了站场阴极保护系统的可靠性。

4 实例验证

某站场使用一套阴极保护系统对四个储罐进行保护，6组深井阳极埋设45 m分别位于储罐的两侧每组阳极由8支研究块组成,共8 m,直径为25 cm，储罐直径为60 m，恒电位仪的输出电位为8 V。对储罐周边测试桩进行检测，测得电位均处于-0.85～-1.2 V(相对铜/硫酸铜参比电极)，达到保护状态，保护效果良好。

通过采用COMSOL软件模拟结果进行比对分析，模拟时，所取土壤电阻率为50 Ω·m，砂垫层200 Ω·m，沥青砂层1 000 Ω·m。采用COMSOL软件模拟得到的储罐保护电位分布云图，通过模拟计算，并将COMSOL软件模拟结果和实测数据进行统计,见表1。

表1 计算结果比较Tab. 1 Comparison of calculated results

COMSOL基于有限元方法进行计算，计算时间较长，一般认为其结果较为精确，但是存在外边界截断误差；而NSSCP是基于边界元法进行开发的，不存在边界截断误差。由表1可见，NSSCP的计算精度接近COMSOL的，证明该软件计算结果可靠，可以作为站场阴极保护系统的检测维护软件。

5 结束语

(1) NSSCP软件在计算结果精度方面接近COMSOL软件，并且计算过程中收敛速率快，计算时间短，提高了使用效率。

(2) NSSCP软件具有较高的实用性，实际应用效果较好，能够用于指导站场阴极保护系统的设计计算和科学维护。

[1] 陈洪源,范志刚,刘玲莉,等. 区域性阴极保护技术在输气站场中的应用[J]. 油气储运,2005,24(5):41-44.

[2] 张丰,陈洪源,李国栋,等. 数值模拟在管道和站场阴极保护中的应用[J]. 油气储运,2011,30(3):208-212.

[3] 陈石,秦朝葵. 边界元法介绍及其在阴极保护领域应用[J]. 上海煤气,2015(3):7-10.

“2017材料物理测试新技术研讨会”征 文 通 知

由上海材料研究所等单位联合主办，《理化检验-物理分册》编辑部、《机械工程材料》编辑部共同承办的“物理测试新技术研讨会”已连续成功举办两届，受到了国内广大材料物理测试专家、学者及从业人员的热烈欢迎和积极参与。在此基础上，“2017材料物理测试新技术研讨会”将继续围绕“新的测试技术与方法”、“新的测试仪器和设备”、“新的试验标准与规范”等方面展开讨论，旨在为广大物理测试工作者以及设备厂商提供一个沟通、交流和展示的平台，以总结经验，促进创新，进而推动我国材料物理测试技术的发展与进步。本次会议将于2017年10月底在上海举行，由上海材料研究所、中国机械工程学会理化检验分会、中国机械工程学会材料分会联合主办，《理化检验-物理分册》编辑部、《机械工程材料》编辑部和上海市工程材料应用与评价重点实验室共同承办。会议将邀请业内权威专家作专题报告，同时热忱欢迎行业专家、学者及从业人员踊跃投稿，积极参会。所征论文将以《理化检验-物理分册》正刊专栏或增刊形式发表，评选出的优秀论文将在会议期间进行口头交流。

1 征文内容

(1) 材料物理测试的新技术、新方法、新标准和新仪器设备的介绍及应用；

(2) 物理测试技术在新材料开发、研究中的应用；

(3) 物理测试技术在产品质量控制过程中的应用；

(4) 物理测试技术在材料表面处理工艺、材料冷热加工工艺、材料改性技术以及先进制造技术方面的应用；

(5) 物理测试技术在产品在役安全监控、寿命估算、可靠性评估等方面的应用；

(6) 机械产品及其构件失效分析的最新思路、技术、方法及典型的失效分析案例介绍；

(7) 试验设备的升级改造技术；

(8) 试验室的管理及在全面质量管理中的成功经验；

(9) 其他与材料物理测试技术相关的论文。

2 征文要求

(1) 未公开发表过的原创性论文。

(2) 来稿要求论点明确、数据可靠、逻辑严密、文字精炼。每篇论文必须包括题目、作者姓名、作者单位、单位所在地及邮政编码、摘要和关键词、正文、参考文献和第一作者及通信作者简介(包括姓名、性别、职称、出生年月、所获学位、目前主要从事的工作和研究方向)，在文稿的首页地脚处注明论文属何项目、何基金(编号)资助，没有的不注明。论文题目、作者姓名、作者单位、单位所在地及邮政编码、摘要和关键词应中英文对照。

(3) 论文摘要尽量写成报道性文摘，包括目的、方法、结果、结论4方面内容(200字左右)，应具有独立性与自含性，关键词选择贴近文义的规范性单词或组合词(3～8个)。

(4) 文稿篇幅(含图表)一般不超过8 000字。文中量和单位的使用请参照中华人民共和国法定计量单位最新标准。外文字符必须分清大、小写，正、斜体，黑、白体，上下角标应区别明显。

(5) 文中的图、表应有自明性，图像要清晰，层次要分明；所有图和表均应标明中英文图题和表题。

(6) 参考文献的著录格式采用顺序编码制，请按文中出现的先后顺序编号。所引文献必须是作者直接阅读参考过的、最主要的、公开出版文献。未公开发表的、且很有必要引用的，请采用脚注方式标明，参考文献不少于5条。

3 征文截止时间

2017年8月31日。

4 论文提交

论文投稿至E-mail: pt@mat-test.com，请注明稿件为“物理测试会议论文”。

5 联系方式

地址： 上海市邯郸路99号，《理化检验-物理分册》编辑部、《机械工程材料》编辑部。

联系人：李玲；金静静。

电话：021-65559079；021-65556775-361。E-mail：pt@mat-test.com。

(会议的具体时间与地点，及会务费用等详情将在第二轮通知中发出。)

Design and Development of Numerical Simulation Software for Cathodic Protection of Stations Based on Boundary Element Method

LIAO Kexi1, MA Xiaolei1, ZHANG Zhongfang2, ZHOU Weijun2, ZHANG Yao2, REN Chiyuan3
(1. College of Petroleum & Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China;2. Tarim Oilfield Company Oil and Gas Distribution Unit, China National Petroleum Corporation, Korla 841000, China;3. College of Science, Southwest Petroleum University, Chengdu 610500, China)

Boundary element method in distributed computing of cathodic protection potential has a distinct advantage. When the boundary element method to scatter a numerical model was used, the calculation nodes were less and calculation speed was fast. The software based on the boundary element method calculated the cathodic protection potential distribution in cathodic station, and the results were visual, shielding and interference position could be found through the results, auxiliary anode position could be optimized in order to resolve the difficulty that the station cathodic corrosion protection system area could not be predicted, so the station cathodic protection system engineering could be more rational and effective.

boundary element method; station; cathodic protection; potential distribution simulation

2016-01-15

廖柯熹(1970-)，教授，博士，从事油气储运系统完整性管理与安全评价工作,liaokxswpi@163.com

10.11973/fsyfh-201707015

TE88

A

1005-748X(2017)07-0551-06