2.1 前言

承接上一篇文章，前面我們已經(jīng)介紹了一維聲波方程有限元求解：

藍(lán)不是藍(lán)：有限單元法(Finite Element Method)實(shí)現(xiàn)聲波方程模擬（Part 1）?zhuanlan.zhihu.com

這一部分將一維問題提升到二維問題。不知道大家有沒有發(fā)現(xiàn)，在一維問題中，我們是通過矩陣相乘的方式來求解，得到的解

是一個(gè)關(guān)于節(jié)點(diǎn)編號(hào)的向量，即對(duì)于一個(gè)研究區(qū)域，將該研究區(qū)域劃分為有限個(gè)單元，兩個(gè)單元間有共用的節(jié)點(diǎn)，假設(shè)一共有 個(gè)節(jié)點(diǎn)并編號(hào)為 ，那么我們實(shí)際上求解得到的是向量 。在二維問題中，也是同樣的求解方法，只不過這里的節(jié)點(diǎn)編號(hào)不再是按照單方向編號(hào)（如一維中沿著 方向），而是按照一定規(guī)律編號(hào)（如按列或按行），編號(hào)規(guī)則可人為設(shè)定。那么這樣得到的一個(gè)向量 就會(huì)在空間上出現(xiàn)跳躍，但是這并不會(huì)影響我們的求解。如果這里沒有理解，先繼續(xù)往后面看。

2.2 二維聲波方程有限元求解

2.2.1 算法推導(dǎo)

對(duì)于二維聲波方程，我們先寫成一個(gè)簡寫形式：

我們同樣使用伽遼金（Galerkin ）方法將該微分方程變成弱形式。并引入一個(gè)二維基函數(shù)

，來實(shí)現(xiàn)對(duì)解的插值近似。這里我們假設(shè) 。

對(duì)上式關(guān)于空間項(xiàng)的積分使用分部積分法則：

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

這里我們?nèi)匀皇褂靡痪S情況下類似的邊界條件：

消掉邊界條件這一項(xiàng)后，再利用基函數(shù)插值近似我們的解（

代表節(jié)點(diǎn)數(shù)）：

最終得到的公式為：

時(shí)間項(xiàng)我們采用二階差分格式，則上式的隱式差分格式為：

簡寫成矩陣形式：

同樣的，簡寫成矩陣形式的顯式差分格式為：

其中的質(zhì)量矩陣

和剛度矩陣 為：

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

和 之后，在利用公式（20.1）或（20.2）就可以得到計(jì)算域中每一個(gè)節(jié)點(diǎn)的值。

2.2.2 網(wǎng)格剖面

在這里我們使用三角形單元構(gòu)建一個(gè)非結(jié)構(gòu)化網(wǎng)格，事實(shí)上如果讓我們自己來構(gòu)建一個(gè)很復(fù)雜的非結(jié)構(gòu)化網(wǎng)格是比較費(fèi)力的，我們可以使用MATLAB自帶的網(wǎng)格生成函數(shù)來幫助我們完成網(wǎng)格生成。

這部分代碼為：

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

其中

完成網(wǎng)格初始化，0.08衡量網(wǎng)格疏密程度，越小網(wǎng)格越密。 完成網(wǎng)格的加密。函數(shù)返回值 的含義見下表：返回值含義

p	在點(diǎn)矩陣p中，第一行和第二行包含網(wǎng)格中節(jié)點(diǎn)編號(hào)對(duì)應(yīng)的x和y坐標(biāo)。
e	在邊界矩陣e中，第一行和第二行包含每一段邊界起始點(diǎn)和結(jié)束點(diǎn)的索引，第三和第四行包含起始和結(jié)束參數(shù)值，第五行包含邊界段編號(hào)，第六和第七行包含左側(cè)和右側(cè)子域編號(hào)。
t	在三角形單元矩陣t中，前三行包含三角形三個(gè)角點(diǎn)的索引，按逆時(shí)針順序給出，第四行包含子域編號(hào)。

