有限元的MATLAB解法

5有限元的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”,再点击6“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”。712.在MATLAB中将编好的计算程序导入，按F5运行。备注：Property（属性）用于画图时选用相应的绘图类型u方程的解abs(grad(u))每个三角形的中心的▽u的绝对值abs(c*grad(u))每个三角形的中心的c·▽u的绝对值-grad(u)u的负梯度-▽u我们也可以用MATLAB程序求解PDE问题，同时显示解的图形；8一个长直接接地金属矩形槽，其侧壁与底面电位均为0，顶盖电位为100V，求槽内的电位分布：100V0V0V0V（1）画出剖分图（尺寸与书上一样）；（2）标出各剖分点坐标值；（3）求出各点电位值（用有限差分）；（4）画出等电位图。解：（1）编写以下程序得：x=0:5y=0:59[X,Y]=meshgrid(x,y)plot(X,Y)holdonplot(Y,X)fori=0:5s=i:5t=0:(5-i)plot(s,t)plot(t,s)end得到剖分图如下：（2）用有限元法编写程序如下：10Nx=6;Ny=6;Xm=5;Ym=15;Np=5;Nq=5;fori=1:Nxforj=1:NyN(i,j)=(i-1)*Ny+j;/i列j行的节点编号/X(N(i,j))=(i-1)*Xm/Np;/节点横坐标/Y(N(i,j))=(j-1)*Ym/Nq;/节点纵坐标/endendfori=1:2*Xmforj=1:Ymifrem(i,2)==1L(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)==0L(i,j)=(i-1)*Ny+j;p(i,j)=(2i-2)*Ny/2+j;q(i,j)=p(i,j)+Ny;11r(i,j)=q(i,j)+1;endendendfori=1:2*Xmforj=1:Ymb(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);12Krq(i,j)=Kqr(i,j);Krr(i,j)=(b(r(i,j))^2+c(r(i,j))^2)/(2*area(i,j));endendfori=1:2*Xmforj=1:YmK(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));endendfori=1:11K(i,:)=0;13K(i,i)=1;endfori=1:11:111K(i,:)=0;K(i,i)=1;endfori=111:121K(i,:)=0;K(i,i)=1;endfori=11:11:121K(i,:)=0;K(i,i)=1;endB=zeros(121,1);fori=11:11:121B(i,1)=100;endU=K\B;14b=1;XX=zeros(11,11)forj=1:11fori=1:11XX(i,j)=U(b,1);b=b+1;endendsubplot(1,2,1),mesh(XX)axis([0,11,0,11,0,100])subplot(1,2,2),contour(XX,15)holdon（3）由上面的程序得到节点电位：V1=0000007.14299.82147.14290018.750025.000018.75000042.857152.678642.857100100.0000100.0000100.0000015（4）由程序得到的电场分布图及等位线图如下：4.用有限元法求矩形波导（b/a=0.45）的：（1）电场分布图；（2）求TE模式下的主模、第一、二高次模的截止波长（5次），画出截至波长图；（3）求TM模式下的主模、第一、二高次模的截止波长（5次），画出截至波长图。解：利用MATLAB中的PDE工具箱：取矩形波导的宽边尺寸为a，窄边尺寸为0.45a。6（1）主模的电场分布图如下：在Neumann边界条件下：在Dirichlet边界条件下：（2）在TE模式下设置边界条件为Neumann条件，使用编制好7的程序计算出主模的截止波长为1.9988a，第一高次模为0.9977a，第二高次模为0.8972a，截止波长图如下：0a2a当λ0.9977a时高次模区当0.9977aλ1.9988a时单模TE10区当1.9988aλ时截止模区λTE10TE20TE01TE11TE21TE30（3）在TM模式下设置边界条件为Dirichlet条件，使用编制好的程序计算出主模的截止波长为0.8179a，第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下：0a2a当λ0.6655a时高次模区当0.6655aλ0.8179a时单模TM11区当0.8179aλ时截止模区λTM11TM21TM12TM31TM135.用时域有限差分求解上述4题中的前两问。解：（1）根据时域有限差分编写的程序可画出主模的电场分布8图如下：在Neumann边界条件下：在Dirichlet边界条件下：（2）根据时域有限差分编写的程序可画出频谱图和场结构图，从左图中可以读出主模截止频率cf值，主模81.499610cfZ=碒，根据01/()ccoflme=，其中1208.8510e-=?,-601.256630610,m=?从而计算出主模截至波长2.0273cal=。9