有限元的MATLAB解法范本模板.docx
《有限元的MATLAB解法范本模板.docx》由会员分享，可在线阅读，更多相关《有限元的MATLAB解法范本模板.docx（13页珍藏版）》请在冰豆网上搜索。
有限元的MATLAB解法范本模板
有限元的MATLAB解法
1.打开MATLAB。
2。
输入“pdetool”再回车，会跳出PDEToolbox的窗口（PDE意为偏微分方程，是partialdifferentialequations的缩写）,需要的话可点击Options菜单下Grid命令，打开栅格。
3.完成平面几何模型：
在PDEToolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型,可画矩形R，椭圆E,圆C,然后在Setformula栏进行编辑并（如双脊波导R1+R2+R3改为RI—R2-R3，设定a、b、s/a、d/b的值从而方便下步设定坐标）
用算术运算符将图形对象名称连接起来，若还需要，可进行储存,形成M文件。
4.用左键双击矩形进行坐标设置:
将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，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中：
点击“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：
5
[X，Y］=meshgrid（x，y）
plot（X，Y）
holdon
plot（Y，X）
fori=0：
5
s=i:
5
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
end
fori=1：
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
end
end
fori=1:
2*Xm
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））=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））;
end
end
fori=1:
2*Xm
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；XX=zeros（11，11）
forj=1：
11
fori=1：
11
XX（i，j）=U（b,1）；
b=b+1;
end
end
subplot（1,2，1），mesh（XX）
axis（[0,11，0,11,0，100］）
subplot（1,2，2）,contour（XX,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）根据时域有限差分编写的程序可画出频谱图和场结构图,从左图中可以读出主模截止频率
值，主模
根据
，其中
从而计算出主模截至波长
。
升级会员 