技术领域
[0001] 本
发明属于三维时域
电磁学数值求解技术领域，针对外加电磁源项的问题，采用一种含源项的时域杂交间断伽辽金数值方法来计算电
磁场。
背景技术
[0002] 麦克斯韦方程组是描述宏观电磁现象普遍规律的基本方程，对各种电磁问题的分析，归根结底是求解不同边界条件和媒质分布情况下Maxwell方程组的解。随着计算机与计算数学的发展，计算电磁学数值方法得到了快速的发展，对应的电磁仿真在
微波、天线、异向介质新材料等多种领域有着广泛的应用。
[0003] 目前的电磁学数值方法主要分为频域和时域两类求解。由于频域计算每次只能获得单个频点下的稳态
电磁场分布，而时域能够包含丰富瞬态信息，经过简单的时频变换可以获得宽频带频域信息，因此采用时域方法能够更加快速、直观地反映出物体在复杂电磁环境下的电磁响应特性。
[0004] 在求解实际的电磁场问题时，需要引入激励源到时域计算电磁学数值方法中。激励源是一种定义在三维物体表面或者二维物体上的已知源，这种激励源可以是电磁场、
电压源、
电流源或者电荷源。因此，对于数值方法分析问题时，一个重要的任务是对激励源的模拟，即选择合适的入射
波形式以及用适当方法将源加入到数值方法
迭代中。现有时域方法有时域有限元方法、时域有限体积方法、时域有限差分法与时域间断伽辽金方法等，这些方法都有各自的优缺点。
[0005] 时域有限差分法采用规则网格剖分，该方法在模拟计算复杂
几何模型往往受到限制，且计算结果
精度有限。通过入射波作为初值条件加入到网格中所有
位置的
电场和磁场分量上，由这些初始场正负和幅度的分布来决定平面波的极化和传播方向。但是该方法加源的主要缺点是当脉冲宽度较大时，计算空间必须扩大到足以容纳完整的脉冲波形，从而造成计算内存不足。尽管时域间断伽辽金方法是目前新型的方法，能支持非结构非共形网格，灵活的时间格式与并行性等优点，但是该方法的单元交界面上的未知量是重复的，导致在获得相同精度时，所需的未知量明显多于经典有限元方法。
[0006] 目前在时域间断伽辽金方法中，还可以通过总场-散射场分离方式加入平面波。对于这种加源方式，平面波从靠近
相位参考点的连接边界出发，传播到内部障碍物的位置时，在与障碍物上的各种媒质发生电磁作用时，所产生的所有误差都将累积到
算法误差中，从而降低了电磁散射建模精度。对于时域间断伽辽金方法，网格
采样密度会因为较小的数值色散误差而比时域有限差分法而稀疏很多，从而导致应用总场-散射场技术时散射场区的
泄漏较大。可见这些常用的加源的方式存在着不同程度的缺点，特别是数值方法本身也有各自的利弊，对于复杂模型的优化仿真，已经不能满足设计者的要求，因此需要一种稳定可靠的数值方法来处理外加电磁源项，进而获得电磁问题的电磁响应特性。
发明内容
[0007] 针对上述存在问题或不足，本发明提供了一种三维时域杂交间断伽辽金方法外加电磁源项的处理方法。
[0008] 一种三维时域杂交间断伽辽金方法外加电磁源项的处理方法，包括以下步骤：
[0009] 步骤A、根据目标
电子器件的物理结构，结合工作环境与边界条件对其仿真建模；
[0010] 本发明基于有源三维时域Maxwell方程组来进行推导说明。首先给出有源三维时域Maxwell方程组如下所示：
[0011]
[0012] 其中，T表示计算电磁学时域中模型仿真计算的最终时间，E(x,y,z,t)和H(x,y,z,t)分别是电场强度和磁场强度矢量。计算区域Ω为三维，即 ε为介质相对
介电常数，μ为介质的相对磁导率。M为磁流密度，J为电流密度，电流密度J等于导体电流密度Jc与外加电流密度Ji之和，即
[0013] J＝Jc+Ji (2)
[0014] 式中Jc＝σeE，σe为电导率。
[0015] 同样地，磁流密度M等于导体磁流密度Mc与外加磁流密度Mi之和，即[0016] M＝Mc+Mi (3)
[0017] 式中Mc＝σmH，σm为磁导率。
[0018] 由于本发明只涉及外加源项，故重写(1)式，则有
[0019]
[0020] 步骤B、采用四面体单元剖分三维求解区域，且面离散和体积离散必须相容；
[0021] 采用四面体单元剖分三维求解区域是一种公知过程，因此本步骤不再详细描述。需要注意的是，面离散和体积离散必须相容。本发明将计算区域Ω划分成Nh个四面体网格的集合 其中每个体单元用τi(i＝1,2,3,…,Nh)表示，即
[0022]
[0023] 定义面集合 是由Nf个三
角形面单元Df组成，即
[0024]
[0025] 步骤C、先考虑只有外加电流密度Ji时，推导时域杂交间断伽辽金方法的半离散形式；
[0026] 由于只有外加电流密度，因此(4)式变为
[0027]
[0028] 对于(7)式。定义有限元函数空间Vh与有限元迹空间Mh后，对于 用基函数Φ将电场和磁场进行展开，得到计算区域解析值(E,H)的逼近解(Eh,Hh)，那么对一个体单元τi有
[0029]
[0030] 根据矢量格林定理，上式等效于：
[0031]
[0032] 其中 是时域间断伽辽金方法上面上的数值通量，n是计算区域边界 上的外法向单位矢量。对于时域杂交间断伽略金方法，用下面的数值迹去替换数值通量[0033]
[0034] 其中 是面上的切向电场与磁场。Λh是整个计算区域上面集合 上的杂交项。τ＞0是局部稳定系数。考虑吸收边界条件Γa的边界方程结合守恒条件：
[0035]
[0036] 通过推导，得到有源时域杂交间断伽辽金方法的半离散形式：
[0037]
[0038] 式中上标inc表示在吸收边界处的入射场值。(12)式的前两个式子中的电场和磁场用杂交量线性表出，是局部线性方程，结合第三个守恒式子，最终形成一个只与杂交量有关的全局线性系统，一旦求出杂交量，则根据局部线性系统求出电磁场值。以上的推导过程与无源时域杂交间断伽辽金方法的半离散形式类似，但是在外加电流密度后，下面给出一种源项处理技术：
[0040] E(r,t)＝E(r)f(t) (13)
[0041] 其中E(r)是空间分布函数，f(t)是时间函数。在总场格式下，将局部辐射源转换为体源来实现激励源的强加。回顾(7)式的第一个方程，在t＝0和t上进行积分可得[0042]
[0043] 为了获得形如(13)式的局部辐射源，通过引入如下体电流源实现：
[0044]
[0045] 这种方法易于实现，还可以通过在源点附近网格加密，以提高计算精度。
[0046] 注意到半离散形式(12)中，关于外加源项还需要进行一个体积分求解，即事实上，由于外加源项的形式较多，往往Ji具有复杂的表达式，无法正常像电场和磁场用基函数展开，因此如何求解空间离散的这个体积分显得尤为重要。针对源项的形式，给出两种处理办法：
[0047] 第一种是能施加形如(13)式的局部辐射源，这种情况下，体电流源用(15)式获得，由于电场能用基函数Φ展开，那么初始的Ji也可以用基函数展开，则在源项体积分求解时，基函数与基函数作用会产生一个关于源项的
质量矩阵 那么
[0048]
[0049] 将(16)式带入(12)式，只需要在时间迭代时增加源项的系数 的迭代。
[0050] 然而对于一些Ji无法用基函数展开表示，用上面的方法会存在一定的误差，因此，这里给出第二种处理办法，即应用高斯数值积分的办法来求解外加源项的体积分。这里给出高斯积分求解四面体体积分的公式为：
[0051]
[0052] 如果在积分公式中涉及到与坐标参量G(x,y,z)有关的量时，则有[0053]
[0054] 其中F是被积函数，(L1i,L2i,L3i,L4i)是被积函数的采样点，wi表示权重，V是体积，(L1,L2,L3,L4)表示四面体
节点基函数。注意当应用高斯数值积分处理时，其积分精度与被积函数的阶数有关，当被积函数的阶数越高，计算时所需要的采样点越多。
[0055] 步骤D、考虑只有外加磁流密度Mi时，推导时域杂交间断伽辽金方法的半离散形式；由于只有外加电流密度，因此(4)式变为
[0056]
[0057] 由于步骤C中已经推导了只有外加电流密度Ji的时域杂交间断伽辽金方法的半离散形式。同样的，对于外加磁流密度Mi也可以用步骤C的方法实现。这里我们从等效原理的角度实施另一种添加外加磁流密度Mi的处理方式。在加源时，等效面磁流与入射电场有着如下瞬时值关系
[0058] Mi＝-n×Einc (20)
[0059] 由于(20)式中，入射电场Einc可以用基函数Φ展开，则实施起来也很容易。将(20)式带入(19)式中，对于 那么对一个体单元τi有
[0060]
[0061] 通过推导得到有源时域杂交间断伽辽金方法的半离散形式：
[0062]
[0063] 对于(22)式中含源项的面积分 将Einc用基函数Φ展开后，基函数与基函数作用会产生一个面质量矩阵 那么
[0064]
[0065] 步骤E、基于步骤C与步骤D得到的半离散形式，推导含源项的局部线性系统和全局线性系统。
[0066] 考虑到杂交量只存在于面单元上，并且保持单值，因此将杂交量看成一个待求常量，只考虑电场和磁场的时间迭代。这里采用二阶隐式的Crank-Nicolson时间格式进行时间离散，由于该过程是一种公知的时间离散格式，故本步骤不再详细描述时间离散过程。下面我们主要讨论时间离散后得到的含源项的局部线性系统和全局线性系统的具体形式。
[0067] 首先考虑局部线性系统。在一个四面体单元上，用待求系数 表示tn+1时刻的电场 和磁场 的待求系数，同理 表示tn时刻的电磁场的待求系数。则根据有源时域杂交间断伽辽金方法的半离散形式(12)、(22)的前两个方程，推导出一个四面体上的局部线性系统为
[0068]
[0069] 其中Λe表示τi上的所有面单元的杂交量。 和 是通过基函数作用后能得到的局部矩阵。而 与Winc是与源项相关的矩阵与向量，是已知的。由于守恒条件并未改变，因此半离散格式(12)和(22)的第三个方程是相同的，因此结合该方程，得到全局线性系统在一个四面体单元内的基本方程
[0070]
[0071] 其中 是通过基函数作用能得到的局部矩阵，be是局部右端项。根据(25)式依次
叠加每个四面体单元后，最终得到一个全局线性系统，即
[0072]
[0073] 矩阵 就是全局线性矩阵，y是右端项向量。一旦通过全局线性系统(26)获得全局杂交量Λ，根据每个四面体上局部杂交量与全局杂交量的映射关系，就可以得到每个四面体上的Λe，从而由局部线性系统(24)式就可以求得每个单元内待求系数 至此，我们已经详述完有源三维杂交时域间断伽辽金数值方法的全部过程。
[0074] 综上所述，本发明提供的三维时域杂交间断伽辽金方法外加电磁源项的处理方法，基于有源的Maxwell方程，推导时域杂交间断伽辽金的构造形式。针对不同的源项进行不同的处理：使其针对不同的问题能够适用于不同的边界条件，例如吸收边界ABC和完美匹配层PML。相比现有的无源时域杂交间断伽辽金方法，这个改进后有源的时域杂交间断伽辽金方法适用性更强。最后在时间离散上，我们选择无条件稳定的隐式时间格式,相比显式时间格式可以适当扩大时间步长，缩减模型的总体仿真时间。相比目前流行的时域间断伽辽金方法，本发明方法能够具有更少的全局未知量以及显著的计算性能。
附图说明
[0076] 图2局部辐射源激励下的金属复合结构的
三维网格剖分图。
[0077] 图3有源三维杂交时域间断伽辽金数值方法仿真金属复合结构所得电场|Ex|分布图，(a)XOZ平面(b)YOZ平面。
[0078] 图4基于二阶基函数金属复合结构某个观察点随时间变化的电场Ex场值变化。
具体实施方式
[0079] 下面结合附图和
实施例进一步的详细说明本发明。
[0080] 参照图1，三维时域杂交间断伽略金方法外加电磁源项的处理方法，包括以下步骤：
[0081] 步骤A、根据目标电子器件的物理结构，结合工作环境与边界条件对其仿真建模；
[0082] 步骤B、采用四面体单元剖分三维求解区域，面离散和体积离散必须相容；
[0083] 步骤C、先考虑只有外加电流密度Ji时，推导时域杂交间断伽辽金方法的半离散形式；
[0084] 由于只有外加电流密度，因此(4)式变为
[0085]
[0086] 通过推导，得到有源时域杂交间断伽辽金方法的半离散形式：
[0087]
[0088] 式中上标inc表示在吸收边界处的入射场值。(12)式的前两个式子中的电场和磁场用杂交量线性表出，是局部线性方程，结合第三个守恒式子，最终形成一个只与杂交量有关的全局线性系统，一旦求出杂交量，则根据局部线性系统求出电磁场值。以上的推导过程与无源时域杂交间断伽辽金方法的半离散形式类似，但是在外加电流密度后，下面给出一种源项处理技术：
[0089] 对于局部辐射源，表示为如下形式
[0090] E(r,t)＝E(r)f(t) (13)
[0091] 其中E(r)是空间分布函数，f(t)是时间函数。在总场格式下，可以将局部辐射源转换为体源来实现激励源的强加。回顾(7)式的第一个方程，我们在t＝0和t上进行积分可得[0092]
[0093] 为了获得形如(13)式的局部辐射源，通过引入如下体电流源实现：
[0094]
[0095] 这种方法易于实现，但为了保证计算精度，还可以通过在源点附近网格加密。
[0096] 注意到半离散形式(12)中，关于外加源项我们还需要进行一个体积分求解，即事实上，由于外加源项的形式较多，往往Ji具有复杂的表达式，我们无法正常像电场和磁场用基函数展开，因此如何求解空间离散的这个体积分显得尤为重要。这里针对源项的形式，给出两种处理办法：
[0097] 第一种是能施加形如(13)式的局部辐射源，这种情况下，体电流源可以用(15)式获得，由于电场能用基函数Φ展开，那么初始的Ji也可以用基函数展开，则在源项体积分求解时，基函数与基函数作用会产生一个关于源项的质量矩阵 那么
[0098]
[0099] 将(16)式带入(12)式，只需要在时间迭代时增加源项的系数 的迭代。
[0100] 然而对于一些Ji无法用基函数展开表示，用上面的方法会存在一定的误差，因此，这里给出第二种处理办法，即应用高斯数值积分的办法来求解外加源项的体积分。这里我们给出高斯积分求解四面体体积分的公式为：
[0101]
[0102] 如果在积分公式中涉及到与坐标参量G(x,y,z)有关的量时，则有[0103]
[0104] 其中F是被积函数，(L1i,L2i,L3i,L4i)是被积函数的采样点，wi表示权重，V是体积，(L1,L2,L3,L4)表示四面体节点基函数。注意当应用高斯数值积分处理时，其积分精度与被积函数的阶数有关，当被积函数的阶数越高，计算时所需要的采样点越多。
[0105] 步骤D、考虑只有外加磁流密度Mi时，推导时域杂交间断伽辽金方法的半离散形式；由于只有外加电流密度，因此(4)式变为
[0106]
[0107] 由于步骤C中我们已经推导了只有外加电流密度Ji的时域杂交间断伽辽金方法的半离散形式。同样的，对于外加磁流密度Mi也可以用步骤C的方法实现。这里我们从等效原理的角度实施另一种添加外加磁流密度Mi的处理方式。在加源时，等效面磁流与入射电场有着如下瞬时值关系
[0108] Mi＝-n×Einc (20)
[0109] 将(20)式带入(19)式中，通过推导，得到有源时域杂交间断伽辽金方法的半离散形式：
[0110]
[0111] 对于(22)式中含源项的面积分 将Einc用基函数Φ展开后，基函数与基函数作用会产生一个面质量矩阵 那么
[0112]
[0113] 步骤E、基于步骤C与步骤D得到的半离散形式，推导含源项的局部线性系统和全局线性系统。
[0114] 考虑到杂交量只存在于面单元上，并且保持单值，因此将杂交量看成一个待求常量，只考虑电场和磁场的时间迭代。这里采用二阶隐式的Crank-Nicolson时间格式进行时间离散，由于该过程是一种公知的时间离散格式，故本步骤不再详细描述时间离散过程。下面我们主要讨论时间离散后得到的含源项的局部线性系统和全局线性系统的具体形式。
[0115] 首先考虑局部线性系统。在一个四面体单元上，用待求系数 表示tn+1时刻的电场 和磁场 的待求系数，同理 表示tn时刻的电磁场的待求系数。则根据有源时域杂交间断伽辽金方法的半离散形式(12)、(22)的前两个方程，可以推导出一个四面体上的局部线性系统为
[0116]
[0117] 其中Λe表示τi上的所有面单元的杂交量。 和 是通过基函数作用后inc能得到的局部矩阵。而 与W 是与源项相关的矩阵与向量，是已知的。由于守恒条件并未改变，因此半离散格式(12)和(22)的第三个方程是相同的，因此结合该方程，最终可以得到一个全局线性系统，即
[0118]
[0119] 矩阵 就是全局线性矩阵，y是右端项向量。一旦通过全局线性系统(26)获得全局杂交量Λ，根据每个四面体上局部杂交量与全局杂交量的映射关系，就可以得到每个四面体上的Λe，从而由局部线性系统(24)式就可以求得每个单元内待求系数 至此，我们已经详述完有源三维杂交时域间断伽辽金数值方法的全部过程。
[0120] 实施例：
[0121] 以下面的一个金属复合球为模型，添加点源来分析有源三维杂交时域间断伽辽金数值方法的局部源辐射问题。
[0122] 图2的模型是由一个完美导体金属球和完美导体金属圆柱体构成的复合结构。其中金属球半径为1m，圆柱体半径和高度分别为0.02m和0.1m,被放置在距离金属球表面0.1m的位置。在金属球与圆柱体间隙的中心rs＝(-1.05,0,0)，设置了一个
频率为f＝300MHz,正z方向的局部辐射源 为了在局部辐射源周围获得较高的计算精度，需要在源点rs周围进行局部网格加密，如图2所示，加密后的网格用阴影立方体(平面形状正方形区域)标记。计算区域为立方体，由吸收边界ABC截断。采用网格离散后，共有90063个四面体，19043个节点，188057个面单元，其中最大和最小网格棱长为3.90E-01m和3.46E-
03m，二者比值为112.7，是一个典型的多尺度问题。
[0123] 为了验证本发明中的有源三维杂交时域间断伽辽金数值方法对源项的处理，我们采用高斯数值积分的办法来求解外加源项的体积分，基于二阶基函数仿真得到电场Ex模值在X0Z平面的分布如图3：
[0124] 图3为有源三维杂交时域间断伽辽金数值方法仿真金属复合结构所得电场|Ex|分布图，(a)X0Z平面(b)Y0Z平面。从两个不同视角的电场分布图，可以看出金属复合结构的仿真结果能够基本符合电磁散射的规律。
[0125] 进一步地，我们还使用显式时域间断伽辽金方法(exDGTD)作为对比算法，计算了两种数值算法下某个观察点上随着时间变化的电场Ex场值变化，如图4：
[0126] 由图4，我们可以看出当采用有源三维杂交时域间断伽辽金数值方法时，基于二阶基函数imHDGTD计算的Ex与exDGTD的结果吻合，说明了结果的可靠性。最后我们进行性能分析，对于exDGTD，其全局自由量未知数为5,403,780，而imHDGTD的全局自由量未知数能够减少很多，仅为2,256,684。同时由于时间步长的增大，总体仿真时间能大量减少,进一步说明该方法能够大幅提高计算性能。
