基于FDM的2D矩形波导模式分析
一、题目分析
用差分法计算$a=b=\frac{11\pi}{7}$的矩形波导中$TM_{11}、TE_{10}$模的场分布、截止频率、特征阻抗（网格划分$ℎ=\frac{\pi}{7}$）
对波导横截面二维标量波动方程定解问题，在波导内：
$$\nabla_{t}^{2}\varphi + K_{c}^{2}\varphi^{2} = \frac{\partial^{2}\varphi}{\partial x^{2}} + \frac{\partial^{2}\varphi}{\partial y^{2}} + K_{c}^{2}\varphi^{2} = 0$$
采用正方形网格划分的有限差分法：
$$\varphi_{i,j} = \frac{1}{4}\left( {\varphi_{i + 1,j} + \varphi_{i,j + 1} + \varphi_{i - 1,j} + \varphi_{i,j - 1} - h^{2}K_{c}^{2}} \right)$$
对每个节点的值，我们可以采用迭代法和直接法求解。
二、迭代法
这里我们采用迭代较快的超松驰迭代法，迭代格式为：
$$\varphi_{i,j}^{({n + 1})} = \varphi_{i,j}^{(n)} + \omega\left\{ {\frac{1}{4 - h^{2}K_{c}^{2}}\left\{ {\varphi_{i + 1,j}^{\{ n\}} + \varphi_{i,j + 1}^{\{ n\}} + \varphi_{i - 1,j}^{\{{n + 1}\}} + \varphi_{i,j - 1}^{\{{n + 1}\}}} \right\} - \varphi_{i,j}^{\{ n\}}} \right\}$$
ω为松弛因子(1<ω<2)，加速求解的迭代；场的初值需介于其最大值和最小边界值，可以减小迭代次数。
2.1 ω的选取
对于正方形场域的第一类边值问题：
$$\omega_{best} = \frac{2}{1 + {\sin\left( \frac{\pi}{m - 1} \right)}}$$
$m$为场域划分的一边节点数
对于正方形场域的第二类边值问题，ω可采用下式进行估算：
$$\omega_{best} = \frac{2}{1 + \sqrt{1 - \lambda}}$$
$$\lambda = {\lim\limits_{n\rightarrow\infty}\frac{\varepsilon^{n + 1}}{\varepsilon^{n}}}$$
$$\varepsilon^{n + 1} = {\max\limits_{}\left| {\varphi^{n + 1} - \varphi^{n}} \right|}$$
$$\varepsilon^{n} = {\max\limits_{}\left| {\varphi^{n} - \varphi^{n - 1}} \right|}$$
2.2 $K_c$的修正与$f_c$
在场值的迭代之前，我们需要对$K_c$设定一个初值，初值设定的好坏直接影响最终$K_c$迭代的精度。在理论上，只要$K_c$足够小或与主模截止波数接近即可。
为使求解的$K_c$精度较高，在对φ迭代3~5次后，对$K_c$进行一次修正。由：
$${\iint_{s}^{}{\varphi\nabla_{t}^{2}}}\varphi ds + {\iint_{s}^{}K_{c}^{2}}\varphi^{2}ds = 0$$
$$K_{c}^{2} = - \frac{{\iint_{s}^{}{\varphi\nabla_{t}^{2}}}\varphi ds}{{\iint_{s}^{}\varphi^{2}}ds}$$
又
$$\nabla_{t}^{2}\varphi \approx \frac{\sum{\varphi - 4\varphi_{0}}}{h^{2}}$$
得
$$K_{c}^{2}h^{2} = - \frac{\sum\limits_{i,j}{\varphi_{i,j}\left( {\varphi_{i + 1,j} + \varphi_{i,j + 1} + \varphi_{i - 1,j} + \varphi_{i,j - 1} - 4\varphi_{i,j}} \right)}}{\sum\limits_{i,j}\varphi_{i,j}^{2}}$$
由波导截止特性知：
$$f_{c} = \frac{cK_{c}}{2\pi} \cdots c = 2.99792458 \times 10^{8}$$
2.3 迭代法模式分析
对于TM模，在波导壁处
$$\left. E_{Z} \right|_{c} = 0$$
在任一网格内点有：
$${E_{Z}}_{i + 1,j} + {E_{Z}}_{i,j + 1} + {E_{Z}}_{i - 1,j} + {E_{Z}}_{i,j - 1} - {4E_{Z}}_{i,j} + h^{2}K_{c}^{2}{E_{Z}}_{i,j} = 0$$
需要注意的是，在对场赋初值的时候，由于波导壁上的$E_Z=0$，若对内点赋零值，零解自动满足亥姆霍兹方程，但无意义，故至少需选择一个内点场值不为零。
由特征阻抗的定义方式，对于TM模求V较为方便：
$$Z_{c} = \frac{V^{2}}{P} = \frac{Z_{TM11}\beta}{\sqrt{K_{c}^{2} + \beta^{2}}}\frac{\left( {E_{Zmax} - E_{Zmin}} \right)^{2}}{K_{c}^{2}{\iint_{s}^{}E_{Z}^{2}}ds}$$
整理得：
$$Z_{c} = Z_{TM11}\sqrt{1 - \left( \frac{f_{c}}{f} \right)^{2}}\frac{- \left( {E_{Zmax} - E_{Zmin}} \right)^{2}}{{\sum\limits_{i,j}{{E_{Z}}_{i,j}\left( {{E_{Z}}_{i + 1,j} + {E_{Z}}_{i,j + 1} + {E_{Z}}_{i - 1,j} + {E_{Z}}_{i,j - 1} - {4E_{Z}}_{i,j}} \right)}}\mathrm{\Delta}S_{i,j}}$$
内点：$∆S_(i,j)=h^2$ 拐点：$∆S_(i,j)=h^2/4$ 边界点：$∆S_(i,j)=h^2/2$
对于TE模，在波导壁处
$$\left. \frac{\partial H_{Z}}{\partial n} \right|_{c} = 0$$
为求得边界的值，需要在边界外侧设置一排虚拟的网格节点
1.jpg如图，在所求场域外再画一排虚拟节点（如虚线所示）
由边界条件知：
$${H_{Z}}_{1,j} = {H_{Z}}_{3,j} \qquad {H_{Z}}_{n,j} = {H_{Z}}_{n-2,j}$$
$${H_{Z}}_{i,j} = {H_{Z}}_{i,3} \qquad {H_{Z}}_{1,n} = {H_{Z}}_{i,n-2}$$
在任一网格内点有：
$${H_{Z}}_{i + 1,j} + {H_{Z}}_{i,j + 1} + {H_{Z}}_{i - 1,j} + {H_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$$
需要注意的是，在对场赋初值的时候，由于波导壁上的$\left. \frac{\partial H_{Z}}{\partial n} \right|_{c} = 0$，若对内点赋相同值，相同值满足亥姆霍兹方程及梯度边界条件，但无意义，故赋值不能全部相等。
由特征阻抗的定义方式，对于TE模求I较为方便：
$$Z_{c} = \frac{P}{I^{2}} = \frac{Z_{TE10}\sqrt{K_{c}^{2} + \beta^{2}}}{\beta}\frac{K_{c}^{2}{\iint_{s}^{}H_{Z}^{2}}ds}{{\left( {H_{Zmax} - H_{Zmin}} \right)^{2}K}_{c}^{2}}$$
整理得：
$$Z_{c} = Z_{TE10}\sqrt{1 - \left( \frac{f_{c}}{f} \right)^{2}}\frac{\sum\limits_{i,j}{{H_{Z}}_{i,j}\left( {{H_{Z}}_{i + 1,j} + {H_{Z}}_{i,j + 1} + {H_{Z}}_{i - 1,j} + {H_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j}} \right)\mathrm{\Delta}S_{i,j}}}{{- \left( {H_{Zmax} - H_{Zmin}} \right)}^{2}}$$
三、直接法
由上述差分格式，针对标量亥姆霍兹方程，可以安装节点分解为求解N个特征方程：
$$\left\lbrack K \right\rbrack\left\lbrack \varphi \right\rbrack + \left( {K_{c}h} \right)^{2}\left\lbrack \varphi \right\rbrack = 0$$
其中K为系数矩阵，K可由不同问题的差分表达式系数合并而来，并且对于不同的边界条件，K内相应值不同。
由此特征值方程可以求出$(K_c h)^2$及φ，从而得到截止波数、截止频率及场值。
3.1 边界条件处理
对于TM模，边界条件为${E_Z |}_c$=0，无需过多处理；对于TE模，在不含拐点边界有$\left. \frac{\partial H_{Z}}{\partial n} \right|_{c} = 0$：
左边界：$2{H_{Z}}_{i + 1,j} + {H_{Z}}_{i,j + 1} + {H_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
右边界：${H_{Z}}_{i,j + 1} + 2{H_{Z}}_{i - 1,j} + {H_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
上边界：${H_{Z}}_{i + 1,j} + {H_{Z}}_{i - 1,j} + {{2H}_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
下边界：${H_{Z}}_{i + 1,j} + {{2H}_{Z}}_{i,j + 1} + {H_{Z}}_{i - 1,j} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
在拐点处：
左下拐点：$2{H_{Z}}_{i + 1,j} + {{2H}_{Z}}_{i,j + 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
右上拐点：$2{H_{Z}}_{i - 1,j} + 2{H_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
左上拐点：${{2H}_{Z}}_{i + 1,j} + {{2H}_{Z}}_{i,j - 1} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
右下拐点：${{2H}_{Z}}_{i,j + 1} + 2{H_{Z}}_{i - 1,j} - {4H_{Z}}_{i,j} + h^{2}K_{c}^{2}{H_{Z}}_{i,j} = 0$
在填充系数矩阵时非拐点边界点需额外理，拐点不参与计算可以不必处理
3.2 系数矩阵填充
内点是标准的5点差分格式，在K矩阵中：
令$i,j$为所有未知数行列值
某未知数自身$i=j$的系数为：$K_{i,j}=-4$
某未知数的上下左右未知数系数为：$K_{i,j}=1$
其余无联系的未知数：$K_{i,j}=0$
则内点的系数矩阵为：
2.png对非拐点的边界上的某一未知数，若有一相邻未知数在同一行/列，且位于边界的内侧，则该相邻未知数的系数为2。
3.png上述$K_{side}$为一般情况的形式，对不同的边界条件，$K_{side}$可能非对称非正定。总体K矩阵为二者的和，在编程的过程中，我们可能需要对其取负以得到正的特征值。
3.3 方程求解及绘图
利用MATLAB的$[V,D]=eig(K)$函数，可以求得方程的特征值矩阵D（即$h^2 K_c^2$和特征向量矩阵V（即φ），需要注意的是，某些特征值十分接近零是无效的，需要去掉。
对所得V矩阵取列向量即为某模对应的列向量值，为了绘制场图，我们需将所得列向量变形为对应网格节点的矩阵形式。
在得到截止波数后，可以求得截止频率和特征阻抗，对于$TE_{10}$模
$$Z_{c} = \frac{a}{b}\frac{Z_{0}}{\sqrt{1 - \left( \frac{f_{c}}{f} \right)^{2}}}$$
四、程序框图
对$TM_{11}$模，由边界条件的易处理性，我们采用迭代法可以轻松求解出截止频率
对$TE_{10}$模，由迭代法为取极限求截止波数、场值，而题设为正方形波导，$TE_{10}$ 、$TE_{01}$模截止波数相同，求解出的场为两者叠加，且边界条件在更新迭代的过程中容易导致场值不收敛，因此此采用直接法求解。
在得到$E_z$或$H_z$的标量值后，利用纵横关系式即可得到对应的$\overset{\rightarrow}{H_{t}}$、$\overset{\rightarrow}{E_{t}}$值，由于得到的场值为离散值，在横纵关系式中，求导可由中心差分格式（内点）、向前/后差分格式（边界点）近似。
4.1 迭代法求解程序流程图
4.png4.2 直接法求解程序流程图
5.png五、仿真结果
5.1 $TM_{11}$模
计算结果（$Z_c$为工作频率50MHz时的值）：
6.png理论$K_c=0.899954$、理论$f_c=42.939916MHz$
$E_z$矩阵及曲面图：
7.png
8.png与HFSS对比(同尺寸比例波导截面)：
$TM_{11}$模$E_z$场分布对比
9.png
10.png$TM_{11}$模$\overset{\rightarrow}{H_{t}}$场分布对比
11.png
12.png5.2 $TE_{10}$模
计算结果（$Z_c$为工作频率40MHz时的值）：
13.png理论$K_c=0.636363$、理论$f_c=30.363106MHz$
$H_z$矩阵及曲面图：
14.png
15.png与HFSS对比(同尺寸比例波导截面)：
$TE_{10}$模$H_z$场分布对比（MATLAB图像蓝色部分为负值表示向内凹陷，HFSS图像左上箭头向里，两者对应。）
16.png
17.png$TE_{10}$模$\overset{\rightarrow}{E_{t}}$场分布对比
18.png
19.png对比可知，两个模式的场分布结果与HFSS十分吻合。
六、MATLAB源代码
6.1 迭代法求解$TM_{11}$模
%基于有限差分迭代法求解矩形波导中TM11模 截面内场分布、截止频率fc和特性阻抗Zc
% if you can't see the Chinese character, please change your coding format
% written by yuzhongzhibi
% https://www.drawrain.com/index.php/archives/160
%理论截止波数 K0 = 0.899954
%理论截止频率 f0 = 42.939916MHz
clear all; clc;
%%%%%%% 初始化 %%%%%%%
%矩形波导截面尺寸
a = 11*pi/7;
b = 11*pi/7;
%正方形网格划分尺寸
h = pi/7;
%网格划分行、列数
m = round(a/h) + 1;
n = round(b/h) + 1;
%最佳收敛因子
w = 2/(1+sin(pi/(m-1)));
%迭代次数
count = 50000;
%误差范围
delta = 10e-5;
%创建电场矩阵
Ez = zeros(m,n);
%创建磁场矩阵
Hx = zeros(m,n);
Hy = zeros(m,n);
% 内点初值
for i = 2:(m - 1)
for j = 2:(n - 1)
Ez(i, j) = 1.15;
end
end
% 边的初值
Ez( 1,:) = 0;
Ez(12,:) = 0;
Ez(:, 1) = 0;
Ez(:,12) = 0;
% Kc初值选取小于理论(Kc*h)^2的值
Kc = 0.8*(2*pi^2/m^2);
%迭代计数
k = 0;
%求解Kc式的分子
sum1 = 0;
%求解Kc式的分母
sum2 = 0;
%求解Zc所需值
sum3 = 0;
%%%%%%% 求解 %%%%%%%
while k < count
%迭代求解电场
error = 0;
for i = 2 : m - 1
for j = 2 : n - 1
temp = Ez(i,j);
Ez(i,j) = Ez(i,j) + w*((Ez(i+1,j)+Ez(i,j+1)+Ez(i-1,j)+Ez(i,j-1))/(4-(Kc*h)^2) -Ez(i,j));
error = error + abs(temp-Ez(i,j));
end
end
k = k + 1;
clear i j;
%迭代五次后求解Kc
if(rem(k,4) == 0)
for i = 2 : m - 1
for j = 2 : n - 1
%由TM边界上Ez为0，则ΔSij可以消掉，这里不予计算
sum1 = sum1 + Ez(i,j)*(Ez(i+1,j)+Ez(i,j+1)+Ez(i-1,j)+Ez(i,j-1)-4*Ez(i,j));
sum2 = sum2 + Ez(i,j)^2;
end
end
Kc = sqrt(-sum1/sum2)/h;
end
%如果总误差小于规定误差，提前跳出循环
if error < delta
break
end
end
clear i j;
%求解sum3
for i = 2 : m - 1
for j = 2 : n - 1
sum3 = sum3 + Ez(i,j)*(Ez(i+1,j)+Ez(i,j+1)+Ez(i-1,j)+Ez(i,j-1)-4*Ez(i,j))*h^2;
end
end
clear i j;
fprintf("迭代次数k = %d\n",k);
fprintf("截止波数Kc = %0.4f\n",Kc);
%绘制TM11模的Ez值热力图
figure(1);
surf(Ez);
title('TM11模 Ez');
colorbar;
%求解截止频率
c = 2.99792458e8;
fc = c*Kc/(2*pi);
fprintf("截止频率fc = %0.4fHz\n",fc);
%绘制TM11模的Ht值热力图
%真空介电常数
Epsilon = 8.854187817e-12;
%场值系数
C = -fc*2*pi*Epsilon/Kc^2;
%中心差分近似计算偏导
for i = 2:m-1
for j = 2:n-1
Hx(i,j) = 0.5*C*(Ez(i+1,j)-Ez(i-1,j))/h;
Hy(i,j) = -0.5*C*(Ez(i,j+1)-Ez(i,j-1))/h;
end
end
%绘图
[X,Y] = meshgrid(1:1:12,1:1:12);
figure(2);
quiver(X,Y,Hx,Hy);
title('TM11模 Ht');
%假设工作频率 Hz
f = 5e7;
%求解特征阻抗
E_max = max(max(Ez));
E_min = min(min(Ez));
Z0 = 120*pi;
Zc = -Z0*sqrt(1-(fc/f)^2)*(E_max-E_min)^2/(sum3*Kc^2);
fprintf("特征阻抗Zc = %0.4f\n",Zc);
6.2 直接法求解$TE_{10}$模
%基于有限差分直接法求解矩形波导中TE10模 截面内场分布、截止频率fc和特性阻抗Zc
% if you can't see the Chinese character, please change your coding format
% written by yuzhongzhibi
% https://www.drawrain.com/index.php/archives/160
%理论截止波数 K0 = 0.636363
%理论截止频率 f0 = 30.363106MHz
clear all;clc;
%波导尺寸
a = 11*pi/7;
b = 11*pi/7;
%剖分大小
h = pi/7;
%节点数量
m = round(a/h)+2;
n = round(b/h)+2;
%求解未知数个数
N = (m-1)*(n-1);
%系数矩阵
K = zeros(N);
%未知数联系矩阵
N_link = zeros(N);
%创建磁场矩阵
Hz = zeros(m-1,n-1);
%创建电场矩阵
Ex = zeros(m-1,n-1);
Ey = zeros(m-1,n-1);
%未知数编号
k = 0;
%计算联系矩阵
for j = 1 : n-1
for i = 1 : m-1
k = k + 1;
N_link(j,i) = k;
end
end
clear i j;
%填充系数矩阵
for i=1:N
for j=1:N
[im,in] = find(N_link==i);
[jm,jn] = find(N_link==j);
if i==j
K(i,j) = 4;
elseif in==jn && im==jm+1
K(i,j) = -1;
elseif in==jn && im==jm-1
K(i,j) = -1;
elseif in==jn+1 && im==jm
K(i,j) = -1;
elseif in==jn-1 && im==jm
K(i,j) = -1;
else
K(i,j)=0;
end
%边界节点处理 拐点不参与计算可省略
%下边界
if in==jn && im==jm-1 && im==1
K(i,j) = -2;
%上边界
elseif in==jn && im==jm+1 && im==(n-1)
K(i,j) = -2;
%左边界
elseif in==jn-1 && im==jm && in==1
K(i,j) = -2;
%右边界
elseif in==jn+1 && im==jm && in==(m-1)
K(i,j) = -2;
end
end
end
clear i j;
%求解特征值和特征向量
[V,D] = eig(K);
%排序处理
Kx = diag(D);
if ~issorted(Kx)
%排序
[Kx,idx] = sort(Kx);
%相同排序
V = V(:, idx);
end
Kc_all = sqrt(Kx)./h;
% 去掉特征值约为零的模
if abs(Kc_all(1)) < 0.0001
fprintf('第一特征值 ≈ %f 去除\n', Kc_all(1));
Kc_all = Kc_all(2:end);
V = V(:, 2:end);
end
%取出主模的Kc、Hz值
Kc = Kc_all(1);
fprintf("截止波数Kc = %0.4f\n",Kc);
Hz_temp = V(:,1);
%对所得特征列向量处理为矩阵
for i = 1:m-1
Hz(i,:) = Hz_temp((i-1)*(n-1)+1:i*(n-1));
end
clear i;
%绘制TE10模的hz值热力图
figure(1);
surf(Hz);
title('TE10模 Hz');
colorbar;
%求解截止频率
c = 2.99792458e8;
fc = c*Kc/(2*pi);
fprintf("截止频率fc = %0.4fHz\n",fc);
%绘制TE10模的Et值热力图
%真空介电常数
miu = pi*4e-7;
%场值系数
C = -fc*2*pi*miu/Kc^2;
%中心差分近似计算偏导
for i = 2:m-2
for j = 2:n-2
Ex(i,j) = 0.5*C*(Hz(i+1,j)-Hz(i-1,j))/h;
Ey(i,j) = -0.5*C*(Hz(i,j+1)-Hz(i,j-1))/h;
end
end
%绘图
[X,Y] = meshgrid(1:1:12,1:1:12);
figure(2);
quiver(X,Y,Ex,Ey);
title('TE10模 Et');
%假设工作频率 40MHz
f = 4e7;
%求解特征阻抗
Z0 = 120*pi;
Zc = a/b*Z0/sqrt(1-(fc/f)^2);
fprintf("特征阻抗Zc = %0.4f\n",Zc);
七、结论
由上述仿真结果与理论值对比，$K_c$、$f_c$与理论值均十分接近，TM模、TE模$f_c$相对误差均在0.34%左右，与理论结果十分契合，圆满完成了上机题目要求目标，充分验证了有限差分法的可行性与准确性。
通过本次上机实验，我深入了解到了有限差分法在求解波导场此类实际问题的应用方法，对理论部分的一些要点体会得更加深刻，也了解到了在实际工程应用中的化简处理步骤，对波导场的两种模式场的分布也有了更深入的理解。
MATLAB是一个强大的数学工具，在实际编程过程中我也遇到了不少困难，一两天都没有好的解决办法，感谢李金艳老师的悉心指导，问题均得解决，祝李老师工作顺利，身体健康。