地震波方程时间域边界元法数值模拟

【摘要】：边界单元法是求解数学物理方程的一种数值计算方法。这种方法是把所研究的问题的微分方程变成边界积分方程，然后将区域的边界划分为有限个单元，也就是把边界积分方程离散化，得到只含有边界上的结点未知量的方程组，然后进行求解。 这种方法在处理问题时与目前广泛使用的有限单元法和有限差分法有某些类似之处，但其出发点完全不同。 有限单元法和有限差分法属于所谓“区域法”。这些方法的出发点是把问题的连续“区域”划分成许多细小的单元或网格，然后把各单元或网格换成简单的等价模型，再把它们联系起来进行全部计算。也就是说，把原来的分布参数系统问题化为集中参数系统问题来求解。这些方法的基本思想是用完全（或局部）满足定义域上边界条件的函数去逼近问题控制微分方程。边界单元法与它们正好相反，它是把定义域的边界划分成一系列单元，用满足控制微分方程的函数去逼近边界条件，在单元上所考虑的函数可以按不同的形式变化，这一做法与有限单元法大致相同。 边界单元法可以分成两种基本类型，即直接法和间接法。间接法是从一个基本解入手，该解在定义域内满足控制微分方程，但却含有某些未知数，这些未知数则通过在许多点（或子域）上施行边界条件来确定。间接边界单元法是用物理意义不一定很明确的变量来表示化成的公式。这种方法曾用于求解由拉普拉斯方程或亥姆霍茨方程所控制的弹性力学问题或其它势问题。直接边界单元法是以格林恒等式为出发点，变量具有明确的物理意义，现在已优先应用在工程科学中，本论文只讨论这种方法。 边界元法是将区域上的控制方程转化为沿区域边界的积分方程，因此它只需要定义边界上的单元，结合边界条件求解，这样就使处理问题的维数降 WP=59 低一维，这一步的关键就是格林公式和基本解。由于这种方法只将边界离散来求解，求解一个问题所建立的方程组阶数低，数据大为减少，因而计算机的内存要求也降低。 但是，边界元法所建立的方程组的系数矩阵是稠密的，一般是非对称的，而且矩阵元素分量的计算量很大，这就抵消了降阶之后矩阵消元所能省出的一部分时间。此外，边界元法在处理非均一、非同质的问题时不够灵活。 边界元法和有限元法都要以一个试探解函数来实现对求解区域的离散，这两种方法在选择这些函数时采用了不同的准则。用边界元法求解边值问题需要找到控制微分方程的一个基本解或控制微分方程在求解区域上的格林函数，这对于某些问题是十分困难的。有限元法的主要缺点之一是不大适合求解无限边界场域边值问题，而只能求解有界问题，因为用有限个单元离散无限域显然是不可能的。用有限元法难于处理的另一类问题是域内具有应力奇异的问题。在固体力学问题中，这类应力奇异通常发生在不规则的凹角或孔洞附近。由于应力奇异可能引起断裂扩展，因此在奇异点附近能否得到一个较为精确的解答，有时就显得十分重要。但是在奇异点处，理论上应力为无限大，用有限元法可能产生毫无意义的分析结果。除了对于狭长形状的解域处，边界元法的求解精度一般高于有限元法。边界元法对两种不同类型变量（即在一些场问题中的势和通量，或在应力分析中的位移和应力）一般都能获得较好的结果。而有限元法一般只是对于控制微分方程中所考虑的变量（即势或位移）的求解精度比较高，但对于这些变量的导数（通量或者是应力）所得结果就不很精确，而且通常使得它们在单元之间成为不连续。如果在媒质中存在高通量或应力集中区域，这个问题将更加突出。而边界元法能够解决以上的问题。边界元法在得出边界近似解之后，虽得不到解析显式，但可以逐点计算域内点的近似解，而有限元法则必须同时对所有域内结点联立求解。因此，当只需对个别点求解时，边界元法较简便。 本文阐述了边界元法研究的国内外现状，介绍了 函数与Green函数的 WP=60 概念与性质，列举了各类微分方程的基本解，研究了二维波动方程时间域边界元解法的数学理论与离散化计算公式，研制了Fortran语言计算机程序。最后，我们用二维波动方程时间域边界元法及其相应的计算机程序进行了地震波场的数值模拟，构造了三个理论地震模型，进行了大量的模拟计算和分析对比，证实了这种模拟方法的有效性，得到如下一些有益的结论： 利用时间域边界元法对地震波场进行数值模拟，其结果是可靠有效的。对于复杂地表情形的地球模型，用时间域边界元法进行数值模拟也是方便有效的，而且地表对模拟结果的影响，可以得到反映。利用第一类边界 作为人工截断边界条件，方法容易实现，但要将计算范围扩大，以减少边界的影响，这也会带来使计算量增大的问题。利用时间域边界元法模拟分块均匀介质的地球模型是可行的，但当分块均匀介质过多时，边界法的优势与有限元法相比就会变得不够明显。利用第三类边界 作为人工截断边界，不必扩大计算区域，可减少计算量，但要准确给出变量 的值还是很困难的，因此计算效果不佳。第三类边界是一种近似的吸收边界，寻求更好的吸收边界条件及相应的Green函数是改进此方法的一个方向。