成长厉程网站,永兴做网站,跨境电商导购网站建设,网站目录管理模板下载2.1 前言承接上一篇文章#xff0c;前面我们已经介绍了一维声波方程有限元求解#xff1a;蓝不是蓝#xff1a;有限单元法(Finite Element Method)实现声波方程模拟#xff08;Part 1#xff09;zhuanlan.zhihu.com这一部分将一维问题提升到二维问题。不知道大家有没有发…2.1 前言承接上一篇文章前面我们已经介绍了一维声波方程有限元求解蓝不是蓝有限单元法(Finite Element Method)实现声波方程模拟Part 1zhuanlan.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如果不是可以平移为0Ly 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*numely1:numnodx*numelynumnodx ...numnodx1: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:xg1,yg:yg1);%这个小矩阵拉直了就是连接矩阵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 m1:3
% for n1: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,:))0abs(p(1,:))50abs(p(2,:))0abs(p(2,:))50);%% 利用递推关系求波场值
for i 2:nt-1U(source_x(1),i) s_t(i);%隐式
% U(:,i1) (Mv^2*dt^2*S)(M*(2*U(:,i)-U(:,i-1)));
% %显式right_vector v^2*dt^2*S*U(:,i);U(:,i1) 2*U(:,i)-U(:,i-1)-Mright_vector;U(boundarynodes,i1) 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 i5: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:文章中如有纰漏欢迎在评论区指出。
