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

warning: 这篇文章距离上次修改已过571天,其中的内容可能已经有所变动。
info:博客公式可右键缩放

一、题目分析

用差分法计算$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.jpg1.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.png2.png

对非拐点的边界上的某一未知数,若有一相邻未知数在同一行/列,且位于边界的内侧,则该相邻未知数的系数为2。

3.png3.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.png

4.2 直接法求解程序流程图

5.png5.png

五、仿真结果

5.1 $TM_{11}$模

计算结果($Z_c$为工作频率50MHz时的值):

6.png6.png

理论$K_c=0.899954$、理论$f_c=42.939916MHz$

$E_z$矩阵及曲面图:

7.png7.png

8.png8.png

与HFSS对比(同尺寸比例波导截面):
$TM_{11}$模$E_z$场分布对比

9.png9.png

10.png10.png

$TM_{11}$模$\overset{\rightarrow}{H_{t}}$场分布对比

11.png11.png

12.png12.png

5.2 $TE_{10}$模

计算结果($Z_c$为工作频率40MHz时的值):

13.png13.png

理论$K_c=0.636363$、理论$f_c=30.363106MHz$

$H_z$矩阵及曲面图:

14.png14.png

15.png15.png

与HFSS对比(同尺寸比例波导截面):
$TE_{10}$模$H_z$场分布对比(MATLAB图像蓝色部分为负值表示向内凹陷,HFSS图像左上箭头向里,两者对应。)

16.png16.png

17.png17.png

$TE_{10}$模$\overset{\rightarrow}{E_{t}}$场分布对比

18.png18.png

19.png19.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是一个强大的数学工具,在实际编程过程中我也遇到了不少困难,一两天都没有好的解决办法,感谢李金艳老师的悉心指导,问题均得解决,祝李老师工作顺利,身体健康。

添加新评论