实用COMSOL后处理二次开发技术＊

时间：2024-05-18

周子杰，刘英伟，张洋，张建能

实用COMSOL后处理二次开发技术＊

周子杰，刘英伟，张洋，张建能

（哈尔滨工程大学材料科学与化学工程学院，黑龙江哈尔滨 150001）

采用COMSOL软件对船体的腐蚀与防护进行了模拟，模拟结果需要通过后处理模块显示出来。一般情况下，COMSOL的后处理模块所包含的功能能满足大部分要求，但有些时候，对于一些特殊的要求却难以满足。比如计算某一等值线所包围的面积以及沿某后天截取的曲线上显示某物理量的分布等。针对这些不足，基于空间几何和线性插值理论，提出一套算法，并根据该算法编制了程序进行计算。结果表明，面积计算非常精确，曲线物理量的计算结果和精确解十分吻合，说明提出的算法正确、程序可靠。该算法也给软件使用者提供了更多的表示形式，同时也大大地拓展了该软件的功能。该算法具有通用性，可以移植到其他软件上。

腐蚀；阴极保护；有限元；插值计算

COMSOL软件全称为COMSOL Multiphysics，是瑞典COMSOL公司于2003年正式推出的全球第一款真正的多物理场耦合分析软件，它以其方便、易用、高效等优点而被广泛应用于各个领域的科学研究和工程计算当中[1-4]。其中的腐蚀模块专门用来模拟和仿真水或电解质溶液中，金属结构的腐蚀及其防护，很多学者使用该模块进行了大量卓有成效的研究[5-7]。该软件功能强大，但也存在不足，比如对于一个防腐设计者而言，除了关心被腐蚀构件的表面电位、电流分布外，还关心构件表面腐蚀或未腐蚀面积的大小，即某一等值线所包围的面积，而COMSOL的后处理模块却没有提供这样的功能；此外，COMSOL的后处理无法将某一曲线上的物理量以一维图形显示出来。对于那些几何建模时，自然形成的曲线（先天曲线，如图1（b）中的I）后处理可以实现上述功能，但对于后天曲线（如图1（b）中的II、III，它是计算完毕后，用某个截面与船体相截形成的曲线），则现有的后处理模块不能实现上述功能。除非在建模开始，就对感兴趣的部分进行切割以形成先天曲线，但这无疑大大增加了建模的复杂性，而且对哪些部位感兴趣一般也无法确定。针对这些不足，本文开发了一套程序，该程序能将计算结果做进一步处理，通过一定的算法，能实现上述功能。这为设计者提供了更多的、更直接信息，这大大有利于设计的优化，同时也大大扩展了该软件的功能。

1 控制方程

带有四个辅助阳极的船体几何模型如图1所示，阳极材料为惰性有色金属铂。当船体浸没在海水中时，除了铂表面外，其他表面将遭受海水的腐蚀，发生如下电化学反应：

2Fe→2Fe2++4e. （1）

O2+2H2O+4e→4OH-. （2）

在反应（2）中，水中的氧不断地从船体表面夺取电子而发生还原反应；而反应（1）需要不断地为反应（2）提供电子；失去了电子的铁原子变成离子状态，在水分子这种极性分子的作用下，很容易溶解到水中，从而导致了船体的腐蚀。为了阻止腐蚀的发生，必须抑制反应（1）的进行。外加电流阴极保护，就是采用外电源向船体表面注入电流，以代替反应（1）。注入的电流使船体表面产生阴极极化，当极化电位达到800 mV时（参比电极为Ag/AgCl/Seawater），表面就会得到保护。

当外电源向船体注入电流时，由于海水是强电解质，具有导电性，因此会在船体周围的一块水体里产生电场，存在电位分布，如果围绕船体周围划分出一块适当大小的水体Ω，如图2所示，则当水体中电位处于稳态时，其分布满足Laplace方程：

式（3）中：为电解液电势；，，∈Ω。

1.1 边界条件

图2中的水体Ω具有封闭的表面，其上存在以下几种边界条件。

1.1.1 电化学反应边界条件

图2中：边界①上发生电化学反应（1）和（2），其电流和电势关系遵循Tafel公式：

表1 电化学参数

参数铁氧平衡电位/V0.760.189 塔菲尔斜率/（V/m）0.410.18 交换电流度/（A/m2）7.7e-77.1e-5 修正后交换电流密度/（A/m2）1.078 7e-79.94e-6 海水电导率/S4