上述代碼生成的網(wǎng)格為：

非結(jié)構(gòu)化網(wǎng)格示意圖

這樣一個(gè)網(wǎng)格包含了4465個(gè)節(jié)點(diǎn)和8728個(gè)單元，已經(jīng)很密集了，再大一些計(jì)算起來就會(huì)很吃力，如果本身計(jì)算機(jī)內(nèi)存不大，算力不夠的話，可以減少單元數(shù)。

使用下面的代碼可以查看網(wǎng)格及節(jié)點(diǎn)、單元編號(hào)，這里先把網(wǎng)格調(diào)稀疏點(diǎn)：

expand = 1000; p = expand*p; %'NodeLabels','on'——可以顯示節(jié)點(diǎn)編號(hào)，'ElementLabels','on'——顯示單元編號(hào) pdemesh(p, e, t,'NodeLabels','on','ElementLabels','on');

這里對(duì)p乘上expand是因?yàn)镸ATLAB生成的網(wǎng)格范圍默認(rèn)是

的，我們需要乘上一個(gè)系數(shù)把網(wǎng)格范圍調(diào)大。最終得到的圖：

網(wǎng)格節(jié)點(diǎn)及單元編號(hào)分布情況

到這里大家應(yīng)該就明白了，為什么向量

會(huì)出現(xiàn)跳躍了，我們是按照節(jié)點(diǎn)編號(hào)順序求解的，得到的值自然會(huì)在計(jì)算域內(nèi)出現(xiàn)跳躍。但是我們會(huì)找到一個(gè)正確的對(duì)應(yīng)關(guān)系，來保證我們計(jì)算的每一個(gè)值是和節(jié)點(diǎn)一一對(duì)應(yīng)的。

2.2.3 求解質(zhì)量矩陣和剛度矩陣

我們已經(jīng)得到了這樣一個(gè)非結(jié)構(gòu)化網(wǎng)絡(luò)，并且每一個(gè)節(jié)點(diǎn)和單元的編號(hào)我們也有，接下來就是求質(zhì)量矩陣和剛度矩陣了。

首先我們求解質(zhì)量矩陣和剛度矩陣并不是一個(gè)節(jié)點(diǎn)一個(gè)節(jié)點(diǎn)的求，而是一個(gè)單元一個(gè)單元的求，在一個(gè)單元里面包含三個(gè)節(jié)點(diǎn)，三個(gè)節(jié)點(diǎn)兩兩組合后構(gòu)成一個(gè)

的矩陣，然后在根據(jù)這三個(gè)節(jié)點(diǎn)對(duì)應(yīng)的編號(hào)組裝到整體的 的矩陣中。這是什么意思呢？給大家一個(gè)直觀的例子，例如我求一個(gè)由編號(hào) 這三個(gè)節(jié)點(diǎn)構(gòu)成的三角形單元 的質(zhì)量矩陣（參照上圖：網(wǎng)格節(jié)點(diǎn)及單元編號(hào)分布情況左下角單元），這個(gè)矩陣長這樣：

單元[7]的質(zhì)量矩陣

可以看到這個(gè)矩陣其實(shí)是對(duì)稱的。之后將這個(gè)小矩陣組裝進(jìn)大矩陣，即總質(zhì)量矩陣，得到大矩陣長這個(gè)樣子的：

單元[7]組裝進(jìn)總質(zhì)量矩陣

不知道大家有沒有注意到，我們之前提到過單元與單元之間存在節(jié)點(diǎn)的重疊，這種重疊就體現(xiàn)在組裝矩陣的時(shí)候，不同單元節(jié)點(diǎn)值的疊加，例如上面的例子中，節(jié)點(diǎn)

同時(shí)被單元 、 、 、 、 、 六個(gè)單元所共享，最終得到的矩陣 這個(gè)元素是這六個(gè)單元的貢獻(xiàn)值之和。而這種疊加現(xiàn)象只出現(xiàn)在對(duì)角元素上。

