CN114528716A - 一种应用于多尺度电磁波问题分析的高效数值仿真方法 - Google Patents

一种应用于多尺度电磁波问题分析的高效数值仿真方法 Download PDF

Info

Publication number
CN114528716A
CN114528716A CN202210199606.7A CN202210199606A CN114528716A CN 114528716 A CN114528716 A CN 114528716A CN 202210199606 A CN202210199606 A CN 202210199606A CN 114528716 A CN114528716 A CN 114528716A
Authority
CN
China
Prior art keywords
grid
coarse
fine
formula
interface
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210199606.7A
Other languages
English (en)
Inventor
康祯
张瑞祥
刘俊峰
姚尧
孙文钊
祁沛垚
韩传高
马晓珑
余俨
王嘉玮
邵景晖
杨方
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Thermal Power Research Institute Co Ltd
Original Assignee
Xian Thermal Power Research Institute Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Thermal Power Research Institute Co Ltd filed Critical Xian Thermal Power Research Institute Co Ltd
Priority to CN202210199606.7A priority Critical patent/CN114528716A/zh
Publication of CN114528716A publication Critical patent/CN114528716A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种应用于多尺度电磁波问题分析的高效数值仿真方法，包括：步骤1，根据多尺度电磁波问题的几何结构，采用局部亚网格剖分技术对计算区域进行空间网格剖分；步骤2，采用FDTD方法对粗网格区域的电磁场量进行求解；步骤3，采用PITD方法对细网格区域的电磁场量进行求解；步骤4，粗网格区域与细网格区域分界面上电磁场量的计算；步骤5，重复上述步骤，完成所有迭代步数。本发明专利在两种算法的分界面上采用了嵌套网格技术，构建了一种稳定的信息交互策略，使得本专利所建立的混合算法在进行长时间迭代后，仍然能够稳定运行，具有很好的长时稳定性。

Description

