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

规则波导的有限元分析

求解原理

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

  1. 实验原理与一维有限元分析原理相同,但是本征值问题并不是形如的求解,而是,除此之外,在区域离散化、基函数选取中也不相同。

  2. 区域离散化
    将二维平面进行离散化,将其划分为若干个三角形,将三角形的顶点坐标、个数、边进行存储;
  3. 选取插值基函数
    插值基函数选取为,通过插值条件判断:
  4. 可得出具体插值基函数为:

  5. 建立并组装刚度矩阵
    由下列方程建立矩阵


    再通过循环进行组装矩阵

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

求解代码

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

输入之后会弹出以下页面

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

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

最后点击mesh,底部有Export mesh

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

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

  1. 通过计算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);
  2.  获取顶点编号,得到矩阵
    ​​​​​​​​​​​​​​
    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
  3. 循环组装刚度矩阵与结果矩阵
    %矩阵组装第二步,组装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
     
  4.  处理边界条件并进行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

  5.  根据理论公式对举行波导进行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 );
    

  6. 最后输出图像
    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波理论值' );
     

不规则波导的求解

求解过程

求解代码与规则波导相同,需要注意的是,我们只需要将第一步画边界时画出自己需要的边界即可,利用区域的加减等方法即可。 

  • 27
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值