式（6）中：，分别是和涂层厚度以及船体所处环境有关的参数，分别为0.02和0.012；是使用年限，这里按10年计算，这样计算得到破损率。

确定了破损率以后，按下式修正交换电流密度：

式（7）中：c是修正后的交换电流密度；是破损率为0时的交换电流密度，具体数据如表1所示。

1.1.2 电流边界条件

电极和海水的接触面②共有4个（图2（a）只标识了2个）。当外电源向船体表面注入电流时，电流也会通过接触面注入到海水中，因此存在如下关系：

式（8）中：为海水电导率；为单位表面法向矢量；_为注入的电流密度。

外加电流密度是通过 DNV RP B401规范确定的：在船体裸露的情况下，要使船体表面达到保护电位，所需的电流密度为120 mA/m2。考虑到破损率，所需电流密度变为=c×120 mA/m2=16.8 mA/m2。船体总面积为11.564 m2，因此所需要的总电流为：=11.564×16.8=194.3 mA=0.194 3 A。因为有4个电极供电，所以每个电极的供电电流为0.194 3/4=0.048 568 8≈0.05 A。

1.1.3 无限远边界条件

所谓无限远不是空间位置的无限远，而是指电场强度很弱，近乎为0的地方，如图2（a）中的表面③，此处表面电流可看作为0，公式为：

图2 边界条件与有限元离散

图3 计算结果

1.2 几何模型离散并求解

设置完边界条件后，对几何模型进行有限元离散，结果如图2（b）所示，之后就可以进行计算。计算结果如图3所示。由图3可知，靠近电极的地方电位较高，电极中间区域电位低一些。总体来看，船体表面电位均小于保护电位0.8 V，可见整个船体都没有得到保护，先前估算的电流偏低，应适度增加电流供应。

2 后处理的二次开发

对设计者而言，有时候还需要了解表面腐蚀面积的大小，即某一等值线所包围的面积（以下简称“面积”），然而COMSOL后处理模块却没有提供这一功能，而这一数据对设计者而言无疑是关键的。另外，计算完毕后，设计者常常想把某物理量沿感兴趣的曲线显示出来，对于那些建模时自然形成的曲线，COMSOL的后处理模块能实现这一功能，但是对于后天通过平面和几何体相交得到的曲线，后处理模块却无法实现这一功能，这大大限制了该软件的功能。为此本文开发了一套程序，通过对计算结果进行二次处理，实现了上述功能。

2.1 数据的导出

由于涉及的计算都和船体表面有关，因此需要把船体表面信息输出到文本文件中，再导入程序中处理。船体表面原来是曲面，用四面体单元离散几何体后，表面就由大量平面三角形组成，因此需将这些三角形单元的信息输出。这些信息包括单元编号及其所包含的节点号、节点的坐标、每个节点的电位等。由于计算域有很多表面，而所处理的表面只是其中一部分，如图3（a）所示，因此需要将这部分表面单独定义。这可以通过COMSOL 软件中的Definitions功能实现。接下来要把这部分单独定义的表面信息提取出来（COMSOL的计算结果遍布整个求解域，并不区分属于哪个表面），方法是在后处理模块的Results节点中，选择Data Sets命令，把先前定义的表面的有关数据从总的结果中提取出来；最后，在后处理模块Results中，选择Export→Data命令，此时右侧会出现界面，在界面里，可以选择要处理的数据集(先前定义的Data Sets会出现)；在界面的Output栏目下，选择数据格式（这里选为Sectionwise）、存储数据的路径和文件名等。需要注意的是，在该界面的Advanced栏目下，不要选Include header一项，因为这会在数据文件中混入英文字符，导致数据处理麻烦。当选择、设置完毕后，按界面上方的Export按钮，就可以把数据输出到文本文件中了，将文本文件打开，可以看到表2所示的数据格式（片段）：前1 828行记录的是1 828个节点的坐标，其编号就是所在的行号，第1 829～5 323行中，每一行代表一个单元所属的三个顶点的编号（例如第1 829行为第1行），据此可知单元数目为5 323－1 829+1=3 495.从5 324行开始到最后，每一行都代表一个节点的电位，其中5 324行和第1行相对应，其余类推到结束。

