二阶偏微分方程组龙格库塔法_有限单元法(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中，前三行包含三角形三个角点的索引，按逆时针顺序给出，第四行包含子域编号。

上述代码生成的网格为：

非结构化网格示意图

这样一个网格包含了4465个节点和8728个单元，已经很密集了，再大一些计算起来就会很吃力，如果本身计算机内存不大，算力不够的话，可以减少单元数。

使用下面的代码可以查看网格及节点、单元编号，这里先把网格调稀疏点：

expand = 1000;
p = expand*p;
%'NodeLabels','on'——可以显示节点编号，'ElementLabels','on'——显示单元编号
pdemesh(p, e, t,'NodeLabels','on','ElementLabels','on');

这里对p乘上expand是因为MATLAB生成的网格范围默认是

的，我们需要乘上一个系数把网格范围调大。最终得到的图：

网格节点及单元编号分布情况

到这里大家应该就明白了，为什么向量

会出现跳跃了，我们是

按照节点编号顺序求解的，得到的值自然会在计算域内出现跳跃。但是我们会找到一个正确的对应关系，来保证我们计算的每一个值是和节点一一对应的。

2.2.3 求解质量矩阵和刚度矩阵

我们已经得到了这样一个非结构化网络，并且每一个节点和单元的编号我们也有，接下来就是求质量矩阵和刚度矩阵了。

首先我们求解质量矩阵和刚度矩阵并不是一个节点一个节点的求，而是一个单元一个单元的求，在一个单元里面包含三个节点，三个节点两两组合后构成一个

的矩阵，然后在根据这三个节点对应的编号组装到整体的

的矩阵中。这是什么意思呢？给大家一个直观的例子，例如我求一个由编号

这三个节点构成的三角形单元

的质量矩阵（参照上图：

网格节点及单元编号分布情况左下角单元），这个矩阵长这样：

单元[7]的质量矩阵

可以看到这个矩阵其实是对称的。之后将这个小矩阵组装进大矩阵，即总质量矩阵，得到大矩阵长这个样子的：

单元[7]组装进总质量矩阵

不知道大家有没有注意到，我们之前提到过单元与单元之间存在节点的重叠，这种重叠就体现在组装矩阵的时候，不同单元节点值的叠加，例如上面的例子中，节点

同时被单元

、

六个单元所共享，最终得到的矩阵

这个元素是这六个单元的贡献值之和。而这种

叠加现象只出现在对角元素上。

那么这个小质量矩阵或刚度矩阵怎么求呢？这里我推荐使用参考单元的方法来求解，因为好理解而且也很好计算。

我们发现，对于上面的非结构化网络，三角单元的节点坐标值不具有规律性，我们如果能用一种规则的、好计算的三角形来求解，例如等腰直角三角形，那我们就好求出每一个小矩阵。将

组成的单元映射到

参考单元：

参考单元及坐标映射

这其实就相当于一种坐标映射，将不规则的三角形映射为规则的三角形，这种映射关系为：

其中的

为坐标映射系数：

令

，则参考单元中三个节点对应的插值基函数为：

有了插值基函数表达式后还没完，我们还要进行积分算子

的坐标转换，对于一个通用的坐标转换例子：

我们从

域变换到

域，使用微分的链式法则：

则有如下关系：

式中的

代表雅可比矩阵行列式，其数值与参与变换的三角形节点坐标有关：

可以发现，通过上面的参考单元，我们就能够找到一种通用的方法来计算每个单元的质量矩阵和刚度矩阵，并且积分区域都统一变成了一样的。

（1）质量矩阵

首先，对于每个单元的质量矩阵，其中的元素计算公式变换为：

这里由于我们已经知道了参考单元每一个节点的基函数表达式，我们可以把雅克比行列式从积分中提出来，将关于基函数的积分计算出来，即：

这个公式说明什么，说明经过参考单元变换后，我们刚刚提到的质量矩阵其实计算起来只需要求

就可以了，形象一点就是：

参考单元下计算单元质量矩阵

细心的你可能会注意到，这里的编号统一都是

，和前面的单元

对应关系为：

，按照逆时针方向一一对应。

（2）刚度矩阵

刚度矩阵要稍微复杂一点点，因为这里面涉及到对基函数的偏导，将这个偏导写成向量形式：

使用偏微分的链式法则：

将公式（31）代入（30）中，则有：

其中的矩阵

代表雅克比矩阵的逆。对参考单元的基函数求导：

最终得到单元每一个元素的刚度矩阵计算公式为：

为了简化计算，我们发现积分公式里面没有任何一项和

有关，即：

这相当于告诉我们：

用这个公式就好计算多了。

2.3 算法实现

（1）首先生成网格，这里封装为函数create_grid()：

function GRID = create_grid()
% This function is used to generate  structured and unstructured mesh
GRID = struct('unstructured',@unstructured,...'structured',@structured,...'connect_mat',@connect_mat);function [p,e,t] = unstructured()[p, e, t] = initmesh('squareg', 'hmax', 1);%单元初始化[p,e,t] = refinemesh('squareg',p,e,t);%单元加密endfunction [p,e,t] = structured()%% 参数设置Lx = 1000; %定义单元右边界（左边界为0，如果不是，可以平移为0）Ly = 1000;%定义单元上边界N = 60;%分割的一个方向的单元数目numelx = N;%定义分割的x方向单元数目（按矩形计算）numely = N;%定义分割的y方向单元数目（按矩形计算）numnodx = numelx + 1; % x方向节点个数比单元个数多1numnody = numely + 1; % y方向节点个数比单元个数多1nel = 3;%每个单元的节点数目,即每个单元上有几个形函数参与作用，单元自由度coordx = linspace(0,Lx,numnodx)'; %等分节点的坐标（为了方便，我这里采用等分的方式，事实上单元长度可以不一致，非均匀网格）coordy = linspace(0,Ly,numnody)'; %等分节点的坐标（为了方便，我这里采用等分的方式，事实上单元长度可以不一致）[X, Y] = meshgrid(coordx,coordy);%张成网格，X和Y分别表示对应位置的横纵坐标X = X';Y = Y';p = [X(:) Y(:)];%把网格一行一行扯开，coord的每一行是对应节点的坐标，按顺序排p = p';connect = connect_mat(numnodx,numnody,nel);%连接矩阵，表示每个单元周围的节点编号，也就是涉及的形函数编号t = connect';e = [1:numnodx numnodx*numely+1:numnodx*numely+numnodx ...numnodx+1:numnodx:numnodx*(numely-1)+1 2*numnodx:numnodx:numely*numnodx]; % 强制性边界点的编号，本例子中是四条边，下上左右边endfunction connect_mat = connect_mat( numnodx,numnody,nel)%输入横纵坐标的节点数目，和单元自由度%输出连接矩阵，每个单元涉及的节点的编号xn = 1:(numnodx*numnody);%拉成一条编号A = reshape(xn,numnodx,numnody);%同形状编号for i = 1:(numnodx-1)*(numnody-1)xg = rem(i,numnodx-1);%xg表示单元为左边界数起第几个if xg == 0xg = numnodx-1;endyg = ceil(i/(numnodx-1));%下边界其数第几个a = A(xg:xg+1,yg:yg+1);%这个小矩阵，拉直了就是连接矩阵a_vec = a(:);connect_mat(2*i-1:2*i,1:nel) = [a_vec([1 4 3])';a_vec([4 1 2])'];endend
end

函数里面除了可以实现非结构化网格，还提供了了等腰直角三角形形式的结构化网络可选。

（2）接着编写计算质量矩阵和刚度矩阵的函数FEM_2D_func()：

function FEM_2D_function = FEM_2D_func()%设置所有需要的函数FEM_2D_function = struct('assemble_matrix_2D',@assemble_matrix_2D,...%组装刚度矩阵方法一'assemble_matrix_2D2',@assemble_matrix_2D2,...%组装刚度矩阵方法二'assemble_vector_2D',@assemble_vector_2D,...%组装右端向量'treat_boundary_condition',@treat_boundary_condition);%处理边界条件%% 实现——组装刚度矩阵方法一function A = assemble_matrix_2D(p, t)fprintf('组装刚度矩阵:n');%获取单元个数和节点个数number_of_elements = length(t);number_of_nodes = length(p);%初始化总刚度矩阵A = zeros(number_of_nodes, number_of_nodes);for n = 1: number_of_elements% 获取单元局部编号到全局编号之间的对应关系local2global = t(1: 3, n);vertices = p(:, local2global);x = vertices(1, :);y = vertices(2, :);a = 0.5*((x(2)*y(3)-x(3)*y(2))-(y(3)-y(2))*x(1)+y(1)*(x(3)-x(2)));%三角形单元面积a = abs(a);a = 1/(2*a);A_local = zeros(3,3);ph = [1,0;0,1;-1,-1];detJ = (x(3)-x(1))*(y(2)-y(3))-(y(3)-y(1))*(x(2)-x(3));J_inv = a*[y(2)-y(3) , x(3)-x(2);y(3)-y(1) , x(1)-x(3)];for i = 1: 3for j = 1: 3A_local(i, j) = -0.5*detJ.*(ph(i,:)*J_inv*(J_inv'*ph(j,:)'));endendA(local2global, local2global) = A(local2global, local2global) + A_local;if(mod(n,100)==0)fprintf('number of  elements:%d/%dn',n,number_of_elements)endendfprintf('刚度矩阵组装完成n')end
%% 组装刚度矩阵方法二——参考单元function A = assemble_matrix_2D2(p, t)fprintf('组装刚度矩阵:n');%获取单元个数和节点个数number_of_elements = length(t);number_of_nodes = length(p);%初始化总刚度矩阵A = zeros(number_of_nodes, number_of_nodes);N{1} = @(xi, eta) 1 - xi - eta; N_xi{1} = -1; N_eta{1} = -1;N{2} = @(xi, eta) xi; N_xi{2} = 1; N_eta{2} = 0;N{3} = @(xi, eta) eta; N_xi{3} = 0; N_eta{3} = 1;ymax = @(xi) 1 - xi;for n = 1: number_of_elements% 获取单元局部编号到全局编号之间的对应关系local2global = t(1: 3, n);vertices = p(:, local2global);xx = vertices(1, :); yy = vertices(2, :);J = [xx(1) * N_xi{1} + xx(2) * N_xi{2} + xx(3) * N_xi{3}, xx(1) * N_eta{1} + xx(2) * N_eta{2} + xx(3) * N_eta{3};yy(1) * N_xi{1} + yy(2) * N_xi{2} + yy(3) * N_xi{3}, yy(1) * N_eta{1} + yy(2) * N_eta{2} + yy(3) * N_eta{3}];detJ = abs(det(J));J_inv_T = inv(J);A_local = zeros(3, 3);for i = 1: 3for j = 1: 3fun = (N_xi{i} * J_inv_T(1, 1) + N_eta{i} * J_inv_T(2, 1))...* (N_xi{j} * J_inv_T(1, 1) + N_eta{j} * J_inv_T(2, 1))...+ (N_xi{i} * J_inv_T(1, 2) + N_eta{i} * J_inv_T(2, 2))...* (N_xi{j} * J_inv_T(1, 2) + N_eta{j} * J_inv_T(2, 2));A_local(i, j) = integral2(@(xi, eta) fun .* eta ./ eta, 0, 1, 0, ymax) * detJ;endendA(local2global, local2global) = A(local2global, local2global) + A_local;if(mod(n,100)==0)fprintf('number of  elements:%d/%dn',n,number_of_elements)endendfprintf('刚度矩阵组装完成n')end
%% 组装质量矩阵M和右端向量Ffunction [M,F] = assemble_vector_2D(p, t)fprintf('组装质量矩阵:n');number_of_elements = length(t);number_of_nodes = length(p);M = zeros(number_of_nodes, number_of_nodes); % 总质量矩阵Me = zeros(3, 3);F = zeros(number_of_nodes, 1);            % 初始化右端向量Fe = zeros(3, 1);                         % 初始单元右端向量% 单元质量矩阵N{1} = @(xi, eta) 1 - xi - eta; N_xi{1} = -1; N_eta{1} = -1;N{2} = @(xi, eta) xi; N_xi{2} = 1; N_eta{2} = 0;N{3} = @(xi, eta) eta; N_xi{3} = 0; N_eta{3} = 1;ymax = @(xi) 1 - xi;for i = 1: number_of_elementslocal2global = t(1: 3, i);%获取当前单元包含的节点编号vertices = p(:, local2global);%获取当前单元所有节点的x,y坐标% 从参考单元到物理单元的映射%     x = @(xi, eta) xx(1) * N{1}(xi, eta) + xx(2) * N{2}(xi, eta) + xx(3) * N{3}(xi, eta);%     y = @(xi, eta) yy(1) * N{1}(xi, eta) + yy(2) * N{2}(xi, eta) + yy(3) * N{3}(xi, eta);xx = vertices(1, :); yy = vertices(2, :);J = [xx(1) * N_xi{1} + xx(2) * N_xi{2} + xx(3) * N_xi{3}, xx(1) * N_eta{1} + xx(2) * N_eta{2} + xx(3) * N_eta{3};yy(1) * N_xi{1} + yy(2) * N_xi{2} + yy(3) * N_xi{3}, yy(1) * N_eta{1} + yy(2) * N_eta{2} + yy(3) * N_eta{3}];detJ = abs(det(J));%积分太慢了
%             for m=1:3
%                 for n=1:3
%                     Me(m,n) = integral2(@(xi,eta)N{m}(xi,eta).*N{n}(xi,eta),0,1,0,ymax)*detJ;
%                 end
%                 Fe(m,1) = integral2(@(xi,eta) N{m}(xi,eta), 0, 1,0,ymax) * detJ;
%             endMe(1,1) = 1/12*detJ;Me(1,2) = 1/24*detJ;Me(1,3) = 1/24*detJ;Me(2,1) = 1/24*detJ;Me(2,2) = 1/12*detJ;Me(2,3) = 1/24*detJ;Me(3,1) = 1/24*detJ;Me(3,2) = 1/24*detJ;Me(3,3) = 1/12*detJ;Fe(:,1) = 1/6*detJ;M(local2global, local2global) = M(local2global, local2global) + Me;F(local2global, 1) = F(local2global, 1) + Fe;if(mod(i,100)==0)fprintf('number of  elements:%d/%dn',i,number_of_elements)endendfprintf('质量矩阵组装完成n')end
%% 边界条件的处理function [M,S,F] = treat_boundary_condition(M,S,F, p, boundarynodes, e)p = p';number_of_boundarynodes = length(boundarynodes);for i = 1: length(e)if p(e(1, i), 1) == -1 && p(e(2, i), 1) == -1f1 = @(x) (p(e(2, i), 2) - x) ./ (p(e(2, i), 2) - p(e(1, i), 2));f2 = @(x) (x - p(e(1, i), 2)) ./ (p(e(2, i), 2) - p(e(1, i), 2));F(e(1, i)) = F(e(1, i)) + integral(f1, p(e(1, i), 2), p(e(2, i), 2));F(e(2, i)) = F(e(2, i)) + integral(f2, p(e(1, i), 2), p(e(2, i), 2));endendfor i = 1: number_of_boundarynodesif p(boundarynodes(i), 2) == -1 || p(boundarynodes(i), 2) == 1M(boundarynodes(i), :) = 0;M(boundarynodes(i), boundarynodes(i)) = 1;F(boundarynodes(i)) = 0;S(boundarynodes(i), :) = 0;S(boundarynodes(i), boundarynodes(i)) = 1;endendendend

这个函数里面还提供了另外一种求解刚度矩阵的方法可选（没有使用到参考单元），并求解右端向量F、设置边界条件(可选)。

（3）接着编写主函数FEM_2D_main():

clear;clc
t_pre = clock;
t_start = clock;
%% 生成单元
GRID = create_grid();
[p,e,t] = GRID.unstructured();
expand = 1000;
p = expand*p;
%'NodeLabels','on'——可以显示节点编号，'ElementLabels','on'——显示单元编号
pdemesh(p, e, t,'NodeLabels','on','ElementLabels','on');%生成结构化三角网格
% GRID = create_grid();
% [p,e,t] = GRID.structured();num_nodes = length(p);
num_elements = length(t);
t_s = sprintf('nodes=%d  elements=%d',num_nodes,num_elements);
title(t_s)
% 在点矩阵p中，第一行和第二行包含网格中点的x和y坐标。
% 在边矩阵e中，第一行和第二行包含起始点和结束点的索引，第三和第四行包含起始和结束参数值，
%第五行包含边缘段编号，第六和第七行包含左侧和右侧子域编号。
% 在三角形矩阵t中，前三行包含角点的索引，按逆时针顺序给出，第四行包含子域编号。%% 参数设置
nt = 150; dt = 0.003;
v = 1500;
T = (1:nt)*dt;
number_of_notes = length(p);
fmain = 40;%% 设置震源
s_t = (1-2*pi^2*fmain*(T-0.2).^2).*exp(-fmain*pi^2*(T-0.2).^2);%雷克子波%% 调用函数计算矩阵
FEM = FEM_2D_func();
S = FEM.assemble_matrix_2D(p, t);%刚度矩阵
[M,F]= FEM.assemble_vector_2D(p, t);%质量矩阵和右端向量
boundarynodes = unique([e(1, :) e(2, :)]);
U = zeros(number_of_notes,nt);
[m,source_x] = find(abs(p(1,:))>=0&abs(p(1,:))<=50&abs(p(2,:))>=0&abs(p(2,:))<=50);%% 利用递推关系求波场值
for i = 2:nt-1U(source_x(1),i) = s_t(i);%隐式
%     U(:,i+1) = (M+v^2*dt^2*S)(M*(2*U(:,i)-U(:,i-1)));
%     %显式right_vector = v^2*dt^2*S*U(:,i);U(:,i+1) = 2*U(:,i)-U(:,i-1)-Mright_vector;U(boundarynodes,i+1) = 0;if(mod(i,10)==0)fprintf('time step=%d total=%dsn',i,nt);end
end
fprintf('time step=%d total=%dsn',nt,nt);
total_time = etime(clock,t_start)/60;
fprintf('total runtime %.2fminutes',total_time)%% 绘图
filename = sprintf('FEM-2D dt=%.3f fm=%.1f v=%d explicit.gif',dt,fmain,v);
pic_num = 1;
for i=5:5:nttrisurf(t(1: 3, :)', p(1, :)', p(2, :)', U(:,i))str = sprintf('time step=%dndt=%.3f fm=%.1f v=%d',i,dt,fmain,v);title(str)colorbarshading interpview([90,90])F = getframe(gcf);I = frame2im(F,256);[I,map]=rgb2ind(I,256);if pic_num == 1imwrite(I,map,filename,'gif','Loopcount',inf,'DelayTime',0.2);elseimwrite(I,map,filename,'gif','WriteMode','append','DelayTime',0.2);endpic_num = pic_num + 1;
end
%% 获取波场快照
trisurf(t(1: 3, :)', p(1, :)', p(2, :)', U(:,150))
str = sprintf('FEM-2D-hemogenous t=%d ms',150*3);
title(str)
shading interp
view([90,90])
colormap('gray')
saveas(gcf, 'FEM-2D-unstructued.png', 'png')

主函数里面可选这使用非结构化和结构化网格，并有隐式差分格式和显式差分格式可选。

2.4 模拟结果

隐式差分格式模拟结果

显式差分格式模拟结果

两种都存在一定的数值耗散和色散现象。

3 声波方程有限单元法数值模拟总结

通过这两篇文章，相信大家对于有限单元法都有了一定的了解，我本来还想向大家介绍一下有关有限元法稳定性分析的，但是本来篇幅就太长了，以后如果感兴趣的人比较多，我再更新这一块吧！

另外，我之前都没有在文章中求赞求关注啥的，不过这次我还是厚脸皮一下，如果你喜欢这篇文章，觉得对你有用的话，请点击一下喜欢，欢迎你关注我，我会在后续继续分享有限体积法以及我所做过的一些比较有意思的数值模拟。

ps:文章中如有纰漏，欢迎在评论区指出。