基于FEM的2D矩形波导模式分析

摘要

  有限元法是里兹法与有限差分法相结合的结果，在理论上以变分为基础，在具体方法构造上又利用了有限差分法网格离散化处理的思想。它通过变分方法，使得误差函数达到最小值并产生稳定解。类比于连接多段微小直线逼近圆的思想，有限元法包含了一切可能的方法，这些方法将许多被称为有限元的小区域上的简单方程联系起来，并用其去估计更大区域上的复杂方程。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的（较简单的）近似解，然后推导求解这个域总的满足条件（如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。其核心是剖分插值，即将连续场分割为有限个单元，然后用较简单的插值函数来表示每个单元的解。
  本文以标准矩形波导BJ-100为例，使用有限元法对TE模与TM模情况下横截面的场分布情况进行了理论分析，利用MATLAB仿真并将结果与HFSS仿真结果进行了对比。

关键词

有限元法；标准矩形波导；横截面；MATLAB仿真

1.引言

  有限元法$(FEM，Finite Element Method)$是一种为求解偏微分方程边值问题近似解的数值技术。它通过变分方法，使得误差函数达到最小值并产生稳定解。。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。也正是因此特性，我们选择使用有限元法来进行有关矩形波导的2D分析。

2.矩形波导2D分析

2.1理论推导

2.1.1边值问题

$$ \begin{cases} \nabla^{2} \phi+k_{c}^{2}=0 \\ \left.\frac{\partial \phi}{\partial n}\right|_{L}=0 \end{cases} \qquad \text {L为波导壁} $$

2.1.2泛函——变分问题（无条件）

$$F \left\lbrack \phi \right\rbrack =\frac{1}{2}{\iint_{D}^{}\left\lbrack \left( \frac{\partial y}{\partial x} \right)^{2} + \left( \frac{\partial y}{\partial x} \right)^{2} - k_{c}^{2}\phi \right\rbrack}dxdy = min$$

2.1.3离散

$$F\left\lbrack \phi \right\rbrack=\sum\limits_{e = 1}^{E}{F_{e}\left\lbrack \phi \right\rbrack}$$
$$F_{e}\left\lbrack \phi \right\rbrack = F_{e}\left\lbrack \phi^{e} \right\lbrack = {\iint_{D}^{}{\frac{1}{2}\left\lbrack \left( \frac{\partial y}{\partial x} \right)^{2} + \left( \frac{\partial y}{\partial x} \right)^{2} - k_{c}^{2}\phi^{e^{2}} \right\rbrack}}dxdy$$
$$N_{s}^{e}\left( {x,y} \right) = \frac{1}{2\delta}\left( {a_{s} + b_{s}x + c_{s}y} \right)~~~s = i,j,m$$

$$ \begin{cases} {\frac{\partial\phi^{e}}{\partial x} = \left( b_{i}\phi_{i} + b_{j}\phi_{j} + b_{m}\phi_{m} \right)/2\delta} \\ {\frac{\partial\phi^{e}}{\partial y} = \left( c_{i}\phi_{i} + b_{j}\phi_{j} + c_{m}\phi_{m} \right)/2\delta} \\ \end{cases} $$

则
$$F_{e}\left\lbrack \phi^{e} \right\lbrack= \frac{1}{2}\phi_{e}^{T}K_{e}\phi_{e} - \frac{k_{c}^{2}}{2}\phi_{e}^{T}T\phi_{e}~~\phi = \left\lbrack N_{e}{\rbrack\left\lbrack \phi \right.}_{e} \right\rbrack$$
$$k_{rs}^{e} = k_{sr}^{e} = \frac{\varepsilon}{4 \bigtriangleup}\left( {b_{r}b_{s} + c_{r}c_{s}} \right)~~\left( r,s = i,j,m \right)$$
$$T_{e} = \frac{\bigtriangleup}{6}~~~~~T_{ij}^{e} = \frac{\bigtriangleup}{12}~~~~T_{im}^{e} = \frac{\bigtriangleup}{12}~~\overset{}{\Rightarrow}~~~T_{rl}^{e} = T_{lr}^{e} = \frac{\bigtriangleup}{12}~\left( 1 + \delta_{rl} \right)$$

