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:文章中如有紕漏，歡迎在評論區指出。

總結

以上是生活随笔為你收集整理的二阶偏微分方程组 龙格库塔法_有限单元法(Finite Element Method)实现声波方程模拟（Part 2）...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。