有限元的MATLAB解法.docx
- 文档编号:10425854
- 上传时间:2023-02-11
- 格式:DOCX
- 页数:8
- 大小:58.14KB
有限元的MATLAB解法.docx
《有限元的MATLAB解法.docx》由会员分享,可在线阅读,更多相关《有限元的MATLAB解法.docx(8页珍藏版)》请在冰豆网上搜索。
有限元的MATLAB解法
有限元的MATLAB解法
1.打开MATLAB。
2.输入“pdetool”再回车,会跳出PDEToolbox的窗口(PDE意为偏微分方程,是partialdifferentialequations的缩写),需要的话可点击Options菜单下Grid命令,打开栅格。
3.完成平面几何模型:
在PDEToolbox的窗口中,点击工具栏下的矩形几何模型进行制作模型,可画矩形R,椭圆E,圆C,然后在Setformula栏进行编辑并(如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标)
用算术运算符将图形对象名称连接起来,若还需要,可进行储存,形成M文件。
4.用左键双击矩形进行坐标设置:
将大的矩形left和bottom都设为0,width是矩形波导的X轴的长度,height是矩形波导的y轴的长度,以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。
5.进行边界设置:
点击“Boundary”中的“BoundaryMode”,再点击“Boundary”中的“SpecifyBoundaryConditions”,选择符合的边界条件,Neumann为诺曼条件,Dirichlet为狄利克雷条件,边界颜色显示为红色。
6.进入PDE模式:
点击"PDE"菜单下“PDEMode”命令,进入PDE模式,单击“PDESpecification”,设置方程类型,“Elliptic”为椭圆型,“Parabolic”为抛物型,“Hyperbolic”为双曲型,“Eigenmodes”为特征值问题。
7.对模型进行剖分:
点击“Mesh”中“InitializeMesh”进行初次剖分,若要剖的更细,再点击“RefineMesh”进行网格加密。
8.进行计算:
点击“Solve”中“SolvePDE”,解偏微分方程并显示图形解,u值即为Hz或者Ez。
9.单击“Plot”菜单下“Parameters”选项,打开“PlotSelection”对话框。
选中Color,Height(3-Dplot)和Showmesh三项,然后单击“Plot”按钮,显示三维图形解。
10.如果要画等值线图和矢量场图,单击“Plot”菜单下“Parameters”选项,打开“PlotSelection”对话框。
选中Contour和Arrows两项,然后单击Plot按钮,可显示解的等值线图和矢量场图。
11.将计算结果条件和边界导入MATLAB中:
点击“ExportSolution”,再点击“Mesh”中“ExportMesh”。
12.在MATLAB中将编好的计算程序导入,按F5运行。
备注:
Property(属性)用于画图时选用相应的绘图类型
u方程的解
abs(grad(u))每个三角形的中心的▽u的绝对值
abs(c*grad(u))每个三角形的中心的c·▽u的绝对值
-grad(u)u的负梯度-▽u
我们也可以用MATLAB程序求解PDE问题,同时显示解的图形;
一个长直接接地金属矩形槽,其侧壁与底面电位均为0,顶盖电位为100V,求槽内的电位分布:
(1)画出剖分图(尺寸与书上一样);
(2)标出各剖分点坐标值;
(3)求出各点电位值(用有限差分);
(4)画出等电位图。
解:
(1)编写以下程序得:
x=0:
5
y=0:
5
[X,Y]=meshgrid(x,y)
plot(X,Y)
holdon
plot(Y,X)
fori=0:
5
s=i:
5
t=0:
(5-i)
plot(s,t)
plot(t,s)
end
得到剖分图如下:
(2)用有限元法编写程序如下:
Nx=6;Ny=6;Xm=5;Ym=15;Np=5;Nq=5;
fori=1:
Nx
forj=1:
Ny
N(i,j)=(i-1)*Ny+j;/i列j行的节点编号/
X(N(i,j))=(i-1)*Xm/Np;/节点横坐标/
Y(N(i,j))=(j-1)*Ym/Nq;/节点纵坐标/
end
end
fori=1:
2*Xm
forj=1:
Ym
ifrem(i,2)==1
L(i,j)=(i-1)*Nq+j;
p(i,j)=2*(i-1)*Ny/2+Ny+j+1;
q(i,j)=p(i,j)-Ny;
r(i,j)=q(i,j)-1;
elserem(i,2)==0
L(i,j)=(i-1)*Ny+j;
p(i,j)=(2i-2)*Ny/2+j;
q(i,j)=p(i,j)+Ny;
r(i,j)=q(i,j)+1;
end
end
end
fori=1:
2*Xm
forj=1:
Ym
b(p(i,j))=Y(q(i,j))-Y(r(i,j));b(q(i,j))=Y(r(i,j))-Y(p(i,j));
b(r(i,j))=Y(p(i,j))-Y(q(i,j));c(p(i,j))=X(r(i,j))-X(q(i,j));
c(q(i,j))=X(p(i,j))-X(r(i,j));c(r(i,j))=X(q(i,j))-X(p(i,j));
area(i,j)=(b(p(i,j))*c(q(i,j))-b(q(i,j))*c(p(i,j)))/2;
K=zeros(Nx*Ny);
Kpp(i,j)=(b(p(i,j))^2+c(p(i,j))^2)/(2*area(i,j));
Kpq(i,j)=(b(p(i,j))*b(q(i,j))+c(p(i,j))*c(q(i,j)))/(2*area(i,j));
Kpr(i,j)=(b(p(i,j))*b(r(i,j))+c(p(i,j))*c(r(i,j)))/(2*area(i,j));
Kqp(i,j)=Kpq(i,j);Kqq(i,j)=(b(q(i,j))^2+c(q(i,j))^2)/(2*area(i,j));
Kqr(i,j)=(b(q(i,j))*b(r(i,j))+c(q(i,j))*c(r(i,j)))/(2*area(i,j));
Krp(i,j)=Kpr(i,j);
Krq(i,j)=Kqr(i,j);Krr(i,j)=(b(r(i,j))^2+c(r(i,j))^2)/(2*area(i,j));
end
end
fori=1:
2*Xm
forj=1:
Ym
K(p(i,j),p(i,j))=Kpp(i,j)+K(p(i,j),p(i,j));
K(p(i,j),q(i,j))=Kpq(i,j)+K(p(i,j),q(i,j));
K(p(i,j),r(i,j))=Kpr(i,j)+K(p(i,j),r(i,j));
K(q(i,j),p(i,j))=Kqp(i,j)+K(q(i,j),p(i,j));
K(q(i,j),q(i,j))=Kqq(i,j)+K(q(i,j),q(i,j));
K(q(i,j),r(i,j))=Kqr(i,j)+K(q(i,j),r(i,j));
K(r(i,j),p(i,j))=Krp(i,j)+K(r(i,j),p(i,j));
K(r(i,j),q(i,j))=Krq(i,j)+K(r(i,j),q(i,j));
K(r(i,j),r(i,j))=Krr(i,j)+K(r(i,j),r(i,j));
end
end
fori=1:
11
K(i,:
)=0;
K(i,i)=1;
end
fori=1:
11:
111
K(i,:
)=0;
K(i,i)=1;
end
fori=111:
121
K(i,:
)=0;
K(i,i)=1;
end
fori=11:
11:
121
K(i,:
)=0;
K(i,i)=1;
end
B=zeros(121,1);
fori=11:
11:
121
B(i,1)=100;
end
U=K\B;
b=1;XX=zeros(11,11)
forj=1:
11
fori=1:
11
XX(i,j)=U(b,1);
b=b+1;
end
end
subplot(1,2,1),mesh(XX)
axis([0,11,0,11,0,100])
subplot(1,2,2),contour(XX,15)
holdon
(3)由上面的程序得到节点电位:
(4)由程序得到的电场分布图及等位线图如下:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 MATLAB 解法