CN115600435A - 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置 - Google Patents

一种介质涂覆导体复合目标电磁散射隐式计算方法及装置 Download PDF

Info

Publication number
CN115600435A
CN115600435A CN202211411927.5A CN202211411927A CN115600435A CN 115600435 A CN115600435 A CN 115600435A CN 202211411927 A CN202211411927 A CN 202211411927A CN 115600435 A CN115600435 A CN 115600435A
Authority
CN
China
Prior art keywords
electromagnetic
grid
medium
iteration
calculation
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
CN202211411927.5A
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.)
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 CN202211411927.5A priority Critical patent/CN115600435A/zh
Publication of CN115600435A publication Critical patent/CN115600435A/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
    • 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

Landscapes

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

Abstract

本申请公开了一种介质涂覆导体复合目标电磁散射隐式计算方法及装置，包括：获取目标数据包括网格数据文件、边界条件文件、目标计算电磁参数、数值计算控制参数，并初始化计算空间电磁场；以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组进行迭代求解：在仿真模型的外层采用物理时间步循环，内层采用虚拟时间步子迭代循环；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量和隐式迭代解计算，更新下一级数值。能够在保障精度的同时提升计算效率，节约成本。

Description

一种介质涂覆导体复合目标电磁散射隐式计算方法及装置
技术领域
本申请涉及电磁技术领域，特别涉及一种介质涂覆导体复合目标电磁散射隐式计算方法及装置。
背景技术
导体/介质复合目标电磁散射问题，有两大类数值模拟方法，一类是求解电磁流积分程组的数值算法—矩量法系列。主要采用表面耦合积分方程(Coupled IntegralEquation，CIE)来分析，在导体表面用电场积分方程(electric field integralequations，EFIE)、磁场积分方程(magnetic field integral equations，MFIE)或组合场积分方程(combined field integral equations，CFIE)来描述，而介质表面则用EFIE或Poggio-Miller-Chang-Harrington-Wu-Tsai(PMCHWT)方程来分析，介质表面的等效电流和等效磁流都需要计算，对于金属介质复合结构目标，面积分方程在导体介质交接面及介质介质交接面需要特别对待以保证较高的精度，面积分方程通常适用于分段均匀介质目标，即使如此随着不均匀介质数目增加，需要考虑不同介质交接面的边界条件，使得面积分方程的效率下降。对于非均匀介质则必须采用解体面混合积分方程算法。综上，尽管矩量法在解决电大尺寸导体目标电磁散射方面取得巨大成功，但对于含介质/导体复合目标，不管采用面电流基函数的面积分方程，还是采用体电流基函数的体面积分方程，未知量较纯金属体散射问题大幅增加，计算效率下降。另一大类是求解电磁学麦克斯韦方程组或赫姆霍兹波动方程，进一步得到电磁场时间、空间分布的微分类算法。在微分方程方法中，传统常用的有有限元法(Finite Element Method，FEM)和有限差分(Finite-Difference Time-Domain，FDTD)方法。FEM和FDTD方法在分析散射体目标时，需要用网格把连续的三维空间区域划分开。对于有限尺寸的空间区域，这两种方法可以有效地对目标分析求解；对于散射体外是无限大空间的情况，需要人为设置边界条件使求解区域具有有限尺寸，这种做法会带来截断误差和网格的色散误差，而且对计算机硬件资源也有较高的要求。另一方面，上世纪90年代以来，与欧拉方程相同的双曲型数学特征促进计算流体力学(Computational FluidDynamics,CFD)技术在电磁场计算中的应用，其中时域有限差分法(Finite DifferenceTime Domain,FDTD)和时域有限体积法(Finite Volume Time Domain,FVTD)最为著名。20世纪60年代K.S.Yee发表先驱性的时域有限差分算法—FDTD，采用笛卡尔直角正交网格，模拟复杂外形容易产生阶梯效应带来误差。微分类电磁场解法不同于积分方程类数值算法的优势在于，介质的引入并不会带来未知量的大量增加，介质主要影响电磁场量之间的本征关系。
介质电磁参数如介电常数、磁化率的变动，电磁波长也相应变化，传统时域有限差分法和时域有限体积法时间推进采用显式格式，时间显式方法有一个最大缺陷是受稳定性限制,整个计算空间必须采用统一的最小全局计算步长,为模拟几何外形剧烈变化生成的贴体加密网格(例如机翼前后缘棱边几何奇异性带来电磁场梯度剧烈变化要求加密网格仔细模拟以及电磁多尺度问题)，会带来很小的全局时间步长，大的网格单元就需要更多时间步在该单元传播信息，从而使得到稳定的时变电磁场需要较长的计算时间,尤其当介质内电磁波长小于真空电磁波长，显式格式稳定性要求使得物理时间步比真空电磁散射问题更小，带来更长计算时间的弊端,结果就是受稳定性限制的小时间步长带来时域电磁场计算量的显著增加，消耗大量计算资源。其次，介质/导体分界面，由于导体内电磁场量处处为0，所以没有布置计算网格，在介质/导体边界，每个时间步都必须特殊处理，给场量插值以及信息交换和并行带来不便，影响计算效率。
综上，如何在保障精度的同时提升介质涂覆导体复合目标电磁散射问题的计算效率，节约成本是目前亟待解决的问题。
发明内容
有鉴于此，本申请的目的在于提供一种介质涂覆导体复合目标电磁散射隐式计算方法及装置，能够在保障精度的同时提升介质涂覆导体复合目标电磁散射问题的计算效率，从而节约成本。其具体方案如下：
第一方面，本申请公开了一种介质涂覆导体复合目标电磁散射隐式计算方法，包括：
获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；
基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
可选的，获取网格数据文件，包括：
采用四边形或六面体结构对仿真模型进行网格剖分，网格在散射壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，得到网格数据文件。
可选的，所述物理时间步循环和虚拟时间步子迭代循环过程为：
增加虚拟时间导数项
Figure BDA0003939007040000031
将待求解的介质时域麦克斯韦方程组修正为：
Figure BDA0003939007040000032
Figure BDA0003939007040000033
Figure BDA0003939007040000041
Figure BDA0003939007040000042
其中，W是电磁场守恒变量，μr是相对磁化率，μ0是真空磁化率，εr是相对介电常数，ε0是真空介电常数，σm是磁导率，σe是电导率。Q是散射电磁场原始变量，τ是虚拟时间，t是归一化物理时间，Fx、Fy、Fz是直角坐标系下电磁通量的X、y、z分量，S是介质磁导率、电导率引发的外加源项，J是外加强迫电流矢量，Jx、Jy、Jz是直角坐标系下外加强迫电流的X、y、z分量，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量，B是实数型散射磁感应强度矢量，D是实数型散射场电位移矢量，Ei是入射场电场强度矢量，Hi是入射场磁场强度矢量，E是散射场电场强度矢量，H是散射场磁场强度矢量，含下标X、y、z标量分别是对应矢量的直角坐标系X、y、z分量，上标t对应的是电磁场总场形式；
当
Figure BDA0003939007040000043
收敛时，该方程组等同于原始方程组，定常虚拟时间τ子迭代表示为：
Figure BDA0003939007040000044
Figure BDA0003939007040000045
其中，W*是W经虚拟时间子迭代后的电磁场守恒变量值，W*是Wn+1的近似，*号上标代表的是虚拟时间子迭代对应守恒变量值，R是通量残差，R*是通量残差R增加了物理时间导数项以及源项后的残差；n是物理时间步数；Δt是物理时间步长，Wn是第n物理时间步时的电磁守恒变量，Wn-1是第n-1物理时间步时的电磁守恒变量；R*(W*)为W*对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度；定常子迭代部分采用隐式算法：
Figure BDA0003939007040000051
其中,ξ、η、ζ分别表示结构网格曲线坐标系的三个不同方向；F,G,H是分别对应曲线坐标系ξ，η，ζ方向电磁通量；
Figure BDA0003939007040000052
是对应曲线坐标系源项；ω是隐式控制参数，取ω＝1的全隐式；下标i，j，k是网格单元编号，
Figure BDA0003939007040000053
是第i，j，k网格单元第n+1虚拟时间步的电磁守恒变量，
Figure BDA0003939007040000054
是第i，j，k网格单元第n虚拟时间步的电磁守恒变量，Wn+1是第n+1物理时间步的电磁守恒变量，Wn是第n物理时间步的电磁守恒变量，Wn-1是第n-1物理时间步的电磁守恒变量，
Figure BDA0003939007040000055
是第i，j，k网格单元第n+1虚拟时间步的空间通量残差，
Figure BDA0003939007040000056
是第i，j，k网格单元第n虚拟时间步的空间通量残差，Δτ是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算，采用不同局部虚拟时间子迭代步长定常计算不同的网格单元。
可选的，所述空间通量计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量。
Figure BDA0003939007040000057
Figure BDA0003939007040000058
Figure BDA0003939007040000059
其中，下标k分别取曲线坐标系ξ，η，ζ方向之一，相应的Fk即为对应ξ，η，ζ方向的电磁通量，
Figure BDA00039390070400000510
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure BDA00039390070400000511
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；S，S-为相似矩阵，Λ+，Λ-分别为正负特征值构成的对角矩阵，QL，QR分别代表分界面处左、右状态变量，采用MUSCL格式：
Figure BDA0003939007040000061
Figure BDA0003939007040000062
其中，φ＝1是限制器，下标i是网格单元编号，
Figure BDA0003939007040000063
对应单元分界面，κ＝1/3是3阶精度格式的控制参数，
Figure BDA00039390070400000614
和Δ分别是后差和前差算符；
Figure BDA0003939007040000064
表示在网格单元
Figure BDA0003939007040000065
分界面处左状态电磁场守恒变量，
Figure BDA0003939007040000066
表示在网格单元
Figure BDA0003939007040000067
分界面处右状态电磁场守恒变量；Qi是第i个网格单元散射电磁场守恒变量，Qi+1是第i+1个网格单元散射电磁场守恒变量。
可选的，在介质/导体边界，使用物理和数值边界条件构造虚拟像点扩充层电磁场值；
其中，使用的物理和数值边界条件为：
Figure BDA0003939007040000068
Figure BDA0003939007040000069
Figure BDA00039390070400000610
Figure BDA00039390070400000611
其中，
Figure BDA00039390070400000612
是导体壁面单位外法矢，Et是导体壁面总电场强度矢量，Bt是导体壁面总磁感应场矢量，Etn是导体壁面总电场强度法向分量，Ht是导体壁面总磁场强度矢量；
构造的虚拟像点扩充层电磁场值为：
Figure BDA00039390070400000613
Figure BDA0003939007040000071
Ex-2 t＝Ex-1 t
Hx-2 t＝Hx-1 t
其中，网格块内单元编号为x，两层扩充层相应编号为x-1和x-2，上标t代表总场。
可选的，对含不同波阻相邻单元根据以下介质之间物理边界条件作修正：
Figure BDA0003939007040000072
Figure BDA0003939007040000073
Figure BDA0003939007040000074
Figure BDA0003939007040000075
其中，上标s，t分别代表散射场和总场，下标L，R分别代表不同波阻分界面左右状态，
Figure BDA0003939007040000076
是介质分界面单位法矢，不同波阻介质分界面MUSCL插值具体修正方式为：分别对相邻电场强度和磁场强度切向分量作MUSCL插值；分别对相邻电位移矢量和磁感应强度矢量的法向分量作MUSCL插值；最后结合两个MUSCL插值结果得到不同波阻介质分界面的
Figure BDA0003939007040000077
可选的，采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure BDA0003939007040000078
其中，A+，A-，B+，B-，C+，C-是分裂后的系数矩阵，RHS是上一个物理时间步计算的空间通量残差，ΔWn是隐式子迭代电磁场差值；将
Figure BDA0003939007040000081
表示为LDU，近似因式分解得到：
Figure BDA0003939007040000082
Figure BDA0003939007040000083
Figure BDA0003939007040000084
Figure BDA0003939007040000085
其中，下标i，j，k是网格单元编号，β是雅可比系数矩阵最大特征值分裂参数，γA，γB，γC是雅克比系数矩阵最大特征值，I是单位对角矩阵，D是如上公式所计算的对角矩阵，D+是对应负特征值分裂后所形成上三角矩阵，D-是对应正特征值分裂后所形成下三角矩阵，ΔW+是分裂所对应上三角部分电磁守恒变量差值，ΔW-分裂所对应下三角部分电磁守恒变量差值；ΔWi-1指网格单元i-1的相邻迭代时间步电磁守恒变量差值；ΔWj-1指网格单元j-1的相邻迭代时间步电磁守恒变量差值；ΔWk-1指网格单元k-1的相邻迭代时间步电磁守恒变量差值；ΔWi+1指网格单元i+1的相邻迭代时间步电磁守恒变量差值；ΔWj+1指网格单元j+1的相邻迭代时间步电磁守恒变量差值；ΔWk+1指网格单元k+1的相邻迭代时间步电磁守恒变量差值；
Figure BDA0003939007040000086
是指相邻i-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000087
是指相邻j-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000088
是指相邻k-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000089
是指相邻i+1网格单元分裂后的系数矩阵；
Figure BDA00039390070400000810
是指相邻j+1网格单元分裂后的系数矩阵；
Figure BDA00039390070400000811
是指相邻k+1网格单元分裂后的系数矩阵；
得到前后向迭代计算电磁场差值ΔW：
Figure BDA0003939007040000091
L＝D+D-
U＝D+D+
其中，D-1是D的逆矩阵，L，U分别是根据D+D-、D+D+计算所得到的上三角矩阵、下三角矩阵；
前向循环：
Figure BDA0003939007040000092
后向循环：
Figure BDA0003939007040000093
其中,
Figure BDA0003939007040000094
是电磁守恒变量差值的中间过渡变量；
最后由介质内电磁场守恒变量减去解析的入射电磁场量得到散射电磁场Q:
Figure BDA0003939007040000095
其中,μr是相对磁化率，εr是相对介电常数，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量。
第二方面，本申请公开了一种介质涂覆导体复合目标电磁散射隐式计算装置，包括：
数据获取及初始化模块，用于获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；
计算模块，用于基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
第三方面，本申请公开了一种电子设备，包括存储器和处理器，其中：
所述存储器，用于保存计算机程序；
所述处理器，用于执行所述计算机程序，以实现前述的介质涂覆导体复合目标电磁散射隐式计算方法。
第四方面，本申请公开了一种计算机可读存储介质，用于保存计算机程序，其中，所述计算机程序被处理器执行时实现前述的介质涂覆导体复合目标电磁散射隐式计算方法。
可见，本申请获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。也即，本申请通过一种介质内守恒电磁场两重(物理和虚拟)时间迭代推进和空间通量残差隐式计算的时域有限体积方法，来替代现有方法中通常采用的时间、空间通量残差都为显式的求解算法。如此一来，使得物理时间步长根据物理问题进行选取而不受稳定性的限制，稳定性则由隐式虚拟时间子迭代满足，从而克服显式方法时间步受稳定性限制且必须采用统一最小全局步长和加密网格带来的大计算量问题，提高计算效率。最后，利用这种双时间步隐式时域有限体积方法高效地求解2维、3维麦克斯韦方程组，获得介质涂覆导体复合目标时变电磁场的时间、空间分布，在保证格式精度的同时提高时间推进效率，节约计算成本。并且，本申请在介质/导体边界，使用四个物理和数值边界条件构造虚拟像点电磁场值，作为网格块端面扩充层，使得在块边界插值处所需电磁场值依然存在，从而使块内所有单元分界面通量计算，即使在网格块端面也不用特殊处理，都保持一致的计算处理方式，有利计算效率。并且，插值遵循了间断电磁场的物理特性，能更准确计算电磁场时空分布。这样，通过本申请提供的方案，能够在保障精度的同时提升介质涂覆导体复合目标电磁散射问题的计算效率，从而节约成本。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本申请公开的一种介质涂覆导体复合目标电磁散射隐式计算方法流程图；
图2为本申请公开的一种具体的介质涂覆导体复合目标电磁散射隐式计算流程图；
图3为本申请公开的一种介质涂覆厚度λd/30，εr＝2.56，金属球E平面双站雷达散射截面分布示意图；
图4为本申请公开的一种介质涂层厚度0.01米，金属球(ka＝3)双站雷达散射截面分布示意图；
图5为本申请公开的一种涂覆金属行波管电磁散射示意图；
图6为本申请公开的一种行波管有无介质涂覆同一空间点电磁场分量比较示意图；
图7为本申请公开的一种行波管无介质涂覆表面诱导电流密度分布示意图；
图8为本申请公开的一种行波管有介质涂覆表面诱导电流密度分布示意图；
图9为本申请公开的一种行波管介质涂覆前后双站雷达散射截面比较示意图；
图10为本申请公开的一种介质涂覆导体复合目标电磁散射隐式计算装置结构示意图；
图11为本申请公开的一种电子设备结构图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
参见图1所示，本申请实施例公开了一种介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，包括：
步骤S11：获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型。
在一种具体的实施方式中，获取网格数据文件的具体过程为采用四边形或六面体结构对仿真模型进行网格剖分，网格在散射壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，得到网格数据文件。
在一种实施方式中，可以用户可以利用仿真软件，根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模，得到仿真模型，然后采用四边形或六面体结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件。本申请获取用户输入的目标计算电磁参数、数值计算控制参数、以及输入网格数据和边界条件信息文件，初始化计算空间电磁场。
并且，网格密度保证每波长15-20个网格点，壁面密度>300点/介质内当地波长，几何奇异处加密到50-100个网格点/介质内当地波长；二维网格时，则在垂直该二维网格所在平面按右手法则推进1层，作为三维问题特例统一计算；网格数据文件包括结构网格块数目及每块3个曲线坐标系下维度。
进一步的，对于输入的目标计算电磁参数、数值计算控制参数：预先设定物理时间步长，二维问题用入射电磁波周期无量纲化0.01量级的物理时间步长来保证足够计算精度，三维问题选择量级0.001的物理时间步长；同时设定子迭代收敛标准判据值为全网格空间相邻子迭代时间步最大幅度差值<0.001，判定为收敛；设定最大子迭代步数为30。
步骤S12：基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数。该介质麦克斯韦方程组为待求解散射场形式、含介质时域麦克斯韦方程组。
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
其中，所述物理时间步循环和虚拟时间步子迭代循环过程为：
增加虚拟时间导数项
Figure BDA0003939007040000131
将待求解散射场形式、含介质时域麦克斯韦方程组修正为：
Figure BDA0003939007040000141
Figure BDA0003939007040000142
Figure BDA0003939007040000143
Figure BDA0003939007040000144
其中，W是电磁场守恒变量，μr是相对磁化率，μ0是真空磁化率，εr是相对介电常数，ε0是真空介电常数，σm是磁导率，σe是电导率。Q是散射电磁场原始变量，τ是虚拟时间，t是归一化物理时间，Fx、Fy、Fz是直角坐标系下电磁通量的X、y、z分量，S是介质磁导率、电导率等引发的外加源项，J是外加强迫电流矢量，Jx、Jy、Jz是直角坐标系下外加强迫电流的X、y、z分量，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量，B是实数型散射磁感应强度矢量，D是实数型散射场电位移矢量，Ei是入射场电场强度矢量，Hi是入射场磁场强度矢量，E是散射场电场强度矢量，H是散射场磁场强度矢量，含下标X、y、z标量分别是对应矢量的直角坐标系X、y、z分量，上标t对应的是电磁场总场形式；
当
Figure BDA0003939007040000145
收敛时，该方程组等同于原始方程组，定常虚拟时间τ子迭代表示为：
Figure BDA0003939007040000146
Figure BDA0003939007040000151
其中，W*是W经虚拟时间子迭代后的电磁场守恒变量值，W*是Wn+1的近似，*号上标代表的是虚拟时间子迭代对应守恒变量值，R是通量残差，R*是通量残差R增加了物理时间导数项以及源项后的残差；n是物理时间步数；Δt是物理时间步长，Wn是第n物理时间步时的电磁守恒变量，Wn-1是第n-1物理时间步时的电磁守恒变量；R*(W*)为W*对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度；定常子迭代部分采用隐式算法：
Figure BDA0003939007040000152
其中,ξ、η、ζ分别表示结构网格曲线坐标系的三个不同方向；F,G,H是分别对应曲线坐标系ξ，η，ζ方向电磁通量；
Figure BDA0003939007040000153
是对应曲线坐标系源项；ω是隐式控制参数，取ω＝1的全隐式，其他参数(即ω不为1)对应显、隐混合格式；下标i，j，k是网格单元编号，
Figure BDA0003939007040000154
是第i，j，k网格单元第n+1虚拟时间步的电磁守恒变量，
Figure BDA0003939007040000155
是第i，j，k网格单元第n虚拟时间步的电磁守恒变量，Wn+1是第n+1物理时间步的电磁守恒变量，Wn是第n物理时间步的电磁守恒变量，Wn-1是第n-1物理时间步的电磁守恒变量，
Figure BDA0003939007040000156
是第i，j，k网格单元第n+1虚拟时间步的空间通量残差，
Figure BDA0003939007040000157
是第i，j，k网格单元第n虚拟时间步的空间通量残差，Δτ是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算，显著区别于显式方法，定常计算的不同的网格单元有不同局部虚拟时间子迭代步长从而加快该单元电磁场收敛。
进一步，所述空间通量计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量。
Figure BDA0003939007040000161
Figure BDA0003939007040000162
Figure BDA0003939007040000163
其中，下标k分别取曲线坐标系ξ，η，ζ方向之一，相应的Fk即为对应ξ，η，ζ方向的电磁通量，
Figure BDA0003939007040000164
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure BDA0003939007040000165
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；S，S-为相似矩阵，Λ+，Λ-分别为正负特征值构成的对角矩阵，QL，QR分别代表分界面处左、右状态变量，采用MUSCL格式而达到最高三阶精度：
Figure BDA0003939007040000166
Figure BDA0003939007040000167
其中，φ＝1是限制器，下标i是网格单元编号，
Figure BDA0003939007040000168
对应单元分界面，κ＝1/3是3阶精度格式的控制参数，
Figure BDA0003939007040000169
和Δ分别是后差和前差算符；
Figure BDA00039390070400001610
表示在网格单元
Figure BDA00039390070400001611
分界面处左状态电磁场守恒变量，
Figure BDA00039390070400001612
表示在网格单元
Figure BDA00039390070400001613
分界面处右状态电磁场守恒变量；Qi是第i个网格单元散射电磁场守恒变量，Qi+1是第i+1个网格单元散射电磁场守恒变量。
Figure BDA00039390070400001614
插值过程中需要相邻单元中心电磁场值，在网格块端面编号会溢出真实区间，因此在每个网格块每个维度两个端面方向，根据网格块之间几何关系、物理边界条件扩充两层，使得在块边界插值处所需Qi-1、Qi+1依然存在，从而使块内所有单元分界面通量计算，即使在网格块端面也不用特殊处理，都保持一致的计算处理方式，有利计算效率。
进一步的，在介质/导体边界，使用物理和数值边界条件构造虚拟像点扩充层电磁场值；
其中，使用的物理和数值边界条件为：
Figure BDA0003939007040000171
Figure BDA0003939007040000172
Figure BDA0003939007040000173
Figure BDA0003939007040000174
其中，
Figure BDA0003939007040000175
是导体壁面单位外法矢，Et是导体壁面总电场强度矢量，Bt是导体壁面总磁感应场矢量，Etn是导体壁面总电场强度法向分量，Ht是导体壁面总磁场强度矢量；
构造的虚拟像点扩充层电磁场值为：
Figure BDA0003939007040000176
Figure BDA0003939007040000177
Ex-2 t＝Ex-1 t
Hx-2 t＝Hx-1 t
其中，网格块内单元编号为x，两层扩充层相应编号为x-1和x-2，上标t代表总场。
例如，结合这四个物理及数值介质/导体边界条件，构造虚拟像点扩充层电磁场值，假设网格块内单元编号为2，那么两层扩充层相应编号为1和0：
Figure BDA0003939007040000178
Figure BDA0003939007040000179
E0 t＝E1 t
H0 t＝H1 t
其中，下标0,1,2分别代表对应编号的场量，上标t代表总场。
上述MUSCL格式计算
Figure BDA0003939007040000181
过程中要求电磁场Qi、Qi+1连续可微才能保持数值精度，但是不同波阻介质单元之间存在间断，非连续可微，按空间方向随意插值会带来大的数值误差。因此本申请实施例对含不同电磁参数(波阻)相邻单元根据以下介质之间物理边界条件作修正：
Figure BDA0003939007040000182
Figure BDA0003939007040000183
Figure BDA0003939007040000184
Figure BDA0003939007040000185
其中，上标s，t分别代表散射场和总场，下标L，R分别代表不同波阻分界面左右状态，
Figure BDA0003939007040000186
是介质分界面单位法矢，不同波阻介质分界面MUSCL插值具体修正方式为：分别对相邻电场强度和磁场强度切向分量作MUSCL插值；分别对相邻电位移矢量和磁感应强度矢量的法向分量作MUSCL插值；最后结合两个MUSCL插值结果得到不同波阻介质分界面的
Figure BDA0003939007040000187
进一步的，隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值的具体过程包括：
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure BDA0003939007040000188
其中，A+，A-，B+，B-，C+，C-是分裂后的系数矩阵，RHS是上一个物理时间步计算的空间通量残差，ΔWn是隐式子迭代电磁场差值；将
Figure BDA0003939007040000191
表示为LDU(即上三角、下三角、对角线)，近似因式分解得到：
Figure BDA0003939007040000192
Figure BDA0003939007040000193
Figure BDA0003939007040000194
Figure BDA0003939007040000195
其中，下标i，j，k是网格单元编号，β是雅可比系数矩阵最大特征值分裂参数，γA，γB，γC是雅克比系数矩阵最大特征值，I是单位对角矩阵，D是如上公式所计算的对角矩阵，D+是对应负特征值分裂后所形成上三角矩阵，D-是对应正特征值分裂后所形成下三角矩阵，ΔW+是分裂所对应上三角部分电磁守恒变量差值，ΔW-分裂所对应下三角部分电磁守恒变量差值；ΔWi-1指网格单元i-1的相邻迭代时间步电磁守恒变量差值；ΔWj-1指网格单元j-1的相邻迭代时间步电磁守恒变量差值；ΔWk-1指网格单元k-1的相邻迭代时间步电磁守恒变量差值；ΔWi+1指网格单元i+1的相邻迭代时间步电磁守恒变量差值；ΔWj+1指网格单元j+1的相邻迭代时间步电磁守恒变量差值；ΔWk+1指网格单元k+1的相邻迭代时间步电磁守恒变量差值；
Figure BDA0003939007040000196
是指相邻i-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000197
是指相邻j-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000198
是指相邻k-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000199
是指相邻i+1网格单元分裂后的系数矩阵；
Figure BDA00039390070400001910
是指相邻j+1网格单元分裂后的系数矩阵；
Figure BDA00039390070400001911
是指相邻k+1网格单元分裂后的系数矩阵；
得到前后向迭代计算电磁场差值ΔW：
Figure BDA0003939007040000201
L＝D+D-
U＝D+D+
其中，D-1是D的逆矩阵，L，U分别是根据D+D-、D+D+计算所得到的上三角矩阵、下三角矩阵；
前向循环：
Figure BDA0003939007040000202
后向循环：
Figure BDA0003939007040000203
其中,
Figure BDA0003939007040000204
是电磁守恒变量差值的中间过渡变量；
最后由介质内电磁场守恒变量减去解析的入射电磁场量得到散射电磁场Q:
Figure BDA0003939007040000205
其中,μr是相对磁化率，εr是相对介电常数，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量。
并且，在具体的实施方式中，子迭代收敛判据为电场和磁场在计算空间最大幅度差的绝对值，且数值计算中设置有最大子迭代步数。
可以理解的是，目前需要一种高效的直接求解时域Maxwell方程组，在保持高数值精度的同时，既放宽很小网格尺度对迭代物理时间步的限制，又能准确模拟介质/导体、介质/介质电磁反射、折射物理特性的数值技术，提升计算性能。为实现高效时域电磁数值方法求解复杂外形、高频、介质涂覆导体复合目标大规模电磁散射问题，在保证时间、空间计算格式高精度的前提下，提升计算效率。本申请提供了一种介质涂覆导体复合目标电磁散射问题的隐式时域数值技术，涉及隐式双时间步时间推进、介质导体物理边界条件虚拟像点处理、电磁场间断所对应介质分界面处，电磁通量根据物理条件分解计算的电磁数值方法。首先通过在控制方程中引入定常的虚拟时间导数项，在每个物理时间段内对物理时间导数作线性化处理，从而使得物理时间推进步长可以根据物理问题进行选取而不受稳定性的限制，计算的稳定性要求由虚拟时间子迭代满足，子迭代的定常计算可以采用局部时间步长加速收敛，极大缩短获取介质内稳定电磁场所需计算时间。其次，通过每个网格块每个维度方向扩充两层，来统一处理虚拟像点介质/导体物理边界，达到不需特殊处理网格块端点通量，与网格内部单元一样计算方式的目的，从而提高程序运行和并行计算效率。最后，根据介质分界面电磁场间断特性构造通量，反映物理特征能更准确计算电磁场时空分布。利用这种双时间步隐式时域有限体积方法高效地求解2维、3维Maxwell方程组，获得介质/导体复合目标电磁散射场时间、空间分布，保证格式精度同时提高时间推进效率，节约计算成本。
与现有技术相比，本申请的有益效果主要包括：
1、本申请通过一种介质内守恒电磁场两重(物理和虚拟)时间迭代推进和空间通量残差隐式计算的时域有限体积方法，来替代现有方法中通常采用的时间、空间通量残差都为显式的求解算法。如此一来，使得物理时间步长根据物理问题进行选取而不受稳定性的限制
2、本申请采用介质内散射场和部分入射场构成的守恒电磁场方式，其中入射场由解析给出，相当于求解散射场形式的Maxwell方程组，由此不用在整个计算网格空间数值计算入射电磁波的传播，避免了入射电磁场的耗散及色散，有利保持数值精度。
3、本申请支持结构网格和多区域分解，采用贴体网格直接可将麦克斯韦方程组守恒律的积分形式应用到离散的曲线坐标系网格单元。
4、本申请不同于传统FDTD，FDTD采用笛卡尔直角正交网格模拟壁面存在阶梯效应影响数值精度，并用电磁场分量时空交叉放置添加人工粘性给2阶中心差分格式，FVTD采用贴体曲线坐标系网格能更好拟合物面和在几何奇异处加密网格，并且电磁场量在网格空间同置于网格单元中心，采用迎风格式来保持人工粘性，更有利保持精度和算法设计。
5、本申请在每个网格块每个维度两个端面方向，根据网格块之间几何关系、物理边界条件扩充两层，使得在块边界插值处所需电磁场值依然存在，从而使块内所有单元分界面通量计算，即使在网格块端面也不用特殊处理，都保持一致的计算处理方式，有利计算效率。需要指出的是如果没有扩充网格，那么在介质/导体边界，每个时间步都必须特殊处理，给场量插值以及信息交换和并行带来不便，影响计算效率。
6、本申请在介质/导体边界，使用四个物理和数值边界条件构造虚拟像点电磁场值，作为网格块端面扩充层。
7、本申请在不同波阻介质/介质网格单元分界面，即含不同电磁参数(波阻)相邻单元需要根据介质之间电磁波反射、折射物理条件对MUSCL插值作修正：分别对相邻电场强度和磁场强度切向分量作MUSCL插值，分别对相邻电位移矢量和磁感应强度矢量的法向分量作MUSCL插值，最后结合两者得到不同波阻介质分界面的左右状态变量，如此插值遵循了间断电磁场的物理特性，能更准确计算电磁场时空分布。需要指出的是，由于不同电磁参数网格单元之间电磁场具有间断特性，计算电磁通量所需分界面电磁场值(由相邻单元中心场值插值得到)，如果不考虑物理特性纯粹像自由空间一样，按空间方向插值，将不能准确模拟介质之间物理特征，带来大的数值误差。
8、本申请采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解,用两次循环替代了稀疏矩阵求逆，工程上简单易用。
9、本申请替代传统Runge-Kutta时间显式推进及空间通量也显式计算的时域有限体积方法，从而放松网格和显式算法对物理时间步极端限制，这种双时间步隐式时域有限体积方法能高效地求解2维、3维Maxwell方程组，获得介质涂覆导体复合目标时变电磁场的时空分布，在保证格式精度的同时提高时间推进效率，提升计算性能。
进一步的，例如参见图2所示，图2为本申请实施例提供的一种具体的介质涂覆导体复合目标电磁散射隐式计算流程图。图2中，整个时域有限体积方法计算电磁场软件按结构可以分为：预处理、电磁场计算和后处理三个部分。预处理主要包括网格数据输入、计算参数数据输入、控制参数输入三个模块，主要用来读入网格数据、计算参数数据输入、控制参数文件，并在此基础上，进行预先处理，为电磁场计算提供计算支撑；电磁场计算包括：空间电磁场MUSCL格式插值、单元分界面通量计算、时间推进、收敛判断模块组成；后处理主要用于输出电磁场的时间、空间分布、目标表面诱导电流密度、雷达散射截面输出。
下面结合要数值模拟的Maxwell方程组两个旋度方程：
Figure BDA0003939007040000231
Figure BDA0003939007040000232
其中，
Figure BDA0003939007040000233
是梯度符号，Bt是实数型总场磁感应强度矢量，Dt是实数型总场电位移矢量，Et是实数型总场电场强度矢量，Ht是实数型总场磁场强度矢量，J是外加强迫电流矢量，Jm是磁化电流矢量，Je是传导电流矢量，t是归一化物理时间。
介绍双时间步时域有限体积法迭代数值计算介质涂覆导体复合目标电磁散射过程。散射场形式、含介质时域麦克斯韦方程组直角坐标系下表示为：
Figure BDA0003939007040000234
Figure BDA0003939007040000235
Figure BDA0003939007040000236
Figure BDA0003939007040000237
其中，W是电磁场守恒变量，μr是相对磁化率，μ0是真空磁化率，εr是相对介电常数，ε0是真空介电常数，σm是磁导率，σe是电导率。Q是散射电磁场原始变量，τ是虚拟时间，t是归一化物理时间，Fx、Fy、Fz是直角坐标系下电磁通量的X、y、z分量，S是介质磁导率、电导率等引发的外加源项，J是外加强迫电流矢量，Jx、Jy、Jz是直角坐标系下外加强迫电流的X、y、z分量，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量，B是实数型散射磁感应强度矢量，D是实数型散射场电位移矢量，Ei是入射场电场强度矢量，Hi是入射场磁场强度矢量，E是散射场电场强度矢量，H是散射场磁场强度矢量，含下标X、y、z标量分别是对应矢量的直角坐标系X、y、z分量，上标t对应的是电磁场总场形式；
对于复杂外形物体，采用的是计算空间贴体多块结构网格，因此均存在坐标变换：
k＝k(x，y，z) k＝ξ，η，ζ
其中，k代表曲线坐标系ξ，η，ζ三个方向，分别取ξ，η，ζ之一。得到所要数值模拟的曲线坐标系下的麦克斯韦方程组守恒形：
Figure BDA0003939007040000241
Figure BDA0003939007040000242
Figure BDA0003939007040000243
Figure BDA0003939007040000244
Figure BDA0003939007040000245
Figure BDA0003939007040000246
其中，V是坐标变换的雅可比矩阵，^上标变量代表曲线坐标系下的值，由坐标变换获取。
Figure BDA0003939007040000247
为曲线坐标系守恒变量，
Figure BDA0003939007040000248
即为当k分别取曲线坐标系方向ξ，η，ζ之一时的曲线坐标系下电磁通量。
为了摆脱传统显式时域有限体积全局最小物理时间步的限制带来的大计算量弊端，为此发明提高介质/导体复合目标时变电磁场时间推进效率的全隐式双时间步计算方法，用双时间步保证时间精度的同时根据物理特征选取物理时间步，结合隐式空间通量残差获得稳定、高效的计算流程，包括以下步骤：
步骤1：根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模。
步骤2：采用四边形(二维)或六面体结构(三维)对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，所述网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件。网格密度保证每波长15-20个网格点，壁面密度>300点/每波长，几何奇异处加密到50-100个网格点/每波长，二维网格在垂直该平面按右手法则推进1层，作为3维问题特例统一计算。网格数据文件包括结构网格块数目及每块3个曲线坐标系下维度。
步骤3：预处理部分，输入目标计算电磁参数、数值计算控制参数。虚拟时间子迭代因是隐式CFL数不受上述显式稳定性要求约束，预先设定物理时间步长，二维问题用入射电磁波周期无量纲化的物理时间步长Δt＝0.01计算精度已可以，三维问题选择Δt≈0.001量级。同时设定子迭代收敛标准判据值(例如全网格空间相邻子迭代时间步最大幅度差值<0.001,判定为收敛)、设定最大子迭代步数(例：isubmax＝30)。
步骤4:输入网格数据和边界条件信息文件，初始化计算空间电磁场。
步骤5：以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解。
步骤5-1：外层物理时间步循环，到计算收敛结束。
步骤5-2：内层虚拟时间步子迭代循环，到子迭代收敛或最大子迭代步数条件满足之一。双时间方法非定常计算的时间精度还受到每一物理时间步上子迭代步数制约。迭代步数多，收敛控制严格，时间精度才有保证。真实时间步小，则子迭代可以取得少一些，数值计算中为防止场变化剧烈处、以及子迭代CFL数取得过大等情况时，残差很有可能无法下降到收，陷入死循环而给出最大子迭代步数，另一方面对子迭代收敛进行判断和控制,在保证一定时间精度前提下尽快结束迭代，子迭代收敛判据使用电场和磁场在计算空间最大幅度差绝对值。
以下介绍步骤5-1和步骤5-2具体算法：采用MUSCL(Monotonic Upstream Schemesfor Conservation Laws,MUSCL)插值。
增加虚拟时间导数项
Figure BDA0003939007040000261
将待求解散射场形式、含介质时域麦克斯韦方程组修正为：
Figure BDA0003939007040000262
明显可见，当
Figure BDA0003939007040000263
收敛时，该方程组等同于原始方程组，定常虚拟时间τ子迭代表示为：
Figure BDA0003939007040000264
Figure BDA0003939007040000265
其中，W*是W经虚拟时间子迭代后的电磁场守恒变量值，W*是Wn+1的近似，*号上标代表的是虚拟时间子迭代对应守恒变量值，R是通量残差，R*是通量残差R增加了物理时间导数项以及源项后的残差；n是物理时间步数；Δt是物理时间步长，Wn是第n物理时间步时的电磁守恒变量，Wn-1是第n-1物理时间步时的电磁守恒变量；R*(W*)为W*对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度。
定常子迭代部分采用隐式算法：
Figure BDA0003939007040000266
其中,ξ为结构网格曲线坐标系方向1，η为结构网格曲线坐标系方向2，ζ为结构网格曲线坐标系方向3；F,G,H是分别对应曲线坐标系ξ，η，ζ方向电磁通量；
Figure BDA0003939007040000267
是对应曲线坐标系源项；ω是隐式控制参数，取ω＝1的全隐式，其他参数(即ω不为1)对应显、隐混合格式；下标i，j，k是网格单元编号，
Figure BDA0003939007040000268
是第i，j，k网格单元第n+1虚拟时间步的电磁守恒变量，
Figure BDA0003939007040000271
是第i，j，k网格单元第n虚拟时间步的电磁守恒变量，Wn+1是第n+1物理时间步的电磁守恒变量，Wn是第n物理时间步的电磁守恒变量，Wn-1是第n-1物理时间步的电磁守恒变量，
Figure BDA0003939007040000272
是第i，j，k网格单元第n+1虚拟时间步的空间通量残差，
Figure BDA0003939007040000273
是第i，j，k网格单元第n虚拟时间步的空间通量残差，Δτ是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算，显著区别于显式方法，定常计算的不同的网格单元有不同局部虚拟时间子迭代步长从而加快该单元电磁场收敛。
步骤5-3：在每个虚拟时间子迭代过程中，逐网格块、逐网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
有限体积法的空间精度体现在能否精确模拟守恒变量W在网格单元分界面处的状态变量，以得到相应精确的分界面通量F,G,H，采用Steger-Warming分裂计算网格单元分界面通量。
Figure BDA0003939007040000274
Figure BDA0003939007040000275
Figure BDA0003939007040000276
其中，下标k分别取曲线坐标系ξ，η，ζ方向之一，相应的Fk即为对应ξ，η，ζ方向的电磁通量F,G,H，
Figure BDA0003939007040000277
代表曲线坐标系ξ，η，ζ，对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure BDA0003939007040000278
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；S，S-为相似矩阵，Λ+，Λ-分别为正负特征值构成的对角矩阵，QL，QR分别代表分界面处左、右状态变量，采用MUSCL格式而达到最高三阶精度：
Figure BDA0003939007040000279
Figure BDA00039390070400002710
其中，φ＝1是限制器，下标i是网格单元编号，
Figure BDA00039390070400002711
对应单元分界面，κ＝1/3是3阶精度格式的控制参数，
Figure BDA0003939007040000281
和Δ分别是后差和前差算符；
Figure BDA0003939007040000282
表示在网格单元
Figure BDA0003939007040000283
分界面处左状态电磁场守恒变量，
Figure BDA0003939007040000284
表示在网格单元
Figure BDA0003939007040000285
分界面处右状态电磁场守恒变量；Qi是第i个网格单元散射电磁场守恒变量，Qi+1是第i+1个网格单元散射电磁场守恒变量。
步骤5-3-1：为MUSCL插值作预备，扩充每个网格块，特别是介质/导体边界构造虚拟像点：
Figure BDA0003939007040000286
插值过程中需要相邻单元中心电磁场值，在网格块端面编号会溢出真实区间，因此在每个网格块每个维度两个端面方向，根据网格块之间几何关系、物理边界条件扩充两层，使得在块边界插值处所需Qi-1、Qi+1依然存在，从而使块内所有单元分界面通量计算，即使在网格块端面也不用特殊处理，都保持一致的计算处理方式，有利计算效率。在介质/导体边界，使用物理和数值边界条件：
Figure BDA0003939007040000287
Figure BDA0003939007040000288
Figure BDA0003939007040000289
Figure BDA00039390070400002810
其中，
Figure BDA00039390070400002811
是导体壁面单位外法矢，Et是导体壁面总电场强度矢量，Bt是导体壁面总磁感应场矢量，Etn是导体壁面总电场强度法向分量，Ht是导体壁面总磁场强度矢量。
结合这四个物理及数值介质/导体边界条件，构造虚拟像点扩充层电磁场值(假设网格块内单元编号为2，那么两层扩充层相应编号为1和0)：
Figure BDA00039390070400002812
Figure BDA0003939007040000291
E0 t＝E1 t
H0 t＝H1 t
其中，下标0,1,2分别对应场量编号，上标t代表总场。
步骤5-3-2：不同波阻介质/介质单元之间，电磁场存在间断，为更好模拟在介质之间电磁波反射、折射特性，在不同波阻介质分界面作特殊MUSCL插值：
上述MUSCL格式计算
Figure BDA0003939007040000292
过程中要求电磁场Qi、Qi+1连续可微才能保持数值精度，但是不同波阻介质单元之间存在间断，非连续可微，按空间方向随意插值会带来大的数值误差。因此含不同电磁参数(波阻)相邻单元需要根据以下介质之间物理边界条件作修正：
Figure BDA0003939007040000293
Figure BDA0003939007040000294
Figure BDA0003939007040000295
Figure BDA0003939007040000296
其中，上标s，t分别代表散射场和总场，下标L，R分别代表不同波阻分界面左右状态，
Figure BDA0003939007040000297
是介质分界面单位法矢，不同波阻介质分界面MUSCL插值具体修正方式为：分别对相邻电场强度和磁场强度切向分量作MUSCL插值；分别对相邻电位移矢量和磁感应强度矢量的法向分量作MUSCL插值；最后结合两者得到不同波阻介质分界面的
Figure BDA0003939007040000298
空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解,用两次循环替代了稀疏矩阵求逆，工程上简单易用。经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到：
Figure BDA0003939007040000299
其中，A+，A-，B+，B-，C+，C-是分裂后的系数矩阵，RHS是上一个物理时间步计算的空间通量残差，ΔWn是隐式子迭代电磁场差值；表示为LDU，近似因式分解得到：
Figure BDA0003939007040000301
Figure BDA0003939007040000302
Figure BDA0003939007040000303
Figure BDA0003939007040000304
其中，下标i，j，k是网格单元编号，β是雅可比系数矩阵最大特征值分裂参数，γA，γB，γC是雅克比系数矩阵最大特征值，I是单位对角矩阵，D是如上公式所计算的对角矩阵，D+是对应负特征值分裂后所形成上三角矩阵，D-是对应正特征值分裂后所形成下三角矩阵，ΔW+是分裂所对应上三角部分电磁守恒变量差值，ΔW-分裂所对应下三角部分电磁守恒变量差值；ΔWi-1指网格单元i-1的相邻迭代时间步电磁守恒变量差值；ΔWj-1指网格单元j-1的相邻迭代时间步电磁守恒变量差值；ΔWk-1指网格单元k-1的相邻迭代时间步电磁守恒变量差值；ΔWi+1指网格单元i+1的相邻迭代时间步电磁守恒变量差值；ΔWj+1指网格单元j+1的相邻迭代时间步电磁守恒变量差值；ΔWk+1指网格单元k+1的相邻迭代时间步电磁守恒变量差值；
Figure BDA0003939007040000305
是指相邻i-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000306
是指相邻j-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000307
是指相邻k-1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000308
是指相邻i+1网格单元分裂后的系数矩阵；
Figure BDA0003939007040000309
是指相邻j+1网格单元分裂后的系数矩阵；
Figure BDA00039390070400003010
是指相邻k+1网格单元分裂后的系数矩阵；
得到前后向迭代计算电磁场差值ΔW：
Figure BDA0003939007040000311
L＝D+D-
U＝D+D+
其中，D-1是D的逆矩阵，L，U分别是根据D+D-、D+D+计算所得到的上三角矩阵、下三角矩阵；
前向循环：
Figure BDA0003939007040000312
后向循环：
Figure BDA0003939007040000313
其中,
Figure BDA0003939007040000314
是电磁守恒变量差值的中间过渡变量。
最后由介质内电磁场守恒变量减去解析的入射电磁场量得到散射电磁场Q:
Figure BDA0003939007040000315
其中,μr是相对磁化率，εr是相对介电常数，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量。
以上是隐式双时间步计算Maxwell方程组控制介质内电磁场迭代过程。
步骤6：收敛判断，后处理过程，输出电磁场的时间、空间分布，输出表面诱导电流和雷达散射截面空间分布数据等。
如图3所示，图3是介质涂覆厚度λd/30，εr＝2.56，金属球E平面双站雷达散射截面分布，计算涂覆一层厚度为λd/30，λd为介质当地波长，相对介电常数εr＝2.56，半径等于入射波长的金属球E平面双站雷达散射截面分布，并与解析的Mie级数解作比较，二者吻合良好。其中选择dt＝0.001，计算网格为4块，介质层内有2个网格块，维度都为16x91x91，介质外自由空间也有2个网格块，维度都为52x91x91。可见介质厚度＝0.0208倍入射波长，其中物面网格选取每波长30个网格点，远场边界在3波长外平均每波长17个网格点，辐向网格壁面加密，壁面网格第一层厚度约等于自由入射电磁波长的1/500，小于当地波长1/300。
如图4所示，图4是介质涂层厚度0.01米，金属球(ka＝3)双站雷达散射截面分布，计算涂覆金属球电磁散射，其中入射电磁波波长0.32米，金属球半径0.1528米，相当于金属球电尺寸ka＝3,涂层厚度0.01米，介质相对介电常数εr＝2+j、相对磁化率μr＝1.5+0.5j，由介质层电磁参数可见不仅有电磁波的折射效应，而且相对介电常数和相对磁化率的虚部，不可避免带来电磁能量的耗损。采用两层计算网格，介质层12x73x37，外部自由空间网格层38x73x37。选择dt＝0.0018，采用介质守恒电磁变量的FVTD方法计算结果，E平面和H平面双站雷达散射截面分布都与解析的Mie级数解吻合很好。
如图5-图9所示，比较金属行波管涂覆介质前后电磁散射特性的变化，计算条件是：X波段入射电磁波f＝9.41GHz,波长31.88mm,HH极化，波矢在XY平面与+X轴成入射角45度入射，介质涂覆厚度d＝0.4mm覆盖于行波管金属壁上表面,电磁参数：相对介电常数εr＝2.56+j4。整个网格计算空间共有35个网格块，460万网格点。图5是涂覆金属行波管电磁散射示意图；图6是行波管有无介质涂覆同一空间点电磁场分量比较，可见涂覆的有耗材料吸收部分入射的电磁能量，减小了散射网格空间的电磁场强；图7是行波管无介质涂覆表面诱导电流密度分布，图8是行波管有介质涂覆表面诱导电流密度分布；图7和图8作比较，可见在同一几何位置，涂覆的有耗材料吸收部分入射的电磁能量引起表面诱导电流减少；图9是行波管介质涂覆前后双站雷达散射截面比较，图中PEC代表完全导电体，RAM为雷达吸波材料英文简写。吸波材料涂覆后，由于电磁波能量的吸收作用，同一空间位置电磁反射波相应减小，XY平面的双站RCS分布图中，主要的反射波瓣φ＝135°处，RCS从完全导电壁5.33dB下降到介质覆盖下的-1.71dB减少了7dB,其他角度RCS也有不同的减缩。
这几个数值算例表明，本申请对应的介质涂覆导体复合目标电磁散射问题的隐式时域数值技术，既能在放松时间步限制的同时，保证数值精度并提升计算效率；本申请中网格块扩充层虚拟像点模拟介质/导体物理、数值边界，不需特殊处理网格块端点通量，提高程序运行和并行效率；本申请根据介质分界面电磁场间断特性构造通量，能准确地计算介质内电磁场时空分布。综合以上技术，本申请能准确获得介质/导体复合目标电磁散射场时间、空间分布，保证格式精度同时提高时间推进效率，节约计算成本。
可见，本申请提供的方案，涉及介质内电磁场隐式双时间步时间推进、介质/导体物理边界条件虚拟像点处理、电磁场间断所对应介质分界面，电磁通量根据物理条件分解计算的电磁数值方法。通过介质内守恒电磁场两重(物理和虚拟)时间迭代推进和空间通量残差隐式计算的时域有限体积方法，来替代现有方法通常采用的时间、空间残差显式求解算法(特别当介质内电磁波长小于真空电磁波长，显式格式稳定性要求使得物理时间步比真空电磁散射问题更小，带来更长计算时间的弊端)，其中物理时间步长根据物理问题进行选取而不受稳定性的限制，稳定性则由隐式虚拟时间子迭代满足，从而克服显式方法时间步受稳定性限制且必须采用统一最小全局步长和加密网格带来的大计算量问题，提高计算效率。其二、通过每个网格块扩充两层、统一处理的介质/导体物理边界的虚拟像点构造方式，从而不需特殊处理网格块端点通量，而是与网格内部单元一样计算方式，从而提高程序运行和并行计算效率。最后，根据介质分界面电磁场间断特性构造通量，反映物理特征能更准确计算电磁场时空分布。利用这种双时间步隐式时域有限体积方法高效地求解Maxwell方程组，获得介质/导体复合目标电磁散射场时间、空间分布，保证格式精度同时提高时间推进效率，节约计算成本。可应用于吸波材料隐身特性、介质电磁特性、介质天线辐射特性计算分析。
参见图10所示，本申请实施例公开了一种介质涂覆导体复合目标电磁散射隐式计算装置，包括：
数据获取及初始化模块11，用于获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；
计算模块12，用于基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
参见图11，本申请实施例公开了一种电子设备20，包括处理器21和存储器22；其中，所述存储器22，用于保存计算机程序；所述处理器21，用于执行所述计算机程序，前述实施例公开的介质涂覆导体复合目标电磁散射隐式计算方法。
关于上述介质涂覆导体复合目标电磁散射隐式计算方法的具体过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
并且，所述存储器22作为资源存储的载体，可以是只读存储器、随机存储器、磁盘或者光盘等，存储方式可以是短暂存储或者永久存储。
另外，所述电子设备20还包括电源23、通信接口24、输入输出接口25和通信总线26；其中，所述电源23用于为所述电子设备20上的各硬件设备提供工作电压；所述通信接口24能够为所述电子设备20创建与外界设备之间的数据传输通道，其所遵循的通信协议是能够适用于本申请技术方案的任意通信协议，在此不对其进行具体限定；所述输入输出接口25，用于获取外界输入数据或向外界输出数据，其具体的接口类型可以根据具体应用需要进行选取，在此不进行具体限定。
进一步的，本申请实施例还公开了一种计算机可读存储介质，用于保存计算机程序，其中，所述计算机程序被处理器执行时实现前述实施例公开的介质涂覆导体复合目标电磁散射隐式计算方法。
关于上述介质涂覆导体复合目标电磁散射隐式计算方法的具体过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。
以上对本申请所提供的一种介质涂覆导体复合目标电磁散射隐式计算方法及装置进行了详细介绍，本文中应用了具体个例对本申请的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本申请的方法及其核心思想；同时，对于本领域的一般技术人员，依据本申请的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本申请的限制。

Claims (10)

1.一种介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，包括：
获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；
基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
2.根据权利要求1所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，获取网格数据文件，包括：
采用四边形或六面体结构对仿真模型进行网格剖分，网格在散射壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，得到网格数据文件。
3.据权利要求1所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，所述物理时间步循环和虚拟时间步子迭代循环过程为：
增加虚拟时间导数项
Figure FDA0003939007030000011
将待求解的介质时域麦克斯韦方程组修正为：
Figure FDA0003939007030000012
Figure FDA0003939007030000021
Figure FDA0003939007030000022
Figure FDA0003939007030000023
其中，W是电磁场守恒变量，μr是相对磁化率，μ0是真空磁化率，εr是相对介电常数，ε0是真空介电常数，σm是磁导率，σe是电导率，Q是散射电磁场原始变量，τ是虚拟时间，t是归一化物理时间，Fx、Fy、Fz是直角坐标系下电磁通量的X、y、z分量，S是介质磁导率、电导率引发的外加源项，J是外加强迫电流矢量，Jx、Jy、Jz是直角坐标系下外加强迫电流的X、y、z分量，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量，B是实数型散射磁感应强度矢量，D是实数型散射场电位移矢量，Ei是入射场电场强度矢量，Hi是入射场磁场强度矢量，E是散射场电场强度矢量，H是散射场磁场强度矢量，含下标X、y、z标量分别是对应矢量的直角坐标系X、y、z分量，上标t对应的是电磁场总场形式；
当
Figure FDA0003939007030000024
收敛时，该方程组等同于原始方程组，定常虚拟时间τ子迭代表示为：
Figure FDA0003939007030000025
Figure FDA0003939007030000026
其中，W*是W经虚拟时间子迭代后的电磁场守恒变量值，W*是Wn+1的近似，*号上标代表的是虚拟时间子迭代对应守恒变量值，R是通量残差，R*是通量残差R增加了物理时间导数项以及源项后的残差；n是物理时间步数；Δt是物理时间步长，Wn是第n物理时间步时的电磁守恒变量，Wn-1是第n-1物理时间步时的电磁守恒变量；R*(W*)为W*对应的中间状态通量残差；采用后向2阶差分的物理时间导数为2阶时间精度；定常子迭代部分采用隐式算法：
Figure FDA0003939007030000031
其中,ξ、η、ζ分别表示结构网格曲线坐标系的三个不同方向；F,G,H是分别对应曲线坐标系ξ，η，ζ方向电磁通量；
Figure FDA0003939007030000032
是对应曲线坐标系源项；ω是隐式控制参数，取ω＝1的全隐式；下标i，j，k是网格单元编号，
Figure FDA0003939007030000033
是第i，j，k网格单元第n+1虚拟时间步的电磁守恒变量，
Figure FDA0003939007030000034
是第i，j，k网格单元第n虚拟时间步的电磁守恒变量，Wn+1是第n+1物理时间步的电磁守恒变量，Wn是第n物理时间步的电磁守恒变量，Wn-1是第n-1物理时间步的电磁守恒变量，
Figure FDA0003939007030000035
是第i，j，k网格单元第n+1虚拟时间步的空间通量残差，
Figure FDA0003939007030000036
是第i，j，k网格单元第n虚拟时间步的空间通量残差，Δτ是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算，采用不同局部虚拟时间子迭代步长定常计算不同的网格单元。
4.根据权利要求1所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，所述空间通量计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量：
Figure FDA0003939007030000037
Figure FDA0003939007030000038
Figure FDA0003939007030000039
其中，下标k分别取曲线坐标系ξ，η，ζ方向之一，相应的Fk即为对应ξ，η，ζ方向的电磁通量，
Figure FDA0003939007030000041
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure FDA0003939007030000042
代表曲线坐标系ξ，η，ζ对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；S，S-为相似矩阵，Λ+，Λ-分别为正负特征值构成的对角矩阵，QL，QR分别代表分界面处左、右状态变量，采用MUSCL格式：
Figure FDA0003939007030000043
Figure FDA0003939007030000044
其中，φ＝1是限制器，下标i是网格单元编号，
Figure FDA0003939007030000045
对应单元分界面，κ＝1/3是3阶精度格式的控制参数，
Figure FDA0003939007030000046
和Δ分别是后差和前差算符；
Figure FDA0003939007030000047
表示在网格单元
Figure FDA0003939007030000048
分界面处左状态电磁场守恒变量，
Figure FDA0003939007030000049
表示在网格单元
Figure FDA00039390070300000410
分界面处右状态电磁场守恒变量；Qi是第i个网格单元散射电磁场守恒变量，Qi+1是第i+1个网格单元散射电磁场守恒变量。
5.根据权利要求4所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，在介质/导体边界，使用物理和数值边界条件构造虚拟像点扩充层电磁场值；
其中，使用的物理和数值边界条件为：
Figure FDA00039390070300000411
Figure FDA00039390070300000412
Figure FDA00039390070300000413
Figure FDA00039390070300000414
其中，
Figure FDA00039390070300000415
是导体壁面单位外法矢，Et是导体壁面总电场强度矢量，Bt是导体壁面总磁感应场矢量，Et n是导体壁面总电场强度法向分量，Ht是导体壁面总磁场强度矢量；
构造的虚拟像点扩充层电磁场值为：
Figure FDA0003939007030000051
Figure FDA0003939007030000052
Ex-2 t＝Ex-1 t
Hx-2 t＝Hx-1 t
其中，网格块内单元编号为x，两层扩充层相应编号为x-1和x-2，上标t代表总场。
6.根据权利要求5所述的所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，对含不同波阻相邻单元根据以下介质之间物理边界条件作修正：
Figure FDA0003939007030000053
Figure FDA0003939007030000054
Figure FDA0003939007030000055
Figure FDA0003939007030000056
其中，上标s，t分别代表散射场和总场，下标L，R分别代表不同波阻分界面左右状态，
Figure FDA0003939007030000057
是介质分界面单位法矢，不同波阻介质分界面MUSCL插值具体修正方式为：分别对相邻电场强度和磁场强度切向分量作MUSCL插值；分别对相邻电位移矢量和磁感应强度矢量的法向分量作MUSCL插值；最后结合两个MUSCL插值结果得到不同波阻介质分界面的
Figure FDA0003939007030000058
7.根据权利要求6所述的所述的介质涂覆导体复合目标电磁散射隐式计算方法，其特征在于，
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure FDA0003939007030000061
其中，A+，A-，B+，B-，C+，C-是分裂后的系数矩阵，RHS是上一个物理时间步计算的空间通量残差，ΔWn是隐式子迭代电磁场差值；将
Figure FDA0003939007030000062
表示为LDU，近似因式分解得到：
Figure FDA0003939007030000063
Figure FDA0003939007030000064
Figure FDA0003939007030000065
Figure FDA0003939007030000066
其中，下标i，j，k是网格单元编号，β是雅可比系数矩阵最大特征值分裂参数，γA，γB，θC是雅克比系数矩阵最大特征值，I是单位对角矩阵，D是如上公式所计算的对角矩阵，D+是对应负特征值分裂后所形成上三角矩阵，D-是对应正特征值分裂后所形成下三角矩阵，ΔW+是分裂所对应上三角部分电磁守恒变量差值，ΔW-分裂所对应下三角部分电磁守恒变量差值；ΔWi-1指网格单元i-1的相邻迭代时间步电磁守恒变量差值；ΔWj-1指网格单元j-1的相邻迭代时间步电磁守恒变量差值；ΔWk-1指网格单元k-1的相邻迭代时间步电磁守恒变量差值；ΔWi+1指网格单元i+1的相邻迭代时间步电磁守恒变量差值；ΔWj+1指网格单元j+1的相邻迭代时间步电磁守恒变量差值；ΔWk+1指网格单元k+1的相邻迭代时间步电磁守恒变量差值；
Figure FDA0003939007030000071
是指相邻i-1网格单元分裂后的系数矩阵；
Figure FDA0003939007030000072
是指相邻j-1网格单元分裂后的系数矩阵；
Figure FDA0003939007030000073
是指相邻k-1网格单元分裂后的系数矩阵；
Figure FDA0003939007030000074
是指相邻i+1网格单元分裂后的系数矩阵；
Figure FDA0003939007030000075
是指相邻j+1网格单元分裂后的系数矩阵；
Figure FDA0003939007030000076
是指相邻k+1网格单元分裂后的系数矩阵；
得到前后向迭代计算电磁场差值ΔW：
Figure FDA0003939007030000077
L＝D+D-
U＝D+D+
其中，D-1是D的逆矩阵，L，U分别是根据D+D-、D+D+计算所得到的上三角矩阵、下三角矩阵；
前向循环：
Figure FDA0003939007030000078
后向循环：
Figure FDA0003939007030000079
其中,
Figure FDA00039390070300000710
是电磁守恒变量差值的中间过渡变量；
最后由介质内电磁场守恒变量减去解析的入射电磁场量得到散射电磁场Q:
Figure FDA00039390070300000711
其中,μr是相对磁化率，εr是相对介电常数，Bi是实数型入射场磁感应强度矢量，Di是实数型入射场电位移矢量。
8.一种介质涂覆导体复合目标电磁散射隐式计算装置，其特征在于，包括：
数据获取及初始化模块，用于获取目标数据，并初始化计算空间电磁场；其中，所述目标数据包括网格数据文件、边界条件文件、目标计算电磁参数以及数值计算控制参数；所述网格数据文件为对仿真模型进行网格剖分得到的网格数据文件，所述仿真模型为根据介质涂覆导体复合目标所仿真电磁问题的物理背景以及边界条件信息进行仿真建模得到的模型；
计算模块，用于基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解，得到所述介质涂覆导体复合目标的电磁场参数；
其中，所述基于所述目标数据，并以基于介质内守恒电磁场时间迭代推进和空间通量残差的隐式双时间步方式，对介质麦克斯韦方程组时变电磁场进行迭代求解包括：在仿真模型的外层采用物理时间步循环，直至计算收敛结束；在仿真模型的内层采用虚拟时间步子迭代循环，直至子迭代收敛结束；在每个虚拟时间子迭代过程中，依次对各个网块格中的各个边缘网格单元根据边界条件进行介质/导体边界虚拟像点扩充，对各个网块格中的任意网格单元根据间断电磁特性进行介质/介质边界分解插值，进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间子迭代步数守恒电磁场数值。
9.一种电子设备，其特征在于，包括存储器和处理器，其中：
所述存储器，用于保存计算机程序；
所述处理器，用于执行所述计算机程序，以实现如权利要求1至7任一项所述的介质涂覆导体复合目标电磁散射隐式计算方法。
10.一种计算机可读存储介质，其特征在于，用于保存计算机程序，其中，所述计算机程序被处理器执行时实现如权利要求1至7任一项所述的介质涂覆导体复合目标电磁散射隐式计算方法。
CN202211411927.5A 2022-11-11 2022-11-11 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置 Pending CN115600435A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211411927.5A CN115600435A (zh) 2022-11-11 2022-11-11 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211411927.5A CN115600435A (zh) 2022-11-11 2022-11-11 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置

Publications (1)

Publication Number Publication Date
CN115600435A true CN115600435A (zh) 2023-01-13

Family

ID=84852399

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211411927.5A Pending CN115600435A (zh) 2022-11-11 2022-11-11 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置

Country Status (1)

Country Link
CN (1) CN115600435A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115995279A (zh) * 2023-03-22 2023-04-21 中国空气动力研究与发展中心计算空气动力研究所 一种材料力学特性评估方法、装置、设备及可读存储介质
CN116052820A (zh) * 2023-03-22 2023-05-02 中国空气动力研究与发展中心计算空气动力研究所 一种材料热性能评估方法、装置、设备及可读存储介质
CN116090262A (zh) * 2023-04-10 2023-05-09 中国空气动力研究与发展中心计算空气动力研究所 一种有限催化分级隐式数值模拟方法、装置、设备及介质

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115995279A (zh) * 2023-03-22 2023-04-21 中国空气动力研究与发展中心计算空气动力研究所 一种材料力学特性评估方法、装置、设备及可读存储介质
CN116052820A (zh) * 2023-03-22 2023-05-02 中国空气动力研究与发展中心计算空气动力研究所 一种材料热性能评估方法、装置、设备及可读存储介质
CN116090262A (zh) * 2023-04-10 2023-05-09 中国空气动力研究与发展中心计算空气动力研究所 一种有限催化分级隐式数值模拟方法、装置、设备及介质

Similar Documents

Publication Publication Date Title
CN112989680B (zh) 减少网格使用量的fvfd远场积分边界条件计算方法
CN115600435A (zh) 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置
Sumithra et al. Review on computational electromagnetics
Tamayo et al. Multilevel adaptive cross approximation (MLACA)
CN113158492B (zh) 一种时变电磁场的全隐式双时间步计算方法
CN113158527B (zh) 一种基于隐式fvfd计算频域电磁场的方法
Huang et al. Multibranch Rao–Wilton–Glisson basis functions for electromagnetic scattering problems
Jia et al. $\mathcal {H} $-Matrices Compressed Multiplicative Schwarz Preconditioner for Nonconformal FEM-BEM-DDM
Li et al. Fast periodic interpolation method for periodic unit cell problems
Negi Memory reduced half hierarchal matrix (H-matrix) for electrodynamic electric field integral equation
Chen et al. A memory-efficient implementation of perfectly matched layer with smoothly varying coefficients in discontinuous Galerkin time-domain method
Rong et al. Fast direct surface integral equation solution for electromagnetic scattering analysis with skeletonization factorization
CN112329204B (zh) 考虑载体平台耦合的重复性结构电磁特性特征模快速分析方法
CN115879276A (zh) 一种目标对象的电磁特性分析方法、装置及设备和介质
Zuo et al. A Massively Parallel Preconditioned FEM-BEM Method for Accurate Analysis of Complex Electromagnetic Field Problems
Zhou et al. Efficient EDM-PO method for the scattering from electrically large objects with the high-order impedance boundary condition
CN115510690A (zh) 基于awe技术的电大不确定外形金属目标的新型电磁特性计算方法
Kim et al. Adaptive integral method for higher order method of moments
Liu et al. A cylindrical MRTD algorithm with PML and quasi-PML
Lai et al. Analysis of nurbs surfaces modeled geometries with higher-order MoM based aim
Tong et al. Multilevel fast multipole algorithm for acoustic wave scattering by truncated ground with trenches
Liu et al. A combined immersed finite element and conservative semi-Lagrangian scheme for plasma-material interactions
Hollander et al. Adaptive multilevel nonuniform grid algorithm for the accelerated analysis of composite metallic–dielectric radomes
Kaufmann The meshless radial point interpolation method for electromagnetics
Tian et al. Accelerated Hybrid Method for Electromagnetic Scattering From Multiple Complex Targets Above a Rough Surface

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