一种用于阴极保护数值模拟的高性能计算方法

(1. 中海石油(中国)有限公司 深圳分公司 深水工程建设中心，深圳 518054; 2. 大连理工大学 船舶工程学院，大连 116024)

阴极保护(Catholic Protection，简称CP)是通过电化学方法降低金属腐蚀速率的一种技术[1]，有牺牲阳极阴极保护(Galvanic Anode Catholic Protection，简称GACP)和外加电流阴极保护(Impressed Current Catholic Protection，简称ICCP)。牺牲阳极阴极保护设计[2]，大多根据经验公式确定需要的阳极数量，然后再根据设计人员的经验确定阳极布置方案，因此该方法的误差较大。采用该方法时一般选取较大的安全系数，但这会导致阳极用量增加和经济成本升高，而且增加了安装施工难度和降低结构强度。外加电流阴极保护系统设计[3]，采用的辅助阳极的单体输出电流较大，辅助阳极的数量较少，结构表面的保护电位分布范围较大，且基于规范经验公式的设计方法，往往会导致结构表面发生欠保护或过保护现象，从而降低结构强度。

为了弥补经验公式的不足，科研人员将有限差分法(Finite Difference Method，简称FDM)[4-5]、边界元法(Boundary Element Method，简称BEM)[6-8]和有限元法(Finite Element Method，简称FEM)[9-11]等数值模拟方法应用到阴极保护设计中。在牺牲阳极阴极保护设计中，通过数值模拟可以确定最优的阳极布置方案，即采用最少的阳极达到最佳的保护效果，从而降低经济成本和施工难度。在外加电流阴极保护设计中，也可通过数值模拟确定最优的辅助阳极数量、布置和电流输出方案，避免过保护问题。

英国Beasy公司针对阴极保护数值模拟，已开发出比较成熟的基于边界元法的商业软件。目前在国内，基于边界元法的阴极保护数值模拟也有了较多的研究[12-15]，但还没有成熟的商业软件。采用边界元法进行阴极保护数值模拟，无需考虑域内电位分布，只需要对边界进行离散，且可以直接得到结构外表面的保护电位分布，因此，边界元法是目前最理想的阴极保护数值模拟方法。

然而，随着海洋结构物的结构形式日益复杂，尺寸越来越大，数值模拟的节点规模也越来越大，大规模模型的数值模拟计算效率很低，严重影响了模拟的进度和效率。当前典型海洋工程结构物边界元模型的节点数通常在104～105数量级，对应的双精度型浮点系数矩阵将占用计算机0.7～74.5 GB内存空间，对计算机硬件资源耗用巨大。通过常规串行计算方式生成如此规模的系数矩阵，计算效率很低。边界元数值模拟中，系数矩阵中每行数据的计算相互独立且没有交叉调用，易于实现并行计算，理论上通过并行算法提高的计算效率倍数为计算机的CPU数，因此在CPU数目较多的服务器上通过并行算法可以大幅提升系数矩阵的生成效率。

C#是一种面向对象的高级程序语言，从C和C++语言演化而来[16]，吸收了C、C++与Java语言的优点，具有Visual Basic编译语言的高效率和C++编译语言的强大功能。C#语言开发的系统具有界面友好、易维护和升级等优点，并能有效地保护算法和数据。因此，C#语言是开发Windows数值模拟应用程序的较优选择。但是C#语言是基于Microsoft .NET Framework(简称.Net)运行环境的，在.Net环境下，一维数组和多维数组作为单个托管对象(由.Net自动分配与释放内存),其存储能力都不超过2 GB，这限制了数值模拟的计算规模(使用双精度浮点数，节点数或矩阵的阶数不得超过16 384)，因此需要研究.Net运行环境下大规模矩阵在内存中的存储方法，以提高数值模拟的计算规模。

对于线性方程组，见式(1)，不同求解方法的计算效率会相差几倍、几十倍，甚至更高，在计算大规模问题时，计算效率会严重影响模拟的进度和效率。因此需要对大规模线性方程组的快速求解方法进行研究，以提高阴极保护数值模拟的计算效率。

Ax=b

(1)

本工作针对边界元阴极保护数值模拟，使用C#语言，基于.Net运行环境和Lapack函数库，以实现数值模拟中系数矩阵并行生成、大规模矩阵存储和大规模线性方程组快速求解，对实际工程应用和开发我国自主的阴极保护数值模拟商业软件具有一定的指导作用。

1 阴极保护数值模拟的基本原理

1.1 控制方程

在电解质中，阴极保护相关电场中的电流密度和电位满足欧姆定律[11]，如式(2)所示。

