2-47 三维 FDTD 算出的近场，怎样外推成完整远场方向图？

MATLAB 3D FDTD + Love（拉夫）等效原理时谐场近远场变换程序

一句话介绍： 这套 MATLAB 程序在 80×80×80 三维 Yee 网格中迭代电磁场，以中心 z 向正弦电场作为激励，在吸收边界内侧建立封闭外推面，提取最后一个周期的切向电磁场并转换为等效电流、等效磁流，再对目标频率做复数累积和六面远场积分，最终得到归一化二维剖面、三维辐射方向图以及方向系数 D。

先看它解决什么问题

三维 FDTD 最擅长的是在有限计算区域里求电磁场随时间怎样传播，但真正关心辐射问题时，往往还想知道：

某个方向上的远区场有多强？整个球面上的辐射分布是什么形状？

水平剖面和垂直剖面分别长什么样？能不能进一步算出方向系数？

如果直接把 FDTD 区域一直扩展到“足够远”的远场位置，三维网格数量和计算量都会迅速增加。

当前程序采用的思路是：

三维 FDTD 时域迭代 → 封闭外推面采样 → Love 等效电流/磁流 → 单频傅里叶提取 → 六面空间积分 → 球坐标远场分量 → 归一化方向图 → 方向系数 D

所以它真正完成的是：

在有限三维 FDTD 区域内先得到近场，再利用封闭等效面把指定频率的时谐场外推到远区。

1. 当前工程到底是什么？

建立三维 Yee 网格并初始化 E、H 六个场分量；

在时域内循环更新电磁场，并处理六个外边界；

在吸收边界内侧的封闭六面体上提取切向场，形成等效电流和等效磁流的目标频率复幅值；

对不同 θ、φ 方向执行六面积分，得到归一化远场方向图和方向系数。

当前工程没有介质模型文件、天线几何文件、测量数据或额外函数，因此所有项目事实都可以直接回到这一份脚本核对。

2. 默认三维 FDTD 网格有多大？

源码使用自由空间参数，激励频率为：f = 100 MHz

由源码参数直接复算：

因为 dx=λ/40，80 个网格在单个方向上的总尺度约为 2λ。

时间步对应的 Courant 数为：c·dt/dx = 0.5

对于三个方向步长相同的三维 FDTD，常用稳定上限为 1/√3 ≈ 0.577，因此当前 0.5 位于这一稳定条件之内。

3. 电磁场在三维网格里怎样循环更新？

程序采用典型的 Yee 交错网格结构：

ex、ey、ez 分别保存三个方向的电场分量；

hx、hy、hz 分别保存三个方向的磁场分量；

E 和 H 在空间位置以及时间上都交错半个步长。

每个时间步先利用当前电场旋度更新 H，再利用新的磁场旋度更新 E。

从代码结构上可以概括为：

Eⁿ → 更新 Hⁿ⁺¹ᐟ² → 更新 Eⁿ⁺¹ → 处理边界 → 提取外推面场值

这种交错关系不仅决定了 FDTD 本身怎样推进，后面做目标频率傅里叶累积时也必须补偿 E、H 之间的半时间步差异。

4. 当前程序实际使用的激励源是什么？

每个时间步更新完电场后，源码都会在计算区域中心执行：

ez(ia,ja,ka) = ez(ia,ja,ka) + sin(2πf·n·dt)

也就是说，当前真正参与运行的是：

位于三维网格中心的 z 向单点正弦电场激励。

频率就是前面设置的 100 MHz。

源码前部还保留了 r0、c1、c2、ll 等细导线相关参数，并写有一整组局部场修正公式，但这些语句目前全部被注释掉，没有参与当前计算。

因此阅读这套程序时需要区分：

当前有效模型 = 中心 z 向正弦 Ez 激励

而不是“已经启用了半径 r0 的完整细导线模型”。

5. 六个外边界怎样减少反射？

为了让有限网格尽量模拟开放空间，程序在 x、y、z 六个方向都设置了外推型吸收边界。

这也是源码为什么同时保留：

ex_1、ex_2、ex_3，ey_1、ey_2、ey_3，ez_1、ez_2、ez_3

以及对应三组历史磁场的原因。

边界更新使用前三个时刻、由内向外相邻网格的场值进行多时刻外推。例如可以概括成：

边界当前值 ≈ 3×前一时刻内层值 - 3×前两时刻更内层值 + 前三时刻再内层值

