★ FEATURED ARTICLE 简介这是一个面向二维横电波TE电磁仿真的交替方向隐式时域有限差分ADI-FDTD实现程序包适合学习时域有限差分FDTD数值方法或从事天线与波导分析等方向的工科学生及研究人员。程序基于C语言编写将隐式时间离散化为上下半步迭代在保证稳定性的同时突破传统库朗Courant稳定性条件限制。程序包共13个文件除核心的C源文件外还包含可直接运行的可执行程序以及VC工程、调试信息和预编译头文件便于按需修改参数并重新编译整体体积仅664KB。目前已有207人学习。通过研读源码与运行示例可掌握Yee网格预处理、半步电场更新及收敛迭代等关键实现细节并以此为基础扩展至复杂结构或激励源设定是高效入门ADI-FDTD算法的实用参考资料。1. 计算二维TE波的ADI-FDTD方法的程序.rar先想清楚三件事再解压计算二维TE波的ADI-FDTD方法的程序.rar这个包名的核心是 ADI-FDTD 和二维TE波。ADI-FDTD 全称 Alternating Direction Implicit Finite-Difference Time-Domain交替方向隐式时域有限差分它把传统FDTD每个时间步拆成两个半步用隐式更新绕过 Courant 稳定条件让时间步长不再被最小网格尺寸卡死。二维TE波则规定了场结构只有 Ex、Ey、Hz 三个分量适合波导、微带、光子晶体这类在横截面上传播的电磁问题。谁适合打开这个包电磁场专业学生、做天线与微波电路的工程师、算超表面和光波导的从业者通常都带着一个具体诉求细网格下FDTD跑不动想用ADI-FDTD把时间步放大。不过拿到rar后第一件事不是解压运行而是想清楚三件事程序里场分量是不是TE、时间步单位和CFL数藏在哪、有没有PML或吸收边界。下面按这个顺序把这个包拆开讲。2. 为什么二维TE波适合用ADI-FDTD场方程与两个半步分裂2.1 二维TE波只有三个分量更新关系比三维好写一个量级相对z轴的二维TE波电场只落在x-y平面磁场只有Hz。程序里要更新的分量就只有三个Maxwell方程组化成$$ \frac{\partial E_x}{\partial t}\frac{1}{\varepsilon}\frac{\partial H_z}{\partial y} $$$$ \frac{\partial E_y}{\partial t}-\frac{1}{\varepsilon}\frac{\partial H_z}{\partial x} $$$$ \frac{\partial H_z}{\partial t}\frac{1}{\mu}\left(\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x}\right) $$这三个方程和三维情况最大的区别是没有z方向导数也没有 Hx、Hy、Ez 的交叉耦合。所以二维TE波的ADI-FDTD程序结构上要比三维程序少一大半代码。rar包里如果只解决TE波那主循环通常只有两个半步每个半步解一个方向的三对角方程复杂度是O(NxNy)而不是三维的O(NxNy*Nz)。拿到包后先看数组定义。我见过的坑里有一类是标题写着TE程序其实是TM波TM波对应 Ez、Hx、Hy不仅场分量不同Yee网格上的空间交错位置也要整体旋转。只看公式不改网格跑起来很稳但算出来的场完全不对。二维TE波在Yee网格上Hz通常在网格中心Ex在竖直边上Ey在水平边上数组尺寸天然差了一格。2.2 ADI分裂的半步一个方向隐式一个方向显式传统FDTD是显式更新稳定条件限制时间步长$$ \Delta t \le \frac{1}{c\sqrt{1/\Delta x^21/\Delta y^2}} $$ADI-FDTD的思路是让每个半步只在一个方向做隐式求解。常见做法是取 τΔt/2第一个半步从n到n1/2把y方向的导数放进隐式$$ E_x^{n1/2}E_x^n\frac{\tau}{\varepsilon}\frac{\partial H_z^{n1/2}}{\partial y} $$$$ E_y^{n1/2}E_y^n-\frac{\tau}{\varepsilon}\frac{\partial H_z^n}{\partial x} $$$$ H_z^{n1/2}H_z^n\frac{\tau}{\mu}\left(\frac{\partial E_x^{n1/2}}{\partial y}-\frac{\partial E_y^n}{\partial x}\right) $$把 Hz 的表达式代进 Ex 的更新式消元后就得到关于 Ex 的三对角方程每一项只沿y方向耦合。这就是标准Thomas算法的用武之地沿着每一行x解一个一维三对角方程组。第二个半步从n1/2到n1换x方向做隐式得到关于Ey的三对角方程。两步合起来才是一个完整的物理时间步Δt。这个分裂格式的好处是每一步都不需要解大型稀疏矩阵只需要解若干条一维三对角链。代价是引入算子分裂误差时间步放得越大这个误差越明显。“无条件稳定”说的是数值上不炸不是结果一定准。2.3 “无条件稳定”不是免死金牌怎么判断自己的问题该不该用它我判断一个电磁问题要不要用ADI-FDTD按下面几步走先看网格尺寸是不是小于介质中最短波长的十分之一。如果网格太粗时间步放得再大也没有意义色散误差已经把波形抹掉了。再看结构里有没有高Q谐振、窄缝、薄介质层。这类结构对时间离散非常敏感ADI-FDTD的大时间步容易让谐振频率偏移。然后用小时间步跑一个基准算例比如CFLN1的FDTD结果再把你期望的ADI时间步放上去对比。如果只是想把时间步放大2到5倍ADI-FDTD很合适如果指望放大20倍通常会在频谱上得到一个“稳定但错误”的结果。程序包自带的readme里如果只夸“无条件稳定”没有给验证算例那拿到手后第一件事不是直接算自己的模型而是先复现它的demo结果。复现不出来后面所有参数都没有参照。3. 解开程序.rar用7-Zip看文件结构再决定是MATLAB还是Fortran3.1 rar里一般装了什么主程序、场更新、源、PML、后处理学术风格的ADI-FDTD程序.rar通常不是可执行exe而是源码加一个readme。我一般解压后先看目录文件类型大致可以分五类模块常见命名职责主程序main.m / main.f90定义网格、时间步、总步数、主循环场更新adi_update.f90 / adi_x.m / adi_y.m两个半步的三对角求解激励源source.m / guass.c加入点源、线源、电流密度边界条件pml.f90 / abc.mPML或吸收边界后处理plot_field.py / output.m导出场快照、探针时域数据别一上来就找.exe正规学术程序很少给编译好的二进制。如果rar里只有readme和一个main.f90那也正常很多程序把子过程写在同一个文件里。先看readme里写的编译器要求Intel Fortran、gfortran还是MATLAB这个信息决定你后面所有安装动作。如果有demo目录先跑demo。我遇到过不少包主程序里写死的路径是C:\Users\admin\Desktop\proj换到Linux什么都不改直接段错误。先把路径和输出文件名改成自己环境下的再谈物理参数。3.2 用7-Zip解压并把文件名乱码处理掉很多人问7-Zip能不能解压rar文件能而且足够干净。WinRAR弹广告已经是老问题7-Zip右键“提取到当前目录”就行。这里有一个容易卡住进度的细节中文文件名乱码。现象是解压后目录树全部正确但主程序、数据文件的名字是一串乱码脚本里引用的文件名对不上。原因不是压缩包损坏而是压缩包内文件名用的是GBK编码7-Zip默认按UTF-8解码。解决方法是换支持强制编码的压缩软件或者直接命令行用ls -b看真实文件名然后mv改成ASCII名字。改完再编译避免后面每一条路径都在踩坑。如果压缩包加了密码我建议直接找作者要。网上一搜“rar密码移除”“强制解压”之类的工具十个里有九个是捆绑下载双击之后不是帮你解密是先给你装全家桶。你手上的科研代码大概率不是加密内容不要在这种地方冒险。解压完用下面命令看一层结构cd adi_te_project find . -maxdepth 2 -type f | sort看到文件列表后心里先给每个文件归类归类完成后再跑程序。这一步能省掉后面很多“找不到模块”的报错。3.3 跑通最小算例先别调物理先让程序不报错很多ADI-FDTD程序是Fortran写的因为三对角求解在Fortran里可以直接用手写Thomas算法性能比MATLAB高一个量级。拿到源程序后的最小操作是# 先用小网格验证编译和运行不要直接跑1000*1000 cat input.dat EOF nx 40 ny 40 dx 0.001 dy 0.001 dt 2e-12 steps 500 EOF ifort -O2 -mkl -o adi_te main.f90 adi_update.f90 source.f90 pml.f90 ./adi_te input.dat代码里input.dat是常见做法我这样写是为了把网格、步长、步数这些参数外置避免每次都要改源码重新编译。如果你拿到的程序没有输入文件最简单的方式是打开main程序找到nx、ny、dt、steps这些赋值行改成最小值。500步是为了快速验证程序能完整跑完一个循环不考虑物理精度。跑完之后检查输出文件。一般会有探针点的时间序列比如Hz_probe.dat用文本编辑器打开看一眼前几十行如果不是全零说明主循环在执行。全零的话检查激励源有没有加对或者程序是不是用了硬源但源点坐标没落在采样范围内。这一步通过之后再调CFL数再上PML。4. 激励源、CFL数与PML三个参数决定结果可不可信4.1 高斯脉冲和调制高斯激励波形决定你能得到宽频还是窄带结果二维TE波程序里最常用的激励是高斯脉冲因为它频谱平滑一次仿真能覆盖宽频段。设最大关注频率为fmax脉冲宽度参数τ一般取$$ \tau \frac{1}{3 f_{max}} $$然后用一个时移形成波形。Python里写出来就是import numpy as np fmax 5e9 # 最大关注频率 tau 1.0 / (3 * fmax) t0 4 * tau t np.linspace(0, 5e-9, 5000) J_base np.exp(-((t - t0) / tau) ** 2) # 如果是窄带问题用调制高斯 f0 2.4e9 J_mod J_base * np.cos(2 * np.pi * f0 * (t - t0))这里J_base是基带高斯适合算宽带响应J_mod是调制高斯适合算某一个频点附近的响应。程序里如果加了硬源也就是直接把电场或磁场赋成这个值波形会带很强的非物理直流分量我一般倾向用软源通过把电流密度加进Maxwell方程实现让波场自己往外辐射避免源点处的寄生反射。看程序里的源实现时注意它加在哪个场分量上。二维TE波如果源应该加在Ez上那是TM波的程序TE波的电流源通常加在Ex或Ey上具体取决于你想要的极化方向。源位置和极化方向错了后面的PML调得再好也是白调。4.2 不要问“CFL能到多少”先算基准再用CFLN2到5对比ADI-FDTD放宽了CFL限制但“能到多少”是个伪问题。正确流程是先把CFLN定义为$$ CFLN\frac{\Delta t_{ADI}}{\Delta t_{FDTD}} $$其中Δt_FDTD就是传统FDTD的稳定性极限。然后做一组对比实验用CFLN1跑一次作为基准结果。用CFLN2、4、8分别跑同一个模型。把探针时域波形叠在一起比较。如果CFLN从2变到4波形早期还能重合后期开始分离说明网格、源、PML都没有大问题分离来自分裂误差。如果从CFLN1到2就明显偏了那就不是时间步的问题而是程序本身有bug。我的习惯是先把程序里的时间步写成CFLN形式不要写裸的Δt值。这样换成不同网格尺寸时稳定性条件会跟着变程序不容易悄悄越过CFL极限。程序包如果给的例子是CFLN50我建议你怀疑它是不是只演示“不炸”没展示“准”。频率响应里峰位偏移超过百分之几对工程来说已经不可接受了。4.3 没有PML的ADI-FDTD只是一半吸收边界怎么接进去二维TE波的开放域问题程序里必须带PML或至少一阶吸收边界。否则波到了截断边界反弹回来探针上看到的是混叠场根本没法分析。PML参数常用多项式渐变$$ \sigma(\rho)\sigma_{max}\left(\frac{\rho}{L}\right)^m $$其中L是PML层厚度对应的物理尺寸m一般取3到4σ_max的经验值约为$$ \sigma_{max}\approx\frac{m1}{150\pi \Delta} $$这里的Δ是靠近PML边界的空间步长。很多程序的PML只做了显式方向ADI-FDTD的隐式方向没把σ放进去。这样造成的现象是边界反射不大但整体波形在后期有一个缓慢的拖尾怎么都消不掉。原因就是隐式求解矩阵里的电导率项缺失PML层在隐式方向上根本没有吸收能力。拿到rar后如果发现程序里只有一层abc.f90而没有PML那要当心。一阶吸收边界在正入射时勉强能用斜入射反射很严重。ADI-FDTD的优势就是可以加大时间步如果边界条件把精度限制死了这个大时间步的优势就名存实亡。5. ADI-FDTD程序落地避坑五个从项目里爬出来的翻车点5.1 解压后中文乱码还没到编译就卡住现象用普通解压工具解开后文件夹里的源码文件名全是乱码Makefile里引用的路径对不上编译器提示找不到文件。原因压缩包在Windows下用GBK编码保存了文件名Linux和macOS下默认按UTF-8显示。rar本身没坏文件字节也没坏。解决不要急着找“rar解压软件”换着试直接用7-Zip解压后在终端里用ls -b看到转义后的真实文件名然后用mv把乱码改成main.f90这类ASCII名字。如果压缩包里所有文件都乱码可以先用convmv或Python的rarfile配合编码转换重命名。切记不要去下载那些号称“强制解压”“rar密码移除”的小工具多数是捆绑木马。5.2 CFLN开大以后波形稳定但频谱偏移现象跑了3000步场值没有发散屏幕上的波形看起来也正常但和CFLN1的基准结果对比谐振峰频率偏移了百分之五甚至出现额外的小峰。原因ADI-FDTD的无条件稳定是针对数值稳定性不是针对精度。时间步越大算子分裂误差和数值色散越严重高频分量被改变得越多。解决把CFLN压在2到5之间然后用CFLN1的结果校准一次。如果程序包只给了一个CFLN20的结果图自己复现时一定要同时保存小时间步基准不要只看“长得很像的波形”。5.3 标题是TE程序却是TM波的场分量现象把程序包里的参数直接改成自己的模型后电场分量输出的量级不对反射系数曲线完全对不上。原因有些作者会把同一个框架压缩成多个版本TE/TM的标记没有改干净。TE波要更新Ex、Ey、HzTM波要更新Ez、Hx、Hy。TE波程序里如果出现Ez数组那很可能主框架是TM波。解决打开程序搜索数组名确认主循环里更新的是哪几个分量。如果是TM波不要只改公式符号Yee网格的磁场位置也要跟着转。二维TE波的Hz在网格中心TM波的Hz不存在这是结构的改变不是参数改变。5.4 Thomas算法越界0基数组和1基数组混在一起现象程序的网格规模超过某个值后出现数组下标越界或者结果在角落出现明显亮斑。原因MATLAB和Fortran的数组默认从1开始Python/C从0开始Yee网格上Ex、Ey、Hz的位置又各偏半格写循环时边界很容易多一格或者少一格。解决先把网格缩到比如5乘5在Thomas算法入口打印矩阵非零元位置确认三对角三条对角线对应的索引没有负值和越界。这是最费时的部分但定位一次之后网格尺寸再大也不会犯同一个错。很多rar里的程序自己就有这个问题小网格没事大网格在PML边界崩。5.5 OpenMP并行时三对角求解出现随机噪声现象加了-qopenmp或-fopenmp编译后程序跑得变快但输出结果每次都有一点点差异后处理时高频噪声反复出现。原因三对角求解本身是按行独立可以并行但如果并行时把内层循环的方向写反多个线程同时读写同一个场数组就会产生竞态条件。解决固定外层并行维度。如果第一步是沿y方向解Ex循环应该写成!$omp parallel do private(j, ...) do i 1, nx call thomas(...) ! 每个i处理一串y方向的点 end do !$omp end parallel do保证每个线程处理不同的i行不会碰同一个Hz数组的不同写入位置。PML区的更新也要放在同样的并行策略下否则边界不反射但边界附近会出现随机涨落。6. 最后验收用能量衰减和Q值证明你的ADI-FDTD程序没白跑程序不报错、波形能出来还不代表能用于工程。我最后一个动作永远是做一个封闭腔或者开放腔的Q值校验。以二维TE波的PEC矩形腔为例在腔内放一个点源记录某个位置的Hz时域信号。稳定后再截掉源让场自然衰减。用Python直接算Q值from scipy.signal import welch fs 1.0 / dt f, P welch(hz_probe, fsfs, nperseg2**12) # 找峰值频率和-3dB带宽 f_r f[np.argmax(P)] half_max np.max(P) / 2 idx_half np.where(P half_max)[0] bandwidth f[idx_half[-1]] - f[idx_half[0]] Q f_r / bandwidth print(ff_r {f_r/1e9:.4f} GHz, Q {Q:.1f})理论Q值由腔体尺寸和材料损耗决定如果程序用的是无耗PEC腔Q值应该主要是数值损耗和PML泄漏的反映。Q值比理论低说明数值损耗太大多半是时间步太大或PML层数不够Q值比理论高说明边界反射把能量留在了腔内PML根本没工作。没有PML时还可以做能量守恒检查把整个域内电磁场能量累加没有激励后能量应当单调下降。ADI-FDTD如果有分裂误差能量曲线会出现小的台阶这个台阶就是每个时间步的数值误差。台阶越大越不能用大CFLN。我自己拿到这种计算二维TE波的ADI-FDTD程序包从来不走“解压后直接跑”的捷径。先把TE/TM分清楚再用小网格验证索引然后用CFLN对比最后用Q值验收。用这套流程替换掉“能跑就行”的直觉在这个方向上的时间投入才算真正落地。希望这些步骤能帮你也少走一圈弯路。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
