基于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），它们的计算形式为：
$$\begin{array}{rl}M& =\sum _{k=1}^N{\int }_{\Omega }{\phi }_k(x,y){\phi }_i(x,y)d\Omega \\ S& =\sum _{k=1}^N{\int }_{\Omega }\nabla {\phi }_k(x,y)\nabla {\phi }_i(x,y)d\Omega \\ F& ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega \end{array}$$

二、网格划分

本文使用的是基于矩形网格划分处的三角单元网格，它们具有结构化网格的特性，便于索引。主要有如下几种形式：

右斜对角的三角网格

左斜对角的三角网格

全对角的三角网格

在上述的三种网格图中，黑色数字字样表示的是节点索引，蓝色数字字样表示单元索引，对于不同网格在计算雅克比矩阵时要注意选取基函数的顺序，谨慎对照参考单元的节点，否则可能会出现假的数值各项异性的现象。

三、数值实验

3.1 求解质量矩阵

首先写出质量矩阵的表达式：
$$\begin{array}{rl}M& ={\int }_{\Omega }{\phi }_k(x,y){\phi }_i(x,y)d\Omega \end{array}$$可以看出质量矩阵其实就是两个基函数乘积的面积分，以全对角的三角网格为例，其基函数如下图所示：

全对角的三角网格的基函数 看起来质量矩阵的求解是要对不同节点的基函数进行求积之后再积分，但是对于两个不相邻的基函数，它们的乘积为零；对于两个有非零重叠部分的基函数，它们乘积的积分也可以划分成若干个小三角单元的面积分，所以最终求取质量矩阵的过程变成了有限个小单元的面积分，通过求取小单元的局部质量矩阵，将它们组装起来得到最终大的质量矩阵。

标准三角形的坐标映射

对于不同形状的三角单元，它们的面积分求取结果不同，通过坐标映射之后，得到以下的问题：
从一个坐标系（x，y）中的已知三点坐标$(x_{i,j,k},y_{i,j,k})$，映射到另一个坐标系（u，v）中的三点坐标$[(1,0),(0,1),(0,0)]$的问题，求解其中的映射关系，讨论不同坐标下的一个微分问题。

简单来说，就是求解一个坐标变换矩阵$M$，将坐标系（x，y）中的三点坐标矩阵 $A$ 转换到坐标系（u，v） $B$ 中去:
$$\begin{array}{rl}B& =M\times A\\ A=\left[\begin{array}{ccc}{x}_i& x_j& x_k\\ y_i& y_j& y_k\\ 1& 1& 1\end{array}\right]& ,B=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 1& 1& 1\end{array}\right]\\ M& =B\cdot A^{-1}\\ & =\left[\begin{array}{ccc}{x}_i-x_k& x_j-x_k& x_k\\ y_i-y_k& y_j-y_k& y_k\\ 0& 0& 1\end{array}\right]\\ x& =(x_i-x_k)u+(x_j-x_k)v+x_k\\ y& =(y_i-y_k)u+(y_j-y_k)v+y_k\end{array}$$
则它们之间的微分关系为：
$$\begin{array}{r}\left|det(J)\right|=\left|\begin{array}{cc}\frac{dx}{du}& \frac{dx}{dv}\\ \frac{dy}{du}& \frac{dy}{dv}\end{array}\right|=\left|\begin{array}{cc}{x}_i-x_k& x_j-x_k\\ y_i-y_k& y_j-y_k\end{array}\right|\end{array}$$

同时他们之间的面积微元关系为：
$$\begin{array}{rl}dx\times dy& =\left(\frac{dx}{du}du+\frac{dx}{dv}dv\right)\times \left(\frac{dy}{du}du+\frac{dy}{dv}dv\right)\\ & =\frac{dx}{dv}dv\times \frac{dy}{du}du+\frac{dx}{du}du\times \frac{dy}{dv}dv\\ & =\left(\frac{dx}{du}\frac{dy}{dv}-\frac{dx}{dv}\frac{dy}{du}\right)du\times dv\\ & =\left|\begin{array}{ll}\frac{dx}{du}& \frac{dx}{dv}\\ \frac{dy}{du}& \frac{dy}{dv}\end{array}\right|du\times dv\\ & =|\det (J)|du\times dv\end{array}$$
对于标准的三角单元，其质量矩阵中的基函数在该面积上的积分为：