一种应用于多尺度电磁波问题分析的高效数值仿真方法
技术领域
本发明属于计算电磁学领域，具体涉及一种应用于多尺度电磁波问题分析的高效数值仿真方法。
背景技术
时域有限差分(Finite-Difference Time-Domain Method，FDTD)方法具有简单、直观和通用等优点，是求解麦克斯韦方程最常用的全波电磁仿真工具之一。然而，由于受到Courant-Friedrich-Levy(CFL)稳定性条件和不断增大的数值色散误差的限制，传统FDTD方法在模拟探地雷达、可穿戴天线或传感器等具有多尺度特征的复杂电磁波问题时，需要采用细网格对整个计算域进行离散，存在内存需求大、迭代步数多及计算效率低等问题。因此，开发精确高效的数值工具来求解多尺度电磁波问题仍然是计算电磁学领域的一个重要研究方向。
采用基于局部亚网格技术的FDTD方法对多尺度问题进行求解，即有选择的使用局部细网格对精细结构及附近区域进行剖分，使用粗网格对其他区域进行剖分，从而达到降低网格数目及内存需求的目的。但是根据CFL稳定性条件的要求，该方法时间步长的选取受到细网格尺寸的限制，即当细网格的尺寸非常小时，时间步长的选取也非常小，从而导致基于亚网格技术的FDTD方法仍然存在迭代步数较多、执行时间较长、累计误差较大等劣势。
为了突破局部亚网格FDTD方法中CFL稳定性条件的限制，提出了多种时域混合数值算法，包括混合FDTD和交替方向隐式FDTD(Alternating-Direction-Implicit FDTD，ADI-FDTD)方法及混合FDTD和局部一维FDTD(Locally One-Dimensional FDTD，LOD-FDTD)方法等。这类混合算法的基本思想是利用ADI-FDTD方法及LOD-FDTD方法等无条件稳定的时域数值算法来求解局部细网格上的场分量，使得时间步长的选取能够由粗网格尺寸来决定。因此，这类混合算法能够在不同尺寸的空间网格上采用统一的较大的时间步长进行求解，从而避免了在时间上进行插值或外推而引起的额外数值误差。但是ADI-FDTD方法和LOD-FDTD方法的数值色散误差会随着时间步长的增大迅速增加，导致这类混合算法在使用大时间步长进行仿真时计算精度有限。
时域精细积分(Precise-Integration Time-Domain，PITD)方法虽然不是无条件稳定的，但其稳定性条件仍然比CFL稳定性条件要宽松得多。与ADI-FDTD方法和LOD-FDTD方法相比，PITD方法的数值色散误差更小，且几乎与时间步长的选取无关。由于PITD方法对内存的需求高于FDTD、ADI-FDTD和LOD-FDTD方法，因此也不适合单独求解多尺度问题。需要说明的是，已公布的发明专利《一种基于亚网格技术的电磁波时域高效数值混合算法》(申请号：201911387294.4)中所提出的混合算法被证明在迭代步数增加时，将不再收敛稳定，即长时稳定性较差。
发明内容
本发明的目的在于提出了一种应用于多尺度电磁波问题分析的高效数值仿真方法，该方法以亚网格剖分技术为基础，结合FDTD方法和PITD方法提出一种应用于多尺度电磁波问题求解的高效混合算法，该算法具有迭代步数少、内存需求低、长时稳定性好等优点。
本发明采用如下技术方案来实现的：
一种应用于多尺度电磁波问题分析的高效数值仿真方法，包括：
步骤1，根据多尺度电磁波问题的几何结构，采用局部亚网格剖分技术对计算区域进行空间网格剖分；
步骤2，采用FDTD方法对粗网格区域的电磁场量进行求解；
步骤3，采用PITD方法对细网格区域的电磁场量进行求解；
步骤4，粗网格区域与细网格区域分界面上电磁场量的计算；
步骤5，重复上述步骤，完成所有迭代步数。
本发明进一步的改进在于，步骤1中，对于二维横电波TE模式，精细结构部分以及附近区域使用局部细网格进行剖分，其他区域则使用粗网格进行剖分。
本发明进一步的改进在于，粗网格的尺寸为所计算问题电磁波波长的1/10到1/8，细网格尺寸小于等于粗网格尺寸的1/2，粗网格尺寸与细网格尺寸的比值为nf。
本发明进一步的改进在于，步骤2中，二维TE波模式下，粗网格区域麦克斯韦方程的展开式为：
Figure BDA0003527025810000031
Figure BDA0003527025810000032
Figure BDA0003527025810000033
式中，Exc、Eyc分别为粗网格上的x方向和y方向的电场分量，Hzc为粗网格上的z方向的磁场分量，ε为计算区域介质的介电常数，σ为计算区域介质的电导率，μ为计算区域介质的磁导率；
FDTD方法在求解电磁场量时，同时对时间偏微分算子和空间偏微分算子进行中心差分离散，从而将偏微分方程组转化为代数方程组，FDTD方法的迭代求解公式如下：
Figure BDA0003527025810000034
Figure BDA0003527025810000035
式中，Mμ、Mε和Mσ分别是包含磁导率、介电常数和电导率等材料参数的对角矩阵，C是离散的旋度算子，E和H分别是电场强度矢量和磁场强度矢量，n是迭代时间步数，Δt为时间步长；
FDTD方法时间步长的选取受到CFL稳定性条件的限制，FL稳定性条件的表达式如下：
Figure BDA0003527025810000036
式中，Δxc为粗网格在x方向的空间步长，Δyc为粗网格在y方向的空间步长，c为电磁波在介质中的传播速度，由介电常数和磁导率决定。
本发明进一步的改进在于，步骤3中，二维TE波模式下，细网格区域麦克斯韦方程的展开式为：
Figure BDA0003527025810000041
Figure BDA0003527025810000042
Figure BDA0003527025810000043
式中，Exf、Eyf分别为细网格上的x方向和y方向的电场分量，Hzf为细网格上的z方向的磁场分量；
PITD方法在求解电磁场量时，将时间偏微分算子和空间偏微分算子分开处理；首先，对式(7)-式(9)的空间偏微分算子进行二阶精度中心差分离散，时间偏微分算子不做处理，从而得到一组常微分方程组如下：
Figure BDA0003527025810000044
Figure BDA0003527025810000045
Figure BDA0003527025810000046
式中，Δxf为细网格在x方向的空间步长，Δyf为细网格在y方向的空间步长；
将式(10)-式(12)统一写成矩阵形式为：
Figure BDA0003527025810000047
式中，X＝(Exf,Eyf,Hzf)T是一个包含空间离散网格中所有电磁场量的一维列向量，M是一个由空间步长和媒质参数所决定并且不随时间变化的系数矩阵，f(t)是由激励源引入的一维列向量；
根据常微分方程理论，并利用三点插值构造高斯积分公式，得到PITD方法的迭代公式为：
Figure BDA0003527025810000051
PITD方法的数值稳定性条件为：
Figure BDA0003527025810000052
式中，l＝2N，N为预取的正整数。
本发明进一步的改进在于，步骤4中，为了推进混合算法迭代求解的进程，对粗网格区域与细网格区域分界面上的电磁场量进行求解，以实现粗网格区域FDTD方法与细网格区域PITD方法之间计算信息的交换；其中以嵌套网格为基础，构建两种算法之间的信息交互策略。
本发明进一步的改进在于，步骤4中，具体实现方法如下
首先，假设粗网格空间步长Δxc＝Δyc＝Δc，细网格空间步长Δxf＝Δyf＝Δf，Δc/Δf＝nf；
步骤4-1，求解x方向分界面上的电磁场量；
对于x方向分界面上的场量Exc和Exf，粗网格场量的节点编号为(i+1/2，j-1/2)，对应的亚网格场量的节点编号为(I+1/2，J+1/2)，...，(I+nf-1/2，J+1/2)；
通过周围的磁场分量Hzc和Hzf，可得x方向分界面上粗网格电场分量Exc的迭代公式如下：
Figure BDA0003527025810000053
式中，
Figure BDA0003527025810000054
Figure BDA0003527025810000055
T是1×nf的单位矩阵；
x方向分界面上细网格电场分量Exf与对应的粗网格电场分量Exc相等，即迭代公式为：
Figure BDA0003527025810000056
步骤4-2，求解y方向分界面上的电磁场量；
对于y方向分界面上的场量Eyc和Eyf，粗网格场量的节点编号为(i-1/2，j+1/2)，对应的亚网格场量的节点编号为(I+1/2，J+1/2)，...，(I+1/2，J+nf-1/2)；
通过周围的磁场分量Hzc和Hzf，可得y方向分界面上粗网格电场分量Eyc的迭代公式如下：
Figure BDA0003527025810000061
式中，
Figure BDA0003527025810000062
y方向分界面上细网格电场分量Eyf的迭代公式为：
Figure BDA0003527025810000063
本发明进一步的改进在于，步骤5中，横磁波TM模式的求解过程及迭代公式可根据TE波模式类推得到，具体操作流程如下：
(1)采用局部亚网格技术对多尺度问题的计算区域进行空间网格剖分；
(2)采用FDTD方法求解粗网格区域的电磁场量Exc、Eyc及Hzc；
(3)采用PITD方法求解细网格区域的电磁场量Exf、Eyf及；
(4)利用式(16)和式(18)插值计算粗网格与细网格分界面上的Exc、Eyc分量；
(5)利用式(17)和式(19)计算粗网格与细网格分界面上的Exf、Eyf分量；
(6)重复流程(1)-流程(5)。
本发明至少具有如下有益的技术效果：
在电磁波多尺度问题中，通常情况下精细结构区域范围较小，其他区域范围较大，如果对全部区域均按照精细结构采用细网格剖分的话，网格剖分数目会很大，导致计算机内存需求也很大，带来计算资源的浪费。本发明以亚网格空间剖分方式为基础，对电磁波多尺度问题中的精细结构及附近区域采用细网格进行剖分，对其他区域采用粗网格进行剖分，并分别采用FDTD方法和PITD方法分别对粗网格区域和细网格区域的场量进行求解。与此同时，针对两种不同算法间计算信息的交换需求，以嵌套网格技术为基础构建合理的插值策略，来计算粗网格与细网格分界面上的电磁场量。该数值算法具有如下优点：
第一：采用局部亚网格技术对多尺度问题进行空间剖分，减少空间网格数目，降低计算机内存需求，避免计算资源浪费；
第二：传统的亚网格FDTD方法由于受到CFL稳定性条件的限制，时间步长的选取需要根据细网格尺寸来选择，即时间步长会很小，导致迭代步数很大，带来较大的累计误差。本发明在细网格区域采用PITD方法进行计算，该方法能够突破CFL稳定性条件的限制，即可以采用与粗网格区域FDTD方法相同的较大时间步长进行仿真。一方面，统一的时间步长能够避免在时间上进行插值或外推而引起的额外数值误差，另一方面，较大的时间步长能够减少时间迭代步数。
第三：传统PITD方法虽然不受CFL稳定性条件的限制，但是其内存需求较大，在整个区域均采用PITD方法进行计算，效率不高。本发明对范围较小的细网格区域采用能够使用较大时间步长的PITD方法进行求解，对范围较大的粗网格区域采用内存需求较小的FDTD方法进行求解，从而最大程度发挥两种算法各自的特点和优势；
第四：本发明以嵌套网格技术为基础构建了计算信息交换策略来求解不同网格分界面的电磁场量，相比于传统的插值方法，该策略能够保证混合算法在进行长时间的迭代仿真后仍然保持收敛，具有更好的长时稳定性；
第五：本发明的分界面上计算信息交换策略适用于自由空间、电介质及良导体等材料中，具有广泛的适用性；
第六：本发明所提出的算法在求解电磁波多尺度问题时具有内存需求小，迭代步数少，长时稳定性好等特点，综上使得该算法具有计算精度高、计算效率高的优点。
附图说明
图1为本发明所提出算法的计算流程示意图。
图2为空间网格剖分示意图。
图3为粗网格区域和细网格区域x方向分界面上场量Exc和Exf分布的示意图。
图4为粗网格区域和细网格区域x方向分界面上场量Eyc和Eyf分布的示意图。
图5为本发明实施例1中谐振腔内观察点处电场分量Eyc随时间的变化曲线图。
图6为本发明实施例2的几何结构及不同亚网格布置方式示意图。
图7为本发明实施例2中物块材料为铜时观察点处电场分量Eyc随时间的变化曲线图。
图8为本发明实施例2中物块材料为电介质时观察点处电场分量Eyc随时间的变化曲线图。
图9为本发明实施例3的几何结构示意图。
图10为本发明实施例3中观察点处电场分量Exc随时间的变化曲线图。
具体实施方式
下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。
本发明提出了一种应用于多尺度电磁波问题分析的高效数值仿真方法，计算流程如附图1所示。以二维TE波模式为例，首先按照附图2所示的空间网格剖分方式，对多尺度问题计算区域分别采用细网格和粗网格进行剖分，其中细网格用于精细结构及附近区域，粗网格用于其他区域。然后采用能够使用大时间步长的PITD方法求解粗网格区域电磁场量(分界面上的电磁场量除外)，采用内存需求较小的FDTD方法求解细网格区域电磁场量(分界面上的电磁场量除外)。接下来，以附图3和附图4所示的嵌套网格技术为基础，建立合理的求解策略计算分界面上细网格和粗网格的电磁场量，从而达到交换不同算法计算信息的目的。最后重复上述过程，直至完成所需迭代时间步数。TM波模式的求解过程及迭代公式可根据TE波模式类推得到。
为了完成上述目标，本发明提出了一种应用于多尺度电磁波问题分析的高效数值仿真方法，包括以下步骤：
步骤1，根据多尺度电磁波问题的几何结构，采用局部亚网格剖分技术对计算区域进行空间网格剖分，以二维TE模式为例，网格空间剖分示意图如附图2所示，精细结构部分以及附近区域使用局部细网格进行剖分，其他区域则使用粗网格进行剖分。其中，粗网格的尺寸通常取为所计算问题电磁波波长的1/10到1/8，细网格尺寸小于等于粗网格尺寸的1/2，粗网格尺寸与细网格尺寸的比值为nf。
步骤2，采用FDTD方法对粗网格区域的电磁场量进行求解。
二维TE波模式下，粗网格区域麦克斯韦方程的展开式为：
Figure BDA0003527025810000091
Figure BDA0003527025810000092
Figure BDA0003527025810000093
式中，Exc、Eyc分别为粗网格上的x方向和y方向的电场分量，Hzc为粗网格上的z方向的磁场分量，ε为计算区域介质的介电常数，σ为计算区域介质的电导率，μ为计算区域介质的磁导率。
FDTD方法在求解电磁场量时，需要同时对时间偏微分算子和空间偏微分算子进行中心差分离散，从而将偏微分方程组转化为代数方程组。FDTD方法的迭代求解公式如下：
Figure BDA0003527025810000094
Figure BDA0003527025810000095
式中，Mu、Mε和Mσ分别是包含磁导率、介电常数和电导率等材料参数的对角矩阵，c是离散的旋度算子，E和H分别是电场强度矢量和磁场强度矢量，n是迭代时间步数，Δt为时间步长。
FDTD方法时间步长的选取受到CFL稳定性条件的限制。CFL稳定性条件的表达式如下：
Figure BDA0003527025810000101
式中，Δxc为粗网格在x方向的空间步长，Δyc为粗网格在y方向的空间步长，c为电磁波在介质中的传播速度，由介电常数和磁导率决定。
步骤3，采用PITD方法对细网格区域的电磁场量进行求解。
二维TE波模式下，细网格区域麦克斯韦方程的展开式为：
Figure BDA0003527025810000102
Figure BDA0003527025810000103
Figure BDA0003527025810000104
式中，Exf、Eyf分别为细网格上的x方向和y方向的电场分量，Hzf，为细网格上的z方向的磁场分量。
PITD方法在求解电磁场量时，将时间偏微分算子和空间偏微分算子分开处理。首先，对式(26)-式(28)的空间偏微分算子进行二阶精度中心差分离散，时间偏微分算子不做处理，从而得到一组常微分方程组如下：
Figure BDA0003527025810000105
Figure BDA0003527025810000106
Figure BDA0003527025810000107
式中，Δxf为细网格在x方向的空间步长，Δyf为细网格在y方向的空间步长。
将式(29)-式(31)统一写成矩阵形式为：
Figure BDA0003527025810000111
式中，X＝(Exf，Eyf，Hzf)T是一个包含空间离散网格中所有电磁场量的一维列向量，M是一个由空间步长和媒质参数所决定并且不随时间变化的系数矩阵，f(t)是由激励源引入的一维列向量。
根据常微分方程理论，并利用三点插值构造高斯积分公式，可得PITD方法的迭代公式为：
Figure BDA0003527025810000112
PITD方法的数值稳定性条件为∶
Figure BDA0003527025810000113
式中，l＝2N，N为预取的正整数。
步骤4，粗网格区域与细网格区域分界面上电磁场量的计算。
为了推进混合算法迭代求解的进程，需要对粗网格区域与细网格区域分界面上的电磁场量进行求解，以实现粗网格区域FDTD方法与细网格区域PITD方法之间计算信息的交换。与此同时，如果计算信息交换策略不适用于混合算法，易出现不收敛、不稳定等现象。本发明以嵌套网格为基础，构建两种算法之间的信息交互策略，既能准确的交换计算信息，也能保证混合算法的稳定性。
首先，假设粗网格空间步长Δxc＝Δyc＝Δc，细网格空间步长Δxf＝Δyf＝Δf，Δc/Δf＝nf。
步骤4-1，求解x方向分界面上的电磁场量。
对于x方向分界面上的场量Exc和Exf，其嵌套网格分布如附图3所示。粗网格场量的节点编号为(i+1/2，j-1/2)，对应的亚网格场量的节点编号为(I+1/2，J+1/2)，...，(I+nf-1/2，J+1/2)。
通过周围的磁场分量Hzc和Hzf，可得x方向分界面上粗网格电场分量Exc的迭代公式如下：
Figure BDA0003527025810000121
式中，
Figure BDA0003527025810000122
Figure BDA0003527025810000123
T是1×nf的单位矩阵。
x方向分界面上细网格电场分量Exf与对应的粗网格电场分量Exc相等，即迭代公式为：
Figure BDA0003527025810000126
步骤4-2，求解y方向分界面上的电磁场量。
对于y方向分界面上的场量Eyc和Eyf，其嵌套网格分布如附图4所示。粗网格场量的节点编号为(i-1/2，j+1/2)，对应的亚网格场量的节点编号为(I+1/2，J+1/2)，...，(I+1/2，J+nf-1/2)。
通过周围的磁场分量Hzc和Hzf，可得y方向分界面上粗网格电场分量Eyc的迭代公式如下：
Figure BDA0003527025810000124
式中，
Figure BDA0003527025810000125
y方向分界面上细网格电场分量Eyf的迭代公式为：
Figure BDA0003527025810000127
步骤5，重复上述步骤，完成所有迭代步数。TM波模式的求解过程及迭代公式可根据TE波模式类推得到。
具体操作流程如下：
(1)采用局部亚网格技术对多尺度问题的计算区域进行空间网格剖分。
(2)采用FDTD方法求解粗网格区域的电磁场量Exc、Eyc及Hzc(不包括粗网格与细网格分界面上的Exc、Eyc分量)；
(3)采用PITD方法求解细网格区域的电磁场量Exf、Eyf及Hzf(不包括粗网格与细网格分界面上的Exf、Eyf分量)；
(4)利用式(35)和式(37)插值计算粗网格与细网格分界面上的Exc、Eyc分量；
(5)利用式(36)和式(38)计算粗网格与细网格分界面上的Exf、Eyf分量；
(6)重复流程(1)-流程(5)。
实施例1：尺寸为60mm×40mm的二维谐振腔。
中心10mm×10mm区域采用局部亚网格进行剖分，其他区域采用粗网格进行剖分。粗网格的尺寸为Δxc＝Δyc＝1mm，粗网格尺寸与细网格尺寸的比例为nf＝4。以FDTD方法应满足的CFL稳定性为依据，仿真时间步长Δt为按照粗网格尺寸所能选取的最大时间步长。整个仿真的迭代时间步数为106，激励源为高斯微分脉冲。
附图5所示为本发明提出的混合方法在计算实施例1时，观察点处电场强度Ey随时间的变化曲线图。经过106个时间迭代步之后，计算结果仍然收敛，即算法具有良好的长时稳定性。
实施例2：如附图6所示，计算区域为40mm×40mm，内部设置一个8mm×8mm的物块，物块的材料按照两种情形设置，分别为铜(电导率为5.6×107s/m)和有耗电介质(电导率为5S/m，相对介电常数为2)。亚网格的配置形式包括包围物块、穿过物块和不包含物块三种。
分别采用全粗网格FDTD方法、全细网格FDTD方法以及三种亚网格配置形式下的混合算法对实施例2进行仿真。粗网格的尺寸为Δxc＝Δyc＝1mm，粗网格尺寸与细网格尺寸的比例为nf＝4。全细网格FDTD方法的时间步长设置为Δt＝0.589ps，全粗网格FDTD方法和混合算法的时间步长设置为4Δt。Δt为细网格尺寸下根据CFL稳定性条件所能选取的最大时间步长。激励源为高斯微分脉冲。整个计算区域采用拉伸坐标系完美匹配层(SC-PML)吸收边界条件进行截断。
附图7所示为采用不同算法在计算实施例2(物块为铜)时，观察点处电场强度Ey随时间的变化曲线图。
附图8所示为采用不同算法在计算实施例2(物块为有耗电介质)时，观察点处电场强度Ey随时间的变化曲线图。
在不同的亚网格布置形式以及不同的分界面材料属性情况下，本发明所提出算法的计算结果与其他算法保持了一致性，具有良好的可行性。
实施例3：探地天线测试模拟，几何结构布置如附图9所示。电介质的材料参数为相对介电常数εr＝6，电导率σ＝0。电介质内设置一条3mm宽的空气气隙。探地雷达发射天线(T点)的激励源设置为中心频率为900MHz的高斯微分脉冲。
分别采用传统全细网格FDTD方法、亚网格FDTD方法及本发明的混合算法对该算例进行模拟。对于亚网格算法，对空气气隙及附近区域采用亚网格进行剖分，亚网格的尺寸设置为Δxf＝Δyf＝1mm，其他区域采用粗网格进行剖分，粗网格的尺寸设置为Δxc＝Δyc＝3mm，即nf＝3。传统FDTD方法按照统一网格尺寸进行剖分，即Δx＝Δy＝1mm。传统全细网格FDTD方法和亚网格FDTD方法的时间步长根据细网格的尺寸来设置，即Δt＝2.35ps，本发明所提出算法的时间步长设置为3Δt。
附图10所示为采用不同算法在计算实施例3时，接收天线(O点)处电场强度Ex随时间的变化曲线图。可以观察到，电介质与空气分界面、空气气隙的两次反射情况，本发明所提出算法的计算结果与其他算法保持了一致性。
附表1所示为采用不同算法计算实施例3的计算资源信息的对比。
表1
Figure BDA0003527025810000141
Figure BDA0003527025810000151
传统FDTD方法采用了统一的细网格进行剖分，主计算区域的网格数目约为360000，SC-PML吸收边界条件的网格剖分数目约为24400，网格剖分数目较大。采用局部亚网格技术进行剖分，主计算区域的网格数目约为40240，SC-PML吸收边界条件的网格剖分数目约为8400，网格数目大大降低。
本发明所提出的混合算法能够使用较大的时间步长进行仿真，迭代步数为其他两种算法的1/3，CPU执行时间最短，约为传统FDTD方法的1/22，亚网格FDTD方法的1/2。
实施例3验证了本发明所提出混合算法能够有效降低计算机内存需求以及CPU执行时间，从而提高在求解电磁波多尺度问题时的计算效率。
虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (8)

