电磁场数值计算上机题报告

电磁场数值计算上机题报告第一题计算长直接地金属槽中的电场分布。金属槽横截面如图1所示，其侧壁与底面电位均为零，顶盖电位相对值为10。槽内电位函数满足拉普拉斯方程。计算槽内电位分布。要求：（1）先用正方形网格粗分，每边取4个网格计算，取不同的松弛因子，比较其收敛速度。取计算精度为千分之一。（2）划分网格加倍，计算电位分布，并与上面计算结果比较。（一）建立问题的数学物理模型首先列出方程及其边界条件槽内的电位满足二维的拉普拉斯方程：仞平仞平八2曜相十节=°其中x,y的范围是：0<x<a,0<y<a边界条件是：里lx=0*L=a=0中=0伸=10y=0y=a图1（二）算法设计及其实现在本题中，因为区域为正方形区域，网格采用正方形网格，每边四个网格，因此，每边要有5个结点，网格数m=n=4，比较少，不能用课本中的公式计算收敛因子，取收敛因子为a=1.4，计算的程序的流程为：选取计算的场域，并划分网格，网格划分如右图2所示：右图表示网格的划分，共16个网格，一共划分了25个结点，每个节点用相应的下标（i,j）来表示，对应的电位为u（i,j）根据题意，边界条件的处理如下:u(i,j=u(，日)un=,i)

u(k,j=10$2k<n—1)

