《有限元的MATLAB解法.docx》由会员分享，可在线阅读，更多相关《有限元的MATLAB解法.docx（9页珍藏版）》请在三一办公上搜索。
1、有限元的MATLAB解法有限元的MATLAB解法 1.打开MATLAB。 2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口，需要的话可点击Options菜单下Grid命令，打开栅格。 3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进行编辑并 用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。 4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角
2、点为原点坐标为参考设置其他矩形坐标。 5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击5 “Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。 6.进入PDE模式：点击PDE菜单下“PDE Mode”命令，进入PDE模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。 7.对模型进行剖
3、分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。 8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。 9.单击“Plot”菜单下“Parameters”选项，打开“Plot Selection”对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单击“Plot”按钮，显示三维图形解。 10.如果要画等值线图和矢量场图，单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Cont
4、our和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。 11.将计算结果条件和边界导入MATLAB中：点击“Export Solution”,再点击“Mesh”中“Export Mesh”。 6 12.在MATLAB中将编好的计算程序导入，按F5运行。 备注： Property用于画图时选用相应的绘图类型 u 方程的解 abs(grad(u) 每个三角形的中心的u的绝对值 abs(c*grad(u) 每个三角形的中心的cu的绝对值 - grad(u) u的负梯度-u 我们也可以用MATLAB程序求解PDE问题，同时显示解的图形； 7 一个长直接接地金属矩形槽，其侧壁与底
5、面电位均为0，顶盖电位为100V，求槽内的电位分布： 100V0V0V0V画出剖分图； 标出各剖分点坐标值； 求出各点电位值； 画出等电位图。 解：编写以下程序得： x=0:5 y=0:5 8 X,Y=meshgrid(x,y) plot(X,Y) hold on plot(Y,X) for i=0:5 s=i:5 t=0:(5-i) plot(s,t) plot(t,s) end 得到剖分图如下： 用有限元法编写程序如下： 9 Nx=6;Ny=6;Xm=5;Ym=15;Np=5;Nq=5; for i=1:Nx for j=1:Ny N(i,j)=(i-1)*Ny+j; /i列j行的节点编号
6、/ X(N(i,j)=(i-1)*Xm/Np;/节点横坐标/ Y(N(i,j)=(j-1)*Ym/Nq;/节点纵坐标/ end end for i=1:2*Xm for j=1:Ym if rem(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; else rem(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; 10 r(i,j)=q(i,j)+1; end end end for i
7、=1:2*Xm for j=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
8、,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); 11 Krq(i,j)=Kqr(i,j); Krr(i,j)=(b(r(i,j)2+
9、c(r(i,j)2)/(2*area(i,j); end end for i=1:2*Xm for j=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)
10、,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 for i=1:11 K(i,:)=0; 12 K(i,i)=1; end for i=1:11:111 K(i,:)=0; K(i,i)=1; end for i=111:121 K(i,:)=0; K(i,i)=1; end for i=11:11:121 K(i,:)=0; K(i,i)=1; end B=
11、zeros(121,1); for i=11:11:121 B(i,1)=100; end U=KB; 13 b=1;XX=zeros(11,11) for j=1:11 for i=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) hold on 由上面的程序得到节点电位： V1=0 0 0 0 00 7.1429 9.8214 7.1429 00 18.7500 25.0000 18.7500 00 42.8571 5
12、2.6786 42.8571 00 100.0000 100.0000 100.0000 014 由程序得到的电场分布图及等位线图如下： 4.用有限元法求矩形波导的： 电场分布图； 求TE模式下的主模、第一、二高次模的截止波长，画出截至波长图； 求TM模式下的主模、第一、二高次模的截止波长，画出截至波长图。 解：利用MATLAB中的PDE工具箱：取矩形波导的宽边尺寸为a，窄边尺寸为0.45a。 15 主模的电场分布图如下： 在Neumann边界条件下： 在Dirichlet边界条件下： 在TE模式下设置边界条件为Neumann条件，使用编制好6 的程序计算出主模的截止波长为1.9988a，第一
13、高次模为0.9977a，第二高次模为0.8972a，截止波长图如下： TE20TE01TE21TE11TE300当0.9977a时高次模区a当0.9977a1.9988a时单模TE10区2a当1.9988a时截止模区TE10在TM模式下设置边界条件为Dirichlet条件，使用编制好的程序计算出主模的截止波长为0.8179a，第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下： TM11TM21TM12TM13TM310当0.6655a时高次模区a当0.6655a0.8179a时单模TM11区2a当0.8179a时截止模区5.用时域有限差分求解上述4题中的前两问。 解：根据时域有限差分编写的程序可画出主模的电场分布7 图如下： 在Neumann边界条件下： 在Dirichlet边界条件下： 根据时域有限差分编写的程序可画出频谱图和场结构图，108从左图中可以读出主模截止频率fc值，主模fc=1.4996碒-12据lc=1/(fcmoe0)，其中e0=8.85?10, Z，根m0=1.2566306?10-6,从而计算出主模截至波长lc=2.0273a。 8 9