那么這個(gè)小質(zhì)量矩陣或剛度矩陣怎么求呢？這里我推薦使用參考單元的方法來求解，因?yàn)楹美斫舛乙埠芎糜?jì)算。

我們發(fā)現(xiàn)，對(duì)于上面的非結(jié)構(gòu)化網(wǎng)絡(luò)，三角單元的節(jié)點(diǎn)坐標(biāo)值不具有規(guī)律性，我們?nèi)绻苡靡环N規(guī)則的、好計(jì)算的三角形來求解，例如等腰直角三角形，那我們就好求出每一個(gè)小矩陣。將

組成的單元映射到 參考單元：

參考單元及坐標(biāo)映射

這其實(shí)就相當(dāng)于一種坐標(biāo)映射，將不規(guī)則的三角形映射為規(guī)則的三角形，這種映射關(guān)系為：

其中的

為坐標(biāo)映射系數(shù)：

令

，則參考單元中三個(gè)節(jié)點(diǎn)對(duì)應(yīng)的插值基函數(shù)為：

有了插值基函數(shù)表達(dá)式后還沒完，我們還要進(jìn)行積分算子

的坐標(biāo)轉(zhuǎn)換，對(duì)于一個(gè)通用的坐標(biāo)轉(zhuǎn)換例子：

我們從

域變換到 域，使用微分的鏈?zhǔn)椒▌t：

則有如下關(guān)系：

式中的

代表雅可比矩陣行列式，其數(shù)值與參與變換的三角形節(jié)點(diǎn)坐標(biāo)有關(guān)：

可以發(fā)現(xiàn)，通過上面的參考單元，我們就能夠找到一種通用的方法來計(jì)算每個(gè)單元的質(zhì)量矩陣和剛度矩陣，并且積分區(qū)域都統(tǒng)一變成了一樣的。

（1）質(zhì)量矩陣

首先，對(duì)于每個(gè)單元的質(zhì)量矩陣，其中的元素計(jì)算公式變換為：

這里由于我們已經(jīng)知道了參考單元每一個(gè)節(jié)點(diǎn)的基函數(shù)表達(dá)式，我們可以把雅克比行列式從積分中提出來，將關(guān)于基函數(shù)的積分計(jì)算出來，即：

這個(gè)公式說明什么，說明經(jīng)過參考單元變換后，我們剛剛提到的質(zhì)量矩陣其實(shí)計(jì)算起來只需要求

就可以了，形象一點(diǎn)就是：

參考單元下計(jì)算單元質(zhì)量矩陣

細(xì)心的你可能會(huì)注意到，這里的編號(hào)統(tǒng)一都是

，和前面的單元 對(duì)應(yīng)關(guān)系為： ，按照逆時(shí)針方向一一對(duì)應(yīng)。

（2）剛度矩陣

剛度矩陣要稍微復(fù)雜一點(diǎn)點(diǎn)，因?yàn)檫@里面涉及到對(duì)基函數(shù)的偏導(dǎo)，將這個(gè)偏導(dǎo)寫成向量形式：

使用偏微分的鏈?zhǔn)椒▌t：

將公式（31）代入（30）中，則有：

其中的矩陣

代表雅克比矩陣的逆。對(duì)參考單元的基函數(shù)求導(dǎo)：

最終得到單元每一個(gè)元素的剛度矩陣計(jì)算公式為：

為了簡化計(jì)算，我們發(fā)現(xiàn)積分公式里面沒有任何一項(xiàng)和

有關(guān)，即：

這相當(dāng)于告訴我們：

用這個(gè)公式就好計(jì)算多了。

2.3 算法實(shí)現(xiàn)

