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

•第一部分用有限元法解亥姆霍兹方程电磁波在定常态下的传播•1.1代数方程的求解•1.2强加边界条件的处理•1.3不均匀介质波导的处理•1.4对称性问题•第二部分matlab仿真实验与结果输出•2.1仿真程序•3.2输出结果10.2.4用有限元法解亥姆霍兹方程在均匀各向相同的波导结构中，电磁波的传播在定常态下服从亥姆霍兹方程。0)(20nk式中：表示电场强度矢量E的分量或磁场强度矢量H的分量；为自由空间中的波传播常数；为介质的折射率；为相对介电常数。r00/2=λkπrε=n当电磁场可以描述为在z方向具有的变化时，方程可以简化为)exp(zj0])[(2202222nkyx0n即如果代表H的分量，则的导体边界处满足第二类边界条件，如果代表的分量E，则在导体边界处满足第一类边界条件。0220[()]0tkn下面是亥姆霍兹方程对应的泛函：meJJ1)()(dyxJ])()[(21)(2220iieJ化成求解广义特征值的代数方程。对上述泛函用有限元求解，令则泛函到达极值时，应该有采用简单的三角形单元分割定义域，则有eeemjimmmjmijmjjjiimijiimmmjmijmjjjiimijiiieieieHKhhhhhhhhhkkkkkkkkkJJJ][)(4122iiiicbk)(4122jjjjcbk)(4122mmmmcbk其中：)(41jijijiijccbbkk)(41mimimiimccbbkk)(41mjmjmjjmccbbkk6mmjjiihhh12Δ==jiijhh12miimhh12mjjmhh0][eeeeHK对各单元矩阵求和，最后得到线性代数方程：0][HK或记作式中：K，H都是是对称矩阵，且H是正定的，下面阐述一下求解中遇到的问题及处理方法。TLLH1.代数方程的求解0=)Φ(•]-)([11IμLKLT对上述代数方程的求解要用到求解广义特征值的方法，求解广义特征值的方法是将正定的对称实矩阵H分解为则上述矩阵方程可以化为该方程与原矩阵方程具有相同的特征值，且转化成为求普通的矩阵特征值问题。0][HK2.强加边界条件的处理当遇到强加边界条件时，上面的代数方程中Φ的某些分量是已知的，这时对方程的求解可以采用两种处理方法。①将矩阵中已知分量所在的行、列的非主对角元素置0，而主对角元素置1，得到新的矩阵P，求解矩阵P的特征向量，将的已知分量替换成已知的数值。][HKipp0][HKi][HKi②的右端减去已知分量所在矩阵的列元素乘以已知分量的数值，将矩阵中已知分量所在的行、列的非主对角元素置0，而主对角元素置1，得到新的矩阵P，然后求解该方程。][HKi3.不均匀介质波导的处理对于存在不同介质的介质波导时，因为，折射率与三角元素所处的位置有关，所以不能直接求解广义特征值方程，而应当将并入系数矩阵中，得到新的代数方程，此时求解得到广义特征值才为介质波导对应的特征值。0][eeeeHK2220enkeneeHnk220eK4.对称性问题当研究的场域具有对称性时，问题将得到简化。例如，当研究的场具有左右对称结构时，此时可以用一半的三角元素去替代另一半的三角元素，但在仿真中发现这种方法在求解高阶模时出现了问题。例如在下面研究的介质波导中，理论上第一个高阶模应当是HE21，而如果用对称条件的话，则不会出现HE21模，而是出现了HE31模，进一步研究还是发现凡是横向为偶数的模都不会出现。这是由于在做对称处理时已经认为对称轴上的电场自动满足第二类边界条件，所以在对称轴上电场永远不会取得0值(否则电场就是恒定为0)。即使此时改边界条件为第一类边界条件也得不到正确解。如果不利用对称条件，按照一般方法得出的解中横向就会出现偶次模，且如理论所预言的第一个高阶模为HE21。所以在对称条件下，所用的三角元素可以减少，但最后的系数矩阵规模仍不能减小。在问题求解的过程中，难点是在于对研究对象的几何形状的描述和边界条件的描述，即怎样自动把物理空间映射到离散的计算空间。在程序中表现为如何标记边界结点或边界元，一种可行的办法是在描述结点或单元的数据结构中增加一个标记是否为边界结点或边界元的标志，另外还需要在结构中增加与该结点或单元相邻的结点的信息。对于规则的几何体，划分后的结点编号可以和结点在几何体上的物理位置直接建立起映射关系，所以可以大大化简运算量。另外，可以将边界描述为满足方程的曲线，如果任何一点满足方程则该点便被标记为边界点。0),(yxf),(iiyx0),(iiyxf左上角的坐标(x,y),边长a，b，最少的单元总数Ne网格单元描述和数据组织结构结点n(xn,yn)，n结点编号，(xn,yn)结点坐标三角形单元nt(it,jt,mt),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))如果nnh,则nt((a+1,b+1),(a+2,b+2),(a+1,b+2)),！上三角如果n=nh-1,则n=n-nh+1nt((a+1,b+1),(a+2,b+1),(a+2,b+2)),！下三角%该程序用有限元法来处理标量亥姆霍兹方程，其实例是脊形介质光波导中的光波%预处理functionFEM(a,h,t,b,Ne)%参数说明%a:波导半宽度，h:脊形部分波导高度，t:脊形层外波导厚度，b:匹配层的边长(考虑为正方形)，Ne三角单元个数%区域的划分：在脊形部分采用细分，在脊形部分之外采用粗分%由于波导是左右对称，所以只考虑右半边部分。将区域分成6个大区域，横向两种划分尺度，纵向三种划分尺度%各区域中三角单元数分别占总数的15%，20%，15%，15%，20%，15%，nc=1;%claddingrefractiveindexnw=3.38;%waveguiderefractiveindexns=3.17;%substraterefractiveindexw=1.55;%wavelength:umkv=2*pi*a/w;%waveconstant: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区水平方向节点数Nv2=round(h/a*Nh1);%3区垂直方向节点数Nh2=ceil((Nh1-1)*ar4/ar3+1);Nh=Nh1+Nh2-1;%对称区水平方向节点数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;%网格步不等分划分因子ifNh21%产生各个尺度的坐标x2，y1，y3s1=(r-1)*b/(r^(Nh2-1)-1)/a;x2=zeros(1,Nh2-1);x2(1)=x1(Nh1)+s1;fori=2:Nh2-1x2(i)=x2(i-1)+s1*r^(i-1);endelsex2=(a+b)/a;endifNv11s1=(1/r-1)*(b+t-h)/((1/r)^(Nv1-1)-1)/a;y1=zeros(1,Nv1-1);y1(1)=(b+t-h)/a;fori=2:Nv1-1y1(i)=y1(i-1)-s1*(1/r)^(i-1);endelsey1=(b+t)/a;endifNv1s1=(r-1)*(b-t)/(r^(Nv3-1)-1)/a;y3=zeros(1,Nv3-1);y3(1)=y2(Nv2)-s1;fori=2:Nv3-1y3(i)=y3(i-1)-s1*r^(i-1);endelsey3=(t-b)/a;end%节点矩阵准备Ni=zeros(2,Nh*Nv);%Ni向量，第一行表示横坐标，第二行表示纵坐标%区域1fori=1:Nv1-1Ni(1,(i-1)*Nh+1:(i-1)*Nh+Nh1)=x1;%横坐标赋值Ni(2,(i-1)*Nh+1:(i-1)*Nh+Nh1)=y1(i);%纵坐标赋值end%区域2fori=1:Nv1-1Ni(1,(i-1)*Nh+Nh1+1:i*Nh)=x2;%横坐标赋值Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y1(i);%纵坐标赋值end%区域3fori=Nv1:Nv1+Nv2-1Ni(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%区域4fori=Nv1:Nv1+Nv2-1Ni(1,(i-1)*Nh+Nh1+1:i*Nh)=x2;%横坐标赋值Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y2(i-Nv1+1);%纵坐标赋值end%区域5fori=Nv1+Nv2:NvNi(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%区域6fori=Nv1+Nv2:NvNi(1,(i-1)*Nh+Nh1+1:i*Nh)=x2;%横坐标赋值Ni(2,(i-1)*Nh+Nh1+1:i*Nh)=y3(i-Nv1-Nv2+1);%纵坐标赋值endti=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)*Nv,(2*Nh-1)*Nv);%常数项矩阵fori=1:2*(Nh-1)*(Nv-1)%由三角形编号得到顶点编号，对于对称结构只准备对称部分的三角单元，由其编号可以推出另一半的编号m=floor((i-1)/2/(Nh-1));n=mod(i-1,2*(Nh-1));ifn=(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);%对称点的标号(m+1)*(2*Nh-1)+(Nh-n-1)delt2=0;ifm=Nv1-2%区域1和区域2neff=nc;elseifm=Nv1+Nv2-2%区域5和区域6折射率neff=ns;elseifm=Nv1+ti-1%波导部分的区域2neff=nw;elseifn=Nh1-1neff=nc;else%区域3neff=nw;endendelse%上三角Nt(1,i)=m