基于核函数理论的快速多极展开及其算例研究

刘 青，李世俊，武 伟,†

(1．西北工业大学航天学院，陕西西安 710072；2．宝鸡市158信箱，陕西宝鸡 721000)

边界元法是近年来受到广泛关注的一种数值方法[1]．它将控制方程转化为边界积分方程，通过求解边界上的未知量获得实际问题的解．它将问题降低一维，形成的方程自由度小，计算量降低，处理方便，求解精度高．因此，近年来边界元法在许多领域得到了广泛应用[2]．然而，常规边界元法得到的系数矩阵为满阵，计算量随自由度数目的增加而呈平方量级增长，因此，难以应用于大规模工作问题中．为解决该问题，近三十年来，人们提出了多种快速算法，如快速多极方法（Fast Multipole Method, FMM）[3-5]、H-矩阵法[5]、快速小波方法[6-7]、PFFT（Precorrected Fast Fourier Transform）[8]、ACA（Adaptive Cross Approximation）[9]等．这些算法可以应用于加速边界元法的求解，形成的快速边界元法可以将边界元法的计算复杂度由O(N2)降低到O(NlogαN)α≥0[10-14]．

快速多极方法因其优点和影响力，已经广泛应用于静电学、声学、电磁学、弹性动力学等许多领域[10-13,15]．它的基本思路是，首先将边界逐次剖分成越来越小的区域，然后对每个小区域，根据核函数的远场低阶解析展开，依次采用S2M（Source to Multipole）、M2M（Multipole to Multipole）、M2L（Multipole to Local）、L2L（Local to Local）、L2T（Local to Target）变换，计算它与远场区域之间的相互作用，最后实现统一．其优点是，可以达到任意指定的精度，而且计算效率高，可以将边界元法的计算量降低到O(N)．但是，缺点是它的各步变换公式是由核函数的解析展开式推导而来，而不同的展开方式对计算效率有较大的影响，且不同物理问题的核函数往往不同，解析展开式也不相同．因此，对于每一个新型物理问题的核函数，都需要重新推导它的低阶解析展开式以及各步变换的计算公式，以构造FMM快速计算格式，这为FMM向不同领域的推广带来了很大的不便．

为解决上述基于核函数解析展开的FMM（下面简称“解析FMM”）的这一缺陷，本文将以文献[16]提出的核无关快速多极方法（Kernel Independent Fast Multipole Method, KIFMM）为例，研究探讨KIFMM与解析FMM之间的关系，进一步深入研究KIFMM中等效源近似与解析FMM中核函数解析展开的关系，从而将KIFMM统一到传统解析FMM的理论框架下，并在此基础上，开展应力计算实例，验证本文理论的正确性．

1 解析FMM和KIFMM原理

基于等效源的KIFMM与解析FMM方法的结构相同，它们都需要构建八叉树（对于二维问题，构建四叉树）结构，在八叉树上进行上行和下行运算．首先，定义第0层的正方体格子C0，使所有源点和目标点均在C0内，然后，将C0均分成8个子格子，这样就得到了八叉树的下一层，即第1层的格子，C0即为这些格子的父格子．一直细分下去，直到八叉树中每个叶子格子中源点个数和目标点个数均不超过预先给定的值s．对于每一个格子C，定义它的近场NC为同层上与它相邻格子以及它自身所在的区域，其余区域称为它的远场FC．对格子与其远场之间的相互作用KIFMM和FMM都是采用上行、下行变换进行快速近似计算，但KIFMM与FMM最显著的不同之处在于，它采用上行等效源代替解析FMM中在源点附近的多极展开，采用下行检测势代替解析FMM中目标点附近的局部展开．本节首先分析解析FMM和KIFMM的加速计算原理．

1.1 解析FMM

解析FMM的远场计算步骤如图1所示．

图1 解析 FMM 中的远场变换计算

假设需要计算某格子内N个源点{yj}点处的源强度{qj}对它某个远场格子内M个目标点{xi}点处位势的贡献{pi}，

其中G(x,y)为核函数．FMM首先将核函数分别在源点附近某点0y处进行低阶多极展开，在目标点附近某点0x处进行低阶局部展开，可得：

其中，展开项数T()ε与选用的精度ε有关．这样，式（1）可以通过如下3步计算：

1）根据源点的源强度qj，计算y0处的多极展开系数ms(y0)．

即FMM中的 S2M 变换．

2）根据y0处的多极展开系数ms(y0)计算x0处的局部展开系数lt(x0)．

即FMM中的M2L变换．

3）根据x0处的局部展开系数lt(x0)计算目标点处的位势pi．

