阴极保护是金属腐蚀防护中常用的电化学保护方法。在阴极保护系统中，系统的电位分布及电流密度是衡量保护效果的重要参数，目前广泛采用数值模拟的方法进行分析计算。常用的数值方法主要用有限差分法[1, 2, 3]、有限单元法[4, 5, 6, 7]及边界元法[8, 9, 10, 11]。其中边界元法具有精确模拟实际边界复杂构形和只在模型边界离散网格的特点，与有限差分法和有限元法等基于体积网格划分的数值方法相比，其在减少方程未知量的同时，可以有效降低网格划分的难度。
然而常规边界元法并不适合处理大规模问题，原因在于边界元法形成的系数矩阵通常是满阵，存储量是O(N2)(其中N是未知量数量)； 使用常规直接或迭代求解器需要O(N2)~O(N3)的计算量。当N增加时，常规边界元法的存储量和计算量将快速上升。为了提高边界元法的求解效率，近年来一些快速数值算法在该领域得到了应用。其中特别值得关注的是快速多极算法(FMM)[12, 13, 14, 15]，该算法可将边界元法的存储量和计算量均降到O(N)，从而可以有效处理大规模问题。快速多极边界元法已经广泛应用于弹性力学和传热分析等领域的大规模计算[16, 17, 18, 19, 20, 21, 22]，关于该算法的综述可见文[23, 24]。 张东东等[25]将快速多极边界元法应用于二维阴极保护系统分析并验证了其有效性。
阴极保护系统中通常包含大量圆柱状管道结构，如输油管道、埋地钢桩等。该结构的主要特征是在一个方向上的尺度远大于其余2个正交方向。用常规三维面单元离散将会产生数量庞大的边界单元，同时单元形状易畸化，导致计算效率和精度均可能下降。考虑到管道结构单一尺度占优的特点，可以预计其任意横截面上电位的变化与长度方向电位的变化相比可以忽略不计。该特征可用于简化管道结构的边界元模型。
本文提出适用于阴极保护系统中圆柱状管道结构的三维管道边界单元模型。该模型将管道结构离散为线单元，同时保留在对应圆柱面上积分的特征，把二维积分转化为一维或解析积分。将快速多极算法用于边界元法的加速求解，以便能够在普通个人电脑上处理大规模的阴极保护系统分析。针对阴极保护系统中的阴极极化边界条件，采用Newton-Raphson法迭代求解由此引入的非线性问题。
1 三维管道边界元模型 1.1 阴极保护分析的边界积分方程阴极保护分析对应的静电场控制方程和边界条件为：
| $ \left\{ \begin{align} & {{\nabla }^{2}}u\left( x \right)=0,\ \ \ \ \ \ \ \ \ \ \ \ x\in V; \\ & u\left( x \right)={{u}_{o}},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\in {{S}_{1}}; \\ & q\left( x \right)=-k\frac{\partial u}{\partial n}={{q}_{0}},\ \ \ x\in {{S}_{2}}; \\ & q\left( x \right)=-k\frac{\partial u}{\partial n}=fc\left( u \right),\ \ \ x\in {{S}_{3}}.\ \ \\ \end{align} \right. $ | (1) |
其中：u是电位；q是电流密度；k是介质电导率；S1、S2和S3分别对应给定电位、电流密度和极化边界。当求解域V是均匀介质时，式(1)对应的边界积分方程为：
| $ \begin{align} & c\left( x \right)u\left( x \right)+\int_{S}{\frac{\partial {{G}^{*}}\left( x,y \right)}{\partial n}}u\left( y \right)dS\left( y \right)= \\ & \int_{S}{{{G}^{*}}}\left( x,y \right)\left[-\frac{1}{k}q\left( y \right) \right]dS\left( y \right). \\ \end{align} $ | (2) |
其中：x和y分别表示边界上的源点和场点；n是场点处的单位外法向量； c(x)是与源点所在边界形状有关的常数项； G*(x,y)是Laplace方程的基本解，即
| $ {{G}^{*}}\left( x,y \right)=\frac{1}{4\pi r}. $ | (3) |
其中r是场点到源点的距离。用常值单元将边界离散为Ne个边界单元，通过配点格式得到式(2)的离散形式：
| $ \begin{align} & \frac{1}{2}{{u}_{i}}+\sum\limits_{j=1}^{{{N}_{e}}}{{{u}_{j}}\int_{{{S}_{j}}}{\frac{\partial {{G}^{*}}\left( {{x}_{i}},y \right)}{\partial n}}}dS\left( y \right)= \\ & -\frac{1}{k}\sum\limits_{j=1}^{{{N}_{e}}}{{{q}_{j}}\int_{{{S}_{j}}}{{{G}^{*}}\left( {{x}_{i}},y \right)}}dS\left( y \right),\ \ \ i=1,...,{{N}_{e}}. \\ \end{align} $ | (4) |
其中ui和qi分别是定义在第i个单元面Si上的电位和电流密度。式(4)的矩阵形式为：
| $ \left( C+H \right)\centerdot u=-\frac{1}{k}G\centerdot q. $ | (5) |
其中：
| $ C=diag\left( \frac{1}{2},\frac{1}{2},\ldots ,\frac{1}{2} \right); $ | (6) |
| $ {{H}_{ij}}=\int_{{{S}_{j}}}{\frac{\partial {{G}^{*}}\left( {{x}_{i}},y \right)}{\partial n}}dS\left( y \right), $ | (7) |
| $ \begin{align} & {{G}_{ij}}=\int_{{{S}_{j}}}{\frac{\partial {{G}^{*}}\left( {{x}_{i}},y \right)}{\partial n}}dS\left( y \right); \\ & \ \ \ u={{\left( {{u}_{1}},{{u}_{2}},\ldots ,{{u}_{{{N}_{e}}}} \right)}^{T}},\\ & \ \ \ q={{\left( {{q}_{1}},{{q}_{2}},\ldots ,{{q}_{{{N}_{e}}}} \right)}^{T}}. \\ \end{align} $ | (8) |
通过引入边界条件、重新排列已知项与未知项，可以得到方程组A·x=B，其中A是系数矩阵，该系数矩阵通常为满阵。
1.2 三维管道模型及线单元对于阴极保护系统中常见的细长管道结构，其直径远小于长度，在其任意横截面内的电位变化也远小于沿轴线方向的变化，因此可假设在任意横截面内电位为常数，电位仅沿轴线方向变化。基于此假设，可将管道结构简化为三维管道模型，在该模型中： 使用轴线代替管道； 划分网格时将轴线离散成线单元，每个线单元代表一段管道；单元积分在线单元对应的管道圆柱面上进行。
使用线单元对管道进行离散后，系数矩阵的计算可在建立在管道轴线上的局部柱坐标系下进行。对于第j个常值线单元Sj，以其对应的管道轴线的中点为坐标原点，轴向为Z向建立柱坐标系Cylinder-1。 单元的法向指向管道轴线，单元长度为Lj，半径为Rj，场点y的柱坐标为y(ρ,θ,z)。坐标系Cylinder-1如图1所示。
|
| 图 1 线单元Sj的柱坐标系Cylinder-1 |
根据源点xi的位置不同，对积分的计算可分成以下3种情况：
1) 当xi在线单元上时，将其定义在轴线中点位置，
此时的积分类似于面单元上的奇异积分，但由于xi不在圆柱面上，而属于求解域外，因此积分本质上是非奇异的。此时应取c(xi)=0，而Gii和Hii可解析积出：
| $ \begin{align} & {{H}_{ij}}=\frac{{{L}_{i}}}{\sqrt{L_{i}^{2}+4R_{i}^{2}}},\\ & {{G}_{ij}}=\frac{1}{2}{{R}_{i}}\ln \left( \frac{\sqrt{L_{i}^{2}+4R_{i}^{2}}+{{L}_{i}}}{\sqrt{L_{i}^{2}+4R_{i}^{2}}-{{L}_{i}}} \right). \\ \end{align} $ | (9) |
2) 当xi不在线单元上、且不在该单元轴线延长线上时，Gij和Hij沿管道轴向(Z向)的积分可解析积出，沿周向(θ方向)的积分需要用数值积分的方法计算，即
| $ \begin{align} & {{H}_{ij}}=\frac{{{R}_{j}}}{4\pi }\int_{0}^{2\pi }{\frac{{{R}_{j}}-{{r}_{s}}\cos \left( \theta -{{\theta }_{s}} \right)}{R_{j}^{2}+r_{s}^{2}-2{{R}_{j}}{{r}_{s}}\cos \left( \theta -{{\theta }_{s}} \right)}}. \\ & \left[\frac{\frac{1}{2}{{L}_{j}}-{{z}_{s}}}{{{R}_{R}}\left( \theta \right)}+\frac{\frac{1}{2}{{L}_{j}}+{{z}_{s}}}{{{R}_{L}}\left( \theta \right)} \right]d\theta ,\\ & {{G}_{ij}}=\frac{{{R}_{j}}}{4\pi }\int_{0}^{2\pi }{\ln }\frac{-{{z}_{s}}+\frac{1}{2}{{L}_{j}}+{{R}_{R}}\left( \theta \right)}{-{{z}_{s}}-\frac{1}{2}{{L}_{j}}+{{R}_{L}}\left( \theta \right)}d\theta . \\ \end{align} $ | (10) |
其中：(Rj,θ,z)、(rs,θs,zs)分别是场点y和源点xi在柱坐标系Cylinder-1下的坐标值。RL(θ)和 RR(θ)分别是图1中xi到2个端面圆周上A点和B点的距离(如图1所示)，与θ有关。从式(10)可以看出，线单元积分为一维，因此与三维面单元相比，积分计算量进一步降低。
3)当xi不在线单元上、但在该单元轴线延长线上时，Gij和Hij可解析积出：
| $ \begin{array}{l} {H_{ij}} = \frac{1}{2}\left[{\frac{{{L_j}/2 - {z_s}}}{{\sqrt {R_R^2{\rm{ + }}R_j^2} }} + \frac{{{L_j}/2 + {z_s}}}{{\sqrt {R_L^2{\rm{ + }}R_j^2} }}} \right]\\ {G_{ij}} = \frac{1}{2}{R_j}\ln \left( {\frac{{ - {z_s} + {L_j}/2 + \sqrt {R_R^2{\rm{ + }}R_j^2} }}{{ - {z_s} - {L_j}/2 + \sqrt {R_L^2{\rm{ + }}R_j^2} }}} \right). \end{array} $ | (11) |
其中：RL和RR分别是xi到单元左右端面圆周上任意一点的距离，与θ无关，有
| $\begin{array}{l} {R_L} = \sqrt {R_j^2 + \left( {{z_s} + {L_j}/2} \right)} ,\\ {R_R} = \sqrt {R_j^2 + \left( {{z_s} - {L_j}/2} \right)} . \end{array} $ | (12) |
需要注意的是：对于管道结构终点的线单元，应将其对应终点的圆端面积分考虑在内。该积分可用二维数值积分求解，由于管道终点数量与总体的线单元数量相比可忽略，该端面积分在整体积分中所占的计算量比例也可忽略。
1.3 非线性极化曲线边界条件的处理采用迭代方法处理极化曲线引入的非线性边界条件。将极化曲线分段线性化，逐段计算斜率a和截距b，从而将极化曲线上u和q的关系表示为：
| $ u = a\left( q \right) \cdot q + b\left( q \right). $ | (13) |
在集成系数矩阵A和右手向量B时，对于第j个具有极化边界的单元Sj，A和B对应元素表示为：
| $ \begin{array}{l} {A_{ij}} = a\left( q \right) \cdot \left( {{c_{ij}} + {H_{ij}}} \right) + \frac{1}{k}{G_{ij}},\\ {B_i} = - b\left( q \right) \cdot \left( {{c_{ij}} + {H_{ij}}} \right). \end{array} $ | (14) |
可以看出此时方程组A·x=B为非线性，迭代求解过程如图2所示。
|
| 图 2 非线性边界条件迭代求解流程图 |
快速多极边界元法的基本思想是：采用自适应树结构作为主要存储和操作对象，在迭代求解方程组的每一个迭代步中，通过树结构操作等效完成系数矩阵与迭代向量的相乘运算，从而避免了系数矩阵的显式存储。由于树结构的存储和操作运算为 O(N)量级，快速多极边界元法是O(N)量级的数值算法，从而能够在普通微机上有效处理大规模问题。
快速多极边界元法主要包含2个关键步骤：多极展开和局部展开。
2.1 多极展开和多极展开系数的传递选择参考点Oy，满足$\left| {\overrightarrow {{O_y}y} } \right| < \left| {\overrightarrow {{O_x}y} } \right|$ 则式(2)中的基本解可展开成2个函数系列的乘积之和，其中一个函数系列只与源点有关，另一个与场点有关，从而分离源点和场点的直接联系。该展开也称作多极展开，本文参考文[24]给出基本解的一个紧凑多极展开格式：
| $ \begin{array}{l} {G^ * }\left( {x,y} \right) = \\ \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {\overline {{S_{nm}}} \left( {\overrightarrow {{O_y}x} } \right){R_{nm}}\left( {\overrightarrow {{O_y}y} } \right)} } . \end{array} $ | (15) |
其中：Snm和Rnm是球谐函数[24]； n是多极展开阶数。将展开代入积分可得：
| $ {G_{ij}} = \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {\overline {{S_{nm}}} \left( {\overrightarrow {{O_y}x} } \right){c_{nm}}\left( {{O_y}} \right).} } $ | (16) |
其中${c_{nm}}\left( {{O_y}} \right) = \int_{{S_j}} {{R_{nm}}} \left( {\overrightarrow {{O_y}y} } \right)q\left( x \right)dS\left( y \right)$称作多极展开系数。可以看出，cnm(Oy)与源点无关，因此一次计算得到cnm(Oy)即可用于多个源点的Gij计算，从而有效减少了计算量。利用Gij与Hij的关系，可以推导得到类似的Hij的多极展开格式和多极展开系数。
选择新的参考点Oy'，满足$\left| {\overrightarrow {{O_y}'y} } \right| < \left| {\overrightarrow {{O_x}'y} } \right|$，则以Oy'为展开点可得到新的多极展开系数cnm(Oy')。该系数可通过cnm(Oy)线性映射得到，即多极展开系数的传递[24]为：
| $ {{c}_{nm}}\left( {{O}_{y}}' \right)=\sum\limits_{n'=0}^{\infty }{\sum\limits_{n'=-n'}^{n'}{{{R}_{n'm'}}\left( \overrightarrow{{{O}_{y}}'{{O}_{y}}} \right){{c}_{n-m',m-m'}}\left( {{O}_{y}} \right).}} $ | (17) |
选择参考点Ox，满足 $\left| \overrightarrow{{{O}_{x}}x} \right|<\left| \overrightarrow{{{O}_{x}}y} \right|$，基本解同样可展开成2个函数系列的乘积之和，称作局部展开。基本解对应积分的局部展开格式为[24]：
| $ {{G}_{ij}}=\frac{1}{4\pi }\sum\limits_{n=0}^{\infty }{\sum\limits_{m=-n}^{n}{{{R}_{nm}}\left( \overrightarrow{{{O}_{x}}x} \right){{d}_{nm}}\left( {{O}_{x}} \right).}} $ | (18) |
其中dnm(Ox)称作局部展开系数，可由多极展开系数线性映射得到，即多极展开向局部展开系数的传递[24]为：
| $ {{d}_{nm}}\left( {{O}_{x}} \right)=\sum\limits_{n'=0}^{\infty }{\sum\limits_{n'=-n'}^{n'}{{{\left( -1 \right)}^{n}}\overline{{{S}_{n+n',n+m'}}}\left( \overrightarrow{{{O}_{y}}{{O}_{x}}} \right){{c}_{n'm'}}\left( {{O}_{y}} \right).}} $ | (19) |
选择新的参考点Ox'，满足$\left| \overrightarrow{{{O}_{x}}'x} \right| < \left| \overrightarrow{{{O}_{x}}'y} \right|$以O'x为展开点的新的局部展开系数dnm(O'x)可通过 dnm(Ox)线性映射得到，即局部展开系数的传递[24]为：
| $ {{d}_{nm}}\left( {{O}_{x}}' \right)=\sum\limits_{n'=n}^{\infty }{\sum\limits_{m'=-n'}^{n'}{{{R}_{n'-n,m'-m}}\left( \overrightarrow{{{O}_{x}}{{O}_{x}}} \right){{d}_{n'm'}}\left( {{O}_{x}} \right).}} $ | (20) |
数值算例使用的边界元计算程序使用C++编写，运行在Intel Core 2 Duo 2 GHz的笔记本电脑上，内存为2 GB。 算例中常规表面边界元划分采用三角形常值单元。该单元格式简洁，易于网格划分并满足工程精度，是边界元分析中常用单元类型之一。
3.1 线单元计算精度验证本算例研究线单元的计算精度。如图3所示，长方体尺寸为10 m×10m×30m，电导率为k=5Ω－1·m－1，中心处布置2根管道，长度为12 m、半径为0.2 m； 管道轴线处在一条直线，间距为 2 m，管道两端到长方体表面的距离也是2 m。 长方体两端(阴影部分)分别给定电位150 V和极化边界条件(极化曲线如图4所示)，管道表面分别给定电流100 A/m2及电位10 V。
|
| 图 3 包含两条管道的长方体模型示意 |
|
| 图 4 极化曲线 |
将长方体表面离散成11 200个三角形常值单元。管道分别用三角形单元和线单元离散，其中采用三角形单元离散时，共有832个单元； 采用线单元离散时，共有40个单元。使用预处理广义极小残量法(GMRES)迭代求解器进行求解，管道1的电位和管道2的电流计算结果见图5，另外采用商用有限元软件Abaqus对相同的模型进行了计算。有限元模型中根据对称性取四分之一建模，单元数量为89 832个，远多于边界单元数量。
|
| 图 5 算例1管道电流和电位分布 |
由图5可见，三角形常值单元和线单元计算结果与有限元结果总体一致，只是在端部线单元精度略有降低。线单元与三角形单元的平均相对误差为1.4%。由于每根管道的线单元离散数量远小于三角形单元离散数量，该算例表明线单元在有效减少离散难度和单元数量的同时，达到较高的计算精度。
3.2 大规模线单元的计算精度验证本算例研究包含大规模管系的模型计算精度。模型是一个24 m×24 m×50 m的长方体区域，包含5×5×4共100根直径0.4 m、 长10 m的管道，分5层排布，每层如图6所示。电导率k=5 Ω－1·m－1，2个正方形表面分别给定电位10 V和150 V，其余表面为绝缘表面。若用有限元离散该模型，将产生普通个人计算机难以处理的数量庞大的体积网格。用边界元离散，当用三角形边界单元离散管道时，共产生41 600个单元，模型单元总数为53 504，使用快速多极算法求解； 而使用线单元时，共产生2 000个线单元，模型单元总数仅为13 904，使用预处理GMRES迭代求解器求解。图7给出了分别用三角形单元和线单元计算得到的模型中心处4根管道的电位分布。从图7中可以看出，2种离散方法得到的计算结果一致。该算例表明，当模型中含有大量细长管道时，相对于采用三角形单元离散，采用线单元离散可以节省大量计算资源，同时大规模的线单元仍然能够保持较高的计算精度。
|
| 图 6 包含5×5×4=100条管道的长方体模型平面示意图 |
|
| 图 7 中心位置管道的电位分布 |
本算例研究基于三维管道模型的快速多极边界元法在阴极保护系统分析中的应用效果，并与常规边界元法进行对比。图8中模型的基本参数为：罐底直径为40 m，给定极化曲线边界条件(见图9)； 将罐底四周直径为90 m、 深度为60 m的圆柱体土壤作为求解域，土壤电阻率为0.10 Ω－1·m－1，界面设为绝缘表面； 在罐底正下方分布有3根长为 70 m、 平行于y轴的管道，深度分别为10 m、30 m和50 m，管道直径为1 m，表面设为绝缘表面； 阳极设置在土壤中段，长20 m，宽7 m，给定电位7 V。
|
| 图 8 阴极保护系统分析模型示意 |
|
| 图 9 罐底极化曲线 |
采用2种网格密度进行划分：使用常规GMRES迭代算法求解时的网格单元总数为8 434，迭代收敛误差为1×10-5； 使用快速多极算法求解的网格单元总数为41 961，多极和局部展开的阶数为10，迭代收敛误差为1×10-5。上述2种方法计算的电位分布云图见图10，管道电位分布见图11。可以看到，这2种方法获得的云图结果相似，管道电位分布保持一致。该算例表明，基于线单元的快速多极边界元法能够高效进行大规模复杂阴极保护系统分析，并能保持较高精度。
|
| 图 10 整体电位分布 |
|
| 图 11 阴极保护系统中间管道电位分布 |
本文将快速多极边界元法应用于阴极保护系统分析，提出了用于离散大规模管系结构的线单元模型。数值算例表明：与采用三角形单元离散管道表面相比，该模型可将管道上的单元数量减少1个数量级，2种单元的平均相对误差约为1.4%； 快速多极算法可求解自由度数为50 000量级的大规模问题。可见，本文模型可降低网格划分难度、减少计算量； 快速多极算法可提高边界元法的求解速度，并降低存储量，从而能够在普通微机上分析大规模阴极保护系统； 迭代算法可有效处理阴极极化行为带来的非线性边界条件。
基于三维管道模型的快速多极边界元法可应用于大规模阴极保护系统分析，为石油储罐罐底等阴极保护系统设计提供依据，在腐蚀防护领域具有广泛的应用前景。
应当看到，本文提出的三维管道模型适用于大区域防腐层良好的多管系阴极保护模拟。对于局部管系问题，如局部管段防腐层漏涂点对邻近管道影响等，则仍需在局部划分三维边界单元。
| [1] | Stromman R, Arild R. Computerized numerical techniquesapplied in design of offshore cathodic protection systems [J]. Materials Performance, 1981, 20(4): 15-20. |
| [2] | 张鸣镝, 杜元龙, 殷正安, 等. 有限差分法计算海底管道阴极保护时的电位分布 [J].中国腐蚀与防护学报, 1994, 14(1): 77-81.ZHANG Mingdi, DU Yuanlong, YIN Zheng'an, et al. Calculating potential distributions of cathodically protedted subsea pipeline with finite difference method [J]. Journal of Chinese Society for Corrosion and Protection, 1994, 14(1): 77-81.(in Chinese) |
| [3] | Kranc S C, Alberto A S. Detailed modeling of corrosion macrocells on steel reinforcing in concrete [J]. Corrosion Science, 2001, 43(7): 1355-1372. |
| [4] | Hassanein A M, Glass G K, Buenfeld N R. Protection current distribution in reinforced concrete cathodic protection systems [J].Cement and Concrete Composites, 2002, 24(1): 159-167 |
| [5] | 邱枫, 徐乃欣. 钢质贮罐底板外侧阴极保护时的电位分布 [J]. 中国腐蚀与防护学报, 1996, 16(1): 29-36.QIU Feng, XU Naixin. Potential distribution on cathodically protected external tank bottom [J]. Journal of Chinese Society for Corrosion and Protection, 1996, 16(1): 29-36. (in Chinese) |
| [6] | 邱枫, 徐乃欣.码头钢管桩阴极保护时的电位分布 [J]. 中国腐蚀与防护学报, 1997, 17(1): 12-18.QIU Feng, XU Naixin. Potential distribution on cathodically protected steel pipe piles [J]. Journal of Chinese Society for Corrosion and Protection, 1997, 17(1): 12-18.(in Chinese) |
| [7] | 邱枫, 徐乃欣. 用带状牺牲阳极对埋地钢管实施阴极保护时的电位和电流分布 [J].中国腐蚀与防护学报, 1997, 17(2): 106-110.QIU Feng, XU Naixin. Potential and current distributions on pipelines cathodically protected with ribbon sacrificial anodes [J]. Journal of Chinese Society for Corrosion and Protection, 1997, 17(2): 106-110.(in Chinese) |
| [8] | Fu J W, Chow J S K. Cathodic protection designs using an intergral equation numercial method [J]. Materials Performance, 1982, 21(10): 9-12. |
| [9] | Degiorgi V G. Evaluation of perfect paint assumptions in modeling of cathodic protection systems [J]. Engineering Analysis with Boundary Elements, 2002, 26(5): 435-445. |
| [10] | Riemer D P, Orazem M E. A mathematical model for the cathodic protection of tank bottoms [J].Corrosion Science, 2005, 47(3): 849-868. |
| [11] | Brichau F, Deconinck J. A numerical model for cathodic protection of buried pipes [J]. Corrosion, 1994, 50(1): 39-49. |
| [12] | Rokhlin V. Rapid solution of integral equations of classical potential theory [J]. Journal of Computational Physics, 1985, 60(2): 187-207. |
| [13] | Greengard L, Rokhlin V. A fast algorithm for particle simulations [J]. Journal of Computational Physics, 1997, 135(2): 280. |
| [14] | Greengard L, Rokhlin V. A new version of the fast multipole method for the Laplace equation in three dimensions [J]. Acta Numerica, 1997, 6(1): 229-269. |
| [15] | Hrycak T, Rokhlin V. An improved fast multipole algorithm for potential fields [J]. SIAM Journal on Scientific Computing, 1998, 19(6): 1804-1826. |
| [16] | Nishimura N, Liu Y J. Thermal analysis of carbon-nanotube composites using a rigid-line inclusion model by the boundary integral equation method [J]. Computational Mechanics, 2004, 35(1): 1-10. |
| [17] | Wang H T, Yao Z H. Application of a new fast multipole bem for simulation of 2D elastic solid with large number of inclusions [J]. Acta Mechanica Sinica, 2004, 20(6): 613-622. |
| [18] | Liu Y, Nishimura N, Otani Y. Large-scale modeling of carbon-nanotube composites by a fast multipole boundary element method [J]. Computational Materials Science, 2005, 34(2): 173-187. |
| [19] | Liu Y J, Nishimura N, Otani Y, et al. A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model [J]. Journal of Applied Mechanics, 2005, 72(1): 115-128. |
| [20] | Wang H T, Yu S Y. Large-scale numerical simulation of mechanical and thermal properties of nuclear graphite using a microstructure-based model [J]. Nuclear Engineering and Design, 2008, 238(12): 3203-3207. |
| [21] | Zhu X, Chen W, Huang Z, et al. Fast multipole boundary element analysis of 2D viscoelastic composites with imperfect interfaces [J]. Science China Technological Sciences, 2010, 53(8): 2160-2171. |
| [22] | Wang H T, Yao Z H. Large-scale thermal analysis of fiber composites using a line-inclusion model by the fast boundary element method [J]. Engineering Analysis with Boundary Elements, 2013, 37(2): 319-326. |
| [23] | Nishimura N. Fast multipole accelerated boundary integral equation methods [J]. Applied Mechanics Reviews, 2002, 55(4): 299-324. |
| [24] | Yoshida K. Applications of Fast Multipole Method to Boundary Integral Equation Method [D]. Kyoto, Japan: Kyoto University, 2001. |
| [25] | 张东东, 杜翠薇, 程学群, 等. 快速多极边界元法在阴极保护求解电位分布中的应用 [C]// 第十二届中国科学技术协会年会, 第二卷.福州: 中国科学技术协会, 2010: 601-608.ZHANG Dongdong, DU Cuiwei, CHEN Xuequn, et al. Fast multipole boundary element methods for cathodic protection potential distribution problems [C]// The Twelfth Annual Conference of CAST, Volume 2. Fuzhou: China Association for Science and Technology, 2010: 601-608.(in Chinese) |