1.一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，包括：
步骤1，根据多尺度电磁波问题的几何结构，采用局部亚网格剖分技术对计算区域进行空间网格剖分；
步骤2，采用FDTD方法对粗网格区域的电磁场量进行求解；
步骤3，采用PITD方法对细网格区域的电磁场量进行求解；
步骤4，粗网格区域与细网格区域分界面上电磁场量的计算；
步骤5，重复上述步骤，完成所有迭代步数。
2.根据权利要求1所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤1中，对于二维横电波TE模式，精细结构部分以及附近区域使用局部细网格进行剖分，其他区域则使用粗网格进行剖分。
3.根据权利要求2所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，粗网格的尺寸为所计算问题电磁波波长的1/10到1/8，细网格尺寸小于等于粗网格尺寸的1/2，粗网格尺寸与细网格尺寸的比值为nf。
4.根据权利要求3所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤2中，二维TE波模式下，粗网格区域麦克斯韦方程的展开式为：
Figure FDA0003527025800000011
Figure FDA0003527025800000012
Figure FDA0003527025800000013
式中，Exc、Eyc分别为粗网格上的x方向和y方向的电场分量，Hzc为粗网格上的z方向的磁场分量，ε为计算区域介质的介电常数，σ为计算区域介质的电导率，μ为计算区域介质的磁导率；
FDTD方法在求解电磁场量时，同时对时间偏微分算子和空间偏微分算子进行中心差分离散，从而将偏微分方程组转化为代数方程组，FDTD方法的迭代求解公式如下：
Figure FDA0003527025800000021
Figure FDA0003527025800000022
式中，Mμ、Mε和Mσ分别是包含磁导率、介电常数和电导率等材料参数的对角矩阵，C是离散的旋度算子，E和H分别是电场强度矢量和磁场强度矢量，n是迭代时间步数，Δt为时间步长；
FDTD方法时间步长的选取受到CFL稳定性条件的限制，FL稳定性条件的表达式如下：
Figure FDA0003527025800000023
式中，Δxc为粗网格在x方向的空间步长，Δyc为粗网格在y方向的空间步长，c为电磁波在介质中的传播速度，由介电常数和磁导率决定。
5.根据权利要求4所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤3中，二维TE波模式下，细网格区域麦克斯韦方程的展开式为：
Figure FDA0003527025800000024
Figure FDA0003527025800000025
Figure FDA0003527025800000026
式中，Exf、Eyf分别为细网格上的x方向和y方向的电场分量，Hzf为细网格上的z方向的磁场分量；
PITD方法在求解电磁场量时，将时间偏微分算子和空间偏微分算子分开处理；首先，对式(7)-式(9)的空间偏微分算子进行二阶精度中心差分离散，时间偏微分算子不做处理，从而得到一组常微分方程组如下：
Figure FDA0003527025800000027
Figure FDA0003527025800000028
Figure FDA0003527025800000031
式中，Δxf为细网格在x方向的空间步长，Δyf为细网格在y方向的空间步长；
将式(10)-式(12)统一写成矩阵形式为：
Figure FDA0003527025800000032
式中，X＝(Exf,Eyf,Hzf)T是一个包含空间离散网格中所有电磁场量的一维列向量，M是一个由空间步长和媒质参数所决定并且不随时间变化的系数矩阵，f(t)是由激励源引入的一维列向量；
根据常微分方程理论，并利用三点插值构造高斯积分公式，得到PITD方法的迭代公式为：
Figure FDA0003527025800000033
PITD方法的数值稳定性条件为：
Figure FDA0003527025800000034
式中，l＝2N，N为预取的正整数。
6.根据权利要求5所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤4中，为了推进混合算法迭代求解的进程，对粗网格区域与细网格区域分界面上的电磁场量进行求解，以实现粗网格区域FDTD方法与细网格区域PITD方法之间计算信息的交换；其中以嵌套网格为基础，构建两种算法之间的信息交互策略。
7.根据权利要求6所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤4中，具体实现方法如下
首先，假设粗网格空间步长Δxc＝Δyc＝Δc，细网格空间步长Δxf＝Δyf＝Δf，Δc/Δf＝nf；
步骤4-1，求解x方向分界面上的电磁场量；
对于x方向分界面上的场量Exc和Exf，粗网格场量的节点编号为(i+1/2,j-1/2)，对应的亚网格场量的节点编号为(I+1/2,J+1/2)，…，(I+nf-1/2,J+1/2)；
通过周围的磁场分量Hzc和Hzf，可得x方向分界面上粗网格电场分量Exc的迭代公式如下：
Figure FDA0003527025800000041
式中，
Figure FDA0003527025800000042
Figure FDA0003527025800000043
T是1×nf的单位矩阵；
x方向分界面上细网格电场分量Exf与对应的粗网格电场分量Exc相等，即迭代公式为：
Exf|I+1/2,J＝…＝Exf|I+nf-1/2,J＝Exc|i+1/2,j (17)
步骤4-2，求解y方向分界面上的电磁场量；
对于y方向分界面上的场量Eyc和Eyf，粗网格场量的节点编号为(i-1/2,j+1/2)，对应的亚网格场量的节点编号为(I+1/2,J+1/2)，…，(I+1/2,J+nf-1/2)；
通过周围的磁场分量Hzc和Hzf，可得y方向分界面上粗网格电场分量Eyc的迭代公式如下：
Figure FDA0003527025800000044
式中，
Figure FDA0003527025800000045
y方向分界面上细网格电场分量Eyf的迭代公式为：
Eyf|I,J+1/2＝…＝Eyf|I,J+nf-1/2＝Eyc|i,j+1/2 (19)。
8.根据权利要求7所述的一种应用于多尺度电磁波问题分析的高效数值仿真方法，其特征在于，步骤5中，横磁波TM模式的求解过程及迭代公式可根据TE波模式类推得到，具体操作流程如下：
(1)采用局部亚网格技术对多尺度问题的计算区域进行空间网格剖分；
(2)采用FDTD方法求解粗网格区域的电磁场量Exc、Eyc及Hzc；
(3)采用PITD方法求解细网格区域的电磁场量Exf、Eyf及；
(4)利用式(16)和式(18)插值计算粗网格与细网格分界面上的Exc、Eyc分量；
(5)利用式(17)和式(19)计算粗网格与细网格分界面上的Exf、Eyf分量；
(6)重复流程(1)-流程(5)。
CN202210199606.7A 2022-03-01 2022-03-01 一种应用于多尺度电磁波问题分析的高效数值仿真方法 Pending CN114528716A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210199606.7A CN114528716A (zh) 2022-03-01 2022-03-01 一种应用于多尺度电磁波问题分析的高效数值仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210199606.7A CN114528716A (zh) 2022-03-01 2022-03-01 一种应用于多尺度电磁波问题分析的高效数值仿真方法

