ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

矩量法金属介质复合结构线性基函数与边元基实现

矩量法金属介质复合结构线性基函数与边元基实现

把金属介质复合结构的矩量法推到代码层面,第一道真正让人卡住好几天的坎,往往不是电场积分方程和磁场积分方程的推导,而是线性基函数怎么定义、怎么选、符号怎么摆。前两步把表面等效电流 J 和等效磁流 M 的积分方程理顺之后,剩下的问题很朴素:面电流和面磁流到底拿什么东西展开?这篇就专门聊这件事,把工程里常见的几种线性基函数——分片常数基、节点线性基、边元基、分段正弦基——逐个拆开,讲清定义、散度为什么是命门、法向连续性为什么省不得,以及在金属面和介质面混排时,怎么用同一套基函数把 J 和 M 一起搭起来。适合已经能跑通纯金属目标矩量法、正准备往介质加载、涂覆、复合天线罩这类结构上扩展的同行,也适合带学生做课题、需要一个能直接抄的基函数实现的人。

1. 先把问题摆正:基函数究竟在为谁服务

1.1 金属面和介质面上的未知量根本不是一回事

很多人做纯金属目标做久了,会形成一种惯性思维:基函数就是拿来展开“电流”的。到了复合结构里,这个惯性会立刻把你带进沟里。金属面上的边界条件是切向电场为零,未知量只有一个切向面电流 J;而介质(或者有耗涂覆层)的分界面上,切向电场和切向磁场都不为零,未知量变成切向电流 J 和切向磁流 M 两套。这意味着同一个几何面片,在金属区贡献一个自由度,在介质区贡献两个自由度,网格在两类区域的交界线上汇合时,自由度编号必须一次规划清楚,否则后面组装矩阵会莫名其妙地漏掉几行几列。

更麻烦的是这两类未知量的物理性质并不对称。J 是真实的面电流密度,单位是安培每米;M 是等效磁流密度,量纲是伏特每米,它在 Maxwell 方程组里的位置和 J 是对偶的,但在实际实现中它并不对应任何真实的电荷运动,因此对它的连续性要求比对 J 宽松得多。J 的法向分量在跨越内边时必须连续,否则就等价于在面上凭空放了一根线电荷,那会直接污染近场和输入阻抗;M 则没有对应的“磁荷守恒”约束那么严格,但也不能随便乱来,因为它参与的是磁流连续性方程,散度处理不当同样会引入伪源。

我在最早实现复合结构求解器的时候,就是因为把金属面和介质面用了同一个数组、同一套编号,结果在金属和介质的接缝处电流基函数被重复计入,金属面上一小块区域的电流密度算出来比周围高出一倍多,查了两天才发现是那条接缝上的边被两个区域各建了一次。这件事之后我给自己立了个规矩:网格可以通用,自由度必须分区管理,金属区只发电流基的编号,介质区同时发电流基和磁流基的编号,接缝上的边归属哪一侧必须在数据结构里写死。

1.2 “随便什么基函数都能收敛”是个危险的错觉

教科书上讲矩量法,往往会花很大篇幅论证基函数和权函数的选取,然后告诉你只要基函数完备、网格足够细,结果总会收敛。这句话在数学上没毛病,在工程上却几乎等于没说。收敛速度、条件数、低频稳定性、近场精度,这四样东西全都和基函数的形状强相关,而它们才是决定你能不能把代码跑出可用结果的关键。

最典型的例子就是分片常数基。用每个三角形上一个常向量去展开电流,理论上网格无限细也能收敛,但实际中你会发现阻抗矩阵的条件数随网格加密迅速恶化,而且面上会出现大量非物理的电荷堆积。原因很简单:相邻两个三角形的电流法向分量在公共边上直接跳变,这个跳变在数学上就是一个面电荷的线密度,物理上根本不存在,但方程里认了。它带来的后果是近场计算严重失真,尤其是贴近表面的区域,电场会抖得很难看。

所以基函数的选择本质上是在做一笔交易:你希望基函数具备哪些“物理上合理”的性质,愿意为此付出多少实现复杂度和计算量。线性基函数这一档里,值得拿出来比较的就那么几种,选型逻辑我后面会展开,但核心判断标准只有一条——这个基函数能不能让电流的散度在面上有一个受控的、可解析的表达。因为散度直接对应电荷,电荷直接对应近场精度和低频行为,这一条过不了关,其他都是白搭。

