在电磁仿真这条路上走得久了你会发现一个很有意思的现象同一套积分方程、同一份网格、同一个求解器换一组基函数结果能从和实测对得上变成差出十几dB。金属介质复合结构的矩量法尤其如此——金属表面只有切向电流介质表面却同时存在等效电流和等效磁流界面两侧的场还要靠边界条件缝合。线性基函数在这个场景里扮演的角色远不是教科书上那句把未知函数展开成一组已知函数的线性组合能概括的。它真正决定的是离散之后电荷账本还能不能平跨介质棱上的法向分量还能不能连续矩阵的条件数会不会在低频段直接爆掉。这篇内容面向已经写过或者正在写矩量法程序的人尤其是卡在能跑出数但结果不对这个阶段的同行。我会从基函数承担的物理职责讲起把脉冲、三角、分段正弦这几类标量基函数的取舍逻辑拆开再进入面元上的矢量线性基函数最后落在金属-介质复合界面上最容易被忽略的拼接规则和几个实打实的坑。1. 基函数不是插值工具它在矩量法里真正管的是电荷账本很多人第一次接触矩量法会把基函数理解成数值分析里的插值基——给定若干采样点用一组基函数拟合未知曲线。这个理解在纯数学层面没错但放到电磁积分方程里就偏了。矩量法求解的是算子方程未知量是表面电流密度或者等效电流密度基函数的作用是把连续分布的电流近似成有限个电流元的叠加。电流元的选取方式直接决定了电荷密度怎么分布因为电荷和电流之间隔着一条连续性方程。1.1 离散链条上基函数所处的位置完整的离散链条大概是这样先由等效原理建立积分方程把未知量写成一组基函数的线性组合再用测试函数做内积把算子方程转成矩阵方程。基函数出现在第一步测试函数出现在第二步两者共同决定了阻抗矩阵的每一个元素。写成式子就是未知电流展开加上加权残差J(r) ≈ Σ a_n f_n(r) w_m, L(f_n) w_m, g这里的f_n是基函数w_m是测试函数L是积分微分算子。基函数一旦选定L(f_n)就固定了后面的所有数值误差都从这里长出来。选脉冲基函数L(f_n)在段边界是不连续的选三角基函数L(f_n)在两个相邻段上是分段线性的。这个差别在做自阻抗项积分的时候会放大成完全不同的奇异积分处理策略。我见过不少人调程序先怀疑网格不够密再怀疑积分点不够多最后才回头检查基函数。实际上网格密度和积分精度都是收敛速度问题基函数选错是收敛到错误答案问题两者性质完全不同。前者加密能救后者加密只会让错误答案更精致。1.2 支撑紧致性与连续性两项绕不开的硬指标挑选基函数绕不开两个硬指标。第一个是支撑紧致性也就是单个基函数的非零区域要尽量小。支撑越小阻抗矩阵越稀疏填充和求解的代价越低。脉冲基函数的支撑就是一个段三角基函数是两个相邻段RWG基函数是两个共享棱的三角形都属于紧支撑。第二个是连续性电流在相邻基函数的交界处必须满足物理上要求的连续条件。金属表面电流的切向分量连续法向分量为零介质表面等效电流的切向分量连续法向分量允许跳变跳变量正好等于表面电荷密度。这两个指标经常打架。支撑最紧致的脉冲基函数恰恰是不连续的连续性最好的分段正弦基函数支撑又偏大。工程上的选择本质是在矩阵规模和电荷守恒之间找平衡点。对于电小尺寸结构矩阵规模不是瓶颈优先保连续性对于电大尺寸结构稀疏性带来的收益更明显可以用不连续基函数配合额外的电荷修正项。1.3 金属介质复合结构对基函数的额外要求纯金属结构的矩量法只需要处理一种未知量——表面电流。金属介质复合结构不一样介质区域要用等效原理替换表面同时出现等效电流J和等效磁流M。这意味着基函数要同时服务于两类未知量而且在金属和介质的交界棱上两类未知量的边界条件必须匹配。具体来说金属表面的电流只有切向分量介质表面的等效电流在切向连续、法向跳变。当一条棱同时属于金属三角形和介质三角形时介质侧的电流法向分量必须为零否则就违背了金属表面的边界条件。这个约束在基函数层面就要处理掉不能留到矩阵方程阶段靠罚函数或者后处理补救。我后面会专门用一节讲这个拼接规则因为这是复合结构程序里最容易出错的地方之一。2. 标量线性基函数的三条路线脉冲、三角与分段正弦先从一维问题切入。细线天线、线栅结构这类问题未知量沿线的电流分布标量基函数就够了。三条常见路线分别是脉冲基函数、三角基函数和分段正弦基函数它们的历史都很长但至今仍在大量程序里出现原因各不相同。2.1 脉冲基函数最省事的方案与它埋下的电荷累积隐患脉冲基函数把细线分成若干段每一段上的电流近似为常数。数学上简单到不能再简单第 n 个基函数在第 n 段上取 1其余段取 0。矩阵元素的积分也容易算因为基函数是常数可以提到积分号外面。def pulse_basis(seg_index, seg_length, n_segments): # 第 seg_index 段上取 1其余为 0 coeff [0.0] * n_segments coeff[seg_index] 1.0 return coeff问题出在段与段之间的边界上。电流在边界处发生跳变而电荷密度和电流的导数直接相关。脉冲基函数在段内导数为零在边界处导数变成冲激函数等效于在每一个段边界上都堆了一坨点电荷。这坨电荷是数值离散人为造出来的物理上并不存在。它会污染近场计算让输入阻抗出现虚部偏差在阵列结构里还会让互耦系数算不准。脉冲基函数并非一无是处。它的矩阵填充速度快在电大尺寸、只需要远场粗略估计的场合仍然可用。但如果你的目标是精确的输入阻抗或者近场分布脉冲基函数带来的伪电荷迟早会来找你麻烦。我个人的经验是用脉冲基函数算出来的电流分布在段边界处会看到明显的锯齿状抖动这个抖动就是伪电荷的直接体现。2.2 三角屋顶基函数用连续性换来的稳定三角基函数把一段电流分摊到两个相邻的段上形成一个三角形的形状在第 n 段上从 0 线性上升到峰值在第 n1 段上从峰值线性下降到 0。这样相邻基函数在公共节点处的值都是零电流在整个线上连续。def triangular_basis(node_index, nodes, n_segments): # 在第 node_index-1 段上线性上升第 node_index 段上线性下降 values [0.0] * (n_segments 1) for i in range(n_segments): if i node_index - 1: values[i] 1.0 # 峰值在 node_index 处 elif i node_index: values[i] 1.0 return values连续性的收益是实打实的。电流连续意味着电荷密度不再出现冲激而是分段常数数值上温和得多。阻抗矩阵的条件数通常比脉冲基函数低一到两个数量级迭代求解的收敛速度也更快。代价是每个基函数覆盖两个段矩阵的非零元素数量大约翻倍填充时间相应增加。三角基函数还有一个容易被忽略的好处它天然满足电流在端点处为零的边界条件。细线天线的端点电流必须为零三角基函数在第一个节点和最后一个节点处的值本来就是零不需要额外施加约束。脉冲基函数做不到这一点端点处的电流必须靠额外的方程强行压到零处理不好会引入反射。2.3 分段正弦基函数谐振长度下的精度红利与积分代价分段正弦基函数在每一段上按正弦规律变化相邻段之间保持连续。它的形状是f_n(s) sin(k(s - s_n)) / sin(k Δ)其中k是波数Δ是段长。这个形式不是随便选的它正好是细线天线在谐振附近的近似解形式。当段长接近半波长的整数分之一时分段正弦基函数能用很少的未知量把电流分布描述得很准这是它最大的卖点。代价在积分上。分段正弦基函数和格林函数做内积的时候被积函数里出现正弦和指数函数的乘积解析积分表达式冗长数值积分又需要更多的积分点才能收敛。在宽带扫频的时候k随频率变化基函数的形状跟着变每一频点都要重新计算基函数的系数填充时间比三角基函数高出一截。我的实际体会是分段正弦基函数适合窄带、高精度、结构简单的场景比如单根偶极子的输入阻抗计算。一旦结构变复杂、频带变宽它的积分代价和实现复杂度就不划算了三角基函数配合足够密的网格能达到接近的精度代码还短得多。2.4 三类标量基函数的横向对比与选型建议把三类基函数放在一张表里对照取舍会更清楚基函数类型连续性每基函数覆盖段数积分复杂度伪电荷风险适用场景脉冲不连续1低高电大尺寸粗估、远场快速扫描三角连续2中低通用细线问题、宽带扫频分段正弦连续2可扩展高低窄带高精度、谐振结构选型的时候我一般遵循一条经验规则先用三角基函数把问题跑通如果发现收敛所需未知量明显多于理论估计再考虑换分段正弦基函数压缩未知量。反过来如果矩阵规模已经大到内存吃紧才退回到脉冲基函数做粗算。这个顺序不要颠倒否则你会在伪电荷带来的各种诡异现象里浪费大量时间。3. 面元上的矢量线性基函数RWG的构造逻辑与散度约束从细线走到面元未知量从标量变成矢量基函数的构造思路要重新设计。面元上的电流是二维分布有方向还要满足法向分量的边界条件。RWG基函数是这一领域最经典的选择它的构造逻辑值得逐层拆开看因为后面处理介质界面的时候所有规则都是从它这里延伸出去的。3.1 从线到面为什么矢量基函数必须盯住法向分量面元上的电流是切向矢量但它跨越公共棱的时候法向分量的连续性才是真正要保证的东西。原因在于边界条件金属表面的电流法向分量为零介质表面的电流法向分量跳变等于表面电荷密度。如果基函数在公共棱上法向分量不连续表面电荷就会凭空产生或者凭空消失散射场的计算结果直接错掉。这就是为什么面元基函数不能用简单的每个三角形上一个常矢量来构造。常矢量在三角形边界处的法向分量取决于方向相邻三角形的常矢量之间没有任何约束法向分量必然跳变。矢量基函数必须让相邻三角形共享同一组自由度并且在公共棱上把法向分量绑死。3.2 RWG基函数的几何定义与归一化系数的来历RWG基函数定义在两个共享一条棱的三角形上公共棱的长度记为l_n两个三角形的面积分别记为A_n^和A_n^-。它的表达式是f_n(r) (l_n / (2 A_n^)) * ρ_n^(r) , r 在 T_n^ 内 f_n(r) (l_n / (2 A_n^-)) * ρ_n^-(r) , r 在 T_n^- 内ρ_n^是从T_n^的不在公共棱上的那个顶点指向观察点r的矢量ρ_n^-则是从观察点指向T_n^-对应顶点的矢量。系数l_n / (2 A_n^)不是拍脑袋定的它的作用是让基函数在公共棱上的法向分量恰好等于 1。推导过程很直接三角形的高h 2A/lρ在垂直于棱的方向上的分量正好等于这个高乘以系数l/(2A)之后法向分量变成l/(2A) * (2A/l) 1。这个归一化是整个RWG构造的核心它把基函数的物理意义固定成沿棱的单位法向电流后面的散度、电荷、连续性全都建立在这个基础上。def rwg_basis(tri_plus, tri_minus, edge_length): # tri_plus, tri_minus 为两个三角形的顶点坐标 A_plus triangle_area(tri_plus) A_minus triangle_area(tri_minus) coef_plus edge_length / (2.0 * A_plus) coef_minus edge_length / (2.0 * A_minus) return coef_plus, coef_minus写这段代码的时候有一个必须检查的点两个三角形的顶点顺序必须保证外法向一致。如果一边的顶点是逆时针、另一边是顺时针ρ的方向就会反基函数的符号跟着反散度的符号也反。这种错误在程序里不会报错只会让结果莫名其妙地偏掉排查起来非常痛苦。3.3 散度与表面电荷密度的对应关系RWG基函数的散度在两个三角形内分别是常数∇·f_n l_n / A_n^ , 在 T_n^ 内 ∇·f_n - l_n / A_n^- , 在 T_n^- 内这个散度直接对应表面电荷密度因为连续性方程给出ρ_s -(1/jω) ∇·J。把整个基函数的散度在两个三角形上积分结果是l_n和-l_n净电荷为零。这就是RWG基函数保证电荷守恒的数学根源电流从正三角形流入公共棱再从公共棱流出到负三角形流入等于流出不在任何地方堆积。这个性质在金属表面是必须的因为金属表面不能存储净电荷除非有外部注入。在介质表面则要小心介质表面允许电荷积累但积累量必须由法向场跳变决定不能是基函数构造带来的伪电荷。RWG基函数的零净电荷性质在这里反而成了一个约束介质表面的真实电荷分布要靠法向场的跳变来体现基函数本身不能自带电荷。3.4 公共棱上法向分量连续性的验证方法写完RWG基函数一定要做一次数值验证在公共棱上取若干采样点计算基函数在两个三角形侧的法向分量看它们是否都等于 1。这个验证只需要几行代码却能提前发现顶点顺序、面积计算、方向定义上的大部分错误。def check_normal_continuity(edge_point, edge_normal, tri_plus, tri_minus): f_plus rwg_value(edge_point, tri_plus, sideplus) f_minus rwg_value(edge_point, tri_minus, sideminus) n_plus dot(f_plus, edge_normal) n_minus dot(f_minus, -edge_normal) assert abs(n_plus - 1.0) 1e-10, 正侧法向分量不等于1 assert abs(n_minus - 1.0) 1e-10, 负侧法向分量不等于1实测下来最容易出问题的地方是退化三角形。当两个三角形中有一个面积趋近于零时归一化系数l/(2A)会爆炸基函数的值变成天文数字矩阵元素跟着失控。网格生成阶段就要过滤掉退化三角形或者在基函数计算时加一个面积下限保护不要指望求解器能容忍这种病态元素。4. 金属-介质复合界面上的基函数拼接J与M的连续性怎么保证到了复合结构前面的规则要叠加使用。金属表面只有等效电流J介质表面同时有J和M两种未知量的基函数展开方式相同但边界条件不同。界面上的拼接规则如果搞错方程组的解会收敛到一个物理上不存在的状态残差很小但结果全错——这是最难排查的一类问题。4.1 表面等效原理下的未知量布局按照表面等效原理介质区域被替换成等效电流和等效磁流它们分布在介质表面两侧。对于闭合介质体外表面和内表面各有一套未知量对于金属介质复合结构金属表面和介质表面的交界处要特别处理。一个常见的布局方式是金属表面只布置J的RWG基函数介质表面同时布置J和M的RWG基函数。金属和介质的交界棱上介质侧的J基函数必须满足法向分量为零的约束因为这条棱同时属于金属表面金属表面不允许电流有法向分量。4.2 公共棱上等效电流与等效磁流的连续性约束介质表面的等效电流切向连续、法向跳变跳变量对应表面电荷等效磁流的情况类似切向连续、法向跳变对应等效磁荷。用RWG基函数展开的时候切向连续性由基函数本身保证法向跳变由散度体现这是RWG构造的自然结果不需要额外处理。真正需要额外处理的是介质和金属的交界棱。这条棱上的介质侧基函数其法向分量在棱上必须为零。实现方式有两种一种是直接把这条棱从介质侧的基函数列表中剔除另一种是在基函数的系数上施加约束。前一种方式更干净但会导致交界处的电流展开不完整后一种方式保留了自由度但需要额外的约束方程。我倾向于第一种方式并且在网格划分阶段就把金属-介质交界棱标记出来让基函数生成器在遍历棱的时候直接跳过这些棱。这样生成的基函数集合天然满足边界条件矩阵方程不需要额外的约束行条件数也更稳定。4.3 跨界面棱的基函数归属规则一条棱可能同时属于金属三角形和介质三角形也可能只属于其中一种。归属规则要写清楚不然会出现同一个自由度被重复计数或者完全遗漏的情况。我的做法是给每条棱打一个类型标签棱的类型两侧三角形J基函数M基函数金属-金属金属 金属保留不生成介质-介质介质 介质保留保留金属-介质金属 介质保留仅金属侧不生成边界棱单侧不生成不生成金属-介质棱上的J基函数只保留金属侧的那一半介质侧不生成对应的J基函数。这样处理之后交界处的电流只有切向分量法向分量自动为零。这个规则看起来简单但如果没有在数据结构层面明确表达代码里很容易出现这条棱到底算谁的的混乱。4.4 体极化电流基函数与表面基函数的耦合如果介质区域不采用表面等效原理而是采用体等效原理介质内部会出现体极化电流需要体基函数来展开。体基函数通常是四面体上的线性基函数或者常矢量基函数它和表面基函数之间的耦合项要单独计算。体基函数和表面基函数的耦合本质上是体积分和面积分的混合。数值实现的时候体积分的积分点要避开表面附近因为表面附近的场有奇异性面积分的积分点要避开棱附近理由相同。我一般会在体积分和面积分之间留一个薄层的缓冲缓冲层内的贡献用解析近似补上这样能显著降低积分点数量同时保证精度。5. 数值落地时的四大坑奇异积分、条件数、低频崩溃与伪电荷基函数的理论再漂亮落到代码里还是会遇到一堆具体问题。这一节把我踩过的、见过的坑集中列一下每个坑都给出可操作的排查思路而不是一句注意处理就带过。5.1 自阻抗项里 1/R 奇异积分的处理手法RWG基函数的自阻抗项也就是同一个基函数和自身做内积的那一项被积函数里有1/R的奇异性。当观察点和源点重合时1/R发散数值积分直接算会得到无穷大或者剧烈振荡的结果。常见的处理手法有三种。第一种是奇异性提取把1/R拆成奇异部分和正则部分奇异部分用解析公式算正则部分用数值积分。第二种是Duffy变换把三角形上的积分区域映射到一个正方形让奇异点落在正方形的角上雅可比行列式在奇异点附近趋近于零正好抵消1/R的奇异性。第三种是局部坐标变换把积分点集中到奇异点附近用更密的采样处理。def duffy_transform(u, v): # 把单位正方形映射到三角形奇异点落在 u0, v0 xi u eta u * v jacobian u return xi, eta, jacobian我实测下来Duffy变换的实现最简单精度也够用缺点是积分点数量要比普通高斯积分多一些。对于自阻抗项我一般用 7 点或者 9 点的Duffy积分邻接项用普通高斯积分远场项用低阶积分这样能在精度和速度之间取得不错的平衡。5.2 基函数阶数、网格密度与矩阵条件数的三角关系基函数的阶数、网格的密度、矩阵的条件数这三者之间有一个不太直观的关系。通常的直觉是网格越密精度越高但条件数会跟着恶化。基函数阶数越高单个基函数覆盖的范围越大条件数反而可能改善因为未知量之间的耦合更平滑。具体到RWG基函数网格密度一般取每波长 8 到 10 个三角形边长对应每波长大约 10 到 20 个基函数。低于这个密度基函数无法分辨电流的快速变化精度掉得很快高于这个密度条件数上升迭代求解的收敛速度变慢而且收益递减。网格密度每波长基函数数相对精度条件数量级建议用途5 到 8低10^2快速预估、优化循环内层10 到 15中10^3通用仿真20 到 30高10^4 到 10^5高精度校核、近场分析大于 40收益递减10^6 以上一般不推荐这张表里的数字是经验值具体问题会有偏差但量级关系是稳定的。如果你发现条件数比表里高出一到两个数量级先检查基函数的顶点顺序和面积计算再检查网格里有没有细长三角形这两个原因占了绝大多数。5.3 低频崩溃现象与基函数选择的关系低频崩溃是电场积分方程的老问题频率趋近于零的时候矩阵的条件数急剧恶化求解器要么不收敛要么收敛到一个错误的结果。根源在于电场积分方程里电流项和电荷项的量级差异频率越低电荷项占比越大矩阵的两个子块量级悬殊数值上接近奇异。从基函数的角度看低频崩溃和基函数的散度性质直接相关。RWG基函数的散度是常数对应电荷项这部分在低频时主导矩阵电流项对应的部分在低频时趋近于零。解决低频崩溃的经典方法是把电流分解成环电流和树电流环电流无散度树电流有散度两者分开处理条件数就能恢复到正常量级。还有一种是增广型公式在方程里额外引入一个标量未知量把电荷项显式分离出来。这种方法实现起来比环树分解简单但矩阵规模会增大适合中等规模的问题。我自己的经验是频率低于第一谐振点十分之一的时候就要考虑低频崩溃了不要等到求解器报错才回头处理。5.4 伪电荷的识别与自检指标伪电荷是基函数不连续带来的数值产物识别它有几个实用的自检指标。第一个指标是电流分布的抖动如果沿一条直线取电流采样发现相邻采样点之间有规律的锯齿状振荡振荡的周期和网格周期一致基本可以确定是伪电荷。第二个指标是电荷密度的积分把散度在整个表面上积分理论上应该为零对于闭合金属表面如果积分结果明显偏离零说明基函数的电荷守恒被破坏了。第三个指标是能量守恒计算辐射功率和输入功率的比值如果比值明显大于 1说明有数值源在注入能量伪电荷是常见嫌疑之一。这三个指标我一般会写成自动化检查脚本每次改完基函数代码就跑一遍比等结果出来再排查省事得多。6. 一套可复现的验证流程从单棱电流到介质加载平板理论讲完最后给一套可以照着做的验证流程。这套流程分四步每一步都有明确的预期结果和判据适合刚写完基函数模块、还不确定对不对的阶段。6.1 单根细线天线的电流分布校验第一步用细线天线验证标量基函数。取一根半波长的细线用三角基函数展开电流计算输入阻抗。预期结果是输入阻抗的实部在 70 到 75 欧姆附近虚部接近零。如果实部偏离超过 10%先检查基函数在端点处是否为零再检查积分点的数量。这一步的意义在于把标量基函数的正确性单独隔离出来。如果细线天线的结果都不对后面的面元基函数和复合结构根本不用往下做先把这里的问题解决干净。6.2 金属平板的表面电流对比第二步用金属平板验证RWG基函数。取一块边长约一个波长的方形平板平面波正入射计算表面电流分布。预期结果是电流在平板中心区域比较均匀边缘处出现绕射引起的增强或者减弱。把计算结果和文献里的经典曲线对比主瓣范围内的相对误差应该控制在 5% 以内。这一步最容易暴露的问题是顶点顺序和法向方向。如果电流分布整体翻转了180度或者出现明显的对称性破缺几乎可以肯定是顶点顺序的问题。我的做法是在网格生成之后、基函数计算之前统一把所有三角形的顶点顺序调整成逆时针从源头上避免这个问题。6.3 介质块加载后的场分布检查第三步在金属平板旁边放一块介质块验证复合界面的基函数拼接。预期结果是介质块内部的场强明显低于外部介质表面的等效电流在金属-介质交界棱处法向分量趋近于零。检查方法是沿交界棱取若干采样点计算电流的法向分量绝对值应该小于切向分量的 1%。如果法向分量明显不为零说明交界棱的基函数归属规则没有正确实施。回头检查第 4.3 节那张表逐条核对棱的类型标签和基函数的生成条件问题一般出在金属-介质棱被错误地当成了介质-介质棱。6.4 收敛性扫描的实操参数设置最后一步做收敛性扫描确认结果的稳定性。扫描三个维度网格密度从每波长 8 个基函数增加到 20 个积分点数量从 3 点增加到 7 点基函数阶数从线性增加到二阶。每一次扫描记录输入阻抗或者散射截面的变化量。判据是网格密度增加一倍结果变化小于 2%积分点增加一倍结果变化小于 0.5%。如果两个判据都满足说明离散误差和积分误差都已经收敛结果可以采信。如果某一个维度怎么加密都不收敛问题大概率不在数值精度而在基函数的构造或者边界条件的实施需要回到前面的章节重新检查。这套流程跑下来大概需要半天到一天的时间取决于你对代码的熟悉程度。我的建议是把它固化成脚本和基函数代码放在同一个仓库里每次修改基函数都跑一遍。基函数是整个矩量法程序的地基地基出问题上面盖什么都是歪的。这一点我在实际项目里体会太深了宁可多花时间在验证上也不要等到结果对不上再从头排查。