其中n其中n为一行对应的节点数，根据题意，这里n=5②用u(i,j)表示节点的电位，设经过第n次迭代之后的结点电位用U(n)来表示，则超松弛迭代法的差分格式(记u—u(i,j))为：n(+)U+n(+1U+n4—：1U)n()+，i1—1，j—，i1，ijU(n+1户Un(书(Un(十Ui,jij4+1,j边界条件前面已经给出给各个节点赋初值，对于非边界上的点(对于边界点的值前面已经赋过初值)，如下U(i,j)="^(j-1)n—1迭代计算，直到已经满足精度条件为止，这里精度为0.0001,最后输出计算的结果，输出结果见生成的数据文件shuju.txt中图3计算框图计算框图如图3所示

用fortran90语言编写计算机程序，计算各点的电位，程序清单如下:图3计算框图programhomeworklimplicitnoneintegeri,j,k,n!定义动态数组!每行的节点数!迭代因子!定义动态数组!每行的节点数!迭代因子!定义动态数组的大小，分配内存reala,b,c,eintegermn=5e=1.4allocate(u(n,n))doi=2,n-1doj=2,n-1u(i,j)=(10.00/(n-1))*(i-1)enddoenddo!给除了边界点以外的点赋初值doi=1,nu(i,1)=0u(i,n)=0enddodoj=2,n-1u(1,j)=0u(n,j)=10enddo!边界条件k=1!k是判断误差是否达标的变量m=0!m是表示迭代次数的变量do!迭代过程if(k==0)exit!当所有节点满足精度要求时，跳出循环，迭代终止doi=2,n-1doj=2,n-1k=0a=u(i,j)b=u(i,j)+e/4*(u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,jT)-4*u(i,j))u(i,j)=bif(abs(a-b)/b>0.0001)thenk=k+1endifenddoenddom=m+1!m为迭代的次数enddo!将计算结果输出，在shuju.txt中查找open(1,file='shuju.txt')write(1,*)'迭代次数为:'，mwrite(1,*)'各节点的点位为:’doi=n,1,-1write(1,'(5f7.4)')(u(i,j),j=1,n)enddoclose(1)!输出数据结果，结果位于生成的shuju.txt文件中endprogramhomework1计算结果及数据分析当网格数目为4（节点数目为5）,加速收敛因子为1.4，计算精度为0.0001时，计算出各结点电位如下所示（对不同的收敛因子的比较后面进行）迭代次数为：10各节点的点位为:0.000010.000010.000010.00000.00000.00004.28625.26794.28580.00000.00001.87502.50061.87500.00000.00000.71390.98240.71460.00000.00000.00000.00000.00000.0000丛生成的数据来看，各节点的电位值是对称的（当然，这里是在做数值计算，不可能完全对称，只能是近似对称的），这点不难从理论上进行分析得到。为了使计算结果更加精确，将网格进一步细分，网格加倍之后（此时节点数为9），计算出的结果为：迭代次数为:13各节点的点位为：0.000010.000010.000010.000010.000010.000010.000010.00000.00000.00004.82836.61407.30917.49607.30776.61184.82640.00000.00002.69894.31905.12705.36815.12364.31372.69440.00000.00001.64542.83443.51123.72583.50552.82561.63790.00000.00001.04271.85782.35472.51712.34681.84561.03250.00000.00000.65921.19181.52711.63761.51781.17730.64740.00000.00000.39240.71190.91590.98240.90680.69780.38040.00000.00000.18020.33400.43060.46140.42470.32540.17650.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000为了便于比较，特地将与第一种情形相同时的结点的电位值标记为红颜色的，通过比较不难发现，当网格数比较少时，相应得到的结果会偏小，这里，当网格数加倍时，对于同一个节点的电位值，最大的偏差有0.1V，相对偏差在2%左右。

为了进一步反应电位的变化情况，将网格进一步细分，共有3232=1024个网格，精确度也取为0.000001，将得到的数据保存，利用matlab软件读取这些数据，并作图，下图是用matlab软件画的反应电位函数①=中（乙y）的三维图形。0010).•”0-..■■■6?：■?，•-0JL-40-图4电位函数①二①^,y）00用一系列的平面O=c（常数），去截取曲面①二①^,y），得到等势线近似的分布如下，为了简便起见，这里只画出10条等势线（间隔△中=1V），知道等势线的分布，过任意点作等势线的垂线便得到电场线的分布。如图5所示77TTT70To"TTTTTTo"o.o图5等势面分布示意图不难发现，电势高的地方，等势面分布密集，电势梯度大，因而电场也较强。相反，电势较低的地方电场强度也比较弱。最后，讨论收敛因子a对迭代次数的影响下表中的数据是当每一行的网格数为4,精确度取为0.000001时不同的a值的情况下的迭代次数。收敛因子1.01.11.21.31.4迭代次数191491112收敛因子1.51.61.71.81.9迭代次数20283850102从上表可以看出，当收敛因子不同时，迭代的次数也不一样，但是，对最后所得到的结果几乎没任何影响（关于这点，调试上述程序进行验证）从上表可以发现，最佳的收敛因子应该在1.1与1.3之间，对这个区间内的数据每隔0.01进行细分，并统计相关次数，得到最佳收敛因子为1.18，此时总共迭代次数为8次。由于收敛因子一般是认为来选择的，怎么选择收敛因子，才能得到较快的收敛速度呢？一般的，这与网格的数目有关系，当一行网格数为8时，仿照上面找出最佳收敛因子的办法，最佳收敛因子为1.46,一般而言，网格数越多，相应的收敛因子可以大一些，当网格数很大时，可以用下式计算收敛因子—（正方形区域）1+sina=〈n2一兀巳二（矩形区域）、\m2n2另外，不同的初始值也会对迭代的次数有影响，但是，初始值完全是人为的给定的，我们在赋初始值时，采用接近实际电位的一组值来赋值，会减少迭代的次数。（四）其他讨论本题也可以利用对称型，取左半区域进行计算，右半区域的电势值与左半区域是对称的甲3,y）=甲（a-x,y），这时，相应的定解条件为：TOC\o"1-5"\h\zd2中d2中苛+及=°x，ye（0,a））中x=0,0<y<aa0y=0,0<x<<2中=10,ay=a,0<x<2=0x=2,0<y<a关于这种算法的程序及执行结果见本文附录中。

第二题（课本52页第三题）（一）物理模型和数学模型的建立如图6所示，根据题意，忽略曲度效应，假设定子和转子的磁导率均为无穷大，关于磁位的定解条件和边界条件为：君+藉=0（（x，y）a）’甲、TO。虬L竺=竺=0、dxABdxEF其中，区域。=ABCDEF各边的长度为：AFAfgIFBCCB'g=8AFAfgIFBCCB'（二）算法的设计和实现划分网格，确定节点，在这里取网格间距为1，仍然采用正方形的网格，网格的划分如右图7所示，由于网格是不规则的图形，但可以明显分为两部分，每个部分都是矩形，因此，分别对这两个网格的节点进行标号，记为（i,j），相应的磁位分别记为u（i,j）。边界条件的处理和以及各节点的五差分格式的确定。节点磁位用表示。结点标号的范围是：1<i<751<j<12a、先讨论边界条件的处理对于AF和BCDE两条边的边界条件，属于第一类边界条件，处理起来比较简单，直接根据题意，可知：u(75,j)=0(8<j<12)u(i,7)=0(9<i<74)'u(9,j)=0(1<j<7)u(1,j)=100(1<j<12)对于位于对称线EF上的节点，属于第二类边界条件，超松弛的五点差分迭代公式为：aU(n+1=un(〜(2u+(+)u+(+1un(4)un)((2<i<75;j=12)(i,j)(j)4—i(j,1—i(j1+)i顶1,)ij(,)同理，对于位于边界线AB上的节点，超松弛的五点差分迭代公式为：aU(n+1=Un(+(2Un(+Un+1)Un—(4Un)(<2<i=；j1)(i,j)(j)4+i(j,1-)i(j+,)i顶1,)ij(,)以上是边界条件的处理。b对于矩形1内和矩形2内的节点，五点差分迭代格式如下：aU(n+1)——U(n)+—(U(n+1)+U(n)+U(n+1)+U(n)—4u(n))(i,j)(i,j)4(i—1,j)(i+1,j)(i,j—1)(i,j+1)(i,j)以上是迭代的算法和及相关边界条件的处理。具体计算时，根据不同节点的属性，选择正确的迭代公式即可。本题中收敛因子取为1.4给节点磁位赋初值给各个节点赋初值，这一步在这里是比较难做的，原因是边界比较复杂，一旦当节点的初值赋的不好，会大大影响程序执行的效率，下面采用半线性半插值的方法赋初值，具体的办法如下：我们把求解域上面那个长为11宽为8的矩形几位矩形1，下面那个长为5宽为66的矩形记为矩形2，对于矩形1内的节点，采用线性赋值的方法，12)u(i,j=10料<伊<j<19—112)对于矩形1内地9行的节点，位于边界上的节点磁位为0，不位于边界上的节点认为其磁位与地八行同一列的磁位值的一半，即：u(9,j=u(8j<8<)12对于矩形2内的节点，也采用线性赋值的方法，边界值分别为0和u(9,j)，这时我们有：U(i,j)=U(9,j)^7^(10<i<74,8<j<12)75—9对于非求解区域内的节点，其位函数值统统认为是0，即u(i,j=01Oi<75<j96编写程序，利用上面讲的迭代算法进行计算，当计算结果达到一定的要求后输出计算的结果。程序计算的框图如图8所示

图8程序框图⑥用fortran90编写的计算机程序清单如下:programhomework2implicitnonereala,b,crealdxrealeintegeri,j,pintegerm,n,g,h,k,lreal,allocatable::u(:,:)dx=1!长度步长g=66/dxh=8/dx+1k=6/dx+1l=5/dxm=g+hn=k+lallocate(u(m,n))doj=1,nu(1,j)=100enddo!边界条件doi=2,h-1doj=1,nu(i,j)=100.0*(h-i)/(h-1)!对求解区域内的点赋初始值enddoenddodoj=1,nif(j<=k)thenu(h,j)=0!边界条件elseu(h,j)=0.5*u(h-1,j)!边界内的点赋初值endifenddodoi=h+1,m!边界条件u(i,k)=0enddodoj=k+1,n!边界条件u(m,j)=0enddodoi=h+1,m-1doj=k+1,nu(i,j)=u(h,j)*(m-i)/(m-h)!求解区域内的点赋初值enddoenddodoi=h+1,mdoj=1,ku(i,j)=0!将求解区域以外的点的磁位赋为0enddoenddoe=1.4!超松弛因子p=1!p是迭代精度是否达标的指标，p=0时，精度达标doif(p==0)exitp=0doi=2,h-1doj=1,na=u(i,j)!迭代过程if(j==1)thenb=e/4*(2*u(i,j+1)+u(i-1,j)+u(i+1,j)-4*u(i,j))elseif(j==n)thenb=e/4*(2*u(i,j-1)+u(i-1,j)+u(i+1,j)-4*u(i,j))elseb=e/4*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)-4*u(i,j))endifb=b+au(i,j)=bif(abs(b-a)/b>0.0001)thenp=p+1endifenddoenddodoi=h,m-1!迭代过程doj=k+1,na=u(i,j)if(j==n)thenb=e/4*(2*u(i,j-1)+u(i-1,j)+u(i+1,j)-4*u(i,j))elseb=e/4*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)-4*u(i,j))