(2)

式中：q为电流密度矢量;ρ为海水电阻率;φ为电位。

取电解质中一立方微元体作为研究对象，则垂直于立方微元体六个面的电流密度如图1所示。

图1 三维微元体Fig. 1 Cubic volume element

任意时刻，图1中的立方微元体中的电流变化量为

(3)

当阴极保护系统相关电场达到平衡状态且无内部场源时，由电荷守恒可知，微元体中的电流通量为0。因此通过式(3)可以得到电场域内控制方程即Laplace方程

(4)

1.2 边界条件

阴极保护边界元数值模拟的两个物理量分别是结构表面的保护电位和保护电流密度。阴极保护问题的边界条件一般有两类：第二类边界条件和第三类边界条件。

(1) 第二类边界条件，表面电流密度为0，例如：海水与空气的交界面和结构防腐蚀涂层完好部分。

(5)

式中：n为边界处的法线方向。

(2) 第三类边界条件，防腐蚀涂层破损部位，裸钢和阳极表面，电位和电流密度的关系满足极化曲线。

(6)

1.3 基于边界元的阴极保护数值模拟

对于Laplace方程，其离散形式的边界积分方程为[17]

(7)

式中：ci为节点处与边界光滑度相关的系数;N为形函数；m为单元数；φi为节点i处的电位；Γi为第j个单元；φi*为基本解；qi*为基本解的法向导数。

(8)

(9)

对每个节点，按式(7)求积分，则可得线性方程组

(10)

其中

C=diag(c1,c2,…,cn)

(11)

(12)

(13)

令

(14)

则式(11)可以化成

Hφ=Gq

(15)

带入1.2中的边界条件，将式(15)变换成线性方程组

Aφ=b

(16)

由于边界条件一般是非线性的，需要反复迭代求解方程组(16)，得到各节点处的保护电位。

对于节点数为N的模型，在求解过程中，需要首先生成H和G矩阵，生成H和G的时间复杂度为O(N2)。通常情况下，对于节点数N=10 000的模型，生成系数矩阵需要10 min左右；而对于节点数N=100 000的模型，生成系数矩阵则需要大概17 h。A是N×N的满阵，在求解线性方程组(14)时，时间复杂度为O(N3)，空间复杂度为O(N2)，使用双精度浮点数(一个双精度浮点数占8个字节的内存)存储矩阵时，占用的内存M可通过下式计算。

(17)

2 并行计算方法

本工作中单机并行计算技术的实现基于微软.NET 4.0平台，利用任务并行库TPL[18](Task Parallel Library)改写计算程序。TPL采用Task类(任务类)，它的声明对象落实到了具体的计算任务，可以自动实现线程的创建和执行。

针对阴极保护数值模拟，从式(7)可知每个节点处的积分计算量基本相同，因此可以依据节点数和线程数划分任务。本研究通过TPL中的Parallel.For函数开发了边界元系数矩阵并行生成算法。

对系数矩阵并行算法进行了测试，测试环境如表1所示，测试结果如图2所示。从图2可以看出，当线程数小于计算机处理器数目时，增加线程数可以提升计算效率。在线程数等于处理器数目时，计算效率最高，进一步增加线程数，会降低计算效率。

表1 测试环境Tab. 1 Testing enviroment

图2 计算时间与线程数的关系Fig. 2 Relationship between computing time and number of threads

在研究中发现,计算机同时进行计算的线程数不会超过处理器核数，当使用的线程数超过处理器核数时，超出的线程在前面的线程执行完毕后才会执行。由于每个线程的计算量是一样的，故在不考虑其他因素的理想情况下，本例中并行计算时间为

(18)

式中：t为计算时间;C为由总任务量决定的常数。

公式(18)说明，当NCPU等于32时，计算效率最高，继续增加NCPU会降低计算效率。

需要注意的是，公式(18)是在理想状态下给出的，不能用于求解并行计算时间，仅用来定性分析计算效率。

另一方面，使用过多的线程数会占用过多的计算机资源，申请和释放线程也会消耗一定的时间，从而也会降低计算效率。

从测试结果和以上分析可知，在并行计算的线程数等于计算机处理器核数时，并行计算的效率最高。

3 大规模矩阵存储方法

在.Net运行环境下，一维数组和多维数组作为单个托管对象(由.Net自动分配与释放内存)的存储能力都不超过2 GB，矩阵阶数不得超过16 384，这限制了数值模拟的计算规模。本工作通过使用非托管内存实现了大规模矩阵的存储。

