CN113158527A - 一种基于隐式fvfd计算频域电磁场的方法 - Google Patents

一种基于隐式fvfd计算频域电磁场的方法 Download PDF

Info

Publication number
CN113158527A
CN113158527A CN202110526876.XA CN202110526876A CN113158527A CN 113158527 A CN113158527 A CN 113158527A CN 202110526876 A CN202110526876 A CN 202110526876A CN 113158527 A CN113158527 A CN 113158527A
Authority
CN
China
Prior art keywords
grid
electromagnetic
matrix
electromagnetic field
frequency domain
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.)
Granted
Application number
CN202110526876.XA
Other languages
English (en)
Other versions
CN113158527B (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 CN202110526876.XA priority Critical patent/CN113158527B/zh
Publication of CN113158527A publication Critical patent/CN113158527A/zh
Application granted granted Critical
Publication of CN113158527B publication Critical patent/CN113158527B/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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Abstract

本发明提供一种基于隐式FVFD计算频域电磁场的方法，包括根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；数值计算区域的四边形或六面体结构网格，在壁面和几何奇异处加密，所述网格逐渐远离散射壁面而逐渐稀疏；输出网格数据文件和设定和输出边界条件文件；输入目标计算电磁参数、数值计算控制参数；输入网格数据和边界条件信息文件，初始化计算空间电磁场；对麦克斯韦方程组频域电磁场进行迭代求解；输出电磁场的实部虚部空间分布，输出表面诱导电流和雷达散射截面空间分布数据。本发明可求解任意复杂外形、高频电大尺寸目标大规模电磁散射问题。

Description

一种基于隐式FVFD计算频域电磁场的方法
技术领域
本发明涉及电磁学的频域数值求解技术领域，尤其是涉及一种基于隐式FVFD（Finite Volume Frequency Domain, FVFD）计算频域电磁场的方法。
背景技术
复杂外形目标电磁散射、复杂电磁环境电磁干扰都需要计算电磁场空间分布，电磁场满足麦克斯韦方程组，随着计算机技术的发展直接求解该方程组成为可能。与欧拉方程相同的双曲型数学特征促进计算流体力学（Computational Fluid Dynamics, CFD）技术在电磁场计算中的应用，其中时域有限差分法(Finite Difference Time Domain, FDTD)和时域有限体积法(Finite Volume Time Domain, FVTD)最为著名。FVTD直接将麦克斯韦方程组守恒律的积分形式应用到离散的网格单元。由于许多有限差分法不能用于不连续函数的计算，而且在非定常计算中要保证一些重要量的守恒性，使得设计差分格式变得复杂，因此转而在物理空间使用守恒律的积分形式，允许计算不连续函数，从而使CFD专家中有限体积法变得非常流行。
时间域计算有利于模拟宽带脉冲电磁波信号辐射、散射，但如果入射波是单频简谐信号，则可以在频域中计算电磁场。传统频域方法主要有解析方法、高频近似方法和全波数值方法。解析方法仅可以求解特殊几何形状简单目标电磁散射，对于实际几何复杂目标无能为力。高频近似方法包括几何光学(Geometrical Optics, GO)、几何绕射理论(Geometrical Theory of Diffraction, GTD)、物理光学(Physical Optics, PO)、物理绕射理论(Physical Theory of Diffraction, PTD)、一致性几何绕射理论(Uniform Theoryof Diffraction, UTD)、一致性渐进理论(Uniform Asymptotic Theory, UAT)和等效边缘电流法(Method of Equivalent Current, MEC)等，基于高频场局部性原理仅考虑部件或细小单元在入射波下产生的散射场，不考虑部件或单元之间的相互耦合，高频方法在分析复杂结构目标的电磁散射时精度较差。全波数值方法直接求解 Maxwell 偏微分方程或电磁流积分方程，不作任何近似具有较高计算精度，在计算机资源允许情况下，全波数值方法可以求解任意频率的电磁问题。高精度全波电磁数值方法主要分为两类：一类是求解以电流为变量积分方程，包括矩量法（Method of Moment, MOM）及后续发展的多极子方法（FastMultipole Method, FMM）、多层快多极子方法(Multi-Level Fast Multipole Algorithm,MLFMA)；另一类求解以电磁场为变量的Maxwell微分或亥姆霍兹波动方程的FDTD方法和有限元方法（Finite Element Method, FEM）。
现有直接计算电磁场的微分类方法中，FDTD采用笛卡尔直角正交网格模拟壁面存在阶梯效应影响数值精度，并用电磁场分量时空交叉放置添加人工粘性给2阶中心差分格式。1992年Huh在AIAA-92-0453：“a compact high-order finite-volume time-domain/frequency-domain method for electromagnetic scattering”中提出一种Huh方法，Huh方法采用紧致差分结合滤波人工粘性构造通量，时间迭代采用点隐4步Runge-Kutta方法，过程繁琐复杂。1998年Bonnet在“Frequency-Domain Finite Volume Method forElectromagnetic Scattering ”中提出Bonnet方法，其Bonnet方法采用的解线性代数方程组的BICGSTAB(1)方法。一方面场计算微分类方法在模拟复杂细节结构、非均匀材料有独特优势和工程需求，另一方面现在频域微分类算法尚不完善。
综上所述，有必要发展一种能模拟任意外形、材料且可叠加各种加速算法的FVFD方法。
发明内容
本发明提供一种基于隐式FVFD计算频域电磁场的方法，其可求解任意复杂外形、高频电大尺寸目标大规模电磁散射问题，或者求解含复杂电子细节结构的多尺度MOM难以处理的电磁问题。
为实现上述目的，本发明提供如下技术方案：一种基于隐式FVFD计算频域电磁场的方法，包括以下步骤：
步骤1，根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；
步骤2，采用四边形或六面体结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件；
步骤3，输入目标计算电磁参数、数值计算控制参数；
步骤4，输入网格数据和边界条件信息文件，初始化计算空间电磁场；
步骤5，基于虚拟时间推进和空间通量残差对麦克斯韦方程组频域电磁场进行迭代求解；
步骤6，输出电磁场的实部虚部空间分布，输出表面诱导电流和雷达散射截面空间分布数据。
优选地，所述步骤5过程如下：对仿真模型做定常虚拟时间步循环，直至收敛结束；在每个虚拟时间迭代过程中，依次对各个网块格、各个网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间迭代步数守恒电磁场数值。
优选地，所述定常虚拟时间步循环，直至收敛结束过程为：
通过虚拟时间
Figure 821174DEST_PATH_IMAGE001
，将待求解的麦克斯韦方程组修正为：
Figure 37130DEST_PATH_IMAGE002
Figure 432339DEST_PATH_IMAGE003
其中，
Figure 184394DEST_PATH_IMAGE004
是入射简谐电磁波频率，
Figure 6857DEST_PATH_IMAGE005
为结构网格曲线坐标系方向1，
Figure 273890DEST_PATH_IMAGE006
为结构网格曲线坐标系方向2，
Figure 410473DEST_PATH_IMAGE007
为结构网格曲线坐标系方向3；
Figure 79352DEST_PATH_IMAGE008
分别对应为
Figure 10399DEST_PATH_IMAGE005
、
Figure 764728DEST_PATH_IMAGE009
、
Figure 937958DEST_PATH_IMAGE007
方向的电磁通量；
Figure 461344DEST_PATH_IMAGE010
是频域复数型电磁场守恒变量，
Figure 828871DEST_PATH_IMAGE001
是虚拟时间，
Figure 70497DEST_PATH_IMAGE011
是直角坐标系下电磁通量的
Figure 283303DEST_PATH_IMAGE012
分量，
Figure 864457DEST_PATH_IMAGE013
是频域复数型磁感应强度矢量，
Figure 199624DEST_PATH_IMAGE014
是频域复数型电位移矢量，
Figure 364764DEST_PATH_IMAGE015
是频域复数型电场强度矢量，
Figure 381261DEST_PATH_IMAGE016
是频域复数型磁场强度矢量，含下标
Figure 879239DEST_PATH_IMAGE012
标量分别是对应矢量的
Figure 650885DEST_PATH_IMAGE012
分量；上标
Figure 539207DEST_PATH_IMAGE017
是虚拟时间迭代步数，下标
Figure 421712DEST_PATH_IMAGE018
是网格单元编号，
Figure 711879DEST_PATH_IMAGE019
是隐式控制参数，取
Figure 654428DEST_PATH_IMAGE020
为全隐式，
Figure 528580DEST_PATH_IMAGE021
是第
Figure 886881DEST_PATH_IMAGE022
网格单元第
Figure 359450DEST_PATH_IMAGE023
虚拟时间迭代步时的电磁守恒变量，
Figure 410583DEST_PATH_IMAGE024
是第
Figure 273497DEST_PATH_IMAGE022
网格单元第
Figure 497805DEST_PATH_IMAGE025
虚拟时间迭代步时的电磁守恒变量，
Figure 995520DEST_PATH_IMAGE026
代表曲线坐标系频域复数型的电磁场守恒变量；
Figure 279871DEST_PATH_IMAGE027
是第
Figure 630081DEST_PATH_IMAGE022
网格单元第
Figure 658079DEST_PATH_IMAGE028
虚拟时间迭代步时的空间通量残差，
Figure 777345DEST_PATH_IMAGE029
是第
Figure 904701DEST_PATH_IMAGE022
网格单元第
Figure 804524DEST_PATH_IMAGE030
虚拟时间迭代步时的空间通量残差；
Figure 370635DEST_PATH_IMAGE031
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算；采用不同局部虚拟时间迭代步长定常虚拟时间步循环计算不同的网格单元，加快相应网格单元的单元电磁场收敛。
优选地，所述空间通量计算和隐式迭代解计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量
Figure 108521DEST_PATH_IMAGE032
Figure 469095DEST_PATH_IMAGE033
；
Figure 793897DEST_PATH_IMAGE034
；
Figure 101382DEST_PATH_IMAGE035
；
式中
Figure 991978DEST_PATH_IMAGE036
分别取曲线坐标系
Figure 789032DEST_PATH_IMAGE037
方向之一，相应的
Figure 335551DEST_PATH_IMAGE038
即为对应
Figure 509044DEST_PATH_IMAGE039
方向的电磁通量；
Figure 955943DEST_PATH_IMAGE040
代表曲线坐标系
Figure 861582DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 957714DEST_PATH_IMAGE041
代表曲线坐标系
Figure 669318DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 472189DEST_PATH_IMAGE042
为相似矩阵，
Figure 611047DEST_PATH_IMAGE043
分别为正负特征值构成的对角矩阵，
Figure 132158DEST_PATH_IMAGE044
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 381874DEST_PATH_IMAGE045
代表自变量为
Figure 803365DEST_PATH_IMAGE046
的
Figure 113124DEST_PATH_IMAGE047
相似矩阵；
Figure 121531DEST_PATH_IMAGE048
代表自变量为
Figure 112621DEST_PATH_IMAGE046
的
Figure 952401DEST_PATH_IMAGE049
对角矩阵；
Figure 370744DEST_PATH_IMAGE050
代表自变量为
Figure 663185DEST_PATH_IMAGE051
的
Figure 520283DEST_PATH_IMAGE052
相似矩阵；
Figure 650788DEST_PATH_IMAGE053
代表自变量为
Figure 240032DEST_PATH_IMAGE054
的
Figure 285348DEST_PATH_IMAGE047
相似矩阵；
Figure 618241DEST_PATH_IMAGE055
代表自变量为
Figure 167034DEST_PATH_IMAGE054
的
Figure 192758DEST_PATH_IMAGE056
对角矩阵；
Figure 459792DEST_PATH_IMAGE057
代表自变量为
Figure 94910DEST_PATH_IMAGE054
的
Figure 498210DEST_PATH_IMAGE058
相似矩阵；
Figure 757153DEST_PATH_IMAGE059
Figure 449165DEST_PATH_IMAGE060
其中
Figure 186177DEST_PATH_IMAGE061
是限制器，下标
Figure 647245DEST_PATH_IMAGE062
是网格单元编号，
Figure 77090DEST_PATH_IMAGE063
对应单元分界面，
Figure 990819DEST_PATH_IMAGE064
是3阶精度格式的控制参数，
Figure 531522DEST_PATH_IMAGE065
和
Figure 611211DEST_PATH_IMAGE066
分别是后差和前差算符；
Figure 946377DEST_PATH_IMAGE067
表示网格单元
Figure 675299DEST_PATH_IMAGE068
分界面处左状态电磁守恒变量，
Figure 691797DEST_PATH_IMAGE069
表示网格单元
Figure 127457DEST_PATH_IMAGE068
分界面处右状态电磁守恒变量；
Figure 633525DEST_PATH_IMAGE070
是第
Figure 787426DEST_PATH_IMAGE062
个网格单元电磁场守恒变量，
Figure 669931DEST_PATH_IMAGE071
是第
Figure 458633DEST_PATH_IMAGE062
+1个网格单元电磁场守恒变量；
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure 401181DEST_PATH_IMAGE072
其中
Figure 776799DEST_PATH_IMAGE073
是分裂后的系数矩阵，
Figure 197416DEST_PATH_IMAGE074
是上一个虚拟时间步计算的空间通量残差，
Figure 669986DEST_PATH_IMAGE075
是隐式子迭代电磁场差值；
将
Figure 455539DEST_PATH_IMAGE076
表示为LDU近似因式分解得：
Figure 380770DEST_PATH_IMAGE077
Figure 542761DEST_PATH_IMAGE078
Figure 306055DEST_PATH_IMAGE079
Figure 590406DEST_PATH_IMAGE080
其中，下标
Figure 2933DEST_PATH_IMAGE081
是网格单元编号，
Figure 703036DEST_PATH_IMAGE082
是雅可比系数矩阵最大特征值分裂参数，
Figure 884618DEST_PATH_IMAGE083
是雅克比系数矩阵最大特征值；
Figure 277553DEST_PATH_IMAGE084
是单位对角矩阵，
Figure 115059DEST_PATH_IMAGE085
是对角矩阵，
Figure 681170DEST_PATH_IMAGE086
为上三角矩阵，
Figure 717259DEST_PATH_IMAGE087
为下三角矩阵，
Figure 779631DEST_PATH_IMAGE088
是对应上三角矩阵的电磁守恒变量差值，
Figure 104433DEST_PATH_IMAGE089
是对应下三角矩阵的电磁守恒变量差值；
Figure 208655DEST_PATH_IMAGE090
指
Figure 302513DEST_PATH_IMAGE091
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 99568DEST_PATH_IMAGE092
指
Figure 646087DEST_PATH_IMAGE093
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 819579DEST_PATH_IMAGE094
指
Figure 289916DEST_PATH_IMAGE095
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 195555DEST_PATH_IMAGE096
指
Figure 291687DEST_PATH_IMAGE097
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 675395DEST_PATH_IMAGE098
指
Figure 806162DEST_PATH_IMAGE099
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 882703DEST_PATH_IMAGE100
指
Figure 466131DEST_PATH_IMAGE101
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 152065DEST_PATH_IMAGE102
是指相邻
Figure 137338DEST_PATH_IMAGE103
网格单元分裂后的系数矩阵；
Figure 447097DEST_PATH_IMAGE104
是指相邻
Figure 455504DEST_PATH_IMAGE105
网格单元分裂后的系数矩阵；
Figure 446594DEST_PATH_IMAGE106
是指相邻
Figure 286374DEST_PATH_IMAGE107
网格单元分裂后的系数矩阵；
Figure 704717DEST_PATH_IMAGE108
是指相邻
Figure 997158DEST_PATH_IMAGE109
网格单元分裂后的系数矩阵；
Figure 290474DEST_PATH_IMAGE110
是指相邻
Figure 984760DEST_PATH_IMAGE111
网格单元分裂后的系数矩阵；
Figure 574005DEST_PATH_IMAGE112
是指相邻
Figure 619321DEST_PATH_IMAGE113
网格单元分裂后的系数矩阵；
最后，得到前后向迭代计算电磁场差值
Figure 952214DEST_PATH_IMAGE114
：
Figure 501007DEST_PATH_IMAGE115
Figure 526731DEST_PATH_IMAGE116
Figure 793765DEST_PATH_IMAGE117
其中，
Figure 163304DEST_PATH_IMAGE118
是对角矩阵
Figure 832183DEST_PATH_IMAGE119
的逆矩阵，
Figure 28809DEST_PATH_IMAGE120
、
Figure 517559DEST_PATH_IMAGE121
分别是根据
Figure 457833DEST_PATH_IMAGE122
、
Figure 981218DEST_PATH_IMAGE123
计算得出的上三角矩阵、下三角矩阵；
前向循环：
Figure 83166DEST_PATH_IMAGE124
后向循环：
Figure 761010DEST_PATH_IMAGE125
其中，
Figure 301713DEST_PATH_IMAGE126
是电磁守恒变量差值的中间过渡变量。
优选地，所述步骤2中: 网格密度保证每波长13-20个网格点，壁面密度>300点/波长，几何奇异处加密到50-100个网格点/波长，二维网格在垂直该二维网格所在平面按右手法则推进一层，作为三维问题特例统一计算；网格数据文件包括结构网格块数目及每块三个曲线坐标系下维度。
优选地，所述步骤3中：如有等离子体外部流畅情况，则还应输入对应流场参数。
优选地，所述定常虚拟时间步循环为隐式，其CFL数不受显式稳定性要求约束。
与现有技术相比，本发明的有益效果是：
1、本发明可求解任意复杂外形、高频电大尺寸目标大规模电磁散射问题，或者求解含复杂电子细节结构的多尺度MOM难以处理的电磁问题。支持曲线坐标系结构网格和多区域分解以及并行算法。对于简谐波单频入射情况（此时时间信号为连续周期信号），为提高计算效率和精度并降低网格量，在Maxwell方程组中对时间直接求导转化到频域，使变量维度从时空4维降到仅余空间的3维，从非定常算法降为定常算法，变量为复数型优美且更有效率，减少时域计算目标电磁特征所需的傅里叶变换环节。
2、本发明的隐式频域有限体积(FVFD)方法不同于FDTD，FDTD采用笛卡尔直角正交网格模拟壁面存在阶梯效应影响数值精度，并用电磁场分量时空交叉放置添加人工粘性给2阶中心差分格式，本发明FVFD采用贴体曲线坐标系网格能更好拟合物面和在几何奇异处加密网格，并且电磁场量在网格空间同置于网格单元中心，采用迎风格式来保持人工粘性，更有利保持精度和算法设计。
3、本发明不同于有限元方法，FEM和FVFD都可采用任意形状网格单元模拟离散计算空间，FEM采用基函数模拟节点或棱边矢量电场、磁场矢量之一，并用变分方法或残值加权构造矩阵形式方程组，得到带状离散全空间矩阵并求解该线性代数方程组，本发明的FVFD直接求解频域Maxwell方程组控制的所有6个电磁场分量，且不需要每步计算稀疏矩阵逆运算。
4、本发明与1992年Huh在AIAA-92-0453：“a compact high-order finite-volumetime-domain/frequency-domain method for electromagnetic scattering”中的FVFD不同，Huh方法采用紧致差分结合滤波人工粘性构造通量，时间迭代采用点隐4步Runge-Kutta方法，过程繁琐复杂，本FVFD采用迎风插值Steger-Warming分裂构造通量及隐式迭代计算电磁场量。
5、本发明与 1998年Bonnet在“Frequency-Domain Finite Volume Method forElectromagnetic Scattering”中的FVFD不同，其一，本发明FVFD通量采用Steger-Warming分裂非简单相位几何关系，其二，本发明FVFD添加虚拟时间步，采用隐式前后向分裂矩阵迭代计算，Bonnet则采用的解线性代数方程组的BICGSTAB(1)方法。
6、本发明提高推进效率的局部时间步定常迭代推进方法，结合隐式空间通量残差获得稳定、高效的计算流程，摆脱传统显式有限体积全局对迭代时间步的限制所带来的大计算量弊端。
附图说明
图1为运用本发明计算频域电磁场的流程示意图；
图2为采用本发明对金属球FVFD计算定常迭代收敛过程(ka=10)示意图；
图3为采用本发明对金属球FVFD计算表面诱导电流分布(ka=10)示意图；
图4为采用本发明对金属球FVFD计算的双站雷达散射截面与解析解比较(ka=10)示意图；
图5为采用本发明对金属橄榄球电磁散射FVFD计算的复数型电磁场量幅值定常迭代收敛过程示意图；
图6为采用本发明对金属橄榄球FVFD计算的表面诱导电流分布(f=1.18GHz)示意图；
图7为采用本发明对金属橄榄球FVFD计算的双站雷达散射截面分布与MOM方法的比较(f=1.18GHz)示意图。
具体实施方式
参照图1至图7对本发明一种基于隐式FVFD计算频域电磁场的方法的实施例进一步说明。
参照附图1，整个频域有限体积方法计算电磁场软件按结构可以分为：预处理、电磁场计算和后处理三个部分。预处理主要包括网格数据输入、计算参数数据输入、控制参数输入三个模块，主要用来读入网格数据、计算参数数据输入、控制参数文件，并在此基础上，进行预先处理，为电磁场计算提供计算支撑；电磁场计算包括：空间电磁场MUSCL格式插值、单元分界面通量计算、时间推进、收敛判断模块组成；后处理主要用于输出电磁场的空间实部虚部分布、目标表面诱导电流密度、雷达散射截面输出。
下面结合要数值模拟的频域Maxwell方程组两个旋度方程（时间因子
Figure 617288DEST_PATH_IMAGE127
），法拉第(Faraday)电磁感应定律：
Figure 218033DEST_PATH_IMAGE128
；安培（Ampere）定理：
Figure 681376DEST_PATH_IMAGE129
，介绍隐式频域有限体积法数值计算过程。其中
Figure 963453DEST_PATH_IMAGE130
是对应复数变量的虚数符号，
Figure 461430DEST_PATH_IMAGE131
是简谐电磁波频率，
Figure 905181DEST_PATH_IMAGE132
是复数型磁感应强度矢量，
Figure 121399DEST_PATH_IMAGE133
是复数型电位移矢量，
Figure 440122DEST_PATH_IMAGE134
是复数型电场强度矢量，
Figure 792606DEST_PATH_IMAGE135
是复数型磁场强度矢量，
Figure 407258DEST_PATH_IMAGE136
是外加强迫电流。
频域麦克斯韦方程组的两个旋度方程无源条件下的直角坐标系守恒形式为：
Figure 845193DEST_PATH_IMAGE137
Figure 469072DEST_PATH_IMAGE138
其中，
Figure 613746DEST_PATH_IMAGE004
是入射简谐电磁波频率，
Figure 727195DEST_PATH_IMAGE139
是频域复数型电磁场守恒变量，
Figure 88644DEST_PATH_IMAGE140
是添加的虚拟时间，
Figure 312952DEST_PATH_IMAGE141
是直角坐标系下电磁通量的
Figure 640028DEST_PATH_IMAGE142
分量，
Figure 862062DEST_PATH_IMAGE143
是频域复数型磁感应强度矢量，
Figure 946693DEST_PATH_IMAGE133
是频域复数型电位移矢量，
Figure 974692DEST_PATH_IMAGE144
是频域复数型电场强度矢量，
Figure 93957DEST_PATH_IMAGE145
是频域复数型磁场强度矢量，含下标
Figure 549209DEST_PATH_IMAGE142
标量分别是对应矢量的
Figure 885251DEST_PATH_IMAGE142
分量。明显可见当
Figure 451361DEST_PATH_IMAGE146
收敛时，该方程组等同于原始方程组。
对于复杂外形物体，采用的是计算空间贴体多块结构网格，因此均存在坐标变换：
Figure 425133DEST_PATH_IMAGE147
Figure 51287DEST_PATH_IMAGE148
其中
Figure 110510DEST_PATH_IMAGE149
代表曲线坐标系
Figure 480311DEST_PATH_IMAGE150
三个方向，得到所要数值模拟的曲线坐标系下的麦克斯韦方程组守恒形：
Figure 636486DEST_PATH_IMAGE151
；
Figure 371224DEST_PATH_IMAGE152
；
Figure 416278DEST_PATH_IMAGE153
；
Figure 324191DEST_PATH_IMAGE154
式中，
Figure 334872DEST_PATH_IMAGE155
是坐标变换的雅可比矩阵， 对应^上标变量代表曲线坐标系下的值，由坐标变换获取。
Figure 240512DEST_PATH_IMAGE156
为曲线坐标系下的电磁守恒变量；
Figure 274327DEST_PATH_IMAGE005
为结构网格曲线坐标系方向1，
Figure 985931DEST_PATH_IMAGE006
为结构网格曲线坐标系方向2，
Figure 54381DEST_PATH_IMAGE007
为结构网格曲线坐标系方向3；
Figure 193238DEST_PATH_IMAGE157
为曲线坐标系下的电磁通量。
Figure 884988DEST_PATH_IMAGE158
分别取
Figure 400283DEST_PATH_IMAGE159
三个曲线坐标系下的方向之一。
为了摆脱传统显式有限体积全局对迭代时间步的限制所带来的大计算量弊端，为此本发明基于FVFD计算频域电磁场的方法利用局部时间步定常迭代推进，结合隐式空间通量残差获得稳定、高效的计算流程，其包括以下步骤：
步骤1：根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模。
步骤2：采用四边形（2维）或六面体（3维）结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格随逐渐远离散射壁面逐渐稀疏。数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件。网格密度保证每波长13-20个网格点，壁面密度>300点每波长，几何奇异处加密到50-100个网格点每波长，二维网格在垂直该二维网格所在平面按右手法则推进一层，作为三维问题特例统一计算。网格数据文件包括结构网格块数目及每块三个曲线坐标系下维度。
步骤3：预处理部分，输入目标计算电磁参数、数值计算控制参数、有等离子体外部流场情况输入对应流场参数文件。虚拟时间迭代因是隐式其CFL数不受显式稳定性要求约束。
步骤4：输入网格数据和边界条件信息文件，初始化计算空间电磁场。
步骤5：基于虚拟时间推进和空间通量残差对麦克斯韦方程组频域电磁场进行迭代求解。
步骤5-1：定常虚拟时间步循环，到计算收敛结束。
Figure 323240DEST_PATH_IMAGE160
其中，
Figure 632999DEST_PATH_IMAGE004
是入射简谐电磁波频率，
Figure 375827DEST_PATH_IMAGE005
为结构网格曲线坐标系方向1，
Figure 429233DEST_PATH_IMAGE006
为结构网格曲线坐标系方向2，
Figure 705232DEST_PATH_IMAGE007
为结构网格曲线坐标系方向3；
Figure 123575DEST_PATH_IMAGE008
分别对应为
Figure 681595DEST_PATH_IMAGE005
、
Figure 210796DEST_PATH_IMAGE009
、
Figure 170662DEST_PATH_IMAGE007
方向的电磁通量；
Figure 759906DEST_PATH_IMAGE010
是频域复数型电磁场守恒变量，
Figure 539644DEST_PATH_IMAGE001
是虚拟时间，
Figure 371071DEST_PATH_IMAGE011
是直角坐标系下电磁通量的
Figure 185443DEST_PATH_IMAGE012
分量，
Figure 7906DEST_PATH_IMAGE013
是频域复数型磁感应强度矢量，
Figure 212622DEST_PATH_IMAGE014
是频域复数型电位移矢量，
Figure 411522DEST_PATH_IMAGE015
是频域复数型电场强度矢量，
Figure 18084DEST_PATH_IMAGE016
是频域复数型磁场强度矢量，含下标
Figure 949131DEST_PATH_IMAGE012
标量分别是对应矢量的
Figure 703461DEST_PATH_IMAGE012
分量；上标
Figure 706052DEST_PATH_IMAGE017
是虚拟时间迭代步数，下标
Figure 665655DEST_PATH_IMAGE018
是网格单元编号，
Figure 767603DEST_PATH_IMAGE019
是隐式控制参数，取
Figure 9229DEST_PATH_IMAGE020
为全隐式，
Figure 487615DEST_PATH_IMAGE021
是第
Figure 865506DEST_PATH_IMAGE022
网格单元第
Figure 466252DEST_PATH_IMAGE023
虚拟时间迭代步时的电磁守恒变量，
Figure 867277DEST_PATH_IMAGE024
是第
Figure 647889DEST_PATH_IMAGE022
网格单元第
Figure 880288DEST_PATH_IMAGE025
虚拟时间迭代步时的电磁守恒变量，
Figure 589618DEST_PATH_IMAGE026
代表曲线坐标系频域复数型的电磁场守恒变量；
Figure 477939DEST_PATH_IMAGE027
是第
Figure 298128DEST_PATH_IMAGE022
网格单元第
Figure 650612DEST_PATH_IMAGE028
虚拟时间迭代步时的空间通量残差，
Figure 29378DEST_PATH_IMAGE029
是第
Figure 467313DEST_PATH_IMAGE022
网格单元第
Figure 91192DEST_PATH_IMAGE030
虚拟时间迭代步时的空间通量残差；
Figure 298182DEST_PATH_IMAGE031
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算。显著区别于显式方法，定常计算不同的网格单元采用不同局部虚拟时间迭代步长从而加快该单元电磁场收敛。
步骤5-2：在每个虚拟时间迭代过程中，逐网格块、逐网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间迭代步数守恒电磁场数值。
有限体积法的空间精度体现在能否精确模拟守恒变量Q在网格单元分界面处的状态变量，以得到相应精确的分界面通量
Figure 349315DEST_PATH_IMAGE161
，采用Steger-Warming分裂计算网格单元分界面通量。
Figure 212229DEST_PATH_IMAGE162
；
Figure 436537DEST_PATH_IMAGE034
；
Figure 199831DEST_PATH_IMAGE035
；
式中
Figure 218603DEST_PATH_IMAGE036
分别取曲线坐标系
Figure 631130DEST_PATH_IMAGE037
方向之一，相应的
Figure 596812DEST_PATH_IMAGE038
即为对应
Figure 716077DEST_PATH_IMAGE039
方向的电磁通量；
Figure 171329DEST_PATH_IMAGE040
代表曲线坐标系
Figure 743256DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 309367DEST_PATH_IMAGE041
代表曲线坐标系
Figure 70691DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 696844DEST_PATH_IMAGE042
为相似矩阵，
Figure 756067DEST_PATH_IMAGE043
分别为正负特征值构成的对角矩阵，
Figure 125869DEST_PATH_IMAGE044
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 219727DEST_PATH_IMAGE045
代表自变量为
Figure 16781DEST_PATH_IMAGE046
的
Figure 563300DEST_PATH_IMAGE047
相似矩阵；
Figure 471213DEST_PATH_IMAGE048
代表自变量为
Figure 918113DEST_PATH_IMAGE046
的
Figure 886069DEST_PATH_IMAGE049
对角矩阵；
Figure 919884DEST_PATH_IMAGE050
代表自变量为
Figure 569171DEST_PATH_IMAGE051
的
Figure 434359DEST_PATH_IMAGE052
相似矩阵；
Figure 510899DEST_PATH_IMAGE053
代表自变量为
Figure 94328DEST_PATH_IMAGE054
的
Figure 45841DEST_PATH_IMAGE047
相似矩阵；
Figure 765535DEST_PATH_IMAGE055
代表自变量为
Figure 12977DEST_PATH_IMAGE054
的
Figure 83701DEST_PATH_IMAGE056
对角矩阵；
Figure 74791DEST_PATH_IMAGE057
代表自变量为
Figure 914571DEST_PATH_IMAGE054
的
Figure 332914DEST_PATH_IMAGE058
相似矩阵。
Figure 890934DEST_PATH_IMAGE163
Figure 918671DEST_PATH_IMAGE164
其中
Figure 612957DEST_PATH_IMAGE061
是限制器，下标
Figure 530098DEST_PATH_IMAGE062
是网格单元编号，
Figure 247518DEST_PATH_IMAGE063
对应单元分界面，
Figure 642727DEST_PATH_IMAGE064
是3阶精度格式的控制参数，
Figure 394782DEST_PATH_IMAGE065
和
Figure 217245DEST_PATH_IMAGE066
分别是后差和前差算符；
Figure 421961DEST_PATH_IMAGE067
表示网格单元
Figure 620862DEST_PATH_IMAGE068
分界面处左状态电磁守恒变量，
Figure 725959DEST_PATH_IMAGE069
表示网格单元
Figure 719322DEST_PATH_IMAGE068
分界面处右状态电磁守恒变量；
Figure 411335DEST_PATH_IMAGE070
是第
Figure 148347DEST_PATH_IMAGE062
个网格单元电磁场守恒变量，
Figure 609415DEST_PATH_IMAGE071
是第
Figure 39259DEST_PATH_IMAGE062
+1个网格单元电磁场守恒变量。
空间通量隐式迭代和雅可比系数矩阵的分裂前后向迭代求解,用两次循环替代了稀疏矩阵求逆，工程上简单易用。经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到计算：
Figure 218568DEST_PATH_IMAGE165
其中
Figure 493691DEST_PATH_IMAGE166
是分裂后的系数矩阵，
Figure 573381DEST_PATH_IMAGE074
是上一个迭代时间步计算的空间通量残差，
Figure 908547DEST_PATH_IMAGE167
是隐式子迭代电磁场差值。该方程表示为LDU近似因式分解
Figure 637469DEST_PATH_IMAGE168
；
Figure 653966DEST_PATH_IMAGE169
；
Figure 89627DEST_PATH_IMAGE079
；
Figure 861274DEST_PATH_IMAGE080
；
其中，下标
Figure 749595DEST_PATH_IMAGE081
是网格单元编号，
Figure 632101DEST_PATH_IMAGE082
是雅可比系数矩阵最大特征值分裂参数，
Figure 984585DEST_PATH_IMAGE083
是雅克比系数矩阵最大特征值；
Figure 363351DEST_PATH_IMAGE084
是单位对角矩阵，
Figure 738969DEST_PATH_IMAGE085
是对角矩阵，
Figure 159586DEST_PATH_IMAGE086
为上三角矩阵，
Figure 632155DEST_PATH_IMAGE087
为下三角矩阵，
Figure 683288DEST_PATH_IMAGE088
是对应上三角矩阵的电磁守恒变量差值，
Figure 608519DEST_PATH_IMAGE089
是对应下三角矩阵的电磁守恒变量差值；
Figure 770510DEST_PATH_IMAGE090
指
Figure 832007DEST_PATH_IMAGE091
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 552576DEST_PATH_IMAGE092
指
Figure 965102DEST_PATH_IMAGE093
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 930784DEST_PATH_IMAGE094
指
Figure 112367DEST_PATH_IMAGE095
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 239723DEST_PATH_IMAGE096
指
Figure 139546DEST_PATH_IMAGE097
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 643340DEST_PATH_IMAGE098
指
Figure 945008DEST_PATH_IMAGE099
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 305582DEST_PATH_IMAGE100
指
Figure 128919DEST_PATH_IMAGE101
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 498721DEST_PATH_IMAGE102
是指相邻
Figure 327000DEST_PATH_IMAGE103
网格单元分裂后的系数矩阵；
Figure 61737DEST_PATH_IMAGE104
是指相邻
Figure 670573DEST_PATH_IMAGE105
网格单元分裂后的系数矩阵；
Figure 781749DEST_PATH_IMAGE106
是指相邻
Figure 792430DEST_PATH_IMAGE107
网格单元分裂后的系数矩阵；
Figure 196604DEST_PATH_IMAGE108
是指相邻
Figure 292736DEST_PATH_IMAGE109
网格单元分裂后的系数矩阵；
Figure 942024DEST_PATH_IMAGE110
是指相邻
Figure 807211DEST_PATH_IMAGE111
网格单元分裂后的系数矩阵；
Figure 946069DEST_PATH_IMAGE112
是指相邻
Figure 467180DEST_PATH_IMAGE113
网格单元分裂后的系数矩阵；
最后得到前后向迭代计算电磁场差值
Figure 654579DEST_PATH_IMAGE170
：
Figure 639852DEST_PATH_IMAGE171
Figure 385829DEST_PATH_IMAGE172
Figure 456553DEST_PATH_IMAGE173
；
其中，
Figure 447643DEST_PATH_IMAGE118
是对角矩阵
Figure 287423DEST_PATH_IMAGE119
的逆矩阵，
Figure 705766DEST_PATH_IMAGE120
、
Figure 998207DEST_PATH_IMAGE121
分别是根据
Figure 855305DEST_PATH_IMAGE122
、
Figure 487274DEST_PATH_IMAGE123
计算得出的上三角矩阵、下三角矩阵；
前向循环：
Figure 138836DEST_PATH_IMAGE174
后向循环
Figure 620370DEST_PATH_IMAGE175
其中，
Figure 953263DEST_PATH_IMAGE126
是电磁守恒变量差值的中间过渡变量。
以上是隐式FVFD计算Maxwell方程组所控制频域电磁场的迭代过程。
步骤6：收敛判断，后处理过程，输出电磁场的实部虚部空间分布，输出表面诱导电流和雷达散射截面空间分布数据等。
如图2至图4所示，是金属球(ka=10)隐式FVFD计算结果，计算参数为：2个网格数据块，维度都为46x97x25,CFL=1000,
Figure 502056DEST_PATH_IMAGE176
为全隐式，其中物面网格选取每波长20个网格点，远场边界在3波长外，辐向网格壁面加密，图2是金属球FVFD计算定常迭代收敛过程(ka=10)，收敛后结束计算标准全空间相邻迭代步电磁场最大幅度插值<0.01，图3是金属球FVFD计算的表面诱导电流分布(ka=10),图4是金属球FVFD计算的双站雷达散射截面与解析解比较(ka=10)，两者吻合良好。而且可见任意CFL数的计算说明隐式FVFD的无条件稳定性。
如图5至图7所示，以金属橄榄球电磁散射为例，计算参数为：1个网格数据块，维度都为40x46x65，CFL=5,
Figure 527780DEST_PATH_IMAGE177
为全隐式，远场边界在3波长外，辐向网格壁面加密，图5是金属橄榄球电磁散射FVFD计算的复数型电磁场量幅值定常迭代收敛过程，明显不同于时域方法简谐波条件收敛后的周期性振荡波形，图6是金属橄榄球FVFD计算的表面诱导电流分布(f=1.18GHz)，图7是金属橄榄球FVFD计算的双站雷达散射截面分布与MOM方法的比较(f=1.18GHz)，可见两者即使在雷达散射截面-40dB=0.0001m2 很小量级依然吻合很好，间接说明该运用本发明基于隐式FVFD计算频域电磁场方法计算出的高数值计算精度非常好。
这几个数值算例表明，该发明对应的隐式频域有限体积FVFD计算方法能在放松传统显式算法对迭代虚拟时间步限制的同时，获得无条件稳定的虚拟时间迭代推进，能保证数值精度并提升计算效率。
以上仅为本发明的较佳实施例而已，并不用以限制本发明，应当指出的是，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于隐式FVFD计算频域电磁场的方法，其特征在于，其包括以下步骤：
步骤1，根据目标所仿真电磁问题的物理背景，结合边界条件信息进行仿真建模；
步骤2，采用四边形或六面体结构对仿真模型进行网格剖分，网格在壁面和几何奇异处加密，网格逐渐远离散射壁面而逐渐稀疏；数值计算对应区域的网格，输出网格数据文件，设定和输出边界条件文件；
步骤3，输入目标计算电磁参数、数值计算控制参数；
步骤4，输入网格数据和边界条件信息文件，初始化计算空间电磁场；
步骤5，基于虚拟时间推进和空间通量残差对麦克斯韦方程组频域电磁场进行迭代求解；
步骤6，输出电磁场的实部虚部空间分布，输出表面诱导电流和雷达散射截面空间分布数据。
2.根据权利要求1所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述步骤5过程如下：对仿真模型做定常虚拟时间步循环，直至收敛结束；在每个虚拟时间迭代过程中，依次对各个网块格、各个网格单元进行空间通量计算和隐式迭代解计算，更新下一级虚拟时间迭代步数守恒电磁场数值。
3.根据权利要求2所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述定常虚拟时间步循环，直至收敛结束过程为：
通过虚拟时间
Figure 505002DEST_PATH_IMAGE001
，将待求解的麦克斯韦方程组修正为：
Figure 837895DEST_PATH_IMAGE002
Figure 386688DEST_PATH_IMAGE003
其中，
Figure 910948DEST_PATH_IMAGE004
是入射简谐电磁波频率，
Figure 115664DEST_PATH_IMAGE005
为结构网格曲线坐标系方向1，
Figure 48985DEST_PATH_IMAGE006
为结构网格曲线坐标系方向2，
Figure 655547DEST_PATH_IMAGE007
为结构网格曲线坐标系方向3；
Figure 914490DEST_PATH_IMAGE008
分别对应为
Figure 340923DEST_PATH_IMAGE005
、
Figure 343514DEST_PATH_IMAGE009
、
Figure 303118DEST_PATH_IMAGE007
方向的电磁通量；
Figure 467383DEST_PATH_IMAGE010
是频域复数型电磁场守恒变量，
Figure 646691DEST_PATH_IMAGE001
是虚拟时间，
Figure 187394DEST_PATH_IMAGE011
是直角坐标系下电磁通量的
Figure 502969DEST_PATH_IMAGE012
分量，
Figure 103715DEST_PATH_IMAGE013
是频域复数型磁感应强度矢量，
Figure 504740DEST_PATH_IMAGE014
是频域复数型电位移矢量，
Figure 849134DEST_PATH_IMAGE015
是频域复数型电场强度矢量，
Figure 783329DEST_PATH_IMAGE016
是频域复数型磁场强度矢量，含下标
Figure 289397DEST_PATH_IMAGE012
标量分别是对应矢量的
Figure 443298DEST_PATH_IMAGE012
分量；上标
Figure 325803DEST_PATH_IMAGE017
是虚拟时间迭代步数，下标
Figure 678287DEST_PATH_IMAGE018
是网格单元编号，
Figure 292939DEST_PATH_IMAGE019
是隐式控制参数，取
Figure 668557DEST_PATH_IMAGE020
为全隐式，
Figure 354753DEST_PATH_IMAGE021
是第
Figure 997962DEST_PATH_IMAGE022
网格单元第
Figure 111411DEST_PATH_IMAGE023
虚拟时间迭代步时的电磁守恒变量，
Figure 974325DEST_PATH_IMAGE024
是第
Figure 136316DEST_PATH_IMAGE022
网格单元第
Figure 463392DEST_PATH_IMAGE025
虚拟时间迭代步时的电磁守恒变量，
Figure 685426DEST_PATH_IMAGE026
代表曲线坐标系频域复数型的电磁场守恒变量；
Figure 832374DEST_PATH_IMAGE027
是第
Figure 296591DEST_PATH_IMAGE022
网格单元第
Figure 478174DEST_PATH_IMAGE028
虚拟时间迭代步时的空间通量残差，
Figure 871109DEST_PATH_IMAGE029
是第
Figure 770932DEST_PATH_IMAGE022
网格单元第
Figure 337042DEST_PATH_IMAGE030
虚拟时间迭代步时的空间通量残差；
Figure 310815DEST_PATH_IMAGE031
是由稳定性控制的虚拟时间步长，由CFL数和局部网格单元几何尺度和特征值计算；采用不同局部虚拟时间迭代步长定常虚拟时间步循环计算不同的网格单元，加快相应网格单元的单元电磁场收敛。
4.根据权利要求3所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述空间通量计算和隐式迭代解计算过程如下：
采用Steger-Warming分裂计算网格单元分界面通量
Figure 936968DEST_PATH_IMAGE032
Figure 996191DEST_PATH_IMAGE033
；
Figure DEST_PATH_IMAGE034
；
Figure 802211DEST_PATH_IMAGE035
；
式中
Figure DEST_PATH_IMAGE036
分别取曲线坐标系
Figure 958385DEST_PATH_IMAGE037
方向之一，相应的
Figure 693123DEST_PATH_IMAGE038
即为对应
Figure 301959DEST_PATH_IMAGE039
方向的电磁通量；
Figure 147555DEST_PATH_IMAGE040
代表曲线坐标系
Figure 158237DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，正特征值进行分裂后得到的电磁通量；
Figure 63876DEST_PATH_IMAGE041
代表曲线坐标系
Figure 160008DEST_PATH_IMAGE037
对应方向Steger-Warming分裂中，负特征值进行分裂后得到的电磁通量；
Figure 331268DEST_PATH_IMAGE042
为相似矩阵，
Figure 462035DEST_PATH_IMAGE043
分别为正负特征值构成的对角矩阵，
Figure 538575DEST_PATH_IMAGE044
分别代表分界面处左状态变量、右状态变量，采用MUSCL格式而达到最高三阶精度；
Figure 856424DEST_PATH_IMAGE045
代表自变量为
Figure 309402DEST_PATH_IMAGE046
的
Figure DEST_PATH_IMAGE047
相似矩阵；
Figure 232359DEST_PATH_IMAGE048
代表自变量为
Figure 542117DEST_PATH_IMAGE046
的
Figure DEST_PATH_IMAGE049
对角矩阵；
Figure 783480DEST_PATH_IMAGE050
代表自变量为
Figure DEST_PATH_IMAGE051
的
Figure 774570DEST_PATH_IMAGE052
相似矩阵；
Figure DEST_PATH_IMAGE053
代表自变量为
Figure 552033DEST_PATH_IMAGE054
的
Figure 32693DEST_PATH_IMAGE047
相似矩阵；
Figure DEST_PATH_IMAGE055
代表自变量为
Figure 528397DEST_PATH_IMAGE054
的
Figure 119915DEST_PATH_IMAGE056
对角矩阵；
Figure DEST_PATH_IMAGE057
代表自变量为
Figure 453682DEST_PATH_IMAGE054
的
Figure 105243DEST_PATH_IMAGE058
相似矩阵；
Figure DEST_PATH_IMAGE059
Figure 822664DEST_PATH_IMAGE060
其中
Figure DEST_PATH_IMAGE061
是限制器，下标
Figure 155556DEST_PATH_IMAGE062
是网格单元编号，
Figure DEST_PATH_IMAGE063
对应单元分界面，
Figure 406146DEST_PATH_IMAGE064
是3阶精度格式的控制参数，
Figure DEST_PATH_IMAGE065
和
Figure 166292DEST_PATH_IMAGE066
分别是后差和前差算符；
Figure DEST_PATH_IMAGE067
表示网格单元
Figure 371008DEST_PATH_IMAGE068
分界面处左状态电磁守恒变量，
Figure DEST_PATH_IMAGE069
表示网格单元
Figure 507592DEST_PATH_IMAGE068
分界面处右状态电磁守恒变量；
Figure 612689DEST_PATH_IMAGE070
是第
Figure 543736DEST_PATH_IMAGE062
个网格单元电磁场守恒变量，
Figure DEST_PATH_IMAGE071
是第
Figure 298065DEST_PATH_IMAGE062
+1个网格单元电磁场守恒变量；
采用空间通量隐式迭代和雅克比系数矩阵的分裂前后向迭代求解，经过通量偏导守恒变量产生的雅克比系数Steger-Warming分裂得到，
Figure 238339DEST_PATH_IMAGE072
其中
Figure DEST_PATH_IMAGE073
是分裂后的系数矩阵，
Figure 699408DEST_PATH_IMAGE074
是上一个虚拟时间步计算的空间通量残差，
Figure 863673DEST_PATH_IMAGE075
是隐式子迭代电磁场差值；
将
Figure 541516DEST_PATH_IMAGE076
表示为LDU近似因式分解得：
Figure 82219DEST_PATH_IMAGE077
Figure 397794DEST_PATH_IMAGE078
Figure 998540DEST_PATH_IMAGE079
Figure 399565DEST_PATH_IMAGE080
其中，下标
Figure 743959DEST_PATH_IMAGE081
是网格单元编号，
Figure 914040DEST_PATH_IMAGE082
是雅可比系数矩阵最大特征值分裂参数，
Figure 685687DEST_PATH_IMAGE083
是雅克比系数矩阵最大特征值；
Figure 72544DEST_PATH_IMAGE084
是单位对角矩阵，
Figure 955049DEST_PATH_IMAGE085
是对角矩阵，
Figure 245216DEST_PATH_IMAGE086
为上三角矩阵，
Figure 125447DEST_PATH_IMAGE087
为下三角矩阵，
Figure 563382DEST_PATH_IMAGE088
是对应上三角矩阵的电磁守恒变量差值，
Figure 249578DEST_PATH_IMAGE089
是对应下三角矩阵的电磁守恒变量差值；
Figure 394252DEST_PATH_IMAGE090
指
Figure 943920DEST_PATH_IMAGE091
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 869150DEST_PATH_IMAGE092
指
Figure 93458DEST_PATH_IMAGE093
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 358217DEST_PATH_IMAGE094
指
Figure 314672DEST_PATH_IMAGE095
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 727199DEST_PATH_IMAGE096
指
Figure 692881DEST_PATH_IMAGE097
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 874463DEST_PATH_IMAGE098
指
Figure 329716DEST_PATH_IMAGE099
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 400178DEST_PATH_IMAGE100
指
Figure 966288DEST_PATH_IMAGE101
网格单元的相邻迭代时间步电磁守恒变量差值；
Figure 205640DEST_PATH_IMAGE102
是指相邻
Figure 769476DEST_PATH_IMAGE103
网格单元分裂后的系数矩阵；
Figure 891016DEST_PATH_IMAGE104
是指相邻
Figure 260817DEST_PATH_IMAGE105
网格单元分裂后的系数矩阵；
Figure 354675DEST_PATH_IMAGE106
是指相邻
Figure 587948DEST_PATH_IMAGE107
网格单元分裂后的系数矩阵；
Figure 196784DEST_PATH_IMAGE108
是指相邻
Figure 42380DEST_PATH_IMAGE109
网格单元分裂后的系数矩阵；
Figure 53062DEST_PATH_IMAGE110
是指相邻
Figure 21018DEST_PATH_IMAGE111
网格单元分裂后的系数矩阵；
Figure 54833DEST_PATH_IMAGE112
是指相邻
Figure 766437DEST_PATH_IMAGE113
网格单元分裂后的系数矩阵；
最后，得到前后向迭代计算电磁场差值
Figure 569308DEST_PATH_IMAGE114
：
Figure 708165DEST_PATH_IMAGE115
Figure 727811DEST_PATH_IMAGE116
Figure 243106DEST_PATH_IMAGE117
其中，
Figure 900484DEST_PATH_IMAGE118
是对角矩阵
Figure 210242DEST_PATH_IMAGE119
的逆矩阵，
Figure 218650DEST_PATH_IMAGE120
、
Figure 272056DEST_PATH_IMAGE121
分别是根据
Figure 49519DEST_PATH_IMAGE122
、
Figure 530179DEST_PATH_IMAGE123
计算得出的上三角矩阵、下三角矩阵；
前向循环：
Figure 524418DEST_PATH_IMAGE124
后向循环：
Figure 115936DEST_PATH_IMAGE125
其中，
Figure 810223DEST_PATH_IMAGE126
是电磁守恒变量差值的中间过渡变量。
5.根据权利要求1所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述步骤2中: 网格密度保证每波长13-20个网格点，壁面密度>300点/波长，几何奇异处加密到50-100个网格点/波长，二维网格在垂直该二维网格所在平面按右手法则推进一层，作为三维问题特例统一计算；网格数据文件包括结构网格块数目及每块三个曲线坐标系下维度。
6.根据权利要求1所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述步骤3中：如有等离子体外部流畅情况，则还应输入对应流场参数。
7.根据权利要求2所述的一种基于隐式FVFD计算频域电磁场的方法，其特征在于：所述定常虚拟时间步循环为隐式，其CFL数不受显式稳定性要求约束。
CN202110526876.XA 2021-05-14 2021-05-14 一种基于隐式fvfd计算频域电磁场的方法 Active CN113158527B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110526876.XA CN113158527B (zh) 2021-05-14 2021-05-14 一种基于隐式fvfd计算频域电磁场的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110526876.XA CN113158527B (zh) 2021-05-14 2021-05-14 一种基于隐式fvfd计算频域电磁场的方法

Publications (2)

Publication Number Publication Date
CN113158527A true CN113158527A (zh) 2021-07-23
CN113158527B CN113158527B (zh) 2021-08-24

Family

ID=76875043

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110526876.XA Active CN113158527B (zh) 2021-05-14 2021-05-14 一种基于隐式fvfd计算频域电磁场的方法

Country Status (1)

Country Link
CN (1) CN113158527B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114494650A (zh) * 2022-04-06 2022-05-13 中国空气动力研究与发展中心计算空气动力研究所 一种分布式非结构网格跨处理器面对接方法及系统
CN116050303A (zh) * 2023-03-06 2023-05-02 中国空气动力研究与发展中心计算空气动力研究所 一种cfd并行计算下的周期性边界条件施加方法
CN117452081A (zh) * 2023-12-26 2024-01-26 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端
CN117610386A (zh) * 2024-01-24 2024-02-27 浙江电驱动创新中心有限公司 基于有限体积法的高频电磁场仿真方法、系统及计算机

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102725654A (zh) * 2010-01-27 2012-10-10 丰田自动车株式会社 雷达系统和方向检测方法
US20150103873A1 (en) * 2013-10-11 2015-04-16 Dell Products L.P. Systems and methods for measurement of electrical channel loss
CN108460241A (zh) * 2018-05-28 2018-08-28 中国民用航空中南地区空中交通管理局 一种仪表着陆系统扰动仿真方法
CN109376485A (zh) * 2018-12-03 2019-02-22 上海无线电设备研究所 基于aca-mlfma加速的区域分解非共形网格的快速仿真建模方法
CN110210129A (zh) * 2019-06-03 2019-09-06 中南大学 自适应有限元gpr频率域正演方法
CN111143991A (zh) * 2019-12-25 2020-05-12 国网辽宁省电力有限公司沈阳供电公司 一种介质包裹导线的横磁波传输模型及其构建方法
CN111899330A (zh) * 2020-07-21 2020-11-06 深圳大学 一种层状结构的全各向异性媒质的仿真方法及相关组件
CN112286886A (zh) * 2020-10-29 2021-01-29 中国空气动力研究与发展中心计算空气动力研究所 一种多块结构网格数据压缩存储格式

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102725654A (zh) * 2010-01-27 2012-10-10 丰田自动车株式会社 雷达系统和方向检测方法
US20150103873A1 (en) * 2013-10-11 2015-04-16 Dell Products L.P. Systems and methods for measurement of electrical channel loss
CN108460241A (zh) * 2018-05-28 2018-08-28 中国民用航空中南地区空中交通管理局 一种仪表着陆系统扰动仿真方法
CN109376485A (zh) * 2018-12-03 2019-02-22 上海无线电设备研究所 基于aca-mlfma加速的区域分解非共形网格的快速仿真建模方法
CN110210129A (zh) * 2019-06-03 2019-09-06 中南大学 自适应有限元gpr频率域正演方法
CN111143991A (zh) * 2019-12-25 2020-05-12 国网辽宁省电力有限公司沈阳供电公司 一种介质包裹导线的横磁波传输模型及其构建方法
CN111899330A (zh) * 2020-07-21 2020-11-06 深圳大学 一种层状结构的全各向异性媒质的仿真方法及相关组件
CN112286886A (zh) * 2020-10-29 2021-01-29 中国空气动力研究与发展中心计算空气动力研究所 一种多块结构网格数据压缩存储格式

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
W. KYLE ANDERSON .ETC: ""Petrov–Galerkin and discontinuous-Galerkin methods for time-domain and frequency-domain electromagnetic simulations"", 《JOURNAL OF COMPUTATIONAL PHYSICS》 *
许勇 等: ""双时间步时域有限体积方法计算时变电磁场"", 《电波科学学报》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114494650A (zh) * 2022-04-06 2022-05-13 中国空气动力研究与发展中心计算空气动力研究所 一种分布式非结构网格跨处理器面对接方法及系统
CN114494650B (zh) * 2022-04-06 2022-06-24 中国空气动力研究与发展中心计算空气动力研究所 一种分布式非结构网格跨处理器面对接方法及系统
CN116050303A (zh) * 2023-03-06 2023-05-02 中国空气动力研究与发展中心计算空气动力研究所 一种cfd并行计算下的周期性边界条件施加方法
CN116050303B (zh) * 2023-03-06 2023-06-27 中国空气动力研究与发展中心计算空气动力研究所 一种cfd并行计算下的周期性边界条件施加方法
CN117452081A (zh) * 2023-12-26 2024-01-26 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端
CN117452081B (zh) * 2023-12-26 2024-04-30 国网天津市电力公司营销服务中心 一种电磁干扰计算方法、装置及存储介质、电子终端
CN117610386A (zh) * 2024-01-24 2024-02-27 浙江电驱动创新中心有限公司 基于有限体积法的高频电磁场仿真方法、系统及计算机

Also Published As

Publication number Publication date
CN113158527B (zh) 2021-08-24

Similar Documents

Publication Publication Date Title
CN113158527B (zh) 一种基于隐式fvfd计算频域电磁场的方法
CN112989680B (zh) 减少网格使用量的fvfd远场积分边界条件计算方法
Özgün et al. MATLAB-based finite element programming in electromagnetic modeling
Natarajan et al. Convergence and accuracy of displacement based finite element formulations over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation
Wang et al. Three-dimensional simple conformal symplectic particle-in-cell methods for simulations of high power microwave devices
CN113158492B (zh) 一种时变电磁场的全隐式双时间步计算方法
US8510091B1 (en) Domain decomposition formulations for simulating electromagnetic fields
CN111079278B (zh) 三维时域杂交间断伽辽金方法外加电磁源项的处理方法
CN102156764A (zh) 一种分析天线辐射和电磁散射的多分辨预条件方法
CN115600435A (zh) 一种介质涂覆导体复合目标电磁散射隐式计算方法及装置
Pan et al. An efficient high order multilevel fast multipole algorithm for electromagnetic scattering analysis
Hejranfar et al. Application of a preconditioned high‐order accurate artificial compressibility‐based incompressible flow solver in wide range of Reynolds numbers
Wang et al. A new family of exponential-based high-order DGTD methods for modeling 3-D transient multiscale electromagnetic problems
Hu et al. A Chebyshev-based high-order-accurate integral equation solver for Maxwell’s equations
CN106777472B (zh) 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法
CN107944113A (zh) 一种计算三维高速平动目标电磁散射场的方法
CN108090296B (zh) 基于高阶辛紧致格式的波导全波分析方法
Gorgizadeh et al. Eigenmode computation of cavities with perturbed geometry using matrix perturbation methods applied on generalized eigenvalue problems
Akçelik et al. Adjoint methods for electromagnetic shape optimization of the low-loss cavity for the international linear collider
Zheng et al. Three dimensional acoustic shape sensitivity analysis by means of adjoint variable method and fast multipole boundary element approach
Gao et al. Element differential method for solving linear and nonlinear electromagnetic problems
CN111144013B (zh) 高精度介质体目标散射的仿真方法
Li et al. A modified dual-level fast multipole boundary element method for large-scale three-dimensional potential problems
CN116401921B (zh) 一种各项异性磁化等离子体媒质处理方法及系统
Garza et al. High-order Chebyshev-based Nyström Methods for electromagnetics

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