《高等边界元法:理论与程序》共9章,第?章为绪论,第2章介绍必要的数学知识,第3~6章介绍与位势问题相关的边界元法,第7~8章介绍线性和非线性力学问题的边界元法,第9章介绍求解多种介质问题的新方法?《高等边界元法:理论与程序》展示了作者多年来的研究成果,如:将任意域积分转换成边界积分的径向积分法?求解大型非对称稀疏矩阵方程的同时消元回代法?计算弱奇异和近奇异积分的单元子分法?计算超奇异积分的幂级数展开法?建立一般非线性问题积分方程的源点隔离法以及用单一积分方程求解多种介质问题的界面积分方程法?。
更多科学出版社服务,请扫码获取。
第1章绪论
1.1数值方法概述随着科学技术的不断发展,需要解决的工程问题也越来越复杂,对于大多数问题,由于求解问题几何形状的复杂性或计算介质的非线性,人们已经很难得到问题的解析答案。另外,随着计算机性能的日益提高,求解工程问题的数值计算方法也不断成熟,现在几乎所有的大型工程问题都要借助数值计算进行分析或评估,为科技人员的工程设计提供依据或参考。
目前,已经发展起来的用得较多的数值计算方法有五大类:有限差分法、有限体积法、有限单元法、无网格法和边界元法。
(1) 有限差分法(finite difference method)[1]。有限差分法是将所考虑的区域划分成网格,用差分近似微分,把微分方程变成差分方程。也就是通过数学上的近似,把求解微分方程的问题变换成求解关于节点未知量的代数方程问题。该方法简单、易懂,便于复杂微分方程的求解,因此在流体力学领域被广泛采用[2]。但当求解问题的几何形状复杂时,按空间坐标的离散变得困难,网格的正交性不易保留,导致计算精度降低。
(2) 有限体积法(finite volume method)[3]。有限体积法是基于物理问题控制方程的积分形式的数值方法。解题思路是:把计算域离散成有限个互不重叠的控制体单元(网格),通过将积分形式的控制方程作用于每一个体单元来建立离散的代数方程组。有限体积法既有有限差分法的特点,又有有限单元法的特点,方法简单、几何适应性好,是现代计算流体力学中占主导地位的数值方法[4]。其缺点是物理量的空间导数相关量(如通量)是由周围体单元中心处的值决定的,因此精度较低,特别是靠近边界的通量值更难计算准确。
(3) 有限单元法(finite element method)[5]。有限单元法是通过变分原理建立含权函数(形函数)的体积分形式方程的方法。解题思路是:首先将计算域离散成有限个互不重叠的有一定规则的节点组成的单元,然后对每个单元积分并通过组集形成总体代数方程组。有限单元法的单元可以按不同的连接方式进行组合,每个单元可以有不同的形状和材料性质,因此有限单元法具有几何适应性强、可灵活处理不同物性参数的优点,在各个领域都得到了广泛应用[6]。有限单元法的缺点是:物理量的空间导数是通过对形函数的求导得到的,因此精度要比物理量本身低一阶;在金属成形、优化计算、渗流问题自由面的确定等运动边界问题中,有限元网格可能产生畸变和重叠,以致计算精度下降或计算中止。
(4) 无网格法(meshless method)[7]。无网格法是基于构造点插值函数的数值方法,因此只需要在计算域内布置一系列的离散点即可,不需要网格单元,具有很强的解题灵活性[8,9]。但无网格法发展得还不够成熟,缺少坚实的理论基础和严格的数学证明,因此计算精度、守恒性等一直没有明确的答案。此外,无网格法是基于点的算法,因此布点数量和方案会影响计算时间和精度。另外,由于不使用单元,对于几何较复杂的运动边界问题,边界变化时要判断重新分布后的点是内部点、外部点还是边界点,这时会存在困难。
(5) 边界元法(boundary element method,BEM)[10]。边界元法是基于格林公式和问题的基本解将控制微分方程转化为边界积分方程的一种数值方法。其主要优点是:①只需要将边界离散成单元,因而准备数据简单、便于复杂几何问题的建模[11];②能够自动满足无限远处的边界条件,因而适合于求解无限与半无限域问题[12];③所求物理量的空间导数的计算公式可以解析地从基本边界积分方程中导出, 因此所求与导数相关的物理量(如通量、应力等)与物理量本身具有同样级别的精度[13];④在求解运动边界问题时,移动边界节点的位移与原边界节点的坐标相加就自然形成了新的边界单元信息,不需要专门重构单元,也不会有网格重叠的问题[1416]。此外,由于在计算域的边界上有单元信息可用,所以通过单元积分很容易判断一点是内部、外部还是边界点。基于这些优点,边界元法在一些领域(如运动边界[1517]、裂纹[18,19]、接触[20,21]、辐射[22,23]、超薄结构[2428]、无限域和半无限域[14,29,30]、声学[31]等)的应用优于有限元法。
1.2边界元法的发展历史
边界元法已发展成为求解工程与科学问题的常用数值分析方法之一。它是一种在经典的积分方程基础上,吸收了有限元法的离散技术而发展起来的数值方法。边界元法通过采用一个满足无限或半无限域场方程的奇异函数——基本解作为权函数,将基于问题控制方程的域积分方程转化为边界积分方程,并将边界离散成边界单元来求解边界未知量的数值解。
边界元法的产生可追溯到19世纪,当时有人提出了积分等式和位势理论的数学问题,把线性偏微分方程的边值问题转化为边界积分方程求解。1905年,Fredholm将积分方程应用于求解弹性力学问题[32]。1953年,Muskhelishvili[33]和Kellogg[34]分别将积分方程法用于求解结构力学和位势问题。1965年,Mikhlin[35]解决了积分方程理论中的奇异性问题,为积分方程法在工程中的应用开辟了道路。
20世纪60年代,电子计算机的发展开创了数值求解积分方程的新时代,积分方程法作为数值计算方法开始应用于实际问题。1963年,Jaswon[36]采用间接边界积分方程法,成功地求解了位势问题和弹性力学问题。这种方法的主要思想是沿边界配置一种虚设的点源密度函数,先确定密度函数,再求边界上的未知物理量。其缺点是待求的点源分布函数是虚构的,不具有明确的物理意义,因此求解过程需要两步完成。但它具有场量方程和场量梯度方程相互独立的优点,因而易于组合求解各种复杂边界条件的边值问题。60年代,一些苏联学者对积分方程尤其是奇异积分方程的理论作了更为深入的研究,为进一步应用边界积分方程方法开辟了道路。与此同时,高速大型计算机的出现及其硬件的迅猛发展使离散求解积分方程成为可能。到了60年代后半期,边界元法的研究更趋活跃,边界元法的直接法被应用于求解各类问题。1967年,Rizzo[37]用直接边界元法求解了二维弹性问题。1969年,Cruse[38]将此法推广到三维弹性力学问题。在直接法中,表述边界积分方程的未知量是真实的物理量,通过求解积分方程可以直接得到边界上所求的未知物理量。1975年,Cruse和Rizzo[39]出版了第一部边界积分方程法的著作。1977年,Banerjee和Butterfield[40]首次采用了boundary element method这一名称。1978年,Brebbia在英国南安普顿召开了第一届国际边界元法会议,出版了专著The Boundary Element Method for Engineers[41]。书中用加权余量法推导出了边界积分方程,并指出加权余量法是最普遍的数值方法,如果以开尔文(Kelvin)解作为权函数,从加权余量法可导出弹性力学问题的边界积分方程,通过将边界离散成单元的方法可数值求解积分方程。至此,边界元法这一名称得到了国际公认。
自1978年第一届国际边界元法会议后,边界元法会议几乎每年在世界各地举办。世界各国已从基本理论与方法的研究向深广领域发展,大量论文和专著先后问世。在此时期,边界元法在数学分析理论和数值积分方法的研究、边界元法的完善和应用范围的拓宽以及边界元法应用软件等方面均得到飞速发展。边界元法的应用遍及固体力学、流体力学、波动学、传热学、电磁学等学科领域。
我国关于边界元法的研究大约开始于1978年,当时杜庆华在国内首先开创了工程中边界元法的研究,开始跟踪国际上这一领域的最新进展,其研究领域主要在固体力学方面[42]。与此同时,王泳嘉[43]、郑颖人等[44]开始了岩土工程中的边界元法研究,杨德全和赵忠生[45]在流体力学边界元法方面开了先河。值得一提的是,岑章志[46]对我国的弹塑性边界元法做了开拓性的研究工作,高效伟和Davies[13]发表了国际上第一个弹塑性边界元程序,结束了三十多年来在非线性力学边界元分析方面只有论文发表、没有程序公布的局面。
近年来,我国学者在快速多极边界元法研究方面取得了一系列的研究成果,代表性工作可见姚振汉与王海涛的著作[11]。此外,我国学者在近奇异积分计算[2428]、多种介质问题[4749]以及与CAD技术结合解决工程问题[50]等方面的研究工作也引人注目。边界元法的研究目前在我国正处于上升阶段[51],近年来在一些回国学者的带领下越来越多的专家学者投身于边界元法的研究中。
1.3边界元法中的难点问题及其研究进展
边界元法在解决接触、断裂力学、运动边界、无限域与半无限域以及超薄结构等问题中具有独特的优势,被大量地用于科学与工程问题的计算分析,在许多应用领域,边界元法被公认为有限元法的一个重要补充。然而,传统的边界元法有下述几方面的弱点:
(1) 奇异性问题。边界元法中所用的基本解是奇异函数,数值计算时首先需要消除积分奇异性才能得到精确的计算结果。多年来已有大量的文献在计算效率与精度方面提出了不少计算奇异积分的方法,如线性单元的解析消除法[10,52]、弱奇异积分的单元子分法[13,53,54]、强奇异积分的间接计算法[13,55]、高阶奇异积分的直接计算法[5658]等。这些方法在计算弱奇异和强奇异积分时非常有效,但在处理高阶奇异积分时的稳定性还需要进一步考查。
(2) 满系数矩阵问题。边界元法在求解问题时所形成的方程组的系数矩阵是满阵,因而占有较大的计算机内存,难以解决大型工程问题。为了解决此问题,研究人员提出了两种有效的解决方法,一种是快速多极法,另一种是区域分解法。快速多极边界元法[11,59],通过将奇异核函数进行级数展开的技术,降低矩阵向量相乘操作的计算量级和存储量级,达到节省内存和提高计算速度的目的。区域分解法[10,60]的基本思想是把总求解域划分成多个子域,对每个子域建立边界元矩阵方程,然后利用子域间公共节点上的面力平衡条件和位移相容性条件把子域方程组集成总体系统方程组。这样形成的系数矩阵是块状稀疏矩阵,可利用现有成熟的稀疏矩阵求解器(如LU分解法和GMRES迭代法(广义最小残量法))对系统方程组进行有效求解。区域分解法是边界元法中被广泛应用的方法,不仅能解决满系数矩阵问题,而且能通过按照材料性质划分子域的手段求解由不同材料组成的复合介质问题[61],还能通过沿裂纹面划分子域的方法求解断裂力学问题[19]。最近,作者基于区域分解法,提出了求解非均匀介质问题的三步变量凝聚技术[49],形成的系统方程组只有公共节点位移作为未知量,可以求解大型工程问题。虽然这些方程组的组建与求解技术能解决相当规模的工程问题,但对于超大型问题,系统方程组的求解仍然是一个极具挑战性的问题。最近,高效伟等提出了求解非对称稀疏方程组的同时消元回代法求解技术[62,63],在计算效率和储存空间方面都有了显著的提高。
(3) 非线性和非均质问题中的域积分问题。传统的边界元法只是在解决线性问题方面比较成熟,对于非线性问题却远非如此。如弹塑性力学问题,虽然从四十多年前就有不少学者开始对该课题进行深入的研究[64,65],但直到2000年,高效伟和Davies[66]才彻底解决了弹塑性边界元法中强奇异域积分的计算问题,并于2002年公开发表了国际上第一个弹塑性边界元程序[13]。我国的其他学者也对弹塑性边界元法的发展做出了重要的贡献[6770]。
用边界元法解决非线性和非均质问题时,由于很难求得控制方程的基本解,所以不得不用对应于线性和均匀介质问题的基本解来建立积分方程,由此导致了域积分的出现。为了计算域积分,传统的方法是将计算域离散成内部网格,采用像有限元法中的方式计算域积分[6470],这样就消除了边界元法只需将边界离散成单元的优点。正是这种求解非线性和非均质问题中需要内部网格的致命弱点,严重影响了边界元法的发展。
为了避免使用内部网格,国际上不少学者致力于无内部网格边界元法的研究,其中最常用的方法是将出现在积分方程中的域积分转换成边界积分。其中应用最广泛的方法是Breb