比较两种计算方法可知，采用式（1）直接计算{pi}，计算量为O(MN)，而若采用单层FMM计算，则计算量为O(T(ε)⋅(M+N+T(ε))．当M、N远大于T(ε)时，采用单层FMM计算即可显著减少计算量．对于边界元法来说，单层FMM最多可将O(N2)的计算量降低到O(N3/2)．对于大规模问题，为获得O(N)的计算复杂度，则需要采用多层FMM算法，即引入M2M变换和L2L变换，其中，M2M变换将低层格子上的多极展开系数变换为高层格子上的多极展开系数，L2L变换将高层格子上的局部展开系数变换为低层格子上的局部展开系数，它们可以采用插值的方法计算[17]．多层FMM的计算步骤如图1所示．从图1中也可以看出，上行运算的S2M、M2M和下行运算的L2T、L2L基本上是对称的．

1.2 KIFMM

KIFMM的算法结构与解析FMM类似，但远场作用采用等效源和检测势来近似计算．单层KIFMM计算过程如图2所示．

图2 单层KIFMM计算步骤（实心点表示等效点，空心圆环表示检测点）

其中，T()ε为上行等效点的个数，T′()ε为上行检测点的个数，它们均与选用的精度有关，写为矩阵形式，得：

对应于图2中的第1步、第2步。

对应于图2中的第5步．

为保持KIFMM中上行运算和下行运算类似于解析FMM的对称性，可以将第1步和第2步称为S2M变换，它将格子C中的源变换为上行等效源，变换矩阵为：

矩阵维数为M×T(ε)．所以，采用这种方法计算量为O(T(ε)⋅(M+N+T(ε))，与解析FMM相同．

对于大规模问题，KIFMM 同样可以通过引入M2M变换和L2L变换，采用多层算法格式以获得O(N)的计算复杂度．其中，M2M变换将低层格子上的上行等效源变换为高层格子上的上行等效源，计算公式与S2M类似，只需要将公式（6）中源点上的源替换为低层格子上的上行等效源即可．L2L变换将高层格子上的下行检测势变换为低层格子上的下行检测势，计算公式与L2T类似，只需要将公式（10）中目标点上的位势替换为低层格子上的下行检测势即可．多层KIFMM的远场计算变换步骤如图3所示．

图3 KIFMM的远场变换计算

2 KIFMM中等效源近似与核函数展开

本节研究KIFMM中等效源近似与解析FMM中核函数解析展开的关系，讨论等效源近似对应的核函数展开形式，将两者统一到核函数展开理论的框架下．

首先观察利用上行等效源计算C内的源在其远场产生的位势场，即图2中第1步和第2步．C内的任意一点y处单位源在远场任意点x处产生的位势场p(x)为：

根据KIFMM，采用上行等效源计算远场处的位势场p(x)，需要进行如下3步计算．

所以，构造上行等效源，实际上相当于采用函数（20）对核函数在上行等效点处进行T()ε项近似低阶展开．

类似地，采用格子D的下行等效源，计算它远场DF中任一点y处单位源在D内任一点x点处产生的位势，可得：

所以，下行计算中，实际上相当于采用函数（23）对核函数在下行等效点处进行T()ε项近似低阶展开．

综合式（21）和（24），C中任一源点y上单位源在它远场格子D中任一目标点x上产生的位势可表示为：

对比式（25）和解析FMM中核函数的解析展开式（2）可知，KIFMM的等效源近似相当于采用（25）式对核函数进行近似低阶展开．

3 算 例

本节以两个弹性场的计算为例来验证本文理论和方法的正确性．

由于弹性动力学核无关快速多极边界元法（Kernel Independent Fast Multipole Boundary Element Method, KIFMBEM）的性能已经在文献[18]中进行过详细讨论，因此本文对此不再重复，仅研究KIFMBEM计算应力应变场的精度．程序采用C++串行编程实现．Nyström方法离散得到的方程用GMRes求解器求解，采用对角块预处理技术以加速收敛．方程组迭代收敛限与选用的算法精度ε相同．

首先采用一个有解析解的算例来验证本文理论和方法的正确性．假设在无限大弹性体中存在一个球形空腔，空腔半径为1 m，内壁面受幅值为100 MPa、频率为2 672 Hz的简谐压强作用，如图4所示．弹性体弹性模量为200 GPa，密度为7 850 kg/m3，泊松比为0.31．计算弹性体内的位移场、应变场和应力场．

图4 受内压的球形空腔

首先在球面上划分网格，求解位移边界积分方程得到边界上的位移和面力．创建外径为3 m、内径为1 m的圆环面，在圆环面上按同样单元大小划分网格，计算网格点上的应力应变，并与解析解对比计算其误差．采用不同大小的单元划分一组网格以观察应力应变场计算结果的收敛性．球面网格单元数分别为48、140、542、1 138、1 920、2 916、11 864和48 192，圆环面上的节点数分别为32、92、356、732、1 032、1 664、7 296和29 186．由于在基于二次元的Nyström离散方法中，每个边界单元上有6个Gauss点，每个Gauss点上有3个自由度，所以边界上的自由度为单元数乘以18．由于应力应变张量的对称性，在求解应力边界积分方程时，圆环面上每个点有6个自由度，所以圆环面上总自由度数为节点数乘以6．

本算例在一台8核（Xeon 5355, 2.66 GHz）、32 GB内存的工作站上计算．在KIFMM中采用不同的上行等效点数p（p= 4,5,6,7），计算了应力应变场，得到的误差与边界上自由度N的关系如图5和图6所示，其中，conventional是用未加速的传统边界元法计算的结果．可以看出，在选用合适的等效点数时，采用本文的KIFMM方法可以有效地保持传统边界元法的精度，表明了本文理论和方法的正确性．

图5 圆环面上的应力误差

图6 圆环面上的应变误差

4 结 论

本文推导证明了基于等效源展开的核无关快速多极方法（KIFMM）等同于采用了一种特殊的核函数展开，从而它与解析FMM可以统一到相同的理论框架下，证明了KIFMM中等效源近似与解析FMM中核函数解析展开相似．在此基础上，将此理论应用于弹性动力学的应力应变场分析，验证了本文理论和方法的正确性．