1.3 电磁场数值计算方法
1.计算电磁学的重要性
在现代科学研究中,科学实验、理论分析、高性能计算已经成为三种重要的研究手段。在电磁学领域中,经典电磁理论只能在11种可分离变量坐标系中求解麦克斯韦方程组或其退化形式,最后得到解析解。解析解的优点如下。
● 可将解表示为已知函数的显式,从而计算出精确的数值结果。
● 可以作为近似解和数值解的检验标准。
● 在解析过程中和在解的显式中,可以观察到问题的内在联系,以及各个参数对数值结果所起的作用。
利用这种方法可以得到问题的准确解,而且效率比较高;但是其适用范围太窄,只能求解具有规则边界的简单问题。当遇到不规则形状或任意形状的边界问题时,需要比较复杂的数学技巧。20世纪60年代以来,随着电子计算机技术的发展,一些电磁场的数值计算方法也迅速发展起来,并在实际工程问题中得到了广泛的应用,从而形成了计算电磁学研究领域,成为现代电磁理论研究的主流。简而言之,计算电磁学是在电磁场与微波技术学科中发展起来的,它建立在电磁场理论的基础上,以高性能计算机技术为工具,运用计算数学方法,专门解决复杂电磁场与微波工程问题。相对于经典电磁理论分析而言,在应用计算电磁学解决电磁学问题时,受边界约束大为减小,可以解决各种类型的复杂问题。从原则上来讲,从直流电到光都属于该学科的研究范围。近几年来,电磁场工程在以电磁能量或信息的传输、转换过程为核心的强电与弱电领域中发挥了重要作用。
2.计算电磁学的分类
(1)时域方法与频域方法。电磁学的数值计算方法可以分为时域(Time Domain,TD)方法和频域(Frequency Domain,FD)方法两大类。
时域方法是指对麦克斯韦方程组按时间步进后求解有关场量。最著名的时域方法是时域有限差分法(Finite Difference Time Domain,FDTD)。这种方法通常适用于求解在外界激励下场的瞬态变化过程。若使用脉冲激励源,则一次求解可以得到一个很宽频带范围内的响应。时域方法具有可靠的精度、更快的计算速度,并且能够真实地反映电磁现象的本质,特别是在诸如短脉冲雷达目标识别、时域测量、宽带无线电通信等研究领域,更是具有不可估量的作用。
频域方法基于时谐微分、积分方程,通过对N个均匀频率采样值进行傅里叶逆变换可得到所需的脉冲响应,即研究在时谐(Time Harmonic)激励条件下,经过无限长时间后的稳态场的分布情况。使用这种方法,每次计算只能求得一个频率点上的响应。在过去,这种方法被大量使用,主要原因是信号、雷达一般工作在窄带上。
当要获取复杂结构时域超宽带响应时,如果采用频域方法,则需要在很大带宽范围内的不同频率点上进行多次计算,然后利用傅里叶变换获得时域响应数据,计算量较大;如果直接采用时域方法,则可以一次性获得时域超宽带响应数据,大大提高计算效率。时域方法还能直接处理非线性媒质和时变媒质问题,具有很大的优越性。时域方法使电磁场的理论与计算从处理稳态问题发展到能够处理瞬态问题,使人们处理电磁现象的范围得到了极大的扩展。
频域方法可以分成基于射线的方法(Ray-based)和基于电流的方法(Current-based)。基于射线的方法包括几何光学法(GO)、几何绕射理论(GTD)和一致性绕射理论(UTD)等;基于电流的方法主要包括矩量法(MoM)和物理光学法(PO)等。基于射线的方法通常用光的传播方式来近似电磁波的行为,考虑射向平面后的反射与经过边缘、尖劈和曲面后的绕射。当然,这些方法都是高频近似方法,主要适用于那些目标表面光滑且其细节对于工作频率而言可以忽略的情况。同时,这些方法对近场的模拟不够精确。基于电流的方法一般通过求解目标在外界激励下的感应电流来求解感应电流产生的散射场,但真实的场为激励场与散射场之和。在基于电流的方法中,最著名的是矩量法。矩量法严格建立在积分方程的基础上,在数字上是精确的。其实,我们并不能判断矢量法是一种低频方法还是一种高频方法,只是该方法所需的存储空间和计算时间随未知元数的快速增长阻止了其在高频情况下的应用,因而被限定在低频至中频的应用中。物理光学法可以认为是矩量法的一种近似,它忽略了各子散射元间的相互耦合作用,这种近似对大而平滑的目标是适用的,但是对于目标上含有边缘、尖劈和拐角等外形的部件,它就失效了。当然,对于形状简单的物体,物理光学法还是一个常用的方法,毕竟其求解过程很迅速,并且所需的存储空间也非常小。
(2)积分方程法与微分方程法。从求解的方程形式上又可以将电磁学的数值计算方法分为积分方程(IE)法和微分方程(DE)法。IE法与DE法相比,其特点如下:①IE法的求解区域维数比DE法的求解区域维数少一维,误差仅限于求解区域的边界,故精度高;②IE法适宜求解无限域问题,而DE法在用于无限域问题的求解时会遇到网格截断问题;③IE法产生的矩阵是满的,阶数小,而DE法产生的矩阵是稀疏的,阶数大;④IE法难以处理非均匀、非线性和时变媒质问题,而DE法则可以直接用于求解这类问题。因此,求解电磁场工程问题的出发点有四种方式:频域积分方程(FDIE)、频域微分方程(FDDE)、时域微分方程(TDDE)和时域积分方程(TDIE)。
计算电磁学也可以分成基于微分方程的方法和基于积分方程的方法。其中,基于微分方程的方法包括时域有限差分法(FDTD)、时域有限体积法(FVTD)、频域有限差分法(FDFD)、有限元法(FEM)。在基于微分方程的方法中,其未知数从理论上讲应定义在整个自由空间内,以满足电磁场在无限远处的辐射条件。但是计算机只有有限的存储量,因此人们引入了吸收边界条件,以此来等效无限远处的辐射条件,使未知数局限于有限空间内。即便如此,其涉及的未知数数目依然庞大(相比于边界积分方程而言)。同时,偏微分方程的局域性使得场在数值网格的传播过程中形成色散误差,而且研究的区域越大,色散误差就越大。数目庞大的未知数和数值耗散问题使得基于微分方程的方法在分析电大尺寸目标时遇到了困难。对于有限元法,早期基于节点(Node-based)的处理方式非常有可能由于插值函数的导数不满足连续性而出现不可预知的伪解问题,使得这种在工程力学中非常成功的方法在电磁学领域无法大展身手,直到一种基于棱边(Edge-based)的处理方式出现后,这个问题才得以解决。
基于积分方程的方法主要包括各类基于边界积分方程与体积积分方程的方法。与基于微分方程的方法不同,其未知数通常定义在源区。例如,对于完全导电体(金属),其未知数仅存在于表面,这显然比基于微分方程的方法少很多未知数。格林函数的引入使得电磁场在无限远处的辐射条件已解析地包含在方程中。场的传播过程可由格林函数来精确描述,因而不存在色散误差的积累效应。
(3)几种主要方法的比较。这里对计算电磁学中几种主要的数值计算方法进行简单的比较,包括时域有限差分法(FDTD)、有限元法(FEM)、矩量法(MoM)、多极子法(MMP)、几何光学绕射法(GTD)、物理光学绕射法(PTD)和传输线法(TLM),如表1-3所示。
表1-3 几种主要的数值计算方法的比较
(4)多种方法的混合使用。由于实际问题的多样性,单独使用以上介绍的方法可能并不能满足需要,如涂敷媒质的目标、印制电路板及微带天线的辐射散射/EMC分析、带复杂腔体和缝隙结构的目标的散射等。因此,工程界常常将各种方法搭配起来使用,形成各种混合方法。常见的混合方法包括边界积分方程与体积积分方程/微分方程方法的混合、高频近似方法与低频精确方法的混合、解析方法与数值方法的混合等。
高频近似方法与低频精确方法的混合方法一般是针对含有复杂细节的电大尺寸目标提出的。由于完全使用低频精确方法处理电大尺寸部分往往会超出目前计算机的能力,而单纯使用高频近似方法又得不到足够精确的近场,所以这种分而治之的折中方案就出现了。常用的混合方法包括弹跳射线法/矩量法混合(SBR/MoM)、物理光学绕射法/矩量法混合(PTD/MoM)、几何光学绕射法/矩量法混合(GTD/MoM)等。虽然引入了高频近似,赢得了速度和空间,但在一定程度上损失了精度。
除了上述几种混合方法,将解析方法和数值方法混合也是一种非常有用的方法。例如,在二维非均匀媒质电磁学问题中,将二维的数值计算转化为径向本征模式展开与纵向解析递推的数值模式匹配法(NMM);对于n维偏微分方程,先使用n-l维数值离散方法将其转化为常微分方程,再用解析方法求其通解的直线法都是很好的例子。
(5)算法的快速求解。快速算法是为了解决矩量法求解过程中存储量和计算量过大的问题而出现的。近年来,许多学者致力于精确方法的快速求解,以满足工程中日益增长的对电大尺寸复杂物体进行精确模拟的需要。由于矩量法产生的是一个满阵,存储量为O(N2),采用直接求解的计算复杂度为O(N3),采用迭代求解的计算复杂度为O(N2)。当未知数数目N增大的时候,存储量和计算量都会快速增加,这极大地限制了其求解能力。而某些基于矩量法的快速算法(如多层快速多极子算法)则可以成功地将存储量和计算复杂度分别降到O(N)和O(NlogN)量级,极大地提升了其求解能力。这些方法主要有基于分组思想的快速多极子方法(FMM)、多层快速多极子算法(MLFMA)、快速非均匀平面波算法(FIPWA)、自适应积分法(AIM)、共轭梯度快速傅里叶变换(CG-FFT)等。
并行计算也称高性能计算,它在现有的算法基础上增加计算资源等硬件设施,把待求解的问题分解为许多小问题,并分别在不同的处理器上求解,通过网络等方式实现进程间的通信,最后得到需要的解,从而实现联合求解。并行计算机从20世纪中期出现以来,出现了多种体系,主要有并行向量机(PVP)、对称多处理机(SMP)、大规模并行处理机(MPP)、集群(Cluster)、分布式共享存储多处理机(DSM)等。
下面对几种主要的计算电磁学数值计算方法进行简单的介绍。
3.有限元法(FEM)
有限元法是在20世纪40年代被提出的,50年代,用于飞机的设计中。后来,这种方法得到发展并广泛地应用于结构分析问题中。目前,作为广泛应用于工程和数学问题的一种通用方法,有限元法已非常著名。
有限元法是以变分原理为基础的一种数值计算方法。它应用变分原理把所要求解的边值问题转化为相应的变分问题,通过对场域进行剖分、插值,将离散化变分问题变为普通多元函数的极值问题,进而得到一组多元的代数方程组,此时只需求解代数方程组就可以得到所求边值问题的数值解。有限元法对微分形式的麦克斯韦方程组在频域进行求解,其求解的未知数是每个小网格的电场与磁场。有限元法一般会对整个求解空间用四面体进行划分,并计算4个格点上的场强,四面体的内部场强分布由4个格点插值得出。通常在待仿真的金属结构和变化比较复杂的部分,网格划分得更密。
有限元法从原理上可以对任意形状的结构进行求解,类似于暴力破解算法,对结构的要求少;但是它消耗的仿真资源多、仿真速度慢。
因为要对整个空间进行网格划分,所以有限元法实际上更适合封闭的空间。但在实际操作中,仿真器通过引入吸收边界条件(Absorbing Boundary Conditions)和完美匹配层(Perfectly Matched Layers)等技术,已经能够极好地解决开放空间的求解问题。因此,有限元法已经不再局限于封闭空间。片上无源空间也是一个开放空间求解问题,我们一般要设置空气盒子的六面为辐射边界,以得到更准确的结果。
有限元法的求解步骤如下。
(1)区域离散化。将场域或物体分为有限个子域,如三角形、四边形、四面体、六面体等。
(2)选择插值函数。选择插值函数的类型(如多项式),用节点(图形定点)的场值求取子域各点的场的近似值。插值函数可以选择为一阶(线性)、二阶(二次)或高阶多项式。尽管高阶多项式的精度高,但通常得到的公式比较复杂。
(3)方程组公式的建立。可以通过里茨方法或伽辽金方法建立方程组公式。
(4)选择合适的代数解法求解代数方程,即可得到待求解边值问题的数值解。
有限元法的特点如下。
● 最终求解的线性代数方程组一般为正定的稀疏系数矩阵。
● 特别适合处理具有复杂几何形状物体和边界的问题。
● 便于处理有多种媒质和非均匀连续媒质的问题。
● 便于计算机实现,可以将其做成标准化的软件包。
4.矩量法(MoM)
矩量法是计算电磁学中最常用的方法之一。自从20世纪60年代Harrington提出矩量法的基本概念以来,它在理论上日臻完善,并广泛应用于工程之中。特别是在电磁辐射与散射及电磁兼容领域,矩量法更显示出其独特的优越性。
矩量法的基本思想是将几何目标剖分离散,并在其上定义合适的基函数,然后建立积分方程,用权函数进行检验,从而产生一个矩阵方程,求解该矩阵方程,即可得到几何目标上的电流分布,从而求得其他近/远场信息。矩量法对积分形式的麦克斯韦方程组在频域求解时,需要求解的未知数为金属的表层电流分布,得到电流分布后,仿真器只需根据格林函数进行数值积分即可得到待求解空间任何点的场分布。在有限元法中,未知数为空间每个点的场分布,其求解矩阵维度大于矩量法的求解矩阵维度。
矩量法的另一大优势是整个无限大的背景结构的信息已经包含在格林函数中,在计算时,它只需对待求解的金属结构划分网格,而不需要对媒质层划分网格,因此,其网格数目小于有限元法的网格数目。因此,对于特定结构(3D层状结构),矩量法的求解速度快、消耗的计算资源少。
对于任意结构或非均匀媒质,矩量法在理论上也可以求解。但是需要使用体积/表面积积分方程对背景环境进行描述,导致未知数数目急剧增加、求解效率急剧下降,反而不如求解微分方程的有限元法高效。
矩量法的求解过程如下。
(1)离散化过程:主要目的是将算子方程转化为代数方程。算子方程为
式中,f为未知等效流或场;g为已知激励源。算子L的定义域适当地选择一组线性无关的基函数(f1,f2,…,fn)(或称展开函数),将未知函数f在算子L的定义域内展开为基函数的线性组合并取有限项近似,即
再将式(1-13)代入算子方程中,利用算子的线性性质,可将算子方程转化为代数方程,即
于是,求解未知函数f的问题就转化为求解未知系数an的问题。
(2)取样检验过程:为了使未知函数f与其近似函数fN之间的误差极小,必须进行取样检验,在抽样点上使加权平均误差为零,从而确定未知系数an。
在算子L的值域内适当地选择一组线性无关的权函数(又称检验函数)Wm,将其与上述代数方程取内积进行抽样检验,即
利用算子的线性和内积性质,将式(1-15)转化为矩阵方程,得
于是,求解代数方程的问题就转化为求解矩阵方程的问题。
(3)矩阵的求逆过程:一旦得到了矩阵方程,通过常规的矩阵求逆或求解线性方程组,就可以得到矩阵方程的解,从而确定展开系数an,得到原算子方程的解。
矩量法的特点如下。
(1)矩量法是基于电磁场积分方程的数值计算方法。积分方程的主要优点在于:一方面,由于格林函数的引入,电磁场在无限远处的辐射条件已经解析地包含在积分方程中,这样可以准确地得到未知数之间的关系,避免数值色散误差;另一方面,它产生的未知数的数目一般比基于微分方程的方法产生的未知数的数目少很多,比较适用于解决电大尺寸的电磁散射问题。
(2)矩量法是一种精确方法,其结果精度仅受计算精度和计算模型精度的限制,因此,它可以实现任意精度下的计算和求解。
(3)矩量法是一种稳定的计算方法,在整个矩量法的求解过程中,不易出现类似于其他计算方法计算过程中出现的伪解问题,同时,它得到的矩阵条件数好,容易求解、求逆。
(4)对于金属表面,矩量法可以利用边界条件直接简化计算,从而导出金属表面的积分方程,而其他计算方法则往往要完全计算整个实体的场分布,这就体现出矩量法在分析金属表面问题时的优越性。
(5)矩量法的全局性使得它产生的矩阵为稠密矩阵,这样,经典矩量法的数据存储量和计算复杂度都很高。因此,快速算法的研究成为矩量法应用研究中的一个热点。
5.时域有限差分算法(FDTD)
从Yee于1966年在解决电磁散射问题的时候提出最初思想到现在,FDTD已经经过了近60年的发展。在此期间,人们不断提出新的思想和方法来克服FDTD的一些缺点。例如,在时间步进算法上,除了传统的Leap-Frog算法,还发展了线性多步时间步进算法(如交错后向差分算法和交错式Adams-Bashforth算法),单步时间步进算法(如Runge-Kutta算法和离散积异卷积法),伪谱算法(如采用Laguerre多项式、交替方向隐式时间步进算法)等;在空间离散上,除了传统的基于Taylor级数展开定理的中心对称有限差分格式,还发展了Discrete Singular Convolution(DSC)格式、非标准的有限差分格式、基于窗函数法的中心对称有限差分格式、最优有限差分格式、FFT等。至此,FDTD已经形成了一个庞大的算法族。
传统电磁场的计算主要是在频域上进行的,但近年来,时域计算方法也越来越受到重视。FDTD是电磁场的一种时域计算方法,它已在很多方面显示出其独特的优越性,在解决有关非均匀媒质、任意形状和复杂结构的散射体及辐射系统的电磁问题中更加突出。FDTD可以直接求解依赖时间变量的麦克斯韦旋度方程,利用二阶精度的中心差分近似把旋度方程中的微分算符直接转换为差分形式,这样实现了在一定体积内和一段时间上对连续电磁场的数据进行取样压缩。电场和磁场分量在空间被交叉放置,这样可以保证在媒质边界切向场分量的连续条件自然得到满足。在笛卡儿坐标系中,每个磁场分量由4个电场分量包围,每个电场分量也由4个磁场分量包围。
这种电磁场的空间放置方法符合法拉第定律和安培定律的自然几何结构。因此,FDTD是计算机在数据存储空间对连续实际电磁波的传播过程在时间进程上进行的数字模拟。在每个网格点上,各场分量的新值均仅依赖该点在同一时间步的值及该点周围邻近点其他场前半个时间步的值,这正是电磁场的感应原理。这些关系构成了FDTD的基本算式,通过逐个时间步对模拟区域各网格点进行计算,在执行到适当的时间步数后,即可获得所需的结果。
FDTD的特点如下。
(1)直接时域计算。FDTD直接把含时间变量的麦克斯韦旋度方程在Yee氏网格空间转换为差分方程。在这种差分格式中,每个网格点上的电场(或磁场)分量仅与它相邻的磁场(或电场)分量及上一时间步该点的场值有关。在每一时间步计算网格空间各点的电场和磁场分量,随着时间步的推进,即能直接模拟电磁波及其与物体的相互作用过程。FDTD把各类问题都作为初值问题来处理,使电磁波的时域特性被直接反映出来。这一特点使它能直接给出非常丰富的电磁场问题的时域信息,给复杂的物理过程描绘出清晰的物理图像。如果需要频域信息,则只需对时域信息进行傅里叶变换即可。如果想获得宽频带信息,则只需在宽频谱的脉冲激励下进行一次计算即可。
(2)广泛的适用性。FDTD的直接出发点是概括电磁场普遍规律的麦克斯韦方程组,这就预示着该方法具有广泛的适用性,近几年的发展完全证实了这一点。从具体的算法来看,在FDTD的差分方程中,被模拟空间电磁性质的参量是按网格空间给出的,因此,只需设定相应的空间点以适应参数,即可模拟各种复杂的电磁结构。媒质的非均匀性、各向异性、色散特性和非线性等能很容易地被精确模拟。由于网格空间中的电场和磁场分量是被交叉放置的,而且在计算中用差分代替了导数,所以媒质边界切向场分量的连续条件能自然得到满足,这就为模拟复杂的电磁结构提供了极大的方便,任何问题只要能正确地对源和结构进行模拟,FDTD就能给出正确的解答。
(3)节省计算机的存储空间和计算时间。很多复杂的电磁场问题都不能计算解决,往往不是因为没有可选用的方法,而是因为计算条件的限制。当代电子计算机的发展方向是运用并行处理技术进一步提高计算速度。并行计算机的发展推动了数值计算中并行处理的研究,适合并行计算的发展将发挥更大的作用。例如,前面指出的FDTD的计算特点是每一网格点上的电场(或磁场)只与其周围相邻点处的磁场(或电场)及其上一时间步的场值有关,这使得它特别适合用于并行计算。施行并行计算可使FDTD所需的存储空间和计算时间减少为只与N1/3成正比。
(4)计算程序的通用性。由于麦克斯韦方程组是FDTD计算任何问题的数学模型,因而它的基本差分方程具有通用性。此外,吸收边界条件和连续条件对很多问题都是通用的,而计算对象的模拟则是通过给网格赋予参数来实现的,与以上各部分没有直接联系,可以独立进行。因此,一个基础的FDTD计算程序对电磁场问题具有通用性,对不同的问题或计算对象只需修改有关部分即可,大部分是相同的。
(5)简单、直观,容易掌握。FDTD直接从麦克斯韦方程组出发,不需要任何导出方程,这样就避免了使用更多的数学工具,使得它成为所有电磁场数值计算方法中最简单的一种。由于FDTD能直接在时域中模拟电磁波的传播及其与物体作用的物理过程,所以该方法又是非常直观的一种方法。由于它既简单又直观,所以掌握它不是件很困难的事情,只要掌握电磁场的基本理论知识(不需要很多数学知识),就可以学习运用这一方法解决很复杂的电磁场问题。因此,这一方法很容易得到推广,并在很广泛的领域发挥作用。
6.弹跳射线法(SBR)
弹跳射线法(SBR)技术是由Hao Ling等于20世纪80年代末提出的高频算法,这是一种结合了几何光学法和物理光学法优点的高频近似方法:首先,将入射的平面电磁波用一定密度的密集射线管来模拟入射波在目标几何结构中的传播情况,根据几何光学法的原理,相邻射线管之间没有能量的耦合和泄漏,当射线照射目标时,按照光学原理发生反射来追踪所有射线管在空间中的行进轨迹;然后,利用几何光学法追踪所有射线管中的场值变化,并在射线管离开目标表面时,利用物理光学法的原理对射线管的远场散射场进行计算;最后,累加所有射线管的远场散射场贡献即可得到目标的总散射场。
7.如何选择合适的电磁场仿真算法
(1)结构的特点。对于层状结构,矩量法能够提供最高效、快速地求解,因此可以优先考虑矩量法。例如,PCB走线、层状倒封装、片上无源器件,都可认为是层状结构。对于转换头、接口、波导、三维天线、BGA等复杂的非层状结构,只能选择有限元法和FDTD。
当然,对于极其简单的结构,如单圈电感等,采用有限元法也可以很快得到结果,速度差别并不明显。
(2)结构的响应类型。由于有限元法和矩量法均从频域求解,所以它们比FDTD更适合那些具备窄带响应或高Q值(品质因数)的结构,如滤波器、谐振腔、波导等。FDTD从时域求解,因而天生更适合分析TDR和EMI、EMC等。对于宽带响应,FDTD也更高效,它在时域求解之后采用傅里叶变换即可得到频率响应;但有限元法和矩量法需要对频带内的频点逐一进行求解(自适应算法可减少求解频点数目)。
(3)结构的复杂程度与端口的数目。对于复杂的3D结构,有限元法的求解效率要高于FDTD的求解效率。另外,有限元法和矩量法的求解时间与端口数目的关系不大,一次求解即可得到所有端口的响应;而FDTD对N个端口需要重复求解N次。因此,像BGA这样的多端口结构,有限元法是更好的选择。表1-4和表1-5分别给出了不同特点的应用及不同的应用领域选择数值计算方法的优/劣势。
表1-4 不同特点的应用选择数值计算方法的优/劣势
注:表中的“+”表示优势,“-”表示劣势。“+”越多表示越适用于该种情况。
表1-5 不同的应用领域选择数值计算方法的方法优/劣势