endifb=b+au(i,j)=bif(abs(b-a)/b>0.0001)thenp=p+1endifenddoenddoenddoopen(1,file='shuju.txt')!将结果保存入文件中doi=1,h-1write(1,'(12f8.4)')(u(i,j),j=1,n)enddodoi=h,mwrite(1,'(48x,6f8.4)')(u(i,j),j=k,n)enddoclose(1)endprogramhomework2(三)计算结果与数据分析程序运行结果如下：100.0000100.0000100.0000100.0000100.0000100.0000100.0000100.0000100.0000100.0000100.0000100.000087.805687.824687.882587.981588.123688.307488.524688.757588.979789.162989.283289.325075.573175.610275.723875.920076.205576.581277.033577.525677.998578.388978.644778.733663.266563.319363.482563.769364.197364.778665.502666.312867.100067.749368.173268.320050.854250.918051.117651.477552.035652.833253.885655.123256.339257.335257.978758.200238.314438.380938.592538.987439.634440.633142.083443.955145.798347.273648.206248.523425.641625.698725.884126.245126.881527.981329.859632.815635.625537.754739.049139.480912.854612.888313.000013.227513.665014.551216.558321.822326.133229.070630.754531.30220.00000.00000.00000.00000.00000.00000.000011.782018.014421.640023.596224.21900.00007.291412.502515.878717.771218.38120.00004.88118.825511.601213.228713.76350.00003.40756.31728.47209.778910.21550.00002.43184.56386.19057.19957.54050.00001.75593.31594.52655.28825.54750.00001.27582.41723.31163.87914.07330.00000.93021.76552.42372.84342.98740.00000.67951.29091.77422.08332.1895