它的目的不是产生新的物理场，而是尽量把传播到计算域边缘的波继续“向外延伸”，降低硬截断带来的反射。

需要注意，当前代码并没有使用 PML，而是这种外推式吸收边界。因此如果把它用于高精度定量辐射计算，边界反射误差仍然需要单独验证。

6. 外推面为什么要放在吸收边界里面？

程序通过：SX=SY=SZ=1+AA

以及：XX=YY=ZZ=81-AA

定义六个封闭外推面。

当前 AA=2，所以这个闭合面位于吸收边界内侧约：

2dx ≈ 0.1499 m的位置。

这样处理有两个目的：

外推面仍位于 FDTD 的有效计算区域里，可以直接取得近场；

又不会直接贴在最外层吸收边界上，避免把边界更新值直接当作等效面场值。

由于 Yee 网格中的 E、H 分量并不位于同一个空间位置，程序还会对外推面附近的场做平均：

电场通常对相邻两个网格值取平均；

磁场通常对相邻四个网格值取平均。

这一处理的核心目的就是：

先把交错存储的电磁场重新对齐到等效面上的同一积分位置。

7. Love 等效原理在这里具体做了什么？

当前程序的近远场变换核心不是直接拿 E、H 去积分，而是先把封闭面上的切向场转换成等效面电流和等效面磁流。

对应的基本关系可以写成：

Jₛ = n × H

Mₛ = -n × E

其中 n 是各个外推面的外法向方向。

源码中的 Jz_fx1、Jy_fx1、Mz_fx1、My_fx1 等变量，就是针对 x、y、z 六个面的不同切向分量，按照面法向分别处理正负号后得到的等效源。

简单理解：

只要封闭面上的切向电磁场已经知道，就可以把这个面看成一层等效辐射源，再由它重建外部远场。

这就是当前程序能够把有限 FDTD 近场区域和远区方向图连接起来的关键。

8. 为什么总共算 5 个周期，却只对最后 1 个周期做傅里叶提取？

源码只有在：n > TT×(SS-1)时才开始累加外推面上的频域量。

当前 TT=80、SS=5，因此：

前 4 个周期主要让时域场逐步建立；

第 5 个周期的 80 个时间步用于提取 100 MHz 的复幅值与相位。

对由 E 构成的等效磁流，程序使用：

exp[-j2π(n-N₀)fdt]

进行复数累积；对由 H 构成的等效电流，则使用：

exp[-j2π(n-1/2-N₀)fdt]

其中半步修正正是为了补偿 Yee 算法中 H 相对 E 的半时间步交错。

所以这里并不是调用 MATLAB 的 fft 对整段波形做宽频谱分析，而是：

只针对已知目标频率 100 MHz，直接累加这一频点的复数幅值和相位。

9. 六个等效面怎样被积分成远区场？

时域循环结束后，程序建立远场观察角：

θ：0～π，共 41 个采样点，步长 4.5°；

φ：0～2π，共 41 个采样点，步长 9°。

对每一个 (θ,φ) 方向，程序都会构造空间相位项：

exp[jk(x sinθ cosφ + y sinθ sinφ + z cosθ)]

再把六个面上的等效电流、等效磁流与这个相位项相乘，并按网格面积 dx² 做离散求和。

于是每个观察方向都会得到来自六个面的综合贡献。

代码随后把这些贡献组合成：

fx、fy、fz：等效电流相关的三个笛卡尔分量；

fmx、fmy、fmz：等效磁流相关的三个笛卡尔分量。

这一步正是“近场到远场”计算量的核心：

六个封闭面上的复数场分布 → 对每个观察方向进行带相位的面积积分 → 得到该方向的远区复场量。

10. 为什么不直接把 FDTD 区域扩展到真正的远场位置？

因为三维 FDTD 的计算量对网格规模非常敏感。

如果每个方向都扩展几倍距离，三维单元总数会按体积快速增长，同时每一个新增网格还要在全部时间步里更新 E、H。

当前方法只需要在一个较小的有限区域里完成时域计算，然后保存一个闭合外推面上的场，就能继续求不同角度的远区分布。

因此这条链路更适合：

近场用 FDTD 精确推进，远场用等效原理和积分外推。

这也是近远场变换比“把整个远区都网格化”更有计算价值的地方。

11. 程序怎样得到归一化方向图和方向系数？

六面积分完成后，程序先把笛卡尔分量投影到球坐标方向，得到：

fθ、fφ、fmθ、fmφ

随后按照当前源码组合远场强度：

