西安交通大学迟明珺、马西奎 等：适用于有源区域电磁波时程精细积分仿真的大时间步长矩阵数值积分方法

阅读提示：本文约 2900 字

时程精细积分法是一种全波电磁数值方法，具有相当宽松的数值稳定性条件，且数值色散误差几乎不受时间步长大小的影响。在有源区域的电磁数值模拟中，时程精细积分法需要计算激励源添加引入的非齐次项所产生的矩阵积分，其计算或存在矩阵求逆运算易不稳定，或难以在大时间步长下保证计算结果的精度，这使得时程精细积分法的推广应用受到一定限制。本文将响应矩阵法引入有源区域的电磁波的时程精细积分法仿真中，不仅避免了矩阵的求逆运算，并且能够在保证结果计算精度的条件下，增大时间步长的选择范围，提高算法的计算效率。

研究背景

时域有限差分（Finite Difference Time Domain, FDTD）法是一种简单直观的全波电磁数值方法，已被广泛应用于电磁学研究的众多领域中。然而，根据Courant Friedrichs Lewy（CFL）稳定性条件，FDTD法的时间步长受限于空间网格尺寸却成为制约提高计算效率的一个瓶颈因素。

时程精细积分（Precise Integration Time Domain, PITD）法突破了CFL稳定性条件的限制，在保证计算精度的前提下，它不仅具有相当宽松的时间步长选取范围，且数值色散误差几乎不受时间步长大小的影响，这使得PITD法在提高计算效率方面具有自然的优势，近年来在计算电磁学领域中受到了人们的广泛关注。

目前，PITD法的数值稳定性分析、数值色散特性分析、吸收边界处理等关键技术已经相当成熟，掣肘其广泛应用的存储占用大的问题也得到了比较明显的缓解，其应用也已推广到运动介质、时变媒质等的电磁波问题分析中。

论文所解决的问题及意义

在有源区域的电磁仿真中，PITD法需要计算添加激励源引起的非齐次项积分——矩阵积分项。目前，PITD法中常用的矩阵积分项计算方法有两类，一类是直接近似积分方法，采用线性函数或三角函数对非齐次项进行近似后直接积分，得到相应的时域递推计算公式。它虽然具有良好的数值计算精度，但存在矩阵的求逆运算。

矩阵求逆不仅计算量大，而且容易数值不稳定，有时还会有逆矩阵不存在的情况。为了确保逆矩阵存在，马西奎等给出了去除系数矩阵线性相关关系的公式，其缺点是在实际计算中会遇到其它问题且编程难度大。

另一类是数值积分方法，例如辛普森积分方法、高斯积分方法等，这类方法虽然不涉及矩阵的求逆运算，但数值计算精度较低，在时间步长较大时尤其明显。如何避免矩阵求逆运算和在大时间步长下保证计算结果的精度，仍然是在PITD法应用中的一个亟待解决的问题。

钟万勰等提出了一种用于非齐次动力方程和非线性动力方程中Duhamel积分项计算的响应矩阵法。本文将响应矩阵法引入到电磁波PITD仿真中，不仅避免了非齐次项矩阵积分的求逆运算，并且保证了在大时间步长下计算结果的高精度，进而提高了算法的计算效率。

论文方法及创新点

PITD法的空间离散与FDTD法相似，沿着x、y和z三个方向将电场强度和磁场强度离散到Yee元胞的各条棱边与各个表面。

综合空间各节点上的常微分方程，可以写为如下矩阵形式：

本文将响应矩阵法引入到PITD法的非齐次项计算中，将非齐次项产生的矩阵积分的计算转化为一系列基本形式函数的响应矩阵的计算，成功地避免了矩阵求逆和实现了在大时间步长下的高精度计算。由于避免了矩阵的求逆运算，所以不受矩阵性态的限制，更便于扩展PITD法的应用范围。

