技能培训专题 工程电磁场讲义工程电磁场导论实验

用有限差分法解静电场边值问题一、目的1．掌握有限差分法的原理与计算步骤；2．理解并掌握求解差分方程组的超松弛迭代法，分析加速收敛因子的作用；3．学会用有限差分法解简单的二维静电场边值问题，并编制计算程序。二、方法原理有限差分法是数值计算中应用得最早而又相当简单、直观的一种方法。应用有限差分法通常所采取的步骤是：⑴采用一定的网格分割方式离散化场域。⑵进行差分离散化处理。用离散的、只含有限个未知数的差分方程组，来近似代替场域内具有连续变量的偏微分方程以及边界上的边界条件（也包括场域内不同媒质分界面上的衔接条件）。⑶结合选定的代数方程组的解法，编制计算机程序，求解由上面所得对应于待求边值问题的差分方程组，所得解答即为该边值问题的数值解。现在，以静电场边值问题)2()()1(02222sfDyxL中在为例，说明有限差分法的应用。f(s)为边界点s的点函数，二位场域D和边界L示于图5.1-1中。0yx02413DLhh图5.1-1有限差分的网格分割1．离散化场域应用有限差分法时，首先需从网格划分着手决定离散点的分布方式。通常采用完全有规律的方式，这样在每个离散点上可得出相同形式的差分方程，有效地提高解题速度。如图5.1-1所示，现采用分别与x，y轴平行的等距（步距为h）网格线把场域D分割成足够多的正方形网格。各个正方形的顶点（也即网格线的交点）称为网格的结点。这样，对于场域内典型的内结点0，它与周围相邻的结点1、2、3和4构成一个所谓对称的星形。2．差分格式造好网格后，需把上述静电场边值问题中的拉普拉斯方程（1）式离散化。设结点0上的电位值为0。结点1、2、3和4上的电位值相应为1、2、3和4，则基于差分原理的应用，拉普拉斯方程（1）式在结点0处可近似表达为1+2+3+4-41=0（3）这就是规则正方形网格内某点的电位所满足的拉普拉斯方程的差分格式，或差分方程。对于场域内的每一个结点，关系式（3）式都成立，都可以列出一个相同形式的差分方程。但是，对于近邻边界的结点，其边界不一定正好落在正方形网格的结点上，而可能如图5.1-2所示。其中1、2为边界线上的结点，p、q为小于1的正数。仿上所述，可推得对这些近邻边界结点的拉普拉斯方程的差分格式为0)11(11)1()1(04321qpqpqqpp（4）式中：1和2分别是给定边界条件函数f(s)在对应边界点处的值，是已知的。01234hhqhph图5.1-2近邻边界的结点3．边界条件的近似处理为了求解给定的边值问题，还必须对边界条件，以及具体问题中可能存在的分界面上的衔接条件，进行差分离散化处理，以构成相应的差分边值问题。这里，我们只考虑正方形网格分割下的边界条件的近似处理。⑴第一类边界条件如果网格结点正好落在边界L上，因此对应于边界条件（2）式的离散化处理，就是把点函数f(s)的值直接赋予对应的边界结点。如果边界L不通过网格分割时所引进的结点（例如图5.1-2中的1、2结点是边界线L与网格线的交点，并不是网格分割时所引进的网格结点），那末在紧邻边界的结点的差分格式应选用（4）式，这时，把点函数f(s)的值直接赋予边界线L与网格线的交点1和2。⑵第二类边界条件应当指出，从实际电场问题的分析出发，如图5.1-3所示，以电力线为边界的第二类齐次边界条件是常见的一种情况。0Ln（5）这时，可沿着场域边界外侧安置一排虚设的网格结点，显然，对于边界结点0，由于该处0n，故必有1=3，因此相应于边界条件（5）式的差分计算格式为21+2+4-40=0（6）D0124h边界s虚设结点3hE10342AAhh图5.1-3第二类齐次边界的一种情况图5.1-4对称线上结点的差分格式同样，在许多工程问题中，常常能够判定待求电场具有某些对称性质，这样只需要计算某一对称部分的场就能完全决定整个场的分布。为此，还必须导出位于场的对称线上的结点所满足的差分计算格式。以对称线与网格结点相重合为例（见图5.1-4），设'AA线为一对称线，对于位于对称线上的任一结点0，由拉普拉斯方程（因对称性，必有1=3）可得相应的差分计算格式是21+2+4-40=0（7）⑶媒质分界面上的衔接条件在此选取两种情况进行差分离散化的处理。分界面与网格线相重合的情况；设分界面L与网格线相重合，如图5.1-5所示，在两种媒质aε和bε中电位都满足拉普拉斯方程。容易导得，两种媒质分界面上衔接条件在结点0的差分格式为04121204321KKK（8）其中baK。02413abhh20413hhab图5.1-5分界面与网格线相重合图5.1-6分界面L对网格呈对角线形态分界面对于网格呈对角线形态的情况：如图5.1-6所示，分界面L对于网格呈对角线形态，在两种媒质a和b中电位都满足拉普拉斯方程。容易导得，两种媒质分界面上衔接条件在结点0的差分格式为04)(12)(1204321KKK（9）其中baK。总之，类似以上的分析处理方法，可以逐个导得各种类型的边界条件和衔接条件差分离散化的计算格式。限于篇幅，在此不再展开。4．差分方程组的求解在对场域D内各个结点（包括所有场域内点和有关的边界结点）逐一列出对应的差分方程，组成差分方程组后，就可选择一定的代数解法，以算出各离散结点上待求的电位值。注意到差分方程组的系数一般是有规律的，且各个方程都很简单，包含的项数不多（最多不超过5项），因此，对于有限差分法，通常都采用逐次近似的迭代方法求解。在迭代法的应用中，为加速迭代解收敛速度，一般采用的是超松弛迭代法。由于编写计算机程序的需要，每一网格结点的位置由双下标（i，j）予以识别，如图5.1-7所示。对于差分方程（3）式，采用超松弛迭代法（规定迭代的运算顺序是：从左下角开始做起，即i小的先做；对固定的i，j小的先做。），则关于结点0迭代到第（n+1）次时的近似值，应由如下迭代公式算得)(),()1()1,()1(),1()()1,()(),1()(),()1(),(4(4njinjinjinjinjinjinji）（10）ij2(i,j+1)0(i,j)3(i-1,j)1(i+1,j)4(i,j-1)图5.1-7结点的双下标（i，j）标号式中：称为加速收敛因子，其取值范围是1≤2，当≥2时，迭代过程将不收敛。加速收敛因子有一个最佳取值问题，但随具体问题而异。对于第一类边值问题，若一正方形场域由正方形网格分割（每边结点数为m+1），则最佳收敛因子0可按下式计算msin120（11）在更一般的情况下，0只能凭借经验取值。应当指出，为加速迭代解收敛速度，在迭代运算前，恰当地给定各内点的初始值（即所谓第0次近似值）也是一个有效的途径。5．迭代解收敛程度的检验在超松弛迭代法的应用中，还必须涉及迭代解收敛程度的检验问题。对此，通常的处理方法是：迭代一直进行到所有内结点上相邻两次迭代解的近似值满足修正条件Wnjinji),()1(),(（12）时，终止迭代。将式(12)作为检查迭代解收敛程度的依据。其中：W是指定的最大允许误差。6．有限差分法的程序框图启动给定边值填写场域内的初始值叠代次数计数n:=0n:=n+1按超松弛迭代法进行一次迭代，求)1(),(nji所有内点相邻二次迭代值的相对误差是否小于W？打印迭代次数n，待求数值解)j,i(停机是否图5.1-8程序框图三、上机作业设有一个长直接地金属矩形槽，（a=2b），如题5.1-1图所示，其侧壁与底面电位均为零，顶盖电位为100V（相对值），求槽内的电位分布。abxy=0=0=0=100Vo题5.1-1图矩形接地金属槽具体要求：⑴编写一个计算机程序（用你熟悉的程序语言）⑵求相邻两次迭代值的指定的最大允许误差小于10-5的迭代收敛解。⑶以步距40ah的正方形网格离散化场域，然后应用有限差分法求电位的数值解。⑷根据场分布的对称性，试以半场域为计算对象，并以步距40ah将该半场域由正方形网格予以分割，然后用有限差分法求电位的数值解。⑸分别取为n个不同的值和最佳值0，求电位的数值解，以此分析加速收敛因子的作用。从迭代收敛时的迭代次数和最终数值解这两方面总结自已的看法。⑹用计算机描绘等位线分布。⑺取中心点)2,2(baP处电位的精确解（解析解）与数值解进行比较，说明误差范围。程序：#includestdio.h#includemath.h#defineM40#defineN20#defineA1.5doublegi[M+1][N+1];//Defineaarraytostorethevaluesoffi[i][j]doubleg(inti,intj)//Defineafunctiontocaculatefi[i][j](n+1){doubleresult;result=gi[i][j]+(gi[i+1][j]+gi[i][j+1]+gi[i][j-1]+gi[i-1][j]-4*gi[i][j])*(A/4.0);returnresult;}main(){doubleg(int,int);intflag=0,i,j,n=0;for(i=0;i=M;i++)//Evaluatetheoriginalvaluesofthetop.gi[i][N]=100;for(i=0;i=M;i++)//Evaluatetheoriginalvaluesofotherareas.for(j=0;jN;j++)gi[i][j]=0;do//Begintocaculate.{n++;flag=0;//Usedasamarktomarkupiftheprocessisfinished.for(i=1;iM;i++)for(j=1;jN;j++){if(fabs(g(i,j)-gi[i][j])10E-5)//Ifthereexistsonepointthatdoesn'tmeettherequirements,thenflagwillnotequals0.flag++;gi[i][j]=g(i,j);}}while(flag);//Ifflagequals0,thenstop.for(j=0;j=N;j++)//Printtheresultsorderbyrows.{for(i=0;i=M;i++)printf(%f,gi[i][j]);printf(\n);}printf(Printed%dtimes.\n,n);//Printthecyclenumber.}运行结果：