2.1.4单元合成

$$F\left( \phi \right) \approx {\sum\limits_{e = 1}^{E}{F_{e}\left( \phi \right) = \frac{1}{2}}}\phi^{T}k\phi - \frac{k_{c}^{2}}{2}\phi^{T}T\phi$$
取极值
$$\frac{\partial F\left( \phi \right)}{\partial\phi} = {\sum\limits_{e = 1}^{E}{\frac{\partial F_{e}\left( \phi \right)}{\partial\phi} = 0}}$$
即$K\phi = k_{c}^{2}T\phi$
在$K\phi = k_{c}^{2}T\phi$中，$K$为对称阵，$T$为对称正定阵
利用平方根法（$cholesky$）令$T=LL^T$，$L$为下三角矩阵

$$ L = \begin{bmatrix} l_{11} & & \\ \vdots & \ddots & \\ l_{n1} & \cdots & l_{nn} \\ \end{bmatrix} $$

$$ L^T = \begin{bmatrix} l_{11} & \cdots & l_{n1} \\ & \ddots & \vdots \\ & & l_{nn} \\ \end{bmatrix} $$

$$ T = \begin{bmatrix} {\begin{matrix} T_{11} & T_{12} & T_{13} \\ T_{21} & T_{22} & T_{23} \\ \vdots & \vdots & \vdots \\ \end{matrix}~~~~\begin{matrix} \ldots & T_{1n} \\ \ldots & T_{2n} \\ \ddots & \vdots \\ \end{matrix}} \\ {\begin{matrix} T_{n1} & T_{n2} & T_{n3} \\ \end{matrix}~~~~\begin{matrix} \ldots & T_{nn} \\ \end{matrix}} \\ \end{bmatrix} $$

则有
$$T_{11} = l_{11} \cdot l_{11}~~~~{T}_{21} = l_{21} \cdot l_{11}$$
$$T_{i1} = l_{i1} \cdot l_{11}~~i = 1,2 \cdot \cdot \cdot \cdot m~~~~~l_{kk} = \sqrt{a_{kk} - {\sum\limits_{p = 1}^{k - 1}l_{kp}^{2}}}$$
$$l_{kk} = a_{ik} - {\sum\limits_{p = 1}^{k - 1}{l_{ip}l_{kp}/l_{kk}}}$$
$$K\phi = k_{c}^{2}T\phi\overset{T = {LL}^{T}}{\rightarrow}L^{- 1}KL^{- T}\left( {L^{T}\phi} \right) = k_{c}^{2}\left( {L^{T}\phi} \right)$$
$$V = L^{- 1}L^{- T}\left( 已知 \right)\overset{z = \phi}{\rightarrow}VZ = k_{c}^{2}z$$

2.1.5解$VZ=k_{c}^{2}z$

  用$householder$矩阵将$V$矩阵上$Hessenberg$化。$H' ▽H=F$，$H$为$householder$矩阵，又因$V$对称，则$F$为三对角矩阵，n个$householder$矩阵相乘，最后得三对角矩阵
$$N_{s}^{e} \bigtriangleup \left( {x,y} \right) = \frac{1}{2 \bigtriangleup}\left( {a_{s} + b_{s}y + c_{s}y} \right)$$
$$\bigtriangledown N_{k}^{e} = \frac{1}{2 \bigtriangleup}\left( {\frac{\partial}{\partial x}{\vec e}_{x} + \frac{\partial}{\partial y}{\vec e}_{y}} \right)N_{k}^{e}\left( {x,y} \right) = \frac{1}{2 \bigtriangleup}\left( {b_{k}{\vec {e}}_{x} + c_{k}{\vec {e}}_{y}} \right)$$
$$N_{l}^{e} \bigtriangledown N_{k}^{e} = \frac{1}{2 \bigtriangleup}\left( {a_{l} + b_{l}y + c_{l}y} \right) \cdot \frac{1}{2 \bigtriangleup}\left( {b_{k}{\vec{e}}_{x} + c_{k}{\vec{e}}_{y}} \right)$$

