CN110580365B - 基于叠层矢量基函数的动态p自适应DG-FETD方法 - Google Patents

Abstract

本发明公开了一种基于叠层矢量基函数的动态p自适应DG‑FETD方法，包括：建立求解模型，使用四面体网格对模型进行离散，得到模型的结构信息；设置仿真参数，读取模型的结构信息；对一阶麦克斯韦旋度方程等式两边采用伽辽金法测试，并用基函数展开，得到最终的迭代公式；通过蛙跳差分格式时域有限元的迭代公式进行时间迭代，每一次时间步迭代开始时计算空间场值波动参数，调整基函数阶数；所有时间步迭代结束得到空间中的电场值和磁场值。本发明实现了不连续伽辽金时域有限元法中高低阶基函数的混合，并且实现了不同离散区域基函数自行选取合适阶数的功能，在保证精度的情况下有效地节省时域有限元法的仿真计算时间，具有很强的实际工程应用价值。

Description

基于叠层矢量基函数的动态p自适应DG-FETD方法

技术领域

本发明属于有限元法基函数自适应领域，特别是一种基于叠层矢量基函数的动态p 自适应DG-FETD方法。

背景技术

在多尺度问题中，在求解某些对精度要求较高的结构时，我们有两种选择：减小剖分尺寸使用更密集的网格进行离散或者增加基函数的阶数。

前一种方法会带来更多的未知量，时间步长的选取也会受到最小剖分尺寸的限制，而这些问题对于高阶基函数影响较小，相比而言后一种方法使用的更加普遍。但在面对电大问题的时候，高阶基函数方法计算时间仍然很长。

发明内容

本发明的目的在于提供一种基于叠层矢量基函数的动态p自适应DG-FETD方法，在高阶基函数的基础上研究了高低阶基函数混合并且阶数实现自适应的方法，能够更为快速高效的解决问题。

实现本发明目的的技术解决方案为：一种基于叠层矢量基函数的动态p自适应DG-FETD方法，步骤如下：

第一步，建立求解模型，使用四面体网格对模型进行离散，得到模型的结构信息，包括四面体的节点信息和单元信息，其中节点信息包括节点序号、节点坐标，单元信息包括单元序号以及该单元内包含的节点序号；

第二步，设置仿真参数，读取模型的结构信息；

第三步，基于以电场强度和磁场强度为未知量的一阶麦克斯韦旋度方程，对方程等式两边采用伽辽金法测试，并用基函数展开，得到最终的迭代公式，根据迭代公式填充计算迭代矩阵；

第四步，通过蛙跳差分格式时域有限元的迭代公式进行时间迭代，每一次时间步迭代开始时计算空间场值波动参数，依据该参数调整基函数阶数，用新的基函数进行计算；所有时间步迭代结束得到空间中的电场值和磁场值；

第五步，根据计算出的场值继续求解需要的物理参数。

本发明与现有技术相比，其显著优点：(1)本发明可以实现时域有限元法中基函数阶数的动态自适应；(2)相比原有的DG-FETD方法，本发明大大提高了计算效率，减少了计算时间。

附图说明

图1是金属球双站RCS示意图。

图2是观察点时域波形示意图。

图3是观察点基函数阶数随时间变化图。

具体实施方式

一种基于叠层矢量基函数的动态p自适应不连续伽辽金时域有限元(DG-FETD)方法，步骤如下：

第一步，建立求解模型，使用四面体网格对模型进行离散，得到模型的结构信息，包括四面体的节点信息以及单元信息；其中节点信息包括节点序号、节点坐标，单元信息包括单元序号以及该单元内包含的节点序号；

第二步，设置仿真参数，读取模型的结构信息；

第三步，根据电场强度E和磁场强度H为未知量的一阶麦克斯韦旋度方程：

上式中ε、μ分别表示离散单元的介电常数和磁导率，等式(1)(2)采用Whitney-I型的叠层矢量基函数进行伽辽金法测试，经由矢量恒等式和散度定理，得：

式(3)中，

测试函数，θV为单元V的面，

为此面上的法向量，式(4)中，

为测试函数。

在不连续伽辽金时域有限元法中，两个单元在单元分界面处是不连续的，那么我们就需要引入中心通量强加电磁场的切向连续性，中心通量如下：

