计算电磁学实验五-矩量法实验

实验五实验五矩量法矩量法实验五实验五矩量法矩量法计算电磁学－矩量法计算电磁学－矩量法计算电磁学－矩量法计算电磁学－矩量法电子科技大学物理电子学院电子科技大学物理电子学院赖生建赖生建[e-mail]@uestc.edu.cn加权余量法加权余量法加权余量法加权余量法„„给定边值问题的场方程给定边值问题的场方程（（微分或积分方程微分或积分方程））及边及边„„给定边值问题的场方程给定边值问题的场方程（（微分或积分方程微分或积分方程））及边及边界条件统一表述为如下的界条件统一表述为如下的算子方程算子方程,LugugV=∈1()sbSuu=r1()sbSsbuq∂=r2()sbSqn∂[e-mail]离散化为矩阵离散化为矩阵„„左边等于左边等于n⎛⎞jjvvWLudVWgdV=∫∫1njiiviWLNudV=⎛⎞⎜⎟⎝⎠∑∫nn右边等于右边等于()()11nnjiiijivviiWLNudVuWLNdV===∑∑∫∫定义内积表达式定义内积表达式()(),jijiWLNdVWLN≡∫„„右边等于右边等于,jjWgdVWg≡∫„„定义内积表达式定义内积表达式()(),jijiv∫,jjvgg∫„„加权余量式加权余量式可简写成可简写成()1,,(1,2,,)nijijiuWLNWgjn===∑1i=[e-mail]离散化为矩阵离散化为矩阵上式为含上式为含个未知数个未知数的的个方程个方程以用以用矩阵的形矩阵的形„„上式为含上式为含nn个未知数个未知数的的nn个方程个方程，，可以用可以用矩阵的形矩阵的形式式来表示来表示iu{}{}{}M{}{}{}Mug=()()()WLNWLNWLN⎧⎫{}1112121222,(),(),(),(),(),()nnWLNWLNWLNWLNWLNWLNM⎧⎫⎪⎪⎪⎪=⎨⎬⎪⎪####12,(),(),()nnnnWLNWLNWLN⎪⎪⎪⎪⎩⎭11,uWg⎧⎫⎧⎫⎪⎪⎪⎪,()jijiMWLN={}{}22,,uWgug⎪⎪⎪⎪⎪⎪⎪⎪==⎨⎬⎨⎬⎪⎪⎪⎪##,nnuWg⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭[e-mail]矩量法编程矩量法编程()MLWN=矩量法编程矩量法编程„„网格剖分网格剖分,()jiijMLWN=格剖分格剖分„„基函数选择基函数选择„„权函数选择权函数选择„„编写基函数选择与方编写基函数选择与方„„编写权函数选择与基编写权函数选择与基编写基函数选择与方编写基函数选择与方程算子结合的计算程算子结合的计算编写权函数选择与基编写权函数选择与基函数内积计算函数函数内积计算函数„„矩阵填充计算矩阵填充计算[e-mail]编写离散网格数据结构编写离散网格数据结构„„编写离散点数据结构编写离散点数据结构ttStPit{„structStruPoint{„doublex;„doubley;„};„„编写网格数据结构编写网格数据结构„structStruElemental{„intid;„doublexs,xe,ys,ye;„StruPoint*pt;„};}[e-mail]编写基函数计算函数块编写基函数计算函数块定义基数类型定义基数类型„„定义基函数类型定义基函数类型„enumBaseFunType{BASFUN_PULSE,BASFUN_TRI}„„编写基函数计算的函数块编写基函数计算的函数块„doubleCmpBaseFun(BaseFunTypetype,doublex)p(ypyp,)„{„...„}[e-mail]编写矩阵元素计算函数块编写矩阵元素计算函数块„„定义权函数类型定义权函数类型„enumWeightFunType{WGT_DELTA,WGT_GALERKIN}„„矩阵元素计算类型矩阵元素计算类型矩阵元素计算类型矩阵元素计算类型„doubleCmpMatrixElement(double(*L)(double),StruElemental*elemi,StruPoint*ptj,BaseFunType)pjypbaseType,WeightFunTypeweightType)„{„doubleval;„if(weightType==WGT_DELTA){„val=L(elemi,ptj,baseType);„}„returnval;„}[e-mail]带电导体棒的电场分布带电导体棒的电场分布„„半径为半径为aa，长度为，长度为LL的细长带电导体棒，给定的细长带电导体棒，给定电位电位，，求此带电导体棒的求此带电导体棒的电荷分布电荷分布0V求此带电导体棒的求此带电导体棒的电荷分布电荷分布0()()4'r'rr'Leeqxdxφπε=−∫22()()rr'''wherexxyy−=−+−导体棒的分段模型导体棒的分段模型„„导体棒的分段模型导体棒的分段模型由边界条件由边界条件„„由边界条件由边界条件000,()rexLyaVφφ===„„基函数展开基函数展开()()''Nennqxax=Π∑1()0x'xnnlxl∈Δ⎧Π=⎨∉Δ⎩[e-mail]1ennn=∑0xnl∉Δ⎩带电导体棒的电场分布带电导体棒的电场分布„„把基函数代入边界积分方程把基函数代入边界积分方程0(1)1114'rr'NndxnndxnVadxπε−==−∑∫2222()()()rr''''xxyyxxa−=−+−=−+1n=()xxa=+„„点匹配法计算（点匹配法计算（权函数作用权函数作用））11NdL(11)001()1,4rrrr'r'NndxLmnnmndxVadxdδδπε−=−−=∫∫∑1N(11)4,1,2,...,1'rr'ndxnNedxmnnadVNxmπε=−−==∑∫„„计算电位系数（计算电位系数（MOMMOM矩阵元素矩阵元素））11''ndxndxdxdxz==∫∫[e-mail]2(1)(1)(')rr'ndxndxmmmndxdxxzxa−−−−+∫∫带电导体棒的电场分布带电导体棒的电场分布01111214NVazzzπε⎧⎫⎧⎫⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪„„矩阵方程为矩阵方程为21222204NzzzaVπε⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪=⎨⎬⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪######1204NNNNNzzzaVπε⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭⎩⎭1d2(1)1(')'ndxmnndxmzdxxxa−=−+∫()()()()22logbmbmxxxxa⎡⎤−+−−⎢⎥=⎢⎥()()2,(1),amambaxxxxawherexndxxndx⎢⎥−+−−⎣⎦==−[e-mail],(),ba带电导体棒编程带电导体棒编程„„定义导体棒的数据结构定义导体棒的数据结构„typedefstructStructWire{„doublelenx;„doublea;„intnx;„StruElemental*pElement;„double(*L)(StruElemental*,StruPoint*);()(,);„}StruMoM,*pStruMoM;[e-mail]带电导体棒编程带电导体棒编程„„导体棒的算子函数导体棒的算子函数„doubleCmpWireAlgo(StruElemental*elemi,StruPoint*ptj)pg(,pj)„{„doubleval,vb2,va2,xa,xb,a2,xm;„doubleval,vb2,va2,xa,xb,a2,xm;„xm=ptj-x;„xa=elemi-xs;„xaelemixs;„xb=elemi-xe;„a2=ptj-y-elemi-pt-y;„a2ptjyelemipty;„a2*=a2;„vb2=(xb-xm)*(xb-xm);„vb2=(xb-xm)(xb-xm);„va2=(xa-xm)*(xa-xm);„val=log((vb2+sqrt(vb2a2))/(va2+sqrt(va2a2)));„val=log((vb2+sqrt(vb2-a2))/(va2+sqrt(va2-a2)));„returnval;„}[e-mail]导电平板的静电场导电平板的静电场平静平静„„设正方形导电板，边长设正方形导电板，边长为为位于位于平面上中平面上中LL，，位于位于ZZ＝＝00平面上，中平面上，中心点如心点如图示，图示，若导电平板若导电平板电位电位试求导电板上的试求导电板上的电位电位V0V0，，试求导电板上的试求导电板上的电荷分布及电荷分布及电容电容/2/2('')LLqxy∫∫/2/2022/2/2(,)4()''(')(')LLeLLqxydxdyxxyyπεφ−−=−+−∫∫[e-mail]导电平板的静电场导电平板的静电场„„11、首先、首先分扳为分扳为NN个均匀小个均匀小块块△△SS，并选，并选基函数基函数为为首先首先分扳为分扳为个均匀小个均匀小块块并并函数函数为为分域脉冲函数分域脉冲函数。。1NS∈Δ⎧r代入边界的积分方程代入边界的积分方程11()(),()0NnennnnnSqaS=∈Δ⎧=ΠΠ=⎨∉Δ⎩∑rr'r'r'r„„代入边界的积分方程代入边界的积分方程1N∑∫∫0004()|4rSVπεφπε=„„22、、选选权函数权函数0022114''(')(')nnSnVadxdyxxyyπε==−+−∑∫∫„„22、、选选权函数权函数()()wxxyyδ=−−()()mmmwxxyyδ[e-mail]导电平板的静电场导电平板的静电场平静平静„„33、求求内积内积求求内积内积,()()()()mnmenmmenSzwLqxxyyLqdxdyδ==−−∫∫1()()''xxyydxdydxdyδ⎡⎤=−−⎢⎥∫∫∫∫1∫∫22()()(')(')nmmSSxxyydxdydxdyxxyyδ⎢⎥⎢−+−⎥⎣⎦∫∫∫∫221''(')(')nSmmdxdyxxyy=−+−∫∫[e-mail]导电平板的静电场导电平板的静电场平静平静44矩阵主角元素计算矩阵主角元素计算„„44、矩阵主角元素计算矩阵主角元素计算„„当当m=nm=n时时„IntegrantisaSingularity221''()()nnSzdxdy=∫∫„IntegrantisaSingularity22(')(')nSmmxxyy−+−∫∫2log(12)a=+πε[e-mail]导电平板的静电场导电平板的静电场平静平静221,''mnSmnzdxdy≠=∫∫22,(')(')nmnSmmyxxyy−+−∫∫SΔ22()()nmnmnSxxyyΔ≈−+−„„55、、矩阵方程矩阵方程[][][]1−[][][]zag=[][][]1nmnnNazg−=[][][]mnnnzag=1Nennqa=Π∑1n=[e-mail]导电平板编程导电平板编程„„定义导体平板的数据结构定义导体平板的数据结构„structStructPlane:{„doublelenx,leny;„intnx,ny;„StruElemental*pElement;p„double(*L)(StruPoint*);„};};[e-mail]带电导体棒编程带电导体棒编程„„导体棒的算子函数导体棒的算子函数„doubleCmpPlaneAlgo(StruElemental*elemi,StruPoint*ptj)pg(,pj)„{doubleds,xn,yn,xm,ym,val,a;„if(ptj==elemi-pt){„if(ptjelemipt){„a=elemi-xe-elemi-xs;„val=2*a*log(1+sqrt(2))/(PI*EPS0);„val2alog(1.+sqrt(2.))/(PIEP