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

《用有限元法解亥姆赫兹方程》由会员分享，可在线阅读，更多相关《用有限元法解亥姆赫兹方程（44页珍藏版）》请在装配图网上搜索。

1、会计学1用有限元法解亥姆赫兹方程用有限元法解亥姆赫兹方程)exp(zj0)(2202222nkyx0n即即 如果如果 代表代表H的分量，则的分量，则 的导体边界处满足的导体边界处满足第二类边界条件第二类边界条件 ，如果，如果 代表的分量代表的分量E，则则 在导体边界处满足第一类边界条在导体边界处满足第一类边界条件件 。 0220()0tk n 第1页/共44页meJJ1)()(dyxJ)()(21)(2220iieJ 化成求解广义特征值的代数方程。对化成求解广义特征值的代数方程。对上述泛函用有限元求解，令上述泛函用有限元求解，令则泛函到达极值时，应该有则泛函到达极值时，应该有第2页/共44页e

2、eemjimmmjmijmjjjiimijiimmmjmijmjjjiimijiiieieieHKhhhhhhhhhkkkkkkkkkJJJ)(4122iiiicbk)(4122jjjjcbk)(4122mmmmcbk其中：其中：第3页/共44页)(41jijijiijccbbkk)(41mimimiimccbbkk)(41mjmjmjjmccbbkk6mmjjiihhh12=jiijhh12miimhh12mjjmhh第4页/共44页0eeeeHK0 HK或记作或记作式中：式中：K，H都是是对称矩阵，且都是是对称矩阵，且H是正定的，是正定的，下面阐述一下求解中遇到的问题及处理方法。下面阐述一

3、下求解中遇到的问题及处理方法。第5页/共44页TLLH0)()(11TTLILKL 对上述代数方程的求解要用到求解广义特征值的对上述代数方程的求解要用到求解广义特征值的方法，求解广义特征值的方法是将正定的对称实方法，求解广义特征值的方法是将正定的对称实矩阵矩阵H分解为分解为则上述矩阵方程可以化为则上述矩阵方程可以化为该方程与原矩阵方程具有相同的特征值，且转化成为该方程与原矩阵方程具有相同的特征值，且转化成为求普通的矩阵特征值问题。求普通的矩阵特征值问题。0 HK第6页/共44页 当遇到强加边界条件时，上面的代数方程中当遇到强加边界条件时，上面的代数方程中的的某些分量是已知的，这时对方程的求解可

4、以采用某些分量是已知的，这时对方程的求解可以采用两种处理方法。两种处理方法。将矩阵将矩阵 中已知分量所在的行、列的非主中已知分量所在的行、列的非主对角元素置对角元素置0，而主对角元素置，而主对角元素置1，得到新的矩阵，得到新的矩阵P，求解矩阵求解矩阵P的特征向量的特征向量 ，将，将 的已知分量替换的已知分量替换成已知的数值。成已知的数值。HKipp第7页/共44页0HKiHKiHKi第8页/共44页2220enkeneeHnk220eK 对于存在不同介质的介质波导时，因对于存在不同介质的介质波导时，因为为 ，折射率，折射率 与三角元素所处的与三角元素所处的位置有关，所以不能直接求解广特征值方程

5、，而应位置有关，所以不能直接求解广特征值方程，而应当将当将 并入系数矩阵并入系数矩阵 中，得到新的代数中，得到新的代数方程，此时求解得到广义特征值才为介质波导对应方程，此时求解得到广义特征值才为介质波导对应的特征值。的特征值。0eeeeHK第9页/共44页第10页/共44页第11页/共44页0),(yxf),(iiyx0),(iiyxf第12页/共44页第13页/共44页% 该程序用有限元法来处理标量亥姆霍兹方程，其实例该程序用有限元法来处理标量亥姆霍兹方程，其实例是脊形介质光波导中的光波是脊形介质光波导中的光波% 预处理预处理function FEM(a,h,t,b,Ne)% 参数说明参数说

6、明%a:波导半宽度，波导半宽度，h:脊形部分波导高度，脊形部分波导高度，t:脊形层外波导脊形层外波导厚度，厚度，b:匹配层的边长匹配层的边长(考虑为正方形考虑为正方形)，Ne三角单元个数三角单元个数%区域的划分：在脊形部分采用细分，在脊形部分之外采区域的划分：在脊形部分采用细分，在脊形部分之外采用粗分用粗分%由于波导是左右对称，所以只考虑右半边部分。将区域由于波导是左右对称，所以只考虑右半边部分。将区域分成分成6个大区域，横向两种划分尺度，纵向三种划分尺度个大区域，横向两种划分尺度，纵向三种划分尺度%各区域中三角单元数分别占总数的各区域中三角单元数分别占总数的15%，20%，15%，15%，2

7、0%，15%，第14页/共44页nc=1; %cladding refractive indexnw=3.38; %waveguide refractive indexns=3.17; %substrate refractive indexw=1.55; %wavelength:umkv=2*pi*a/w; %wave constant:1/umar1=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区水平方向节点数区水平方

8、向节点数Nv2=round(h/a*Nh1); %3区垂直方向节点数区垂直方向节点数Nh2=ceil(Nh1-1)*ar4/ar3+1);Nh=Nh1+Nh2-1; %对称区水平方向节点数对称区水平方向节点数第15页/共44页Nv1=ceil(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

9、; %网格步不等分划分因子网格步不等分划分因子if Nh21 %产生各个尺度的坐标产生各个尺度的坐标x2，y1，y3 s1=(r-1)*b/(r(Nh2-1)-1)/a; x2=zeros(1,Nh2-1); x2(1)=x1(Nh1)+s1; for i=2:Nh2-1 x2(i)=x2(i-1)+s1*r(i-1); endelse x2=(a+b)/a;end第16页/共44页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-

10、1)-s1*(1/r)(i-1); endelse y1=(b+t)/a;end第17页/共44页if Nv1 s1=(r-1)*(b-t)/(r(Nv3-1)-1)/a; y3=zeros(1,Nv3-1); y3(1)=y2(Nv2)-s1; for i=2:Nv3-1 y3(i)=y3(i-1)-s1*r(i-1); endelse y3=(t-b)/a;end第18页/共44页%节点矩阵准备节点矩阵准备 Ni=zeros(2,Nh*Nv); %Ni向量，第一行表示横坐标，第向量，第一行表示横坐标，第二行表示纵坐标二行表示纵坐标%区域区域1for i=1:Nv1-1 Ni(1,(i-1)

11、*Nh+1:(i-1)*Nh+Nh1)=x1; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+1:(i-1)*Nh+Nh1)=y1(i); %纵坐标赋值纵坐标赋值end%区域区域2for i=1:Nv1-1 Ni(1,(i-1)*Nh+Nh1+1:i*Nh)=x2; %横坐标赋值横坐标赋值 Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y1(i); %纵坐标赋值纵坐标赋值end第19页/共44页%区域区域3for i=Nv1:Nv1+Nv2-1Ni(1,(i-1)*Nh+1:(i-1)*Nh+Nh1)=x1; %横坐标赋值横坐标赋值Ni(2,(i-1)*Nh+1:(i-1)*Nh+

12、Nh1)=y2(i-Nv1+1); %纵坐标赋值纵坐标赋值end%区域区域4for 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(i-Nv1+1); %纵坐标赋值纵坐标赋值end第20页/共44页%区域区域5for 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%区域区域6fo

