边界元通用大规模快速算法研究.doc

边界元通用大规模快速算法研究王鹏 王海涛 于溯源（清华大学核能与新能源技术研究院， 北京， 100084）摘 要：边界元法已经广泛应用于声场、电磁场和温度场等领域的数值计算。受计算效率的限制，传统边界元法不适合求解大规模问题。本文将ACA（Adaptive Cross Approximation）算法用于边界元法的大规模快速求解，分析使用ACA算法的边界元求解计算复杂度和求解流程，数值研究ACA算法的计算精度和适用的求解范围，并与其他算法比较。结果表明，ACA算法与传统求解算法相比，在求解效率上有数量级的提高，同时可以控制精度，能够在单台普通微机上完成大规模复杂结构的边界元数值仿真。关键词：边界元；Adaptive Cross Approximation；大规模；声场；电磁场Research of Generally-purposed Large-scale Fast Algorithm for the Boundary Element MethodPeng WANG Haitao WANG Suyuan YU(Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing, 100084, China)Abstract：The boundary element (BEM) has been widely used in the numerical simulation of the acoustic, electromagnetic and temperature fields. The conventional BEM is not suitable for large-scale problems because of the poor computational efficiency. In this paper, the ACA (Adaptive Cross Approximation) algorithm is introduced for the large-scale fast solution of the BEM. The computational complexity and performance procedure of the BEM solution by use of the ACA algorithm are analyzed, and the corresponding accuracy and computational scope is studied numerically in comparison with other algorithms. The results show that the efficiency of ACA is at least one order of magnitude higher than that of the conventional algorithms and that the accuracy is controllable. The BEM simulation of large-scale complex structures can be performed on one desktop computer by use of this new algorithm.Key words：boundary element；Adaptive Cross Approximation；large-scale; acoustic field；electro-magnetic field60 引言与有限元法相比，边界元法只需在结构表面离散网格1，而且易于处理无限域，因此适用于声场、电磁场、温度场等多种物理场的数值计算。例如在结构振动噪声仿真中，边界元法可用于内域和外域声场计算，并可与有限元结构计算耦合分析2, 3。传统边界元法形成的线性代数方程组AX = B中的系数矩阵A通常是满阵，对某些物理场问题如弹性力学，A也是非对称的。满阵的存储量级是，其中N是边界元方程的未知量个数。使用直接算法如高斯消去、LU分解等求解该方程组所需计算量级为；使用迭代算法如共轭梯度(CG)、广义极小残差法(GMRES)等求解计算量为，其中k是迭代收敛步数，通常近似看做是远小于N的常数。可以看出，当N增加时，传统边界元法无论使用直接解法还是迭代解法，其存储和计算量均呈2次或更高的增长速度，将很快达到计算机的求解能力上限。因此传统边界元法不适合处理大规模问题。为提高边界元法的求解效率，将边界元法与求解具有类似特征物理场的快速算法结合成为近年来的边界元主要研究方向之一，其中最有代表性的是快速多极边界元法。边界单元之间的相互作用类似于有势场中粒子之间的相互作用，如天体之间的万有引力、带电粒子之间的静电力等，这种相互作用的特点是与作用双方的距离呈反比。N个粒子的相互作用的直接计算量为。1986年，Barnes和Hut提出树代码算法(Tree Code)4，将多粒子场的计算量降至。1987年，Greengard和Rokhlin提出快速多极算法(Fast Multipole Method, FMM)5，将此类问题的计算量进一步降至。1997年，Greengard和Rokhlin在1987版快速多极算法的基础上进一步提出新型快速多极算法，引入了指数展开的概念，将该算法在三维问题的求解速度大幅度提高6。从上世纪90年代末至今，将快速多极算法和边界元法结合的快速多极边界元法的研究有很大的进展，并用于声场、电磁场、弹性场和温度场等多个物理领域的大规模数值仿真7-14。快速多极边界元法的求解效率与边界元的未知量呈线性关系，因此使边界元法能够求解大规模问题。但该算法的难点是与特定物理问题相关，需要针对不同物理场的基本解寻找特定的展开格式。2003年，Bebendorf和Rjasanow提出了ACA(Adaptive Cross Approximation)算法15，用于对秩很小的矩阵进行快速向量内积分解和存储。考虑到使用树结构递归分解求解域后每两个相距较远的树节点包含的边界单元集合形成的系数子矩阵的秩很小，因此可以将ACA算法用于递归的子矩阵的快速近似和存储，从而可以用于边界元大规模快速求解，而且与物理背景无关。本文介绍声场、电场和温度场等多类物理场的边界积分方程描述格式，将ACA算法用于边界元法的大规模快速求解，分析使用ACA算法的边界元求解计算复杂度和求解流程，数值研究ACA算法的计算精度和适用的求解范围，并与其他算法比较。结果表明，ACA算法与传统GMRES迭代法相比，在计算速度上显著提高，同时所需的存储量大大减少，而且可以控制计算精度。与快速多极算法相比，ACA算法在中等规模以下问题的求解速度也具备优势。使用ACA算法能够在单台普通微机上实现大规模复杂结构的边界元数值仿真。1 多类物理场的边界积分方程描述作为物理场的微分方程描述的两种等价形式，变分方程和积分方程分别构成有限元和边界元的数学基础。后者通常只描述求解域的边界物理量，因此作为其离散形式的边界元只需要对边界生成网格，降低了有限元法的域内离散网格难度，减少了求解未知量个数。忽略体积分的边界积分方程描述为：(1)其中S是求解域的边界。p和q是边界物理量，在声场分析中分别表示边界声压和声压法向导数；在静电场分析中分别表示边界电位和法向电流；在温度场分析中分别表示边界温度和热流。x和y分别代表边界的源点和场点，是与x所在边界几何形状有关的函数，对光滑边界，；和称作边界积分方程的基本解，与具体物理问题有关。在声场分析中，简谐结构振动产生的辐射声压p满足波动方程：(2)其中是波数，w是圆频率，c是流体介质声速。声场问题对应的边界积分方程基本解为：(3)式中r是x和y的距离，i是纯虚数。静电场和稳态温度场对应的物理方程描述相同。在静电场/稳态温度场分析中，无内源的电位/温度p满足Laplace方程：(4)静电场/稳态温度场问题对应的边界积分方程基本解为：(5)式中k在静电场中对应电导率，在稳态温度场分析中对应热导率。将边界S离散成M个单元，每个单元上有L个结点，边界单元上物理量的插值公式分别为：(6)其中是定义在二维局部坐标上的插值函数。将式(6)代入式(1)，得到离散形式的边界积分方程：(7)利用配点格式的加权余量法，即将源点x设为边界单元上的结点，并通过边界单元上的数值积分，引入边界条件，就得到最终的线性代数方程组AX = B。2 边界元的ACA算法传统边界元法的系数矩阵A是满阵，A中的每一个元素都需要显式存储和参与运算。ACA算法的基本原理是：通过树结构将矩阵A递归分解为一系列子矩阵的集合，每个子矩阵对应一组源点集合（共个）和一组场点集合（共个）；当源点集合与场点集合的距离足够远时，的秩很低，可以用一组向量内积的和进行低秩近似（Low-rank approximation）15：(8)其中和是向量，长度分别是和。k是向量个数，与低秩近似控制残差有关，通常。可以看出，通过ACA算法的低秩近似，原子矩阵的存储量降低至向量组的存储量，而且子矩阵-向量相乘中涉及的运算量也进一步降低：(9)使用ACA算法进行边界元求解的步骤如下。步骤一：使用树结构对求解域进行递归分解。以二维问题为例，首先使用一个正方体涵盖求解域边界上的所有边界单元，该正方体称为树结构的根节点（root）；将正方体分解为4个相同的子正方体，每个子正方体均包含一定数量的边界单元；对每个子正方体继续递归分解，直至分解后的子正方体中包含的边界单元数量小于预先设定的下限，就停止分解，不再分解的子正方体称为叶子节点（leaf）。分解后的树结构见图1。图1 求解域边界的树结构递归分解步骤二：矩阵分解。生成树结构后，利用树结构对矩阵进行分解。首先定义树结构节点之间的关系，对树结构中的某个节点ND，其邻居节点定义为与ND大小相同、和ND共享至少一个角点的节点；其相互作用节点定义为与ND大小相同、不为邻居、但其父节点和ND的父节点是邻居的节点。节点关系的示意见图2。遍历整个树结构，对每个节点，将其中包含的边界单元作为源点，其邻居节点包含的边界单元作为场点，二者作用生成的子矩阵定义为；其相互作用节点包含的边界单元作为场点，二者作用生成的子矩阵定义为。对所有节点计算得到的和即构成定义边界元整体系数矩阵A所需的所有子矩阵。图2 节点关系定义步骤三：子矩阵ACA近似。对所有的子矩阵，根据式(8)进行低秩近似。步骤四：迭代计算。使用迭代算法如GMRES进行边界元求解的过程中，在每个迭代步，系数矩阵A-迭代向量的相乘运算由A的一系列子矩阵和与迭代向量对应的局部子向量相乘运算代替。其中与子向量的相乘运算根据式(9)由向量内积运算近似。3 数值算例数值算例用于验证边界元的ACA求解算法的计算精度和计算复杂度，并与传统的GMRES迭代法和快速多极边界元法比较。边界元程序使用标准C+编写，求解方程以Laplace方程为例，可用于静电场和稳态温度场的模拟。算例运行环境为1台配备Intel Core2 Duo E8400处理器的个人桌面电脑，内存为2GB。算例模型是一个内部均布8个立方体空洞的大立方体，大立方体的一个端面给定均匀电位/温度，相对的另一个端面给定均匀电流/热流，其余表面绝缘/绝热。模型的三维透视图见图3。边界单元类型选择三角形常值面单元。图3 计算模型三维透视简图3.1 ACA计算精度验证为了研究边界元法ACA求解的计算精度，分别使用传统的直接矩阵-向量相乘法和ACA法计算边界元形成的方程组中的右手向量B。使用BD表示直接矩阵-向量相乘运算得到的结果向量，BA表示ACA法近似计算得到的结果向量。二者的差别用相对残差表示为：(10)其中2是向量的2范数运算符。ACA算法通过设置低秩近似控制残差控制矩阵的近似程度。图4给出了模型边界元离散自由度为7200时随ACA低秩近似控制残差的变化关系。可以看出，随着的减小，也迅速减小，与呈很好的对数线性关系。因此使用ACA算法进行边界元求解的精度可控。图4 ACA计算精度3.2 ACA计算复杂度研究为了研究边界元法ACA求解的计算复杂度，包括计算量和存储量，对模型划分不同尺寸的网格，对应100050000的自由度范围；分别使用ACA法、传统的GMRES迭代法和新型快速多极算法进行计算，比较不同算法的计算速度和存储量随自由度的变化。图5是计算速度随自由度的变化曲线，图6是存储量随自由度的变化曲线。图5 计算量与自由度的关系曲线图6 存储量与自由度的关系曲线从图5可以看出，在算例求解的自由度范围内，ACA法的求解速度均高于传统的GMRES法。当自由度超过10000时，受存储量的限制，传统GMRES法不能计算更大规模，但根据ACA法和GMRES法的计算量曲线趋势可以看出，随着自由度的增加，ACA法的求解速度与GMRES法相比有量级上的提高。与快速多极算法相比，ACA法在30000自由度以下的中低规模问题上求解速度占有优势，两类算法的平衡交叉点约为30000自由度。当自由度超过30000时，由于快速多极算法的计算效率是，其计算速度将超过ACA法。从图6可以看出，ACA算法的存储量介于GMRES法和快速多极算法之间，而且根据ACA法和GMRES法的存储量曲线趋势可以看出，随着自由度的增加，ACA法的存储量与GMRES法相比有量级上的减少。快速多极边界元法的存储量与ACA法相比又有显著减少。根据各算法存储量曲线的趋势，对于2GB的内存配置，传统GMRES迭代法的求解自由度上限约15000；ACA法约为55000，边界单元数量达到该自由度的问题已属大规模问题；而快速多极算法可达几十万。由于与物理问题类型无关，ACA法可方便用于其他物理问题的边界元求解。4 结论边界元法只在结构表面离散单元，易于处理无限域，在模拟声场、电磁场等物理问题时具备独特优势。传统边界元法受系数矩阵为满阵的限制，求解效率低下；快速多极算法的出现使边界元法可以在单台普通微机上处理大规模问题，但该算法与特定物理背景有关；ACA法基于矩阵的近似计算，与物理背景无关，具备作为边界元通用求解器的潜力。本文回顾不同物理场的边界积分方程描述，将ACA算法应用于边界元法的求解，分析使用ACA算法的边界元求解计算复杂度和求解流程，数值研究ACA算法的计算精度和适用的求解范围，并与传统的GMRES迭代算法和快速多极边界元法比较。结果表明，ACA算法与传统求解算法相比，在计算速度上显著提高，同时所需的存储量大大减少，而且可以控制计算精度。与快速多极算法相比，ACA算法在中等规模以下问题的求解速度也具备优势。使用ACA算法能够在单台普通微机上完成约50000边界自由度的大规模复杂结构的边界元数值仿真。参考文献:1 Yao ZH, Kong FZ, Wang HT, Wang PB. 2D Simulation of composite materials using BEMJ. Engineering Analysis with Boundary Elements, 2004, 28: 927935.2 Gaul L, Wenzel W. A coupled symmetric BE-FE method for acoustic uid-structure interactionJ. Engineering Analysis with Boundary Elements, 2002, 26(7): 629636.3 Citarella R, Federico L, Cicatiello A. Modal acoustic transfer vector approach in a FEM-BEM vibro-acoustic analysisJ. Engineering Analysis with Boundary Elements, 2007, 31(3): 248258.4 Barnes J, Hut P. A hierarchical O(NlogN) force calculation algorithmJ. Nature, 1986, 324: 446-449.5 Greengard L, Rokhlin V. A fast algorithm for particle simulationsJ. J. Comput. Phys., 1987, 73: 325348.6 Greengard L, Rokhlin V. A new version of the fast multipole method for the Laplace equation in three dimensionsJ. Acta Numer., 1997, 6: 229270.7 Fischer M, Gauger U, Gaul L. A multipole Galerkin boundary element method for acousticsJ. Engineering Analysis with Boundary Elements, 2004, 28(2):155162.8 Fischer M, Gaul L. Fast BEM-FEM mortar coupling for acoustic-structure interactionJ. International Journal for Numerical Methods in Engineering, 2005, 62(12): 16771690.9 Nemitz N, Bonnet M. Topological sensitivity and FMM-accelerated BEM applied to 3D acoustic inverse scatteringJ. Engineering Analysis with Boundary Elements, 2008, 32(11): 957970.10 Aoki S, Amaya K, Urago M, Nakayama A. Fast multipole boundary element analysis of corrosion problems. CMESJ. C