采用PITD法分析电磁问题时，常用的激励源按随时间变化的形式分为两类，一类是周期性变化的时谐源；另一类是脉冲形式变化的脉冲源，例如，高斯脉冲源、升余弦脉冲源、微分高斯脉冲源、截断三余弦脉冲源、双指数脉冲源和调制高斯脉冲源。通常时谐源、升余弦脉冲源、截断三余弦脉冲源和双指数脉冲源都能写成e指数函数展开形式；而高斯脉冲源、微分高斯脉冲源、调制高斯脉冲源可以采用Taylor展开写成多项式函数展开形式。

因此，本文给出了这两种基本形式函数的响应矩阵的计算公式。值得指出，在非齐次项就是e指数函数或多项式函数形式的情况下，响应矩阵法可以得到矩阵积分项在计算机上的精确解。

这里，以由皮肤-脂肪-肌肉组成的二维生物组织模型为例，比较响应矩阵法和数值积分方法的存储占用和计算效率。如图1所示，分析由源点发出的400 MHz正弦波对生物组织的影响。

图1 由皮肤-脂肪-肌肉组成的二维生物组织模型

在保证计算误差控制在5 %内的前提下，采用辛普森积分方法、三点高斯积分方法，S最大可取3；采用响应矩阵法，S最大可取到60。S是△t 的倍数，取时间步长△t 为12.73 ps。如表1所示，响应矩阵法的运行时间为57.28 s，分别为高斯积分方法的14.51 %和辛普森积分方法的49.26 %。响应矩阵法的存储占用为451.65 MB，为高斯积分方法的73.71 %。显然，在保证计算精度的前提下，响应矩阵法的时间步长选取范围更加宽松、计算效率更高和存储占用也相对较低。

表1 响应矩阵法与数值积分方法的存储占用及运行时间

结论

本文将响应矩阵法引入到电磁波PITD仿真中，不仅避免了非齐次项矩阵积分的求逆运算，并且保证了在大时间步长下计算结果的高精度，进而提高了算法的计算效率。此外，为了不失一般性，本文给出了多项式函数展开形式、e指数函数展开形式和升余弦函数展开形式的PITD法递推公式。根据递推公式可知，响应矩阵法避免了矩阵求逆运算，使得PITD法不再掣制于矩阵性态，不仅提高了算法的数值稳定性，更有利于扩展其应用范围。

团队介绍

• 马西奎，教授，博士生导师，主要研究方向为电气工程中的多种物理场的耦合理论及其数值分析和软件技术，通信与电子系统中的电磁场与电磁波，非线性电路与系统中的分歧与混沌，混沌的控制、同步及其应用等。

• 迟明珺，博士研究生，研究方向为电磁场数值计算。

• 马亮，博士研究生，研究方向为电磁场数值计算。

• 朱晓杰，博士研究生，研究方向为电磁场数值计算。

本工作成果发表在2024年第21期《电工技术学报》，论文标题为“适用于有源区域电磁波时程精细积分仿真的大时间步长矩阵数值积分方法“。本课题为国家自然科学基金资助项目。

引用本文

迟明珺, 马西奎, 马亮, 朱晓杰. 适用于有源区域电磁波时程精细积分仿真的大时间步长矩阵数值积分方法[J]. 电工技术学报, 2024, 39(21): 6604-6613. Chi Mingjun, Ma Xikui, Ma Liang, Zhu Xiaojie. A Large Time Step Size Matrix Numerical Integral Method for Precise Integration Time Domain Simulation of Electromagnetic Waves in Active Regions. Transactions of China Electrotechnical Society, 2024, 39(21): 6604-6613.

阅读全文

中国电工技术学会

新媒体平台

学会官方微信

电工技术学报

CES电气

学会官方B站

CES TEMS

今日头条号

学会科普微信

大赛官方微信

☎️ 《电工技术学报》：010-63256949/6981；邮箱：dgjsxb@vip.126.com ☎️ 《电气技术》：010-63256943；邮箱：dianqijishu@126.com ☎️ 《中国电工技术学会电机与系统学报（英文）》：电话：010-63256823；邮箱：cestems@126.com ☎️ 编务：010-63256994 ☎️ 订阅：010-63256817 ☎️ “电工技术学报”微信号运营编辑：13121222619（微信同号） ☎️ “CES电气”微信号运营编辑：18500877291（微信同号）