有限元地MATLAB解法

《有限元地MATLAB解法》由会员分享，可在线阅读，更多相关《有限元地MATLAB解法（17页珍藏版）》请在装配图网上搜索。

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

2、起来，若还需要，可进行储存， 形成M文件。III 3|I12 dIII I bII s II c1=c2Ic1_|_c2_ |4. 用左键双击矩形进行坐标设置：将大的矩形left和bottom都设 为0， width是矩形波导的X轴的长度，height是矩形波导的y轴 的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标5. 进行边界设置：点击“Boundary”中的“Boundary Mode”，再点 击 “Boundary” 中的“Specify Boundary Conditions”，选择符合 的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边 界颜色显示

3、为红色。6. 进入PDE模式：点击PDE菜单下“PDEMode ”命令，进入PDE模 式，单击“PDESpecification”，设置方程类型，“Elliptic”为椭 圆型，“ Parabolic为抛物型，“ Hyperbolic为双曲型，“ Eige nmodes” 为特征值问题。7. 对模型进行剖分：点击“Mesh”中“InitializeMesh”进行初次 剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。8. 进行计算：点击“Solve”中“Solve PDE，解偏微分方程并显示 图形解， u 值即为 Hz 或者 Ez。9. 单击 “Plot” 菜单下 “Param

4、eters” 选项，打开 “Plot Selection” 对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单 击“Plot”按钮，显示三维图形解。10. 如果要画等值线图和矢量场图，单击Plot菜单下“Parameters” 选项，打开“Plot Selection”对话框。选中Contour和Arrows两 项，然后单击Plot按钮，可显示解的等值线图和矢量场图。11. 将计算结果条件和边界导入MATLAB中：点击“Export Solution”, 再点击 “Mesh”中 “Export Mesh”。12.在 MATLAB 中将编好的计算程序导入，按

5、 F5 运行。备注：Property （属性）用于画图时选用相应的绘图类型方程的解abs(grad(u)每个三角形的中心的u的绝对值abs(c*grad(u)每个三角形的中心的C u的绝对值- grad(u)u 的负梯度 - u我们也可以用 MATLAB 程序求解 PDE 问题，同时显示解的图形；4 4 二=inima2h （ f circlr1 , r, 1 ；初卅化网珞error 二;err-1:亡rr0.021,r t =ref inemesh （ J circl&gp ,p, c . r. U 功口犠网格：匚 $ H亠尸T . d二：遇求解：）,-2-p（2;: j m 牛定义将确解F

6、r二nnTTT （* inf J :；cm二；ci rDr err-；君h （p, j r. j举绘制同旃图 passer； （R匕” 鲁绘制耕前曲面因 pdes2rt （pT , u-exoct） P象制许拦图形一个长直接接地金属矩形槽，其侧壁与底面电位均为0，顶盖电 位为1OOV，求槽内的电位分布：（1 ）画出剖分图（尺寸与书上一样）；（2）标出各剖分点坐标值；（3）求出各点电位值（用有限差分）（4）画出等电位图。解：（1 ）编写以下程序得：x=0:5y=0:5X,Y=meshgrid(x,y)plot(X,Y)hold onplot(Y,X)for i=0:5s=i:5t=0:(5-i)

7、plot(s,t)plot(t,s)end得到剖分图如下：2）用有限元法编写程序如下：Nx=6;Ny=6;Xm=5;Ym=15;Np=5;Nq=5;for i=1:Nxfor j=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;/ 节点纵坐标 /endendfor i=1:2*Xmfor j=1:Ymif rem(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

8、;r(i,j)=q(i,j)-1;else rem(i,2)=0L(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;endendendfor i=1:2*Xmfor 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(

9、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

10、);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);endendfor i=1:2*Xmfor 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(

11、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);endendfor i=1:11K(i,:)=0;K(i,i)=1;endfor i=1:11:111K(i,:)=

12、0;K(i,i)=1;endfor i=111:121K(i,:)=0;K(i,i)=1;endfor i=11:11:121K(i,:)=0;K(i,i)=1;endB=zeros(121,1);for i=11:11:121B(i,1)=100;endU=KB;b=1;XX=zeros(11,11)for j=1:11for i=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)V1=hold on0000007.14299.8214

13、7.14290018.750025.000018.75000042.857152.678642.857100100.0000100.0000100.000003）由上面的程序得到节点电位：4）由程序得到的电场分布图及等位线图如下：4. 用有限元法求矩形波导（b/a=0.45）的：（ 1 ）电场分布图；（2）求 TE 模式下的主模、第一、二高次模的截止波长（ 5 次），画出截至波长图；（3）求 TM 模式下的主模、第一、二高次模的截止波长（ 5 次），画出截至波长图。解：利用MATLAB中的PDE工具箱：取矩形波导的宽边尺寸为a, 窄边尺寸为 0.45a。（ 1）主模的电场分布图如下：在 Neu

14、mann 边界条件下：Lambda(2)=9.8811 Contour: u0.80.70.60.50.40.30.20.1-0.1-0.20.20.40.81.20.060.040.020-0.02-0.04-0.06在 Dirichlet 边界条件下：-0.01-0.02-0.03-0.04-0.05-0.06-0.07-0.08-0.09-0.1（2）在 TE 模式下设置边界条件为 Neumann 条件，使用编制好 的程序计算出主模的截止波长为1.9988a ，第一高次模为0.9977a ， 第二高 次模为 0.8972a ， 截 止 波长图 如 下：TETETEoJTEiiTE30峙1

15、12a当乎0.997时L当0.9977aX1.9988am 当 1-9988a 时高次模区单模TE 0区截止模区（3）在 TM 模式下设置边界条件为 Dirichlet 条件，使用编制好的程序计算出主模的截止波长为0.8179a，第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下:0a2a当九v0.6655a时高次模区当0.8179av九时截止模区当 0.6655av 九 v0.8179a 时单模TM 1区5. 用时域有限差分求解上述 4 题中的前两问。解：（1 ）根据时域有限差分编写的程序可画出主模的电场分布图如下：在 Dirichlet 边界条件下：在 Neumann 边界条件下： (2)根据时域有限差分编写的程序可画出频谱图和场结构图， 从左图中可以读出主模截止频率 f 值，主模 f = 1.4996 碒108 ，根 c c Z据l = 1/( f Jm e )，其中 e = 8.85? 10-12,m = 1.2566306? 10-6,从而计算出主模截至波长/ = 2.0273a。X 1040 000.511.522.5频率（Hz）x1QsxIO-65 2 5 1 5 2.10.sAScnv