其中，E⁺和H⁺分别表示和此单元共用边界θV的相邻单元上的电场和磁场。

将式(5)(6)代入式(3)(4)，并用Whitney-I型的叠层矢量基函数将E和H展开，得到：

得矩阵方程：

其中，

[T_hh]_ij表示矩阵[T_hh]中第i行第j列的元素；

在填充[T_hh]和[T_ee]时，按照四面体单元的编码进行顺序填充，因为每一个单元与其他单元没有相互作用，故可得到具有块对角特性的[T_hh]和[T_ee]；

第四步，将式(9)和式(10)用蛙跳差分格式展开得到：

[T_ee]·eⁿ⁺¹＝Δt([P_eh]+[S_eh]+[Ss_eh])·h^n+1/2+[T_ee]·eⁿ (11)

[T_ee]·eⁿ⁺¹＝Δt([P_eh]+[S_eh]+[Ss_eh])·h^n+1/2+[T_ee]·eⁿ (12)

其中e和h分别代表本单元中的电场值和磁场值，角标代表时刻。

上式为最后的迭代公式。

整个计算区域内所有的单元在每个时间步开始，

会计算一个关于空间场值波动的参数：

其中m_a表示四面体单元第a条棱边的中点，c表示四面体单元的重心，e表示该点的电场值，d_a表示第a条棱边中点到重心的距离。

在程序最初设定了P^L和P^H，其中P^L表示整个计算区域内基函数的最低阶数，反之P^H表示整个计算区域内基函数的最高阶数，本实施例中P^L取1，P^H取5；然后在程序迭代的最初我们把所有单元的基函数阶数设置为P^k＝P^L。在时间步t时刻，计算得到一个给定单元每条棱边的u_a，只要此单元中任意一个u_a超过了阈值ζ(本实施例中设置为 10^-4，可自行调节)，那么在下一个时间步中，此单元的基函数阶数增加1，反之，如果此单元中没有一个u_a超过ζ，那么单元的基函数阶数保持不变。

一旦一个单元阶数发生了改变，那么就称作p自适应，在随后的时间步里p自适应需要遵守以下的规则：(1)在某单元里若所有u_a满足u_a≤ζ，且在P^k＝P^L的情况下，或者若存在u_a满足u_a＞ζ，且在P^k＝P^H的情况下，那么此单元不发生基函数阶数上的变化；(2)在某单元里若存在u_a满足u_a＞ζ，且在P^k＜P^H的情况下，那么此单元基函数阶数加1；(3)在某单元里若所有u_a满足u_a≤ζ，且在P^k＞P^L的情况下，那么此单元基函数阶数减1。

得到新的基函数阶数后，接着按照前面的方法进行电场值和磁通量更新，空间中每点均可得到正确的场值，至此一个时间步计算全部完成。重复上述步骤直至时间迭代结束。

第五步，利用得到的场值信息计算双站RCS，因为与DG-FETD方法中后处理一样，此处不再赘述。

为了验证本发明的正确性与有效性，下面分析了一个金属球的散射特性。

计算一个半径为0.5m的金属球，球的中心点取在坐标原点。加入的调制高斯源频带范围从10MHz到600MHz，中心频率取300MHz，脉宽相关量取2.3ns，平面波沿着 Z轴正方向入射，电场为X方向极化，取观察点坐标(1.1m,1.1m,1.1m)。用加入p自适应的程序计算出300MHz的双站RCS如图1，观察点处的时域波形以及阶数随时间的变化如图2、图3。且保证精度相同的情况下，p自适应时域有限元法程序耗时923s，传统DG-FETD程序耗时3102s。

Claims (1)

1.一种基于叠层矢量基函数的动态p自适应DG-FETD方法，其特征在于，步骤如下：

第一步，建立求解模型，使用四面体网格对模型进行离散，得到模型的结构信息，包括四面体的节点信息和单元信息，其中节点信息包括节点序号、节点坐标，单元信息包括单元序号以及该单元内包含的节点序号；

第二步，设置仿真参数，读取模型的结构信息；

第三步，基于以电场强度和磁场强度为未知量的一阶麦克斯韦旋度方程，用基函数展开电场强度和磁场强度，根据伽辽金法对方程等式两边采用测试基函数测试，得到最终的迭代公式，根据迭代公式填充计算迭代矩阵；

