边界元法的若干进展和它在固体力学中的应用

1、边界元法的若干进展和它在固体力学中的应用,清华大学工程力学系 姚振汉,引 言,弹性力学的三种提法,微分提法,变分提法,积分提法,偏微分方程边值问题,泛函极值问题,边界积分方程问题,求解析解,差分法求数值解,李兹法求近似解,有限元法求数值解,边界元法求数值解,归结为求解线性代数方程组,常化为常微分方程,三种提法是完全等价的,边界元发展历史回顾,边界积分方程边界元法 有限元法（1955，56）之后发展起来的一种精确高效的工程数值分析方法 在固体力学领域有限元法最重要的补充 边界元法间接法 位势问题（Smith & Pierce, 1958) 弹性力学（Massonet，1965） 边界元法直接法

2、位势问题（Jaswon，1963） 弹性力学（Rizzo，1967）,边界元发展历史回顾,19942003被SCI收录的论文 与边界元法有关的有3904篇 与有限元法有关的为16823篇 19902002被EI收录的论文 与边界元法有关的有19968篇 与有限元法有关的为75184篇 与断裂力学有关的为23647篇 在工程应用方面 在应用最多的部门也从未超过有限元法的十分之一,研究组边界元研究历史回顾,我们研究组边界元法研究开始于1979年 基于弹性力学问题Rizzo型边界积分方程边界元法研究了弹性应力集中问题和薄板弯曲问题 研究了边界元有限元耦合方法 研究了边界元法在形状优化缺陷识别等逆问题

3、中的应用 Galerkin对称边界元法用于结构极限与安定分析等问题 研究了精确高效的计算方法，提出了边界元法误差的一种直接估计。 对于弹性接触问题，提出了单元与单元间协调的接触方案，研究了二维、三维移动、滚动接触。,研究组边界元研究近年工作,2000前后针对复合材料，对于含随机分布大量夹杂的二维弹性固体提出了一种重复相似子域边界元法，计算了100多个夹杂、近万自由度的问题。 研究了在微机机群上的并行算法。2000年在由8台微机组成的机群上最大计算规模45,000自由度。 近年来多极快速算法在边界元法中的应用给边界元法解决复杂工程与科学问题展示了广阔的前景。 用于含随机分布夹杂二维、三维弹性体数

4、值模拟，一台微机可计算数十万自由度的问题，在微机机群并行系统最大的二维算例有800万自由度。 计算了含16384条随机分布裂纹的二维无限弹性体，1,572,864自由度，研究了相应的裂纹扩展问题。,鉴于 边界元法始终是计算力学中一个重要的研究领域 快速多极边界元法是近年受到特别关注的一个研究方向 本人已经从事边界元法研究25年，几乎没有间断 近年研究组研究重点放在快速多极边界元法的研究方面 因此 今天就以此为题，希望能对各位同学的学习有一定的帮助。,弹性力学边界积分方程,弹性静力学的边界积分方程 其中,对于同一弹性体的两种变形状态,通常一个状态是待求状态 另一状态是已知的辅助状态,弹性力学的边

5、界积分方程可以由Betti定理出发导出,以无限弹性空间中一点作用单位集中力的Kelvin问题的解为辅助解,由Betti定理可得Somigliana等式（无体积力情况）,再将集中力作用点（源点）P 趋于边界点 p，考虑到源点是积分的奇异点，作适当处理，即可得到前面给出的边界积分方程。,弹性力学边界元法,将边界分成边界单元 在每个单元上将边界变量插值离散 可得,将弱解代入方程得到误差 对于加权余量法的配点格式，权函数采用Dirac-delta函数，要求 即得,弹性力学边界元法,将方程写成矩阵形式，得 其中H，G矩阵元素由核函数与形函数在单元上的积分求得。 将边界变量列矢量 U，T 按未知量与给定量

6、重新排列，可得边界元法的求解代数方程组 即,弹性力学边界元法,核函数与形函数乘积在单元上的积分 矩阵元素都是核函数与形函数乘积在单元上的积分，矩阵是满阵。 主要计算量就是计算这些积分，以及求解满阵代数方程组。 对于规模不太大的问题，计算积分的工作量是主要计算量。 非奇异积分采用等精度高斯积分格式求积，高斯点数由精度要求对不同情况自动确定。 当核函数的源点落入积分单元时出现奇异积分，包括弱奇异积分和柯西主值积分。,弹性力学边界元法,边界元法的优缺点,边界元法与有限元法及其它数值方法相比较的优缺点 优点：高精度（由于采用了解析基本解） 降低了维数，便于模拟复杂边界形状 对于高梯度、甚至有奇异性的问

