一、操作步骤
1、网格剖分
首先,打开PDE工具箱(命令pdetool),画个矩形,在按钮“三角形”上单击或从菜单Mesh中选择生成初始网格,结果如下(所谓的无组织网格)。要生成直角三角形构成的组织化网格,应进入菜单Mesh选择Parameters,然后键入最大边大小inf,网格细化几次后结果如图。
第二步是把网格导入到主工作区,可以从菜单Mesh中选择Export Mesh项,然后将数组p变为更一般的三维形式:
p(3,:)=0
然后,用下面命令将数组p和t保存到二进制文件*.mat中
save 文件名 p t
这里画的矩形是1m*1m的矩形pec,如下:
------p:3*289 记录所有点的三维坐标,一共289个点
------t:4*512 记录所有三角形的三个顶点编号
2、得到边元
运行rwg1.m,输入参数为p、t
程序主要算法是循环每个三角形,如果两个三角形有一个公共边,则这两个三角形组成一个边元。
输出参数为mesh1
------Area 1*512 每个三角形的面积
------Center 3*512 所有三角形中心的坐标
------TriangleTotal
512 三角形的个数
------EdgesTotal 736 边元的个数
------TriangleMinus 1*736
每个边元的负三角形的编号
------TrianglePlus 1*736 每个边元的正三角形的编号
------Edge_ 2*736 所有边元的公共边两端点的编号
------Edgelength 1*736 所有边元公共边的长度
3、进一步求出边元的一些参数
运行rwg.m,输入参数为mesh1
输出参数为mesh2
------RHO_Plus 3*736 每个边元正三角形自由顶点->中心点的矢量
------RHO_Minus 3*736 每个边元负三角形中心点->自由顶点的矢量
------Center_ 3*9*512
每个三角形被分成9个小三角形后,每个小三角形的中心点坐标 ------RHO__Plus 3*9*736 每个边元正三角形被分成9个三角形后,边元正三角形顶点->小三角形中心点的矢量
------RHO__Minus 3*9*736 每个边元负三角形被分成9个三角形后,小三角形中心点->边元负三角形顶点的矢量
3、计算阻抗举证Z
阻抗矩阵与问题的类型(辐射或散射)无关,但是与频率有关
运行rwg3.m,输入参数为mesh2
输出为impedance
------K = jk
------FactorA 736*1 (j * omega * mu_ * EdgeLength)/(144*pi)的转置 A前面的系数
------FactorFi 736*1 EdgeLength/(36*pi*j*omega
*epsilon_)的转置 Fi前面的系数
------Plus 表示第一个三角形作为正边元的边元号(第一次循环)个数可能是0-3
------Minus 表示第一个三角形作为负边元的边元号(第一次循环)个数可能是0-3
------D 表示第一个三角形中心到其它三角形被分成9个三角形后中心的矢量(第一次循环)
------R D的模值表示第一个三角形中心到其它三角形被分成9个三角形后中心的距离(第一次循环)
------g 格林函数
exp(-KR)/R 表示第一个三角形对其它三角形被分成9个三角形后的作用(第一次循环)
------gP 表示第一个三角形对每个边元正三角形分成的9个小三角形的格林函数(第一次循环)
------gM 表示第一个三角形对每个边元负三角形分成的9个小三角形的格林函数(第一次循环)
这里我把书中的源代码略加修改,因为感觉书中是错的,
4、计算V向量,解MoM矩阵
运行rwg4.m,输入mesh2和impedance
输出为current
5、电流可视化
运行rwg5.m,输入mesh2和current
输出为可视化图形
二、附
MATLAB程序
rwg1.m
%RWG1 Geometry calculations - all
Chapters
% Uses the structure mesh file,
e.g. platefine.mat,
% as an input.
%
% Creates the RWG edge element
for every inner edge of
% the structure. The total
number of elements is EdgesTotal.
% Outputs the following
arrays:
%
% Edge first node
number Edge_(1,1:EdgesTotal)
% Edge second node
number Edge_(2,1:EdgesTotal)
% Plus triangle
number TrianglePlus(1:EdgesTotal)
% Minus triangle
number TriangleMinus(1:EdgesTotal)
% Edge
length EdgeLength(1:EdgesTotal)
% Edge element
indicator EdgeIndicator(1:EdgesTotal)
%
% Also outputs areas and
midpoints of separate triangles:
% Triangle