技术领域
[0001] 本
发明属于三维时域电磁学数值求解技术领域，涉及一种三维时域电磁学杂交时域间断伽辽金数值方法，基于高阶叠层基函数。
背景技术
[0002] 在复杂时域电磁环境下，往往需要求解多尺度装备的瞬态电
磁场。而这些多尺度装备包含着各种非结构的网格，既有尺寸远小于
波长的精细网格，也有类比于甚至远大于波长的均匀网格。如军用飞机、舰船等，其整体结构较大，却又包含很多像孔缝、天线这样的精细结构。而这种存在精细网格的局部加密网格，对时域数值
算法的时间步长与计算性能有着严重的影响。比如显式时间
迭代格式需要非常小的时间步长来保证它的条件
稳定性，但是大大地增加了计算时间；虽然隐式时间迭代格式是无条件稳定的，可以增大时间步长，但是它需要求解一个全局线性系统，导致内存的消耗。特别地，对于高阶问题，随着模型网格数量的增大，全局矩阵的维数会更大，并且全局矩阵很有可能是高度病态的，这对矩阵求解造成了很大的难度。因此，如何选用一种高效的时域电磁数值方法来求解这种局部加密网格是至关重要的。
[0003] 然而，传统的数值方法，比如时域有限差分法、时域有限元法与时域有限体积法等往往会因为网格尺寸与收敛
精度而受到限制。近年来，在时域计算电磁学数值领域中出现了一种时域间断伽辽金法，该方法由于支持非结构网格，有着高精度、天然的高性能并行技术等而被广泛应用，目前已成为一个较为热
门的研究领域。因此，我们可以看出时域间断伽辽金法相对于传统的数值方法有着许多优势。
[0004] 然而，时域间断伽辽金法的所有优势都是基于这样一个条件：虽然在单元交界面上场是不连续的，但是每个单元需要维持自己的基函数。而这个条件导致了它的一个致命的缺点：单元交界面上的未知量是重复的，从而导致在获得相同精度下所需的全局未知量明显多于经典时域有限元方法。对于尺寸较大的均匀网格，采用显式时间迭代格式不会存在太大的问题，但是对于网格数量很多的大规模问题而言，则时域间断伽辽金法的未知量会变得很大。特别是求解局部加密网格，采用显式时间格式则需要更小的时间步长，进而严重增加计算时间与内存消耗。因此对于这种局部加密网格，现有的方法大都结合了隐式时间迭代格式。
[0005] 然而，正如前文所述如何快速有效地求解全局矩阵仍是隐式时间迭代格式的一个很大的难点。并且随着分析的结构复杂化，低阶基函数的精度已经不能满足设计者的要求，而高阶基函数的构造在无形中增加了全局未知量的数目及矩阵求解的困难，因此，迫切需要一种高效的时域电磁数值法来求解这种非结构局部加密网格的问题，在能保证高阶精度的前提下,提升计算性能。
发明内容
[0006] 针对上述存在问题或不足，为实现高效时域电磁数值法求解非结构局部加密网格时，在保证高阶精度的前提下，提升计算性能。本发明提出一种三维时域电磁学杂交时域间断伽辽金数值方法。利用该方法可以高效地求解非结构局部加密网格的三维时域Maxwell方程组，且具有较少的全局未知量、较高的精度以及显著的计算性能。
[0007] 具体包括以下步骤：
[0008] 步骤A、根据目标
电子器件的物理结构，结合工作环境与边界条件对其仿真建模；
[0009] 步骤B、采用四面体单元剖分三维求解区域，面离散和体积离散必须相容；
[0010] 步骤C、选择有限元标量叠层基函数，将
电磁场与杂交量用基函数展开；
[0011] 三维时域电磁学杂交时域间断伽辽金数值方法需要求解即先构成只与杂交量Λh有关的全局线性系统，求出杂交量后，根据局部系统来获得电磁场Eh和Hh， 为间断有限元函数空间， 为有限元迹空间；
[0012] 步骤D、在空间上，建立三维杂交时域间断伽辽金数值方法的半离散形式；
[0013] 引入杂交量来替换时域间断伽辽金法的数值迹，并增加第三个守恒方程来保证全局系统的稳定性，三维杂交时域间断伽辽金数值方法空间上的半离散形式：
[0014]
[0015] ε为计算区域Ω中介质的相对
介电常数，μ为计算区域Ω中介质的相对磁导率，n是计算区域边界 上的外法向单位矢量，τ＞0是局部稳定系数。ginc表示吸收边界条件(ABC)的边界方程， 和 是测试函数。 是切向磁场。
[0016] 步骤E、在时间上，视杂交量为常量，只考虑电磁场的时间离散，结合D步骤形成全离散方程形式，从而得到全局线性系统；
[0017] 定义时间步长Δt，将总仿真时间[0,T]离散为等间隔的时间步tn＝nΔt, 其中n的最大值为Nt，即时间迭代总次数。令tn时刻的电磁场 tn+1时刻的电磁场 采用隐式Crank-Nicolson时间格式，对于(1)式
中的时间偏导项有：
[0018]
[0019] 而不含时间偏导项的近似为：
[0020]
[0021] 在时间上，视杂交量为常量，只考虑电磁场的时间离散(2)与(3)，结合D步骤产生的空间半离散形式，通过进一步推导后有
[0022]
[0023] 其中(4)式的等式左边都关于tn+1时刻，而y1和y2是关于tn时刻，只有y3含有tn和tn+1时刻。具体的右端项形式为：
[0024]
[0025] 首先考虑局部线性系统:由于一个四面体 上含有4个三
角形面单元，我们定义τi上的所有面单元的杂交量为Λe。根据(4)式的前两个方程，推导出在一个四面体上的局部线性系统为
[0026]
[0027] 是tn+1时刻的电磁场， 是tn时刻的电磁场。 和 是通过基函数作用能得到的局部矩阵。结合(4)式的最后一个方程，得到形成全局线性系统的基本方程[0028]
[0029] 其中 是通过基函数作用能得到的局部矩阵，be是局部右端项。根据(7)式,依次
叠加每个四面体单元后，得到一个全局线性系统，即
[0030]
[0031] 矩阵 为全局线性矩阵，y是右端项。
[0032] 步骤F、求解全局线性系统杂交量。
[0033] 进一步的,所述步骤E完成后对其构造p型多重网格进行预处理，再进行步骤F。具体如下：
[0034] 从(8)式出发，根据叠层基函数的性质，重新构造全局矩阵 如下[0035]
[0036] 子矩阵 表示叠层基函数的低阶部分的自作用矩阵，且可以看做是粗网格矩阵。表示叠层基函数的高阶部分的自作用矩阵。而 和 则是低阶与高阶基函数相互作用的耦合矩阵。整个 可以看成是一个细网格矩阵。由于 中包含不同的基函数，而不同阶数的基函数性质不同。因此在矩阵求解前，将矩阵 的对角元进行归一化。然后运用Schur分解，我们可以得到
[0037]
[0038] 其中 是Schur补矩阵。先对Schur补矩阵做一个近似，即然后对两个子矩阵 和 进行不完全Choleski分解，大量的数值例子表明对于和 在近似时所需抛弃因子的
阈值分别为10-5和10-2，进一步得到 和
结合(10)式，得到最终的p型多重网格的预处理矩阵
[0039]
[0040] 在时间迭代开始前，将全局矩阵 分解一次即可用于之后每次的时间迭代求解；根据(11)式，我们可知在整个时间迭代开始前，只需要对矩阵 和 进行一次不完全Choleski分解；将本步骤提出的p型多重网格的预处理技术运用到E步骤的全局线性系统的求解部分，通过全局线性系统(8)获得Λ，根据每个四面体上局部杂交量与全局杂交量的映射关系，得到每个四面体上的Λe，从而根据局部线性系统(6)式来获得每个单元的电磁场待求系数
[0041] 本发明将频域的杂交间断伽辽金法引入到时域中，并以一种更为简单的方式来形成全离散形式，以提升计算性能；构造高阶插值叠层基函数来获得高精度的数值模拟结果，然后构造p型多重网格的预处理来
加速全局线性矩阵的求解，以提升计算性能。
[0042] 综上所述，本发明求解非结构局部加密网格时，在保证高阶精度的前提下，提升了计算性能。
附图说明
[0044] 图2是间断有限元函数空间 与有限元迹空间 的示意图。
具体实施方式
[0045] 下面结合附图和
实施例进一步的详细说明本发明。
[0046] 参照图1，一种基于函数逼近自适应误差分析的三维
微波管输入输出窗模型降阶的数值方法，包括以下步骤：
[0047] 步骤A.将目标电子器件结构结合材料特性进行仿真建模；
[0048] 根据目标电子器件的物理结构，结合工作环境与边界条件对其仿真建模。
[0049] 步骤B.采用四面体网格离散求解域；
[0050] 采用四面体单元剖分三维求解区域是一种公知过程，因此本步骤不再详细描述。需要注意的是，面离散和体积离散必须相容。本实施例中的计算区域 被划分成Nh个四面体网格的集合 其中每个体单元用τi(i=1,2,3,…,Nh)表示，即
[0051]
[0052] 并且面集合 是由Nf个三角形面单元Df组成，即
[0053]
[0054] 步骤C.选择标量叠层基函数，将电磁场与杂交量用基函数展开；
[0055] 首先，在三维时域电磁学杂交时域间断伽辽金数值方法中，所用到的间断有限元函数空间 与有限元迹空间 是一种公知的空间，因此其具体形式不再详细描述。这里为了更好说明二者的差异，在图2中给出了这两种空间不同阶数的示意图，其中 和分别表示这二个空间的
自由度。
[0056] 一般的时域间断伽辽金方法，只需要求解电磁场 然而，三维时域电磁学杂交时域间断伽辽金数值方法需要求解 即先构成只
与杂交量Λh有关的全局线性系统，求出杂交量后，根据局部系统来获得电磁场(Eh,Hh)。
[0057] 选择标量叠层基函数，将电磁场与杂交量(Eh,Hh,Λh)用基函数展开。对于每个四面体 定义局部电磁场为 且
[0058]
[0059] 为了实现高阶精度和p型多重网格的预处理技术，这里采用标量叠层基函数来展开(14)式中的各个方向上的场分量，即
[0060]
[0061] 其中 和 是电磁场的待求系数， 的维数为di＝(p+1)(p+2)(p+3)/6。在本实施例中，以二阶叠层基函数为例进行说明。故di＝10，根据二阶叠层基函数的性质，我们可知 是
节点类型，而 是棱边类型。对于一个四面体上的二阶叠层基函数具体形式如下：
[0062]
[0063] 定义该四面体的电磁场为We＝[Ee,He]T，其对应的待求系数为We，且未知量个数为6di。
[0064] 杂交量Λh的基函数，由于每个面上的杂交量具有单值性，对于面集合上的杂交量可以表示为：
[0065]
[0066] 这里f∈[1,Nf]是一个面单元 的全局面编号，对应的局部杂交量为Λf，具体形式为：
[0067]
[0068] 其中uf和wf是面坐标，且该面上外法向量单位矢量n＝uf×wf。这里采用标量叠层基函数 来展开(18)式中分量，即
[0069]
[0070] 其中 和 是局部杂交量Λf的待求系数， 的维数为df＝(p+1)(p+2)/2。在本实施例中，以二阶叠层基函数为例进行说明。故df＝6，根据二阶叠层基函数的性质，我们可知对于一个面上的二阶叠层基函数具体形式如下：
[0071]
[0072] 定义局部杂交量Λf的待求系数Λf，且未知量个数为2df。
[0073] 综上，我们已经详述完关于电磁场与杂交量的叠层基函数。最后，从未知量个数的角度来说明杂交时域间断伽辽金法相比时域间断伽辽金法的优势。由于时域间断伽辽金法的未知量个数是与电磁场相关，如前所述一个四面体上的未知量个数为6di，因此在整个求解区域上Eh和Hh的未知量个数为6diNh。对于杂交时域间断伽辽金法，其未知量个数只与杂交量相关，由于一个面单元上的未知量个数为2df，因此在整个面集合上Λh的未知量个数为2dfNf。注意到di/df＝(p+3)/3，在实际问题中，离散后的四面体个数与面个数往往不同，但一般而言，我们可以估计二者的关系如下：Nh＝Nf/2。
[0074] 可知时域间断伽辽金法与杂交时域间断伽辽金法的未知量个数之比为(p+3)/2。而这说明了随着阶数的增加，杂交时域间断伽辽金法相比时域间断伽辽金法需要更少的未知量，特别是对于复杂的含有许多网格的模型，杂交时域间断伽辽金法的这一优势就更加突出了。
[0075] 步骤D.在空间上，建立三维杂交时域间断伽辽金数值方法的半离散形式；
[0076] 与传统的时域间断伽辽金法在空间上的离散不同，杂交时域间断伽辽金法需要引入杂交量来替换时域间断伽辽金法的数值迹，并需要增加第三个守恒方程来保证全局系统的稳定性。其具体的公式推导与频域杂交间断伽辽金法类似，是一种公知过程，这里不再阐述。下面只给出三维杂交时域间断伽辽金数值方法的半离散形式：
[0077]
[0078] 这里的ε为计算区域Ω中介质的相对介电常数，μ为计算区域Ω中介质的相对磁导率，n是计算区域边界 上的外法向单位矢量，τ＞0是局部稳定系数。 表示吸收边界条件(ABC)的边界方程， 和 是测试函数， 是切向磁场。
[0079] (21)式只是空间上的半离散形式，还需要处理其中的时间偏导项，即下面的步骤E的时间离散过程，从而得到全离散形式，形成全局线性系统。
[0080] 步骤E.在时间上，视杂交量为常量，只考虑电磁场的时间离散，结合步骤D形成全离散方程形式，从而得到全局线性系统；
[0081] 定义时间步长Δt，将总仿真时间[0,T]离散为等间隔的时间步 其中n的最大值为Nt，即时间迭代总次数。令tn时刻的电磁场 tn+1时刻的电磁场 采用二阶隐式的Crank-Nicolson时间格式，对
于(21)式中的时间偏导项有：
[0082]
[0083] 而不含时间偏导项的近似为：
[0084]
[0085] 考虑到杂交量只存在于面单元上，并且保持单值，因此本发明将杂交量看成一个待求常量，进而提出一种更为简单的方式来构造全局线性系统。即在时间上，视杂交量为常量，只考虑电磁场的时间离散(22)与(23)，结合D步骤产生的空间半离散形式，通过进一步推导后有：
[0086]
[0087] 其中(24)式的等式左边都关于tn+1时刻，而y1和y2是关于tn时刻，只有y3含有tn和tn+1时刻。具体的右端项形式为：
[0088]
[0089] 从上面的离散过程，我们可以看出杂交时域间断伽辽金法的全离散形式包括两个系统：局部线性系统和全局线性系统。事实上，(24)的前两个方程可以构成局部线性系统，而第三个方程就是杂交时域间断伽辽金法的守恒条件，根据这一守恒条件我们可以得到只含杂交量的全局线性系统。一旦解得杂交量，那么每个单元的电磁场就可以通过局部线性系统获得。下面我们从这两个系统的角度进行进一步的分析。
[0090] 首先考虑局部线性系统。由于一个四面体 上含有4个三角形面单元，我们定义τi上的所有面单元的杂交量为Λe。根据(24)式的前两个方程，我们可以推导出在一个四面体上的局部线性系统为
[0091]
[0092] 这里 是tn+1时刻的电磁场， 是tn时刻的电磁场。 和 是通过基函数作用能得到的局部矩阵。结合(24)式的最后一个方程，我们可以得到形成全局线性系统的基本方程
[0093]
[0094] 其中 是通过基函数作用能得到的局部矩阵，be是局部右端项。根据(27)式依次叠加每个四面体单元后，我们可以得到一个全局线性系统，即
[0095]
[0096] 矩阵 就是全局线性矩阵，y是右端项。一旦我们通过全局线性系统(28)获得Λ，根据每个四面体上局部杂交量与全局杂交量的映射关系，就可以得到每个四面体上的Λe，从而根据局部线性系统(26)式来获得每个单元的电磁场待求系数 至此，我们已经详述完三维杂交时域间断伽辽金数值方法的全部过程。
[0097] 步骤F.构造p型多重网格的预处理技术来加速E步骤的全局线性系统的求解[0098] 虽然步骤E已经展示了三维杂交时域间断伽辽金数值方法的全离散形式，但是如何有效地加速全局线性系统(28)的求解是本发明的一个关键之处。正如前文所述，对于这种局部加密网格，虽然隐式时间迭代格式是无条件稳定的，可以增大时间步长，但是它需要求解一个全局线性系统，导致内存的消耗。特别地，对于高阶问题，随着模型网格数量的增大，全局矩阵的维数会更大，并且全局矩阵很有可能是高度病态的，这对矩阵求解造成了很大的难度。而对于一个矩阵方程的求解，一个好的预处理矩阵是加速数值迭代方法收敛和稳定的关键。然而现有的一些预处理技术往往会因为矩阵维数较大或者矩阵的负定性和高度病态性的本质原因，使得迭代法收敛变得非常缓慢，甚至是不收敛的，从而使得高效求解大型稀疏矩阵变得异常困难。比如直接法对负定病态矩阵进行LU分解然后再回代求解后，会导致稀疏性受到限制。
[0099] 由于多重网格法的收敛速度几乎可以和矩阵的维数无关，这里采用p型多重网格的预处理技术来加速E步骤的全局线性系统的求解。相比传统的h型多重网格需要受到网格尺寸的限制，这种p型多重网格是以基函数的叠层性来构造粗细两套网格，灵活性更强。下面从(28)式出发，根据二阶叠层基函数的性质，重新构造全局矩阵 如下[0100]
[0101] 这里的子矩阵 表示一阶的节点基函数的自作用矩阵，且可以看做是粗网格矩阵。 表示二阶叠层型基函数的高阶部分的自作用矩阵。而 和 则是低阶与高阶基函数相互作用的耦合矩阵。整个 可以看成是一个细网格矩阵。由于 中包含不同的基函数，而不同阶数的基函数性质不同。因此在矩阵求解前，将矩阵 的对角元进行归一化。然后运用Schur分解，我们可以得到：
[0102]
[0103] 其中 是Schur补矩阵。众所周知，预处理矩阵越接近于 的逆矩阵，则预处理效率更好，但是计算花销也随之增大。因此，为了减少计算开销，我们先对Schur补矩阵做一个近似，即 然后对两个子矩阵 和 进行不完全Choleski分解，大量的数值例子表明对于 和 在近似时所需抛弃因子的阈值分别为10-5和10-2，进一步得到 和 结合(30)式，得到最终的p型多重网格的预处理矩
阵：
[0104]
[0105] 由于本发明提出的三维杂交时域间断伽辽金数值方法在每次的时间迭代中全局矩阵 是不变的，且是稀疏对称矩阵。因此只需要在时间迭代开始前，将全局矩阵 分解一次即可用于之后每次的时间迭代求解。根据(31)式，我们可知在整个时间迭代开始前，只需要对矩阵 和 进行一次不完全Choleski分解。通过将本步骤提出的p型多重网格的预处理技术运用到E步骤的全局线性系统的求解部分，可以进一步减少计算时间与内存消耗，从而显著增加全局线性系统的求解速度，提高整个三维杂交时域间断伽辽金数值方法的计算性能。