在C#语言中，Marshal.AllocHGlobal(IntPtr cb)函数可以申请连续的非托管内存(通过代码申请与释放的内存)，IntPtr为指针类型的数据，可以表示内存中的所有地址，即AllocHGlobal函数可以申请计算机空闲的最大连续内存空间。内存申请完毕后，可以通过非安全指针进行访问。

AllocHGlobal申请的内存空间是连续的，将矩阵存储在其中后，可以直接使用Lapack函数库的矩阵运算函数进行计算，因此在.Net运行环境下，使用非托管内存进行大规模矩阵的存储是最理想的选择。

4 大规模线性方程组快速求解方法

目前，求解线性方程组主要借助Lapack函数库，Intel和AMD公司分别开发了针对自己处理器的Lapack函数库：MKL(Math Kernel Library)和ACML(AMD Core Math Library)。使用其中的矩阵运算函数，比自行编写的矩阵运算函数效率高很多。

求解线性方程组(16)，一般使用Lapack函数库中dgetrf和dgetrs两个函数，其中dgetrf函数对矩阵A进行三角分解，dgetrs利用分解后的矩阵进行求解。

目前Lapack函数库都封装在DLL(Dynamic Link Library)中，DLL分为两种：一种是直接由机器语言组成(由C、C++和Fortran等语言编译的DLL)；另一种是由.Net中间语言组成，在被调用时，再由.Net翻译为机器语言。

目前MKL和ACML提供的Lapack函数库都是采用Fortran开发的，因此其中的函数由机器语言组成。使用C、C++和Fortran语言都可以通过Lapack函数库对应的头文件直接调用其中的函数。使用C#语言可以直接调用基于.Net的动态链接库(DLL文件)中的函数，但不能够直接调用普通DLL的函数。在C#语言中，可以通过DllImport[19]调用基于机器语言的动态链接库函数。

采用表1中的测试环境，对MKL和ACML中dgetrf和dgetrs函数的求解性能进行测试，测试结果如图3所示。从图3中可以看出：在本测试环境下，ACML的求解效率是MKL的两倍左右，因此在实际使用中应根据处理器类型采用合适的函数库。

图3 矩阵求解时间与节点数的关系Fig. 3 Matrix solving time vs number of nodes

5 求解案例

采用两个实际工程案例对求解性能进行了分析，测试环境如表1所示，并行计算时使用了32个线程。

5.1 张力腿平台阴极保护数值模拟

模型包括张力腿平台的水下主体部分、钻井立管、张力筋腱和牺牲阳极，共有63 024个节点，115 336个三角形单元，如图4所示。求解方程组时，共迭代4次，模拟结果如图5所示，求解耗时如表2所示。从表2可以看出：采用并行生成系数矩阵和AMD Lapack进行求解，可以大幅提升计算效率。

图4 张力腿平台边界元网格模型Fig. 4 Boundary element model of tension-leg platform

图5 张力腿平台牺牲阳极阴极保护的模拟结果Fig. 5 Simulated result of sacrificial anode cathodic protection for tension-leg platform

计算方案生成系数矩阵时长方程组求解时长总时长串行生成&MKL求解566.2113.6679.8并行生成&AMD求解18.057.475.4

5.2 半潜平台阴极保护数值模拟

模型包括半潜平台的水下主体部分、钻井立管和张力筋腱，共有56 594个节点，111 204个三角形单元，如图6所示。求解方程组时，共迭代3次，模拟结果如图7所示，求解耗时如表3所示。从表3可以看出：采用并行生成系数矩阵和AMD Lapack进行求解，可以大幅提升计算效率。

图6 半潜平台边界元网格模型Fig. 6 Boundary element model of semi-submersible platform

图7 半潜平台外加电流阴极保护的模拟结果Fig. 7 Simulated result of impressed current cathodic protection for semi-submersible platform

计算方案生成系数矩阵时长方程组求解时长总时长串行生成&MKL求解456.961.7518.6并行生成&AMD求解14.531.345.8

6 结论

本工作使用C#编程语言，在.Net运行环境下，针对阴极保护边界元数值模拟，从系数矩阵生成、存储和大规模线性方程组求解三个方面对高性能计算方法进行了研究，并得到了以下结论。

(1) 使用单机并行计算方法可以显著提升边界元数值模拟中系数矩阵的生成效率，线程数等于处理器数目时，计算效率最高；

(2) 在.Net环境下，使用非托管内存可以使用计算机空闲的最大连续内存空间，是数值模拟中存储大规模矩阵的理想选择；

(3) 在.Net环境下，可以使用Lapack函数库中的高性能运算函数对线性方程组进行求解，不同处理器应采用相应的函数库。