参考标准三角形

$$\begin{array}{r}{\int }_0^1{\int }_0^{1-x}{\phi }_i(x,y){\phi }_k(x,y)dxdy\end{array}$$
基函数在节点上为分段线性插值函数，这些基函数这个标准三角形区域（这里应该为uv坐标，但是为了计算写成了xy坐标）对应的表达式为
$$\begin{array}{rl}& {\phi }_1(x,y)=x\\ & {\phi }_2(x,y)=y\\ & {\phi }_3(x,y)=1-x-y\end{array}$$
对于$i$等于$k$的情况下，该面积积分有
$$\begin{array}{rl}1.i=k=1& {\int }_0^1{\int }_0^{1-x}x\cdot xdxdy\\ & ={\int }_0^1(1-x)x^{2}dx\\ & =(\frac{1}{3}{x}^2-\frac{1}{4}{x}^4){|}_0^1\\ & =\frac{1}{12}\\ 2.i=k=2& {\int }_0^1{\int }_0^{1-x}y\cdot ydxdy\\ & ={\int }_0^1\frac{1}{3}(1-x{)}^{3}dx\\ & =\left[-\frac{1}{12}(1-x{)}^4\right]{|}_0^1\\ & =\frac{1}{12}\\ 3.i=k=3& {\int }_0^1{\int }_0^{1-x}(1-x-y{)}^{2}dxdy\\ & ={\int }_0^1\frac{1}{3}(1-x{)}^{3}dx\\ & =\left[-\frac{1}{12}(1-x{)}^4\right]{|}_0^1\\ & =\frac{1}{12}\end{array}$$

当$i$不等于$k$时，有，

