A study of fast DGTD method for metal sheet slits
-
摘要:
某些精细结构(如周期结构)仿真中常含有金属片狭缝结构,此类结构精确仿真需要适应狭缝网格尺寸,消耗计算资源较大。为解决此问题,本文借鉴本小组前期细导线研究结果,基于电磁学理论中的对偶原理,即细导线可等效为电流、金属狭缝可等效为磁流,在研究中将金属片狭缝等效为磁流,从而避开复杂的前处理建模过程,可以降低建模及计算难度。算例表明,本文方案精确度较高,且能大幅降低计算资源,对其进一步优化研究可为相关多尺度问题提供技术储备。
Abstract:Metal sheets in periodic structure simulation often contain slit structures, and accurate simulation of such structures consumes large computational resources. In order to solve this problem, this paper draws on the thin wire fast scheme, which equates the slit of the metal sheet to the magnetic flow, thus avoiding the modeling part of it, which can reduce the difficulty of modeling and computation. Examples show that the proposed scheme has high accuracy and can significantly reduce the computational resources, and its further optimization study can provide technical reserves for related multi-scale problems.
-
0 引 言
周期结构常用金属狭缝材料作为基础,两侧覆盖基底材料固定,此类目标中金属狭缝与整体模型相比尺度较小,可归于典型的多尺度问题。此外在金属机箱的设计中,为散热及电磁兼容的设计中也含有金属狭缝几何构造。常用的处理方法是将金属狭缝精确建模,导致狭缝决定的网格尺度远小于待分析电磁波长。虽然时域非连续伽辽金法[1-2](discontinuous Galerkin time domain method, DGTD)等非结构网格算法[3-5]可以将网格尺寸调整为狭缝附近较小、周围较大的方案,但这一方案会造成海量的资源消耗,局域时间步[6]加速效果也有限。
从电磁理论角度来看,金属狭缝与金属细导线具有对偶性,换言之,在数值仿真中可借鉴金属细导线的处理方案来处理金属狭缝。1981年Holland等提出等效电流元法[7],将细导线等效为若干段电流元,联合电流方程与Maxwell方程组进行求解。该方案中细导线通过电流元与空间电磁场相互作用,可以实现对任意取向的细导线仿真,其优势是通用性好,可用于多种算法[7-11]。本小组对其在DGTD中的应用进行了深入研究,包括复杂结构构建、交叉点问题等[2]。
基于对偶性,可将金属狭缝视为缝隙天线,该天线又可等效为线磁流激励。狭缝为一片状模型,而细导线为圆柱体,算法中的若干参数选取也不同。本文基于DGTD研究了金属狭缝的仿真问题,给出了无需建立金属狭缝精细模型的处理方案,并验证了方案的有效性。
1 金属片狭缝仿真方案
从电磁波理论角度来看,金属导线与金属狭缝有一定的对偶性,导线可等效为电流元构成,狭缝可等效为磁流元构成。虽然圆柱形导线与狭缝的几何形态不一致,但在电磁波波长较大的情况下可借鉴细导线方案,忽略狭缝的几何尺寸,将其视作磁流元从而给出一种可行的工程方案。本文讨论的金属片狭缝暂忽略金属片厚度,在算法中将其视作片状无厚度模型处理,如图1所示,其中 \hat {\boldsymbol{\xi}} 为细导线 (PEC 材料)或金属片狭缝的方向。
首先简述细导线方案,基于文献[7]理论,相关方程组为:
\left\{ \begin{aligned} & L\frac{{\partial I}}{{\partial t}} + \frac{{\partial V}}{{\partial \xi }} = {E_\xi } + {{\tilde V}^{{\text{source}}}} \\ & \frac{{\varepsilon \mu }}{L}\frac{{\partial V}}{{\partial t}} + \frac{{\partial I}}{{\partial \xi }} = 0 \\ & L = \frac{\mu }{{2{\text{π }}}}\ln \left(\frac{{{\rho _0} + a}}{{2a}}\right) \end{aligned}\right. (1) 式中: L为单位长度电感;V 及I为电压和电流; {E_\xi } 为与导线距离R处的电场; {\tilde V^{{\text{source}}}} 为电压源,该项仅在有电压激励时不为0(例如同轴激励); \varepsilon 、 \mu 分别为介电系数和磁导率; {\rho _0} 为虚拟柱的半径; a 为导线半径。
基于对偶性,金属片狭缝采用相似方程组
\left\{ \begin{aligned} & {L_{\rm{m}}}\frac{{\partial {I_{\rm{m}}}}}{{\partial t}} + \frac{{\partial {V_{\rm{m}}}}}{{\partial \xi }} = {H_\xi } + {{\tilde V}_{\rm{m}}}^{{\text{source}}} \\ & \frac{{\varepsilon \mu }}{{{L_{\rm{m}}}}}\frac{{\partial {V_{\rm{m}}}}}{{\partial t}} + \frac{{\partial {I_{\rm{m}}}}}{{\partial \xi }} = 0 \\ & {L_{\rm{m}}} = \frac{\varepsilon }{{2{\text{π }}}}\ln \left(\frac{{{\rho _0} + a}}{{2a}}\right) \end{aligned}\right. (2) 式中: {V_{\mathrm{m}}} 及 {I_{\mathrm{m}}} 分别为磁压及磁流; {H_\xi } 为与导线距离 {\rho _0} 处的磁场; {\tilde V_{\mathrm{m}}}^{{\text{source}}} 为磁压源;本文中 {\rho _0} 为平均棱边长度的1.5倍。
式(2)将狭缝等效为图1(右)中的磁流元,金属片狭缝的几何特性为片状而非圆柱体,因此在式(2)中有相应的磁导线半径a需等效处理。经本小组前期算例测试,采用如下近似表达式:
a = w/50 (3) 式中, w 为狭缝宽度。
式(2)可采用一维差分形式处理,其迭代方程为
\left\{ \begin{aligned} {I_{{\rm{m}}k}^{n + \frac{1}{2}}} =& {I_{{\rm{m}}k}^{n + \frac{1}{2}}} +\\ & \frac{{\Delta t}}{{{L_{\rm{m}}}}}\left[ {\left( {{H_{\xi k}^n } + {{\tilde V}_{\rm{m}}}^{{\mathrm{source}}}} \right) - \frac{{{V_{{\rm{m}}{k + \frac{1}{2}}}^n}- {V_{{\rm{m}}{k - \frac{1}{2}}}^n}}}{{\Delta \xi }}} \right] \\ {V_{{\rm{m}}{k + \frac{1}{2}}}^{n + 1}} = & {V_{{\rm{m}}{k + \frac{1}{2}}}^n} - \frac{{{L_{\rm{m}}}\Delta t}}{{\varepsilon \mu \Delta \xi }}\left( {{I_{{\rm{m}}{k + 1}}^{n + \frac{1}{2}}} - {I_{{\rm{m}}k}^{n + \frac{1}{2}}}} \right) \end{aligned}\right. (4) 式中, {I_{{\mathrm{m}}k}^{n + \frac{1}{2}}} 为k位置处 n + \dfrac{1}{2} 时间步磁流变量。
式(4)为一维差分时域迭代,其相关磁流、磁压变量设置如图2,W_1^1与W_1^2代表1号狭缝的两个端点坐标,两个端点构成了W_1^{}狭缝,两端分配给磁压 {V_{\mathrm{m}}} ,磁流与磁压交替分布。
狭缝算法需与DGTD算法配合,三维情形DGTD最终的线性方程组如下:
\left\{ \begin{aligned} & \varepsilon [{{\boldsymbol{M}}^{\text{e}}}]\frac{\partial }{{\partial t}}\{ {{\boldsymbol{E}}^{\text{e}}}\} - [{{\boldsymbol{S}}^{\text{e}}}]\{ {{\boldsymbol{H}}^{\text{e}}}\} {\text{ }} - \{ {\boldsymbol{F}}{{\boldsymbol{h}}^{\text{e}}}\} + \sigma [{{\boldsymbol{M}}^{\text{e}}}]\{ {{\boldsymbol{E}}^{\text{e}}}\} + \{ {{\boldsymbol{J}}^{\text{e}}}\} = 0 \\ & \mu [{{\boldsymbol{M}}^{\text{e}}}]\frac{\partial }{{\partial t}}\{ {{\boldsymbol{H}}^{\text{e}}}\} + [{{\boldsymbol{S}}^{\text{e}}}]\{ {{\boldsymbol{E}}^{\text{e}}}\} {\text{ }} + \{ {\boldsymbol{F}}{{\boldsymbol{e}}^{\text{e}}}\} + {\sigma _{\text{m}}}[{{\boldsymbol{M}}^{\text{e}}}]\{ {{\boldsymbol{H}}^{\text{e}}}\} + \{ {\boldsymbol{J}}_{\mathrm{m}}^{\text{e}}\} = 0 \end{aligned}\right. (5) 式中: \{ {{\boldsymbol{J}}^{\text{e}}}\} 为电流项,在细导线方案中为导线(电流元),本文金属狭缝可令 \{ {{\boldsymbol{J}}^{\mathrm{e}}}\} 为0; \{ {{\boldsymbol{J}}_{\mathrm{m}}^{\mathrm{e}}}\} 为磁流项,与式(4)中的磁流数据关联,表达式为
{J_{{\mathrm{m}}i}^{\text{e}}} = \iiint\limits_{{\Omega ^{\text{e}}}} {{\boldsymbol{N}}_i^{\text{e}} \cdot {{\boldsymbol{J}}_{\mathrm{m}}}{\text{d}}\varOmega } = \left( {{\boldsymbol{N}}_i^{\text{e}}\left| {_{x',y', {\textit{z}}'}} \right. \cdot \hat {\boldsymbol{l}}} \right){I_{\mathrm{m}}}l (6) 式中: \hat {\boldsymbol{l}} 为磁偶极子方向单位矢量; ( {x}^{\prime },{y}^{\prime },{ {\textit{z}}}^{\prime } )为磁偶极子位置; l 为磁流元长度(即空间步长)。由于DGTD采用非结构网格,磁流方程得到的磁流元可能在四面体的任意位置,因此我们给出的表达式可以处理任意位置、任意取向的激励,将模型与狭缝分离。
2 磁场与磁流元的耦合
在狭缝(磁导线)附近定义一个柱状区域[10],用于处理磁导线及附近磁场的耦合,如图3所示,该柱状域环绕第k个磁流元。式 (4)中的磁场采用如下近似加权计算[9]:
{H_{\xi k}^n} = \frac{{\displaystyle \int_{}^{} {\left[ {{\boldsymbol{H}}(r) \cdot \hat {\boldsymbol{\xi}} } \right]} g(r)\varPhi (\xi ){\text{d}}V}}{{\Delta \xi }} (7) 式中: g(r) 为权函数;定义域为 a < r < {\rho _0} ; \varPhi (\xi ) 为一维基函数[9-11]。 式(7)用于计算环绕某个磁流元 {I_{\mathrm{m}}} 的磁场,如图3所示,计算中每个磁流元均设置柱状区域用于处理磁场的加权计算。
式(7)中的加权计算需遵循以下步骤:
1)对第k个单元,搜索含有此单元柱状区域的四面体;
2)对式(7)进行数值积分,挑选若干积分采样点,采样点的电磁场采用四面体基函数合成计算得到;
3)利用第2步的采样点磁场结果计算数值积分。
此外,磁流元在DGTD中的激励也采用加权形式:
{{\boldsymbol{J}}_{{\text{Thin wire}}}} = \sum\limits_{}^{} {\left[ {{I_{}}_k^{n + \frac{1}{2}}(\xi )g(r)\varPhi (\xi )} \right]\hat {\boldsymbol{\xi}} } (8) 与数值积分类似,挑选柱状区域中的若干激励采样点,将磁流元按权值进行激励。
3 数值算例
3.1 直狭缝算例
采用直狭缝周期结构模型进行测试,如图4所示。该直线狭缝长0.5 m、宽1 cm,位于金属片中心,金属片为正方形,其边长0.75 m,吸收边界为完全匹配层(perfectly matched layer,PML)。由于该模型狭缝两侧为断点,因此将两侧的磁流元设为0。
入射波电场垂直于狭缝,即电场y极化。采用精确建模方案及等效近似方案进行仿真,如图5所示为纵向一维坐标模型,z轴为相关边界坐标。图6为直狭缝的入射波、透射波、反射波电场y分量时域波形。可以看出,透射波幅值较小,大部分被反射。
两种方案的计算内存需求和计算耗时如表1所示。可以看出,精确建模计算时间是等效近似的89倍。 本文DGTD算法采用1.5 阶叠层矢量基函数(自由度为20),计算中采用蛙跳时间步进方案,计算平台为Intel(R) Core(TM) i7-10700 CPU。
表 1 直狭缝两种方案计算参数Tab. 1 Line slilt calculation parameters of 2 schemes四面体单元数量 计算方法 内存需求/MB 计算时间/min 11 467 等效近似 97.69 1.00 54 042 精确建模 337.49 89.25 图7给出了傅里叶变换后的反射及透射系数。本例以精确建模方案数据为参考值,近似方案频域数据中,反射系数均方根误差为0.052 4,透射系数均方根误差为0.030 0,计算结果精确性很好。
3.2 六边形狭缝算例
六边形狭缝是一种较为常见的周期结构形式,其精确模型和近似模型如图8所示。六边形狭缝边长0.2 m,狭缝宽度1 cm,位于金属片中心,金属片为正方形,其边长0.75 m,吸收边界也为PML。
入射波电场取y极化方向。采用精确建模方案及等效近似方案进行仿真,相关边界坐标与图5相同。图9为六边形狭缝的入射波、透射波、反射波电场y分量时域波形。六边形狭缝结果与直狭缝结果相似,即透射波幅值较小,大部分被反射。
精确建模及等效近似方案的计算内存需求和计算耗时如表2所示。可以看出,本例中近似方案同样计算时间远小于精确方案。
图10给出了傅里叶变换后的反射及透射系数,同样以精确方案数据为参考,近似方案反射系数均方根误差为0.040 9,透射系数均方根误差为0.029 7,计算结果精确性很好。
表 2 六边形狭缝两种方案计算参数Tab. 2 Hexagonal slit calculation parameters of 2 schemes四面体单元数量 计算方法 内存需求/MB 计算时间/min 11 467 等效近似 97.69 1.02 78 620 精确建模 445.92 156.33 4 结 论
本文基于DGTD算法,讨论了对于金属片狭缝的快速仿真方案,即将狭缝等效为磁流处理,参考了与其具有对偶性的细导线方案。仿真结果表明本文的等效近似方案可以达到很好的精度,并大幅提高了计算效率,无需建立精确狭缝模型,具备工程应用前景。下一步工作将聚焦于研究此类近似模型的理论推导工作及复杂情形的适配工作,使其具备在多尺度分析中的实战能力。
-
表 1 直狭缝两种方案计算参数
Tab. 1 Line slilt calculation parameters of 2 schemes
四面体单元数量 计算方法 内存需求/MB 计算时间/min 11 467 等效近似 97.69 1.00 54 042 精确建模 337.49 89.25 表 2 六边形狭缝两种方案计算参数
Tab. 2 Hexagonal slit calculation parameters of 2 schemes
四面体单元数量 计算方法 内存需求/MB 计算时间/min 11 467 等效近似 97.69 1.02 78 620 精确建模 445.92 156.33 -
[1] ANGULO L D,ALVAREZ J,PANTOJA M F,et al. Discontinuous Galerkin time domain methods in computational electrodynamics:state of the art[C]//Forum for Electromagnetic Research Methods and Application Technologies,2015,10:1-24.
[2] 魏兵, 陈娟, 何欣波, 等. 电磁场时域数值算法的新进展[J]. 电波科学学报,2020,35(1):55-68. doi: 10.13443/j.cjors.2019090206 WEI B, CHEN J, HE X B, et al. New development of time-domain numerical algorithms for electromagnetic field[J]. Chinese journal of radio science,2020,35(1):55-68. (in Chinese) doi: 10.13443/j.cjors.2019090206
[3] 宋凤丽,吴语茂. 基于自适应网格技术的快速物理光学方法[J]. 电波科学学报,2023,38(2):187-194. SONG F L, WU Y M. The fast physical optics method based on the adaptive mesh technique[J]. Chinese journal of radio science,2023,38(2):187-194. (in Chinese)
[4] 张慧雯,陈宇浩,黄晓伟,等. 基于多区域面积分方程的电大非均匀等离子体电磁特性高效计算方法[J]. 电波科学学报,2024,39(1):26-38. doi: 10.12265/j.cjors.2023313 ZHANG H W, CHEN Y H, HUANG X W, et al. A multi-region surface integral equation method for electromagnetic analysis of electrically large inhomogeneous plasma sheath[J]. Chinese journal of radio science,2024,39(1):26-38. (in Chinese) doi: 10.12265/j.cjors.2023313
[5] DOLEAN V,FAHS H,FEZOUI L,et al. Locally implicit discontinuous Galerkin method for time domain electromagnetics[J]. Journal of computational physics,2009,229(2):512-526.
[6] HOLLAND R,SIMPSON L. Finite-difference analysis of EMP coupling to thin struts and wires[J]. IEEE transactions on electromagnetic compatibility,1981(2):88-97.
[7] GUIFFAUT C,REINEIX A,PECQUEUX B. New oblique thin wire formalism in the FDTD method with multiwire junctions[J]. IEEE transactions on antennas and propagation,2011,60(3):1458-1466.
[8] LI P,SHI Y,JIANG L J,et al. Transient analysis of lumped circuit networks-loaded thin wires by DGTD method[J]. IEEE transactions on antennas and propagation,2016,64(6):2358-2369. doi: 10.1109/TAP.2016.2543803
[9] EDELVIK F. A new technique for accurate and stable modeling of arbitrarily oriented thin wires in the FDTD method[J]. IEEE transactions on electromagnetic compatibility,2003,45(2):416-423. doi: 10.1109/TEMC.2003.811294
[10] EDELVIK F,LEDFELT G,LOTSTEDT P,et al. An unconditionally stable subcell model for arbitrarily oriented thin wires in the FETD method[J]. IEEE transactions on antennas & propagation,2003,51(8):1797-1805.
[11] YANG Q,WEI B,LI L,et al. Analysis of the shielding effectiveness of thin wire cages by a hybrid DGTD-FDTD method[J]. IEEE transactions on electromagnetic compatibility,2022,64(6):2198-2206. doi: 10.1109/TEMC.2022.3201385