1.3 本篇讨论的边界在哪里

需要先说清楚范围,免得读者按图索骥却找不到想要的东西。这里讨论的全部是平面三角形面元上的一阶(线性)基函数,几何是直的、不弯曲,网格是三角形剖分。不涉及高阶基函数、不涉及曲面面元上的参数化基函数、不涉及多层快速算法里那种基于聚合树的基函数分组策略,也不涉及体基函数。之所以限定在三角形一阶这一档,是因为它是工程上性价比最高的一段:曲面拟合能力够用,未知量数量可控,实现复杂度在一个人的能力范围内,而且几乎所有的加速算法都是在这个基础上做的二次开发。

还有一个前提要交代:下面涉及的所有符号约定,都建立在“网格是一致定向的”这个基础上,也就是每个三角形的外法向可以通过顶点环序统一确定,相邻三角形在公共边上的环序必然相反。如果你的网格是从外部导入的,一定要先做一次定向一致性检查,否则后面所有的法向、符号、散度全是错的。这个检查本身不难,从一个种子三角形出发广度优先遍历,遇到相邻三角形就根据公共边的方向判断需不需要翻转,遍历完再抽查几个点的法向是否朝外。我吃过这个亏,一个从别的渠道拿来的网格文件,局部有几个三角形的环序是反的,肉眼完全看不出来,但算出来的散射方向图在那个区域有个诡异的凹陷。

2. 四种典型线性基函数的选型逻辑

2.1 分片常数基:最朴素,但电荷表达上先天吃亏

分片常数基的思路极其直白:把每个三角形上的电流当成一个常向量,一个三角形贡献两个未知量(面内的两个切向分量,或者等价的 x、y 两个分量)。它的优点是没有之一——实现简单,编号简单,连求积都省了,因为基函数是常数,积分核和基函数的乘积可以直接提出来。做一维问题、做教学演示、或者写个对照程序验证主程序的正确性,它非常好用。

但它的短板也同样明显。首先是未知量数量翻倍:一个三角形两个未知量,总未知量大约是三角形数的两倍,而三角形数大约是内边数的两倍,所以分片常数基的未知量规模差不多是边元基的四倍。其次是法向不连续,这一点我在上一节已经说过后果了。第三是它无法自然表达电流沿表面的流动趋势,因为相邻三角形的电流值之间没有任何约束,求解器可以自由地让它们在公共边上打架,这种自由度在数学上是冗余的,在数值上就表现为矩阵条件数变差。

我至今保留着分片常数基的代码,但只在一个场景下用:新写一个积分核,先用它跑个最粗糙的结果,确认量级对不对、奇偶性对不对,然后再换正经基函数。它的价值在于“快”,不在于“准”。

2.2 节点线性基(屋顶基):连续性换来了什么,又付出了什么

节点线性基是把面内的电流分量表示成顶点上的分片线性函数的线性组合。想象一个三角形,三个顶点各有一个标量权重,三角形内部的值由重心坐标插值得到,这就是标量屋顶基;电流有两个面内分量,所以每个顶点需要两个未知量。它的支撑域是共顶点的所有三角形,看起来比边元基更“平滑”。

它确实带来了一个好处:在三角形内部,电流是线性变化的,比常数基更能拟合真实分布,尤其是在电流变化剧烈的边缘附近。但它的问题在于跨边的连续性。两个相邻三角形共享一条边,节点基在公共边上的值由两端的顶点值决定,因此它在边上的切向分量是连续的;但法向分量呢?节点基的向量场在跨过公共边时,法向分量一般不连续,除非你特意构造。这个法向跳变又变成了伪线电荷。

还有一个容易被忽略的坑:节点基在结构边缘(开放边界)上的处理非常别扭。顶点上的未知量在边界处没有“对面的三角形”来配对,需要额外构想边界条件,通常的做法是强行置零或者用半基函数,两种做法都会在边缘附近引入误差。金属板的边缘电流本来就是奇异的,这个误差会被放大。

所以节点线性基在实际的金属介质复合结构求解器里用得不多,它更多出现在流体、热传导那类标量场的有限元里。电磁这边,大家更青睐边元。