$$\begin{array}{rl}1.i=1,k=2& {\int }_0^1{\int }_0^{1-x}x\cdot ydxdy\\ & ={\int }_0^1\frac{1}{2}(x^3-2x^2+x)dx\\ & =\frac{1}{2}(\frac{1}{4}{x}^4-\frac{2}{3}{x}^3+\frac{1}{2}{x}^2){|}_0^1\\ & =\frac{1}{24}\\ 2.i=1,k=3& {\int }_0^1{\int }_0^{1-x}x\cdot (1-x-y)dxdy\\ & ={\int }_0^1\frac{1}{2}(x^3-2x^2+x)dx\\ & =\frac{1}{2}(\frac{1}{4}{x}^4-\frac{2}{3}{x}^3+\frac{1}{2}{x}^2){|}_0^1\\ & =\frac{1}{24}\\ 3.i=2,k=3& {\int }_0^1{\int }_0^{1-x}y\cdot (1-x-y)dxdy\\ & ={\int }_0^1\frac{1}{6}(1-x{)}^{3}dx\\ & =\frac{1}{6}(-\frac{1}{4}(1-x{)}^4){|}_0^1\\ & =\frac{1}{24}\end{array}$$
综上所述，可得在参考的标准单元上的质量矩阵为：
$${\int }_{\Omega }{\phi }_i(x,y){\phi }_k(x,y)dxdy=\{\begin{array}{l}\frac{1}{12},当i=k时\\ \frac{1}{24},当i\ne k时\end{array}$$最终得到质量矩阵的计算式为：
$$\begin{array}{rl}M& ={\int }_{\Omega }{\phi }_k(x,y){\phi }_i(x,y)d\Omega \\ & ={\int }_0^1{\int }_0^{1-u}{\phi }_k(u,v){\phi }_i(u,v)dudv\times \left|det(J)\right|\end{array}$$

3.2 刚度矩阵的计算

同样的，写出刚度矩阵的计算式：
$$\begin{array}{rl}S& ={\int }_{\Omega }\nabla {\phi }_k(x,y)\nabla {\phi }_i(x,y)d\Omega \end{array}$$跟质量矩阵一样，刚度矩阵也是分单元计算的，只不过这里计算的是基函数导数乘积的面积分，同时在经过坐标映射之后，基函数求导也会发生变化，利用微分的链式法则，我们得到：
$$\begin{array}{l}\frac{\partial \phi }{\partial x}=\frac{\partial \phi }{\partial u}\frac{\partial u}{\partial x}+\frac{\partial \phi }{\partial v}\frac{\partial v}{\partial x}\\ \frac{\partial \phi }{\partial y}=\frac{\partial \phi }{\partial u}\frac{\partial u}{\partial y}+\frac{\partial \phi }{\partial v}\frac{\partial v}{\partial y}\end{array}\}$$
$$\left(\frac{\partial \phi }{\partial x},\frac{\partial \phi }{\partial y}\right)=\left(\frac{\partial \phi }{\partial u},\frac{\partial \phi }{\partial v}\right)\left(\begin{array}{ll}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right)$$
其中$\begin{array}{r}{J}^{-1}=\left(\begin{array}{ll}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right),\nabla \phi (x,y)=\nabla \phi (u,v)J^{-1}\end{array}$，由标准三角形的基函数可知：
$$\begin{array}{rlrl}& {\phi }_1(u,v)=u& \nabla {\phi }_1(u,v)& =(1,0)\\ & {\phi }_2(u,v)=v& \nabla {\phi }_2(u,v)& =(0,1)\\ & {\phi }_3(u,v)=1-u-v& \nabla {\phi }_3(u,v)& =(-1,-1)\end{array}$$最终，得到了刚度矩阵的计算式：
$$\begin{array}{rl}S& ={\int }_{\Omega }\nabla {\phi }_k(x,y)\nabla {\phi }_i(x,y)d\Omega \\ & ={\int }_0^1{\int }_0^{1-u}\left[\nabla {\phi }_k(u,v)J^{-1}\right]{\left[\nabla {\phi }_i(u,v)J^{-1}\right]}^{T}dudv\times \left|det(J)\right|\\ & =\frac{1}{2}\left[\nabla {\phi }_k(u,v)J^{-1}\right]{\left[\nabla {\phi }_i(u,v)J^{-1}\right]}^T\left|det(J)\right|({\int }_0^1{\int }_0^{1-u}dudv=\frac{1}{2})\end{array}$$

3.3 向量F的计算

写出 $F$ 的计算式：
$$\begin{array}{rl}F& ={\int }_{\Omega }f{\phi }_i(x,y)d\Omega \\ & =f{\int }_0^1{\int }_0^{1-u}{\phi }_i(x,y)dudv\left|det(J)\right|\\ & =\frac{f\left|det(J)\right|}{6}({\int }_0^1{\int }_0^{1-u}{\phi }_i(x,y)dudv=\frac{1}{6}=四面体的体积1\times 1\times \frac{1}{2}\times \frac{1}{3})\end{array}$$

四、数值实验

本次实验选取的矩形网格大小为100x100，速度为3500m/s，矩形长宽均为400m，其他具体数据如下：

本次实验震源使用的是高斯震源，震源时间函数如下所示：

高斯震源 质量矩阵、刚度矩阵和向量F为

（九点震源）非全对角网格的刚度矩阵、质量矩阵、向量F

（九点震源）全对角网格的刚度矩阵、质量矩阵、向量F

实验结果

左三角网格

右三角网格

全对角网格

• 从图中可以看出，网格的不均匀性会导致波出现频散不一致现象，三角网格中沿着对角线频散现象较严重，在实际设计网格要避免网格出现局部区域的大部变化。
• 全对角网格模拟的网格节点数和单元数都显著增多，计算的时间明显增大。
• 上面的结果说明了，有限元能够处理二维声波的传播问题。
• 同时本文还实现了128x128的全对角三角网格的结果，其计算时间明显增大，期间计算机还出现重启的现象，说明有限元对于大型问题所需要的计算资源是极其庞大的。

其他网格结果