13、r 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第21页/共44页ti=find(y2=t/a);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)

14、*Nv,(2*Nh-1)*Nv); %常数项矩阵常数项矩阵for i=1:2*(Nh-1)*(Nv-1)%由三角形编号得到顶点编号，对于对称结构只准备对称由三角形编号得到顶点编号，对于对称结构只准备对称部分的三角单元，由其编号可以退出另一半的编号部分的三角单元，由其编号可以退出另一半的编号 m=floor(i-1)/2/(Nh-1); n=mod(i-1,2*(Nh-1);第22页/共44页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*(

15、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); %对称点的标号对称点的标号(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 第24页/共44页else %上三角上三角

16、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+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;第25页/共44页if m=Nv1+Nv2-2 neff=ns; else if m=ti+Nv1-2 n

17、eff=nw; %波导部分的区域波导部分的区域2 elseif n=Nh1-1 neff=nc; else neff=nw; end end end第26页/共44页 bi=Ni(2,Nt(2,i)-Ni(2,Nt(3,i); %bi bj=Ni(2,Nt(3,i)-Ni(2,Nt(1,i); %bj bm=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

18、*cj-bj*ci); if d_wguide=1 neff=0; end 第27页/共44页%生成有限元系数矩阵生成有限元系数矩阵 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

19、)/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,NNt3)=K(NNt2,NNt3)+(bj*bm+cj*cm)/s-(kv*neff)2*s/4; %Kjm=Kmj K(NNt3,NNt2)=K(NNt2,NNt3);第28页/共44页 H(NNt1,NNt1)=H(NNt1,NNt1)+2*s/3; %Hii H(NNt2,NN

20、t2)=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; %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; N

21、Nt2=NNt2-n-n-delt2; NNt3=NNt3-n-n-2;第29页/共44页 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*s/4; %

22、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,NNt3)=K(NNt2,NNt3)+(bj*bm+cj*cm)/s-(kv*neff)2*s/4; %Kjm=Kmj K(NNt3,NNt2)=K(NNt2,NNt3); 第30页/共44页 H(NNt1,NNt1)=H(NNt1,NNt1)+2*s/3; %Hii H(NNt2,NNt2)=H(NNt2,NNt2)+2*s

23、/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; %Him=Hmi H(NNt3,NNt1)=H(NNt1,NNt3); H(NNt2,NNt3)=H(NNt2,NNt3)+s/4; %Hjm=Hm H(NNt3,NNt2)=H(NNt2,NNt3);end第31页/共44页%求解方程求解方程if d_wguide=1 EVC,EL=eig(K,-H,cho

24、l);else EVC,EL=eig(K,H,chol);end%边界处理边界处理clear Nt; indexi,indexj=find(EL(kv*ns)2); %寻找导模寻找导模if length(indexj)0 j=1; s=input(noshou next eigenmode?,s); while(s=y|s=Y)&(jlength(indexj) fprintf(nall modes is showed) else s=input(show next eigenmode?,s); end endelse fprintf(no guide mode exists,maybe there is some error of your problem!) endfprintf(nend of this program)% copy 到运行空间到运行空间%function FEM%FEM(1.5,1.5,0.2,2,600); 第36页/共44页第37页/共44页第38页/共44页第39页/共44页第40页/共44页第41页/共44页第42页/共44页第43页/共44页