二阶偏微分方程组 龙格库塔法_有限单元法(Finite Element Method)实现声波方程模拟（Part 2）...

2.1 前言

承接上一篇文章，前面我们已经介绍了一维声波方程有限元求解：

蓝不是蓝：有限单元法(Finite Element Method)实现声波方程模拟（Part 1）​zhuanlan.zhihu.com

这一部分将一维问题提升到二维问题。不知道大家有没有发现，在一维问题中，我们是通过矩阵相乘的方式来求解，得到的解

是一个关于节点编号的向量，即对于一个研究区域，将该研究区域划分为有限个单元，两个单元间有共用的节点，假设一共有

个节点并编号为

，那么我们实际上求解得到的是向量

。在二维问题中，也是同样的求解方法，只不过这里的节点编号不再是按照单方向编号（如一维中沿着

方向），而是按照一定规律编号（如按列或按行），编号规则可人为设定。那么这样得到的一个向量

就会在空间上出现跳跃，但是这并不会影响我们的求解。如果这里没有理解，先继续往后面看。

2.2 二维声波方程有限元求解

2.2.1 算法推导

对于二维声波方程，我们先写成一个简写形式：

我们同样使用伽辽金（Galerkin ）方法将该微分方程变成弱形式。并引入一个二维基函数

，来实现对解的插值近似。这里我们假设

。

对上式关于空间项的积分使用分部积分法则：

将上式代入公式（14）后得到：

这里我们仍然使用一维情况下类似的边界条件：

消掉边界条件这一项后，再利用基函数插值近似我们的解（

代表节点数

）：

最终得到的公式为：

时间项我们采用二阶差分格式，则上式的隐式差分格式为：

简写成矩阵形式：

同样的，简写成矩阵形式的显式差分格式为：

其中的质量矩阵

和刚度矩阵

为：

那么，我们先利用公式（21）和（22）求出

和

之后，在利用公式（20.1）或（20.2）就可以得到计算域中每一个节点的值。

2.2.2 网格剖面

在这里我们使用三角形单元构建一个非结构化网格，事实上如果让我们自己来构建一个很复杂的非结构化网格是比较费力的，我们可以使用MATLAB自带的网格生成函数来帮助我们完成网格生成。

这部分代码为：

function [p,e,t] = unstructured()
        [p, e, t] = initmesh('squareg', 'hmax', 0.08);%单元初始化
        [p,e,t] = refinemesh('squareg',p,e,t);%单元加密
end

其中

完成网格初始化，0.08衡量网格疏密程度，越小网格越密。

完成网格的加密。函数返回值

的含义见下表：

返回值	含义
p	在点矩阵p中，第一行和第二行包含网格中节点编号对应的x和y坐标。
e	在边界矩阵e中，第一行和第二行包含每一段边界起始点和结束点的索引，第三和第四行包含起始和结束参数值，第五行包含边界段编号，第六和第七行包含左侧和右侧子域编号。
t	在三角形单元矩阵t中，前三行包含三角形三个角点的索引，按逆时针顺序给出，第四行包含子域编号。

上述代码生成的网格为：