（1）首先生成網(wǎng)格，這里封裝為函數(shù)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()%% 參數(shù)設(shè)置Lx = 1000; %定義單元右邊界（左邊界為0，如果不是，可以平移為0）Ly = 1000;%定義單元上邊界N = 60;%分割的一個(gè)方向的單元數(shù)目numelx = N;%定義分割的x方向單元數(shù)目（按矩形計(jì)算）numely = N;%定義分割的y方向單元數(shù)目（按矩形計(jì)算）numnodx = numelx + 1; % x方向節(jié)點(diǎn)個(gè)數(shù)比單元個(gè)數(shù)多1numnody = numely + 1; % y方向節(jié)點(diǎn)個(gè)數(shù)比單元個(gè)數(shù)多1nel = 3;%每個(gè)單元的節(jié)點(diǎn)數(shù)目,即每個(gè)單元上有幾個(gè)形函數(shù)參與作用，單元自由度coordx = linspace(0,Lx,numnodx)'; %等分節(jié)點(diǎn)的坐標(biāo)（為了方便，我這里采用等分的方式，事實(shí)上單元長度可以不一致，非均勻網(wǎng)格）coordy = linspace(0,Ly,numnody)'; %等分節(jié)點(diǎn)的坐標(biāo)（為了方便，我這里采用等分的方式，事實(shí)上單元長度可以不一致）[X, Y] = meshgrid(coordx,coordy);%張成網(wǎng)格，X和Y分別表示對(duì)應(yīng)位置的橫縱坐標(biāo)X = X';Y = Y';p = [X(:) Y(:)];%把網(wǎng)格一行一行扯開，coord的每一行是對(duì)應(yīng)節(jié)點(diǎn)的坐標(biāo)，按順序排p = p';connect = connect_mat(numnodx,numnody,nel);%連接矩陣，表示每個(gè)單元周圍的節(jié)點(diǎn)編號(hào)，也就是涉及的形函數(shù)編號(hào)t = connect';e = [1:numnodx numnodx*numely+1:numnodx*numely+numnodx ...numnodx+1:numnodx:numnodx*(numely-1)+1 2*numnodx:numnodx:numely*numnodx]; % 強(qiáng)制性邊界點(diǎn)的編號(hào)，本例子中是四條邊，下上左右邊endfunction connect_mat = connect_mat( numnodx,numnody,nel)%輸入橫縱坐標(biāo)的節(jié)點(diǎn)數(shù)目，和單元自由度%輸出連接矩陣，每個(gè)單元涉及的節(jié)點(diǎn)的編號(hào)xn = 1:(numnodx*numnody);%拉成一條編號(hào)A = reshape(xn,numnodx,numnody);%同形狀編號(hào)for i = 1:(numnodx-1)*(numnody-1)xg = rem(i,numnodx-1);%xg表示單元為左邊界數(shù)起第幾個(gè)if xg == 0xg = numnodx-1;endyg = ceil(i/(numnodx-1));%下邊界其數(shù)第幾個(gè)a = A(xg:xg+1,yg:yg+1);%這個(gè)小矩陣，拉直了就是連接矩陣a_vec = a(:);connect_mat(2*i-1:2*i,1:nel) = [a_vec([1 4 3])';a_vec([4 1 2])'];endend end

函數(shù)里面除了可以實(shí)現(xiàn)非結(jié)構(gòu)化網(wǎng)格，還提供了了等腰直角三角形形式的結(jié)構(gòu)化網(wǎng)絡(luò)可選。

（2）接著編寫計(jì)算質(zhì)量矩陣和剛度矩陣的函數(shù)FEM_2D_func()：

