用有限元法解亥姆赫兹方程课件

1、 第一部分第一部分 用有限元法解亥姆霍兹方程用有限元法解亥姆霍兹方程 电磁波在定常态下的传播电磁波在定常态下的传播 1.1 代数方程的求解代数方程的求解 1.2 强加边界条件的处理强加边界条件的处理 1.3 不均匀介质波导的处理不均匀介质波导的处理 1.4 对称性问题对称性问题 第二部分第二部分 matlab仿真实验与结果输出仿真实验与结果输出 2.1 仿真程序仿真程序 3.2 输出结果输出结果 10.2.4 用有限元法解亥姆霍兹方程用有限元法解亥姆霍兹方程 在均匀各向相同的波导结构中，电磁波的传播在均匀各向相同的波导结构中，电磁波的传播 在定常态下服从亥姆霍兹方程。在定常态下服从亥姆霍兹方程

2、。 0)( 2 0 nk 式中：式中： 表示电场强度矢量表示电场强度矢量E的分量或磁场强的分量或磁场强 度矢量度矢量H的分量；的分量； 为自由空间中的波为自由空间中的波 传播常数；传播常数； 为介质的折射率；为介质的折射率； 为相对为相对 介电常数。介电常数。 r 00 /2=k r =n 当电磁场可以描述为在当电磁场可以描述为在z方向具有方向具有 的变化时，方程可以简化为的变化时，方程可以简化为 )exp(zj 0)( 22 0 2 2 2 2 nk yx 0 n 即即 如果如果 代表代表H的分量，则的分量，则 的导体边界处满的导体边界处满 足第二类边界条件足第二类边界条件 ，如果，如果 代

3、表的分量代表的分量 E，则，则 在导体边界处满足第一类边界条在导体边界处满足第一类边界条 件件 。 0 22 0 ()0 t k n 下面是亥姆霍兹方程对应的泛函：下面是亥姆霍兹方程对应的泛函： m e JJ 1 )()( d yx J)()( 2 1 )( 222 0 i i e J 化成求解广义特征值的代数方程。对上化成求解广义特征值的代数方程。对上 述泛函用有限元求解，令述泛函用有限元求解，令 则泛函到达极值时，应该有则泛函到达极值时，应该有 采用简单的三角形单元分割定义域，则有采用简单的三角形单元分割定义域，则有 eee m j i mmmjmi jmjjji imijii mmmjm

4、i jmjjji imijii i e i e i e HK hhh hhh hhh kkk kkk kkk J J J )( 4 1 22 iiii cbk )( 4 1 22 jjjj cbk )( 4 1 22 mmmm cbk 其中：其中： )( 4 1 jijijiij ccbbkk )( 4 1 mimimiim ccbbkk )( 4 1 mjmjmjjm ccbbkk 6 mmjjii hhh 12 = jiij hh 12 miim hh 12 mjjm hh 0 e eee HK 对各单元矩阵求和，最后得到线性代数方程：对各单元矩阵求和，最后得到线性代数方程： 0 HK 或

5、记作或记作 式中：式中：K，H都是是对称矩阵，且都是是对称矩阵，且H是正定的，是正定的， 下面阐述一下求解中遇到的问题及处理方法。下面阐述一下求解中遇到的问题及处理方法。 T LLH 1.代数方程的求解代数方程的求解 0=)( -)( 11 ILKL T 对上述代数方程的求解要用到求解对上述代数方程的求解要用到求解广义广义特征值特征值 的方法，求解广义特征值的方法是将正定的对称的方法，求解广义特征值的方法是将正定的对称 实矩阵实矩阵H分解为分解为 则上述矩阵方程可以化为则上述矩阵方程可以化为 该方程与原矩阵方程具有相同的特征值，且转化成为该方程与原矩阵方程具有相同的特征值，且转化成为 求求普通

6、普通的矩阵特征值问题。的矩阵特征值问题。 0 HK 2.强加边界条件的处理强加边界条件的处理 当遇到强加边界条件时，上面的代数方程中当遇到强加边界条件时，上面的代数方程中 的某些分量是已知的，这时对方程的求解可以采的某些分量是已知的，这时对方程的求解可以采 用两种处理方法。用两种处理方法。 将将矩阵矩阵 中已知分量所在的行、列的非主中已知分量所在的行、列的非主 对角元素置对角元素置0，而主对角元素置，而主对角元素置1，得到新的矩阵，得到新的矩阵P， 求解矩阵求解矩阵P的特征向量的特征向量 ，将，将 的已知分量替换的已知分量替换 成已知的数值。成已知的数值。 HK i p p 0HK i HK

7、i 的右端减去已知分量所在的右端减去已知分量所在 矩阵矩阵 的列元素乘以已知分量的数值，将矩的列元素乘以已知分量的数值，将矩 阵阵 中已知分量所在的行、列的非主对角元中已知分量所在的行、列的非主对角元 素置素置0，而主对角元素置，而主对角元素置1，得到新的矩阵，得到新的矩阵P，然后，然后 求解该方程。求解该方程。 HK i 3.不均匀介质波导的处理不均匀介质波导的处理 对于存在不同介质的介质波导时，因对于存在不同介质的介质波导时，因 为为 ，折射率，折射率 与三角元素所处的与三角元素所处的 位置有关，所以不能直接求解广义特征值方程，而位置有关，所以不能直接求解广义特征值方程，而 应当将应当将

8、并入系数矩阵并入系数矩阵 中，得到新的代中，得到新的代 数方程，此时求解得到广义特征值才为介质波导对数方程，此时求解得到广义特征值才为介质波导对 应的特征值。应的特征值。 0 e eee HK 222 0 e nke n ee Hnk 22 0 e K 4.对称性问题对称性问题 当研究的场域具有对称性时，问题将得到简化。当研究的场域具有对称性时，问题将得到简化。 例如，当研究的场具有左右对称结构时，此时可以用例如，当研究的场具有左右对称结构时，此时可以用 一半的三角元素去替代另一半的三角元素，但在仿真一半的三角元素去替代另一半的三角元素，但在仿真 中发现这种方法在求解高阶模时出现了问题。中发现

9、这种方法在求解高阶模时出现了问题。 例如在下面研究的介质波导中，理论上第一个高例如在下面研究的介质波导中，理论上第一个高 阶模应当是阶模应当是HE21，而如果用对称条件的话，则不会出，而如果用对称条件的话，则不会出 现现HE21模，而是出现了模，而是出现了HE31模，进一步研究还是发模，进一步研究还是发 现凡是现凡是横向为偶数的模横向为偶数的模都不会出现。这是由于在做对都不会出现。这是由于在做对 称处理时已经认为对称轴上的电场自动满足第二类边称处理时已经认为对称轴上的电场自动满足第二类边 界条件，所以在对称轴上电场永远不会取得界条件，所以在对称轴上电场永远不会取得0值值(否则否则 电场就是恒定

10、为电场就是恒定为0)。 即使此时改边界条件为第一类边界条件也得不到正确解。如即使此时改边界条件为第一类边界条件也得不到正确解。如 果不利用对称条件，按照一般方法得出的解中横向就会出现偶次果不利用对称条件，按照一般方法得出的解中横向就会出现偶次 模，且如理论所预言的第一个高阶模为模，且如理论所预言的第一个高阶模为HE21。所以在对称条件下，。所以在对称条件下， 所用的三角元素可以减少，但最后的所用的三角元素可以减少，但最后的系数矩阵规模仍不能减小系数矩阵规模仍不能减小。 在问题求解的过程中，难点是在于对研究对象的在问题求解的过程中，难点是在于对研究对象的几何形状的几何形状的 描述和边界条件的描述

11、描述和边界条件的描述，即怎样自动把，即怎样自动把物理空间映射到离散的计物理空间映射到离散的计 算空间算空间。在程序中表现为如何标记。在程序中表现为如何标记边界结点或边界元边界结点或边界元，一种可行，一种可行 的办法是在描述结点或单元的的办法是在描述结点或单元的数据结构数据结构中增加一个中增加一个标记标记是否为边是否为边 界结点或边界元的标志，另外还需要在结构中增加与该结点或单界结点或边界元的标志，另外还需要在结构中增加与该结点或单 元相邻的结点的信息。元相邻的结点的信息。 对于规则的几何体，划分后的结点编号可以和对于规则的几何体，划分后的结点编号可以和 结点在几何体上的物理位置直接建立起映射关

12、系，结点在几何体上的物理位置直接建立起映射关系， 所以可以大大化简运算量。另外，可以将边界描述所以可以大大化简运算量。另外，可以将边界描述 为满足方程为满足方程 的曲线，如果任何一点的曲线，如果任何一点 满足方程满足方程 则该点便被标记为则该点便被标记为边界点边界点。 0),(yxf ),( ii yx 0),( ii yxf 左上角的坐标左上角的坐标(x,y),边长边长a，b，最少的单元总数，最少的单元总数Ne 网格单元描述和数据组织结构网格单元描述和数据组织结构 结点结点n(xn,yn)，n 结点编号结点编号，(xn,yn)结点坐标结点坐标 三角形单元三角形单元nt(it,jt,mt),

13、nt三角形编号三角形编号，(it,jt,mt)定点编定点编 号号。左上角的。左上角的顶点顶点映射为映射为1号结点，按先行后列反的顺序号结点，按先行后列反的顺序 编号，编号，三角形的编号三角形的编号顺序也与之相同。顺序也与之相同。 于是于是结点编号结点编号和和行列号行列号的关系为的关系为 n=(nh-1) *(nh+nv)， nh为一行内的节点数，为一行内的节点数， nv为一列内的节点数，设三角形为一列内的节点数，设三角形 的编号为的编号为nt，则它对应的结点编号为，则它对应的结点编号为 m=floor(nt-1)/2*(nh-1), n=mod(nt-1/2*(nh-1) 如果如果n=nh-1

14、,则则n=n-nh+1 nt(a+1,b+1), (a+2,b+1), (a+2,b+2), ！下三角！下三角 % 该程序用有限元法来处理标量亥姆霍兹方程，其实例该程序用有限元法来处理标量亥姆霍兹方程，其实例 是是脊形介质光波导脊形介质光波导中的光波中的光波 % 预处理预处理 function FEM(a,h,t,b,Ne) % 参数说明参数说明 %a:波导半宽度，波导半宽度，h:脊形部分波导高度，脊形部分波导高度，t:脊形层外波导脊形层外波导 厚度，厚度，b:匹配层的边长匹配层的边长(考虑为正方形考虑为正方形)，Ne三角单元个数三角单元个数 %区域的划分：在脊形部分采用细分，在脊形部分之外采

15、区域的划分：在脊形部分采用细分，在脊形部分之外采 用粗分用粗分 %由于波导是左右对称，所以只考虑右半边部分。将区域由于波导是左右对称，所以只考虑右半边部分。将区域 分成分成6个大区域，横向两种划分尺度，纵向三种划分尺度个大区域，横向两种划分尺度，纵向三种划分尺度 %各区域中三角单元数分别占总数的各区域中三角单元数分别占总数的15%，20%，15%， 15%，20%，15%， nc=1; %cladding refractive index nw=3.38; %waveguide refractive index ns=3.17; %substrate refractive index w=1.

16、55; %wavelength:um kv=2*pi*a/w; %wave constant:1/um ar1=0.15; ar3=0.15; ar4=0.15; ar5=0.20; d_wguide=1; Nh1=(1+a/h)+(1+a/h)2+4*(Ne*ar3/2-1)*a/h)0.5)*0.5; Nh1=round(Nh1); %3区水平方向节点数区水平方向节点数 Nv2=round(h/a*Nh1); %3区垂直方向节点数区垂直方向节点数 Nh2=ceil(Nh1-1)*ar4/ar3+1); Nh=Nh1+Nh2-1; %对称区水平方向节点数对称区水平方向节点数 Nv1=ceil