Publications (1)

Publication Number Publication Date
CN114528716A true CN114528716A (zh) 2022-05-24

Family

ID=81626159

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210199606.7A Pending CN114528716A (zh) 2022-03-01 2022-03-01 一种应用于多尺度电磁波问题分析的高效数值仿真方法

Country Status (1)

Country Link
CN (1) CN114528716A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115542315A (zh) * 2022-09-23 2022-12-30 安徽大学 基于adi-fdtd的探地雷达正演模拟方法及系统
CN117195650A (zh) * 2023-09-19 2023-12-08 安徽大学 基于高阶矩阵指数完美匹配层的fdtd计算方法及系统
CN117390935A (zh) * 2023-12-11 2024-01-12 芯瑞微(上海)电子科技有限公司 一种计算fdtd电磁仿真收敛检测触发时刻的算法
CN117452081A (zh) * 2023-12-26 2024-01-26 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115542315A (zh) * 2022-09-23 2022-12-30 安徽大学 基于adi-fdtd的探地雷达正演模拟方法及系统
CN115542315B (zh) * 2022-09-23 2023-08-18 安徽大学 基于adi-fdtd的探地雷达正演模拟方法及系统
CN117195650A (zh) * 2023-09-19 2023-12-08 安徽大学 基于高阶矩阵指数完美匹配层的fdtd计算方法及系统
CN117195650B (zh) * 2023-09-19 2024-04-05 安徽大学 基于高阶矩阵指数完美匹配层的fdtd计算方法及系统
CN117390935A (zh) * 2023-12-11 2024-01-12 芯瑞微(上海)电子科技有限公司 一种计算fdtd电磁仿真收敛检测触发时刻的算法
CN117390935B (zh) * 2023-12-11 2024-03-01 芯瑞微(上海)电子科技有限公司 一种计算fdtd电磁仿真收敛检测触发时刻的方法
CN117452081A (zh) * 2023-12-26 2024-01-26 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端
CN117452081B (zh) * 2023-12-26 2024-04-30 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端