2.2 三角形单元被保护面积的计算

将上述数据导入到程序中，并设计合适的变量将其储存起来，然后就可以利用这些数据计算保护面积的大小。由于几何体采用四面体单元离散，单元形状函数为线性[8]，因此四面体的每一个三角形的单元的形状函数也为线性，因此单元边上的电位分布为线性（如图4（b）所示）。另外需要注意的是，三角形顶点坐标是三维的。

表2 数据结构

1﹣5.320.867 461 8310.097 393 2﹣5.320.963 158 033﹣0.069 61 3﹣5.320.818 393 3240.089 777 ………… 1 828﹣0.412 73﹣0.014 350.3 1 829123 1 830142 1 831564 ………… 5 3232342531 827 5 3240.801 651 5 3250.802 263 ………… 7 1520.812 09

2.2.1 特殊情形的处理

如果设保护电位PROTECT=0.8 V的话，那么三角形3个顶点电位和保护电位相比存在几种大小关系，如表3所示。情况1下，由于3个顶点电位都大于0.8 V，因此，根据插值函数的性质，整个三角形内任意点的电位均大于0.8 V，因此整个单元均得到保护，保护面积就等于三角形面积，其大小可根据相关公式[9]求得；情况2的情形和情况1相反，3个顶点的电位均小于0.8 V，同样根据插值函数的性质，三角形内任一点的电位均低于0.8 V，因此整个三角形都没有被保护，保护面积为零。情况3比情况4复杂些，需进一步讨论。

图4 三角形单元

2.2.2 一般情况的处理

以情况3为例，图4中三角形123的边12，可用三维空间直线方程描述：

式（10）中：为参数且∈[0，1]。

因为在12边上，电位分布为线性，因此梯度为：

这样，边12上任意一点的电位可通过线性插值确定：

式（12）中：1a为点1、a间距离，待求。

由于1＞0.8且2＜0.8，因此在1、2点之间，一定存在一点a，使得a=PROTECT=0.8 V；同理边23上同样存在，使b=PROTECT=0.8 V。这样一来，根据插值函数的特点，三角形Δ2内任意点的电位均低于0.8 V，因此三角形Δ2未得到保护。

将a=0.8代入式（12），就可以求得1a，再根据式（10），得到1a的另一表达式：

二者相等，就确定了a，将a代入式（10），就求出了点的坐标。同理，可求出23边上b点的坐标。这样三角形2三个顶点的坐标已知，三角形的面积Δa2b也就确定了，这样一来，三角形123中只有部分面积得到了保护，其大小为PROTECTED=Δ123－Δa2b.情况4的处理与此类似，只不过受保护面积为Δa2b而不是Δ123－Δa2b.这样，一个三角形单元内，受保护部分的面积就确定下来了，将所有三角形做如此处理，最后将总的保护面积累加起来就得到整个表面被保护面积的大小。

2.3 三角形截线上的物理量的计算

这里仅考虑垂直于轴、轴的特殊截面和船体表面交截后形成的曲线，如图1（b）中的曲线II、III。曲线是由截面和船体相交形成的，由于船体离散后，船体表面实际上是由平面三角形构成，因此截面和船体表面相交，实际上就是与大量的三角形相交，因而会截得大量直线段，这条曲线，实际上是这些直线段构成的折线。求出线段端点（交点）的空间坐标和物理量（电位）后，那么沿着这条折线物理量的分布就得到了。下面以曲线III（是由截面=0=0截得）为例，介绍如何求出上述信息。根据三角形单元和截面的位置关系，将二者相交情况分为5类，如表4所示。

2.3.1 特殊情况

在各种相交关系中，情况1因为单元与截面无交点，因此不予考虑；情况2、3中，单元顶点和截面重合，因此不用计算，直接将单元顶点信息存储到数据结构Point_Cut中，以情况2为例：

count_point=count_point+1

Point_Cut（count_point，1）=x2

Point_Cut（count_point，2）=y2=y0

Point_Cut（count_point，3）=z2

Point_Cut（count_point，4）=V2

其中，Point_Cut为数据结构，存储了三角形与截面相交的点的信息；count_point为结构当前存储的交点数目。情况4、5中的三角形的边都和截面相交，须通过一定的算法求得交点信息。截面与三角形位置关系分类如表4所示。

表3 三角形顶点电位大小比较关系

