有限元的MATLAB解法.docx
《有限元的MATLAB解法.docx》由会员分享，可在线阅读，更多相关《有限元的MATLAB解法.docx（8页珍藏版）》请在冰豆网上搜索。
有限元的MATLAB解法
有限元的MATLAB解法
1.翻开MATLAB。
2.输入“pdetool〞再回车，会跳出PDEToolbo*的窗口〔PDE意为偏微分方程，是partialdifferentialequations的缩写〕，需要的话可点击Options菜单下Grid命令，翻开栅格。
3.完成平面几何模型：
在PDEToolbo*的窗口中，点击工具栏下的矩形几何模型进展制作模型，可画矩形R，椭圆E，圆C，然后在Setformula栏进展编辑并〔如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标〕
用算术运算符将图形对象名称连接起来，假设还需要，可进展储存，形成M文件。
4.用左键双击矩形进展坐标设置：
将大的矩形left和bottom都设为0，width是矩形波导的*轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。
5.进展边界设置：
点击“Boundary〞中的“BoundaryMode〞,再点击“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中：
点击“E*portSolution〞,再点击“Mesh〞中“E*portMesh〞。
12.在MATLAB中将编好的计算程序导入，按F5运行。
备注：
Property〔属性〕用于画图时选用相应的绘图类型
u方程的解
abs（grad（u））每个三角形的中心的▽u的绝对值
abs（c*grad（u））每个三角形的中心的c·▽u的绝对值
-grad（u）u的负梯度-▽u
我们也可以用MATLAB程序求解PDE问题，同时显示解的图形；
一个长直接接地金属矩形槽，其侧壁与底面电位均为0，顶盖电位为100V，求槽的电位分布：
〔1〕画出剖分图〔尺寸与书上一样〕；
〔2〕标出各剖分点坐标值；
〔3〕求出各点电位值〔用有限差分〕；
〔4〕画出等电位图。
解：
〔1〕编写以下程序得：
*=0:
5
y=0:
5
[*,Y]=meshgrid（*,y）
plot（*,Y）
holdon
plot（Y,*）
fori=0:
5
s=i:
5
t=0:
（5-i）
plot（s,t）
plot（t,s）
end
得到剖分图如下：
〔2〕用有限元法编写程序如下：
N*=6;Ny=6;*m=5;Ym=15;Np=5;Nq=5;
fori=1:
N*
forj=1:
Ny
N（i,j）=（i-1）*Ny+j;/i列j行的节点编号/
*（N（i,j））=（i-1）**m/Np;/节点横坐标/
Y（N（i,j））=（j-1）*Ym/Nq;/节点纵坐标/
end
end
fori=1:
2**m
forj=1:
Ym
ifrem（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;
elserem（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;
r（i,j）=q（i,j）+1;
end
end
end
fori=1:
2**m
forj=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））=*（r（i,j））-*（q（i,j））;
c（q（i,j））=*（p（i,j））-*（r（i,j））;c（r（i,j））=*（q（i,j））-*（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（N**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））;
end
end
fori=1:
2**m
forj=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）,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
fori=1:
11
K（i,:
）=0;
K（i,i）=1;
end
fori=1:
11:
111
K（i,:
）=0;
K（i,i）=1;
end
fori=111:
121
K（i,:
）=0;
K（i,i）=1;
end
fori=11:
11:
121
K（i,:
）=0;
K（i,i）=1;
end
B=zeros（121,1）;
fori=11:
11:
121
B（i,1）=100;
end
U=K\B;
b=1;**=zeros（11,11）
forj=1:
11
fori=1:
11
**（i,j）=U（b,1）;
b=b+1;
end
end
subplot（1,2,1）,mesh（**）
a*is（[0,11,0,11,0,100]）
subplot（1,2,2）,contour（**,15）
holdon
〔3〕由上面的程序得到节点电位：
〔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，截止波长图如下：
〔3〕在TM模式下设置边界条件为Dirichlet条件，使用编制好的程序计算出主模的截止波长为0.8179a，第一高次模为0.6655a，第二高次模为0.5315a，截止波长图如下：
5.用时域有限差分求解上述4题中的前两问。
解：
〔1〕根据时域有限差分编写的程序可画出主模的电场分布图如下：
在Neumann边界条件下：
在Dirichlet边界条件下：
〔2〕根据时域有限差分编写的程序可画出频谱图和场构造图，从左图中可以读出主模截止频率
值，主模
，根据
，其中
从而计算出主模截至波长
。
升级会员 