17、(Nh2-1)*ar1/ar3+1); %1、2区垂直方向节点数区垂直方向节点数 Nv3=ceil(Nh1-1)*ar5/ar3+1); %5、6区垂直方向节点数区垂直方向节点数 Nv=Nv1+Nv2+Nv3-2; %垂直方向的总节点数垂直方向的总节点数 x1=linspace(0,a,Nh1)/a; y2=linspace(h,0,Nv2)/a; r=1.0001; %网格步不等分划分因子网格步不等分划分因子 if Nh21 %产生各个尺度的坐标产生各个尺度的坐标x2，y1，y3 s1=(r-1)*b/(r(Nh2-1)-1)/a; x2=zeros(1,Nh2-1); x2(1)=x1(N

18、h1)+s1; for i=2:Nh2-1 x2(i)=x2(i-1)+s1*r(i-1); end else x2=(a+b)/a; end if Nv11 s1=(1/r-1)*(b+t-h)/(1/r)(Nv1-1)-1)/a; y1=zeros(1,Nv1-1); y1(1)=(b+t-h)/a; for i=2:Nv1-1 y1(i)=y1(i-1)-s1*(1/r)(i-1); end else y1=(b+t)/a; end if Nv1 s1=(r-1)*(b-t)/(r(Nv3-1)-1)/a; y3=zeros(1,Nv3-1); y3(1)=y2(Nv2)-s1; for