$$ N_{l}^{e} \bigtriangledown N_{k}^{e} = \frac{1}{2 \bigtriangleup}\left( {a_{l} + b_{l}y + c_{l}y} \right) \cdot \frac{1}{2 \bigtriangleup}\left( {b_{k}{\vec{e}}_{x} + c_{k}{\vec{e}}_{y}} \right) = \frac{1}{4 \bigtriangleup^{2}}\left( {a_{l} + b_{l}y + c_{l}y} \right)\left( {b_{k}{\vec{e}}_{x} + c_{k}{\vec{e}}_{y}} \right) $$

同理
$$N_{k}^{e} \bigtriangledown N_{l}^{e} = \frac{1}{4 \bigtriangleup^{2}}\left( {a_{k} + b_{k}y + c_{k}y} \right)\left( {b_{l}{\vec{e}}_{x} + c_{l}{\vec{e}}_{y}} \right)$$
$$N_{lk}^{e}\left( {x,y} \right) = \left\lbrack \left( {a_{l}b_{k} - {a_{k}b}_{l} + c_{l}b_{k}y - c_{k}b_{l}y} \right){\vec {e}}_{x} - \left( {a_{l}c_{k} - {a_{k}c}_{l} + b_{l}c_{k}x - b_{k}c_{l}x} \right){\vec {e}}_{y} \right\rbrack l_{lk}^{e}$$
$$n\left( {l,k} \right) = i~~~~~\left( l,k点确定棱边i~ \right)$$
记
$$p = a_{l}b_{k} - {a_{k}b}_{l}~~~~~~~~q = a_{l}c_{k} - {a_{k}c}_{l}$$