0.00000.49690.94451.29901.52601.60410.00000.36360.69130.95111.11761.17500.00000.26610.50610.69640.81850.86060.00000.19480.37050.50990.59940.63020.00000.14270.27130.37340.43890.46150.00000.10450.19870.27340.32140.33800.00000.07650.14550.20020.23540.24750.00000.05600.10650.14660.17240.18120.00000.04100.07800.10740.12620.13270.00000.03000.05710.07860.09240.09720.00000.02200.04180.05760.06770.07120.00000.01610.03060.04220.04960.05210.00000.01180.02240.03090.03630.03820.00000.00860.01640.02260.02660.02790.00000.00630.01200.01660.01950.02050.00000.00460.00880.01210.01420.01500.00000.00340.00640.00890.01040.01100.00000.00250.00470.00650.00760.00800.00000.00180.00350.00480.00560.00590.00000.00130.00250.00350.00410.00430.00000.00100.00190.00260.00300.00320.00000.00070.00140.00190.00220.00230.00000.00050.00100.00140.00160.00170.00000.00040.00070.00100.00120.00120.00000.00030.00050.00070.00090.00090.00000.00020.00040.00050.00060.00070.00000.00020.00030.00040.00050.00050.00000.00010.00020.00030.00030.00040.00000.00010.00020.00020.00020.00030.00000.00010.00010.00020.00020.00020.00000.00000.00010.00010.00010.00010.00000.00000.00010.00010.00010.00010.00000.00000.00000.00010.00010.00010.00000.00000.00000.00000.00010.00010.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000用matlab对上述数据作图，由于生成的数据后面几行的值很小，为了提高生成图像的区分度，将后面几行小于0.001的数据全部去掉，另外，将求解区域以外的点的磁位全部设定为0，得到节点磁位与空间坐标的曲面①=①3,y）的图像，如图9所示。图9通过结果可以发现，区域内的磁位主要集中在上方的矩形内，从上方往下看，磁位在一定的区域内近似是线性递减的，在这些区域内磁场强度几乎是定值，在线性区域外的区域内，磁位很快的衰减到0，因此，下方的矩形区域内，磁位是很快衰减的，主要集中在靠近上极板的区域内，在这个区域内，磁场强度也是比较小的。对于某些地方，甚至可以认为磁场强度为0。为了更好的反应磁位的变化曲线，将那些磁位的数值很小的点再次删除掉，得到如图10所示的图形图10附录（一）两点说明1报告会附上编写的源程序一块叫上，有一点说明的是：第一题会交上两个程序，分别对应于不同的方法，报告中会提到这两种方法（其中对于利用对称性和第二类边界条件的方法报告中不会具体的讨论）2对于第一题第一种方法交上的源程序和报告中提供的源程序清单不大一样，因为要分析不同的加速收敛因子和不同的网格的结果，提交的源程序中加速收敛因子和网格数目都是自己输入进去的，对于每输入一次加速收敛因子和网格数目，程序都会自己计算出结果，这样便于不同情况下的比较，只要编译运行一遍程序即可，不必每次都要在源程序上都输入一遍网格数目以及加速收敛因子，然后在编译运行。（二）上机题第一题的另一种方法及其源程序代码这种方法也比较简单，主要是根据求解区域的对称性，并且对称线上的电位的法向导数为0，取一半区域划分网格进行编程，此时得到的微分方程及其边界条件前面列出过，这儿不再重复，程序清单及运行结果如下：用fortran90编写的程序清单如下：programhomeworklimplicitnoneinteger::i,j,m,n,k,real,allocatable::reala,b,crealdxrealecount!变量定义区u（:,:）!加速收敛因子e=1.4dx=1!dx可以认为是步长，用以调整网格的数目m=8/dx+1n=4/dx+1!每边的节点数allocate(u(m,n))doi=1,m!边界条件u(i,1)=0enddodoj=2,nu(1,j)=10u(m,j)=0enddodoi=2,m-1doj=2,nu(i,j)=10.0*(i-1)/(m-1)!边界内的结点赋初值enddoenddok=1count=0!count表示迭代次数do!k是判断迭代是否达到精度要求的变量，if(