19、 i=2:Nv3-1 y3(i)=y3(i-1)-s1*r(i-1); end else y3=(t-b)/a; end %节点矩阵准备节点矩阵准备 Ni=zeros(2,Nh*Nv); %Ni向量，第一行表示横坐标，第向量，第一行表示横坐标，第 二行表示纵坐标二行表示纵坐标 %区域区域1 for i=1:Nv1-1 Ni(1,(i-1)*Nh+1:(i-1)*Nh+Nh1)=x1; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+1:(i-1)*Nh+Nh1)=y1(i); %纵坐标赋值纵坐标赋值 end %区域区域2 for i=1:Nv1-1 Ni(1,(i-1)*Nh+Nh1+1:

20、i*Nh)=x2; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y1(i); %纵坐标赋值纵坐标赋值 end %区域区域3 for i=Nv1:Nv1+Nv2-1 Ni(1,(i-1)*Nh+1:(i-1)*Nh+Nh1)=x1; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+1:(i-1)*Nh+Nh1)=y2(i-Nv1+1); %纵坐标赋值纵坐标赋值 end %区域区域4 for i=Nv1:Nv1+Nv2-1 Ni(1,(i-1)*Nh+Nh1+1:i*Nh)=x2; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y2

21、(i-Nv1+1); %纵坐标赋值纵坐标赋值 end %区域区域5 for i=Nv1+Nv2:Nv Ni(1,(i-1)*Nh+1:(i-1)*Nh+Nh1)=x1; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+1:(i-1)*Nh+Nh1)=y3(i-Nv1-Nv2+1); %纵坐标赋值纵坐标赋值 end %区域区域6 for i=Nv1+Nv2:Nv Ni(1,(i-1)*Nh+Nh1+1:i*Nh)=x2; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y3(i-Nv1-Nv2+1); %纵坐标赋值纵坐标赋值 end ti=find(y2=t/a)