情况分类保护面积情况1V1＞0.8，V2＞0.8，V3＞0.8SPROTECTED=SΔ123 情况2V1＜0.8，V2＜0.8，V3＜0.8SPROTECTED=0 情况3V1＞0.8，V2＜0.8，V3＞0.8SPROTECTED=SΔ123-SΔa2b 情况4V1＜0.8，V2＞0.8，V3＜0.8SPROTECTED=SΔa2b

2.3.2 一般情况

以情况5中的三角形①为例，截面（=0）与三角形123相交，在边12和13上分别截取交点、，如图5所示。

三角形123的边12，用形如公式（10）的三维直线方程描述：

由于n=0，代入式（13）求得：

将1代入式（13），就可确定n、n。同理可确定点坐标（m，m，m）。接下来求、点的物理量——电位，这可以根据式（12）确定：

其中：

1、2分别为线段12和13上电位梯度：

表4 截面与三角形位置关系分类

三角形单元与截面的位置关系判据情况1Cut plane（y=y0）①（y1－y0）＜0∩（y2－y0）＜0∩（y3－y0）＜0；②（y1－y0）＞0∩（y2－y0）＞0∩（y3－y0）＞0 情况2Cut plane（y=y0）①（y1－y0）＜0∩（y2－y0）=0∩（y3－y0）＜0 ；②（y1－y0）＞0∩（y2－y0）=0∩（y3－y0）＞0 情况3Cut plane（y=y0）①（y1－y0）=0∩（y2－y0）=0∩（y3－y0）＜0；②（y1－y0）=0∩（y2－y0）=0∩（y3－y0）＞0 情况4Cut plane（y=y0）（y1－y0）=0∩（y2－y0）＞0∩（y3－y0）＜0 情况5Cut plane（y=y0）①（y1－y0）>0∩（y2－y0）<0∩（y3－y0）<0；②（y1－y0）>0∩（y2－y0）>0∩（y3－y0）<0

至此，三角形与某一截面相交，交点的坐标以及电位均可求得，将这些结果加入到结构Point_Cut中。每个三角形都这样处理，最终就得到了一个由截面和船体表面相截而得到的曲线上若干交点的信息，经过进一步处理（去掉冗余、排序等），就可以把这些物理量沿曲线的分布绘制出来，如图6（d）、6（f）所示。

2.3.3 曲线II的求解

曲线II和曲线III的不同在于二者所垂直的坐标轴不同，但在算法上，十分类似，只需对2.3.2中的算法略作修改即可。此时曲线II是由垂直于轴的平面所截，设此平面方程为=0，将它代入式（13）中第一式，就可确定参数，进而可确定节点坐标，接下来，参照2.3.2中的余下步骤就可确定节点电位，这里不再赘述。

2.4 算例

由图3（a）可知，船体表面未得到有效保护，必须增大电流密度，现将电流密度增加到0.1 A/m2，结果如图6（a）、6（b）所示，此时，表面电位有所增大。用开发的程序计分别算了等值线为0.795 2 V、0.8 V（保护电位）和0 V所包围的面积，结果如表5所示。可见，尽管0.795 2 V和0.8 V相差并不大，可所包围的面积却相差很大，不过由于0.795 2 V和0.8 V相差不大，因此可以认为船体表面大部分得到了保护。另外，等值线0所包围的面积就是船体底面总面积，计算结果为12.193 9 m2。

图5 交点信息的求解

表5 不同电位等值线包围面积

电流密度/（A/m2）0.1 等值线值/V0.795 20.80 包围面积/m27.727 12.359 312.193 9

COMSOL后处理具有计算某一表面的总面积的功能，是在Results→Derived Values菜单里实现的。为了验证算法的正确性，本文利用COMSOL提供的这一功能计算了船体表面积，如图3（a）所示，结果为12.194 m2，而本程序计的结果为12.193 9 m2，二者基本一致（考虑舍入的话，二者几乎相等），这证明了算法正确性以及程序的可靠性。

图6（c）、6（e）是在建模时，在相关位置进行切割得到曲线(先天曲线)，经程序计算后，得到的物理量沿曲线精确分布；6（d）、6（f）分别为没有经过切割，直接采用本文开发的程序，在相同部位曲线上计算得到的物理量的分布。由图可以看出，二者符合得非常好，说明算法的正确性和程序的可靠性。