把金属介质复合结构的矩量法推到代码层面第一道真正让人卡住好几天的坎往往不是电场积分方程和磁场积分方程的推导而是线性基函数怎么定义、怎么选、符号怎么摆。前两步把表面等效电流 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所以 v0p0、v1p2。在 T_a 的环序 p0→p1→p2→p0 中边 p0→p2 并不作为一个有向边出现出现的是 p2→p0所以 T_a 是反向的那一侧。在 T_b 的环序 p0→p2→p3→p0 中p0→p2 恰好出现方向与有向对一致。因此 T T_bT− T_a。T 的自由顶点是 p3T− 的自由顶点是 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.82843T− 内是 −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():