当前位置:文档之家› 有限元的MATLAB解法

有限元的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改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标)用算术运算符将图形对象名称连接起来,若还需要,可进行储存,形成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 Specification”,设置方程类型,“Elliptic”为椭圆型,“Parabolic”为抛物型,“Hyperbolic”为双曲型,“Eigenmodes”为特征值问题。

7.对模型进行剖分:点击“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”对话框。

选中Contour和Arrows两项,然后单击Plot按钮,可显示解的等值线图和矢量场图。

11.将计算结果条件和边界导入MATLAB中:点击“Export Solution”,再点击“Mesh”中“Export Mesh”。

12.在MATLAB中将编好的计算程序导入,按F5运行。

备注:Property(属性)用于画图时选用相应的绘图类型u 方程的解abs(grad(u)) 每个三角形的中心的▽u的绝对值abs(c*grad(u)) 每个三角形的中心的c·▽u的绝对值- grad(u) u的负梯度-▽u我们也可以用MATLAB程序求解PDE问题,同时显示解的图形;一个长直接接地金属矩形槽,其侧壁与底面电位均为0,顶盖电位为100V,求槽内的电位分布:100V0V0V0V(1)画出剖分图(尺寸与书上一样);(2)标出各剖分点坐标值;(3)求出各点电位值(用有限差分);(4)画出等电位图。

解:(1)编写以下程序得:x=0:5y=0:5[X,Y]=meshgrid(x,y)plot(X,Y)hold onplot(Y,X)for i=0:5s=i:5t=0:(5-i)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;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(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)); 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(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,:)=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=K\B;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)hold on(3)由上面的程序得到节点电位:V1=0 0 0 0 00 7.1429 9.8214 7.1429 00 18.7500 25.0000 18.7500 00 42.8571 52.6786 42.8571 00 100.0000 100.0000 100.0000 0(4)由程序得到的电场分布图及等位线图如下:4.用有限元法求矩形波导(b/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,截止波长图如下:单模TE10区λ时截止模区λ(3)在TM模式下设置边界条件为Dirichlet条件,使用编制好的程序计算出主模的截止波长为0.8179a,第一高次模为0.6655a ,第二高次模为0.5315a ,截止波长图如下:高次模区单模TM 11区截止模区λ5.用时域有限差分求解上述4题中的前两问。

解:(1)根据时域有限差分编写的程序可画出主模的电场分布图如下:在Neumann 边界条件下:在Dirichlet 边界条件下:(2)根据时域有限差分编写的程序可画出频谱图和场结构图,从左图中可以读出主模截止频率c f 值,主模81.499610c f Z=碒,根据1/(c f l =,其中1208.8510e -=?,-60 1.256630610,m =?从而计算出主模截至波长 2.0273c a l =。

.。

相关主题