从上图中，可以看出由网格不均匀性导致出现反射的假象

示例代码

需要先封装好creat_grid.py。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri

def creat_mesh(x, y, nx, ny, e_type):
    '''
    x:x方向上的距离
    x:y方向上的距离
    nx:x方向上的element的数量
    ny:y方向上的element的数量
    e_type:SQ表示是矩形单元,TR表示二分三角单元,TR2另一种二分三角单元,FE表示一个矩形剖分为四个三角单元

    '''
    # 矩形单元的四角坐标
    q = np.array([[0., 0.], [x, 0.], [0, y], [x, y]])
    # node的数量
    numN = (nx+1)*(ny+1)
    # element的数量
    numE = nx*ny
    # 矩形element的角
    NofE = 4
    # 二维坐标
    D = 2
    # nodes 坐标
    NC = np.zeros([numN, D])
    # dx,dy的计算
    dx = q[1, 0]/nx
    dy = q[2, 1]/ny
    # nodes 坐标计算
    n = 0
    for i in range(1, ny+2):
        for j in range(1, nx+2):
            NC[n, 0] = q[0, 0] + (j-1)*dx
            NC[n, 1] = q[0, 1] + (i-1)*dy

            n += 1
    # element 索引，一个element由四个角的节点进行索引
    EI = np.zeros([numE, NofE])

    for i in range(1, ny+1):
        for j in range(1, nx+1):
            # 从底层开始类推
            if j == 1:
                EI[(i-1)*nx+j-1, 0] = (i-1)*(nx+1) + 1
                EI[(i-1)*nx+j-1, 1] = EI[(i-1)*nx+j-1, 0] + 1
                EI[(i-1)*nx+j-1, 3] = EI[(i-1)*nx+j-1, 0] + (nx+1)
                EI[(i-1)*nx+j-1, 2] = EI[(i-1)*nx+j-1, 3] + 1
            else:
                EI[(i-1)*nx+j-1, 0] = EI[(i-1)*nx+j-2, 1]
                EI[(i-1)*nx+j-1, 3] = EI[(i-1)*nx+j-2, 2]
                EI[(i-1)*nx+j-1, 1] = EI[(i-1)*nx+j-1, 0] + 1
                EI[(i-1)*nx+j-1, 2] = EI[(i-1)*nx+j-1, 3] + 1
    # 至此完成了矩形单元的划分工作
    # 三角形单元需要将每一个矩形单元进行拆分，即一分二成两个三角形
    if e_type == 'TR':
        # 三角形的三个角
        NofE_new = 3
        # 单元数量
        numE_new = numE * 2
        # 新的三角单元索引
        EI_new = np.zeros([numE_new, NofE_new])

        # 对矩形单元进行逐个剖分
        for i in range(1, numE+1):
            EI_new[2*(i-1), 0] = EI[i-1, 0]
            EI_new[2*(i-1), 1] = EI[i-1, 1]
            EI_new[2*(i-1), 2] = EI[i-1, 3]

            EI_new[2*(i-1)+1, 0] = EI[i-1, 1]
            EI_new[2*(i-1)+1, 1] = EI[i-1, 2]
            EI_new[2*(i-1)+1, 2] = EI[i-1, 3]
    # 至此完成了矩形单元的划分工作
        EI = EI_new
    if e_type == 'TR2':
        # 三角形的三个角
        NofE_new = 3
        # 单元数量
        numE_new = numE * 2
        # 新的三角单元索引
        EI_new = np.zeros([numE_new, NofE_new])

        # 对矩形单元进行逐个剖分
        for i in range(1, numE+1):
            EI_new[2*(i-1), 0] = EI[i-1, 0]
            EI_new[2*(i-1), 1] = EI[i-1, 1]
            EI_new[2*(i-1), 2] = EI[i-1, 2]

            EI_new[2*(i-1)+1, 0] = EI[i-1, 0]
            EI_new[2*(i-1)+1, 1] = EI[i-1, 2]
            EI_new[2*(i-1)+1, 2] = EI[i-1, 3]
    # 至此完成了矩形单元的划分工作
        EI = EI_new
    # 现在要做的是将一个矩形单元划分成四个小的三角单元
    # 这种情况下，单元总数将会是原来的四倍，同时网格节点的数量也会增加
    # 对于 (m+1) x (n+1) 大小的网格，网格节点会增加(m) x (n)
    if e_type == 'FE':
        
        # 三角形的三个角
        NofE_new = 3
        # 单元数量
        numE_new = numE*4
        # 节点数量
        numN_new = (nx)*(ny)+(nx+1)*(ny+1)
        # 新的三角单元索引
        EI_new = np.zeros([numE_new, NofE_new])
        # 新的网格节点索引
        NC_new = np.zeros([numN_new,2])
        # 新网格节点由原矩形网格节点和新生成的矩形内部网格节点组成
        NC_new[0:numN,:]=NC
        count=0
        for j in range(ny):
            for i in range(nx):
                NC_new[numN+count,0]=NC[i+(nx+1)*j,0]+dx/2
                NC_new[numN+count,1]=NC[i+(nx+1)*j,1]+dy/2
                count += 1
        for i in range(numE):
            EI_new[4*i,0]=EI[i,0]
            EI_new[4*i,1]=EI[i,1]
            EI_new[4*i,2]=numN+i+1
            
            EI_new[4*i+1,0]=EI[i,1]
            EI_new[4*i+1,1]=EI[i,2]
            EI_new[4*i+1,2]=numN+i+1
            
            EI_new[4*i+2,0]=EI[i,2]
            EI_new[4*i+2,1]=EI[i,3]
            EI_new[4*i+2,2]=numN+i+1
            
            EI_new[4*i+3,0]=EI[i,3]
            EI_new[4*i+3,1]=EI[i,0]
            EI_new[4*i+3,2]=numN+i+1
            
                    
        NC = NC_new
        EI = EI_new
    EI = EI.astype(int)
    
    return NC, EI

