无条件稳定和有条件稳定混合的时域谱元电磁分析方法

技术领域

本发明属于有条件稳定和无条件稳定混合的时域方法，特别是一种针对多尺度电磁 问题的快速分析技术。

背景技术

1864年，麦克斯韦（Maxwell）在前人的研究基础上，提出了能够对宏观电磁场基 本归路进行概括的数学方程组——著名的麦克斯韦方程，由此奠定了电磁学理论研究的 基础。电磁场理论的研究渐渐渗透到地学、生命科学、医学、材料科学和信息科学等众 多技术科学领域，大大促进了科学技术的发展和人类生活的变化。

早期很长一段时间，电磁场理论的研究致力于得到一些问题的解析解，然而完全用 解析方法求解的问题是十分有限的，不能解决什么问题。于是，为了解决科学技术中的 电磁场问题，又发展了一些近似方法和数值方法。但是，限于当时的计算条件，无法充 分发挥这些方法的作用，令某些问题得不到实质性的解决。随着电子计算机技术的飞速 发展，以高性能计算机技术为手段，结合电磁场理论和计算数学提供的各种数值方法， 应运而生了一门交叉学科——计算电磁学。

使用计算电磁学对电磁现象进行分析时，首先根据分析对象的特点建立相应的电磁、 数学模型。然后，选择合适的算法并在计算机上实现。当前计算电磁学方法按解域分可 分为频域方法和时域方法。频域方法主要有：以电磁场问题的积分方程为基础的矩量法 （MOM）和基于变分原理的有限元法（FEM）等；时域方法主要有：时域有限差分法 （FDTD）、时域有限元法（FETD）、时域积分法（TDIE）和时域伪谱方法(PSTD)等。

时域谱元法（Joon-Ho Lee and Qing Huo Liu,“A3-D Spectral-Element Time-Domain  Method for Electromagnetic Simulation,”IEEE Transactions on Microwave Theory and  Techniques.,vol.55,no.5,pp.983-991,May2007）可认为是一种特殊的时域有限元法，由于 时域谱元法采用的差分方式为中心差分，系数矩阵只含有质量矩阵，又由于该方法中所 选用的基函数为正交基函数，所以系数矩阵是块对角的，矩阵求逆会变的很容易，与时 域有限元方法比，这将大大减少计算时间。时域谱元法对网格采用的离散方式为曲六面 体离散，这能很好地拟合各种复杂的电磁结构，又因为时域谱元法网格离散的尺寸可以 很大，与时域有限差分方法相比，这将大大减少计算的未知量。

不连续迦辽金方法（Discontinuous Galerkin,DG）的发展已有了显著的进展，这些方 法适用于具有复杂结构和不均匀煤质的大尺寸问题。这些方法很大程度上得益于80年代 上半时有里德和希尔提出的中子传输方程的求解。这些方法的最重要的特点是允许基函 数（因此，数值解）在不同单元的交界面上不连续。在每个单元上引入一套局部的基函 数而不是在整个计算区域，且不同类型的单元例如六面体，棱柱体或者四面体可以在模 型中共存，同时可以允许两边的方程使用不同的差分格式，这样计算起来更加灵活方便。

现有的分析多尺度电磁问题的时域谱元法主要存在以下两个问题：

(1)采用六面体跟四面体混合剖分，但是六面体部分采用时域谱元法，四面体部分 采用时域有限元方法，而时域有限元方法生成的质量矩阵为稀疏症，求逆花费大量时间 再加上时间步长设置的很小，所以求解速度较慢。

(2)如果不同区域均采用六面体网格剖分，使用时域谱元法的时候，为保证算法稳 定性，时间步长必须按照最小的网格尺寸设定，这样整体的时间迭代步数就会很多，导 致求解时间慢。

发明内容

本发明的目的在于提供一种无条件稳定和有条件稳定混合的时域谱元电磁分析方 法，从而快速分析复杂多尺度电磁问题。

实现本发明目的的技术解决方案为：一种无条件稳定和有条件稳定混合的时域谱元 电磁分析方法，步骤如下：

第一步，对所要分析的电磁问题进行几何建模，将整体模型采用曲六面体进行剖分， 剖分之后得到各个体单元的顶点编号、坐标以及体的编号；