r = |Zfθ + fmφ|² + |-Zfφ + fmθ|²

再执行：

r = r / max(r)

所以最终用于绘图和方向系数计算的 r 是一个最大值归一化为 1 的相对辐射强度方向图，不是带绝对单位的 V/m 远场值。

方向系数采用：

D = 4π / Σ[r(θ,φ) sinθ Δθ Δφ]

其中当前 Δθ=π/40，Δφ=2π/40。

这里有一个会影响定量结果的源码细节：phi=(0:n)/n*2*pi 同时包含 0 和 2π，而后面的积分循环又把 41 个 φ 点全部计入，因此周期端点会被重复积分一次。若要把 D 用于更严格的定量比较，建议积分时去掉重复的 2π 端点，例如让 φ 只取 0～2π-Δφ。

12. 当前程序最终能输出什么？

程序计算结束后会得到三类可视化结果：

一个固定方位下的 θ 向极坐标剖面；

一个 θ=90° 平面内随 φ 变化的极坐标剖面；

由 r(θ,φ) 转换得到的三维归一化辐射方向图。

三维方向图使用：

x = r sinθ cosφ

y = r sinθ sinφ

z = r cosθ

把每个角度的归一化强度直接映射成球坐标半径，再通过 surf 绘制空间方向性。

除此之外，方向系数 D 保留在 MATLAB 工作区，toc 会输出本次脚本运行耗时。

当前附件没有提供配套运行截图或保存下来的数值结果，因此这里不虚构具体的 D 数值；文章中的网格尺度、时间步、角度分辨率等数值均可由当前源码直接复算。

13. 它适合用在哪里？

三维 FDTD 近远场变换教学

这套程序把 Yee 网格更新、吸收边界、等效面采样、Love 等效电流/磁流、单频复幅值提取和远场积分放在同一份脚本里，非常适合理解“时域近场怎样一步步变成频域远场方向图”。

时谐辐射方向图算法验证

当前激励频率固定且已知，程序只提取这一频点的复幅值，因此适合验证单频近远场外推流程、角度采样和方向系数计算。

天线或散射问题的前期原型

如果后续把中心单点激励替换为真实天线、散射体或介质模型，并保证外推面完整包围目标区域，当前六面等效积分框架可以继续作为近远场变换基础。

不过需要明确：当前版本本身并不是完整工业电磁仿真器。它使用外推式吸收边界而非 PML，实际启用的是中心单点 Ez 激励，而且方向系数积分还存在 φ 周期端点重复计数问题，因此更适合教学、算法理解和原型验证，而不是未经校核直接用于高精度工程定值。

14. 哪些参数最值得修改？

当前脚本里真正会明显改变计算行为的参数主要有下面几个。

尤其要注意最后一点：当前程序把 dx 直接绑定到波长。如果要做“固定物理尺寸模型”的频率扫描，就不能只循环修改 f，还需要重新设计网格、模型尺寸与材料参数之间的关系。

15. 怎么运行？

把 threeD_sourceE_NearToFar.m 放在 MATLAB 当前工作目录中，直接运行该脚本即可。

这是完整三维时域循环，主网格为 80×80×80，同时还保存多组历史 E/H 数组，内存和运行时间都明显高于二维 FDTD；

代码每个时间步都执行 pause(.01)，它会额外增加等待时间，但不参与电磁场数值计算；

如果要把当前程序改成新的天线、散射体或介质问题，不能只改绘图部分，必须同时保证激励、材料/几何模型、外推面位置和吸收边界仍然满足新的物理问题。

对于第一次运行，最合适的做法是先保持原始参数不变，确认二维剖面、三维方向图和工作区变量 D 都能正常生成，再逐项修改 AA、SS、角度离散数或网格规模做对比。

16. 一句话看懂这个项目

这是一个三维 FDTD 时谐场近远场变换程序：它在 80×80×80 Yee 网格中以 100 MHz 的中心 z 向正弦 Ez 源进行 5 个周期、共 400 个时间步的电磁场迭代，在距离吸收边界 2 个网格的封闭六面体上对交错 E/H 场进行空间对齐，按照 Love 等效原理构造面电流和面磁流，只对最后一个周期提取目标频率的复幅值，再对 41×41 个 θ、φ 方向执行六面带相位积分，最终得到归一化二维剖面、三维辐射方向图和方向系数 D；当前版本最适合用于 3D FDTD、单频近远场外推与等效原理的教学和原型验证。