Fast electromagnetic simulation of large-scale dielectric-metal periodic structures based on generalized PMCHWT-EFIE and sub-entire-domain basis function method
-
摘要:
在计算大规模介质-金属复合周期结构的电磁散射时，传统积分方程方法存在未知量大、存储占用多和计算时间长等问题。本文采用广义Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT)-电场积分方程(electric field integral equation, EFIE)方法计算均匀介质金属复合结构的电磁响应。该方法通过在分界面处设置区域连接模型(contact-region modeling, CRM)来保证边界处的连续性。为加速子阵列阻抗矩阵填充，本文采用快速偶极子方法(fast dipole method, FDM)来提高计算效率并降低内存占用。结合子全域(sub-entire-domain, SED)基函数方法，子阵列的电流分布特征可被推广到大规模介质金属复合周期结构的电磁场计算中。数值算例表明，该方法能够在保证计算精度的同时大幅度降低计算代价，内存占用降低至商业软件Altair FEKO(使用快速多层多级子方法)的1/10，计算误差在2.6 dB以内。
-
关键词:
- 子全域 (SED) 基函数 /
- 区域连接模型 (CRM) /
- 快速偶极子方法 (FDM) /
- 大规模介质-金属复合周期结构
Abstract:Traditional integral methods may encounter some challenges, including massive unknowns, high storage requirements, and long computation times when calculating large-scale dielectric-metallic composite periodic structures. In this paper, the generalized Poggio-Miller-Chang-Harrington-Wu-Tsai(PMCHWT)-electric field integral equation (EFIE) method is used to calculate the dielectric-metallic structure with the contact-region modeling (CRM) technique to guarantee the continuity on the contact surface. To accelerate the filling of subarray impedance matrix, The fast dipole method (FDM) is introduced in this paper to reduce the computation time and storage requirements. Finally, the current distribution characteristics of subarrays are generalized to the composite dielectric-metallic periodic structures by the sub-entire-domain (SED) basis function method. Numerical results show that this method can significantly reduce computation time and memory usage while ensuring calculation accuracy. Compared with multilevel fast multipole method (MLFMM) of commercial software Altair FEKO, the computational memory is reduced by 1/10 and the calculation error is within 2.6 dB.
-
0 引 言
介质-金属复合周期结构如相控阵天线、频率选择表面被广泛应用于在雷达隐身、波束赋形等领域，因此周期结构的电磁特性研究极为重要。常用的电磁计算方法包括基于Floquet原理[1]的谱域法和全波方法。但在实际应用中，周期结构存在边界截断效应[2]，不满足周期边界条件，谱域法难以得到真实准确的计算结果。矩量法(method of moments, MoM)[3]、时域有限差分(finite-difference time-domain, FDTD)[4]等传统全波方法可以精确处理单元间耦合效应和周期结构边缘效应，但难以处理大规模介质金属复合结构带来的庞大未知量问题。因此，亟需研究面向大规模介质-金属复合周期结构的MoM[5-6]加速技术。
宏基函数方法是用于计算大规模有限周期结构的全波方法加速技术，包括特征基函数方法[7-8]、综合基函数方法[9]和子全域(sub-entire-domain, SED)基函数方法[10-12]等。宏基函数方法首先将待求解区域分解为数个子区域，分析子区域中阵列的电流分布特征，再利用周期性将其拓展至整个周期阵列。其中，SED基函数方法利用周期特性构造缩减矩阵，无需计算每一对单元间的互阻抗矩阵，极大地减少计算未知量，同时保留了对截断效应、耦合效应分析的精确度。
使用SED基函数方法计算介质-金属复合周期结构电磁响应的难点在于如何处理单元之间耦合以及介质-金属分界面的边界连续性。为解决这一问题，文献[13]使用基于SED基函数方法的介质表面积分方程(surface integral equation, SIE) (Poggio-Miller-Chang-Harrington-Wu-Tsai, PMCHWT)分析了耦合效应较弱情形下单元不相连时的石墨烯贴片阵列，可以计算介质-金属相连结构，但单元之间有间距。文献[14]使用多重平面波入射(multiple plane waves, MPWs)技术和奇异值分解(singular value decomposition, SVD)技术对基于SIE的SED基函数方法进行改进，在边界处采用半RWG基函数处理网格，使其能够分析金属周期单元连接情形下的周期阵列。文献[15]引入间断伽辽金(discontinuous Galerkin, DG)边界条件，使用快速偶极子方法(fast dipole method, FDM)和SED基函数方法加速体-面积分方程(volume surface integral equation, VSIE)的MoM方程求解过程，并分析了单元相邻情形下的介质-金属复合周期结构的辐射特性。
针对单元之间接触面以及介质-金属分界面的边界连续性问题，本文基于采用区域连接模型(contact-region modeling, CRM)[16-17]的广义PMCHWT-电场积分方程(electric field integral equation, EFIE)和SED基函数方法，提出了CRM-SIE-SED方法来计算介质-金属复合周期结构的电磁响应。首先，本文推导了基于广义PMCHWT-EFIE的SED基函数方程，并使用CRM方法处理介质-介质、介质-金属之间的边界连续性问题；然后采用FDM[18-19]提高阻抗矩阵填充效率，推导了基于广义PMCHWT-EFIE的FDM，给出广义PMCHWT-EFIE的聚集函数、转移函数以及发散函数的表达式；并针对单元间耦合效应较强的特点，使用多重平面波技术拓展SED基函数的解空间，并采用SVD技术筛选出SED基函数中占主导作用的特征向量。算例表明该方法具有较高的计算精度与计算效率。
1 基于广义PMCHWT-EFIE的CRM-SIE-SED算法
1.1 广义PMCHWT-EFIE
在CRM-SIE-SED方法中，为保证介质-金属分界面以及单元接触面上电流的连续性，采用基于CRM技术的广义PMCHWT-EFIE来计算介质-金属复合体。图1所示为介质-金属连接结构。设自由空间有一均匀介质体，介电常数为ε1，磁导率为μ1，介质体表面的法向向量为{{{\boldsymbol{n}}}_{\text{D}}}。自由空间的介电常数为{\varepsilon _0}，磁导率为{\mu _0}。入射电场为{{{\boldsymbol{E}}}_{{\text{inc}}}}，磁场为{{{\boldsymbol{H}}}_{{\text{inc}}}}。金属表面的法向向量为{{{\boldsymbol{n}}}_{\text{M}}}。入射波照射在介质-金属连接结构后，在自由空间产生电场{{{\boldsymbol{E}}}_1}和磁场{{{\boldsymbol{H}}}_1}，在介质体内产生电场{{{\boldsymbol{E}}}_2}和磁场{{{\boldsymbol{H}}}_2}。设介质体表面的等效面电流为{{{\boldsymbol{J}}}_{{\text{DS}}}}，等效面磁流为{{{\boldsymbol{M}}}_{{\text{DS}}}}，金属表面电流为{{{\boldsymbol{J}}}_{{\text{MS}}}}。
假设介质体与金属相隔距离\delta 的情形下，定义广义PMCHWT-EFIE[18]的算子L和算子K为
\begin{split} &{L}_{t}{{\boldsymbol{X}}}\left({\boldsymbol{r}}\right)=-\mathrm{i}{k}_{t}{\displaystyle {\iint }_{{S}^{\prime }}\left\{{\boldsymbol{X}}\left({{\boldsymbol{r}} '}\right)+\frac{1}{{k}_{t}^{2}}\left[\nabla \cdot {\boldsymbol{X}}\left({{\boldsymbol{r}} '}\right)\right]\nabla \right\}\frac{{\mathrm{e}}^{\mathrm{i}{k}_{t}R}}{4\text{π}R}\text{d}{S}^{\prime }}\\ &{K}_{t}^{\pm \text{，}0}{\boldsymbol{X}}\left({\boldsymbol{r}}\right)=\pm \frac{1}{2}{\boldsymbol{X}}\left({\boldsymbol{r}}\right)\times {{\boldsymbol{n}}}_{t}+\text{P}\text{.V}\text{.}{\displaystyle {\iint }_{{S}^{\prime }}\nabla \frac{{\mathrm{e}}^{\mathrm{i}{k}_{t}R}}{4\text{π}R}\times {\boldsymbol{X}}\left({{\boldsymbol{r}} '}\right)\text{d}{S}^{\prime }} \end{split} (1) 式中：{k_t}为媒质{{t}}中的波数，t = 0表示媒质为自由空间，t = 1表示媒质{{t}}为介质体； {{\boldsymbol{X}}} 为表面电流{{{\boldsymbol{J}}}_{{\text{DS}}}}、{{{\boldsymbol{J}}}_{{\text{MS}}}}或表面磁流{{{\boldsymbol{M}}}_{{\text{DS}}}}；{{\boldsymbol{r}}}为场点，{{\boldsymbol{r}}'}为源点；R为源点与场点的距离；{{\boldsymbol{n}}}为场点所在表面的法向向量； {\text{P}}{\text{.V}}{\text{.}} 为柯西主值积分。算子K的留数项的正负取决于场点位置[16]2。在计算形如\left\langle {{{{\boldsymbol{X}}}_w},K_t^{ \pm ,0}{{{\boldsymbol{X}}}_t}} \right\rangle 的阻抗元素时，当{{{\boldsymbol{X}}}_w}所在表面与{{{\boldsymbol{X}}}_t}所在表面不接触时，算子K的留数项取0。当{{{\boldsymbol{X}}}_w}所在表面与{{{\boldsymbol{X}}}_t}所在表面接触且{{{\boldsymbol{X}}}_w}所在表面指向{{{\boldsymbol{X}}}_t}所在表面与{{{\boldsymbol{X}}}_t}所在表面法线方向相反时，算子K的留数项取正值。当{{{\boldsymbol{X}}}_w}所在表面与{{{\boldsymbol{X}}}_t}所在表面接触且{{{\boldsymbol{X}}}_w}所在表面指向{{{\boldsymbol{X}}}_t}所在表面与{{{\boldsymbol{X}}}_t}所在表面法线方向相同时，算子K的留数项取负值。
介质-金属相连结构的表面积分方程[17]3为：
\begin{split} & {{{\boldsymbol{n}}}_{\text{D}}} \times \left[ \begin{gathered} {\eta _0}{L_0}\left( {{{{\boldsymbol{J}}}_{{\text{DS}}}}} \right) + {\eta _1}{L_1}\left( {{{{\boldsymbol{J}}}_{{\text{DS}}}}} \right) \\ + K_0^ - \left( {{{{\boldsymbol{M}}}_{{\text{DS}}}}} \right) + K_1^ + \left( {{{{\boldsymbol{M}}}_{{\text{DS}}}}} \right) \\ + {\eta _0}{L_0}\left( {{{{\boldsymbol{J}}}_{{\text{MS}}}}} \right) \\ \end{gathered} \right] = - {{{\boldsymbol{n}}}_{\text{D}}} \times {{{\boldsymbol{E}}}_{{\text{inc}}}} \\ &{{{\boldsymbol{n}}}_{\text{D}}} \times \left[ \begin{gathered} - K_0^ - \left( {{{{\boldsymbol{J}}}_{{\text{DS}}}}} \right) - K_1^ + \left( {{{{\boldsymbol{J}}}_{{\text{DS}}}}} \right) \\ + \frac{{{L_0}\left( {{{{\boldsymbol{M}}}_{{\text{DS}}}}} \right)}}{{{\eta _0}}} + \frac{{{L_1}\left( {{{{\boldsymbol{M}}}_{{\text{DS}}}}} \right)}}{{{\eta _1}}} \\ - K_0^ + \left( {{{{\boldsymbol{J}}}_{{\text{MS}}}}} \right) \\ \end{gathered} \right] = - {{{\boldsymbol{n}}}_{\text{D}}} \times {{{\boldsymbol{H}}}_{{\text{inc}}}} \\ & {{{\boldsymbol{n}}}_{\text{M}}} \times \left[ \begin{gathered} {\eta _0}{L_0}\left( {{{{\boldsymbol{J}}}_{{\text{DS}}}}} \right) + K_0^ + \left( {{{{\boldsymbol{M}}}_{{\text{DS}}}}} \right) \\ + {\eta _0}{L_0}\left( {{{{\boldsymbol{J}}}_{{\text{MS}}}}} \right) \\ \end{gathered} \right] = - {{{\boldsymbol{n}}}_{\text{M}}} \times {{{\boldsymbol{E}}}_{{\text{inc}}}} \end{split} (2) 广义PMCHWT-EFIE的连续性条件已包含在CRM技术中[16]3，不需要引入半RWG基函数。用RWG基函数对方程组(2)左右两边做内积，可以得到矩阵方程：
{{\boldsymbol{Z}}} \cdot {{\boldsymbol{I}}}{\text{ = }}{{\boldsymbol{V}}} (3) 其中，阻抗矩阵Z可以表示为子矩阵的形式：
\left[ {\begin{array}{*{20}{c}} {{\eta _0}{\boldsymbol{P}}_0^E + {\eta _1}{\boldsymbol{P}}_1^E}&{{\eta _0}\left( {{\boldsymbol{Q}}_0^{E - } + {\boldsymbol{Q}}_1^{E + }} \right)}&{{\eta _0}{\boldsymbol{P}}_0^E} \\ { - {\eta _0}\left( {{\boldsymbol{Q}}_0^{H - } + {\boldsymbol{Q}}_1^{H + }} \right)}&{{\eta _0}{\boldsymbol{P}}_0^H + \frac{{\eta _0^2}}{{{\eta _1}}}{\boldsymbol{P}}_1^H}&{ - {\eta _0}{\boldsymbol{P}}_0^{H + }} \\ {{\eta _0}{\boldsymbol{P}}_0^E}&{{\eta _0}{\boldsymbol{Q}}_0^{E + }}&{{\eta _0}{\boldsymbol{P}}_0^E} \end{array}} \right] (4) 其中，
\begin{split} &{\boldsymbol{P}}_{tmn}^E = \iint_S {{{\boldsymbol{f}}_m}} \cdot {L_t}\left( {{{\boldsymbol{f}}_n}} \right)dS = {\boldsymbol{Q}}_{tmn}^H \\ &{\boldsymbol{Q}}_{tmn}^{E \pm ,0} = \iint_S {{{\boldsymbol{f}}_m}} \cdot K_t^{ \pm ,0}\left( {{{\boldsymbol{f}}_n}} \right)dS = {\boldsymbol{P}}_{tmn}^{H \pm ,0} \end{split} (5) 式中：下标{\text{m}}、n分别代表场点{\text{m}}以及源点n；{{\boldsymbol{f}}}为RWG基函数。
激励向量计算表达式为：
\left\{ \begin{gathered} {{\boldsymbol{V}}}_{{\text{DS}}}^E = - \iint_S {{\boldsymbol{f}}} \cdot {{{\boldsymbol{E}}}_{{\text{inc}}}}{\text{d}}S \\ {{\boldsymbol{V}}}_{{\text{DS}}}^H = - {\eta _0}\iint_S {{\boldsymbol{f}}} \cdot {{{\boldsymbol{H}}}_{{\text{inc}}}}{\text{d}}S \\ {{\boldsymbol{V}}}_{{\text{MS}}}^E = - \iint_S {{\boldsymbol{f}}} \cdot {{{\boldsymbol{E}}}_{{\text{inc}}}}{\text{d}}S \\ \end{gathered} \right. (6) 1.2 SED方法原理
SED方法通过子阵列的电流性质来模拟阵列中所有单元的电流分布特征。如图2所示，假设整个阵列共有N \times N个单元。原问题分解为两个子问题：求解 3 \times 3 子阵列的RWG基函数展开系数和求解N \times N缩减矩阵展开系数。根据位置不同，所有单元的局域基函数都可以表示为 3 \times 3 子阵列某个单元的SED基函数。再用各个单元的局域基函数来构建各个单元间互阻抗矩阵的缩减矩阵。相比于其他宏基函数方法，SED基函数方法构造局域基函数的过程更加节省计算时间。
![]()
图 2 3 \times 3 子阵列及其9种SED基函数位置分布Fig. 2 The schematic view of the 3 \times 3 subarray and distributions of 9 SED basis functions将每个单元中所有面电流的贡献看成一个基函数，即单元p的SED基函数如下[15]：
{{\boldsymbol{I}}}_p^{{SED} } = \sum\limits_{n = 1}^{{N_a}} {{I_{p,n}}} {f_{p,n}}\left( {{\boldsymbol{r}}} \right) (7) 式中： {I_{p,n}} 为单元p的RWG基函数展开系数； {N_a} 为该单元的RWG基函数未知量。
基于SED基函数阻抗矩阵的缩减矩阵 {{\boldsymbol{Z}}}_{pq}^{{SED} } 、激励向量{{{\boldsymbol{V}}}^{{SED} }}与展开系数{{\boldsymbol{\alpha}} }[15]为：
\begin{gathered} {{\boldsymbol{Z}}}_{pq}^{\rm{SED} } = {\left[ {{{\boldsymbol{I}}}_p^{\rm{SED} }} \right]^{{\mathrm{T}}} } \cdot {{{\boldsymbol{Z}}}_{pq}} \cdot {{\boldsymbol{I}}}_q^{\rm{SED} } \\ {{\boldsymbol{V}}}_p^{\rm{SED} } = {\left[ {{{\boldsymbol{I}}}_p^{\rm{SED} }} \right]^{{\mathrm{T}}} } \cdot {{\boldsymbol{V}}} \\ {{\boldsymbol{\alpha}} = }{\left( {{{{\boldsymbol{Z}}}^{\rm{SED} }}} \right)^{ - 1}} \cdot {{{\boldsymbol{V}}}^{\rm{SED} }} \\ \end{gathered} (8) 式中：p与q为阵列单元的编号； {{{\boldsymbol{Z}}}_{pq}} 为场点在p单元、源点在q单元的互阻抗矩阵，维度为{{{N}}^a} \times {N^a}； {{\boldsymbol{Z}}}_{pq}^{{{\mathrm{SED}}} } 为1 \times 1的缩减矩阵； {{\boldsymbol{V}}}_p^{{{\mathrm{SED}}} } 为p单元的激励函数，维度为1 \times 1。求解得到展开系数{{\boldsymbol{\alpha}} }后，可以得到整个阵列电流系数的表达式：
{{{\boldsymbol{I}}}_n} = \sum\limits_{p = 1}^N {{\alpha _{p,k}}I_{p,k}^{{SED} }} (9) 1.3 MPWs法构建SED基函数
计算单元相连的周期结构时，仅仅一维的SED基函数并不能反映来自非相邻单元的耦合影响。为了更精准地计算周期阵列的散射性质，引入特征基函数中的多重平面波激励方法去扩展SED基函数的维度，将子阵列的SED基函数从一维扩展到 {N_{{{\mathrm{MPWs}}} }} 维。对于 3 \times 3 子阵列，采用不同极化方向不同入射角度的平面波入射于该子阵列，子阵列RWG基函数总个数记为{N_{{\text{sub}}}}，激励函数和电流系数扩展到 {N_{{{\mathrm{MPWs}}} }} \times {N_{{\text{sub}}}} 的矩阵，其矩阵方程表达式如下：
{{{\boldsymbol{Z}}}_{{N_{{\text{sub}}}} \times {N_{{\text{sub}}}}}} \cdot {{{\boldsymbol{I}}}_{{N_{{\text{sub}}}} \times N}} = {{{\boldsymbol{V}}}_{{N_{{\text{sub}}}} \times N}} (10) 虽然MPWs大大增加电流系数的解空间，但得到的SED基函数有自由度冗余。为加快求解SED基函数的缩减矩阵，引入SVD技术，将所有电流系数表示为一系列不相干的特征向量的线性加和：
{{{\boldsymbol{I}}}_{{N_{{\text{sub}}}} \times w}} = {{\boldsymbol{U\varSigma V}}} (11) 式中： {\varSigma } 为对角矩阵，有w个特征值，从上到下依次为 {\sigma _1},{\sigma _2}, \cdots {\sigma _w} ，并且 {\sigma _1} > {\sigma _2} > \cdots > {\sigma _w} 。矩阵 {U} 中第 d 个列向量 {u_d} 对应矩阵 {\varSigma } 中第 d 个特征值 {\sigma _d} 。选取前 \tau 个特征值，使得 {\sigma _\tau }/{\sigma _1} > 1\text{%} 并且 {\sigma _{\tau + 1}}/{\sigma _1} < 1\text{%} 。最终的电流 {{{\boldsymbol{I}}}^{{{\mathrm{SVD}}} }} 可以表示为矩阵 {U} 的列向量线性组合：
{{{\boldsymbol{I}}}^{{{\mathrm{SVD}}} }} = \left[ {\begin{array}{*{20}{c}} {{u_1}}&{{u_2}}& \cdots &{{u_\tau }} \end{array}} \right] (12) 类似于公式(8)，拓展后的缩减矩阵 {{\boldsymbol{Z}}}_{pq}^{{{\mathrm{SVD}}} } 、激励向量 {{\boldsymbol{V}}}_p^{{{\mathrm{SVD}}} } 与展开系数 {{{\boldsymbol{\alpha}} }^{{{\mathrm{SVD}}} }} 表达式为
\begin{gathered} {{\boldsymbol{Z}}}_{pq}^{\rm{SVD} } = {\left[ {{{{\boldsymbol{I}}}^{\rm{SVD} }}} \right]^{T} } \cdot {{{\boldsymbol{Z}}}_{pq}} \cdot {{{\boldsymbol{I}}}^{\rm{SVD} }} \\ {{\boldsymbol{V}}}_p^{\rm{SVD} } = {\left[ {{{{\boldsymbol{I}}}^{\rm{SVD} }}} \right]^{T} } \cdot {{\boldsymbol{V}}} \\ {{{\boldsymbol{\alpha}} }^{\rm{SVD} }}{ = }{\left( {{{{\boldsymbol{Z}}}^{\rm{SVD} }}} \right)^{ - 1}} \cdot {{{\boldsymbol{V}}}^{\rm{SVD} }} \\ \end{gathered} (13) 式中： {{\boldsymbol{Z}}}_{pq}^{{{\mathrm{SVD}}} } 为w \times w阶的p单元与q单元的缩减矩阵； {{\boldsymbol{V}}}_p^{{{\mathrm{SVD}}} } 为w \times 1阶的激励函数； {{{\boldsymbol{\alpha}} }^{{{\mathrm{SVD}}} }} 为w \times 1阶的展开系数。
经过MPWs扩充后的SED基函数能够处理单元间耦合效应较强的情形，得到精确的阵列散射性质，同时SVD只保留了部分贡献较大的基函数，节省了运算时间。
1.4 CRM-SIE-SED算法流程
CRM-SIE-SED算法流程如下：首先计算 3 \times 3 子阵列在多平面波入射下得到的扩展的子阵列电流系数 {{{\boldsymbol{I}}}_{{N_{{\text{sub}}}} \times N}} ，在有表面相互接触时采用CRM方法处理算子K的留数项；然后使用SVD技术挑选出对于阵列散射特性影响较大的基函数向量组 {{{\boldsymbol{I}}}^{{{\mathrm{SVD}}} }} ，计算 N \times N 阵列各单元的互阻抗矩阵，采用式(13)计算缩减矩阵 {{\boldsymbol{Z}}}_{pq}^{{{\mathrm{SVD}}} } 、激励向量 {{\boldsymbol{V}}}_p^{{{\mathrm{SVD}}} } 与展开系数 {{{\boldsymbol{\alpha}} }^{{{\mathrm{SVD}}} }} ；最后由式(9)计算整个阵列的电流系数 {{{\boldsymbol{I}}}_n} 。
2 FDM加速矩阵填充
由于SED基函数方法生成的缩减矩阵规模小，该算法主要计算时间是在阻抗矩阵填充上。采用FDM加速 3 \times 3 子阵列阻抗矩阵以及单元间的互阻抗矩阵的填充过程。对于边界间距大于0.15\lambda 的SED基函数，RWG基函数等效于偶极子，点v处RWG基函数的偶极矩表达式[20]为
{{{\boldsymbol{m}}}_v} = \int_{T_{\text{v}}^ + + T_v^ - } {{{{\boldsymbol{f}}}_v}} {\text{d}}S = {l_v}\left( {{{\boldsymbol{r}}}_{\text{v}}^ + - {{\boldsymbol{r}}}_v^ - } \right) (14) 式中： {{{l}}_v} 为点v处RWG基函数的公共边边长； {{\boldsymbol{r}}}_{\text{v}}^ + 、 {{\boldsymbol{r}}}_{\text{v}}^ - 分别为点v处RWG基函数正负三角形的几何中心。
自由空间中任意一点{{{\boldsymbol{r}}}_n}的电偶极子在{{{\boldsymbol{r}}}_m}处的远辐射场可以表示为[16]：
\begin{split} & {{{ {\boldsymbol{H}}}}^{{\text{sca}}}}({{\boldsymbol{r}}}) = \frac{{ - {{\mathrm{i}}} {k_t}}}{{4{\text{π} }}}\left( {{{{\boldsymbol{m}}}_n} \times {{\boldsymbol{r}}}} \right)C{{e} ^{{{\mathrm{i}}} {k_t}R}} \\ & {{{\boldsymbol{E}}}^{{\text{sca}}}}({{\boldsymbol{r}}}) = \frac{{\eta {{e} ^{{{\mathrm{i}}} {k_t}R}}}}{{4{\text{π} }}}\left[ \begin{gathered} \left( {{{{\boldsymbol{M}}}_n} - {{{\boldsymbol{m}}}_n}} \right)\left( {C - \frac{{{i} {k_t}}}{R}} \right) \\ + 2{{{\boldsymbol{M}}}_n}C \\ \end{gathered} \right] \end{split} (15) 其中，
\begin{split} &R = |{{\boldsymbol{r}}|} = \left| {{{{\boldsymbol{r}}}_m} - {{{\boldsymbol{r}}}_n}} \right| \\ &C = \frac{1}{{{R^2}}}\left( {1 + \frac{{{\mathrm{i}}} }{{kR}}} \right) \\ &{M_n} = \left( {{\hat {\boldsymbol{r}}} \cdot {{{\boldsymbol{m}}}_n}} \right){\hat {\boldsymbol{r}}} \end{split} (16) 阻抗矩阵可以表示为电场与磁场的表达式：
\begin{gathered} {{\boldsymbol{P}}}_{tmn}^E = - {{{\boldsymbol{m}}}_m} \cdot {{{\boldsymbol{E}}}^{{\text{sca}}}}\left( {{{{\boldsymbol{r}}}_m} - {{{\boldsymbol{r}}}_n}} \right) \\ {{\boldsymbol{Q}}}_{tmn}^E = {{{\boldsymbol{m}}}_m} \cdot {{{\boldsymbol{H}}}^{{\text{sca}}}}\left( {{{{\boldsymbol{r}}}_m} - {{{\boldsymbol{r}}}_n}} \right) \\ {{\boldsymbol{Q}}}_{tmn}^{E \pm } = {{\boldsymbol{Q}}}_{tmn}^E \pm \frac{1}{2}{{{\boldsymbol{m}}}_m} \cdot \left( {{{{\boldsymbol{m}}}_n} \times {{{\boldsymbol{n}}}_n}} \right) \\ \end{gathered} (17) 将式(17)代入到式(15)中，可以得到：
\begin{gathered} {{\boldsymbol{P}}}_{tmn}^E = \frac{{\eta {{e} ^{{{\mathrm{i}}} {k_t}R}}}}{{4{\text{π} }}}\left[ \begin{gathered} {{{\boldsymbol{m}}}_m} \cdot {{{\boldsymbol{m}}}_n}\left( {C - \frac{{{{\mathrm{i}}} {k_t}}}{R}} \right) \\ - \left( {{{{\boldsymbol{m}}}_m} \cdot {\hat {\boldsymbol{r}}}} \right)\left( {{\hat {\boldsymbol{r}}} \cdot {{{\boldsymbol{m}}}_n}} \right) \times \left( {3C - \frac{{{{\mathrm{i}}} {k_t}}}{R}} \right) \\ \end{gathered} \right] \\ {{\boldsymbol{Q}}}_{tmn}^E = \frac{{ - {{\mathrm{i}}} {k_t}C{{{\mathrm{e}}} ^{{{\mathrm{i}}} {k_t}R}}}}{{4{\text{π} }}}{{{\boldsymbol{m}}}_m} \cdot \left( {{{{\boldsymbol{m}}}_n} \times {{\boldsymbol{r}}}} \right) \\ {{\boldsymbol{Q}}}_{tmn}^{E \pm } = {{\boldsymbol{Q}}}_{tmn}^E \pm \frac{1}{2}{{{\boldsymbol{m}}}_m} \cdot \left( {{{{\boldsymbol{m}}}_n} \times {{{\boldsymbol{n}}}_n}} \right) \\ \end{gathered} (18) FDM先将网格依据几何位置分组到数个正方体中，根据场点所在的正方体中心点与源点所在正方体中心点位置来区分远场区以及近场区。近场区仍然采用常规数值方法进行计算。远场区则类似于多层快速多极子方法(multilevel fast multipole method, MLFMM)，将远场区之间阻抗矩阵写为聚集-转移-发散的向量矩阵相乘形式[16]：
\begin{gathered} {{\boldsymbol{P}}}_{tmn}^E = {{{\boldsymbol{M}}}_m} \cdot {{{\boldsymbol{T}}}_E} \cdot {{{\boldsymbol{M}}}_m} \\ {{\boldsymbol{Q}}}_{tmn}^E = {{{\boldsymbol{M}}}_m} \cdot \left( {{{{\boldsymbol{T}}}_H} \times {{{\boldsymbol{M}}}_m}} \right) \\ {{\boldsymbol{Q}}}_{tmn}^{E \pm } = {{\boldsymbol{Q}}}_{tmn}^E \pm \frac{1}{2}{{{\boldsymbol{m}}}_m} \cdot \left( {{{{\boldsymbol{m}}}_n} \times {{{\boldsymbol{n}}}_n}} \right) \\ \end{gathered} (19) 其中，聚集函数与转移函数表示为：
\begin{gathered} {{{\boldsymbol{M}}}_{m,n}} = {{{\boldsymbol{m}}}_{m,n}}{{{\mathrm{e}}} ^{{{\mathrm{i}}} k{{h}_{m,n}}}} \\ {{{\boldsymbol{h}}}_m} = {{{\boldsymbol{\hat r}}}_{ji}} \cdot {{\boldsymbol{r}}_{mj}} + \frac{{r_{mj}^2 - {{\left( {{{{\boldsymbol{\hat r}}}_{ji}} \cdot {{\boldsymbol{r}}_{mj}}} \right)}^2}}}{{2{r_{ji}}}} \\ {{{\boldsymbol{h}}}_n} = {{{\boldsymbol{\hat r}}}_{ij}} \cdot {{\boldsymbol{r}}_{ni}} + \frac{{r_{ni}^2 - {{\left( {{{{\boldsymbol{\hat r}}}_{ij}} \cdot {{\boldsymbol{r}}_{ni}}} \right)}^2}}}{{2{r_{ji}}}} \\ \end{gathered} (20) 式中： {{\hat {\boldsymbol{r}}}_{ij}} 为远场组i中心与远场组j中心的单位位置矢量； {{\hat {\boldsymbol{r}}}_{mj}} 为场点到远场组j中心的位置矢量； {{\hat {\boldsymbol{r}}}_{ni}} 为源点到远场组i中心的位置矢量。
由于算子K在留数项不为0时不能直接表示为聚集函数、转移函数、发散函数相乘的形式，需要先计算留数项不为0时的 Q_{tmn}^E ，再加上留数项的偶极子表达式。转移函数表示为：
\begin{gathered} {{{\boldsymbol{T}}}_E} = \frac{{\eta {{\text{e}}^{{{\mathrm{i}}} k{r_{ji}}}}}}{{4{\text{π} }}}\left[ {\overline {{\boldsymbol{\bar I}}} \left( {C - \frac{{{{\mathrm{i}}} k}}{{{{{r}}_{ji}}}}} \right) - {{{\hat {\boldsymbol{r}}}}_{ji}}{{{\hat {\boldsymbol{r}}}}_{ji}}\left( {3C - \frac{{{{\mathrm{i}}} k}}{{{r_{ji}}}}} \right)} \right] \\ {{{\boldsymbol{T}}}_H} = \frac{{ - {{\mathrm{i}}} kC{{\text{e}}^{{{\mathrm{i}}} k{r_{ji}}}}}}{{4{\text{π} }}}{{\boldsymbol{r}}} \\ \end{gathered} (21) 式中： {{\hat {\boldsymbol{r}}}_{ji}}{{\hat {\boldsymbol{r}}}_{ji}} 为二阶并矢， \overline {{\boldsymbol{\bar I}}} 为三阶单位并矢。
聚集函数与转移函数是对称的，可以复用、存储；并且 FDM方法避免了格林函数中双重面积分计算，节省了阻抗矩阵的计算时间。同时，由于周期结构中各单元之间互阻抗矩阵的格林函数具有一定的周期性和对称性，SED方法并不需要计算每一对单元间的互阻抗，提升了缩减矩阵阻抗填充效率。
3 数值算例与分析
本节围绕大规模复合材质周期阵列电磁散射特性计算，展现了结合CRM技术的SIE-SED方法的高效性与准确性。利用MPWs方法拓展基函数的自由度时，采用LU分解法处理3 \times 3子阵列的阻抗矩阵。入射波在水平角以间隔10°采样，俯仰角为间隔15°采样，每个方向的入射波都有水平极化与垂直极化两种极化方式，总入射波的数量为468； FDM中近场组判断标准为立方体中心距离小于3倍的立方体边长。
计算平台为Intel-Core i7-13700H @ 2.4 GHz。对比算法为FEKO的全波算法，FEKO版本为2021-
2100 (x64) 。3.1 圆形微带贴片结构算例
首先选择圆形微带贴片结构来验证CRM技术的精确度。圆形微带贴片结构的衬底以及尺寸如图3所示。衬底相对介电常数为2，厚度为0.2\lambda 。工作频率为2.4 GHz。剖分面元数为575，未知量为1 609。入射波方向为-z方向，极化方向沿-x方向。
图4是 CRM-SIE方法与FEKO在yoz平面与xoz平面的RCS计算结果，可以看出，在yoz平面RCS均值误差为0.20 dB，xoz平面为0.23 dB，说明CRM-SIE与FEKO的计算精度是一致的。
3.2 领结形周期阵列算例
算例2为15 \times 15领结形周期阵列，其单元分布和单元尺寸参数如图5所示。衬底相对介电常数为2，衬底厚度为0.2\lambda 。入射波为频率3.5 GHz的正入射平面波。剖分面元数为95 040，未知量为277 875。
![]()
图 5 15 \times 15 领结形周期阵列分布与天线单元尺寸图Fig. 5 15 \times 15 bowtie antenna array and the geometry of the cell图6为利用CRM-SIE-SED和FEKO的MLFMM计算得到的15 \times 15领结形周期阵列RCS结果，可以看出，yoz平面RCS均值误差为1.07 dB，xoz平面为1.06 dB，两种计算结果基本吻合。
![]()
图 6 15 \times 15 领结形周期阵列RCSFig. 6 Bistatic RCS of the 15 \times 15 bowtie antenna array3.3 圆形周期阵列算例
算例3为20×20圆形周期阵列，其单元分布和单元尺寸参数参照图7。衬底相对介电常数为2，厚度为0.2\lambda 。入射波为频率2.4 GHz的正入射平面波。剖分面元数为230 000，未知量为643 600。CRM-SIE-SED与FEKO的MLFMM的RCS计算结果如图8所示，yoz平面RCS均值误差为2.52 dB，xoz平面为1.79 dB。在大部分散射角度CRM-SIE-SED的计算精度都与MLFMM吻合，仅在计算yoz平面\theta = 160{\text{°}} 和\theta = 200{\text{°}} 时有偏差。推测原因为FDM近似精度不够。
![]()
图 7 20 \times 20 圆形周期阵列分布与天线单元尺寸图Fig. 7 20 \times 20 circular patch antenna array and the geometry of the cell![]()
图 8 20 \times 20 圆形周期阵列RCSFig. 8 Bistatic RCS of the 20 \times 20 circular patch antenna array3.4 计算时间对比
表1为CRM-SIE-SED和MLFMM的内存占用和计算时间对比。
表 1 天线阵列计算时间与内存Tab. 1 CPU time and memory of antenna arrays模型 计算方法 内存占用/GB 计算时间/min 领结形周期阵列 MLFMM(FEKO) 16.6 20.4 CRM-SIE-SED 1.4 15.8 圆形周期阵列 MLFMM(FEKO) 58.8 97.4 CRM-SIE-SED 5.1 60.1 在领结形周期阵列算例中，CRM-SIE-SED内存占用仅为MLFMM的8.4%，原因是该算例中CRM-SIE-SED最大内存消耗在计算3 \times 3子阵列上，需要存储9个单元构成的阻抗矩阵。CRM-SIE-SED比MLFMM的计算速度快了1.29倍。
在20×20圆形周期阵列算例中，CRM-SIE-SED方法在计算缩减矩阵时的内存超过了计算3 \times 3子阵列时需要的内存，但本算例中CRM-SIE-SED比FEKO的MLFMM仍有明显优势，内存占用仅为MLFMM的8.7%，计算速度比MLFMM快了1.61倍。
4 结 论
为求解大规模介质-金属复合周期结构电磁散射问题，本文提出基于广义PMCHWT-EFIE的CRM-SIE-SED方法，有效地降低了介质-金属复合目标的网格数目。该方法无需处理复杂的边界条件，边界之间电流连续性自动得到满足。本文给出了广义PMCHWT-EFIE的偶极子近似表达式，采用FDM方法对广义PMCHWT-EFIE的矩阵填充进行加速，采用MPWs法拓展子域基函数的自由度，使其能够计算单元相连的情形。算例结果表明，CRM-SIE-SED方法具有计算精度高和效率高的特点，与商业软件FEKO的MLFMM方法比较，内存占用降低到1/10，计算速度快了1.3倍。该方法具备计算大规模介质-金属复合周期结构散射问题的能力。
-
![]()
图 2 3 \times 3 子阵列及其9种SED基函数位置分布
Fig. 2 The schematic view of the 3 \times 3 subarray and distributions of 9 SED basis functions
![]()
图 5 15 \times 15 领结形周期阵列分布与天线单元尺寸图
Fig. 5 15 \times 15 bowtie antenna array and the geometry of the cell
![]()
图 6 15 \times 15 领结形周期阵列RCS
Fig. 6 Bistatic RCS of the 15 \times 15 bowtie antenna array
![]()
图 7 20 \times 20 圆形周期阵列分布与天线单元尺寸图
Fig. 7 20 \times 20 circular patch antenna array and the geometry of the cell
![]()
图 8 20 \times 20 圆形周期阵列RCS
Fig. 8 Bistatic RCS of the 20 \times 20 circular patch antenna array
表 1 天线阵列计算时间与内存
Tab. 1 CPU time and memory of antenna arrays
模型 计算方法 内存占用/GB 计算时间/min 领结形周期阵列 MLFMM(FEKO) 16.6 20.4 CRM-SIE-SED 1.4 15.8 圆形周期阵列 MLFMM(FEKO) 58.8 97.4 CRM-SIE-SED 5.1 60.1 
下载: 导出CSV
-
[1] STEVANOVIC I,CRESPO-VALERO P,MOSIG J R. An integral-equation technique for solving thick irises in rectangular waveguides[J]. IEEE transactions on microwave theory and techniques,2006,54(1):189-197. doi: 10.1109/TMTT.2005.860305
[2] CWIK T,MITTRA R. The effects of the truncation and curvature of periodic surfaces:a strip grating[J]. IEEE transactions on antennas and propagation,1998,36(5):612-622.
[3] 吴安雯,吴语茂,杨杨,等. 矩量法-物理光学混合算法计算多尺度复合目标电磁散射场[J]. 电波科学学报,2019,34(1):83-90. WU A W. WU Y M,YANG Y,et al. The MoM-PO hybrid method for calculating the scattered field of multi-scale composite targets[J]. Chinese journal of radio science,2019,34(1):83-90. (in Chinese)
[4] 周永金. 三维人工表面等离子体结构及其FDTD分析的软件实现[D]. 南京:东南大学,2011. ZHOU Y J. Three-dimensional designer surface plasmons structures and the FDTD implementation[D]. Nanjing:Southeast University,2011. (in Chinese)
[5] BACZEWSKI A D,DAULT D L,SHANKER B. Accelerated cartesian expansions for the rapid solution of periodic multiscale problems[J]. IEEE transactions on antennas and propagation,2012,60(9):4281-4290. doi: 10.1109/TAP.2012.2207037
[6] WATANABE Y,IGARASHI H. Accelerated FDTD analysis of antennas loaded by electric circuits[J]. IEEE transactions on antennas and propagation,2012,60(2):958-963. doi: 10.1109/TAP.2011.2173148
[7] PRAKASH V,MITTRA R. Characteristic basis function method:a new technique for efficient solution of method of moments matrix equations[J]. Microwave and optical technology letters,2003,36(2):95-100. doi: 10.1002/mop.10685
[8] YEO J,KOKSOY S,MITTRA R,et al. Efficient generation of method of moments matrices using the characteristic function method[J]. IEEE transactions on antennas and propagation,2004,52(12):3405-3410. doi: 10.1109/TAP.2004.836418
[9] LADISLAU M,VALERIU L,VECCHI G. Analysis of large complex structures with the synthetic-functions approach antennas and propagation[J]. IEEE transactions on antennas and propagation,2007,55(9):2509-2521. doi: 10.1109/TAP.2007.904073
[10] CUI T J,LU W B,QIAN Z G,et al. Accurate analysis of large-scale periodic structures using an efficient sub-entire-domain basis function method[J]. IEEE transactions on antennas and propagation,2004,4(11):4459-4462.
[11] LU W B,CUI T J,ZHAO H. Acceleration of fast multipole method for large-scale periodic structures with finite sizes using sub-entire-domain basis functions[J]. IEEE transactions on antennas and propagation,2007,55(2):414-421. doi: 10.1109/TAP.2006.889805
[12] WANG Q Q,ZHU H B,CHEN R S,et al. Analysis of finite frequency selective surfaces backed by dielectric substrate using sub-entire-domain basis function method[J]. COMPEL the international journal for computation and mathematics in electrical and electronic engineering,2015,34(4):1144-1155. doi: 10.1108/COMPEL-08-2014-0205
[13] YANG W,DING Y,LU W B,et al. Analysis of finite periodic structures of graphene with dielectric substrate using EFIE-PMCHW-SED method[C]//IEEE International Conference on Computational Electromagnetics,2020:31-33.
[14] XIONG T,LU W B,YANG W,et al. A new sub-entire-domain basis function for the accurate analysis of large-scale finite-sized periodic structures with connected-cell[C]//The 12th European Conference on Antennas and Propagation,2018:1-3.
[15] XIANG W,YANG W,LU W B. Fast subentire-domain basis functions method for analysis of composite finite periodic structures with dielectric-conductor cells[J]. IEEE antennas and wireless propagation letters,2023,22(2):233-237. doi: 10.1109/LAWP.2022.3207508
[16] CHU Y H,WENG C W,ZHAO J S,et al. A surface integral equation formulation for low-frequency scattering from a composite object[J]. IEEE transactions on antennas and propagation,2003,51(10):2837-2844. doi: 10.1109/TAP.2003.817999
[17] LIANG H F,SHAO H R,HU J. Fast frequency and material parameters sweep for the calculation of array structures[J]. IEEE antennas and wireless propagation letters,2022,21(7):1328-1332. doi: 10.1109/LAWP.2022.3166958
[18] CHEN X L,GU C Q,NIU Z Y,et,al. A hybrid fast dipole method and adaptive modified characteristic basis function method for electromagnetic scattering from perfect electric conducting targets[J]. Journal of electromagnetic waves and applications,2012,60(2):1186-1191.
[19] 陈新蕾. 基于特征基函数的电磁散射高效算法研究[D]. 南京:南京航空航天大学,2014. CHEN X L. Characteristic basis function (CBF)-based efficient algorithms for electromagnetic scattering[D]. Nanjing:Nanjing University of Aeronautics and Astronautics,2014. (in Chinese)
[20] YUAN J,GU C Q,HAN G. Efficient generation of method of moments matrices using equivalent dipole-moment method[J]. IEEE antennas and wireless propagation letters,2009,8:716-719. doi: 10.1109/LAWP.2009.2024337
计量
- 文章访问数: 58
- HTML全文浏览量: 5
- PDF下载量: 15
