CN113158492B - 一种时变电磁场的全隐式双时间步计算方法 - Google Patents

一种时变电磁场的全隐式双时间步计算方法 Download PDF

Info

Publication number
CN113158492B
CN113158492B CN202110527205.5A CN202110527205A CN113158492B CN 113158492 B CN113158492 B CN 113158492B CN 202110527205 A CN202110527205 A CN 202110527205A CN 113158492 B CN113158492 B CN 113158492B
Authority
CN
China
Prior art keywords
time
grid
electromagnetic
iteration
matrix
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.)
Active
Application number
CN202110527205.5A
Other languages
English (en)
Other versions
CN113158492A (zh
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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202110527205.5A priority Critical patent/CN113158492B/zh
Publication of CN113158492A publication Critical patent/CN113158492A/zh
Application granted granted Critical
Publication of CN113158492B publication Critical patent/CN113158492B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/18—Network design, e.g. design based on topological or interconnect aspects of utility systems, piping, heating ventilation air conditioning [HVAC] or cabling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/12—Timing analysis or timing optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种时变电磁场的全隐式双时间步计算方法，其包括根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；数值计算区域的四边形或六面体结构网格，在壁面和几何奇异处加密，所述网格逐渐远离散射壁面而逐渐稀疏；输出网格数据文件和设定和输出边界条件文件；输入目标计算电磁参数、数值计算控制参数；输入网格数据和边界条件信息文件，初始化计算空间电磁场；以基于时间迭代推进和空间通量残差的隐式双时间步方式对麦克斯韦方程组时变电磁场进行迭代求解。通过在控制方程中引入定常的虚拟时间导数项，从而使得物理时间推进步长可以根据物理问题进行选取而不受稳定性的限制，在保持高数值精度的同时，提升计算性能。

Description

一种时变电磁场的全隐式双时间步计算方法
技术领域
本发明涉及电磁学的时域数值求解技术领域，尤其是涉及一种可提高时变电磁场时间推进效率的全隐式双时间步计算方法。
背景技术
复杂外形目标电磁散射、复杂电磁环境电磁干扰，应用时域方法能能相容地模拟散射、多重散射、孔穿透、腔激励等复杂现象，并且能准确模拟时间历程更直观，不像传统高频渐进方法对特殊部件、特殊电磁现象例如棱边衍射提出特殊处理方式。
时域中的时变电磁场满足时域麦克斯韦方程组，随着计算机技术的发展，直接求解该方程组成为可能。与欧拉方程相同的双曲型数学特征促进计算流体力学（Computational Fluid Dynamics, CFD）技术在电磁场计算中的应用，其中时域有限差分法(Finite Difference Time Domain, FDTD)和时域有限体积法(Finite Volume TimeDomain, FVTD)最为著名。20世纪60年代K.S.Yee发表先驱性的时域有限差分算法，直接差分计算时变麦克斯韦微分方程组，成功地模拟电磁脉冲与理想导电体作用的时域响应，开创一种新的电磁场时域计算方法。在Yee算法中，首先在感兴趣区域（目标及其周围一定空间）生成笛卡尔直角正交网格，电场和磁场各分量在网格空间的取值点被交叉放置，使得在每个坐标平面上每个电场分量的四周由磁场分量环绕，同时每个磁场分量的四周由电场分量所环绕，这样的电磁场配置符合法拉第感应定律和安培环路定律的要求，这种网格通常称为Yee氏网格。
传统时域有限差分法和时域有限体积法时间推进采用2阶中心差分或时空耦合Lax-Wendroff格式以及多步Runge-Kutta法，其共同点是时间计算的显式格式。以Runge-Kutta方法为代表的显式方法既有编程简便也有易实现时间高精度的优点，对时域电磁场计算是可靠的时间离散方法。但是时间显式方法有一个最大缺陷，其时间步长受稳定性的限制, 整个计算空间必须采用统一的最小全局计算步长,为模拟几何外形剧烈变化生成的贴体加密网格（例如机翼前后缘棱边几何奇异性带来电磁场梯度剧烈变化要求加密网格仔细模拟以及电磁多尺度问题），会带来很小的全局时间步长，大的网格单元就需要更多时间步在该单元传播信息，从而使得到稳定的时变电磁场需要较长的计算时间, 特别是在求解高频、电大尺寸目标电磁散射时域问题时, 受稳定性限制的小时间步长带来时域电磁场计算量的显著增加，消耗大量计算资源。另一方面，隐式计算方法能够放宽计算步长的稳定性限制, 但伴随而来是时间精度下降和加密导致系数矩阵维度增加提高矩阵求逆运算难度。
综上而言，迫切需要一种高效的直接求解时域麦克斯韦方程组，在保持高数值精度的同时，放宽很小网格尺度对迭代物理时间步的限制，提升计算性能。
发明内容
本发明提供一种基于时间迭代和空间通量残差隐式的时变电磁场的全隐式双时间步计算方法，其物理时间步长根据物理问题进行选取而不受稳定限制，稳定性由隐式虚拟时间子迭代满足，从而可克服现有技术中显示时间步计算方法受稳定性限制且必须采用统一最小全局步长和加密网络进行计算带来计算量大的问题，大大提高了计算效率。
为实现上述目的，本发明提供如下技术方案：一种时变电磁场的全隐式双时间步计算方法，包括以下步骤：
步骤1，根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；
步骤2，采用四边形或六面体结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件；
步骤3，输入目标计算电磁参数、数值计算控制参数；
步骤4，输入网格数据和边界条件信息文件，初始化计算空间电磁场；
步骤5，以基于时间迭代推进和空间通量残差的隐式双时间步方式对麦克斯韦方程组时变电磁场进行迭代求解。
优选地，所述步骤5过程如下：仿真模型的外层为物理时间步循环，直至计算收敛结束；仿真模型的内层为虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格、各个网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
优选地，所述物理时间步循环和虚拟时间步子迭代循环过程为：
通过虚拟时间导数项
Figure 942632DEST_PATH_IMAGE001
，将待求解的麦克斯韦方程组修正为：
Figure 440609DEST_PATH_IMAGE002
其中，
Figure 149939DEST_PATH_IMAGE003
是实数型的电磁场守恒变量，
Figure 100578DEST_PATH_IMAGE004
是归一化物理时间，
Figure 983083DEST_PATH_IMAGE005
是直角坐标系下电磁通量的
Figure 522518DEST_PATH_IMAGE006
分量；
Figure 465066DEST_PATH_IMAGE007
是实数型磁感应强度矢量，
Figure 840684DEST_PATH_IMAGE008
是实数型电位移矢量，
Figure 261301DEST_PATH_IMAGE009
是电场强度矢量，
Figure 422286DEST_PATH_IMAGE010
是磁场强度矢量，
Figure 535735DEST_PATH_IMAGE011
、
Figure 398649DEST_PATH_IMAGE012
、
Figure 622957DEST_PATH_IMAGE013
分别是
Figure 871405DEST_PATH_IMAGE007
的
Figure 155756DEST_PATH_IMAGE014
分量；
Figure 505966DEST_PATH_IMAGE015
、
Figure 533964DEST_PATH_IMAGE016
、
Figure 715547DEST_PATH_IMAGE017
分别是
Figure 590706DEST_PATH_IMAGE008
的
Figure 490529DEST_PATH_IMAGE006
分量；
Figure 994322DEST_PATH_IMAGE018
、
Figure 295991DEST_PATH_IMAGE019
、
Figure 843516DEST_PATH_IMAGE020
分别是
Figure 230635DEST_PATH_IMAGE009
的
Figure 538119DEST_PATH_IMAGE014
分量；
Figure 428715DEST_PATH_IMAGE021
、
Figure 225769DEST_PATH_IMAGE022
、
Figure 523021DEST_PATH_IMAGE023
分别是
Figure 696513DEST_PATH_IMAGE010
的
Figure 644878DEST_PATH_IMAGE006
分量；当
Figure 612834DEST_PATH_IMAGE024
收敛时，该方程组等同于原始方程组，定常虚拟时间
Figure 895916DEST_PATH_IMAGE025
子迭代表示为：
Figure 607520DEST_PATH_IMAGE026
其中，
Figure 410391DEST_PATH_IMAGE027
是
Figure 549249DEST_PATH_IMAGE003
经虚拟时间子迭代后的电磁场守恒变量值，
Figure 132677DEST_PATH_IMAGE027
是
Figure 67878DEST_PATH_IMAGE028
的近似；
Figure 53152DEST_PATH_IMAGE029
是通量残差，
Figure 300593DEST_PATH_IMAGE030
是通量残差
Figure 371318DEST_PATH_IMAGE029
增加了物理时间导数项后的残差；
Figure 424724DEST_PATH_IMAGE031
是物理时间步数；
Figure 451455DEST_PATH_IMAGE032
是物理时间步长，
Figure 932115DEST_PATH_IMAGE033
是第
Figure 162239DEST_PATH_IMAGE031
物理时间步时的电磁守恒变量，
Figure 19337DEST_PATH_IMAGE034
是第
Figure 402039DEST_PATH_IMAGE035
物理时间步时的电磁守恒变量；
Figure 53600DEST_PATH_IMAGE036
为
Figure 98916DEST_PATH_IMAGE027
对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度；定常子迭代部分采用隐式算法：
Figure 431809DEST_PATH_IMAGE037
其中，
Figure 980602DEST_PATH_IMAGE038
为结构网格曲线坐标系方向1，
Figure 255594DEST_PATH_IMAGE039
为结构网格曲线坐标系方向2，
Figure 522627DEST_PATH_IMAGE040
为结构网格曲线坐标系方向3；
Figure 659211DEST_PATH_IMAGE041
分别对应为
Figure 62510DEST_PATH_IMAGE038
、
Figure 321453DEST_PATH_IMAGE039
、
Figure 784706DEST_PATH_IMAGE040
方向的电磁通量；
Figure 521718DEST_PATH_IMAGE042
是隐式控制参数，取
Figure 982786DEST_PATH_IMAGE043
的全隐式，其他参数对应显、隐混合格式；下标
Figure 412630DEST_PATH_IMAGE044
是网格单元编号，
Figure 575627DEST_PATH_IMAGE045
是第
Figure 116330DEST_PATH_IMAGE044
网格单元第
Figure 697484DEST_PATH_IMAGE046
虚拟时间迭代步时的电磁守恒变量，
Figure 32651DEST_PATH_IMAGE047
是第
Figure 449988DEST_PATH_IMAGE044
网格单元第
Figure 528802DEST_PATH_IMAGE048
虚拟时间迭代步时的电磁守恒变量，
Figure 964463DEST_PATH_IMAGE049
是第
Figure 470530DEST_PATH_IMAGE050
物理时间步时的电磁守恒变量，
Figure 873699DEST_PATH_IMAGE033
是第
Figure 756204DEST_PATH_IMAGE031
物理时间步时的电磁守恒变量，
Figure 46371DEST_PATH_IMAGE034
是第
Figure 988919DEST_PATH_IMAGE035
物理时间步时的电磁守恒变量，
Figure 112340DEST_PATH_IMAGE051
是第
Figure 532957DEST_PATH_IMAGE044
网格单元第
Figure 5526DEST_PATH_IMAGE052
虚拟时间迭代步时的空间通量残差，
Figure 791080DEST_PATH_IMAGE053
是第
Figure 716311DEST_PATH_IMAGE044
网格单元第
Figure 127569DEST_PATH_IMAGE031
虚拟时间迭代步时的空间通量残差；
Figure 454645DEST_PATH_IMAGE054
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算；采用不同局部虚拟时间子迭代步长定常计算不同的网格单元，从而加快相应网格单元的单元电磁场收敛。
优选地，所述空间通量计算和隐式迭代解计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量
Figure 676679DEST_PATH_IMAGE055
Figure 777622DEST_PATH_IMAGE056
Figure DEST_PATH_IMAGE057
；
Figure 477724DEST_PATH_IMAGE058
；
式中
Figure DEST_PATH_IMAGE059
分别取曲线坐标系
Figure 846258DEST_PATH_IMAGE060
方向之一，相应的
Figure 301510DEST_PATH_IMAGE061
即为对应
Figure 139016DEST_PATH_IMAGE062
方向的电磁通量；
Figure 705126DEST_PATH_IMAGE063
代表曲线坐标系
Figure 426701DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 52855DEST_PATH_IMAGE064
代表曲线坐标系
Figure 377657DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 481879DEST_PATH_IMAGE065
为相似矩阵，
Figure 638054DEST_PATH_IMAGE066
分别为正负特征值构成的对角矩阵，
Figure 622059DEST_PATH_IMAGE067
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 230895DEST_PATH_IMAGE068
代表自变量为
Figure 342071DEST_PATH_IMAGE069
的
Figure DEST_PATH_IMAGE070
相似矩阵；
Figure 41167DEST_PATH_IMAGE071
代表自变量为
Figure 9123DEST_PATH_IMAGE069
的
Figure DEST_PATH_IMAGE072
对角矩阵；
Figure 42939DEST_PATH_IMAGE073
代表自变量为
Figure 488963DEST_PATH_IMAGE069
的
Figure DEST_PATH_IMAGE074
相似矩阵；
Figure 806681DEST_PATH_IMAGE075
代表自变量为
Figure DEST_PATH_IMAGE076
的
Figure 883222DEST_PATH_IMAGE070
相似矩阵；
Figure 466650DEST_PATH_IMAGE077
代表自变量为
Figure 401851DEST_PATH_IMAGE076
的
Figure DEST_PATH_IMAGE078
对角矩阵；
Figure 324808DEST_PATH_IMAGE079
代表自变量为
Figure 634566DEST_PATH_IMAGE076
的
Figure 705291DEST_PATH_IMAGE074
相似矩阵；
Figure DEST_PATH_IMAGE080
；
Figure 945648DEST_PATH_IMAGE081
；
其中
Figure DEST_PATH_IMAGE082
是限制器，下标
Figure 723111DEST_PATH_IMAGE083
是网格单元编号，
Figure DEST_PATH_IMAGE084
对应单元分界面，
Figure 892187DEST_PATH_IMAGE085
是3阶精度格式的控制参数，
Figure DEST_PATH_IMAGE086
和
Figure 122311DEST_PATH_IMAGE087
分别是后差和前差算符；
Figure DEST_PATH_IMAGE088
表示网格单元
Figure 166359DEST_PATH_IMAGE089
分界面处左状态电磁守恒变量，
Figure DEST_PATH_IMAGE090
表示网格单元
Figure 860645DEST_PATH_IMAGE089
分界面处右状态电磁守恒变量；
Figure 449890DEST_PATH_IMAGE091
是第
Figure 495206DEST_PATH_IMAGE083
个网格单元电磁场守恒变量，
Figure DEST_PATH_IMAGE092
是第
Figure 575901DEST_PATH_IMAGE083
+1个网格单元电磁场守恒变量；
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure 62377DEST_PATH_IMAGE093
其中，
Figure 150419DEST_PATH_IMAGE094
是分裂后的系数矩阵，
Figure 417452DEST_PATH_IMAGE095
是上一个物理时间步计算的空间通量残差，
Figure 537724DEST_PATH_IMAGE096
是隐式子迭代电磁场差值；
将
Figure 206603DEST_PATH_IMAGE093
表示为LDU，近似因式分解得：
Figure 403229DEST_PATH_IMAGE097
Figure 891979DEST_PATH_IMAGE098
Figure 894570DEST_PATH_IMAGE099
Figure 106371DEST_PATH_IMAGE100
其中，下标
Figure 270636DEST_PATH_IMAGE101
是网格单元编号，
Figure 449945DEST_PATH_IMAGE102
是雅可比系数矩阵最大特征值分裂参数，
Figure 990647DEST_PATH_IMAGE103
是雅克比系数矩阵最大特征值；
Figure 555490DEST_PATH_IMAGE104
是单位对角矩阵，
Figure 156235DEST_PATH_IMAGE105
是对角矩阵，
Figure 619578DEST_PATH_IMAGE106
为上三角矩阵，
Figure 901654DEST_PATH_IMAGE107
为下三角矩阵，
Figure 399632DEST_PATH_IMAGE108
是对应上三角矩阵的电磁守恒变量差值，
Figure 591185DEST_PATH_IMAGE109
是对应下三角矩阵的电磁守恒变量差值；
Figure 807403DEST_PATH_IMAGE110
指
Figure 627592DEST_PATH_IMAGE111
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 980076DEST_PATH_IMAGE112
指
Figure 657045DEST_PATH_IMAGE113
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 281930DEST_PATH_IMAGE114
指
Figure 968126DEST_PATH_IMAGE115
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 112800DEST_PATH_IMAGE116
指
Figure 226249DEST_PATH_IMAGE117
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 839895DEST_PATH_IMAGE118
指
Figure 64203DEST_PATH_IMAGE119
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 391279DEST_PATH_IMAGE120
指
Figure 613313DEST_PATH_IMAGE121
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 760261DEST_PATH_IMAGE122
是指相邻
Figure 975210DEST_PATH_IMAGE123
网格单元分裂后的系数矩阵；
Figure 156793DEST_PATH_IMAGE124
是指相邻
Figure 549728DEST_PATH_IMAGE125
网格单元分裂后的系数矩阵；
Figure 449551DEST_PATH_IMAGE126
是指相邻
Figure 701148DEST_PATH_IMAGE127
网格单元分裂后的系数矩阵；
Figure 737237DEST_PATH_IMAGE128
是指相邻
Figure 363390DEST_PATH_IMAGE129
网格单元分裂后的系数矩阵；
Figure 422613DEST_PATH_IMAGE130
是指相邻
Figure 792414DEST_PATH_IMAGE131
网格单元分裂后的系数矩阵；
Figure 135540DEST_PATH_IMAGE132
是指相邻
Figure 932595DEST_PATH_IMAGE133
网格单元分裂后的系数矩阵；
最后，得到前后向迭代计算电磁场差值
Figure 479114DEST_PATH_IMAGE134
：
Figure 387027DEST_PATH_IMAGE135
Figure 86124DEST_PATH_IMAGE136
Figure 54080DEST_PATH_IMAGE137
其中，
Figure 150212DEST_PATH_IMAGE138
是对角矩阵
Figure 799499DEST_PATH_IMAGE139
的逆矩阵，
Figure 930266DEST_PATH_IMAGE140
、
Figure 256074DEST_PATH_IMAGE141
分别是根据
Figure 573923DEST_PATH_IMAGE142
、
Figure 26901DEST_PATH_IMAGE143
计算得出的上三角矩阵、下三角矩阵；
前向循环：
Figure 12174DEST_PATH_IMAGE144
后向循环：
Figure 321933DEST_PATH_IMAGE145
其中，
Figure 812564DEST_PATH_IMAGE146
是电磁守恒变量差值的中间过渡变量。
优选地，所述步骤2中：网格密度保证每波长13-20个网格点，壁面密度>300点/波长，几何奇异处加密到50-100个网格点/波长；二维网格时，则在垂直该二维网格所在平面按右手法则推进1层，作为三维问题特例统一计算；网格数据文件包括结构网格块数目及每块3个曲线坐标系下维度。
优选地，所述步骤3中：如有等离子体外部流场情况，则还应输入对应流场参数。
优选地，所述步骤3中，预先设定物理时间步长，二维问题用入射电磁波周期无量纲化的物理时间步长
Figure 865970DEST_PATH_IMAGE147
计算精度，三维问题选择
Figure 643434DEST_PATH_IMAGE148
量级；同时设定子迭代收敛标准判据值为全网格空间相邻子迭代时间步最大幅度差值<0.001，判定为收敛；设定最大子迭代步数为50。
优选地，所述步骤5中，子迭代收敛判据为电场和磁场在计算空间最大幅度差的绝对值，且数值计算中设置有最大子迭代步数。
与现有技术相比，本发明的有益效果是：
1、本发明通过一种新的两重（物理和虚拟）时间迭代推进和空间通量残差隐式计算的时域有限体积方法，来替代现有方法中通常采用的时间、空间通量残差都为显式的求解算法。如此一来，使得物理时间步长根据物理问题进行选取而不受稳定性的限制，稳定性则由隐式虚拟时间子迭代满足，从而克服显式方法时间步受稳定性限制且必须采用统一最小全局步长和加密网格带来的大计算量问题，提高计算效率。最后，利用这种双时间步隐式时域有限体积方法高效地求解2维、3维麦克斯韦方程组，获得时变电磁场的时间、空间分布，在保证格式精度的同时提高时间推进效率，节约计算成本。
2、时间推进的双时间步迭代和空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解。时间双时间步和空间通量隐式相结合的时域有限体积方法能稳定、高精度、高效地获得电磁场时间、空间分布，适用于数值模拟电磁波传播、反射等现象，特别适合电大尺寸目标大规模、大运算量电磁场计算和散射隐身特性计算。
3、本发明支持结构网格和多区域分解，采用贴体网格直接可将麦克斯韦方程组守恒律的积分形式应用到离散的曲线坐标系网格单元。
4、本发明在保持高数值精度的同时，放宽很小网格尺度对迭代物理时间步的限制，提升计算性能。
5、本发明用以计算时变电磁场和相应目标电磁特性,通过在控制方程中引入定常的虚拟时间导数项，在每个物理时间段内对物理时间导数作线性化处理，从而使得物理时间推进步长可以根据物理问题进行选取而不受稳定性的限制，计算的稳定性要求由虚拟时间子迭代满足，子迭代的定常计算可以采用局部时间步长加速收敛，极大缩短获取稳定电磁场所需计算时间。
6、在求解时变电磁场时通过一些因素包括物理时间迭代精度、物理时间步长、子迭代收敛标准以及最大子迭代步数的选择，通过它们的协同作用保证时间推进效率和计算精度。
7、本发明的
Figure 124093DEST_PATH_IMAGE149
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算。显著区别于显式方法，定常计算的不同的网格单元采用不同局部虚拟时间子迭代步长从而加快该单元电磁场收敛。
8、本发明采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解,用两次循环替代了稀疏矩阵求逆，工程上简单易用。
9、本发明替代传统的Runge-Kutta时间显式推进及空间通量也显式计算的时域有限体积方法，从而放松了网格和显式算法对物理时间步极端限制，这种双时间步隐式时域有限体积方法能高效地求解2维、3维Maxwell方程组，获得时变电磁场的时空分布，在保证格式精度的同时提高时间推进效率，提升计算性能。
附图说明
图1为运用本发明计算电磁场的流程示意图；
图2为TM波物理时间步长对圆柱(ka=2)雷达截面的影响示意图；
图3为TE波物理时间步长对圆柱(ka=2)雷达截面的影响示意图；
图4为时变散射场振荡历程比较(金属球,ka=2)示意图；
图5为双时间步计算金属球双站RCS分布比较(金属球,ka=2)示意图；
图6为电磁场频域收敛历程比较(金属球,ka=2)示意图。
具体实施方式
参照图1至图6对本发明一种时变电磁场的全隐式双时间步计算方法的实施例进一步说明。
参照附图1，整个时域有限体积方法计算电磁场软件按结构可以分为：预处理、电磁场计算和后处理三个部分。预处理主要包括网格数据输入、计算参数数据输入、控制参数输入三个模块，主要用来读入网格数据、计算参数数据输入、控制参数文件，并在此基础上，进行预先处理，为电磁场计算提供计算支撑；电磁场计算包括：空间电磁场MUSCL格式插值、单元分界面通量计算、时间推进、收敛判断模块组成；后处理主要用于输出电磁场的时间、空间分布、目标表面诱导电流密度、雷达散射截面输出。
下面结合要数值模拟的Maxwell方程组两个旋度方程，
Figure 869064DEST_PATH_IMAGE150
，
Figure 460583DEST_PATH_IMAGE151
，其中，
Figure 358132DEST_PATH_IMAGE152
是梯度符号，
Figure 9693DEST_PATH_IMAGE153
是实数型磁感应强度矢量，
Figure 477845DEST_PATH_IMAGE154
是实数型电位移矢量，
Figure 873055DEST_PATH_IMAGE155
是电场强度矢量，
Figure 687427DEST_PATH_IMAGE156
是磁场强度矢量，
Figure 447573DEST_PATH_IMAGE157
是归一化物理时间。
介绍时域有限体积法迭代数值计算过程。时域麦克斯韦方程组的两个旋度方程无源条件下的直角坐标系守恒形式为：
Figure 714606DEST_PATH_IMAGE158
其中，
Figure 100457DEST_PATH_IMAGE003
是电磁场守恒变量，
Figure 769335DEST_PATH_IMAGE004
是归一化物理时间，
Figure 700382DEST_PATH_IMAGE005
是直角坐标系下电磁通量的
Figure 454712DEST_PATH_IMAGE006
分量；
Figure 978726DEST_PATH_IMAGE159
是实数型磁感应强度矢量，
Figure 502111DEST_PATH_IMAGE160
是实数型电位移矢量，
Figure 666376DEST_PATH_IMAGE161
是电场强度矢量，
Figure 845685DEST_PATH_IMAGE162
是磁场强度矢量，
Figure 386388DEST_PATH_IMAGE011
、
Figure 951230DEST_PATH_IMAGE012
Figure 551976DEST_PATH_IMAGE013
为标量，其分别是
Figure 953001DEST_PATH_IMAGE163
的
Figure 297395DEST_PATH_IMAGE014
分量；
Figure 529793DEST_PATH_IMAGE015
、
Figure 989856DEST_PATH_IMAGE016
、
Figure 940494DEST_PATH_IMAGE017
是标量，其分别是
Figure 760683DEST_PATH_IMAGE164
的
Figure 113167DEST_PATH_IMAGE006
分量；
Figure 242665DEST_PATH_IMAGE018
、
Figure 680600DEST_PATH_IMAGE019
、
Figure 366796DEST_PATH_IMAGE020
为标量，其分别是
Figure 511470DEST_PATH_IMAGE165
的
Figure 624919DEST_PATH_IMAGE014
分量；
Figure 235636DEST_PATH_IMAGE021
、
Figure 459944DEST_PATH_IMAGE022
、
Figure 724703DEST_PATH_IMAGE023
是标量，其分别是
Figure 743475DEST_PATH_IMAGE166
的
Figure 156001DEST_PATH_IMAGE006
分量。
对于复杂外形物体，采用的是计算空间贴体多块结构网格，因此均存在坐标变换：
Figure 370951DEST_PATH_IMAGE167
其中
Figure 552534DEST_PATH_IMAGE168
代表曲线坐标系
Figure 945469DEST_PATH_IMAGE169
三个方向，其分别取
Figure 579712DEST_PATH_IMAGE170
之一。得到所要数值模拟的曲线坐标系下的麦克斯韦方程组守恒形：
Figure 834239DEST_PATH_IMAGE171
Figure 135907DEST_PATH_IMAGE172
Figure 762060DEST_PATH_IMAGE173
Figure 821283DEST_PATH_IMAGE174
式中，V是坐标变换的雅可比矩阵， ^上标变量代表曲线坐标系下的对应值，由坐标变换获取。
Figure 191085DEST_PATH_IMAGE175
为曲线坐标系守恒变量，
Figure 534210DEST_PATH_IMAGE176
即为当
Figure 331265DEST_PATH_IMAGE177
分别取曲线坐标系方向
Figure 877784DEST_PATH_IMAGE178
之一时的曲线坐标系下电磁通量。
为了摆脱传统显式时域有限体积全局最小物理时间步的限制带来的大计算量弊端，为此发明提高时变电磁场时间推进效率的全隐式双时间步计算方法，用双时间步保证时间精度的同时根据物理特征选取物理时间步，结合隐式空间通量残差获得稳定、高效的计算流程，包括以下步骤：
步骤1：根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模。
步骤2：采用四边形（二维）或六面体结构（三维）对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，所述网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件。网格密度保证每波长13-20个网格点，壁面密度>300点/每波长，几何奇异处加密到50-100个网格点/每波长，二维网格在垂直该平面按右手法则推进1层，作为3维问题特例统一计算。网格数据文件包括结构网格块数目及每块3个曲线坐标系下维度。
步骤3：预处理部分，输入目标计算电磁参数、数值计算控制参数、有等离子体外部流场情况输入对应流场参数文件。虚拟时间子迭代因是隐式CFL数不受上述显式稳定性要求约束，预先设定物理时间步长，二维问题用入射电磁波周期无量纲化的物理时间步长
Figure 785697DEST_PATH_IMAGE179
计算精度已可以，三维问题选择
Figure 796378DEST_PATH_IMAGE180
量级。同时设定子迭代收敛标准判据值(例如全网格空间相邻子迭代时间步最大幅度差值<0.001,判定为收敛)、设定最大子迭代步数（例：isubmax=50）。
步骤4：输入网格数据和边界条件信息文件，初始化计算空间电磁场。
步骤5：以基于时间迭代推进和空间通量残差的隐式双时间步方式对麦克斯韦方程组时变电磁场进行迭代求解。
步骤5-1：外层物理时间步循环，到计算收敛结束。
步骤5-2：内层虚拟时间步子迭代循环，到子迭代收敛或最大子迭代步数条件满足之一。双时间方法非定常计算的时间精度还受到每一物理时间步上子迭代步数制约。迭代步数多，收敛控制严格，时间精度才有保证。真实时间步小，则子迭代可以取得少一些，数值计算中为防止场变化剧烈处、以及子迭代 CFL 数取得过大等情况时，残差很有可能无法下降到收，陷入死循环而给出最大子迭代步数，另一方面对子迭代收敛进行判断和控制, 在保证一定时间精度前提下尽快结束迭代，子迭代收敛判据使用电场和磁场在计算空间最大幅度差绝对值。
以下介绍步骤5-1和步骤5-2具体算法：采用MUSCL (Monotonic UpstreamSchemes for Conservation Laws,MUSCL)插值,最大CFL数上限为1.745。
通过增加一项虚拟时间
Figure 449820DEST_PATH_IMAGE181
导数项，双时间步法Maxwell控制方程组修改为：
Figure 545952DEST_PATH_IMAGE182
明显可见，当
Figure 195239DEST_PATH_IMAGE183
收敛时，该方程组等同于原始方程组，定常虚拟时间
Figure 60427DEST_PATH_IMAGE184
子迭代表示为：
Figure 386235DEST_PATH_IMAGE185
其中，
Figure 969663DEST_PATH_IMAGE027
是
Figure 484958DEST_PATH_IMAGE003
经虚拟时间子迭代后的电磁场守恒变量值，
Figure 142336DEST_PATH_IMAGE027
是
Figure 452094DEST_PATH_IMAGE028
的近似；
Figure 211234DEST_PATH_IMAGE029
是通量残差，
Figure 264641DEST_PATH_IMAGE030
是通量残差
Figure 42104DEST_PATH_IMAGE029
增加了物理时间导数项后的残差；
Figure 522764DEST_PATH_IMAGE031
是物理时间步数；
Figure 80784DEST_PATH_IMAGE032
是物理时间步长；
Figure 859253DEST_PATH_IMAGE186
为
Figure 553539DEST_PATH_IMAGE027
对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度。
定常子迭代部分采用隐式算法：
Figure 408363DEST_PATH_IMAGE187
其中，
Figure 188100DEST_PATH_IMAGE038
为结构网格曲线坐标系方向1，
Figure 268795DEST_PATH_IMAGE039
为结构网格曲线坐标系方向2，
Figure 83167DEST_PATH_IMAGE040
为结构网格曲线坐标系方向3；
Figure 905630DEST_PATH_IMAGE041
分别对应为
Figure 110346DEST_PATH_IMAGE038
、
Figure 309247DEST_PATH_IMAGE039
、
Figure 165076DEST_PATH_IMAGE040
方向的电磁通量；
Figure 158440DEST_PATH_IMAGE042
是隐式控制参数，取
Figure 850452DEST_PATH_IMAGE043
的全隐式，其他参数对应显、隐混合格式；下标
Figure 587464DEST_PATH_IMAGE044
是网格单元编号，
Figure 110849DEST_PATH_IMAGE045
是第
Figure 229109DEST_PATH_IMAGE044
网格单元第
Figure 470735DEST_PATH_IMAGE046
虚拟时间迭代步时的电磁守恒变量，
Figure 683541DEST_PATH_IMAGE047
是第
Figure 327012DEST_PATH_IMAGE044
网格单元第
Figure 849129DEST_PATH_IMAGE048
虚拟时间迭代步时的电磁守恒变量，
Figure 578051DEST_PATH_IMAGE049
是第
Figure 594548DEST_PATH_IMAGE050
物理时间步时的电磁守恒变量，
Figure 92526DEST_PATH_IMAGE033
是第
Figure 549659DEST_PATH_IMAGE031
物理时间步时的电磁守恒变量，
Figure 500297DEST_PATH_IMAGE034
是第
Figure 382802DEST_PATH_IMAGE035
物理时间步时的电磁守恒变量，
Figure 672970DEST_PATH_IMAGE051
是第
Figure 615518DEST_PATH_IMAGE044
网格单元第
Figure 240403DEST_PATH_IMAGE052
虚拟时间迭代步时的空间通量残差，
Figure 661020DEST_PATH_IMAGE053
是第
Figure 71273DEST_PATH_IMAGE044
网格单元第
Figure 184722DEST_PATH_IMAGE031
虚拟时间迭代步时的空间通量残差；
Figure 798369DEST_PATH_IMAGE054
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算；采用不同局部虚拟时间子迭代步长定常计算不同的网格单元，从而加快相应网格单元的单元电磁场收敛。
步骤5-3：在每个虚拟时间子迭代过程中，逐网格块、逐网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
有限体积法的空间精度体现在能否精确模拟守恒变量Q在网格单元分界面处的状态变量，以得到相应精确的分界面通量
Figure 22676DEST_PATH_IMAGE188
，采用Steger-Warming分裂计算网格单元分界面通量。
Figure 84173DEST_PATH_IMAGE189
；
Figure 306207DEST_PATH_IMAGE057
；
Figure 718734DEST_PATH_IMAGE058
；
式中
Figure 933684DEST_PATH_IMAGE059
分别取曲线坐标系
Figure 115266DEST_PATH_IMAGE060
方向之一，相应的
Figure 242622DEST_PATH_IMAGE061
即为对应
Figure 142445DEST_PATH_IMAGE062
方向的电磁通量；
Figure 394042DEST_PATH_IMAGE063
代表曲线坐标系
Figure 695710DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 993967DEST_PATH_IMAGE064
代表曲线坐标系
Figure 381086DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 750888DEST_PATH_IMAGE065
为相似矩阵，
Figure 828434DEST_PATH_IMAGE066
分别为正负特征值构成的对角矩阵，
Figure 625489DEST_PATH_IMAGE067
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 172008DEST_PATH_IMAGE068
代表自变量为
Figure 345500DEST_PATH_IMAGE069
的
Figure 44597DEST_PATH_IMAGE070
相似矩阵；
Figure 12553DEST_PATH_IMAGE071
代表自变量为
Figure 46368DEST_PATH_IMAGE069
的
Figure 757972DEST_PATH_IMAGE072
对角矩阵；
Figure 623160DEST_PATH_IMAGE073
代表自变量为
Figure 948968DEST_PATH_IMAGE069
的
Figure 532396DEST_PATH_IMAGE074
相似矩阵；
Figure 719795DEST_PATH_IMAGE075
代表自变量为
Figure 705068DEST_PATH_IMAGE076
的
Figure 14827DEST_PATH_IMAGE070
相似矩阵；
Figure 771037DEST_PATH_IMAGE077
代表自变量为
Figure 824444DEST_PATH_IMAGE076
的
Figure 601907DEST_PATH_IMAGE078
对角矩阵；
Figure 82567DEST_PATH_IMAGE079
代表自变量为
Figure 561958DEST_PATH_IMAGE076
的
Figure 419056DEST_PATH_IMAGE074
相似矩阵。
Figure 51026DEST_PATH_IMAGE190
Figure 702587DEST_PATH_IMAGE191
其中，
Figure 747903DEST_PATH_IMAGE082
是限制器，下标
Figure 831528DEST_PATH_IMAGE083
是网格单元编号，
Figure 380321DEST_PATH_IMAGE084
对应单元分界面，
Figure 406046DEST_PATH_IMAGE085
是3阶精度格式的控制参数，
Figure 673079DEST_PATH_IMAGE086
和
Figure 871979DEST_PATH_IMAGE087
分别是后差和前差算符；
Figure 462229DEST_PATH_IMAGE088
表示网格单元
Figure 721172DEST_PATH_IMAGE089
分界面处左状态电磁守恒变量，
Figure 413185DEST_PATH_IMAGE090
表示网格单元
Figure 150197DEST_PATH_IMAGE089
分界面处右状态电磁守恒变量；
Figure 359068DEST_PATH_IMAGE091
是第
Figure 788912DEST_PATH_IMAGE083
个网格单元电磁场守恒变量，
Figure 764958DEST_PATH_IMAGE092
是第
Figure 243344DEST_PATH_IMAGE083
+1个网格单元电磁场守恒变量。
空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解,用两次循环替代了稀疏矩阵求逆，工程上简单易用。经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到计
Figure 886815DEST_PATH_IMAGE192
其中，
Figure 408932DEST_PATH_IMAGE094
是分裂后的系数矩阵，
Figure 137854DEST_PATH_IMAGE095
是上一个物理时间步计算的空间通量残差，
Figure 154351DEST_PATH_IMAGE193
是隐式子迭代电磁场差值；将该方程表示为LDU，近似因式分解得：
Figure 652329DEST_PATH_IMAGE097
Figure 158397DEST_PATH_IMAGE098
Figure 63030DEST_PATH_IMAGE099
Figure 945535DEST_PATH_IMAGE100
其中，下标
Figure 235702DEST_PATH_IMAGE101
是网格单元编号，
Figure 178250DEST_PATH_IMAGE102
是雅可比系数矩阵最大特征值分裂参数，
Figure 803136DEST_PATH_IMAGE103
是雅克比系数矩阵最大特征值；
Figure 223753DEST_PATH_IMAGE104
是单位对角矩阵，
Figure 696322DEST_PATH_IMAGE105
是对角矩阵，
Figure 481876DEST_PATH_IMAGE106
为上三角矩阵，
Figure 407106DEST_PATH_IMAGE107
为下三角矩阵，
Figure 340338DEST_PATH_IMAGE108
是对应上三角矩阵的电磁守恒变量差值，
Figure 667414DEST_PATH_IMAGE109
是对应下三角矩阵的电磁守恒变量差值；
Figure 889448DEST_PATH_IMAGE110
指
Figure 301975DEST_PATH_IMAGE111
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 64394DEST_PATH_IMAGE112
指
Figure 432928DEST_PATH_IMAGE113
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 888180DEST_PATH_IMAGE114
指
Figure 725686DEST_PATH_IMAGE115
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 291796DEST_PATH_IMAGE116
指
Figure 16301DEST_PATH_IMAGE117
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 642454DEST_PATH_IMAGE118
指
Figure 29573DEST_PATH_IMAGE119
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 71479DEST_PATH_IMAGE120
指
Figure 227653DEST_PATH_IMAGE121
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 211659DEST_PATH_IMAGE122
是指相邻
Figure 820495DEST_PATH_IMAGE123
网格单元分裂后的系数矩阵；
Figure 993987DEST_PATH_IMAGE124
是指相邻
Figure 942352DEST_PATH_IMAGE125
网格单元分裂后的系数矩阵；
Figure 910308DEST_PATH_IMAGE126
是指相邻
Figure 691925DEST_PATH_IMAGE127
网格单元分裂后的系数矩阵；
Figure 137950DEST_PATH_IMAGE128
是指相邻
Figure 206400DEST_PATH_IMAGE129
网格单元分裂后的系数矩阵；
Figure 345258DEST_PATH_IMAGE130
是指相邻
Figure 928686DEST_PATH_IMAGE131
网格单元分裂后的系数矩阵；
Figure 365352DEST_PATH_IMAGE132
是指相邻
Figure 350626DEST_PATH_IMAGE133
网格单元分裂后的系数矩阵。
最后得到前后向迭代计算电磁场差值
Figure 598067DEST_PATH_IMAGE194
：
Figure 668792DEST_PATH_IMAGE195
Figure 410614DEST_PATH_IMAGE196
其中，
Figure 250394DEST_PATH_IMAGE138
是对角矩阵
Figure 668737DEST_PATH_IMAGE139
的逆矩阵，
Figure 961178DEST_PATH_IMAGE140
、
Figure 5226DEST_PATH_IMAGE141
分别是根据
Figure 699513DEST_PATH_IMAGE142
，
Figure 351074DEST_PATH_IMAGE143
计算得出的上三角矩阵、下三角矩阵。
前向循环：
Figure 334073DEST_PATH_IMAGE197
后向循环
Figure 729283DEST_PATH_IMAGE198
，其中，
Figure 963561DEST_PATH_IMAGE146
是电磁守恒变量差值的中间过渡变量。
以上是隐式双时间步计算Maxwell方程组控制电磁场迭代过程。
步骤6：收敛判断，后处理过程，输出电磁场的时间、空间分布，输出表面诱导电流和雷达散射截面空间分布数据等。
如图2和图3所示，比较了TM波和TE波不同物理时间步长对圆柱RCS结果的影响，k为波数，a为圆柱半径，计算网格为51x46，其中物面网格选取每波长20个网格点，远场边界在3波长外，辐向网格壁面加密，平均约每波长15个网格点，随着物理时间间隔的增大, 计算结果与参考值差别逐渐增大, 在过大的物理时间步长条件下, 双时间步方法的计算精度变低。
如图4-图6所示，以金属球电磁散射为例，电尺寸为ka=2,a为球半径，计算网格为61x31x46,同一空间点分别用显式Runge-Kutta法和双时间步方法计算，其散射电磁场量的时间振荡历程，空间两者几无区别，说明物理时间后向差分、2阶精度的双时间方法精度满足傅里叶变换需要。图5比较了显式和双时间步计算的雷达截面双站分布比较，可见后者在前向（
Figure 51603DEST_PATH_IMAGE199
）与解析解吻合更好。图6则是这两种方法逐周期进行傅里叶变换所得频域值收敛历程比较，在达到同样的收敛标准（0.001），双时间步收敛略快，这对大网格量和存在很小体积网格单元的情况，能大大提高计算效率，因为此时显式
Figure 256320DEST_PATH_IMAGE200
量级甚至更小，相反本发明双时间步方法则根据时间精度选择
Figure 189641DEST_PATH_IMAGE201
。
这几个数值算例表明，该本发明对应的隐式双时间步时域有限体积计算方法能在放松时间步限制的同时，保证数值精度并提升计算效率。
以上仅为本发明的较佳实施例而已，并不用以限制本发明，应当指出的是，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种时变电磁场的全隐式双时间步计算方法，其特征在于，其包括以下步骤：
步骤1，根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；
步骤2，采用四边形或六面体结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件；
步骤3，输入目标计算电磁参数、数值计算控制参数；
步骤4，输入网格数据和边界条件信息文件，初始化计算空间电磁场；
步骤5，以基于时间迭代推进和空间通量残差的隐式双时间步方式对麦克斯韦方程组时变电磁场进行迭代求解：仿真模型的外层为物理时间步循环，直至计算收敛结束；仿真模型的内层为虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格、各个网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值；
所述物理时间步循环和虚拟时间步子迭代循环过程为：
通过虚拟时间导数项
Figure 951680DEST_PATH_IMAGE001
，将待求解的麦克斯韦方程组修正为：
Figure 517790DEST_PATH_IMAGE002
其中，
Figure 757142DEST_PATH_IMAGE003
是实数型的电磁场守恒变量，
Figure 304667DEST_PATH_IMAGE004
是归一化物理时间，
Figure 691786DEST_PATH_IMAGE005
是直角坐标系下电磁通量的
Figure 186221DEST_PATH_IMAGE006
分量；
Figure 14500DEST_PATH_IMAGE007
是实数型磁感应强度矢量，
Figure 811554DEST_PATH_IMAGE008
是实数型电位移矢量，
Figure 607341DEST_PATH_IMAGE009
是电场强度矢量，
Figure 780833DEST_PATH_IMAGE010
是磁场强度矢量，
Figure 729198DEST_PATH_IMAGE011
、
Figure 884104DEST_PATH_IMAGE012
、
Figure 917919DEST_PATH_IMAGE013
分别是
Figure 629524DEST_PATH_IMAGE007
的
Figure 494711DEST_PATH_IMAGE014
分量；
Figure 820519DEST_PATH_IMAGE015
、
Figure 341631DEST_PATH_IMAGE016
、
Figure 591346DEST_PATH_IMAGE017
分别是
Figure 763570DEST_PATH_IMAGE008
的
Figure 73329DEST_PATH_IMAGE006
分量；
Figure 81736DEST_PATH_IMAGE018
、
Figure 322094DEST_PATH_IMAGE019
、
Figure 161874DEST_PATH_IMAGE020
分别是
Figure 642534DEST_PATH_IMAGE009
的
Figure 872658DEST_PATH_IMAGE014
分量；
Figure 729755DEST_PATH_IMAGE021
、
Figure 610993DEST_PATH_IMAGE022
、
Figure 262554DEST_PATH_IMAGE023
分别是
Figure 245553DEST_PATH_IMAGE010
的
Figure 827713DEST_PATH_IMAGE006
分量；当
Figure 376506DEST_PATH_IMAGE024
收敛时，该方程组等同于原始方程组，定常虚拟时间
Figure 464548DEST_PATH_IMAGE025
子迭代表示为：
Figure 669264DEST_PATH_IMAGE026
其中，
Figure 868165DEST_PATH_IMAGE027
是
Figure 458415DEST_PATH_IMAGE003
经虚拟时间子迭代后的电磁场守恒变量值，
Figure 717358DEST_PATH_IMAGE027
是
Figure 409370DEST_PATH_IMAGE028
的近似；
Figure 146382DEST_PATH_IMAGE029
是通量残差，
Figure 856718DEST_PATH_IMAGE030
是通量残差
Figure 286562DEST_PATH_IMAGE029
增加了物理时间导数项后的残差；
Figure 200292DEST_PATH_IMAGE031
是物理时间步数；
Figure 740994DEST_PATH_IMAGE032
是物理时间步长，
Figure 571416DEST_PATH_IMAGE033
是第
Figure 844266DEST_PATH_IMAGE031
物理时间步时的电磁守恒变量，
Figure 573187DEST_PATH_IMAGE034
是第
Figure 838952DEST_PATH_IMAGE035
物理时间步时的电磁守恒变量；
Figure 336930DEST_PATH_IMAGE036
为
Figure 780681DEST_PATH_IMAGE027
对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度；定常子迭代部分采用隐式算法：
Figure DEST_PATH_IMAGE037
其中，
Figure 183849DEST_PATH_IMAGE038
为结构网格曲线坐标系方向1，
Figure 4037DEST_PATH_IMAGE039
为结构网格曲线坐标系方向2，
Figure 356521DEST_PATH_IMAGE040
为结构网格曲线坐标系方向3；
Figure 486020DEST_PATH_IMAGE041
分别对应为
Figure 861638DEST_PATH_IMAGE038
、
Figure 282255DEST_PATH_IMAGE039
、
Figure 941775DEST_PATH_IMAGE040
方向的电磁通量；
Figure 789646DEST_PATH_IMAGE042
是隐式控制参数，取
Figure 652559DEST_PATH_IMAGE043
的全隐式，其他参数对应显、隐混合格式；下标
Figure 87255DEST_PATH_IMAGE044
是网格单元编号，
Figure 414332DEST_PATH_IMAGE045
是第
Figure 636366DEST_PATH_IMAGE044
网格单元第
Figure 48892DEST_PATH_IMAGE046
虚拟时间迭代步时的电磁守恒变量，
Figure 998263DEST_PATH_IMAGE047
是第
Figure 117528DEST_PATH_IMAGE044
网格单元第
Figure 572780DEST_PATH_IMAGE048
虚拟时间迭代步时的电磁守恒变量，
Figure 659554DEST_PATH_IMAGE049
是第
Figure 225665DEST_PATH_IMAGE050
物理时间步时的电磁守恒变量，
Figure 199437DEST_PATH_IMAGE033
是第
Figure 825590DEST_PATH_IMAGE031
物理时间步时的电磁守恒变量，
Figure 399660DEST_PATH_IMAGE034
是第
Figure 503882DEST_PATH_IMAGE035
物理时间步时的电磁守恒变量，
Figure 597740DEST_PATH_IMAGE051
是第
Figure 394795DEST_PATH_IMAGE044
网格单元第
Figure 190581DEST_PATH_IMAGE052
虚拟时间迭代步时的空间通量残差，
Figure 301757DEST_PATH_IMAGE053
是第
Figure 312438DEST_PATH_IMAGE044
网格单元第
Figure 467345DEST_PATH_IMAGE031
虚拟时间迭代步时的空间通量残差；
Figure 563477DEST_PATH_IMAGE054
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算；采用不同局部虚拟时间子迭代步长定常计算不同的网格单元，加快相应网格单元的单元电磁场收敛。
2.根据权利要求1所述的一种时变电磁场的全隐式双时间步计算方法，其特征在于：所述空间通量计算和隐式迭代解计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量
Figure 947185DEST_PATH_IMAGE055
Figure 77952DEST_PATH_IMAGE056
；
Figure 403760DEST_PATH_IMAGE057
；
Figure 987188DEST_PATH_IMAGE058
；
式中
Figure 174587DEST_PATH_IMAGE059
分别取曲线坐标系
Figure 159860DEST_PATH_IMAGE060
方向之一，相应的
Figure 656570DEST_PATH_IMAGE061
即为对应
Figure 727294DEST_PATH_IMAGE062
方向的电磁通量；
Figure 718384DEST_PATH_IMAGE063
代表曲线坐标系
Figure 745114DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 225774DEST_PATH_IMAGE064
代表曲线坐标系
Figure 518215DEST_PATH_IMAGE060
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 312996DEST_PATH_IMAGE065
，
Figure 194233DEST_PATH_IMAGE066
为相似矩阵，
Figure 845794DEST_PATH_IMAGE067
分别为正负特征值构成的对角矩阵，
Figure 828794DEST_PATH_IMAGE068
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 224003DEST_PATH_IMAGE069
代表自变量为
Figure 772796DEST_PATH_IMAGE070
的
Figure 47789DEST_PATH_IMAGE065
相似矩阵；
Figure 252505DEST_PATH_IMAGE071
代表自变量为
Figure 185826DEST_PATH_IMAGE070
的
Figure 854705DEST_PATH_IMAGE072
对角矩阵；
Figure 300598DEST_PATH_IMAGE073
代表自变量为
Figure 789348DEST_PATH_IMAGE070
的
Figure 729623DEST_PATH_IMAGE066
相似矩阵；
Figure 253008DEST_PATH_IMAGE074
代表自变量为
Figure 604224DEST_PATH_IMAGE075
的
Figure 845849DEST_PATH_IMAGE065
相似矩阵；
Figure 324235DEST_PATH_IMAGE076
代表自变量为
Figure 702127DEST_PATH_IMAGE075
的
Figure 489823DEST_PATH_IMAGE077
对角矩阵；
Figure 953165DEST_PATH_IMAGE078
代表自变量为
Figure 235242DEST_PATH_IMAGE075
的
Figure 733220DEST_PATH_IMAGE066
相似矩阵；
Figure 426238DEST_PATH_IMAGE079
；
Figure 642456DEST_PATH_IMAGE080
；
其中
Figure 462644DEST_PATH_IMAGE081
是限制器，下标
Figure 815128DEST_PATH_IMAGE082
是网格单元编号，
Figure 679048DEST_PATH_IMAGE083
对应单元分界面，
Figure 116982DEST_PATH_IMAGE084
是3阶精度格式的控制参数，
Figure 740862DEST_PATH_IMAGE085
和
Figure 947852DEST_PATH_IMAGE086
分别是后差和前差算符；
Figure 248252DEST_PATH_IMAGE087
表示网格单元
Figure 173483DEST_PATH_IMAGE088
分界面处左状态电磁守恒变量，
Figure 335474DEST_PATH_IMAGE089
表示网格单元
Figure 662550DEST_PATH_IMAGE088
分界面处右状态电磁守恒变量；
Figure 133852DEST_PATH_IMAGE090
是第
Figure 280799DEST_PATH_IMAGE082
个网格单元电磁场守恒变量，
Figure 246481DEST_PATH_IMAGE091
是第
Figure 428064DEST_PATH_IMAGE082
+1个网格单元电磁场守恒变量；
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure 70267DEST_PATH_IMAGE092
其中，
Figure 970089DEST_PATH_IMAGE093
是分裂后的系数矩阵，
Figure 473883DEST_PATH_IMAGE094
是上一个物理时间步计算的空间通量残差，
Figure 509972DEST_PATH_IMAGE095
是隐式子迭代电磁场差值；
将
Figure 323076DEST_PATH_IMAGE092
表示为LDU，近似因式分解得：
Figure 444616DEST_PATH_IMAGE096
Figure 752101DEST_PATH_IMAGE097
Figure 908276DEST_PATH_IMAGE098
Figure 892281DEST_PATH_IMAGE099
其中，下标
Figure 501117DEST_PATH_IMAGE100
是网格单元编号，
Figure 346713DEST_PATH_IMAGE101
是雅可比系数矩阵最大特征值分裂参数，
Figure 357394DEST_PATH_IMAGE102
是雅克比系数矩阵最大特征值；
Figure 512301DEST_PATH_IMAGE103
是单位对角矩阵，
Figure 608433DEST_PATH_IMAGE104
是对角矩阵，
Figure 257720DEST_PATH_IMAGE105
为上三角矩阵，
Figure 388487DEST_PATH_IMAGE106
为下三角矩阵，
Figure 714295DEST_PATH_IMAGE107
是对应上三角矩阵的电磁守恒变量差值，
Figure 32144DEST_PATH_IMAGE108
是对应下三角矩阵的电磁守恒变量差值；
Figure 485122DEST_PATH_IMAGE109
指
Figure 470396DEST_PATH_IMAGE110
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 967105DEST_PATH_IMAGE111
指
Figure 772250DEST_PATH_IMAGE112
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 763340DEST_PATH_IMAGE113
指
Figure 603120DEST_PATH_IMAGE114
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 270730DEST_PATH_IMAGE115
指
Figure 766434DEST_PATH_IMAGE116
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 357952DEST_PATH_IMAGE117
指
Figure 317818DEST_PATH_IMAGE118
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 156330DEST_PATH_IMAGE119
指
Figure 873750DEST_PATH_IMAGE120
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 268959DEST_PATH_IMAGE121
是指相邻
Figure 270282DEST_PATH_IMAGE122
网格单元分裂后的系数矩阵；
Figure 92745DEST_PATH_IMAGE123
是指相邻
Figure 297461DEST_PATH_IMAGE124
网格单元分裂后的系数矩阵；
Figure 496361DEST_PATH_IMAGE125
是指相邻
Figure 352191DEST_PATH_IMAGE126
网格单元分裂后的系数矩阵；
Figure 345554DEST_PATH_IMAGE127
是指相邻
Figure 37567DEST_PATH_IMAGE128
网格单元分裂后的系数矩阵；
Figure 40158DEST_PATH_IMAGE129
是指相邻
Figure 750494DEST_PATH_IMAGE130
网格单元分裂后的系数矩阵；
Figure 914759DEST_PATH_IMAGE131
是指相邻
Figure 94068DEST_PATH_IMAGE132
网格单元分裂后的系数矩阵；
最后，得到前后向迭代计算电磁场差值
Figure 634770DEST_PATH_IMAGE133
：
Figure 199613DEST_PATH_IMAGE134
Figure 800358DEST_PATH_IMAGE135
Figure 201384DEST_PATH_IMAGE136
其中，
Figure 545778DEST_PATH_IMAGE137
是对角矩阵
Figure 965126DEST_PATH_IMAGE138
的逆矩阵，
Figure 736773DEST_PATH_IMAGE139
、
Figure 625095DEST_PATH_IMAGE140
分别是根据
Figure 507600DEST_PATH_IMAGE141
、
Figure 47035DEST_PATH_IMAGE142
计算得出的上三角矩阵、下三角矩阵；
前向循环：
Figure 989583DEST_PATH_IMAGE143
后向循环：
Figure 365201DEST_PATH_IMAGE144
其中，
Figure 51397DEST_PATH_IMAGE145
是电磁守恒变量差值的中间过渡变量。
3.根据权利要求1所述的一种时变电磁场的全隐式双时间步计算方法，其特征在于：所述步骤2中：网格密度保证每波长13-20个网格点，壁面密度>300点/波长，几何奇异处加密到50-100个网格点/波长；二维网格时，则在垂直该二维网格所在平面按右手法则推进1层，作为三维问题特例统一计算；网格数据文件包括结构网格块数目及每块3个曲线坐标系下维度。
4.根据权利要求1所述的一种时变电磁场的全隐式双时间步计算方法，其特征在于：所述步骤3中：如有等离子体外部流场情况，则还应输入对应流场参数。
5.根据权利要求1所述的一种时变电磁场的全隐式双时间步计算方法，其特征在于：所述步骤3中，预先设定物理时间步长，二维问题设定物理时间步长
Figure 468776DEST_PATH_IMAGE146
，三维问题设定
Figure 582225DEST_PATH_IMAGE147
；同时设定子迭代收敛标准判据值为全网格空间相邻子迭代时间步最大幅度差值<0.001，判定为收敛；设定最大子迭代步数为50。
6.根据权利要求1所述的一种时变电磁场的全隐式双时间步计算方法，其特征在于：所述步骤5中，子迭代收敛判据为电场和磁场在计算空间最大幅度差的绝对值，且数值计算中设置有最大子迭代步数。
CN202110527205.5A 2021-05-14 2021-05-14 一种时变电磁场的全隐式双时间步计算方法 Active CN113158492B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110527205.5A CN113158492B (zh) 2021-05-14 2021-05-14 一种时变电磁场的全隐式双时间步计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110527205.5A CN113158492B (zh) 2021-05-14 2021-05-14 一种时变电磁场的全隐式双时间步计算方法

Publications (2)

Publication Number Publication Date
CN113158492A CN113158492A (zh) 2021-07-23
CN113158492B true CN113158492B (zh) 2021-08-20

Family

ID=76875078

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110527205.5A Active CN113158492B (zh) 2021-05-14 2021-05-14 一种时变电磁场的全隐式双时间步计算方法

Country Status (1)

Country Link
CN (1) CN113158492B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114357904A (zh) * 2021-12-29 2022-04-15 中国航天空气动力技术研究院 一种基于金属流体流场的电磁主动控制方法及设备
CN114492251B (zh) * 2022-04-18 2022-07-15 国家超级计算天津中心 超算环境的低速流场发散处理方法、装置、设备及介质
CN116401935A (zh) * 2023-02-21 2023-07-07 哈尔滨工业大学 一种建筑动态热负荷神经网络预测方法及系统
CN116738891B (zh) * 2023-08-08 2023-10-13 中国空气动力研究与发展中心计算空气动力研究所 一种增强飞行器流场模拟稳定性的lu-sgs改进方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8608826B2 (en) * 2011-04-11 2013-12-17 King Fahd University Of Petroleum And Minerals Method of modeling fly ash collection efficiency in wire-duct electrostatic precipitators
CN103970717A (zh) * 2014-05-08 2014-08-06 中国人民解放军理工大学 基于Associated Hermite 正交函数的无条件稳定FDTD算法
CN109948293A (zh) * 2019-04-02 2019-06-28 安徽大学 一种随机混合显隐式时域有限差分方法
CN111639447A (zh) * 2020-04-30 2020-09-08 南京理工大学 多级局部时间步进技术的任意高阶混合网格时域不连续伽辽金方法
CN112307644A (zh) * 2020-11-20 2021-02-02 金陵科技学院 一种电大尺寸目标的rcs计算方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107766288A (zh) * 2017-10-19 2018-03-06 中国空气动力研究与发展中心计算空气动力研究所 针对高精度格式的鲁棒高效隐式时间推进方法
CN112347667A (zh) * 2020-09-28 2021-02-09 中国民用航空中南地区空中交通管理局 一种用于仪表着陆系统的电磁仿真方法及电子设备

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8608826B2 (en) * 2011-04-11 2013-12-17 King Fahd University Of Petroleum And Minerals Method of modeling fly ash collection efficiency in wire-duct electrostatic precipitators
CN103970717A (zh) * 2014-05-08 2014-08-06 中国人民解放军理工大学 基于Associated Hermite 正交函数的无条件稳定FDTD算法
CN109948293A (zh) * 2019-04-02 2019-06-28 安徽大学 一种随机混合显隐式时域有限差分方法
CN111639447A (zh) * 2020-04-30 2020-09-08 南京理工大学 多级局部时间步进技术的任意高阶混合网格时域不连续伽辽金方法
CN112307644A (zh) * 2020-11-20 2021-02-02 金陵科技学院 一种电大尺寸目标的rcs计算方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"A scalable, fully implicit algorithm for the reduced two-field low-beta extended MHD model";Chacon, L.等;《JOURNAL OF COMPUTATIONAL PHYSICS》;20161201;第326卷;第763-772页 *
"双时间步时域有限体积方法计算时变电磁场";许勇 等;《电波科学学报》;20150815;第30卷(第4期);第647-652页 *
"用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术";孙建国 等;《石油地球物理勘探》;20090815;第44卷(第4期);第494-500页 *

Also Published As

Publication number Publication date
CN113158492A (zh) 2021-07-23

Similar Documents

Publication Publication Date Title
CN113158492B (zh) 一种时变电磁场的全隐式双时间步计算方法
CN112989680B (zh) 减少网格使用量的fvfd远场积分边界条件计算方法
Zhuang et al. Fracture modeling using meshless methods and level sets in 3D: framework and modeling
CN113158527B (zh) 一种基于隐式fvfd计算频域电磁场的方法
Zheng et al. A wideband fast multipole boundary element method for three dimensional acoustic shape sensitivity analysis based on direct differentiation method
Özgün et al. MATLAB-based finite element programming in electromagnetic modeling
CN115600435A (zh) 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置
Liu et al. Efficient Kriging-based aerodynamic design of transonic airfoils: some key issues
CN116738891B (zh) 一种增强飞行器流场模拟稳定性的lu-sgs改进方法
Barbarino et al. A BEM–FMM approach applied to the combined convected Helmholtz integral formulation for the solution of aeroacoustic problems
Lu et al. NNW-GridStar: interactive structured mesh generation software for aircrafts
Wang et al. An adaptive fast multipole boundary face method with higher order elements for acoustic problems in three-dimension
CN111079326B (zh) 二维各向异性网格单元度量张量场光滑化方法
KR20120057274A (ko) Ｉｄ－ｆｄｔｄ 방법을 이용한 전자기파 수치해석 방법
Zheng et al. Three dimensional acoustic shape sensitivity analysis by means of adjoint variable method and fast multipole boundary element approach
CN108090296B (zh) 基于高阶辛紧致格式的波导全波分析方法
CN115879276A (zh) 一种目标对象的电磁特性分析方法、装置及设备和介质
CN116151135B (zh) 一种电大尺寸目标的电磁仿真方法及系统
Ning et al. A grid generator for 3-D explosion simulations using the staircase boundary approach in Cartesian coordinates based on STL models
Nicolas et al. Improved adaptive mesh refinement for conformal hexahedral meshes
CN105589980A (zh) 一种阻抗匹配层的截断边界
Tong et al. Multilevel fast multipole algorithm for acoustic wave scattering by truncated ground with trenches
Huang et al. A novel fdtd cells generation technology and matlab-gui implementation
Hollander et al. Adaptive multilevel nonuniform grid algorithm for the accelerated analysis of composite metallic–dielectric radomes
Dziekonski et al. GPU-accelerated finite element method

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
GR01 Patent grant
GR01 Patent grant