2.3 边元基:散度友好,法向严格连续,工程上的默认选项

边元基(业内习惯用三位作者名字的首字母缩写来称呼,这里就用这个通用缩写)是我认为在三角形网格上做面电流展开的最优解,没有之一。它的构造思想可以这么理解:把一个未知量绑定在网格的一条内边上,这个未知量代表“跨过这条边的总电流”,基函数在两个相邻三角形上分别是一个线性向量场,方向沿着从自由顶点指向该边的方向,大小按面积归一化。

这个构造的妙处在于三件事同时成立。第一,法向分量在公共边上严格连续,用统一的法向去度量,两侧的值都等于 1,这意味着基函数本身不带任何净线电荷。第二,散度在两个三角形上分别是常数,一个正一个负,面积加权后总和为零,所以每个基函数的净电荷严格为零,电荷只以面密度的形式分布在两个三角形内部。第三,未知量的数量等于内边数,对于一个大致的三角剖分,内边数约为顶点数的 1.5 倍,比节点基少,比分片常数基少得多。

正是这三点,让边元基成了几乎所有工程级矩量法代码的默认选择。它的散度是分段常数这一点,在后面积分计算那一节会显示巨大的价值:所有涉及散度的项都能从双重面积分里提出来,只剩下一个纯标量的格林函数积分,而这个积分是可以预先算好、四个矩阵块共用的。

它当然也有局限。边元基天然只适合描述切向电流,对于开放结构(比如一块孤立的薄金属板)的边界边,它没有定义,必须特殊处理;对于细线结构,三角形剖分本身就不合适,得用线段上的基函数。另外,边元基在低频时会出现所谓的“低频崩溃”,这不是基函数本身的错,而是电场积分方程在低频下的算子性质问题,但基函数的选择会影响崩溃的严重程度。

2.4 分段正弦基:窄带里的老牌选手

分段正弦基是在每条边上定义一段正弦分布的电流,跨接两个三角形,相位和幅度由边两端的值确定。它在细线天线时代是绝对的主流,因为它非常贴近真实电流在细导线上的驻波分布,在谐振频率附近的精度极高,往往用很少的分段就能得到很好的结果。

但它的缺点在宽带和面结构上暴露得很彻底。正弦函数的积分核比线性函数复杂得多,矩阵元素的计算量明显上升;而且正弦基的“波长适应性”是双刃剑,一旦频率偏离设计点,基函数的形状和实际电流分布就对不上了,精度反而下降得比线性基更快。在三角形面元上,分段正弦基的定义也没有边元基那么自然,需要额外约定沿边的参数化方式。

我的建议很直接:如果你的问题是细线天线、窄带、而且对精度要求极高,分段正弦基值得考虑;如果是面结构、宽带、或者复合结构,老老实实用边元基,别折腾。

2.5 四者的横向对比

把上面这些整理成一张表,选型的时候可以直接对照。需要说明的是,“未知量数量”这一列是相对量级,具体数值取决于网格的疏密和形状。

基函数类型支撑域未知量数量级法向连续性散度表达主要优势典型陷阱
分片常数基单个三角形约 2 倍三角形数不连续三角形内为零,边上为奇异线源实现极简,适合做对照验证未知量爆炸,伪电荷严重,近场失真
节点线性基共顶点三角形集合约 2 倍顶点数切向连续,法向一般不连续分片常数,非零净电荷面内线性变化,拟合能力强边界顶点处理别扭,边缘误差放大
边元基一对相邻三角形约等于内边数严格连续(统一法向度量)分片常数,净电荷为零未知量少,电荷受控,散度可提出积分边界边无定义,需特殊处理
分段正弦基一对相邻三角形约等于内边数连续正弦型,需数值积分窄带精度高,贴合驻波分布宽带性能差,核函数复杂

表格只能给出定性的判断,真正的选型还要结合你的频率范围、结构类型和精度要求。我个人的经验是:只要结构里有面,就用边元基;只有在纯细线结构、且带宽很窄的情况下,才考虑分段正弦基。分片常数基和节点基,留在工具箱里当调试工具就好。

3. 边元基的完整定义与实现细节

3.1 定义式与四个几何量

对网格中的第 n 条内边,它有长度 l_n,两侧各挂一个三角形,记作 T_n^+ 和 T_n^-,面积分别是 A_n^+ 和 A_n^-。每个三角形都有一个“自由顶点”,也就是不在公共边上的那个顶点,记作 r_n^{+free} 和 r_n^{-free}。基函数定义为:

在 T_n^+ 内,f_n(r) = (l_n / (2A_n^+)) · (r − r_n^{+free}); 在 T_n^- 内,f_n(r) = −(l_n / (2A_n^-)) · (r − r_n^{-free}); 在其他区域,f_n(r) = 0。

这个式子里藏着一个很容易被忽略的几何事实:在三角形内,从自由顶点指向任意一点的向量 r − r_free,它在垂直于公共边方向上的投影,最大值恰好等于该边上的高 h。而三角形面积 A = l·h/2,所以系数 l/(2A) 乘以 h 正好等于 1。这就是为什么基函数在公共边上的法向分量恒等于 1——系数不是随便凑的,是被面积归一化“设计”出来的。

四个几何量里,边长和面积都好算,自由顶点也好找,真正需要小心的是“哪个三角形是 T+、哪个是 T−”。这个选择不影响最终结果的正确性,但影响所有基函数的符号,一旦在不同的模块里用了不同的约定,矩阵就会不对称。后面 3.3 节专门讲这个。

3.2 散度:为什么它能决定电荷,进而决定精度

对基函数取面散度,结果干净得让人感动:在 T_n^+ 内,∇·f_n = l_n / A_n^+;在 T_n^- 内,∇·f_n = −l_n / A_n^-;其他地方为零。这个结果来自一个恒等式:面内向量场 (r − r_0) 的面散度恒等于 2,因为它在两个面内方向上的偏导各自贡献 1。

把这个散度在两个三角形上分别积分,得到 l_n 和 −l_n,加起来正好是零。这就是“每个基函数净电荷为零”的数学表达。要注意,净电荷为零不等于没有电荷——电荷是以面密度的形式存在的,正电荷集中在一个三角形里,负电荷集中在另一个三角形里,总量抵消。面电荷密度由连续性方程给出,等于散度除以 jω(带符号),这个量在计算近场、计算电荷分布、做低频修正的时候都要用到。

散度是分段常数这件事,价值远不止“算起来方便”。它意味着在组装阻抗矩阵的时候,所有包含散度的项——比如电场积分方程里那个 −1/(jωε) 乘以散度再乘格林函数的项——可以把散度当作常数提到面积分外面,剩下的只是一个标量的双重面积分。而这个标量积分,对于 J-J、J-M、M-J、M-M 四个矩阵块来说是完全一样的,可以算一次存起来反复用。这一招能让组装时间直接降一个量级,是我认为边元基最被低估的一个优点,很多讲基函数的资料只提到法向连续,压根不提这一层。

还有一个细节:散度在公共边上是不连续的,从 l/A+ 跳到 −l/A−。这个跳变看起来像是问题,其实正是它保证了电荷守恒。如果散度在边上连续,那么两侧的散度同号,积分下来净电荷就不为零了,那才是真正的错误。

3.3 符号约定:这一节值得你多看两遍

如果只能给一条忠告,那就是:把 T+ 和 T− 的判定规则写成代码,别靠脑子记。

我见过的符号错误大概有这么几种表现:阻抗矩阵不再对称(对无耗结构,理论上应该对称),而且不对称只出现在某些行某些列;散射参数的相位整体翻转;远场方向图出现镜像对称的错误;最隐蔽的一种是结果看起来“差不多对”,但和商业软件比始终差几个百分点,怎么调网格都不收敛。

避免这些的唯一办法是建立一套完全自动化的判定逻辑,并且在整个代码里只认这一个逻辑。我用的规则是这样的:对每条内边,用它的两个端点构造一个有向对 (v0, v1),规定 v0 的全局编号小于 v1。然后遍历这两个三角形,检查这条边在三角形顶点环序中的出现方向。在环序中方向恰好是 v0→v1 的那个三角形,定义为 T+;另一个定义为 T−。这套规则的前提是网格一致定向,也就是所有三角形的外法向可以通过环序统一确定。

