[发明专利]复杂介质弹性波传播模拟的多网格切比雪夫并行谱元法

说明书

技术领域

本发明涉及地球物理，工程地震学，计算声学等领域的高性能计算。

背景技术

复杂介质中弹性波的传播模拟在地球物理，工程地震学和计算声学等领域具有重要地位。如何找到一种更精确、更快速、更适宜并行化的数值模拟方法也一直是国、内外相关领域研究者的工作目标和重点。

当前在弹性波方程模拟中普遍使用的数值方法，主要有有限差分法、有限元法、伪谱法和谱元法等。它们有各自的优缺点并在不同的领域得到应用。

有限差分法将方程中的导数用差分近似，从而将偏微分方程转换为代数方程来求解，其误差由离散化配置点数目和Taylor级数的截断误差决定。有限差分法由于编程简捷和计算速度较快被广泛应用于计算地球物理和工程地震中。但低阶有限差分法需要大量的网格点来保证精度，而高阶有限差分法又很难高效地处理自由界面和复杂结构问题。

有限元方法基于波动方程的弱形式，将计算区域划分为有限个互不重叠的单元，在每个单元内，波场由低阶多项式(如分段线性函数)来近似。有限元法适合处理边界和内部复杂结构，在结构分析和瞬态模拟方面获得了广泛的应用。但有限元法在大型的弹性波方程模拟领域应用相对较少，其原因一是有限元法需要耗用大量计算资源，二是低阶有限元存在频散现象，而通常的高阶有限元存在伪波。

伪谱法利用无限可微的全局函数（如富里叶Fourier级数），将问题从物理空间转化为波数域进行处理，具有良好的处敛性。伪谱法只需要很少的空间网格点数即可达到很高的精度。为了处理边界条件，在空间展开时经常使用切比雪夫Chebyshev或拉格朗日Legendre正交多项式来取代Fourier级数。但由于正交多项式的配置点分布很不均匀，给时间步长的选择带来很大限制。同时，如同有限差分法一样，伪谱法不能自如地处理弯曲的自由界面和复杂结构的问题。为解决这类问题有人提出了弯曲坐标和区域分解等方法，但代价是计算量的增加。

谱元法由Patera在1984年提出，早期主要应用于流体力学。Seriani等在1991年首次将Chebychev谱元法引入到声波方程的数值模拟中，各国研究者在其后做了大量的相关研究。它把有限元法和伪谱法相结合，兼具了有限元处理边界和结构的灵活性和伪谱法的快速收敛特性。在达到相同的精度前提下，可以采用比传统有限元更稀疏的单元划分，减少了计算时间和内存需求。

谱元法的基本方法是在每一个单元上用高阶谱展开。选取以截断的正交多项式表示的基函数，在各个单元上利用配置点插值，以提高解的收敛速度。其主要步骤是：(1)首先把计算的区域分成许多子域（单元），每个子域由若干节点（配置点）组成;(2)在每个子域中把近似解表示成截断的正交多项式展开;(3)用Galerkin方法求解正交问题的变分格式，得到全局的近似解。

谱元法可以用较稀疏的网格和单元获得较高精度，但传统的谱元法中，每个单元内只能有单一均匀介质，某些情况下严重降低计算效率。比如当介质结构复杂，变化尺度很小甚至小于波长时，必须按介质变化的小尺度划分单元和求解，这样将会造成极大的计算规模。

此外，传统的谱元法求解波动方程时，需要形成全局刚度矩阵和质量矩阵，考虑衰减时还需引进阻尼矩阵，需要耗费大量的存贮空间，在算法上限制了其并行化的效率。注意到在波动方程的时间迭代过程中，需要做多次矩阵和向量乘积，如Sp_k，这里S和p_k分别为全局刚度（或质量）矩阵和全局向量。考虑到全局矩阵的稀疏特性，全局的矩阵和相量相乘耗费了大量的计算时间。

以上两点都极大地阻碍了谱元法在地球物理，工程地震学，计算声学等领域的大规模实际应用。

发明内容

本发明的目的在于提供一种能够克服上述缺点的算法。

为此，本发明提供一种多网格Chebyshev并行谱元方法。该方法包括：

a.将计算区域划分为大而规则的单元；

b.在单元内定义主网格和辅助网格，主网格为Chebyshev配置点；

c.利用主网格将单元内波场作Chebyshev截断展开，利用辅助网格将单元内介质参数和外力作截断展开；

d.将波场、介质参数和外力的截断展开式代入波动方程，得到单元质量矩阵、单元刚度矩阵和单元外力向量；

e.基于各单元的单元质量矩阵、单元刚度矩阵和单元外力向量，在主网格上求解弹性波方程。

优选地，利用预条件共轭梯度法求解弹性波方程。