第二步，采用比较边长的方法找出剖分得到的网格中尺寸小于设定值的六面体，将 其标记为小尺寸区域，剩余网格标为大尺寸区域；

第三步，将电场值定义在整体模型剖分后的网格中的每个点上，并用时域谱元法中 的GLL多项式作为矢量基函数对电场在XYZ三个方向进行展开，代入时域波动方程， 并采用伽辽金测试，即测试基函数与展开基函数相同，得到矩阵方程。

第四步，将方程中的时间项用时间差分展开，在标记出来的小尺寸的区域采用具有 无条件稳定的Newmark-β差分格式，其他区域采用有条件稳定的中心差分格式，在进行 时间迭代时整体时间步长按中心差分区域设定，根据总的时间步数每一步直接求解，最 终求得时域电场值。

步骤一中剖分采用的曲六面体单元边长为1/10λ，λ为电磁波波长。

步骤二中设定值的选取与剖分后得到网格的尺寸有关，既要保证小于此设定值的网 格的数目占整体网格数目的比例尽可能小，又要保证大尺寸区域中最小网格的尺寸跟整 个区域最小网格的尺寸之比尽可能大。

步骤三中使用的GLL基函数形式如下：

${\Phi }_{\mathrm{rst}}^{\xi }=\hat{\xi }{\phi }_r^{\left(N_{\xi }\right)}\left(\xi \right){\phi }_s^{\left(N_{\eta }\right)}\left(\eta \right){\phi }_t^{\left(N_{\zeta }\right)}\left(\zeta \right)$

${\Phi }_{\mathrm{rst}}^{\eta }=\hat{\eta }{\phi }_r^{\left(N_{\xi }\right)}\left(\xi \right){\phi }_s^{\left(N_{\eta }\right)}\left(\eta \right){\phi }_t^{\left(N_{\zeta }\right)}\left(\zeta \right)$

${\Phi }_{\mathrm{rst}}^{\zeta }=\hat{\zeta }{\phi }_r^{\left(N_{\xi }\right)}\left(\xi \right){\phi }_s^{\left(N_{\eta }\right)}\left(\eta \right){\phi }_t^{\left(N_{\zeta }\right)}\left(\zeta \right)$

其中，${\Phi }_j^{\left(N\right)}\left(\xi \right)=\frac{1}{N(N+1)L_N\left({\xi }_j\right)}\frac{(1-{\xi }^2)L_N^{\prime }\left(\xi \right)}{\xi -{\xi }_j},$j＝0,1,LN，L_N(ξ)是N阶Legendre 多项式，将ξ∈[-1,1]内的节点{ξ_j,j＝0,1,LN}作为GLL积分点，它们是方程式 $(1-{\xi }_j^2)L_N^{\prime }\left({\xi }_j\right)=0$的(N+1)个根；

将电场用基函数展开，代入到时域波动方程

采用伽辽金法测试，即测试函数与基函数相同，获得总系数矩阵，求解方程

$\left[T\right]\frac{d^{2}e}{{\mathrm{dt}}^2}+\left[S\right]e=0$

步骤四中，在小尺寸的区域采用Newmark-β差分格式后方程为：

([T]+Δt²β[S])eⁿ⁺¹＝(2[T]-Δt²(1-2β)[S])eⁿ-([T]+Δt²β[S])e^n-1

在大尺寸区域采用中心差分格式后方程变为：

[T]eⁿ⁺¹＝(2[T]-Δt²[S])eⁿ-[T]e^n-1

求解方程，在每一步时间迭代中，先求解中心差分区域处的电场，再求解Newmark-β差 分区域的电场，最终得到总的时域电场值。

本发明与现有技术相比，其显著优点：(1)将模型用曲六面体剖分，可以很好的拟合 复杂物体的外形。(2)整体的时间步长不再受最小剖分网格尺寸的限制，可以设置的较大， 大大加快求解速度。

附图说明

图1是介质圆环谐振腔的结构示意图。

图2是本发明方法跟传统的中心差分法算得的频谱对比图。

具体实施方式

本发明基于一种无条件稳定和有条件稳定混合的时域谱元电磁分析方法，步骤如下：

