简介这是一份基于C实现的时域有限差分FDTD微带天线仿真代码库适合电磁场与天线方向的研究生、工程师以及FDTD初学者研读。源码围绕主程序展开并附带可执行文件与配套数据结果覆盖网格定义、场迭代更新、PML边界条件、S参数和辐射方向图输出等核心环节用户可修改基板介电常数、厚度或馈线尺寸研究不同参数对天线增益、带宽及方向性的影响。压缩包共12个文件以cpp源码、exe可执行程序、txt结果记录以及ild/ilc/tds/csm等辅助数据文件为主整体仅363KB结构精简。目前已有264人学习下载。代码包含可运行程序与输出文件既便于入门者对照学习FDTD实现细节也可作为微带天线仿真实验的参考基座适用于课程设计、毕业课题或工程预研中的快速建模与验证。1. 为什么微带天线仿真要自写 FDTD C 而不是直接开软件微带天线设计者拿到一个 2.45 GHz 矩形贴片需求通常会先用近似公式算长宽再用商业全波仿真软件验证。近似公式在薄基板、低介电常数时够用但遇到探针馈电、缝隙耦合、短路针或多层介质结构就明显失真商业软件能处理这些几何细节麻烦的是批量扫参、联合仿真以及把电磁内核嵌入自建工具链时授权和接口都会限制发挥空间。FDTD 把麦克斯韦方程组在时域显式推进一次宽频脉冲激励就能覆盖整个频段而微带天线恰好是窄带结构一次运行拿到全频带响应再取谐振点比逐频点求解更适合设计初期的快速评估。用 C 实现 FDTD 的优势在于更新循环全是数组运算没有隐式矩阵求解内存峰值和单步耗时都可以预估代价是激励源、吸收边界、材料填充和端口提取都需要自己写。这篇文章按网格设计、更新循环、微带建模、结果验证的顺序讲重点放在最容易算错或跑偏的环节。2. FDTD C 的网格与内存设计Yee 布局、CFL 条件和 std::vector 组织2.1 Yee 网格的场分量布置金属边界要落在哪条棱边上FDTD 的第一步不是写方程而是决定电场和磁场在空间上的存放位置。Yee 网格把电场放棱边中点、磁场放面心E_x 位于 (i1/2, j, k)E_y 位于 (i, j1/2, k)E_z 位于 (i, j, k1/2)H_x 位于 (i, j1/2, k1/2)H_y 位于 (i1/2, j, k1/2)H_z 位于 (i1/2, j1/2, k)。这个空间错位让麦克斯韦旋度方程中的空间差分变成中心差分每个差商的中点正好落在被更新场的位置上时间上再错开半步就得到不需要求解线性方程组的显式递推格式。这个错位在微带天线建模里直接决定金属掩膜怎么写。理想导体的边界条件是切向电场为零因此贴片和接地板必须用“棱边电场位置”来定义而不是用一个布尔矩阵标记整个网格单元。如果只用一个单元掩膜把所有电场分量在金属单元内清零金属边界就等于向单元内部移动了半个网格。以一个 2.45 GHz、基板介电常数 4.4、Δx 1 mm 的贴片为例半个网格对应约 0.5 mm贴片等效电长度变化约 1.5%算出来的谐振频率会偏 2% 到 4%这个量级足够让 -10 dB 带宽判断失真。所以在代码里我会把 PEC 掩膜拆成 metalEx、metalEy、metalEz 三个数组分别标记三条棱边是否位于金属内部。更新电场时按分量跳过金属内部的电场始终为零。这样金属的几何位置和 Yee 场分量的实际位置严格对应不会凭空多出或丢失半个网格的材料。2.2 网格分辨率、库朗时间步长与微带基板的波长缩短网格大小直接影响微带天线结果的可靠性。天线仿真里通常要求在最高工作频率对应的介质波长方向上放至少 20 到 40 个网格。微带基板会缩短波长所以不能直接用自由空间波长来取网格。以 FR4ε_r ≈ 4.4ε_eff ≈ 3.3为例2.45 GHz 的自由空间波长约 122 mm介质波长约 67 mm若按 λ/20 取网格步长约 3.3 mm这对 1.6 mm 厚的基板来说太粗厚度方向只有不到一格谐振频率误差会非常大。经验上要保证基板厚度方向至少有 3 到 5 个网格同时水平方向满足 λ/30 以上仿真结果才算基本可信。工作频率FR4εr≈4.4推荐网格Rogers RT5880εr≈2.2推荐网格典型时间步长Δz0.5mm2.45 GHz0.81.2 mm1.01.5 mm约 1.2 ps5.8 GHz0.40.6 mm0.50.8 mm约 0.6 ps10 GHz0.250.4 mm0.30.5 mm约 0.3 ps时间步长由库朗条件决定。三维均匀网格的稳定性条件是Δt ≤ 1 / (c_max × √(1/Δx² 1/Δy² 1/Δz²))其中 c_max 必须取计算域中的最大相速度通常就是真空光速 c0因为 PML 和空气区域的存在不允许按介质中较慢的波速放大时间步。代码实现里我会用一个显式的系数 CFLN 0.9 来控制余量double c0 2.99792458e8; double dt_max 1.0 / (c0 * std::sqrt( 1.0 / (dx * dx) 1.0 / (dy * dy) 1.0 / (dz * dz))); double dt 0.9 * dt_max;这段代码的作用是先算出满足稳定性的理论上限 dt_max再乘 0.9 留出安全余量。CFLN 不取 1.0 是因为实际网格中材料不均匀、PML 内部参数渐变接近临界值时容易出现晚期不稳定性尤其在高介电常数材料与空气交界处。取 0.9 后单步时间只会增加约 10%对总步数影响可以忽略。2.3 一维数组、介电常数掩膜和 PEC 棱边掩膜的 C 数据结构三维场最直接的做法是vectorvectorvectordouble但更新循环的索引会变成三层解引用可读性和性能都不理想。我会把所有场分量和材料参数压成一维std::vectordouble用一个内联函数把 (i, j, k) 映射到线性地址。这样循环访问最内层维度时内存是连续的缓存命中率更高而且所有掩膜数组的索引方式统一不容易把坐标搞混。#include vector #include cstddef #include cmath class FDTDGrid { public: int nx, ny, nz; double dx, dy, dz, dt; std::vectordouble Ex, Ey, Ez; std::vectordouble Hx, Hy, Hz; std::vectordouble epsR, muR; std::vectorunsigned char metalEx, metalEy, metalEz; FDTDGrid(int nx_, int ny_, int nz_, double dx_, double dy_, double dz_) : nx(nx_), ny(ny_), nz(nz_), dx(dx_), dy(dy_), dz(dz_), Ex(nx_ * ny_ * nz_, 0.0), Ey(nx_ * ny_ * nz_, 0.0), Ez(nx_ * ny_ * nz_, 0.0), Hx(nx_ * ny_ * nz_, 0.0), Hy(nx_ * ny_ * nz_, 0.0), Hz(nx_ * ny_ * nz_, 0.0), epsR(nx_ * ny_ * nz_, 1.0), muR(nx_ * ny_ * nz_, 1.0), metalEx(nx_ * ny_ * nz_, 0), metalEy(nx_ * ny_ * nz_, 0), metalEz(nx_ * ny_ * nz_, 0) {} inline size_t idx(int i, int j, int k) const { return (static_castsize_t(i) * ny j) * nz k; } void setSubstrate(int k0, int k1, double er) { for (int i 0; i nx; i) for (int j 0; j ny; j) for (int k k0; k k1; k) epsR[idx(i, j, k)] er; } };代码逻辑说明构造函数给所有场分量分配内存并初始化为零介电常数默认 1.0掩膜默认全部是空气。setSubstrate按 z 方向分层填充介质基板k0 和 k1 是基板上下表面所在的网格层号。金属掩膜按电场分量分开是因为 Yee 网格中同一个网格点上 Ex、Ey、Ez 并不在同一个物理位置一个整体单元掩膜无法准确表达金属边界。需要强调的是使用size_t作为idx的返回类型是为了避免网格数超过 int 上限时索引溢出。微带天线仿真规模通常在 100 万到 1000 万网格以内int 暂时够用但写成size_t不会损失性能还能防止未来加密网格或把代码改造成多线程时出现隐蔽的负下标问题。3. FDTD C 更新循环磁场推进、宽频激励与 S11 提取3.1 显式更新循环的最小实现先由 E 推 H再由 H 推 E写好网格布局后核心任务是把麦克斯韦旋度方程改写成显式时间递推。以一个完整时间步为例次序是先由当前时刻的电场推导下一半时刻的磁场再由更新后的磁场推导下一时刻的电场。下面是三维 FDTD 更新循环的骨架代码#include cmath #include vector const double MU0 4e-7 * M_PI; const double EPS0 8.8541878128e-12; void stepH(FDTDGrid g, double dt) { for (int i 0; i g.nx - 1; i) for (int j 0; j g.ny - 1; j) for (int k 0; k g.nz - 1; k) { size_t p g.idx(i, j, k); double invMu dt / MU0; g.Hx[p] invMu * ( (g.Ey[g.idx(i, j, k 1)] - g.Ey[p]) / g.dz - (g.Ez[g.idx(i, j 1, k)] - g.Ez[p]) / g.dy); g.Hy[p] invMu * ( (g.Ez[g.idx(i 1, j, k)] - g.Ez[p]) / g.dx - (g.Ex[g.idx(i, j, k 1)] - g.Ex[p]) / g.dz); g.Hz[p] invMu * ( (g.Ex[g.idx(i, j 1, k)] - g.Ex[p]) / g.dy - (g.Ey[g.idx(i 1, j, k)] - g.Ey[p]) / g.dx); } } void stepE(FDTDGrid g, double dt) { for (int i 1; i g.nx - 1; i) for (int j 1; j g.ny - 1; j) for (int k 1; k g.nz - 1; k) { size_t p g.idx(i, j, k); if (g.metalEx[p]) { g.Ex[p] 0.0; continue; } double invEps dt / (EPS0 * g.epsR[p]); g.Ex[p] invEps * ( (g.Hz[g.idx(i, j, k)] - g.Hz[g.idx(i, j - 1, k)]) / g.dy - (g.Hy[g.idx(i, j, k)] - g.Hy[g.idx(i, j, k - 1)]) / g.dz); // Ey、Ez 的更新模式与 Ex 轴对称一致此处省略正文不展开 } }这段代码要注意几个细节磁场更新循环不检查金属掩膜因为 PEC 内部磁场并不为零只是电场被强制为零电场更新时才需要按 metalEx、metalEy、metalEz 分别跳过。另一个细节是用而不要用覆盖旧值因为 FDTD 的时间递推本身就是旧值加增量如果覆盖会把历史场信息丢掉。更新结束后应该检查总能量确认没有出现指数发散这个验证方法会在最后一章给出具体实现。3.2 宽频激励源微带天线用软源电压缝隙而不是硬点源根据不同仿真目的激励源大概有三类选择激励方式实现写法适用对象常见问题硬源Ez[p] pulse;只关心稳态场的简单结构对反射波全反射S11 完全错误软源Ez[p] pulse;天线端口、微带馈线激励源点位置要离不连续处有足够距离总场散射场平面波在边界上对 E 和 H 同时修正入射波、方向图、RCS 计算只适用于线性无源媒质微带天线最常用的是软源电压缝隙在微带线上截一个 2 到 3 个网格的小缝沿 z 方向施加 Ez 电场相当于在端口处加了一个集总电压源。用把脉冲加到现有场上电磁波从缝隙向两侧传播反射回来的波可以自由穿过源区不会像硬源那样被直接弹回。// 高斯脉冲源参数t0 为脉冲峰值时刻tau 控制脉冲宽度 double t n * dt; double t0 4.0 * tau; double pulse std::exp(-((t - t0) * (t - t0)) / (tau * tau)); // 在缝隙位置si, sj, sk沿 z 方向施加软源 for (int k sk; k sk 1; k) { size_t p grid.idx(si, sj, k); grid.Ez[p] pulse; }脉冲宽度 tau 要由最高仿真频率 f_max 决定一般取 tau 小于等于 1 / (π × f_max)否则高频段频谱能量下降太多S11 曲线在高频端会因信噪比不足而出现抖动。t0 取 4 倍 tau 是为了让脉冲从零平滑上升避免初始时刻出现阶跃旁瓣。3.3 从时域波形到 S11 参数端口电压电流和离散傅里叶变换有了激励源和更新循环接下来要把仿真结果变成天线工程师熟悉的 S11 曲线。常见做法是在馈电缝隙处记录时域电压 V(t) 和电流 I(t)电压通过对沿缝隙高度的 Ez 字段做线积分得到电流通过计算环绕源点的磁场环路得到。得到 V(t) 和 I(t) 后输入阻抗在频域的表达式是Z_in(ω) V(ω) / I(ω)S11(ω) (Z_in(ω) − Z_0) / (Z_in(ω) Z_0)离散傅里叶变换不需要在整个仿真跑完后才做可以在时间循环里递推累积。下面代码展示了在旁边存下每个频点的复电压和复电流std::vectordouble freqList; // 要评估的频点比如 2.0 GHz 到 3.0 GHz std::vectorstd::complexdouble Vf(freqList.size()), If(freqList.size()); for (int n 0; n totalSteps; n) { // 更新磁场、电场施加激励源 // 记录当前时刻的 V 和 I double V 0.0, I 0.0; // ... 从场分量计算端口电压电流 ... for (size_t f 0; f freqList.size(); f) { double phase 2.0 * M_PI * freqList[f] * n * dt; std::complexdouble rot(cos(phase), -sin(phase)); Vf[f] V * rot * dt; If[f] I * rot * dt; } } for (size_t f 0; f freqList.size(); f) { std::complexdouble Zin Vf[f] / If[f]; std::complexdouble s11 (Zin - 50.0) / (Zin 50.0); double s11dB 20.0 * log10(std::abs(s11)); // 输出频率、S11 幅度和相位 }这段代码的要点是前面频进的 DFT 结果必须除以仿真总时长归一化吗不只需将每步增量累积即可相位因子中-sin(phase)对应的是复指数e^{-jωt}的展开。Z_0 取 50 还是其他值取决于微带馈线的特性阻抗而微带线本身是否严格是 50Ω 又取决于线宽和网格取整这点在第 4 章展开讲。4. 微带天线 FDTD 建模的三个硬细节馈线宽度、PML 参数与介质边界4.1 用整数网格近似 50Ω 微带馈线宽度取整和端口阻抗校准微带天线的馈线若直接连接 50Ω 微带线就需要在网格中表示出特定宽度。以厚度 h 1.6 mm 的 FR4 为例50Ω 微带线的宽度约 3.0 mm。若 Δx 1 mm则线宽正好取 3 个网格若 Δx 0.8 mm则是 3.75 个网格取 4 格后线宽变成 3.2 mm特性阻抗大约会降到 47Ω。这个偏差在 S11 计算中会表现为端口参考阻抗失配导致测量出的 S11 与仿真不一致。处理方式有两种一是调整全局网格让线宽尽量接近整格数二是先按取整后的线宽重新算实际 Z_0并把 S11 公式里的参考阻抗改掉。我一般会先用下面的近似公式估算 50Ω 对应的线宽再决定网格大小// Hammerstad 近似适用于 w/h 1 的宽微带线 double h 1.6e-3; double er 4.4; double Z0_target 50.0; double A Z0_target / 60.0 * sqrt((er 1.0) / 2.0) (er - 1.0) / (er 1.0) * (0.23 0.11 / er); double w h * 8.0 * exp(A) / (exp(2.0 * A) - 2.0); // 先从窄带公式起步实际手算时可以直接查微带线设计表或者用开源工具算完再填进代码。关键步骤是先定网格步长再看取整后的线宽偏离目标值多少最后选择是否微调基板厚度或介电常数来补偿。如果偏差超过 3%推荐改网格而不是改材料因为材料参数一变贴片谐振频率也会跟着漂移。4.2 PML 吸收边界的层数、电导率渐变和放置距离开放空间的微带天线必须在计算域外围设置吸收边界否则辐射波会在边界处反弹回来污染 S11。常见的做法是使用 PML 吸收层。PML 参数好坏直接影响仿真动态范围PML 参数常用取值参数过大参数过小层数 npml812内存和步数浪费低频反射明显S11 出现波纹渐变阶数 m34相邻层电导率跳变大产生数值反射吸收效率不足反射系数 R01e-31e-4反射偏大需要更厚的 PML 或更复杂渐变PML 的实现方式有很多最传统的是 Berenger 分裂场 PML现代 FDTD 代码里更常用单轴各向异性 PML。无论哪种电导率分布都采用从内到外逐渐增大的曲线典型形式是σ(z) σ_max × (z / d)^m其中 z 是 PML 内部到内边界的距离d 是 PML 厚度σ_max 按理论反射系数 R0 计算。对于 C 代码我会在一个数组里预先算好每层 PML 的电导率避免更新循环中重复计算指数或幂函数const int npml 10; const double R0 1e-4; const double m 3.0; const double d npml * dz; const double sigma_max -(m 1) * log(R0) / (2.0 * 377.0 * d); std::vectordouble sigmaPML(npml); for (int k 0; k npml; k) { double ratio (k 0.5) / npml; sigmaPML[k] sigma_max * std::pow(ratio, m); }代码中的0.5偏移是为了让电导率取样在 PML 层中心而不是边界上。σ_max 公式来自平面波垂直入射的理论反射系数实际使用中还要注意 PML 外侧还需要一层 PEC 边界并且 PML 内部场分量需要做裂场或辅助变量处理。PML 离天线太近会把天线的表面波在 PML 入口处反射回来工程经验是 PML 内边界距天线至少留 0.5 个自由空间波长。4.3 介质基板表面与贴片边缘介电常数平均和棱边掩膜微带天线中基板表面是电介质与空气的交界。Yee 网格更新 Ex、Ey 时需要用到电场所在位置的介电常数。如果网格边界刚好落在基板表面而代码只用当前网格位置的 epsR那么介电常数会在空气和 FR4 之间跳变造成场分布误差。一个简单的处理是在介质边界处对介电常数做平均// 在k层和k1层之间取平均介电常数作为边界棱边位置的等效值 double epsAvg 0.5 * (epsR[idx(i, j, k)] epsR[idx(i, j, k 1)]);这种做法虽然粗糙但对矩形贴片天线已经足够。若需要更精确的结果推荐采用亚网格或者保角技术不过那会让代码复杂度上升一个量级。金属贴片边缘的处理原则是沿贴片边界的切向电场必须位于金属棱边掩膜上如果贴片边缘跑到了网格中间需要移动金属掩膜或者对掩膜做细网格处理。最终检查时看一眼电场幅度分布图若贴片边缘出现不正常的场尖峰通常就说明金属边界和电介质边界没有对齐。5. FDTD C 微带天线仿真的三个快速验证技巧5.1 网格收敛性检查写完代码后不能直接信一组仿真的 S11。最快速的验证是网格收敛测试用当前网格 Δx 跑一次再把网格加密到 Δx / 2 跑一次比较两次 S11 最低点对应的谐振频率。二者偏差小于 1% 才认为网格基本收敛偏差在 1% 到 3% 时可以接受但需要标注不确定性超过 3% 就要重新加密网格。收敛测试不需要做全频带只取谐振频点附近窄带即可能节省一半以上的计算时间。5.2 用解析公式对照贴片尺寸和谐振频率对矩形微带贴片主模谐振频率可以由半波长公式估算f_r c / (2 × (L 2ΔL) × √ε_eff)其中 L 是贴片物理长度ΔL 是由边缘场引起的等效伸长ε_eff 是有效介电常数。把用 FDTD 算出的谐振频率与这个公式对比如果偏差在 2% 以内说明贴片尺寸、介质填充和端口设置基本正确偏差偏大时优先检查基板介电常数是否赋对、贴片长度对应的是不是最长边方向的主模以及端口源是否放在了贴片的谐振点附近。5.3 用能量监视器发现 PML 不稳定和库朗条件破坏PML 参数设置错误或 dt 略超库朗条件时仿真往往在前几百步看着正常几千步后开始指数发散。这类问题很难从单点波形发现但把计算域总能量随时间画出来就非常明显。总能量定义为E_total Σ (0.5 × ε × |E|² 0.5 × μ × |H|²) × ΔV在源停止注入后总能量应单调下降因为 PML 在吸收辐射波。如果能量曲线出现上升或反复震荡不衰减说明 PML 不稳定或库朗余量不够。下面代码用于在更新循环里递推计算总能量double energy 0.0; for (size_t p 0; p grid.Ex.size(); p) { energy 0.5 * EPS0 * grid.epsR[p] * (grid.Ex[p] * grid.Ex[p] grid.Ey[p] * grid.Ey[p] grid.Ez[p] * grid.Ez[p]); energy 0.5 * MU0 * (grid.Hx[p] * grid.Hx[p] grid.Hy[p] * grid.Hy[p] grid.Hz[p] * grid.Hz[p]); } energy * grid.dx * grid.dy * grid.dz;把这行代码放进每个时间步末尾最后把 energy 序列写到文件里。判断标准是源结束后的前几百步能量应平滑下降下降速率由 PML 反射率决定如果能量在 20 步内突然上升直接检查 dt 系数是否超过 0.9或者把 PML 层数从 10 加到 12 再跑一次。另一个把能量曲线做对数坐标观察直线衰减不出折点通常说明 PML 工作正常。本文还有配套的精品资源点击获取
简介这是一份基于C实现的时域有限差分FDTD微带天线仿真代码库适合电磁场与天线方向的研究生、工程师以及FDTD初学者研读。源码围绕主程序展开并附带可执行文件与配套数据结果覆盖网格定义、场迭代更新、PML边界条件、S参数和辐射方向图输出等核心环节用户可修改基板介电常数、厚度或馈线尺寸研究不同参数对天线增益、带宽及方向性的影响。压缩包共12个文件以cpp源码、exe可执行程序、txt结果记录以及ild/ilc/tds/csm等辅助数据文件为主整体仅363KB结构精简。目前已有264人学习下载。代码包含可运行程序与输出文件既便于入门者对照学习FDTD实现细节也可作为微带天线仿真实验的参考基座适用于课程设计、毕业课题或工程预研中的快速建模与验证。1. 为什么微带天线仿真要自写 FDTD C 而不是直接开软件微带天线设计者拿到一个 2.45 GHz 矩形贴片需求通常会先用近似公式算长宽再用商业全波仿真软件验证。近似公式在薄基板、低介电常数时够用但遇到探针馈电、缝隙耦合、短路针或多层介质结构就明显失真商业软件能处理这些几何细节麻烦的是批量扫参、联合仿真以及把电磁内核嵌入自建工具链时授权和接口都会限制发挥空间。FDTD 把麦克斯韦方程组在时域显式推进一次宽频脉冲激励就能覆盖整个频段而微带天线恰好是窄带结构一次运行拿到全频带响应再取谐振点比逐频点求解更适合设计初期的快速评估。用 C 实现 FDTD 的优势在于更新循环全是数组运算没有隐式矩阵求解内存峰值和单步耗时都可以预估代价是激励源、吸收边界、材料填充和端口提取都需要自己写。这篇文章按网格设计、更新循环、微带建模、结果验证的顺序讲重点放在最容易算错或跑偏的环节。2. FDTD C 的网格与内存设计Yee 布局、CFL 条件和 std::vector 组织2.1 Yee 网格的场分量布置金属边界要落在哪条棱边上FDTD 的第一步不是写方程而是决定电场和磁场在空间上的存放位置。Yee 网格把电场放棱边中点、磁场放面心E_x 位于 (i1/2, j, k)E_y 位于 (i, j1/2, k)E_z 位于 (i, j, k1/2)H_x 位于 (i, j1/2, k1/2)H_y 位于 (i1/2, j, k1/2)H_z 位于 (i1/2, j1/2, k)。这个空间错位让麦克斯韦旋度方程中的空间差分变成中心差分每个差商的中点正好落在被更新场的位置上时间上再错开半步就得到不需要求解线性方程组的显式递推格式。这个错位在微带天线建模里直接决定金属掩膜怎么写。理想导体的边界条件是切向电场为零因此贴片和接地板必须用“棱边电场位置”来定义而不是用一个布尔矩阵标记整个网格单元。如果只用一个单元掩膜把所有电场分量在金属单元内清零金属边界就等于向单元内部移动了半个网格。以一个 2.45 GHz、基板介电常数 4.4、Δx 1 mm 的贴片为例半个网格对应约 0.5 mm贴片等效电长度变化约 1.5%算出来的谐振频率会偏 2% 到 4%这个量级足够让 -10 dB 带宽判断失真。所以在代码里我会把 PEC 掩膜拆成 metalEx、metalEy、metalEz 三个数组分别标记三条棱边是否位于金属内部。更新电场时按分量跳过金属内部的电场始终为零。这样金属的几何位置和 Yee 场分量的实际位置严格对应不会凭空多出或丢失半个网格的材料。2.2 网格分辨率、库朗时间步长与微带基板的波长缩短网格大小直接影响微带天线结果的可靠性。天线仿真里通常要求在最高工作频率对应的介质波长方向上放至少 20 到 40 个网格。微带基板会缩短波长所以不能直接用自由空间波长来取网格。以 FR4ε_r ≈ 4.4ε_eff ≈ 3.3为例2.45 GHz 的自由空间波长约 122 mm介质波长约 67 mm若按 λ/20 取网格步长约 3.3 mm这对 1.6 mm 厚的基板来说太粗厚度方向只有不到一格谐振频率误差会非常大。经验上要保证基板厚度方向至少有 3 到 5 个网格同时水平方向满足 λ/30 以上仿真结果才算基本可信。工作频率FR4εr≈4.4推荐网格Rogers RT5880εr≈2.2推荐网格典型时间步长Δz0.5mm2.45 GHz0.81.2 mm1.01.5 mm约 1.2 ps5.8 GHz0.40.6 mm0.50.8 mm约 0.6 ps10 GHz0.250.4 mm0.30.5 mm约 0.3 ps时间步长由库朗条件决定。三维均匀网格的稳定性条件是Δt ≤ 1 / (c_max × √(1/Δx² 1/Δy² 1/Δz²))其中 c_max 必须取计算域中的最大相速度通常就是真空光速 c0因为 PML 和空气区域的存在不允许按介质中较慢的波速放大时间步。代码实现里我会用一个显式的系数 CFLN 0.9 来控制余量double c0 2.99792458e8; double dt_max 1.0 / (c0 * std::sqrt( 1.0 / (dx * dx) 1.0 / (dy * dy) 1.0 / (dz * dz))); double dt 0.9 * dt_max;这段代码的作用是先算出满足稳定性的理论上限 dt_max再乘 0.9 留出安全余量。CFLN 不取 1.0 是因为实际网格中材料不均匀、PML 内部参数渐变接近临界值时容易出现晚期不稳定性尤其在高介电常数材料与空气交界处。取 0.9 后单步时间只会增加约 10%对总步数影响可以忽略。2.3 一维数组、介电常数掩膜和 PEC 棱边掩膜的 C 数据结构三维场最直接的做法是vectorvectorvectordouble但更新循环的索引会变成三层解引用可读性和性能都不理想。我会把所有场分量和材料参数压成一维std::vectordouble用一个内联函数把 (i, j, k) 映射到线性地址。这样循环访问最内层维度时内存是连续的缓存命中率更高而且所有掩膜数组的索引方式统一不容易把坐标搞混。#include vector #include cstddef #include cmath class FDTDGrid { public: int nx, ny, nz; double dx, dy, dz, dt; std::vectordouble Ex, Ey, Ez; std::vectordouble Hx, Hy, Hz; std::vectordouble epsR, muR; std::vectorunsigned char metalEx, metalEy, metalEz; FDTDGrid(int nx_, int ny_, int nz_, double dx_, double dy_, double dz_) : nx(nx_), ny(ny_), nz(nz_), dx(dx_), dy(dy_), dz(dz_), Ex(nx_ * ny_ * nz_, 0.0), Ey(nx_ * ny_ * nz_, 0.0), Ez(nx_ * ny_ * nz_, 0.0), Hx(nx_ * ny_ * nz_, 0.0), Hy(nx_ * ny_ * nz_, 0.0), Hz(nx_ * ny_ * nz_, 0.0), epsR(nx_ * ny_ * nz_, 1.0), muR(nx_ * ny_ * nz_, 1.0), metalEx(nx_ * ny_ * nz_, 0), metalEy(nx_ * ny_ * nz_, 0), metalEz(nx_ * ny_ * nz_, 0) {} inline size_t idx(int i, int j, int k) const { return (static_castsize_t(i) * ny j) * nz k; } void setSubstrate(int k0, int k1, double er) { for (int i 0; i nx; i) for (int j 0; j ny; j) for (int k k0; k k1; k) epsR[idx(i, j, k)] er; } };代码逻辑说明构造函数给所有场分量分配内存并初始化为零介电常数默认 1.0掩膜默认全部是空气。setSubstrate按 z 方向分层填充介质基板k0 和 k1 是基板上下表面所在的网格层号。金属掩膜按电场分量分开是因为 Yee 网格中同一个网格点上 Ex、Ey、Ez 并不在同一个物理位置一个整体单元掩膜无法准确表达金属边界。需要强调的是使用size_t作为idx的返回类型是为了避免网格数超过 int 上限时索引溢出。微带天线仿真规模通常在 100 万到 1000 万网格以内int 暂时够用但写成size_t不会损失性能还能防止未来加密网格或把代码改造成多线程时出现隐蔽的负下标问题。3. FDTD C 更新循环磁场推进、宽频激励与 S11 提取3.1 显式更新循环的最小实现先由 E 推 H再由 H 推 E写好网格布局后核心任务是把麦克斯韦旋度方程改写成显式时间递推。以一个完整时间步为例次序是先由当前时刻的电场推导下一半时刻的磁场再由更新后的磁场推导下一时刻的电场。下面是三维 FDTD 更新循环的骨架代码#include cmath #include vector const double MU0 4e-7 * M_PI; const double EPS0 8.8541878128e-12; void stepH(FDTDGrid g, double dt) { for (int i 0; i g.nx - 1; i) for (int j 0; j g.ny - 1; j) for (int k 0; k g.nz - 1; k) { size_t p g.idx(i, j, k); double invMu dt / MU0; g.Hx[p] invMu * ( (g.Ey[g.idx(i, j, k 1)] - g.Ey[p]) / g.dz - (g.Ez[g.idx(i, j 1, k)] - g.Ez[p]) / g.dy); g.Hy[p] invMu * ( (g.Ez[g.idx(i 1, j, k)] - g.Ez[p]) / g.dx - (g.Ex[g.idx(i, j, k 1)] - g.Ex[p]) / g.dz); g.Hz[p] invMu * ( (g.Ex[g.idx(i, j 1, k)] - g.Ex[p]) / g.dy - (g.Ey[g.idx(i 1, j, k)] - g.Ey[p]) / g.dx); } } void stepE(FDTDGrid g, double dt) { for (int i 1; i g.nx - 1; i) for (int j 1; j g.ny - 1; j) for (int k 1; k g.nz - 1; k) { size_t p g.idx(i, j, k); if (g.metalEx[p]) { g.Ex[p] 0.0; continue; } double invEps dt / (EPS0 * g.epsR[p]); g.Ex[p] invEps * ( (g.Hz[g.idx(i, j, k)] - g.Hz[g.idx(i, j - 1, k)]) / g.dy - (g.Hy[g.idx(i, j, k)] - g.Hy[g.idx(i, j, k - 1)]) / g.dz); // Ey、Ez 的更新模式与 Ex 轴对称一致此处省略正文不展开 } }这段代码要注意几个细节磁场更新循环不检查金属掩膜因为 PEC 内部磁场并不为零只是电场被强制为零电场更新时才需要按 metalEx、metalEy、metalEz 分别跳过。另一个细节是用而不要用覆盖旧值因为 FDTD 的时间递推本身就是旧值加增量如果覆盖会把历史场信息丢掉。更新结束后应该检查总能量确认没有出现指数发散这个验证方法会在最后一章给出具体实现。3.2 宽频激励源微带天线用软源电压缝隙而不是硬点源根据不同仿真目的激励源大概有三类选择激励方式实现写法适用对象常见问题硬源Ez[p] pulse;只关心稳态场的简单结构对反射波全反射S11 完全错误软源Ez[p] pulse;天线端口、微带馈线激励源点位置要离不连续处有足够距离总场散射场平面波在边界上对 E 和 H 同时修正入射波、方向图、RCS 计算只适用于线性无源媒质微带天线最常用的是软源电压缝隙在微带线上截一个 2 到 3 个网格的小缝沿 z 方向施加 Ez 电场相当于在端口处加了一个集总电压源。用把脉冲加到现有场上电磁波从缝隙向两侧传播反射回来的波可以自由穿过源区不会像硬源那样被直接弹回。// 高斯脉冲源参数t0 为脉冲峰值时刻tau 控制脉冲宽度 double t n * dt; double t0 4.0 * tau; double pulse std::exp(-((t - t0) * (t - t0)) / (tau * tau)); // 在缝隙位置si, sj, sk沿 z 方向施加软源 for (int k sk; k sk 1; k) { size_t p grid.idx(si, sj, k); grid.Ez[p] pulse; }脉冲宽度 tau 要由最高仿真频率 f_max 决定一般取 tau 小于等于 1 / (π × f_max)否则高频段频谱能量下降太多S11 曲线在高频端会因信噪比不足而出现抖动。t0 取 4 倍 tau 是为了让脉冲从零平滑上升避免初始时刻出现阶跃旁瓣。3.3 从时域波形到 S11 参数端口电压电流和离散傅里叶变换有了激励源和更新循环接下来要把仿真结果变成天线工程师熟悉的 S11 曲线。常见做法是在馈电缝隙处记录时域电压 V(t) 和电流 I(t)电压通过对沿缝隙高度的 Ez 字段做线积分得到电流通过计算环绕源点的磁场环路得到。得到 V(t) 和 I(t) 后输入阻抗在频域的表达式是Z_in(ω) V(ω) / I(ω)S11(ω) (Z_in(ω) − Z_0) / (Z_in(ω) Z_0)离散傅里叶变换不需要在整个仿真跑完后才做可以在时间循环里递推累积。下面代码展示了在旁边存下每个频点的复电压和复电流std::vectordouble freqList; // 要评估的频点比如 2.0 GHz 到 3.0 GHz std::vectorstd::complexdouble Vf(freqList.size()), If(freqList.size()); for (int n 0; n totalSteps; n) { // 更新磁场、电场施加激励源 // 记录当前时刻的 V 和 I double V 0.0, I 0.0; // ... 从场分量计算端口电压电流 ... for (size_t f 0; f freqList.size(); f) { double phase 2.0 * M_PI * freqList[f] * n * dt; std::complexdouble rot(cos(phase), -sin(phase)); Vf[f] V * rot * dt; If[f] I * rot * dt; } } for (size_t f 0; f freqList.size(); f) { std::complexdouble Zin Vf[f] / If[f]; std::complexdouble s11 (Zin - 50.0) / (Zin 50.0); double s11dB 20.0 * log10(std::abs(s11)); // 输出频率、S11 幅度和相位 }这段代码的要点是前面频进的 DFT 结果必须除以仿真总时长归一化吗不只需将每步增量累积即可相位因子中-sin(phase)对应的是复指数e^{-jωt}的展开。Z_0 取 50 还是其他值取决于微带馈线的特性阻抗而微带线本身是否严格是 50Ω 又取决于线宽和网格取整这点在第 4 章展开讲。4. 微带天线 FDTD 建模的三个硬细节馈线宽度、PML 参数与介质边界4.1 用整数网格近似 50Ω 微带馈线宽度取整和端口阻抗校准微带天线的馈线若直接连接 50Ω 微带线就需要在网格中表示出特定宽度。以厚度 h 1.6 mm 的 FR4 为例50Ω 微带线的宽度约 3.0 mm。若 Δx 1 mm则线宽正好取 3 个网格若 Δx 0.8 mm则是 3.75 个网格取 4 格后线宽变成 3.2 mm特性阻抗大约会降到 47Ω。这个偏差在 S11 计算中会表现为端口参考阻抗失配导致测量出的 S11 与仿真不一致。处理方式有两种一是调整全局网格让线宽尽量接近整格数二是先按取整后的线宽重新算实际 Z_0并把 S11 公式里的参考阻抗改掉。我一般会先用下面的近似公式估算 50Ω 对应的线宽再决定网格大小// Hammerstad 近似适用于 w/h 1 的宽微带线 double h 1.6e-3; double er 4.4; double Z0_target 50.0; double A Z0_target / 60.0 * sqrt((er 1.0) / 2.0) (er - 1.0) / (er 1.0) * (0.23 0.11 / er); double w h * 8.0 * exp(A) / (exp(2.0 * A) - 2.0); // 先从窄带公式起步实际手算时可以直接查微带线设计表或者用开源工具算完再填进代码。关键步骤是先定网格步长再看取整后的线宽偏离目标值多少最后选择是否微调基板厚度或介电常数来补偿。如果偏差超过 3%推荐改网格而不是改材料因为材料参数一变贴片谐振频率也会跟着漂移。4.2 PML 吸收边界的层数、电导率渐变和放置距离开放空间的微带天线必须在计算域外围设置吸收边界否则辐射波会在边界处反弹回来污染 S11。常见的做法是使用 PML 吸收层。PML 参数好坏直接影响仿真动态范围PML 参数常用取值参数过大参数过小层数 npml812内存和步数浪费低频反射明显S11 出现波纹渐变阶数 m34相邻层电导率跳变大产生数值反射吸收效率不足反射系数 R01e-31e-4反射偏大需要更厚的 PML 或更复杂渐变PML 的实现方式有很多最传统的是 Berenger 分裂场 PML现代 FDTD 代码里更常用单轴各向异性 PML。无论哪种电导率分布都采用从内到外逐渐增大的曲线典型形式是σ(z) σ_max × (z / d)^m其中 z 是 PML 内部到内边界的距离d 是 PML 厚度σ_max 按理论反射系数 R0 计算。对于 C 代码我会在一个数组里预先算好每层 PML 的电导率避免更新循环中重复计算指数或幂函数const int npml 10; const double R0 1e-4; const double m 3.0; const double d npml * dz; const double sigma_max -(m 1) * log(R0) / (2.0 * 377.0 * d); std::vectordouble sigmaPML(npml); for (int k 0; k npml; k) { double ratio (k 0.5) / npml; sigmaPML[k] sigma_max * std::pow(ratio, m); }代码中的0.5偏移是为了让电导率取样在 PML 层中心而不是边界上。σ_max 公式来自平面波垂直入射的理论反射系数实际使用中还要注意 PML 外侧还需要一层 PEC 边界并且 PML 内部场分量需要做裂场或辅助变量处理。PML 离天线太近会把天线的表面波在 PML 入口处反射回来工程经验是 PML 内边界距天线至少留 0.5 个自由空间波长。4.3 介质基板表面与贴片边缘介电常数平均和棱边掩膜微带天线中基板表面是电介质与空气的交界。Yee 网格更新 Ex、Ey 时需要用到电场所在位置的介电常数。如果网格边界刚好落在基板表面而代码只用当前网格位置的 epsR那么介电常数会在空气和 FR4 之间跳变造成场分布误差。一个简单的处理是在介质边界处对介电常数做平均// 在k层和k1层之间取平均介电常数作为边界棱边位置的等效值 double epsAvg 0.5 * (epsR[idx(i, j, k)] epsR[idx(i, j, k 1)]);这种做法虽然粗糙但对矩形贴片天线已经足够。若需要更精确的结果推荐采用亚网格或者保角技术不过那会让代码复杂度上升一个量级。金属贴片边缘的处理原则是沿贴片边界的切向电场必须位于金属棱边掩膜上如果贴片边缘跑到了网格中间需要移动金属掩膜或者对掩膜做细网格处理。最终检查时看一眼电场幅度分布图若贴片边缘出现不正常的场尖峰通常就说明金属边界和电介质边界没有对齐。5. FDTD C 微带天线仿真的三个快速验证技巧5.1 网格收敛性检查写完代码后不能直接信一组仿真的 S11。最快速的验证是网格收敛测试用当前网格 Δx 跑一次再把网格加密到 Δx / 2 跑一次比较两次 S11 最低点对应的谐振频率。二者偏差小于 1% 才认为网格基本收敛偏差在 1% 到 3% 时可以接受但需要标注不确定性超过 3% 就要重新加密网格。收敛测试不需要做全频带只取谐振频点附近窄带即可能节省一半以上的计算时间。5.2 用解析公式对照贴片尺寸和谐振频率对矩形微带贴片主模谐振频率可以由半波长公式估算f_r c / (2 × (L 2ΔL) × √ε_eff)其中 L 是贴片物理长度ΔL 是由边缘场引起的等效伸长ε_eff 是有效介电常数。把用 FDTD 算出的谐振频率与这个公式对比如果偏差在 2% 以内说明贴片尺寸、介质填充和端口设置基本正确偏差偏大时优先检查基板介电常数是否赋对、贴片长度对应的是不是最长边方向的主模以及端口源是否放在了贴片的谐振点附近。5.3 用能量监视器发现 PML 不稳定和库朗条件破坏PML 参数设置错误或 dt 略超库朗条件时仿真往往在前几百步看着正常几千步后开始指数发散。这类问题很难从单点波形发现但把计算域总能量随时间画出来就非常明显。总能量定义为E_total Σ (0.5 × ε × |E|² 0.5 × μ × |H|²) × ΔV在源停止注入后总能量应单调下降因为 PML 在吸收辐射波。如果能量曲线出现上升或反复震荡不衰减说明 PML 不稳定或库朗余量不够。下面代码用于在更新循环里递推计算总能量double energy 0.0; for (size_t p 0; p grid.Ex.size(); p) { energy 0.5 * EPS0 * grid.epsR[p] * (grid.Ex[p] * grid.Ex[p] grid.Ey[p] * grid.Ey[p] grid.Ez[p] * grid.Ez[p]); energy 0.5 * MU0 * (grid.Hx[p] * grid.Hx[p] grid.Hy[p] * grid.Hy[p] grid.Hz[p] * grid.Hz[p]); } energy * grid.dx * grid.dy * grid.dz;把这行代码放进每个时间步末尾最后把 energy 序列写到文件里。判断标准是源结束后的前几百步能量应平滑下降下降速率由 PML 反射率决定如果能量在 20 步内突然上升直接检查 dt 系数是否超过 0.9或者把 PML 层数从 10 加到 12 再跑一次。另一个把能量曲线做对数坐标观察直线衰减不出折点通常说明 PML 工作正常。本文还有配套的精品资源点击获取
)
)
