《计算电磁学》第三讲ppt课件.ppt
11/16/2022,第三讲 边界条件及有限差分法应用,Dr. Ping DU (杜平),E-mail: ,School of Electronic Science and Applied Physics， Hefei University of Technology (HFUT),11/16/2022,2,边界条件及其处理,积分形式的麦克斯韦方程,微分形式的麦克斯韦方程,(3.1a),(3.1b),(3.1c),(3.1d),(3.2a),(3.2b),(3.2c),(3.2d),其中, H为磁场强度;B为磁感应强度; D为电通量密度；,为电荷密度;,为电流密度,11/16/2022,3,对线性、均匀、各向同性媒质，有 ， 。,其中 、 分别为介质的介电常数和磁导率。,1不同介质分界面上的处理方法 在实际问题中，常遇到所分析的场域存在不同介质。在不同介质分界面上，电通量是连续的，有,其中， 为电位，,-,(3.3),在不同介质分界面处，电位也是连续的。,切向电场强度也是连续的,11/16/2022,4,图3-1 直线形介质分界面处的差分格式,对式(3-3)进行面积分，并利用二维Gauss定理，得,式中，是垂直于区域S围线l的外法线矢量。将S区域各边上的 用其所在边中心点处的两点差分表示，可得（3-4）式左边的积分值。如，对a-a边，沿线的积分为,(3.4),(3.5),11/16/2022,5,对其他三个边类似处理，可得,经过整理，可得,从式(3-7)可以看出，在分界面上的等效相对介电常数为 ，即取平均值。,对于具有角点的介质交界面情形,(3.7),(3.6),11/16/2022,6,图3-2 含角点的介质分界面,角点处的电位 为,(3.8),11/16/2022,7,2. 边界条件的处理,三类边界条件,第一类边界条件：,当网格节点位于边界C上时，则取所在位置的值。,若节点不位于边界C上时，有三种处理办法。,(2)线性插值法,(3.10),(3.9),(1)直接转移法,11/16/2022,8,图 3-3 第一类边界条件的差分网格,(3.11),11/16/2022,9,(3)双向插值法,若 ， ， 代入Poisson方程，则有,第二类边界条件：,(a),(b),图3-4 第二类边界条件的差分网格,(3.13),(3.12),11/16/2022,10,若网格点和边界点重合，,否则，可令 ，再利用不等距差分公式计算。,当 为0时，为齐次边界条件。,第三类边界条件：,当 时，降为第二类边界条件。,(3.15),(3.14),11/16/2022,11,图3-5 第二类和第三类边界条件的差分网格,处理办法：过点O向边界作垂线PQ，与边界交于Q点。令OP、PR、VP的长度分别为ah，bh和ch。对点O有，,点P的值由点V和R的插值得到，,(3.16),(3.17),11/16/2022,12,代入（3-16）,且由于,有,由式(3-15),有,由式(3-19)和(3-20)，得点O的差分格式为,(3.21),(3.18),(3.19),(3.20),11/16/2022,13,有限差分法的应用 (Application of the Finite Differential Method ),差分方程组的建立,分析二维Poisson方程的第一类边界问题为例。设场域D为正方形：,假设x方向和y方向的步长相等 。, 。,以这样的网格离散该区域，如图3-6所示。,图3-6 正方形区域的差分网格,11/16/2022,14,五点差分格式为：,提示：用一个例子加以说明；把 写成一维列向量。在确定了网格格式后，就要据此建立线性方程组。用矩阵符号可写成，,其中，K为系数矩阵， 为未知量,F为已知量,引入x方向的层向量,(3.23),(3.22),11/16/2022,15,一般地，,先确定矩阵K。假设一共有33个内节点(N=4)，并以此为例。,图3-7 含33个内节点的区域,(3.24),11/16/2022,16,对节点(1,1)，其差分格式为,图3-7 含33个内节点的区域,对节点(2,1)，其差分格式为,对节点(3,1)，其差分格式为,(3.27),(3.25),(3.26),11/16/2022,17,对节点(1,2)，其差分格式为,对节点(2,2)，其差分格式为,对节点(3,2)，其差分格式为,对节点(1,3)，其差分格式为,对节点(2,3)，其差分格式为,对节点(3,3)，其差分格式为,(3.28),(3.31),(3.29),(3.30),(3.32),(3.33),11/16/2022,18,我们将所有的内节点 ( , )写成一列向量 ，其为,系数矩阵K为,(3.34),(3.35),11/16/2022,19,F为,差分方程具有如下特征：,系数矩阵K是大型稀疏矩阵；,矩阵K往往是对称正定的，且其前主子式都大于零；但当边界和网格节点不重合时，K的对称性将遭到破坏；（如果具有对称性，则利用这一特性可减少约一半的存储量）,(3.36),11/16/2022,20, K通常不可约，因而方程组不能有其中的一部分单独求解。,差分方程组的求解：,(1) 直接法，如高斯消元法，LU分解。在MATLAB中，可以用phi=KF求得。,(2) 迭代法，Jacobi法，Gauss-Seidel法，SSOR法.,若用Jacobi法，第次的近似值可由第次的近似值得到，其公式为,若用Gauss-Seidel法，第n+1次的迭代中，部分值是第n次迭代得到的；有些是刚更新的，其公式为,(3.37),(3.38),11/16/2022,21,观察这两个公式可以看出，前者需要存储第n、n+1这两次迭代的近似值；后者只需要存储第n+1次的近似值。,另外，由数值分析知道，Jacobi法的收敛速度慢于Gauss-Seidel法。,为了进一步加快迭代速度，我们引入加速因子 。由公式(b)构造，可构造新的迭代公式,这就是所谓的超松弛迭代法(SSOR) .,(3.39),11/16/2022,22,加速因子 满足： 。当 时，公式(3.39)变为(3.38)。当 时，迭代公式会发散。,的值对收敛速度有很大影响。对正方形场域的第一类边值问题，最佳的 为,其中，l为每边的节点数。对于用矩形网格分割的矩形区域，假设每边的节点数分别为l+1，m+1，则,(3.40),(3.41),11/16/2022,23,两个算例(Two numerical examples)：,算例一：二维区域中的电位分布 (The potential distribution in a 2-D domain);,算例二：矩形波导的截止波长(The cutoff wavelength of the rectangular waveguides).,首先分析第一个例子。,例1 一无限长接地金属槽，其侧壁及底面电位均为0，顶的电位为100，如图3-8所示。计算槽内的电位分布。,图3-8 无限长接地金属槽,11/16/2022,24,Analytical solution:,(1),(矩形区域，x、y方向长度分别为a，b；上边界的电位为 ，其他三个边界电位均为0)；,(2),(矩形区域，x、y方向长度分别为a，b；右边界的电位为 ，其他三个边界电位均为0)；,步骤：,(1) 离散场域。我们采用正方形网格离散该区域，每边的节点为l+1=5。,(3.42),(3.43),11/16/2022,25,图3-9 无限长接地金属槽的差分网格,(2) 推出采用SSOR法的差分方程形式。在区域内，由于无自由电荷，则f=0。 加速因子 .,(3) 给出边界条件。该问题中，是第一类边界条件，则直接赋值即可。,(4)选初值。令除边界节点外的其他节点电位为0。,11/16/2022,26,(5)给定迭代收敛的条件。如可令每个节点第n和n+1次近似值的误差绝对值小于 。,(6) 程序流程图。,(7) 编写程序。,(8) 计算结果。,(To be contd),11/16/2022,27,图 3-10 程序流程图,(continued),11/16/2022,28,例2 用有限差分法求解矩形金属波导的截止波长和场分布。,我们假设： 波导壁为完纯导体(PEC); 波导内的介质线性、均匀、各向同性； 波导内无自由电荷和传导电流。 波导工作在匹配状态，具有均匀的截面。波导内部不存在反射波。,由电磁学(电磁场与电磁波，电动力学或工程电磁学等)知道，矩形波导内传播的波可分为TE波和TM波。,对于这两类波，可以归结为求解相应的纵向分量 或 所描述的问题。,以 来标记纵向分量，波导场的分析是定义在波导横截面平面内的二维标量波动方程的定解问题，即,(在波导内部),边界条件：,(1)对TE波( ), ;,(2)对TM波( ), .,(3.45),(3.44),11/16/2022,29,分析TE波模式下的截止波长。假设在x和y方向的步长相等，其差分格式为,当节点位于边界上时，在边界外侧设置一排虚设的网格节点。其边界条件是,。,则差分格式分四种情况 (参考图3-11)：,(1) 位于左边界时，差分格式为,(2) 位于右边界时，差分格式为,(3) 位于上边界时，差分格式为,(3.46),(3.47),(3.48),(3.49),11/16/2022,30,图3-11 矩形金属波导的差分网格,11/16/2022,31, 位于下边界时，差分格式为,拐角：,左下拐角：,右下拐角：,左上拐角：,右上拐角：,(3.50),(3.51),(3.52),(3.53),(3.54),11/16/2022,32,将以上个差分格式应用于网格节点，可得到网格节点上的未知量 的n个差分方程。这样，就构成了矩阵方程,观察该矩阵方程，可知求解截止波长变为求解矩阵K的特征值。截止波长公式为,(3.56),(3.55),* 矩形波导截止波长的解析解,矩形波导中可能存在TE波和TM波。截止波长为,(3.57),( ; ),11/16/2022,33,补充：如何求矩阵的特征值。对于方阵A来说，其特征值就是满足,的所有的值。其中,v为相应的本征向量。,求取矩阵的特征值可以用直接法和迭代法。其中迭代法有Lanczos 法、Davidson 法等。这些方法可以在数值分析或计算方法等书中找到，在此不再讨论。,科学计算软件MATLAB里面提供了求取矩阵本征值(特征值)和本征向量(特征向量)的内置函数eig。,例子 求矩阵A的特征值和特征向量。其中A为,用MATLAB可以得到，其命令为：, A=3 7 5;4 8 3;5 0 9; eig(A) % 求取矩阵A的特征值,11/16/2022,34,输出结果为ans = 15.0000 -1.2749 6.2749,求取矩阵A的特征值和特征向量。, v,d=eig(a) % 求取矩阵A的特征值(d)和特征向量(v)；,v = -0.5774 -0.8802 -0.2653 -0.5774 0.2489 -0.6395 -0.5774 0.4041 0.7215d = 15.0000 0 0 0 -1.2749 0 0 0 6.2749,11/16/2022,35,Homework:,Write the code for the No.1 example and analyze the results.,Write the code for calculating the cutoff wavelength of a rectangular waveguide (1.5m1m).,Analyze the relation between the step and the accuracy.,Analyze the relation between the step and the accuracy.,You can choose one of them. A complete report (Title, author, abstract, key terms, introduction, formulas, results &analysis, and conclusion) needs to be turned in. In addition, the code is also required.,11/16/2022,36,Thank you!,
