有限元的MATLAB解法范本模板.docx

1、有限元的MATLAB解法范本模板有限元的MATLAB解法1.打开MATLAB。2。输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写）,需要的话可点击Options菜单下Grid命令，打开栅格。3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型,可画矩形R，椭圆E,圆C,然后在Set formula栏进行编辑并（如双脊波导R1+R2+R3改为RIR2-R3，设定a、b、s/a、d/b的值从而方便下步设定坐标）用算术运算符将图形对象名称连接起来，若

2、还需要，可进行储存,形成M文件。4.用左键双击矩形进行坐标设置:将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。6.进入PDE模式：点击”PDE菜单下“PDE Mode命令,进入PDE模式，单击“PDE Specif

3、ication”,设置方程类型，“Elliptic为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。7。对模型进行剖分：点击“Mesh中“Initialize Mesh进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。8.进行计算:点击“Solve中“Solve PDE”，解偏微分方程并显示图形解，u值即为Hz或者Ez。9。单击“Plot”菜单下“Parameters”选项，打开“Plot Selection”对话框.选中Color,Height（3D plot)和Show mesh三项,然后单击“Plot

4、按钮，显示三维图形解.10。如果要画等值线图和矢量场图，单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Contour和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。11。将计算结果条件和边界导入MATLAB中：点击“Export Solution,再点击“Mesh中“Export Mesh”。12。在MATLAB中将编好的计算程序导入，按F5运行。备注：Property（属性）用于画图时选用相应的绘图类型u 方程的解abs（grad（u) 每个三角形的中心的u的绝对值abs（cgrad(u) 每个三角形的中心的cu的

5、绝对值 grad（u) u的负梯度u我们也可以用MATLAB程序求解PDE问题，同时显示解的图形；一个长直接接地金属矩形槽,其侧壁与底面电位均为0，顶盖电位为100V，求槽内的电位分布：（1）画出剖分图（尺寸与书上一样）;（2）标出各剖分点坐标值；（3）求出各点电位值(用有限差分）;（4）画出等电位图。解:(1）编写以下程序得：x=0:5y=0：5X，Y=meshgrid(x，y)plot(X，Y)hold onplot(Y，X)for i=0：5 s=i:5 t=0:（5-i) plot（s，t） plot(t,s）end得到剖分图如下：(2)用有限元法编写程序如下：Nx=6;Ny=6；Xm

6、=5；Ym=15;Np=5；Nq=5;for i=1:Nx for j=1：Ny N(i,j）=（i1）Ny+j； /i列j行的节点编号/ X（N（i,j)）=（i-1）Xm/Np；/节点横坐标/ Y(N(i,j)）=（j1）*Ym/Nq；/节点纵坐标/ endendfor i=1：2Xm 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)=（2i2)*Ny/

7、2+j; q（i，j）=p(i,j)+Ny; r(i，j）=q(i，j）+1； end endendfor i=1:2*Xm for j=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=zer

8、os（NxNy); 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）)）/（2area(

9、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）); endendfor i=1:2*Xm for j=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(

10、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）; endendfor i=1:11 K（i，:)=0； K(i，i)=1;endfor i=1：11:111 K（i，：）=0； K（i,i）=1；endfor i=111：121 K(i,:)=0； K(

11、i，i）=1；endfor i=11:11:121 K(i，:)=0; K(i,i)=1;endB=zeros(121,1）；for i=11:11:121 B(i，1)=100;endU=KB;b=1；XX=zeros(11，11)for j=1：11 for i=1：11 XX（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)hold on（3）由上面的程序得到节点电位:(4)由程序得到的电场分布图及等位线图如下：4.用有限元法求矩形波导（b

12、/a=0.45）的： （1）电场分布图； (2）求TE模式下的主模、第一、二高次模的截止波长（5次），画出截至波长图； （3)求TM模式下的主模、第一、二高次模的截止波长(5次)，画出截至波长图。解:利用MATLAB中的PDE工具箱：取矩形波导的宽边尺寸为a，窄边尺寸为0。45a。（1）主模的电场分布图如下： 在Neumann边界条件下：在Dirichlet边界条件下：（2）在TE模式下设置边界条件为Neumann条件，使用编制好的程序计算出主模的截止波长为1.9988a,第一高次模为0.9977a，第二高次模为0。8972a，截止波长图如下: (3）在TM模式下设置边界条件为Dirichlet条件，使用编制好的程序计算出主模的截止波长为0.8179a,第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下：5.用时域有限差分求解上述4题中的前两问.解：（1）根据时域有限差分编写的程序可画出主模的电场分布图如下： 在Neumann边界条件下:在Dirichlet边界条件下：（2)根据时域有限差分编写的程序可画出频谱图和场结构图,从左图中可以读出主模截止频率值，主模,根据，其中, 从而计算出主模截至波长。