规则波导的有限元分析
求解原理
上一部分我们使用一维当中的平行电容求解,现在我们开始解决规则波导的求解
-
实验原理与一维有限元分析原理相同,但是本征值问题并不是形如
的求解,而是
,除此之外,在区域离散化、基函数选取中也不相同。
- 区域离散化
将二维平面进行离散化,将其划分为若干个三角形,将三角形的顶点坐标、个数、边进行存储; - 选取插值基函数
插值基函数选取为,通过插值条件判断:
-
可得出具体插值基函数为:
- 建立并组装刚度矩阵
由下列方程建立矩阵
再通过循环进行组装矩阵 - 处理边界条件并求解
由于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波理论值' );
不规则波导的求解
求解过程
求解代码与规则波导相同,需要注意的是,我们只需要将第一步画边界时画出自己需要的边界即可,利用区域的加减等方法即可。