基于Python的二维有限元声波方程正演计算

基于Python的二维有限元声波方程正演计算

一、基础理论与相关公式的导出

1. 什么是有限元方法？

• 有限元是计算复杂数学问题近似解的工具。当数学方程过于复杂，无法用正常的方法求解，并且一定程度的误差是可以容忍的时候，通常使用这种方法。工程师们通常使用有限元法，因为他们关心的是为实际应用设计产品，而不需要近乎完美的解。有限元法可以适应不同的精度要求，可以减少设计过程中对物理原型的需要。

2. 有限元方法的优点：

1. 可以说，有限元相比于其他的数值计算方法，它对复杂几何图形或非规则形状的建模要简单的多，这或许就是有限最大的优点，任意复杂的区域可以被有限个单元所离散。
2. 优良的可移植性，不同形式的材料特性可以被轻松应用到模型中去。
3. 有限元单元形式多种多样，如果需要可以使用高阶的有限单元，同时有限元方法可以处理多种物理实际问题，譬如静力学、动力学等问题都可以进行建模。
4. 有限元方法是相较简单的，它几乎适定任何问题，在实际工程中应用的最为广泛，受广大工程师的喜爱。
5. 人们为有限元开发了大量软件和相应算法，使得有限元方法对于新手来说都是具有强大的解决问题的能力。

3. 有限元方法的缺点：

1. 跟其他数值方法一样，有限元方法逼近是一种近似计算方法，存在不可避免的误差。
2. 有限元方法的基础理论简单，但是已经很多年没有在这方面有所突破。
3. 有限元的计算量可以说是非常巨大，本人算一个128x128大小的都把计算机算到重启！同时其需要庞大的内存来储存如刚度矩阵这样的大型矩阵。
4. 有限元计算的结果还需工程师进行后续的判断。

4. 二维声波方程的伽辽金方法（Galerkin Method）

• 伽辽金方法采用微分方程对应的弱形式，其原理为通过选取有限多项式函数（又称基函数或形函数），将它们叠加，再要求结果在求解域内及边界上的加权积分（权函数为试函数本身）满足原方程，便可以得到一组易于求解的线性代数方程，且自然边界条件能够自动满足。
  写出我们要求解的方程如下：

$$\begin{array}{r}{\partial }_t^{2}u=v^2{\nabla }^{2}u+f\end{array}$$其中，$\nabla =\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$为二维拉普拉斯算子，$f$为震源扰动项。

• 伽辽金方法简单来说分为几步，一是将解用插值函数的和来表示，二是偏微分方程的弱形式，三是将空间项的偏微分进行分部积分：

1. 解用二维分段线性插值函数来表示。

所谓的分段线性插值函数就是在某个节点处值为1，其余处值均为零的函数，并且值为1的节点与附近0值节点之间为线性关系。假设这个二维分段线性插值函数为基函数${\phi }_i(x,y)$，则我们可以得到插值后的近似解为：
$$\begin{array}{r}u(x,y,t)=\sum _{k=1}^N{u}_k(t){\phi }_k(x,y)\end{array}$$其中，$N$ 表示节点数，对于一个局部的三角形单元 $N=3$。
2. 微分方程的弱形式
即加权积分满足原方程：
$$\begin{array}{r}{\int }_{\Omega }\left[{\nabla }_t^{2}u{\phi }_i(x,y)-v^2{\nabla }^{2}u{\phi }_i(x,y)\right]d\Omega ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega \end{array}$$
对上式关于空间项的积分使用分部积分法则:
$$\begin{array}{r}{\int }_{\Omega }{\nabla }^{2}u{\phi }_i(x,y)d\Omega ={\nabla u{\phi }_i(x,y)|}_{l_1}^{l_2}-{\int }_{\Omega }\nabla u\nabla {\phi }_i(x,y)d\Omega \end{array}$$
将上式代入得到:
$$\begin{array}{r}{\int }_{\Omega }\left[{\nabla }_t^{2}u{\phi }_i(x,y)+v^2\nabla u\nabla {\phi }_i(x,y)\right]d\Omega ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega +{\nabla u{\phi }_i(x,y)|}_{l_1}^{l_2}\end{array}$$如果令上式中的${\nabla u{\phi }_i(x,y)|}_{l_1}^{l_2}$等于零的话，那么就正好满足了自由边界条件。
最终，我们得到的公式为：
$$\begin{array}{r}\sum _{k=1}^N{\nabla }_t^2{u}_k{\int }_{\Omega }{\phi }_k(x,y){\phi }_i(x,y)d\Omega +v^2\sum _{k=1}^N{u}_k{\int }_{\Omega }\nabla {\phi }_k(x,y)\nabla {\phi }_i(x,y)d\Omega ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega \end{array}$$
时间项我们采用二阶差分格式，则上式的隐式差分格式（所谓隐式差分就是上述方程的第二项采用了未知时刻的波场或者说下一时刻的波场值$u_k^{n+1}$）为:
$$\begin{array}{r}\sum _{k=1}^N\frac{u_k^{n+1}-2u_k^n+u_k^{n-1}}{\Delta t^2}{\int }_{\Omega }{\phi }_k(x,y){\phi }_i(x,y)d\Omega +v^2\sum _{k=1}^N{u}_k^{n+1}{\int }_{\Omega }\nabla {\phi }_k(x,y)\nabla {\phi }_i(x,y)d\Omega ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega \end{array}$$
简写成矩阵形式:
$$\begin{array}{r}\left(M+\Delta t^2{v}^{2}S\right)U^{n+1}=\left(2U^n-U^{n-1}\right)M+\Delta t^{2}F\end{array}$$
同样的，如果使用当前时刻的波场值（显式差分），则矩阵形式的显式差分格式为:
$$\begin{array}{r}{U}^{n+1}=2U^n-U^{n-1}+\left[M^{-1}\left(F-v^{2}S{U}^n\right)\right]\Delta t^2\end{array}$$由于历史原因，有限元开发最初的目的就是为了解决刚体或者说是结构力学中的受力问题，这里的$ M $称之为质量矩阵（mass matrix）， $S$ 称之为刚度矩阵（stiffness matrix），它们的计算形式为：
M = ∑ k = 1 N ∫ Ω φ k ( x , y ) φ i ( x , y ) d Ω S = ∑ k = 1 N ∫ Ω ∇ φ k ( x , y ) ∇ φ i ( x , y ) d Ω F = ∫ Ω f φ i ( x , y ) d Ω \begin{aligned} M &= \sum_{k=1}^N\int_{\Omega} \varphi_k(x, y) \varphi_i(x, y) d \Omega\\ S&=\sum_{k=1}^N\int_{\Omega} \nabla \varphi_k(x, y) \nabla \varphi_i(x, y) d \Omega\\ F&=\int_{\Omega} f \varphi_i(x, y) d \Omega \end{aligned} MSF​=k=1∑N​∫Ω​φk​(x,y)<