$$ \bigtriangledown \times N_{i}^{e} = \left| \begin{matrix} {\vec{e}}_{x} & {\vec{e}}_{y} & {\vec{e}}_{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ {\left. p + \left( b \right._{k}c_{l} - b_{l}c_{k} \right)y} & {q + \left( b_{l}c_{k} - b_{k}c_{l} \right)} & 0 \\ \end{matrix} \right| = \frac{2l_{lk}^{e}}{4 \bigtriangleup^{2}}\left( b_{l}c_{k} - b_{k}c_{l} \right){\vec{e}}_{z} $$

2.2MATLAB仿真（TE模）

2.2.1建立模型与网格划分

利用MATLAB的PDE model生成有限元模型（图2-1-1）。

图2-2-1 有限元模型与网格划分（MaxElementSize=1.006）

图2-2-2利用generateMesh函数产生网格划分

图2-2-3 利用meshToPet函数获取全局编号

2.2.2解特征值方程

由标量亥姆霍兹方程经变分计算导出的特征值方程为：$[K][\phi]=k_c^2 [T][\phi]$，针对TE模为Hz。
利用自定义的Get_K_T函数，通过对单元系数矩阵的计算及全局拓展获取特征值方程的$K$矩阵（图2-2-4~图2-2-5）和$T$矩阵（图2-2-6~图2-2-7），求解出特征值矩阵特征值向量矩阵。
Get_K_T.m

function [K,T] = Get_K_T(Elements,Nodes)
%GET_K_T 计算系数矩阵K,T

%相对介电常数
Epsilon = 1;

%获取元素数量和节点数
E_num = size(Elements,1);
N_num = size(Nodes,1);

%获取节点坐标
x = Nodes(:,1);
y = Nodes(:,2);

%创建节点及其abc值矩阵
a = zeros(N_num,3);
b = zeros(N_num,3);
c = zeros(N_num,3);

%创建三角元面积矩阵
Area = zeros(E_num,1);

%创建临时矩阵
temp = zeros(2,3);

%创建三角元系数矩阵
K_all = zeros(N_num);
T_all = zeros(N_num);

for num = 1:E_num
    %获取全局节点编号
    [i,j,m] = deal(Elements(num,1),Elements(num,2),Elements(num,3));
    %获取编号对应的坐标值
    [xi,xj,xm] = deal(x(i),x(j),x(m));
    [yi,yj,ym] = deal(y(i),y(j),y(m));
    %获取a(i,j,m)、b(i,j,m)、c(i,j,m)值
    [a(num,1),a(num,2),a(num,3)] = deal(xj*ym - xm*yj,xm*yi - xi*ym,xi*yj - xj*yi);
    [b(num,1),b(num,2),b(num,3)] = deal(yj-ym,ym-yi,yi-yj);
    [c(num,1),c(num,2),c(num,3)] = deal(xm-xj,xi-xm,xj-xi);
    %获取每个三角元的面积
    Area(num) = 0.5*(b(num,1)*c(num,2)- b(num,2)*c(num,1));
end

%清除变量，后面仍会使用
clear num i j m;

for num = 1:E_num
    %获取bi、bj、bm;ci、cj、cm
    temp(1,:) = b(num,:);
    temp(2,:) = c(num,:);
    %求每个三角元系数矩阵
    Ke = ((temp')*temp).*(0.25*Epsilon/Area(num));
    %获取全局节点编号
    [i,j,m] = deal(Elements(num,1),Elements(num,2),Elements(num,3));
    %拓展为全局系数矩阵
    K_all(i,i) = K_all(i,i) + Ke(1,1);
    K_all(i,j) = K_all(i,j) + Ke(1,2);
    K_all(j,i) = K_all(j,i) + Ke(2,1);
    K_all(i,m) = K_all(i,m) + Ke(1,3);
    K_all(m,i) = K_all(m,i) + Ke(3,1);
    K_all(j,j) = K_all(j,j) + Ke(2,2);
    K_all(j,m) = K_all(j,m) + Ke(2,3);
    K_all(m,j) = K_all(m,j) + Ke(3,2);
    K_all(m,m) = K_all(m,m) + Ke(3,3);
    %求解T矩阵，并拓展为全局编号
    Te = [1/6,1/12,1/12;1/12,1/6,1/12;1/12,1/12,1/6].*Area(num);
    T_all(i,i) = T_all(i,i) + Te(1,1);
    T_all(i,j) = T_all(i,j) + Te(1,2);
    T_all(j,i) = T_all(j,i) + Te(2,1);
    T_all(i,m) = T_all(i,m) + Te(1,3);
    T_all(m,i) = T_all(m,i) + Te(3,1);
    T_all(j,j) = T_all(j,j) + Te(2,2);
    T_all(j,m) = T_all(j,m) + Te(2,3);
    T_all(m,j) = T_all(m,j) + Te(3,2);
    T_all(m,m) = T_all(m,m) + Te(3,3);
end

K = K_all;
T = T_all;

end

图2-2-4 特征值方程的K矩阵（TE；MaxElementSize=1.006）

图2-2-5 K矩阵的散点图（TE；MaxElementSize=1.006）

图2-2-6 特征值方程的T矩阵（TE；MaxElementSize=1.006）

图2-2-7 T矩阵的散点图（TE；MaxElementSize=1.006）

2.2.3 MaxElementSize=1.006时TE模的仿真结果

  MaxElementSize设定为1.0006（默认值）时，$k_c$的实际值与$k_c$的理论值对比（图2-2-8）与前六个TE模仿真结果如下（图2-2-9）。

图2-2-8 kc的实际值（左）与kc的理论值（右）对比（TE；MaxElementSize=1.006）

图2-2-9 MaxElementSize设定为1.0006时前六个TE模的仿真结果

2.2.4 MaxElementSize=2时TE模的仿真结果

  MaxElementSize设定为2时，$k_c$的实际值与$k_c$的理论值对比（图2-2-10）与前六个TE模仿真结果如下（图2-2-11）。

图2-2-10 kc的实际值（左）与kc的理论值（右）对比（TE；MaxElementSize=2）

图2-2-11 MaxElementSize设定为2时前六个TE模的仿真结果

2.2.5 MaxElementSize=0.5时TE模的仿真结果

  MaxElementSize设定为0.5时，$k_c$的实际值与$k_c$的理论值对比（图2-2-12）与前六个TE模仿真结果如下（图2-2-13）。

图2-2-12 kc的实际值（左）与kc的理论值（右）对比（TE；MaxElementSize=0.5）

图 2-2-13 MaxElementSize设定为0.5时前六个TE模的仿真结果

2.2.6 主程序代码（TE）

%基于有限元方法的BJ-100矩形波导正面TE模磁场的2D分析程序
%If you can't see the Chinese character,please change your character encoding format

%清除所有数据和命令行窗口
clear all; clc;

%创建2D矩形模型
BJ100_2D = createpde(1);
R1 = [3,4,-11.43,11.43,11.43,-11.43,-5.08,-5.08,5.08,5.08]';    %矩形大小
g = decsg(R1);  %分解成最小几何
geometryFromEdges(BJ100_2D,g);  %由边界生成矩形模型
%画出确定的矩形
figure(1);
pdegplot(BJ100_2D,'EdgeLabels','on');
xlim([-12,12]);
ylim([-6,6]);
axis equal;

%确定neumann边界条件 TE模只有neumann条件
applyBoundaryCondition(BJ100_2D,'neumann','Edge',1,'g',0,'q',0);
applyBoundaryCondition(BJ100_2D,'neumann','Edge',2,'g',0,'q',0);
applyBoundaryCondition(BJ100_2D,'neumann','Edge',3,'g',0,'q',0);
applyBoundaryCondition(BJ100_2D,'neumann','Edge',4,'g',0,'q',0);

%线性网格剖分
mesh = generateMesh(BJ100_2D,'GeometricOrder','linear');
% %画出网格剖分
% figure(2);
% pdeplot(BJ100_2D);

%获得P E T值
[p,e,t] = meshToPet(mesh);
Elements = t(1:3,:)';   %获取所有三角形元素
Edges = e(1:2,:)';   %获取所有边
Nodes = p';      %获取所有节点

%获取节点编号
N_num = size(Nodes,1);

%获取系数矩阵的值
[K,T] = Get_K_T(Elements,Nodes);

%画出系数矩阵K的散点图
figure(3);
subplot(1,2,1);imagesc(full(K)~=0);
title(sprintf('K(%d,%d)',N_num,N_num));colorbar;axis('image');
subplot(1,2,2);imagesc(full(T)~=0);
title(sprintf('T(%d,%d)',N_num,N_num));colorbar;axis('image');

%特征值求解设置
opts.isreal = 1;%为实矩阵
%求解特征值和特征向量
[V,D,flag] = eigs(K, T, 11, 'smallestabs', opts); %求解10个特征值 第一个特征值为0 忽略
fprintf('收敛标志:%d\n', flag);
    
Kc2 = diag(D); %获取对角线上的特征值
if ~issorted(Kc2)
       [Kc2,idx] = sort(Kc2);
       V = V(:, idx);  %获取对应特征向量
end
Kc = sqrt(Kc2);

%去掉特征值约为零的模
fprintf('第一特征值 = %f', Kc(1));
if abs(Kc(1)) < 0.0001
       Kc2 = Kc2(2:end);
       Kc = Kc(2:end);
       V = V(:, 2:end);
end

%画出前六个TE模
figure(4);
for eigen = 1:6
   subplot(2,3,eigen);
       Hz_all = V(:,eigen);
       trisurf(Elements,Nodes(:,1),Nodes(:,2),Hz_all);
       view(2); %平面视图
       xlabel('x'); ylabel('y'); zlabel('H_z');axis('equal');
       title(sprintf('Kc=%0.4f',Kc(eigen)));
       colorbar('Location','westoutside');
end 
suptitle('前六个TE模的Hz');

2.3 MATLAB仿真（TM模）

2.3.1 MaxElementSize=1.006时TM模的仿真结果

  与TE模的仿真过程类似，我们采用了建立模型与网格划分（图2-3-1）、解特征值方程（图2-3-2~图2-3-5）与查看仿真结果（图2-3-6~图2-3-7）这三个步骤对TM模进行了仿真（MaxElementSize=1.0006）。

图2-3-1 有限元模型与网格划分（TM；MaxElementSize=1.006）

图2-3-2 特征值方程的K矩阵（TM；MaxElementSize=1.006）

图2-3-3 K矩阵的散点图（TM；MaxElementSize=1.006）

图2-3-4 特征值方程的T矩阵（TM；MaxElementSize=1.006）

图2-3-5 T矩阵的散点图（TM；MaxElementSize=1.006）

图2-3-6 kc的实际值（左）与kc的理论值（右）对比（TM；MaxElementSize=1.006）

图2-3-7 MaxElementSize设定为1.0006时前六个TM模的仿真结果

2.3.2 MaxElementSize=2时TM模的仿真结果

  MaxElementSize设定为2时，$k_c$的实际值与$k_c$的理论值对比（图2-3-8）与前六个TM模仿真结果如下（图2-3-9）。

图2-3-8 kc的实际值（左）与kc的理论值（右）对比（TM；MaxElementSize=2）

图2-3-9 MaxElementSize设定为2时前六个TM模的仿真结果

2.3.3 MaxElementSize=0.5时TM模的仿真结果

  MaxElementSize设定为0.5时，$k_c$的实际值与$k_c$的理论值对比（图2-3-10）与前六个TM模仿真结果如下（图2-3-11）。

图2-3-10 kc的实际值（左）与kc的理论值（右）对比（TM；MaxElementSize=0.5）

图2-3-11 MaxElementSize设定为0.5时前六个TM模的仿真结果

2.3.4 主程序代码（TM）

%基于有限元方法的BJ-100矩形波导正面TM模电场的2D分析程序
%If you can't see the Chinese character,please change your character encoding format

%清除所有数据和命令行窗口
clear all; clc;

%创建2D矩形模型
BJ100_2D = createpde(1);
R1 = [3,4,-11.43,11.43,11.43,-11.43,-5.08,-5.08,5.08,5.08]';    %矩形大小
g = decsg(R1);  %分解成最小几何
geometryFromEdges(BJ100_2D,g);  %由边界生成矩形模型
% %画出确定的矩形
% figure(1);
% pdegplot(BJ100_2D,'EdgeLabels','on');
% xlim([-12,12]);
% ylim([-6,6]);
% axis equal;

%确定dirichlet边界条件 TM模只有dirichlet条件
applyBoundaryCondition(BJ100_2D,'dirichlet','Edge',1,'h',1,'r',0);
applyBoundaryCondition(BJ100_2D,'dirichlet','Edge',2,'h',1,'r',0);
applyBoundaryCondition(BJ100_2D,'dirichlet','Edge',3,'h',1,'r',0);
applyBoundaryCondition(BJ100_2D,'dirichlet','Edge',4,'h',1,'r',0);

%线性网格剖分
mesh = generateMesh(BJ100_2D,'GeometricOrder','linear');
% %画出网格剖分
% figure(2);
% pdeplot(BJ100_2D);

%获得P E T值
[p,e,t] = meshToPet(mesh);
Elements = t(1:3,:)';   %获取所有三角形元素
Edges = e(1:2,:)';   %获取所有边
Nodes = p';      %获取所有节点

%获取节点编号
E_num = size(t,1);
%边界编号
B_num = size(Edges,1);
%节点编号
N_num = size(Nodes,1);

%获取系数矩阵的值
[K,T] = Get_K_T(Elements,Nodes);

%画出系数矩阵K的散点图
figure(3);
subplot(1,2,1);imagesc(K~=0);
title(sprintf('K(%d,%d)',N_num,N_num));colorbar;axis('image');
subplot(1,2,2);imagesc(T~=0);
title(sprintf('T(%d,%d)',N_num,N_num));colorbar;axis('image');

%所有编号
all_n = (1:N_num).';
%unique去掉重复行
edge_n = unique(Edges(:));
%setdiff(A,B)函数返回在向量A中却不在向量B中的元素，并且C中不包含重复元素，并且从小到大排序。
int_n = setdiff(all_n, edge_n);
%内部元素编号
in_num = size(int_n,1);

%表示内部或边缘节点的逻辑数组
int_or_edge = ones(N_num, 1);
int_or_edge(edge_n) = 0;
int_or_edge = logical(int_or_edge);

%特征值求解设置
opts.issym = 1;%对称阵
opts.isreal = 1;%实矩阵
%求解特征值和特征向量
[V,D,flag] = eigs(K(int_n, int_n), T(int_n, int_n), 10, 'smallestabs', opts); %求解10个特征值 第一个特征值为0 忽略
fprintf('收敛标志:%d\n', flag);
Kc2 = diag(D);
if ~issorted(Kc2)
       [Kc2,idx] = sort(Kc2);
       V = V(:, idx);
end
Kt = sqrt(Kc2);

%画出前六个TM模
figure(4);
for eigen=1:6
       subplot(2,3,eigen);
       Ez_all = zeros(N_num, 1);
       Ez_all(int_n) = V(:,eigen);
       trisurf(Elements,Nodes(:,1),Nodes(:,2),Ez_all); 
       view(2);
       xlabel('x'); ylabel('y'); zlabel('E_z');axis('equal');
       title(sprintf('Kc=%0.4f',Kt(eigen)));
       colorbar('Location','westoutside');
end
suptitle('前六个TM模的Ez');

3. HFSS仿真结果

3.1 仿真结果对比（TE）

  继MATLAB仿真后，我们又将模型放入HFSS软件中进行仿真测试，所得TE模式仿真结果对比图如下（图3-1-1~图3-1-6）。

图 3-1-1 TE10 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

图 3-1-2 TE20 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

图 3-1-3 TE01 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

图 3-1-4 TE11 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

图 3-1-5 TE30 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

图 3-1-6 TE21 模 MATLAB 仿真（左）与 HFSS 仿真（右）结果对比

3.2 仿真结果对比（TM）

  与TE模类似，继MATLAB仿真后，我们又将模型放入HFSS软件中进行仿真测试，所得TM模式仿真结果对比图如下（图3-2-1~图3-2-6）。

图3-2-1 TM11模MATLAB仿真（左）与HFSS仿真（右）结果对比

图3-2-2 TM21模MATLAB仿真（左）与HFSS仿真（右）结果对比

图3-2-3 TM31模MATLAB仿真（左）与HFSS仿真（右）结果对比

图3-2-4 TM41模MATLAB仿真（左）与HFSS仿真（右）结果对比

图3-2-5 TM12模MATLAB仿真（左）与HFSS仿真（右）结果对比

图3-2-6 TM22模MATLAB仿真（左）与HFSS仿真（右）结果对比

4.结语

  我们利用有限元法，以标准矩形波导BJ-100为例，对TE模与TM模情况下横截面的场分布情况进行了理论分析，利用MATLAB仿真并将结果与HFSS仿真结果进行了对比，实现了对算法的学习利用，并求解了实际问题，准确度也较高，成功完成了设计目标。

参考文献

• 赵 霞. 基于有限元法的波导主模特性分析, 科技信息（学术研究）, 2007,5.
• 吴雅祥, 刘仲武, 王开祥. 有限元求解电磁场问题的误差分析, 磁性材料及器件, 2019,1.
• 邓素芬，杨显清，陈友臸. 基于有限元法的脊波导特征值分析, 电子对抗技术, 2005,(20).
• 文 庆. 微波无源器件有限元—模式匹配法研究, 电子科技大学, 2016,11.
• G.E. Fasshauer，J.G. Zhang. Numer. Algorithms, 2007.
• P. Silvester. High-Order Finite Element Waveguide Analysis (Program Descriptions), 1969.