技术领域
[0001] 本
发明属于时域计算电磁学领域，涉及一种基于混合网格及时间步长的电磁学多尺度计 算方法。
背景技术
[0002] 在现代化战场上，军用飞机，舰艇，战车等武装设备都携带有大量的
电子武器，这些武 器有的可以检测地方军情，有的可以定向发射高功率电子脉冲进行电子对抗作战。因此，在 一定的战场空间内存在着非常复杂的
电磁干扰以及高频脉冲等复杂电磁现象。信息化武装设 备对于复杂电磁环境是非常敏感的，为了保证这些设备可以在复杂的战场电磁环境下依旧发 挥其效用，需要在设计时能充分模拟战场的复杂多尺度电磁环境，这样才能为电磁武装设备 的分析和设计提供充分的理论
支撑。
[0003] 另外，随着电子制造工艺的进步，大规模集成
电路以及现代天线的结构越来越复杂，且 具有复杂多尺度结构电功能符合材料问世，如
光子晶体、超
磁性材料、电热材料和光操纵材 料等，对其进行电磁模拟的难度随之上升。为了很好地模拟这些复杂电路、新型天线以及复 杂材料组成的电子系统，需要对其所有微观结构的影响都加以考虑，因其结构复杂，尺寸相 差很大，对其进行研究会面临较大的困难。
[0004] 传统的数值方法在处理多尺度问题时遇到了较大的困难，例如上述复杂电磁环境下的电 磁干扰，新型纳米结构电磁材料性能分析等。这些电磁问题具有典型多尺度特征，有的具有 细节远小于
波长的电精细结构，而有的具有远大于波长的电大结构。对于这种几何尺度跨度 很大的问题，若采用传统的方法进行网格离散，为了使网格能精确计算求解域的几何、材料 特征，需要划分足够精细的网格才能得到满意的
精度。但是这样会导致网格数目非常巨大， 在建模过程中会产生较大规模的未知量，特别是时域问题，随着网格尺寸的减小，需要进行 的时间
迭代次数也随之增加，这会导致巨量计算资源的花费，有时甚至会导致无法计算。因 此，如何通过数值方法解决大规模多尺度时域电磁问题是目前计算电磁学研究的一大难点。
发明内容
[0005] 针对上述存在问题或不足，为解决现有数值方法求解大规模多尺度时域电磁问题时计 算量巨大的问题，本发明提供了一种基于混合网格及时间步长的电磁学多尺度计算方法。
[0006] 一种基于混合网格及时间步长的电磁学多尺度计算方法，包括以下步骤：
[0007] 步骤1、根据实际电磁问题建立计算模型；
[0008] 步骤2、对于步骤1所建计算模型，根据实际电磁问题的需要，忽略所有的微观细节部分， 在宏观尺度上进行网格剖分，将整个计算模型的求解域划分为K个互不重叠的宏观单元。宏 观单元的尺寸大于求解域内的最小特征尺寸。
[0009] 步骤3、考虑所有微观结构及微观信息，根据实际电磁问题的需要使用混合网格对宏观单 元进行第二次网格剖分，得到微观细网格。
[0010] 微观单元的剖分在宏观单元的离散
基础上进行，微观单元的大小不大于宏观单元内部微 观结构的特征尺度。对于宏观单元内部没有微观结构的部分，二维采用四边形网格剖分、三 维采用六面体网格；对于内部具有微观结构的部分，二维采用三
角形网格剖分、三维用四边 体网格。
[0011] 步骤4、通过伽辽金方法构造麦克斯韦方程组的杂交弱形式，并对电
磁场进行线性分解。
[0012] 步骤4-1、建立求解域的电磁问题。
[0013] 设求解域为Ω∈R3， 为求解域边界，其中ΓPEC为理想电壁边界，ΓABC为吸收边界， 则求解域内的电磁问题写为：
[0014]
[0015] 步骤4-2、通过伽辽金方法建立麦克斯韦方程组的杂交弱形式。
[0016] 求解域Ω进行宏观网格离散后区域为TH，区域内的离散后的单元K∈TH。假定Pm(K) 为单元内维数不小于m的多项式空间，用V来表示空间离散，则有 V:＝{v∈L2(Ω):v|K∈[Pm(K)]3for_all(K∈TH)}，麦克斯韦方程组的伽辽金弱形式可以表示 为：
[0017]
[0018] 离散后的所有宏观单元边界面集合记为εH，在单元边界面 存在切向分量函数 空间Λ，可以表示为 引入杂交变量λ∈Λ，并将式(2)变为杂交 形式：
[0019]
[0020] 步骤4-3、根据宏观问题与微观问题的需要对
电磁场进行线性分解
[0021] 将
电场与磁场进行线性分解，分解为两个部分：
[0022]
[0023] 其中， 与 分量只与该时刻杂交变量λ及宏观边界条件有关，与单元内部的电磁场 及
电流作用 无关，而 与Hf分量只与内部的电磁场及电流作用 有关，与杂交变量λ及 宏观边界条件无关。对电磁场进行线性分解的目的是将宏观尺度边界条件的作用与单元内部 电磁场的作用分离，从而得以构建全局宏观问题与微观局部问题，并且便于实现宏观问题与 微观问题之间的耦合。
[0024] 步骤5、引入多尺度基函数，并进行基函数离散，之后在宏观尺度进行分析，构造宏观“全 局问题”。
[0025] 步骤5-1、引入多尺度基函数，并进行基函数离散。
[0026] 宏观问题建立在宏观单元的交界面上，首先令ψ为宏观粗网格面上的基函数，并且令 dim(ΛH)表示宏观单元交界面总的
自由度维数。因为与宏观网格相关的量主要有λ、及 所以在宏观尺度上对杂交变量λ及电磁场 离散后可表示为：
[0027]
[0028] 其中αi为宏观尺度的待定系数，需要通过宏观问题求解得到，而为多 尺度基函数，需要通过数值方法在微观单元中进行求解才能得到其具体形式。
[0029] 步骤5-2、通过宏观尺度进行分析，构造宏观“全局问题”。
[0030] 宏观问题的本质目的是建立整个求解域的连接性，使每一个宏观单元可以彼此耦合，若 区域内不存在
能量损耗(如只有PEC边界)，则全局问题存在以下形式：
[0031]
[0032] 若宏观区域存在吸收边界及激励，则需要将边界信息引入全局问题，引入吸收边界及激 励后全局问题为：
[0033]
[0034] 在宏观尺度进行基函数离散后，全局问题为以下形式：
[0035]
[0036] 步骤6、在微观单元上进行独立的时间和空间离散，得到微观单元的基函数以及微观局部 时间步长，并结合实际的网格特征，集成微观单元矩阵，构造微观“局部问题”，用以多尺度 基函数的求解。
[0037] 步骤6-1、在局部微观单元上进行独立的时间和空间离散。
[0038] 将每个宏观单元K∈TH做为独立的有限元空间进行离散，离散后的区域记为 每一个 微观单元的离散区域记为k，则有 令φ表示微观细单元上的基函数， 为宏观单元 K内局部微观网格的总自由度，则多尺度基函数在微观尺度离散为：
[0039]
[0040] 在时间离散方面，首先选取基于宏观网格的宏观时间间隔△t:＝[tn-1,tn]，因为每一个局部 单元的微观网格对时间步长的要求不同，因此需要根据微观单元的需要设置满足其要求的局 部时间间隔，局部时间间隔由宏观时间间隔进一步离散得到。令τK为宏观单元K子域内的局 部时间步长，为宏观时间间隔分为s个间隔，则有tn-1+sτK＝tn。
[0041] 步骤6-2、构造微观“局部问题”。
[0042] 采用时域间断伽辽金法DGTD构造微观“局部问题”，二次离散后宏观单元K内的微观单 元面集合记为 令 表示微观单元的面。对于不包含宏观单元K的面边界的面集合记 为 结合式(3)与间断伽辽金法可得：
[0043]
[0044] 通过线性分解，将式(10)进行分解，可以得到局部问题：
[0045]
[0046]
[0047] 其中{e}与{h}表示数值通量，因为微观问题本质上是通过DGTD方法构建的，因此数值 通量可以使用中心通量、迎
风通量或罚项通量。在微观单元与宏观单元交界处{e}与{h}的取 值为实际电磁场的值e与h。而[|φ|]＝nk×φk+nk'×φk'，其中k'表示与k相邻的微观单元编号。
[0048] 步骤6-3、对各微观单元进行相应求解，集成各微观单元矩阵。
[0049] 步骤7、在每一个微观单元采用其自己独立的微观时间步长进行局部问题求解的时间迭 代，在所有局部问题迭代求解完成后，将得到的多尺度基函数代入宏观全局问题，用宏观时 间步长进行全局问题的求解。
[0050] 步骤8、更新电场和磁场，并继续通过时间迭代得到所需的时域电磁场的解。
[0051] 进一步的，所述步骤7中宏观时间步长与网格最大的微观单元时间步长一致，以使得整 个求解域的时间迭代效果最佳。
[0052] 综上所述，本发明通过将多尺度时域电磁问题分解为宏观尺度与微观尺度两部分，实现 了尺度分离，有效降低了原问题的求解规模，并且在微观单元采用适合其特点的混合网格进 行离散，进一步降低了其求解自由度。在时间迭代方面，对于不同的微观单元采用了独立的 时间步长进行时间离散，有效避免了局部微观网格尺寸过小导致时间迭代次数过多的缺点。 在通过微观问题求解出微观信息后，实现宏观问题与微观问题之间的耦合，将其代入宏观问 题并采用宏观时间步长进行迭代求解。最终达到高效求解多尺度时域电磁问题的技术效果。
附图说明
[0055] 图3为实施例的局部微观单元剖分图。
[0056] 图4为实施例的等参坐标变换示意图。
[0057] 图5为实施例的平面波电场分布结果。
[0058] 图6为实施例的观察点A代码计算结果与理论值对比。
[0059] 图7为实施例的L2范数误差。
具体实施方式
[0060] 下面结合附图对本发明的实施做进一步的说明。
[0061] 参阅图1，本发明提供了一种基于混合网格及混合时间步长的时域电磁学多尺度计算方 法，具体流程如下：
[0062] 步骤1、本实施例中，以边长为0.5m的正方形区域传播例进行分析，根据要解决的实际 电磁场的类型，选择计算模型进行计算。
[0063] 步骤2、对于步骤1所建计算模型，根据实际电磁问题的需要，忽略所有的微观细节部分， 在宏观尺度上进行网格剖分，将整个求解域(计算模型)划分为K＝64个互不重叠的宏观单 元，采用的宏观网格划分如图2所示，且每个宏观单元边被进一步分为两段。宏观单元的尺 寸大于求解域内的最小特征尺寸。
[0064] 步骤3、考虑所有微观结构及微观信息，根据问题的需要使用混合网格对宏观单元进行第 二次网格剖分，得到微观细网格。
[0065] 微观单元的剖分在宏观单元离散基础上进行，微观单元的大小不大于宏观单元内部微观 结构的特征尺度。在二次剖分时，采用了混合网格，需要根据实际情况选择适合的微观细网 格进行剖分，对于内部没有微观结构的部分，采用四边形网格剖分，对于内部具有微观结构 的部分采用三角形网格剖分。最终采用的微观网格划分如图3所示。
[0066] 步骤4、通过伽辽金方法构造麦克斯韦方程组的杂交弱形式，并对电磁场进行线性分解。
[0067] 步骤4-1、建立求解域的电磁问题。
[0068] 设求解域为Ω∈R3, 为求解域边界，其中ΓPEC为理想电壁边界，ΓABC为吸收边界， 则求解域内的电磁问题如式(1)。
[0069] 步骤4-2、通过伽辽金方法建立麦克斯韦方程组的杂交弱形式。
[0070] 求解域Ω进行宏观网格离散后区域为TH，区域内的离散后的单元K∈TH。假定Pm(K) 为单元内维数不小于m的多项式空间，用V来表示空间离散，则有V:＝{v∈L2(Ω):v|K∈[Pm(K)]3for_all(K∈TH)}，麦克斯韦方程组的伽辽金弱形式可以表示 为式(2)。
[0071] 离散后的所有宏观单元边界面集合记为εH，在单元边界面 存在切向分量函数 空间Λ，可以表示为 引入杂交变量λ∈Λ，并将式(2)变为杂交 形式(3)。
[0072] 步骤4-3、根据宏观问题与微观问题的需要对电磁场进行线性分解；
[0073] 将电场与磁场进行线性分解，分解为两个部分，如式(4)所示。
[0074] 步骤5、引入多尺度基函数，并进行基函数离散，之后在宏观尺度进行分析，构造宏观“全 局问题”。
[0075] 步骤5-1、引入多尺度基函数，并进行基函数离散。
[0076] 宏观问题主要建立在宏观单元交界面上，首先令ψ为宏观粗网格面上的基函数，并且令 dim(ΛH)表示宏观单元交界面总的自由度维数。因为与宏观网格相关的量主要有及 所以在宏观尺度上对杂交变量λ及电磁场 离散后可表示为式(5)。
[0077] 步骤5-2、通过宏观尺度进行分析，构造宏观“全局问题”。
[0078] 宏观问题的本质目的是建立整个求解域的连接性，使每一个局部单元可以彼此耦合，因 此，局部单元内部在边界处的通量应该与整个宏观部分的通量，因为该求解区域只存在吸收 边界与平面波激励，因此需要吸收边界及激励引入全局问题，引入吸收边界及激励后全局问 题如式(7)所示，在宏观尺度进行基函数离散后，全局问题可以写成以式(8)所示形式。
[0079] 步骤6、在微观单元上进行独立的时间和空间离散，得到微观单元的基函数以及微观局部 时间步长，并结合实际的网格特征，集成微观单元矩阵，构造微观“局部问题”，用以多尺度 基函数的求解。
[0080] 步骤6-1、在局部微观单元上进行独立的时间和空间离散。
[0081] 将每个单元K∈TH做为独立的的有限元空间进行离散，离散后的区域记为 每一个小 的离散区域记为k，则有 令φ表示微观细单元上的基函数，并令 为宏观单元K内 局部微观网格的总自由度，则多尺度基函数可以在微观尺度离散如式(9)所示。
[0082] 在时间离散方面，将宏观时间间隔与矩形微观网格单元区域的时间间隔保持一致，三角 形微观网格单元区域的时间间隔选取为矩形部分的1/3。
[0083] 步骤6-2、采用时域间断伽辽金法DGTD构造微观“局部问题”，二次离散后宏观单元K 内的微观单元面集合记为 令 表示微观单元的面。对于不包含宏观单元K的面边界 的面集合记为 结合式(3)与间断伽辽金法可得式(10)。
[0084] 通过线性分解，将式(10)进行分解，可以得到如式(11)(12)所示局部问题：
[0085] 其中{e}与{h}表示数值通量，因为微观问题本质上是通过DGTD方法构建的，实施例中 数值通量采用迎风通量。在单元边界处{e}与{h}的取值为实际电磁场的值e与h。而 [|φ|]＝nk×φk+nk'×φk'，其中k'表示与k相邻的微观单元编号。
[0086] 步骤6-3、集成局部微观单元矩阵；
[0087] 在本实施例中，微观单元网格由三角形网格与四边形网格构成，对于三角形网格，其质 量矩阵 与
刚度矩阵 都有解析解，因此我们可以得到：
[0088]
[0089]
[0090] 其中，Sk为单元的面积，δ(i,j)为一个函数，当i＝j时值为1，否则值为0。大小三 角形边的长度，方向指向三角形内部。
[0091] 对于四边形网格，需要用到等参变换进行求解，等参变换如图4所示，首先引入参数坐 标(ξ,η)，其与直角坐标(x,y)的坐标转换关系为：
[0092]
[0093] 其中， (ξ1,η1)＝(-1,-1),(ξ2,η2)＝(1,-1)， (ξ1,η1)＝(1,1),(ξ1,η1)＝(-1,1)。
[0094] 因此四边形网格的
质量矩阵 可以写成参数坐标形式：
[0095]
[0096] 其中[J]为雅克比矩阵，在二维情况下为：
[0097]
[0098] 对于四边形网格刚度矩阵 根据偏微分理论，有：
[0099]
[0100] 因此，刚度矩阵 可以写为：
[0101]
[0102] 在这里，质量矩阵 与刚度矩阵 都没有解析解，因此需要通过数值积分进行求解。
[0103] 对于通量系数矩阵Fh(仅为系数，不包括向量部分)，在二维情况下是一个线积分，对于 线积分基函数可以表示为：
[0104]
[0105] 其中ξ为坐标点在线段单元上的相对
位置，因此，通量系数矩阵F可以表示为：
[0106]
[0107] 系数m和n表示单元编号，若m＝n，则为同一单元面积分，若m≠n，则为相邻单元。 之后通过解析方法或数值积分方法进行求解即可。三角形单元与矩形单元的通量系数矩阵都 可以通过该方法得到。
[0108] 另外，局部问题中还有一个与杂交变量相关的矩阵 这个矩阵在二维情况下为 线积分，为一个宏观单元线段基函数ψ与微观单元线段基函数φ的积分，因为这两个基函数 的尺度不同，因此同样需要数值积分进行计算。
[0109] 最终可以得到，微观局部问题(11)与(12)的最终的矩阵形式如下：
[0110]
[0111]
[0112] 其中k表示计算单元，k'表示相邻单元，l表示所在单元与相邻单元的集合。
[0113] 步骤7、在每一个微观单元采用其自己独立的微观时间步长进行局部问题求解的时间迭 代，在所有局部问题迭代求解完成后，将得到的多尺度基函数代入宏观全局问题，用宏观时 间步长进行全局问题的求解。
[0114] 步骤7-1、在每一个局部微观区域采用自己独立的微观时间步长进行局部问题求解的时间 迭代。
[0115] 在进行时间迭代时，首先对局部问题进行时间迭代，对于每个局部单元K，用下标n表示 第n个宏观时间步，用上标s表示第s个局部时间步。对于矩形微观单元区域，s＝1；对于三 角形网格区域，s＝3，其时间迭代形式如下：
[0116]
[0117]
[0118] 在实际计算中，式(24)只需要在第一个宏观时间步内进行计算，在之后的时间迭代中只需 要计算式(25)即可，局部时间迭代完成后，需要将求解得到的多尺度基函数数值解eλ,ef,hλ,hf代入宏观问题进行求解。
[0119] 步骤7-2、用宏观时间步长进行全局问题的求解。
[0120] 求解宏观问题时，首先需要通过微观单元矩阵组成宏观单元矩阵，将局部小矩阵Mh、Sh组合成全局大矩阵形式MH＝{Mh}、SH＝{Sh}，并将局部问题求解得到的结果 也组装为 全局形式 并在宏观边界所对应的单元求出吸收边界及激励所对应的 全局系数大矩阵SA＝{sa}(吸收边界系数全局大矩阵)Sinc＝{sinc}(激励系数全局大矩阵) 以及激励Ginc＝{ginc}，结合式(8)可得宏观全局问题的矩阵形式：
[0121]
[0122] 求解得到宏观变量 后，即可进行新时刻的电场与磁场的 更新。
[0123] 步骤8、更新电场和磁场，并继续通过时间迭代得到所需的时域电磁场的解。
[0124] 在完成一次整个求解域的时间迭代后，需要结合宏观变量与微观多尺度基函数对电场与 磁场进行还原更新，新时刻电场与磁场的更新方法为：
[0125]
[0126] 其中，sn表示每次局部时间迭代的最终时刻， 与 只需要在第一个宏观时间步内 计算即可。
[0127] 在更新完所有的电场与磁场后，需要进行下一时刻的计算，从而得到每一时刻的时域电 磁场的值。
[0128] 计算得到电场分布结果如图5所示，计算域内观察点A的电场理论值与代码计算值对比 如图6所示，L2误差范数误差图7所示。可以看出，计算得到的电场的分布符合实际平面波 电场分布规律，理论值与代码计算值也稳定的很好，且当稳定时L2范数误差稳定在一个很小 的值，结果正确。
[0129] 综上可见，本发明采用的
算法是完全可行的，虽然该实施例较为简单，但其实现过程涵 盖了本发明的所有步骤，该算法可用于更为复杂的电磁学问题的计算，最终达到高效求解多 尺度时域电磁问题的技术效果。