# -*- coding: utf-8 -*-
'''
该程序属于本人计算地球物理作业的一部分，用于计算二维有限元声波方程的正演计算
@作者：fm
'''

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import matplotlib.tri as mtri
from creat_grid import creat_mesh
import time
import os
import shutil
T1 = time.time()
# ----------------------------------------------------------
# 参数设置

eps = 0.2  # Stability limit
x = y = 800
nx = ny = 64
element_type = 'FE'
nt = 300  # Number of time steps
v = 3500
dx = x/nx
dt = eps*dx/v
# 生成网格
NC, EI = creat_mesh(x, y, nx, ny, element_type)
numN = np.size(NC, 0)
u = np.zeros([numN, 1])
uold = np.zeros([numN, 1])
unew = np.zeros([numN, 1])
# initialize time axis
t = np.arange(1, nt+1)*dt

# Initialization of the source time function
# ---------------------------------------------------------------
pt = 20*dt     # Gaussian width
t0 = 3*pt      # Time shift
src = -2/pt**2 * (t-t0) * np.exp(-1/pt**2 * (t-t0)**2)
plt.figure()
plt.plot(t, src, color='b', lw=2, label='Source time function')
plt.ylabel('Amplitude', size=16)
plt.xlabel('time', size=16)
plt.legend()
plt.grid(True)
# plt.show()
# 高斯滤波


def gauss(x, y, sigma=1.5):
    G = 1./(2.*np.pi*sigma**2)*np.exp(-(x*x+y*y)/(2*sigma**2))
    return G


def guass_matrix():
    '''
    滤波矩阵系数 3x3
    '''
    c = np.zeros([3, 3])
    for j, jj in zip(range(3), range(-1, 2, 1)):
        for i, ii in zip(range(3), range(1, -2, -1)):
            # print(j,i)
            c[j, i] = gauss(ii, jj)
    c = c/np.sum(c)
    return c