22、; ti=ti(length(ti); % 三角形单元矩阵准备三角形单元矩阵准备 Nt=zeros(3,(Nh-1)*(Nv-1); %前三行表示逆时针方前三行表示逆时针方 向排列的三个顶点的标号向排列的三个顶点的标号 K=zeros(2*Nh-1)*Nv,(2*Nh-1)*Nv); %刚度矩阵刚度矩阵 H=zeros(2*Nh-1)*Nv,(2*Nh-1)*Nv); %常数项矩阵常数项矩阵 for i=1:2*(Nh-1)*(Nv-1) %由三角形编号得到顶点编号，对于对称结构只准备对称由三角形编号得到顶点编号，对于对称结构只准备对称 部分的三角单元，由其编号可以部分的三角单元，由其编号可以

23、推推出另一半的编号出另一半的编号 m=floor(i-1)/2/(Nh-1); n=mod(i-1,2*(Nh-1); if n=(Nh-1); %下三角形下三角形 n=n-Nh+1; Nt(1,i)=m*Nh+n+1; Nt(2,i)=(m+1)*Nh+n+1; Nt(3,i)=(m+1)*Nh+n+2; NNt1=m*(2*Nh-1)+(Nh+n); %对称点的标号对称点的标号m*(2*Nh-1)+(Nh-n) NNt2=(m+1)*(2*Nh-1)+(Nh+n); %对称点的标号对称点的标号(m+1)*(2*Nh-1)+(Nh-n) NNt3=(m+1)*(2*Nh-1)+(Nh+n+1

24、); %对称点的标号对称点的标号(m+1)*(2*Nh-1)+(Nh-n-1) delt2=0; if m=Nv1+Nv2-2 %区域区域5和区域和区域6 折射率折射率 neff=ns; else if m=Nh1-1 neff=nc; else %区域区域3 neff=nw; end end else %上三角上三角 Nt(1,i)=m*Nh+n+1; Nt(2,i)=(m+1)*Nh+n+2; Nt(3,i)=m*Nh+n+2; NNt1=m*(2*Nh-1)+(Nh+n); %对称点的标号对称点的标号 m*(2*Nh-1)+(Nh-n) NNt2=(m+1)*(2*Nh-1)+(Nh+n

25、+1); %对称点的标号对称点的标号 (m+1)*(2*Nh-1)+(Nh-n-1) NNt3=m*(2*Nh-1)+(Nh+n+1); %对称点的标号对称点的标号 m*(2*Nh-1)+(Nh-n-1) delt2=2; if m=Nv1+Nv2-2 neff=ns; else if m=ti+Nv1-2 neff=nw; %波导部分的区域波导部分的区域2 elseif n=Nh1-1 neff=nc; else neff=nw; end end end bi=Ni(2,Nt(2,i)-Ni(2,Nt(3,i); %bi bj=Ni(2,Nt(3,i)-Ni(2,Nt(1,i); %bj b

26、m=Ni(2,Nt(1,i)-Ni(2,Nt(2,i); %bm ci=Ni(1,Nt(2,i)-Ni(1,Nt(3,i); %ci cj=Ni(1,Nt(3,i)-Ni(1,Nt(1,i); %cj cm=Ni(1,Nt(1,i)-Ni(1,Nt(2,i); %cm s=0.5*abs(bi*cj-bj*ci); if d_wguide=1 neff=0; end %生成有限元系数矩阵生成有限元系数矩阵 K(NNt1,NNt1)=K(NNt1,NNt1)+(bi*bi+ci*ci)/s- 2*(kv*neff)2*s/3; %Kii K(NNt2,NNt2)=K(NNt2,NNt2)+(bj

27、*bj+cj*cj)/s- 2*(kv*neff)2*s/3; %Kjj K(NNt3,NNt3)=K(NNt3,NNt3)+(bm*bm+cm*cm)/s- 2*(kv*neff)2*s/3; %Kmm K(NNt1,NNt2)=K(NNt1,NNt2)+(bj*bj+cj*cj)/s- (kv*neff)2*s/4; %Kij=Kji K(NNt2,NNt1)=K(NNt1,NNt2); K(NNt1,NNt3)=K(NNt1,NNt3)+(bi*bm+ci*cm)/s- (kv*neff)2*s/4; %Kim=Kmi K(NNt3,NNt1)=K(NNt1,NNt3); K(NNt2,

28、NNt3)=K(NNt2,NNt3)+(bj*bm+cj*cm)/s- (kv*neff)2*s/4; %Kjm=Kmj K(NNt3,NNt2)=K(NNt2,NNt3); H(NNt1,NNt1)=H(NNt1,NNt1)+2*s/3; %Hii H(NNt2,NNt2)=H(NNt2,NNt2)+2*s/3; %Hjj H(NNt3,NNt3)=H(NNt3,NNt3)+2*s/3; %Hmm H(NNt1,NNt2)=H(NNt1,NNt2)+s/4; %Hij=Hji H(NNt2,NNt1)=H(NNt1,NNt2); H(NNt1,NNt3)=H(NNt1,NNt3)+s/4;

29、%Him=Hmi H(NNt3,NNt1)=H(NNt1,NNt3); H(NNt2,NNt3)=H(NNt2,NNt3)+s/4; %Hjm=Hmj H(NNt3,NNt2)=H(NNt2,NNt3); %对对称点的操作对对称点的操作 NNt1=NNt1-n-n; NNt2=NNt2-n-n-delt2; NNt3=NNt3-n-n-2; K(NNt1,NNt1)=K(NNt1,NNt1)+(bi*bi+ci*ci)/s- 2*(kv*neff)2*s/3; %Kii K(NNt2,NNt2)=K(NNt2,NNt2)+(bj*bj+cj*cj)/s- 2*(kv*neff)2*s/3; %Kjj K(NNt3,NNt3)=K(NNt3,NNt3)+(bm*bm+cm*cm)/s- 2*(kv*neff)2*s/3; %Kmm K(NNt1,NNt2)=K(NNt1,NNt2)+(bj*bj+cj*cj)/s- (kv*neff)2