规则与不规则波导本征值问题有限元求解（利用Matlab pdetool）

bit_ljf

已于 2024-04-05 12:03:17 修改

阅读量1.2k

点赞数 27

分类专栏：电子信息大二下 # 电磁场理论计算与应用一文章标签：算法 matlab

于 2024-04-04 21:29:35 首次发布

本文链接：https://blog.csdn.net/m0_62655963/article/details/137380097

版权

电子信息大二下同时被 2 个专栏收录

2 篇文章 8 订阅

订阅专栏

电磁场理论计算与应用一

2 篇文章 1 订阅

订阅专栏

规则波导的有限元分析

求解原理

上一部分我们使用一维当中的平行电容求解，现在我们开始解决规则波导的求解

实验原理与一维有限元分析原理相同，但是本征值问题并不是形如的求解，而是，除此之外，在区域离散化、基函数选取中也不相同。
区域离散化
将二维平面进行离散化，将其划分为若干个三角形，将三角形的顶点坐标、个数、边进行存储；
选取插值基函数
插值基函数选取为，通过插值条件判断：
可得出具体插值基函数为：
建立并组装刚度矩阵
由下列方程建立矩阵

再通过循环进行组装矩阵
处理边界条件并求解
由于TM波边界条件为0，需要在刚度矩阵设立边界条件，即在边界处设为0，TE波的边界条件已经被包含在泛函当中，无需再次设立。
得到设立边界条件后的刚度矩阵与结果矩阵后可通过求解出，即该方程的特征值；
求解精确解
因为TE波与TM波的矩形空波导计算公式为可画出TE波与TM波前十个本征模的传播常数不进行对比
可以通过数值解与精确解的一一比对，来确保数值解的精确性。

求解代码

在使用matlab求解时，使用pdetool进行网格划分

输入之后会弹出以下页面

我们可以按照波导的形状进行绘图

在绘图之后点击三角进行网格划分

最后点击mesh，底部有Export mesh

导出三角形顶点坐标p、每个单元节点的节点编号e、与三角形单元编号的trig三个矩阵。

点击OK，可以关闭pdetool页面，将文件保存为.m格式

通过计算p与trig矩阵大小来确定刚度矩阵与结果矩阵的大小
其中在矩阵创建时可以使用稀疏矩阵，但是稀疏矩阵带来的缺点是搜索时效率下降，而稠密矩阵中存储太多0使空间利用率较差
```
Np=size(p);%p is mesh of points
Nt=size(trig);%trig is number of triangles
%e is the edge of triangle
pn=Np(2);
tn=Nt(2);
%矩阵组装中的第一步，将Ae与Be求出后组装
b=zeros(3,1);
c=zeros(3,1);
Ae=zeros(tn,3,3);
Be=zeros(tn,3,3);
```

获取顶点编号，得到矩阵

for k=1:tn
    i1=trig(1,k); i2=trig(2,k); i3=trig(3,k);
    x1=p(1,i1);   x2=p(1,i2);   x3=p(1,i3);
    y1=p(2,i1);   y2=p(2,i2);   y3=p(2,i3);
    %求得三角形面积
    S = x3 * (y1 - y2) + x1 * (y2 - y3) + x2 * (y3 - y1);
    %从数学角度来讲是要求梯度
    b(1)=y2-y3;b(2)=y3-y1;b(3)=y1-y2;
    c(1)=x3-x2;c(2)=x1-x3;c(3)=x2-x1;
    %求梯度的内积
    for i=1:3
        for j=1:3
            Ae(k,i,j)=(b(i)*b(j)+c(i)*c(j))/(S);
            Be(k,i,j)=(1+~(i-j))*S/12;
        end
    end
end

循环组装刚度矩阵与结果矩阵

%矩阵组装第二步，组装A与B
A = zeros( pn , pn );
B = zeros( pn , pn );
for k=1:tn
    for i=1:3
        for j=1:3
            %求刚度矩阵A
            A(trig(i,k),trig(j,k)) = A(trig(i,k),trig(j,k)) + Ae(k,i,j);
            B(trig(i,k),trig(j,k)) = B(trig(i,k),trig(j,k)) + Be(k,i,j);
        end
    end
end

处理边界条件并进行TE波与TM波求解

%求解TE波，TE波不需要进行边界条件处理
[Hz,kc_TE]=eigs(A,B,10,'smallestreal');
kc_te=sqrt(kc_TE);
kc_2_te =  sort( eig( B\A ) );
kc_Te = sqrt( kc_2_te );
%求解TM波，TM波需要将边界置为0
A_tm=A;
B_tm=B;
for i=1:length(e)
    A_tm(e(1,i), :) = zeros(1, pn);% 将对应行全部设为0
    A_tm(:, e(1,i)) = zeros(pn, 1);% 将对应列全部设为0
    B_tm(e(1,i), :) = zeros(1, pn);% 将对应行全部设为0
    B_tm(:, e(1,i)) = zeros(pn, 1);% 将对应列全部设为0
end
[Ez , D] = eig( pinv(B_tm) * A_tm );%求特征值和特征向量
[kc_2_tm0 , I] = sort( diag(D) );%kc_2_tm是kc方的升序，其中可能含有0；I为索引
Ez_tm = Ez(:,I); %对特征向量排序，使之对应,即为Ez
j = 1;
for i = 1:pn
    if(kc_2_tm0(i) == 0)
        continue;
    else
        kc_2_tm(j,1) = kc_2_tm0(i);    %kc_2_tm中存的是去掉0的kc方的值
        j = j+1;
    end
end

根据理论公式对举行波导进行TE波与TM波计算

%计算TE波的理论值
k=1;
for m = 0:5
    for n = 0:5
        if( m==0 && n==0 )
            continue;
        else
            kc_th_2_te( k ) = ( m*pi / 0.4572).^2 + ( n*pi / 0.2286).^2;
            k = k+1;
        end
    end     
end
kc_th_2_te = sort( kc_th_2_te );

%计算TM波的理论值
k=1;
for m = 1:5
    for n = 1:5
        kc_th_2_tm( k )=( m*pi / 0.4572).^2 + ( n*pi / 0.2286).^2;
        k = k+1;
    end     
end
kc_th_2_tm = sort( kc_th_2_tm );

最后输出图像

x = 1:10;
plot( x, kc_2_te(2:11) );
title('前10个本征模Kc平方的对比结果');
hold on;
plot( x, kc_2_tm(1:10) );
hold on;
plot( x, kc_th_2_te(1:10) );
hold on;
plot( x, kc_th_2_tm(1:10) );
legend('TE波程序计算值', 'TM波程序计算值', 'TE波理论值', 'TM波理论值' );