7、题，不仅有较高的精度，而且在同等精度条件下有较高效率 适合于处理无限域、半无限域问题 适合于处理弹性接触等边界条件非线性问题 缺点：适用范围远没有有限元法广泛 解题规模受限制（方程系数矩阵为满阵） 对于域内方程非线性问题优势减弱或丧失,边界元法的特点,早期有的文章通过简单算例高估了边界元法的计算精度，声称能达到0.01%，甚至更高。 这里要分两类问题，一类是没有离散误差的问题，称为简单问题，离散插值的边界变量能精确满足给定边界条件，此时能达到很高精度，例如106，甚至108。 对于一般问题，有离散误差，只有合理划分足够多的边界单元，才能达到要求的计算精度。 薄板梁纯弯曲问题，只要采用二次单元，

8、就是一个简单问题。 通常认为细长薄板梁要划分较多单元才能达到满意精度，我们只用4个二次元，就达104精度。,有限元法中单元边长比不能太大，相邻单元尺寸也不能相差太大。 边界元法则并不受此限制，上述100:1纯弯薄板梁相邻单元长度比为100，只要保证积分精度等运算精度，还是可以得到高精度的结果。 有限元法中高斯积分通常只用12，或22个高斯积分点。 边界元法常采用等精度高斯积分，根据给定积分精度要求来确定用多少高斯点。对于上述细长薄板梁的算例，源点在边长为1的短边中点、在长边二次单元上的积分，需要5060个高斯点进行高斯积分。,边界元法的特点,边界元法在解出未知的边界变量之后，通过域内变量的边界

9、积分公式可以求得域内任意点的位移、应变和应力。 这样得到的域内位移场精确满足Navier方程，但是当它趋于边界时的极限一般说来和边界给定量或边界元法解出的边界变量并不一致。这种差别是和离散误差密切相关的。对于没有离散误差的简单问题，域内解趋于边界的极限和边界变量是一致的，只有很小的运算误差。 这也导致在有离散误差的情况下，在近边界区求得的内点变量有较大的误差，即通常称为“边界效应”。,边界元法的特点,边界元法模拟光弹实验,为了便于将得到的结果和工程界熟悉的光弹性实验结果比较，域内应力分布也画出了等色线和等倾线图。 而且我们一再表明，在数值模拟已经很成熟的今天，进行光弹性实验这样的物理模拟已经没

10、有多大意义。因为它也不是实测，光弹性模型和实际问题相比引进了许多误差，得到的结果并不比数值模拟结果更接近真实情况；而且数值模拟比光弹性模拟更经济便捷。 下图是带孔板单向拉伸的等色线图（条纹等级反映主应力差）。,边界元法模拟光弹实验,边界元法并行计算的工程应用,边界元法结合并行计算，在8台微机组成的机群并行系统上计算了一些规模稍大的工程实际问题。例如： 石油钻杆的偏心钻挺，49818自由度，8台微机，计算43小时41分钟。 渤海石油公司“滨海109”海洋铺管船的锚机滚筒轴，24444自由度，7台微机，计算436分钟。 这是2000年博士生尹欣得到的结果，当时文献中最大计算规模10万自由度。,子域

11、 1的网格划分,子域 25的网格划分,偏心钻挺模型（5万自由度）,有限元模型（16000结点，13000千单元）,钻挺顶部Mises应力图（边界元法，单位：MPa）,钻挺顶部Mises应力图（边界元法，单位：MPa）,钻挺顶部Mises应力图（有限元法，单位：MPa）,钻挺顶部Mises应力图（有限元法，单位：MPa）,滚筒轴边界元模型：5个子域，8184结点，1721单元,滚筒轴有限元模型：13000结点，64000单元,边界元解：左端局部Mises应力（单位：MPa）,有限元解：左端局部Mises应力（单位：MPa）,边界元解：右端局部Mises应力（单位：MPa）,有限元解：右端局部Mi

12、ses应力（单位：MPa）,在复合材料模拟中的应用,针对复合材料数值模拟，研究组对于含随机分布孔洞或夹杂的二维弹性体进行了深入的研究。提出了重复相似子域边界元法，提高了计算效率。分析计算了含100个圆孔、圆形夹杂、椭圆形夹杂、多种形状夹杂的二维弹性体，得到了等效材料常数和体积比、夹杂与基底材料常数比等的关系。考虑了夹杂和基底之间有界面层的情况。 对此类问题，引进快速算法之后，使计算规模大大提高，初步工作已经从100多个夹杂增加到了1600个夹杂。,在复合材料模拟中的应用,重复相似子域边界元法,在复合材料模拟中的应用,含100个随机分布圆形夹杂的二维弹性体,在复合材料模拟中的应用,含100个随机

