等几何边界元法

第1章绪论 　　1.1引言 　　实际中的各种工程问题都需要用解析或数值方法来预测其力学行为。由于解析方法适用范围有限，所以有必要发展数值方法来获得问题的近似解。常用的数值方法有很多，比如有限元法和边界元法。有限元法已被证明是一种可行的方法，被广泛应用于解决各种工程问题，其发展已趋于成熟。边界元法在一些领域里具有独*的优势，比如弹性力学[1]、断裂力学[2]、接触力学[3]、形状优化[4]和声学[5]等。然而，在进行计算机辅助工程(CAE)分析之前，有限元法和边界元法需要将计算机辅助设计(CAD)连续参数的几何模型通过网格生成器离散为分析时的计算模型，导致CAD几何模型和CAE分析模型不一致，其几何离散误差使得计算精度降低[6]。为此，Hughes等[6]于2005年提出了等几何分析方法，其核心思想是利用非均匀有理B样条(NURBS)基函数同时表达CAD几何模型和CAE分析模型，实现了CAD和CAE的无缝连接。等几何分析方法完全消除了网格划分的概念，其有助于实现精确的几何实体造型，并可以得到更好的数值解，也节省了大量的计算时间。 　　本章首先简要介绍等几何有限元法的一些工作，然后着重介绍等几何边界元法及其应用，*后简单介绍本书的章节内容。 　　1.2等几何有限元法简介 　　等几何分析技术[6]的出现引起了许多研究者的关注，并应用于各种问题[7]，例如：①接触问题的等几何分析能够得到更精确的结果，原因在于其基于NURBS基函数固有的高阶连续性，实现了接触表面的光滑表示[8-10]；②NURBS基函数的高阶连续性使其能够以非常简单和统一的方式构建C1或更高阶的近似模型，特别适合于薄板和薄壳的有限元分析，而且与标准有限元板壳单元相比，等几何板壳单元表现出更不明显的剪切锁定[11，12]；③结构振动问题的等几何分析比传统有限元方法更有优势，已有结果表明，在波传播问题中，k细化法比高阶有限元p法具有更强的鲁棒性和更准确的频谱[13，14]；④在优化问题中，由于设计模型直接用于等几何分析，实现了设计模型与分析模型之间的紧密耦合[15，16]。NURBS在CAD和计算机图形学中是普遍存在的，但从计算几何的角度来看，它存在很多问题。 　　也许遇到的*大困难是NURBS模型通常是由很多片组成，这些片之间并不是无缝拼接的，这通常会使网格生成复杂化。从分析的角度来看，NURBS的张量积结构是低效的，无法进行局部加密，这会导致低效的误差估计和自适应算法。为了克服NURBS模型的缺点，一些学者相继在等几何分析的研究工作中引入了T样条[17，18]、PHT样条[19，20]和LR样条[21]等。这些样条具有局部加细的特点，已应用于多个领域，如弹性力学[22]、形状优化[23]、断裂力学[24-27]、动力学分析[28-31]等。 　　以上只是对等几何有限元法做了简单介绍，它具有广泛的应用前景。有兴趣的读者可以关注计算力学的一些知名期刊，例如：Computer Methods in Applied Mechanics and Engineering、International Journal for Numerical Methods in Engineering及Computational Mechanics。这些期刊时常发表一些新的等几何分析的研究成果。另外，一些有关等几何分析的专著[32-35]及评述性的论文[36-38]也很有参考价值。 　　1.3等几何边界元法简介 　　边界元模型只需要边界几何信息，从而使CAD与边界元法实现真正的无缝对接，而有限元法所要求的体积表示与CAD中几何模型的边界表示使得等几何有限元法(IGABEM)在分析三维实体模型时面临挑战。鉴于此，许多学者对等几何边界元法产生了兴趣。Politis等[39]于2009年首次将等几何概念引入到边界元法中，并用之解决外场黎曼问题。Li和Qian[40]提出了一种基于边界积分的等几何分析和形状优化方法，利用NURBS基函数来表示边界形状和近似分析中的物理场，该方法已经成功地应用于解决弹性和位势问题。Simpson等[41]将等几何边界元法应用于求解二维弹性力学问题，详细介绍了等几何边界元法具体实施中涉及的配点的选取方法和积分方程中奇异积分的求解方法。Beer等[42]在其等几何边界元法专著中介绍了稳态势问题、弹性力学、非均质问题、材料非线性、黏性流动及与时间相关的问题等。张见明研究组[43]采用等几何边界元法研究了三维位势问题，数值实验表明，该方法在精度和收敛性方面都有较好的性能。与基于NURBS基函数的等几何边界元法不同，张见明提出的边界面法[44]也是一种等几何边界元法，它用于积分计算的几何量和物理量都是通过CAD模型边界表征的参数曲面得到的，已被用于弹性力学[45]、稳态和瞬态热传导[46，47]、声学[48]和弹性接触[49]等问题的解决中。本书着重介绍基于NURBS基函数的等几何边界元法的基本理论及其应用研究。 　　1.3.1奇异积分及拟奇异积分的计算 　　由于等几何边界元法将等几何分析中的NURBS基函数引入到传统的边界元 　　法中，使得该方法既继承了边界元法和等几何分析的优点，同时又继承了它们的弱点。各类奇异积分的存在是影响边界元法和等几何边界元法的计算精度和计算效率的重要因素，处理这些奇异积分需要特殊的计算方法。在这些特殊的方法中，Guiggiani等[50]提出的局部规则化方法受到众多边界元法领域里学者的青睐，已被用于处理各种奇异积分，但这种方法需要将被积函数中的所有量都进行泰勒级数展开，而且只保留至二阶精度，数值计算繁琐。高效伟[51]提出的基于径向积分法的高阶奇异曲面积分的直接计算方法具有一些优点，例如：不需要对被积函数中的每个量进行泰勒级数展开；具有可以采用级数展开的高阶项以及便于程序编制。 　　对于等几何边界元法，Simpson等[41]讨论了弹性问题等几何边界元法中奇异积分的计算，其中的弱奇异积分通过Telles变换法[52]进行处理，强奇异积分则是通过Guiggiani等[53]提出的奇异性分离法解决。上述奇异积分求解技术也应用于声学等几何边界元法[54，55]和液体夹杂等几何边界元法[56]分析。公颜鹏等[57]采用高效伟提出的幂级数展开法[51]计算了势问题等几何边界元法中的奇异积分，并与传统边界元法进行了比较，结果表明等几何边界元法具有较高的计算精度。Simpson等[41]认为传统边界元法中的刚体位移法不适用于计算等几何边界元法中的强奇异积分，但事实并非如此。徐闯等[58，59]通过简单变换法运用常温法和刚体位移法间接求解了热传导和弹性动力学等几何边界元法中的强奇异积分。另一种避免求解奇异积分的方法是基于规则化等几何边界积分方程的方法。Heltai等[60]提出了一个三维Stokes流动问题的非奇异等几何边界积分方程，并对其收敛性进行了数值验证。在数值实现中，标准高斯求积法足以对规则化等几何边界积分方程进行积分。Simpson等[54]针对声学低频问题提出了一个基于规则化的Burton-Miller公式，其中所有积分都是弱奇异性的，该公式保证了声学等几何边界元法对所有波数的稳定性。 　　等几何边界元法在处理薄体或涂层结构时会遇到拟奇异积分，其精确计算有待进一步研究。公颜鹏等[61]将指数变换法[62，63]推广到等几何边界元法中，成功地解决了等几何边界元法的边界层问题。随后，公颜鹏和董春迎[64]又将自适应积分法[65]引入到三维势问题的等几何边界元法中。由于自适应积分方法是基于单元细分的思想，随着拟奇异性的增强，通过单元细分得到的子单元数量会急剧增加。因此，对于具有强拟奇异性的积分，该方法的计算效率将大大降低。针对这一问题以及二维、三维薄体和涂层结构，公颜鹏等[66，67]研究了拟奇异积分的计算精度和效率，提出了一种混合积分计算方法，保证了计算精度和效率之间的平衡。Keuchel等[55]采用sinh变换法计算了声学问题等几何边界元法中出现的拟奇异积分。Han等[68]利用减法技术将等几何边界元法中的拟奇异积分分离为非奇异部分和奇异部分。奇异部分的积分核用泰勒级数多项式表示，其中不同阶导数用NURBS插值。通过一系列分部积分，导出了具有近似核的奇异部分的解析表达式，而非奇异部分则采用高斯积分求解。该半解析方法能准确地计算出更靠近边界的内点的位移和应力。Han等[69]也利用半解析方法研究了二维势问题等几何边界元法中存在的拟奇异积分，并与指数法和sinh变换法进行了比较。结果表明，该方法具有竞争力，特别是在势流密度模拟过程中的近强奇异积分和高阶奇异积分计算方面。 　　1.3.2等几何边界元法的快速计算 　　等几何边界元法在继承传统边界元法优点的同时也继承了它的缺点。不适合模拟大规模问题是边界元法的一个突出弱点，这是因为由边界积分方程离散得到的线性方程组系数矩阵通常是非对称满秩矩阵，对于N自由度问题，系数矩阵的存储需要.(N2)量级。如果使用直接求解技术，如高斯消去法，计算量会达到.(N3)量级，即使使用迭代算法也需要.(N2)的计算量，这样就导致了边界元法无法解决大规模工程问题。 　　在过去30多年里，快速多极算法(FastMultipoleMethod，FMM)[70-72]、小波变换法[73]、基于傅里叶变换的算法(Fourier Transformation Based Method)[74]和分层矩阵法(Hierarchical Matrix，H-matrix)[75]等快速算法得到了广泛的发展。借助于迭代方法求解线性方程组，快速多极算法有效地加速了边界元方程的求解。在使用快速多极算法后，边界元法的计算复杂度可接近于的量级。快速多极边界元法的难点在于边界积分方程中基本解的多极、局部、指数展开格式的实现和展开系数相互之间的传递关系。不同问题的边界积分方程的基本解是不一样的，即使同一问题的二维和三维的基本解也不相同，如拉普拉斯(Laplace)方程基本解和Helmholtz方程基本解等。因此，对于不同问题，相应的快速算法实施并不是一件简单的事情。Hackbusch等[75，76]提出分层矩阵的概念，将一个稠密矩阵逐层地分成一些子矩阵，并指出其中的一些子矩阵可以被低秩矩阵很好地近似。在分层矩阵的基础上，Bebendorf[77]提出自适应交叉近似算法(Adaptive Cross Approximation，ACA)用于对秩很小的矩阵进行快速向量内积分解和存储，并将其应用于边界元法中。由于ACA是一种纯代数的矩阵压缩方法，不需要将边界积分方程的核函数解析展开，而且与物理背景无关，因此，基于ACA的分层矩阵法得到广泛的应用[78，79]。值得注意的是，快速多极边界元法和基于ACA的分层矩阵边界元法都使用迭代法求解线性方程组，如广义极小残差法。当迭代收敛的速度很快，迭代步足够小时，两种方法能够达到很高的求解效率。然而，迭代的收敛速度依赖于系数矩阵的条件数，系数矩阵的条件数越小，收敛速度越快。因此，边界元迭代算法的预处理技术经常被用于加速迭代的收敛。然而，对于许多问题即便使用了预处理技术，迭代求解的收敛速度仍然缓慢。 　　为了避免迭代算法收敛性的问题，一些学者开始研究稠密矩阵的直接快速算法[80-83]。算法的主要思想是先将矩阵逐层分解，并构造其低秩子矩阵的低秩近似，然后递归地更新方程组的解。Martinsson和Rokhlin[80]提出了基于可分离的分层块矩阵的快速直接算法，并用之对二维势问题进行了求解。Kong等[81]给出了基于分层非对角低秩(Hierarchically Off-Diagonal Low-Rank，HODLR)矩阵的快速直接算法来求解二维势问题。Lai等[82]基于HODLR矩阵结合ACA算法求解了二维高频散射问题。Huang和Liu[83]基于HODLR矩阵概念提出了一个边界元直接快速算法，并对三维势问题进行了快速求解。HODLR矩阵是分层矩阵的一种特殊形式，它的主要特点是所有的非对角子矩阵都用低秩矩阵近似。