Similar Documents

Publication Publication Date Title
CN114528716A (zh) 一种应用于多尺度电磁波问题分析的高效数值仿真方法
Notaros Higher order frequency-domain computational electromagnetics
CN108959772B (zh) 一种大规模有限周期阵列结构特征模式分析方法
Nair et al. Generalized method of moments: A novel discretization technique for integral equations
Zhang et al. Performance of a massively parallel higher-order method of moments code using thousands of CPUs and its applications
Cai et al. Volume surface integral equation method based on higher order hierarchical vector basis functions for EM scattering and radiation from composite metallic and dielectric structures
Pan et al. Hierarchical interpolative decomposition multilevel fast multipole algorithm for dynamic electromagnetic simulations
CN114519287A (zh) 电大多尺度复杂目标的三维电磁场求解方法
Álvarez González A discontinuous Galerkin finite element method for the time-domain solution of Maxwell equations
Killian et al. Electromagnetic scattering from electrically large arbitrarily-shaped conductors using the method of moments and a new null-field generation technique
Zhao et al. A hybrid 2-D/3-D multilevel Green’s function interpolation method for electrically large multilayered problems
Ding et al. A novel hierarchical two-level spectral preconditioning technique for electromagnetic wave scattering
Carpinteri et al. The partition of unity quadrature in element-free crack modelling
Amor-Martin et al. Study of accuracy of a non-conformal finite element domain decomposition method
Garza et al. High-order Chebyshev-based Nyström methods for electromagnetics
Shi et al. An OpenMP parallelized multilevel Green's function interpolation method accelerated by fast Fourier transform technique
Koziel et al. On improved-reliability design optimization of high-frequency structures using local search algorithms
Hussein Efficient near-field computation for radiation and scattering from conducting surfaces of arbitrary shape
Khalef et al. An unconditionally stable radial point interpolation meshless method based on the Crank–Nicolson scheme solution of wave equation
Xin et al. An Efficient and Effective Preconditioner for the Discontinuous Galerkin Domain Decomposition Method of Surface Integral Equation
Kwon et al. Impedance matrix generation by using the fast matrix generation technique
Kaklamani et al. Extension of method of moments for electrically large structures based on parallel computations
Önol et al. Efficient three-layer iterative solutions of electromagnetic problems using the multilevel fast multipole algorithm
JP3927485B2 (ja) 電磁波解析装置及び方法
Lee et al. Numerical methods in antenna modeling

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination