“作为一个CAE设计师,有必要了解计算电磁学吗?”
电磁从业人员按照分工大致可以分为两类:一类从事CEM(计算电磁学),一类从事CAE(计算机辅助设计)。CEM工程师日常的工作内容主要就是“码代码”,目的在于设计“速度更快,精度更高,内存占用更小”的仿真工具,而CAE工程师的日常工作内容主要就是“跑仿真”,目的是利用各种仿真工具,完成诸如天线/滤波器/功分器等元器件的设计。他们的分工与“铸剑师”和“剑客”的关系是一样一样的。
该篇为“缘起收敛性”系列的第二篇文章,也是"电磁CAEer"基础理论部分体系性的一篇文章,旨在为电磁CAE工程师提供一份CEM几种重要算法以及对应的CAE仿真工具的使用说明,并使用最为通俗易懂的语言对这些算法的底层原理进行说明,以期让大家对CEM的底层理论获得最为直观的认知。
1.说在前面
上一篇文章,我们说到,基于Maxwell方程的电磁求解是电磁学研究最为重要的组成,“一张纸,一支笔就可以推演出电磁世界的万千种种”的愿望虽然美好,然而现实却总是那么“冰冷”,解析求解这种“炫技”只适用于自由空间辐射或者少数形状十分规整的散射体。直到“离散”的思想被引入,“计算电磁学的春天”来了……
介绍一下“计算电磁学”的发展沿革以及现状,重点介绍时域有限差分算法(FDTD),有限元算法(FEM),矩量法(MoM)以及以MoM为基础发展而来的多层快速多级子算法(MLFMM)四种主流算法的发展过程以及不同使用环境下的优劣势对比;
使用最为通俗易懂的语言介绍FDTD,FEM,MoM以及MLFMM的核心思想以及底层原理,尤其是针对MoM算法进行详细,将会以你意想不到的丝滑方式将它镌刻至你脑海里。
2.历史与现状
电磁计算算法大致可以按照“时域/频域”和“微分/积分”两个标准,通过两两组合划分为四类。三大算法中,按照第一个标准,FDTD为时域算法,而FEM和MoM同为频域算法;按照第二个标准,FDTD和FEM同为微分方法,而MoM为积分方法。
FDTD算法:采用差分直接离散时域Maxwell方程,电磁场的求解基于时间步的迭代,无需存储全空间的电磁场信息,内存消耗较小,同时采用立方体网格和差分算法,网格形式和算法均十分简单,计算速度也贼快,基于时域算法,特别适合“宽带问题”的求解。但是,简单的立方体方体网格带来的弊端就是模型拟合精度较低,对于含有精细结构的模型,计算精度较低,同时基于“微分方程”,计算区域需要设置截断。FDTD比较适合于不含有较多精细结构的电大尺寸模型的电性能计算以及宽带问题的计算;
FEM算法:采用四面体网格对目标进行离散,拟合精度比FDTD算法更高,计算精度也要明显优于FDTD算法。但是,FEM基于频域/微分算法,需要同时对整个区域内的电磁场信息进行求解和存储,内存消耗大,计算速度慢,计算模型的电尺寸也相对较小。FEM主要适合于微波电路器件,天线等目标“辐射问题”的精确计算;
MoM算法:通过“场-源关系”,将“场”的求解问题转化为“源”求解问题,采用的基函数“格林函数”天然满足辐射条件,无需设置截断,计算精度高,同时矩阵的计算采用直接计算,不存在收敛性的问题,同时由于网格的剖分仅存在于目标体表面或内部,未知量数目大幅降低,矩阵规模小于FDTD和FEM,但是由于“源”之间均存在耦合,因此矩阵为“稠密”矩阵,计算复杂度大,计算速度慢。MoM主要适合于含有精细结构的电小尺寸目标“散射问题”的精确计算;
MLFMM算法:以MoM算法为基础发展而来的加速算法MLFMM,显著提高了MoM的计算速度,大幅降低了内存消耗,处理电大尺寸目标的计算得心应手,但是迭代求解的收敛性较差,对于一些病态矩阵,时常出现计算不收敛的情况,同时在对MoM加速的过程中,采取了一些近似处理,计算略精度低于MoM。MLFMM适合“电大尺寸目标”散射问题的精确求解,对于包含有波导缝隙天线,大规模超材料,磁性媒质等目标的计算问题,则经常会出现计算不收敛的情况。
3.底层原理
如上一篇文章所述,方程中包含了对电场 和磁场 分布函数的时间微分运算以及空间微分运算,以电场 为例,其中 表示 对时间的微分运算, 和 表示 对空间的微分运算,其中 和 最终均可以拆解为 , , 的组合,这四种微分运算可以用 这种通用的表达式进行描述,其可近似为一种简单差分形式:
利用这个近似运算,可以将电场 对时间的微分运算转换成前一时刻的电场 与后一时刻电场 之间的运算关系。同样地电场E对空间微分运算转换成前一位置处的电场 与后一位置处的电场 之间的运算关系 ,因而只要给定了电磁场的初始值 和 (初始条件)以及边界值 和 (边界条件),即可以基于差分关系式,通过不断的循环迭代,求解出任意时刻,任意位置处的电磁场分布 和 和。其优缺点十分明显:
优点:由于采用差分方程以及立方体网格,理论基础和网格划分都相对简单,对于电大尺寸的目标求解速度快,内存消耗小,同时采用基于时域算法,因此适合于宽带问题的求解;
缺点:由于采用矩形单元对目标进行离散,对于复杂形状的目标拟合精度相对较差,从而导致含有精细结构的目标求解精度差。
3.2有限元(FEM)
有限元法基于频域Maxwell方程,其求解的对象是时谐电磁场,即电磁场在时间维度上是周期性分布,循环往复,无始无终,时间变量自然也就失去了意义,电磁场只是空间变量的分布函数:
其采用了拟合效果更好的四面体网格对求解区域体进行剖分。
其中 为第j条棱边的切向分量,为待求的切向分量,而 为相应棱边上对应的基函数,一旦将所有未知量 求解出来,则整个空间中电场分布就完成了求解。这类似于傅里叶级数中使用三角函数展开任意形式的周期函数,所要做的就是求解每个基函数前面系数,然而对于形状不规则的电磁问题,这种基函数的寻找是及其困难甚至不可能的,有限元法的做法是将目标离散成小的单元(三角形,四面体),然后使用非常简单的线性函数或二次函数来近似这个单元上的未知解,这些简单的基函数是一种子域基函数,其与上文中傅里叶级数展开中的全局基函数有着很大的不同。利用有限元将目标离散,并依据电场E在空间Ω满足的波动方程和在边界Γ上满足的边界条件条件建立子域基函数的系数所满足的方程组:
该方程未知量为子域基函数的系数,完成所有未知量的求解,整个空间的电场分布既可以表示为子域函数的叠加。
优点:有限元算法的优点在于其使用了模型拟合度更好的四面体网格,求解精度更高,适合于含有精细结构目标的电磁计算,同时基于Maxwell方程的微分形式,系数矩阵为稀疏矩阵,迭代计算的收敛性较好。
缺点:有限元算法基于“场”的统治方程,有“场”分布的地方就需要进行剖分,然而“场”的分布是充满整个空间的,因此需要设置空气盒子对求解区域进行有限截断,对于电大尺寸目标的计算问题,计算区域十分庞大,剖分的单元和未知量的数目十分庞大,对计算时间和内存的需求也格外高。
3.3矩量法(MoM)
1)格林函数
格林函数为: ,其中 表示场点位置r与源点位置 之间的距离,其为波动方程:
代表空间点源的分布形式,而 就是这个点源在空间中产生的电磁场分布。
2)源-场关系
如上所述,格林函数给定了一个点源在空间中辐射的电磁波的分布形式,同样的我们用光的传播来进行类比,格林函数给出了一个功率为1W的LED小灯泡在空间各个地方产生的光照强度的计算公式,但是如果空间中不只一个LED灯泡,而是有N个LED小灯泡组成的阵列,而且这些每个小灯泡的功率还不都是1W(姑且假设第i个小灯泡位置为 的功率为 ),那么这个LED阵列在空间中任意位置r处产生的光照强度 就可表示为每个小灯泡在r处产生的光照强度的叠加,而每一个小灯泡在r处产生光照强又应该是1w小灯泡的 倍,这应该是显而易见的,即为:
好的,我们现在把视线拨回“电磁场”的求解,通过对Maxwell方程的变换,可以得到电场“E”和磁场“H”所满足的波动方程,方程的形式为:
当你看到这个方程时,不知是否有立马原地爆炸的冲动,但其实这个方程与格林函数满足的波动方程本质上并无不同,方程的左边与
的左边完全相同,这确定了方程的类型以及解的形式,方程的右边只是由原来的点源分布 变成了更为复杂的 ,在电场“E”满足的方程中:
在磁场“H”满足的方程中:
其中“J”和“K”分别为目标体上的电流与磁流,其仅存在于目标体上,在自由空间中没有分布。此处的 (电流源,磁流源的体或面分布)与 (点源)的区别就如同包含无数LED小灯泡组成的灯列阵与一个LED小灯泡的区别,那么”E”和“H“的解,可以类比的表示为:
所不同的是,前面LED灯阵列为离散分布,而此处的电流“J”和磁流“K”的分布为连续分布,因此求和“∑”就变成了求积分“∫”。至此,空间中的场分布“E”和“H”就和目标体上的源分布“J”和“K”建立了联系,这就是“源-场关系”。
3)等效原理
体等效原理
体等效原理就说:以上两种物理模型,对于2区域的电磁场分布是等效的,目标体1对 空间中电磁场分布的影响可以用分布在目标体内的等效源(体电流 和体磁流 )在空间中产生的散射场进行等效。等效源 和 产生的散射场 和 可以依据上文中的源-场关系进行求解:
至于等效源 和 具体是多少,那就是后续需要求解的未知量啦。
面等效原理
面等效原理就说:以上两种物理模型,对于2区域的电磁场分布是等效的,目标体1对 空间中电磁场分布的影响可以用分布在目标体表面的等效源(面电流 和面磁流 )在空间中产生的散射场进行等效。同样地,等效源 和 产生的散射场 和 可以依据源-场关系进行求解:
而等效源和会作为未知量进行求解,这里需要补充说明以下:我们在FEKO中利用MoM或者MLFMM算法进行计算的时候,绝大部分情况下都是利用面三角形网格对目标体的表面进行剖分,而并不对目标体内部进行剖分,其主要原因就是对于均匀介质或者金属,仅存在于目标体表面的面等效源即可等效替代整个目标体的作用,因此考虑目标体内部的情况。
4)矩量法计算模式
Step1:建立支配方程
在理想导体表面,由于电场的切向分量始终为零,所以只有面电流 没有面磁流 ,因此等效电流产生的散射场就可以简化为:
同时空间中的总电场 ,即总电场可以表示为入射场 与散射场 之和,从而有:
其中散射场 可以由源-场关系转化为电流源 的表达式,略加化简,可得:
Step2:离散
同样地,台阶折线的表达式,也是由许许多多的乐高积木块组合而成,这里的乐高积木块叫做“子域基函数”,它的表达式为:
“基函数”就表示它是一个基本单元,“子域”表示它在特定区域内有值,区域外均为0。如图所示的黄色曲线就是 的形式,那么位于x=14处的台阶就可以表示为:
其中 为基函数的系数,也是台阶的高度。整个台阶的分布函数 可以表示为所有台阶之和,即:
整个台阶的分布函数可以表示为所有台阶之和,即:
其中 , , 为棱边 在 上对应顶点坐标, 为棱边 在 上对应顶点坐标, 和 分别为三角形 和 的面积,如此这般,三角形上的每个棱边都可以定义一个矢量形式的基函数 ,那么分布于该三角面元上的电流 可以用定义于三个棱边上的矢量基函数 , 以及 通过线性叠加的方式进行表示:
以上关系对于目标体上的任意三角面元均成立,即处于任意三角面元上的电流分布总是可以表示为公共棱边上的基函数的线性叠加,所以对于目标体表面电流 有:
Step3:匹配
也通过目标体的离散,建立了表面电流 的“假想解”:
其实这个“假想解”还是一个未知解,基函数 前面的系数 统统都是未知量,我们最终的目的就是要求解出所有的 ,从而就完成了表面电流 的求解,这个过程就是“匹配”,顾名思义,用“假想解”去匹配“真实解”。
其中, 依据入射电场求得,为已知量, 为电流分布的系数,是待求的系数, 的具体表达形式较为复杂,我就不展开了,但其物理含义代表了m位置处的电流源 和n处的电流源 之间的耦合,可以通过积分进行求解。通过解方程组,可以求得所有位置处的电流源 ,则整个目标体上的电流分布为:
再利用“源-场关系”可以求得等效源的散射场“ ”和“ ”的分布:
空间中总场分布就是散射场与入射场之和。
3.4快速多级子(MLFMM)
“矩量法”的模式
“快速多级子”的模式
好了,如同”裹脚布“一般的这篇文章终于结束了。首次尝试这种万字长文,从年前写到年后,拖更了快两个月,在经历了一个多月的“学习-消化-理解“这样一个过程后,终于以一种我认为最通俗易懂的方式完成了枯燥深奥的计算电磁学理论的讲解。希望各位读者大人能分三个层次掌握计算电磁学的基本理论。第一层次:对文中有关计算电磁学基本理论的类比图形成深刻的印象;第二层次:能掌握文中关于计算电磁学理论的几个重要的数学推演过程;第三层次:能通过文末补充的专业学术文献和课程笔记全面掌握计算电磁学的所有数学推演过程与结论(当然,这也就成神了)。
最后,我想回到文章最开始的疑问:“作为一个CAE设计师,有必要了解计算电磁学吗?”
看完这篇文章,明天你的仿真该花多长时间可能还是花多少时间 ,该遇到的bug可能一个也不会少,该不收敛的时候可能还是不收敛......,它可能无法解决你仿真过程中遇到的棘手问题,但却可以像武林绝学中的“心法“一样,帮助你在电磁CAE的修罗场里,能够始终多一份笃定与从容......