13、分布圆形夹杂的二维弹性体变形图,在复合材料模拟中的应用,含100个随机分布圆形夹杂的二维弹性体应力分布,在复合材料模拟中的应用,含100个随机分布椭圆形夹杂的二维弹性体,在复合材料模拟中的应用,含100个随机分布椭圆形夹杂二维弹性体的变形图,在复合材料模拟中的应用,含100个随机分布椭圆形夹杂二维弹性体的应力分布,在复合材料模拟中的应用,含100个随机分布不同形状夹杂的二维弹性体,在复合材料模拟中的应用,含100个随机分布不同形状夹杂二维弹性体的变形图,在复合材料模拟中的应用,含100个随机分布不同形状夹杂二维弹性体的应力分布,在复合材料模拟中的应用,含100个随机分布圆形夹杂（有界面层）的二

14、维弹性体,在复合材料模拟中的应用,含圆形夹杂二维弹性体等效弹性模量和模量比的关系,在复合材料模拟中的应用,含圆形夹杂二维弹性体等效体积模量和模量比的关系,在复合材料模拟中的应用,含圆形夹杂二维弹性体等效剪切模量和模量比的关系,在复合材料模拟中的应用,含圆形夹杂二维弹性体等效弹性模量和模量比的关系,快速多极算法,快速多极算法 (fast multipole method, FMM) 最初针对问题：多粒子相互作用的有势场 存储和计算复杂度：O(N2) Barnes, Hut, (1986) 树结构算法 存储和计算复杂度：O(NlogN) Greengard, Rokhlin, (1986) 快速多

15、极算法 存储和计算复杂度：O(N) Hrycak T, Rokhlin V, (1998) 新型快速多极算法 存储和计算复杂度：O(N),快速多极边界元法,常规边界元法由于方程组的系数矩阵为非对称满阵，对于N自由度系统，其运算量为N 3量级，对于存储的要求为N 2量级。 边界元快速算法通过快速多极展开，使运算量和存储量均减少到NlogN 量级。最近引入了局部展开的思想，进一步使运算量和存储量要求减少到了N 量级。 研究组初步的工作已经把解题规模明显扩大，由原来8台微机并行最多计算49,818自由度问题，到用边界元快速算法一台微机计算544,000自由度和1,572,864自由度。,边界元快速算

16、法和常规边界元法的比较,快速多极边界元法基本思想,快速多极算法的基本思想： 如果每个人要自己把信送到每个收信人，是非常费时的，当收信人的范围很大时是不可能的。 但是现代的邮政系统，使每个人可以方便地给世界各地的收信人发信。 传统边界元法形成满阵方程组的系数要一个个地独立计算，就像每个人亲自送信一样。 快速多极算法就像现代邮政系统一样大大提高了效率。,树型存储结构,将求解区域与边界有关部分划分为四叉树（二维问题）或八叉树（三维问题）,近场与远场,对源点 x，近场按常规边界元法计算方程系数，远场用快速多极展开，不再计算一个个系数。,快速多极算法基本步骤,对于远场的快速多极算法有四个基本步骤：,核函

17、数的多极展开 多极展开系数（多极矩）g (y, k) 对所有核函数源点只要计算一次。,核函数的多极展开（平面问题常值元为例）,基层的邮局可以方便地将任何地方转来的邮件送给管辖区域任何住户,展成复数泰勒级数,多极展开系数的转移（M2M）,上级邮局利用下级邮局可以方便地将任何地方转来的邮件送给管辖区域任何住户,多极展开系数向局部展开系数的转移（M2L）,发信方上级邮局利用收信方邮局可以将任何地方转来的邮件送达收信的任何住户,局部展开系数的转移（L2L）,发信方将邮件交给最近的基层邮局就可以利用邮局系统送达任何收信人,快速多极边界元法采用迭代解法,常规边界元法解题规模较小，自由度数通常小于一万，常用

18、Gauss消去法解非对称满阵方程组。,快速多极边界元法，解题规模达到数百万、上千万，并不形成满阵方程组，不存储方程组系数，一定配合使用迭代解法。,为了保证迭代尽快收敛，采用适当的预条件处理技术也十分重要。,快速多极边界元法考题,计算考例（和边界元法比较）,多极快速边界元法考题,快速多极边界元法算例,含1600个随机分别圆形夹杂的二维弹性体,Eb =200.0MPa,快速多极边界元法算例,夹杂体积比、模量比对等效体积模量的影响,快速多极边界元法算例,快速算法老方案和常规边界元法计算时间比较,快速多极边界元法算例,快速算法老方案和常规边界元法存储需求比较,快速多极边界元法算例,快速算法新方案和老方案计算时间比较,快速多极边界元法算例,含300个随机分布球形夹杂的三维弹性体,快速多极边界元法算例,300个随机分布球形夹杂的界面正应力,快速多极边界元法算例,300个随机分布球形夹杂的界面剪应力,快速多极边界元法算例,等效弹性模量与夹杂模量比的关系,快速多极边界元法算例,等效体积模量与夹杂模量比的关系,快速多极边界元法算例,等效弹性模量与夹杂模量比的关系,快速