第一步，模型剖分，将模型统一采用曲六面体网格进行剖分。

第二步，自适应找出网格中尺寸小于设定值的六面体，将其标记出来。

第三步，将电场用基函数进行展开，代入矢量波动方程，并采用伽辽金测试。

第四步，求解矩阵方程，将时间项用时间差分展开，在小尺寸的区域采用具有无条 件稳定的Newmark-β差分格式，其他区域采用中心差分，先进行Newmark-β区域的求 解，再进行中心差分区域的求解。

下面结合附图对本发明做进一步说明。

第一步，对分析的复杂结构进行建模，得到所需的几何信息。然后将模型采用曲六 面体网格进行剖分，剖分尺寸六面体单元边长为1/10λ(λ为电磁波波长)。剖分之后得 到各个体单元的顶点编号和坐标以及体的编号等。

第二步，设定出一个阈值，自适应找出网格中尺寸小于设定值的六面体，将其标记 出来，记为小尺寸网格，剩余的记为大尺寸网格。

第三步，将电场在结点处用GLL基函数进行展开，

在1-D标准参考单元ξ∈[-1,1]中，我们定义N阶GLL(Gauss-Lobatto-Legendre，高 斯-洛巴托-勒让德)基函数为：

${\phi }_j^{\left(N\right)}\left(\xi \right)=\frac{-1}{N(N+1)L_N\left({\xi }_j\right)}\frac{(1-{\xi }^2){L_N}^{\prime }\left(\xi \right)}{(\xi -{\xi }_j)}$

其中，j＝0,1,LN，L_N(ξ)是N阶勒让德多项式，L_N'(ξ)是它的导数。将ξ∈[-1,1] 内的网格点{ξ_j,j＝0,1,LN}作为GLL积分点，它们是方程式的(N+1) 个根，基函数满足φ_j(ξ_i)＝δ_ij的特性。

求解矢量波动方程

$▿\times ▿\times \overrightarrow{E}=-\mathrm{\mu ϵ}\frac{{\partial }^2\overrightarrow{E}}{{\partial t}^2}$

采用伽辽金测试，即测试基函数与展开基函数相同，获得紧凑格式

$\left[S\right]e+\left[T\right]\frac{d^{2}e}{{\mathrm{dt}}^2}=0$

其中，${\left[S\right]}_{\mathrm{ij}}^e=\underset{V^e}{\int \int \int }▿\times {\stackrel{\Gamma }{N}}_u^e·▿\times {\stackrel{\Gamma }{N}}_j^e\mathrm{dxdydz}$

${\left[T\right]}_{\mathrm{ij}}^e=\mathrm{\mu ϵ}\underset{V^e}{\int \int \int }{\stackrel{\Gamma }{N}}_i^e·{\stackrel{\Gamma }{N}}_j^e\mathrm{dxdydz}$

在大尺寸的区域采用中心差分格式

$\left[S\right]e^n+\left[T\right]\frac{e^{n+1}-2e^n+e^{n-1}}{\Delta t^2}=0$

[T]eⁿ⁺¹＝(2[T]-Δt²[S])eⁿ-[T]e^n-1，

在小尺寸的区域采用Newmark-β差分格式

$\left[S\right](\beta e^{n+2}+(1-2\beta )e^n+\beta e^{n-1})+\left[T\right]\frac{e^{n+1}-2e^n+e^{n-1}}{\Delta t^2}=0$

([T]+Δt²β[S])eⁿ⁺¹＝(2[T]-Δt²(1-2β)[S])eⁿ-([T]+Δt²β[S])e^n-1

先进行Newmark-β区域的求解，再进行中心差分区域的求解，最终可得到各点的电 场值。

为了验证本发明方法的有效性，下面分析了一个介质谐振腔的典型算例。

如图1所示的谐振腔中放有一介电常数9.8的介质环，a1=207.25mm,a2=440.75 mm,b=242mm,c=43mm,r1=9.0mm,r2=10.0mm,h=14.0mm，整体采用中心差分和 采用本发明混合方法计算得到的频谱如图2所示，可见两者结果吻合的很好。计算耗时 按表1所示，分别采用32,16,8,4个进程运行程序，可以看出本发明方法比传统的方法计 算时间大大节省。

表1