def creat_stiffness_matrix(NC, EI):
    '''
    使用参考单元法，建立刚度矩阵
    NC:节点位置列表
    EI:单元索引列表
    '''
    EI = EI-1
    # 计算节点数量
    numN = np.size(NC, 0)

    # 刚度矩阵，大小为 numN x numN
    S = np.zeros([numN, numN])
    # 循环遍历所有单元
    for ele in EI:
        # 获取三角形单元三个顶点的坐标
        x0, y0 = NC[ele[0], 0], NC[ele[0], 1]
        x1, y1 = NC[ele[1], 0], NC[ele[1], 1]
        x2, y2 = NC[ele[2], 0], NC[ele[2], 1]
        # 计算雅克比矩阵
        J = np.array([[x0-x2, x1-x2],
                      [y0-y2, y1-y2]])
        # 计算雅克比矩阵的逆
        inv_J = np.linalg.inv(J)
        # 基函数的偏导数矩阵
        partial_phi = np.array([[1, 0],
                                [0, 1],
                                [-1, -1]])
        # 初始化局部三角单元的刚度矩阵（3x3）
        S_local = np.zeros([3, 3])
        # 计算局部刚度矩阵
        # \partial \phi x inv_J
        for i in range(np.size(S_local, 1)):
            for j in range(np.size(S_local, 0)):
                S_local[j, i] = np.dot(np.dot(partial_phi[j, :].reshape(1, 2), inv_J), np.dot(
                    partial_phi[i, :].reshape(1, 2), inv_J).T)*np.linalg.det(J)/2.
                # print(inv_J,partial_phi[i, :],inv_J.shape,partial_phi[i, :].shape)
                # 组装总的刚度矩阵
        S[np.ix_(ele, ele)] += S_local
    print('刚度矩阵组装完成\nsize={}x{}'.format(numN, numN))
    return S


def creat_mass_matrix(NC, EI):
    '''
    使用参考单元法，建立质量矩阵
    NC:节点位置列表
    EI:单元索引列表
    '''
    EI = EI-1
    # 计算节点数量
    numN = np.size(NC, 0)
    numE = np.size(EI, 0)
    # 初始化总质量矩阵
    M = np.zeros([numN, numN])
    # 初始化局部质量矩阵
    M_local = np.zeros([3, 3])
    # 循环遍历所有单元
    for ele in EI:
        # 获取三角形单元三个顶点的坐标
        x0, y0 = NC[ele[0], 0], NC[ele[0], 1]
        x1, y1 = NC[ele[1], 0], NC[ele[1], 1]
        x2, y2 = NC[ele[2], 0], NC[ele[2], 1]

        # 计算雅克比矩阵
        J = np.array([[x0-x2, x1-x2],
                      [y0-y2, y1-y2]])
        # 计算局部质量矩阵 3x3
        M_local = 1/24.*np.linalg.det(J)*np.ones_like(M_local)
        M_local[0, 0] = M_local[1, 1] = M_local[2, 2] = 1/12.*np.linalg.det(J)
        # 组装总质量矩阵
        M[np.ix_(ele, ele)] += M_local
    print('组装质量矩阵完成\nsize={}x{}'.format(numN, numN))
    return M


def creat_right_vector(NC, EI, filter):
    '''
    建立右端向量F函数
    F是为波场添加一个扰动的项
    建立F时假设全域均有相同的单位扰动，
    后续使用时可以利用适当的滤波器或delta函数来施加有限区域的扰动
    filter: 表示源加载的方案
        one_point: 表示源加载在一个点上
        9 points 表示源加载在高斯滤波后的九点上
    '''
    EI = EI-1
    # 计算节点数量
    numN = np.size(NC, 0)
    # 初始化总的F
    F = np.zeros([numN, 1])
    # 初始化局部F
    F_local = np.zeros([3, 1])
    # 遍历filter 找到非零的扰动项
    for ele in EI:
        # 获取三角形单元三个顶点的坐标
        x0, y0 = NC[ele[0], 0], NC[ele[0], 1]
        x1, y1 = NC[ele[1], 0], NC[ele[1], 1]
        x2, y2 = NC[ele[2], 0], NC[ele[2], 1]
        # 计算雅克比矩阵
        J = np.array([[x0-x2, x1-x2],
                      [y0-y2, y1-y2]])
        # 计算局部F
        F_local = filter[ele]/6*np.linalg.det(J)
        # 组装总的F
        F[ele] += F_local

    print('F向量组装完成\nF.size={}x{}'.format(numN, 1))
    return F