function FEM_2D_function = FEM_2D_func()%設(shè)置所有需要的函數(shù)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);%處理邊界條件%% 實(shí)現(xiàn)——組裝剛度矩陣方法一function A = assemble_matrix_2D(p, t)fprintf('組裝剛度矩陣:n');%獲取單元個(gè)數(shù)和節(jié)點(diǎn)個(gè)數(shù)number_of_elements = length(t);number_of_nodes = length(p);%初始化總剛度矩陣A = zeros(number_of_nodes, number_of_nodes);for n = 1: number_of_elements% 獲取單元局部編號(hào)到全局編號(hào)之間的對(duì)應(yīng)關(guān)系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');%獲取單元個(gè)數(shù)和節(jié)點(diǎn)個(gè)數(shù)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% 獲取單元局部編號(hào)到全局編號(hào)之間的對(duì)應(yīng)關(guān)系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 %% 組裝質(zhì)量矩陣M和右端向量Ffunction [M,F] = assemble_vector_2D(p, t)fprintf('組裝質(zhì)量矩陣:n');number_of_elements = length(t);number_of_nodes = length(p);M = zeros(number_of_nodes, number_of_nodes); % 總質(zhì)量矩陣Me = zeros(3, 3);F = zeros(number_of_nodes, 1); % 初始化右端向量Fe = zeros(3, 1); % 初始單元右端向量% 單元質(zhì)量矩陣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);%獲取當(dāng)前單元包含的節(jié)點(diǎn)編號(hào)vertices = p(:, local2global);%獲取當(dāng)前單元所有節(jié)點(diǎn)的x,y坐標(biāo)% 從參考單元到物理單元的映射% 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('質(zhì)量矩陣組裝完成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

這個(gè)函數(shù)里面還提供了另外一種求解剛度矩陣的方法可選（沒有使用到參考單元），并求解右端向量F、設(shè)置邊界條件(可選)。

（3）接著編寫主函數(shù)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'——可以顯示節(jié)點(diǎn)編號(hào)，'ElementLabels','on'——顯示單元編號(hào) pdemesh(p, e, t,'NodeLabels','on','ElementLabels','on');%生成結(jié)構(gòu)化三角網(wǎng)格 % 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) % 在點(diǎn)矩陣p中，第一行和第二行包含網(wǎng)格中點(diǎn)的x和y坐標(biāo)。 % 在邊矩陣e中，第一行和第二行包含起始點(diǎn)和結(jié)束點(diǎn)的索引，第三和第四行包含起始和結(jié)束參數(shù)值， %第五行包含邊緣段編號(hào)，第六和第七行包含左側(cè)和右側(cè)子域編號(hào)。 % 在三角形矩陣t中，前三行包含角點(diǎn)的索引，按逆時(shí)針順序給出，第四行包含子域編號(hào)。%% 參數(shù)設(shè)置 nt = 150; dt = 0.003; v = 1500; T = (1:nt)*dt; number_of_notes = length(p); fmain = 40;%% 設(shè)置震源 s_t = (1-2*pi^2*fmain*(T-0.2).^2).*exp(-fmain*pi^2*(T-0.2).^2);%雷克子波%% 調(diào)用函數(shù)計(jì)算矩陣 FEM = FEM_2D_func(); S = FEM.assemble_matrix_2D(p, t);%剛度矩陣 [M,F]= FEM.assemble_vector_2D(p, t);%質(zhì)量矩陣和右端向量 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);%% 利用遞推關(guān)系求波場值 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')

主函數(shù)里面可選這使用非結(jié)構(gòu)化和結(jié)構(gòu)化網(wǎng)格，并有隱式差分格式和顯式差分格式可選。

2.4 模擬結(jié)果

隱式差分格式模擬結(jié)果

顯式差分格式模擬結(jié)果

兩種都存在一定的數(shù)值耗散和色散現(xiàn)象。

3 聲波方程有限單元法數(shù)值模擬總結(jié)

通過這兩篇文章，相信大家對(duì)于有限單元法都有了一定的了解，我本來還想向大家介紹一下有關(guān)有限元法穩(wěn)定性分析的，但是本來篇幅就太長了，以后如果感興趣的人比較多，我再更新這一塊吧！

另外，我之前都沒有在文章中求贊求關(guān)注啥的，不過這次我還是厚臉皮一下，如果你喜歡這篇文章，覺得對(duì)你有用的話，請(qǐng)點(diǎn)擊一下喜歡，歡迎你關(guān)注我，我會(huì)在后續(xù)繼續(xù)分享有限體積法以及我所做過的一些比較有意思的數(shù)值模擬。

ps:文章中如有紕漏，歡迎在評(píng)論區(qū)指出。

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。