有限元的MATLAB解法Word下载.docx
《有限元的MATLAB解法Word下载.docx》由会员分享，可在线阅读，更多相关《有限元的MATLAB解法Word下载.docx（11页珍藏版）》请在冰豆网上搜索。
点击“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”。
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）编写以下程序得：
x=0:
5
y=0:
[X,Y]=meshgrid（x,y）
plot（X,Y）
holdon
plot（Y,X）
fori=0:
s=i:
t=0:
（5-i）
plot（s,t）
plot（t,s）
end
得到剖分图如下：
（2）用有限元法编写程序如下：
Nx=6;
Ny=6;
Xm=5;
Ym=15;
Np=5;
Nq=5;
fori=1:
Nx
forj=1:
Ny
N（i,j）=（i-1）*Ny+j;
/i列j行的节点编号/
X（N（i,j））=（i-1）*Xm/Np;
/节点横坐标/
Y（N（i,j））=（j-1）*Ym/Nq;
/节点纵坐标/
end
2*Xm
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
2*Xm
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））=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））;
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））;
11
K（i,:
）=0;
K（i,i）=1;
11:
111
fori=111:
121
fori=11:
B=zeros（121,1）;
B（i,1）=100;
U=K\B;
b=1;
XX=zeros（11,11）
forj=1:
fori=1:
XX（i,j）=U（b,1）;
b=b+1;
subplot（1,2,1）,mesh（XX）
axis（[0,11,0,11,0,100]）
subplot（1,2,2）,contour（XX,15）
（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）根据时域有限差分编写的程序可画出主模的电场分布图如下：
（2）根据时域有限差分编写的程序可画出频谱图和场结构图，从左图中可以读出主模截止频率
值，主模
，根据
，其中
从而计算出主模截至波长
。
升级会员 