def give_a_source(NC, EI, points=1, nx=nx, ny=ny):
    '''
    在区域中心施加一个震源，有九点震源可供选择
    '''
    numN = np.size(NC, 0)
    filter = np.zeros([numN, 1])
    mid = int(numN/2)
    filter[mid, 0] = 1
    if points == 9:
        c = guass_matrix()
        filter[mid, 0] = c[1, 1]
        filter[mid+1, 0] = filter[mid-1, 0] = c[0, 1]
        filter[mid+nx+1, 0] = filter[mid-nx-1, 0] = c[1, 0]
        filter[mid+nx+1+1, 0] = filter[mid+nx+1-1, 0]\
            = filter[mid-nx-1+1, 0] = filter[mid-nx-1-1, 0] = c[0, 0]

    if element_type == 'FE':
        filter_new = np.zeros([numN, 1])
        mid = int((nx+1)*(ny+1)/2)
        filter_new[mid, 0] = 1
        filter = filter_new
    return filter
# 计算边界节点索引


def caculate_BC_index(nx=nx, ny=ny):
    '''
    计算矩形区域周围的节点索引
    '''
    # 大矩形的四条边
    l1 = range(0, nx+1)
    l2 = range(0, (nx+1)*ny, nx+1)
    l3 = range((nx+1)*(ny+1)-1, (nx+1)*(ny)-1, -1)
    l4 = range(nx, (nx+1)*(ny+1)-1, nx+1)
    # 返回四条边的索引的唯一值
    return np.unique(np.concatenate([l1, l2, l3, l4], axis=0))
def RemoveDir(filepath):
    '''
    如果文件夹不存在就创建，如果文件存在就清空！
    
    '''
    if not os.path.exists(filepath):
        os.mkdir(filepath)
    else:
        shutil.rmtree(filepath)
        os.mkdir(filepath)

S = creat_stiffness_matrix(NC=NC, EI=EI)
M = creat_mass_matrix(NC=NC, EI=EI)
F = creat_right_vector(
    NC=NC, EI=EI, filter=give_a_source(NC=NC, EI=EI, points=1))
plt.figure()
plt.subplot(221)
plt.imshow(S[0:10, 0:10])
plt.title('Finite-Element Stiffness Matrix $\mathbf{S}$')
plt.axis("off")
plt.subplot(222)
plt.imshow(M[0:10, 0:10])
plt.title('Finite-Element Mass Matrix $\mathbf{M}$')
plt.axis("off")
plt.subplot(212)
plt.title('Finite-Element F vector $\mathbf{F}$')
plt.plot(F)

Minv = np.linalg.inv(M)

bci = caculate_BC_index()
plt.figure()
x = np.squeeze(NC[:, 0])
y = np.squeeze(NC[:, 1])
tri = EI-1
triang = mtri.Triangulation(x, y, tri)
unewp = plt.tripcolor(triang, np.squeeze(unew))
filepath='images'
RemoveDir(filepath)
# 利用欧拉时间步进计算波场值
for it in range(nt):
    unew = (dt**2)*Minv @ (F*src[it]-v**2*S @ u) + 2*u - uold

    # boundary condition
    unew[bci] = 0
    uold, u = u, unew
    if not it % 20:
        unewp = plt.tripcolor(triang, np.array(np.squeeze(unew)))
        plt.axis('equal')
        plt.title('t={:.2f}s'.format(it*dt))
        plt.xlabel('distance/m')
        plt.ylabel('distance/m')
        plt.savefig(filepath+'/{}.png'.format(it))
        print('{}%{}'.format(it, nt))

np.save('unew.npy', unew)

T2 = time.time()
print('程序运行时间:%s毫秒' % ((T2 - T1)*1000))
plt.show()

感谢

感谢互联网其他作者给予的灵感和帮助！