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

warning: 这篇文章距离上次修改已过559天,其中的内容可能已经有所变动。
info:博客公式可右键缩放
冉春霖,周傲杰,马悦敏,刘庚,贠泽宇,周国立 1
(1. 电子科技大学,成都)

摘要

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

利用$generateMesh$函数产生网格划分,默认尺寸为$1mm$(图2-2-2)。
图2-2-2利用generateMesh函数产生网格划分图2-2-2利用generateMesh函数产生网格划分

利用$meshToPet$函数获取元素、点和边的全局编号(图2-2-3)。
图2-2-3  利用meshToPet函数获取全局编号图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-4 特征值方程的K矩阵(TE;MaxElementSize=1.006)

图2-2-5  K矩阵的散点图(TE;MaxElementSize=1.006)图2-2-5 K矩阵的散点图(TE;MaxElementSize=1.006)

图2-2-6  特征值方程的T矩阵(TE;MaxElementSize=1.006)图2-2-6 特征值方程的T矩阵(TE;MaxElementSize=1.006)

图2-2-7  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-8 kc的实际值(左)与kc的理论值(右)对比(TE;MaxElementSize=1.006)

图2-2-9  MaxElementSize设定为1.0006时前六个TE模的仿真结果图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-10 kc的实际值(左)与kc的理论值(右)对比(TE;MaxElementSize=2)

图2-2-11  MaxElementSize设定为2时前六个TE模的仿真结果图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-12 kc的实际值(左)与kc的理论值(右)对比(TE;MaxElementSize=0.5)

图 2-2-13 MaxElementSize设定为0.5时前六个TE模的仿真结果图 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-1 有限元模型与网格划分(TM;MaxElementSize=1.006)

图2-3-2  特征值方程的K矩阵(TM;MaxElementSize=1.006)图2-3-2 特征值方程的K矩阵(TM;MaxElementSize=1.006)

图2-3-3  K矩阵的散点图(TM;MaxElementSize=1.006)图2-3-3 K矩阵的散点图(TM;MaxElementSize=1.006)

图2-3-4  特征值方程的T矩阵(TM;MaxElementSize=1.006)图2-3-4 特征值方程的T矩阵(TM;MaxElementSize=1.006)

图2-3-5  T矩阵的散点图(TM;MaxElementSize=1.006)图2-3-5 T矩阵的散点图(TM;MaxElementSize=1.006)

图2-3-6  kc的实际值(左)与kc的理论值(右)对比(TM;MaxElementSize=1.006)图2-3-6 kc的实际值(左)与kc的理论值(右)对比(TM;MaxElementSize=1.006)

图2-3-7  MaxElementSize设定为1.0006时前六个TM模的仿真结果图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-8 kc的实际值(左)与kc的理论值(右)对比(TM;MaxElementSize=2)

图2-3-9  MaxElementSize设定为2时前六个TM模的仿真结果图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-10 kc的实际值(左)与kc的理论值(右)对比(TM;MaxElementSize=0.5)

图2-3-11  MaxElementSize设定为0.5时前六个TM模的仿真结果图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-1 TE10 模 MATLAB 仿真(左)与 HFSS 仿真(右)结果对比

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

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

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

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

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

由图3-1-1~图3-1-6不难看出,TE模式下MATLAB仿真与HFSS仿真结果契合度很高。

3.2 仿真结果对比(TM)

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

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

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

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

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

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

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

由图3-2-1~图3-2-6不难看出,TM模式下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.

添加新评论