具体实施步骤：首先基于电场强度E和磁场强度H的一阶麦克斯韦旋度方程(1)、(2)：

上式中ε、μ分别表示离散单元的介电常数和磁导率，E和H分别表示电场强度和磁场强度；

接着对方程等式进行伽辽金法测试，并引入中心通量强加相邻单元之间电磁场连续性，再将E和H用基函数展开，得：

其中e和h分别代表本单元中的电场值和磁场值，[T_hh]、[P_he]、[S_he]、[Ss_he]、[T_ee]、[P_eh]、[S_eh]、[Ss_eh]分别为形成的矩阵，角标有关矩阵的第一、二维长度；

等式(1)(2)采用Whitney-I型的叠层矢量基函数进行伽辽金法测试，经由矢量恒等式和散度定理，得：

式(5)中，

为测试函数，

为单元V的面，

为此面上的法向量，式(6)中，

为测试基函数；

在不连续伽辽金时域有限元法中，两个单元在单元分界面处是不连续的，引入中心通量强加电磁场的切向连续性，中心通量如下：

其中，E⁺和H⁺分别表示和此单元共用边界

的相邻单元上的电场和磁场；

将式(7)(8)代入式(5)(6)，并用Whitney-I型的叠层矢量基函数将E和H展开，得到：

最终得紧凑格式的矩阵方程：

其中，

[T_hh]_ij表示矩阵[T_hh]中第i行第j列的元素；

第四步，通过蛙跳差分格式时域有限元的迭代公式进行时间迭代，每一次时间步迭代开始时计算空间场值波动参数，依据该参数调整基函数阶数，用新的基函数进行计算；所有时间步迭代结束得到空间中的电场值和磁场值；

具体实施步骤：首先通过时域有限元的迭代公式进行时间迭代，使用蛙跳差分格式，将式(3)和(4)分别按照蛙跳差分格式展开得到：

[T_hh]·h^n+1/2＝-Δt([P_he]+[S_he]+[Ss_he])·eⁿ+[T_hh]·h^n-1/2 (13)

[T_ee]·eⁿ⁺¹＝Δt([P_eh]+[S_eh]+[Ss_eh])·h^n+1/2+[T_ee]·eⁿ (14)

其中e和h分别代表本单元中的电场值和磁场值，角标代表时刻；

整个计算区域内所有的单元在每个时间步开始，计算一个关于空间场值波动的参数：

其中m_a表示四面体单元第a条边的中点，c表示四面体单元的重心，e表示该点的电场值，d_a表示第a条棱边中点到重心的距离；

依据空间场值波动参数的大小，调整各自单元内基函数的阶数，实现单元基函数阶数的动态自适应，然后按照前面的方法进行电场值和磁通量更新，空间中每点均可得到正确的场值，至此一个时间步计算全部完成；重复上述步骤直至时间迭代结束；

然后依据空间场值波动参数的大小，调整各自单元内基函数的阶数，具体为:

初始设定P^L和P^H，其中P^L表示整个计算区域内基函数的最低阶数，P^H表示整个计算区域内基函数的最高阶数；然后在迭代的最初把所有单元的基函数阶数设置为P^k＝P^L；在时间步t时刻，计算得到一个给定单元每条棱边的u_a，若单元中任意一个u_a超过阈值ζ，那么在下一个时间步中，此单元的基函数阶数增加1，反之，如果此单元中没有一个u_a超过ζ，则单元的基函数阶数保持不变；一旦一个单元阶数发生改变，称作p自适应，在随后的时间步里p自适应需要遵守以下的规则：(1)在某单元里若所有u_a满足u_a≤ζ，且在P^k＝P^L的情况下，或者若存在u_a满足u_a＞ζ，且在P^k＝P^H的情况下，那么此单元不发生基函数阶数上的变化；(2)在某单元里若存在u_a满足u_a＞ζ，且在P^k＜P^H的情况下，那么此单元基函数阶数加1；(3)在某单元里若所有u_a满足u_a≤ζ，且在P^k＞P^L的情况下，那么此单元基函数阶数减1；

第五步，根据计算出的场值继续求解需要的物理参数。

Images