这套规则的好处是它只依赖拓扑,不依赖坐标,也不会因为网格旋转、平移而改变。而且它自动继承了“基函数在公共边上的法向分量,用从 T+ 指向 T− 的单位法向去度量,两侧都等于 +1”这个性质。这里要特别强调“用统一的法向去度量”这句话,因为很多人第一次验证连续性的时候,用的是各自三角形的外法向,结果发现一侧是 +1,另一侧是 −1,就以为自己算错了。其实没错,两个外法向互为反向,用它们各自度量当然会差一个负号,这恰恰说明电流向量本身是同一个,跨边连续。

我的建议是把这条性质做成单元测试:随便挑一条内边,在边上取中点,分别用 T+ 和 T− 的基函数表达式算出向量,这两个向量必须完全相等(浮点误差范围内),并且和统一法向的点积必须等于 1。这个测试跑通了,后面八十 Prozent 的符号问题都不会发生。

3.4 手算一个最小算例:两个三角形拼成的正方形

光看公式容易飘,拿一个能手算的例子走一遍,很多疑问会自己消失。取一个单位正方形,顶点依次为 p0=(0,0,0)、p1=(1,0,0)、p2=(1,1,0)、p3=(0,1,0),沿对角线 p0-p2 剖分成两个三角形,T_a=(p0,p1,p2),T_b=(p0,p2,p3)。两个三角形都按逆时针排列,外法向都是 +z。

公共边就是 p0-p2,长度 l = √2 ≈ 1.41421。两个三角形的面积都是 0.5,这是显然的,对角线把正方形平分成两半。

现在按 3.3 节的规则判定:边的有向对是 (p0, p2),索引 0 小于 2,所以 v0=p0、v1=p2。在 T_a 的环序 p0→p1→p2→p0 中,边 p0→p2 并不作为一个有向边出现,出现的是 p2→p0,所以 T_a 是反向的那一侧。在 T_b 的环序 p0→p2→p3→p0 中,p0→p2 恰好出现,方向与有向对一致。因此 T+ = T_b,T− = T_a。T+ 的自由顶点是 p3,T− 的自由顶点是 p1。

写出基函数:在 T+ 内 f = (1.41421/(2×0.5))·(r − p3) = 1.41421·(r − p3);在 T− 内 f = −1.41421·(r − p1)。

取公共边中点 m = (0.5, 0.5, 0)。T+ 中,f = 1.41421·((0.5,0.5,0) − (0,1,0)) = 1.41421·(0.5, −0.5, 0) = (0.70711, −0.70711, 0)。T− 中,f = −1.41421·((0.5,0.5,0) − (1,0,0)) = −1.41421·(−0.5, 0.5, 0) = (0.70711, −0.70711, 0)。两个向量完全相等,跨边连续,验证通过。

再算统一法向。从 T+ 指向 T− 的方向,也就是从 p3 一侧指向 p1 一侧,单位向量是 (0.70711, −0.70711, 0)。用这个法向点乘 f,结果是 0.5 + 0.5 = 1,两侧都是 1。

散度:T+ 内是 l/A+ = 1.41421/0.5 = 2.82843;T− 内是 −2.82843。面积分验证:2.82843 × 0.5 = 1.41421 = l,负侧同理,和为零。

这个例子虽然简单,但把定义、符号、连续性、散度、电荷守恒全串了一遍。我建议每个刚接手代码的人都手算一遍,然后拿代码验证,对上了再往下走。

3.5 纯 Python 实现:边表构建与基函数求值

数据结构这一层有两个常见选择:一是直接遍历所有三角形的三条边,用哈希表按无向边归并;二是先在网格生成阶段就输出边表。前者对导入的网格更友好,我一般用前者。

import math def sub(a, b): return (a[0]-b[0], a[1]-b[1], a[2]-b[2]) def cross(a, b): return (a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]) def dot(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] def norm(a): return math.sqrt(dot(a, a)) def tri_area(nodes, tri): a, b, c = nodes[tri[0]], nodes[tri[1]], nodes[tri[2]] return 0.5 * norm(cross(sub(b, a), sub(c, a))) def build_edges(tris): """返回内边表:(v0, v1, t_plus, t_minus),v0 < v1""" table = {} for ti, (a, b, c) in enumerate(tris): for u, v in ((a, b), (b, c), (c, a)): key = (u, v) if u < v else (v, u) table.setdefault(key, []).append((ti, u, v)) edges = [] for key, owners in table.items():
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进