CN116227185A - 一种高精度快速计算随机磁化等离子体电磁统计特性的方法 - Google Patents

一种高精度快速计算随机磁化等离子体电磁统计特性的方法 Download PDF

Info

Publication number
CN116227185A
CN116227185A CN202310143038.3A CN202310143038A CN116227185A CN 116227185 A CN116227185 A CN 116227185A CN 202310143038 A CN202310143038 A CN 202310143038A CN 116227185 A CN116227185 A CN 116227185A
Authority
CN
China
Prior art keywords
calculating
electromagnetic
updating
field component
directions
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
CN202310143038.3A
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.)
Xiangtan University
Original Assignee
Xiangtan University
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 Xiangtan University filed Critical Xiangtan University
Priority to CN202310143038.3A priority Critical patent/CN116227185A/zh
Publication of CN116227185A publication Critical patent/CN116227185A/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/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本发明公开了一种基于相关系数法的随机加权拉盖尔时域差分有限方法，具体按照以下步骤实施：输入模型文件；初始化参数和设置参数；添加场源到y方向上的电场分量系数中；更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure DDA0004088166810000011
更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure DDA0004088166810000012
更新计算整个计算区域的x,y,z方向上电流密度分量系数均值
Figure DDA0004088166810000013
更新计算整个计算区域的电磁场均值分量系数的辅助变量；使用MC‑CC方法计算所需的相关系数；更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure DDA0004088166810000014
更新计算整个计算区域的x,y,z方向上磁场分量系数标准差
Figure DDA0004088166810000015
更新计算整个计算区域的x,y,z方向上电流密度分量系数标准差
Figure DDA0004088166810000016
更新计算整个计算区域的电磁场分量系数标准差的辅助变量；更新计算观测点处的电磁场分量均值及标准差；判断拉盖尔多项式的阶数q是否达到预设值。本发明的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，计算速度快、精度高，且更适用于随机媒质中的电磁统计特性分析。

Description

一种高精度快速计算随机磁化等离子体电磁统计特性的方法
技术领域
本发明属于计算电磁学技术领域，具体涉及一种高精度快速计算随机磁化等离子体电磁统计特性的方法。
背景技术
众所周知，时域有限差分(Finite-difference time-domain，FDTD)方法在处理随机媒质时，需要考虑媒质中参数的随机性，难度较大。在分析随机媒质中电磁统计特性时，传统方法是蒙特卡洛(Monte Carlo，MC)方法，该方法需要大量的模拟，从而耗时很长，效率较低。因此，研究者们基于FDTD方法，提出了用于分析随机媒质电磁统计特性的随机时域有限差分方法(Stochastic Finite-difference time-domain，S-FDTD)，这种方法只需要一次计算便可以得到所需的电磁场均值和方差，在实际应用中效率很高。但是这种方法是显式迭代的方法，其时间步长受限于Courant-Friedrich-Levy(CFL)稳定性条件，导致仿真时间较长，在精细结构上，特别是磁化等离子体中，CFL约束条件可能更严格。
因此学者们提出了一种S-WLP-FDTD方法，这种方法不需要处理时间步长，同时还可以对随机媒质进行电磁统计特性分析，这使得它的计算效率比传统的S-FDTD方法在解决具有精细结构的随机问题方面效率更高。但是目前所提出的S-WLP-FDTD方法仅将等离子体的电子密度作为随机各向同性冷等离子体的随机变量，并未涉及磁化等离子体随机媒质统计特性的计算，而且没有提供确定不同建模场景的互相关值的最佳方法，这将影响算法的计算精度，导致算法计算精度较低。
发明内容
本发明的目的是提供一种高精度快速计算随机磁化等离子体电磁统计特性的方法，这种方法是将基于蒙特卡洛相关系数(Monte Carlo Correlation Coefficient,MC-CC)法的相关系数近似值引入到S-WLP-FDTD算法中，使得算法计算速度较快、精度高，然后使用这种算法对随机磁化等离子体中统计特性进行分析。
本发明所采用的技术方案是，一种高精度快速计算随机磁化等离子体电磁统计特性的方法，具体按照以下步骤实施：
步骤1：输入模型文件；
步骤2：初始化参数和设置参数；
步骤3：添加场源到y方向上的电场分量系数中；
步骤4：更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure BDA0004088166790000021
步骤5：更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure BDA0004088166790000022
步骤6：更新计算整个计算区域的x,y,z方向上电流密度分量系数均值
Figure BDA0004088166790000023
步骤7：更新计算整个计算区域的电磁场分量系数均值的辅助变量；
步骤8：使用MC-CC方法计算所需的相关系数；
步骤9：更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure BDA0004088166790000031
步骤10：更新计算整个计算区域的x,y,z方向上磁场分量系数标准差
Figure BDA0004088166790000032
步骤11：更新计算整个计算区域的x,y,z方向上电流密度分量系数标准差
Figure BDA0004088166790000033
步骤12：更新计算整个计算区域的电磁场分量系数标准差的辅助变量；
步骤13：更新计算观测点处的电磁场分量均值及标准差；
步骤14：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。
本发明的特点还在于：
步骤1输入模型文件，具体为：
计算区域大小Nx×Ny×Nz，其中Nx为x方向的网格数，Ny为y方向的网格数；Nz为z方向的网格数；空间步长Δζ，ζ＝x,y,z，x为横坐标，y为纵坐标，z为竖坐标；时间步长Δt；真空中的电导率σ，磁导率μ0，介电常数ε0；吸收边界SC-PML与相关参数κζmax，σζmax，αζmax；其中，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax/σopt取值范围为(0，12]；仿真计算时长Tf；迭代次数k，k≥0且为整数；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[109，1013]；观测点；场源参数。
步骤2初始化参数和设置参数具体为：
初始化的参数包括：
将整个计算区域的电磁场系数均值
Figure BDA0004088166790000041
整个计算区域的电磁场系数标准差
Figure BDA0004088166790000042
整个计算区域的电磁场分量系数的和
Figure BDA0004088166790000043
整个计算区域的辅助变量
Figure BDA0004088166790000044
拉盖尔多项式
Figure BDA0004088166790000045
均初始化为零，其中Fη＝Eζ,Hζ，ζ＝x,y,z，
Figure BDA0004088166790000046
初始化PML系数(c1ζ,c2ζ)，具体为：
c1ζ＝1/(1+0.5ε0s)
c2ζ＝0
其中，ζ＝x,y,z，ε0是空气中的介电常数，s为时间尺度因子，s取值范围为[109，1013]；
设置的参数具体为：
设置CFS-PML吸收边界的参数，具体为：
σζ＝σζmax|ζ-ζ0|m/dm
κζ＝1+(κζmax-1)|ζ-ζ0|m/dm
αζ＝αζmaxζ0/d
式中ζ＝x,y,z，ζ0为PML层与非PML截面位置，d是PML吸收边界的厚度，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax根据σopt来设置，σζmax/σopt取值范围为(0，12]；σopt＝(m+1)/150πΔζ，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δζ取值范围
Figure BDA0004088166790000047
λ为源的波长；
设置PML系数，具体为：
Figure BDA0004088166790000053
步骤3中所添加的场源的表达式为：
Ey(t)＝exp(-(t-t0)2/τ2)
其中，t0,τ为场源参数。
步骤4更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure BDA0004088166790000051
ζ＝x,y,z；具体更新公式为：
Figure BDA0004088166790000052
Figure BDA0004088166790000061
Figure BDA0004088166790000071
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000072
Figure BDA0004088166790000081
Figure BDA0004088166790000082
Figure BDA0004088166790000083
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤5更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure BDA0004088166790000084
ζ＝x,y,z，具体更新公式为：
Figure BDA0004088166790000085
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000086
Figure BDA0004088166790000091
Figure BDA0004088166790000092
Figure BDA0004088166790000093
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤6更新计算整个计算区域的x,y,z方向上的电流密度分量系数均值
Figure BDA0004088166790000094
ζ＝x,y,z，具体更新公式为：
Figure BDA0004088166790000095
Figure BDA0004088166790000096
Figure BDA0004088166790000097
其中，p1＝2ωb/(s+2v),
Figure BDA0004088166790000098
i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
s是时间尺度因子，v是等离子体碰撞频率，ωp是等离子体频率，ε0是相对介电常数。
步骤7更新计算整个计算区域的电磁场分量系数均值的辅助变量，具体更新公式为：
Figure BDA0004088166790000101
Figure BDA0004088166790000102
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
步骤8使用MC-CC方法计算所需的相关系数，具体公式如下：
Figure BDA0004088166790000103
步骤9更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure BDA0004088166790000104
ζ＝x,y,z；具体更新公式为：
Figure BDA0004088166790000111
Figure BDA0004088166790000121
Figure BDA0004088166790000131
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000141
Figure BDA0004088166790000142
Figure BDA0004088166790000143
Figure BDA0004088166790000144
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤10更新计算整个计算区域的x,y,z方向上磁场分量系数标准差
Figure BDA0004088166790000145
ζ＝x,y,z；具体更新公式为：
Figure BDA0004088166790000146
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000151
Figure BDA0004088166790000152
Figure BDA0004088166790000153
Figure BDA0004088166790000154
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤11更新计算整个计算区域的x,y,z方向上电流密度分量系数标准差
Figure BDA0004088166790000155
ζ＝x,y,z,具体更新公式为：
Figure BDA0004088166790000161
Figure BDA0004088166790000162
Figure BDA0004088166790000163
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
s是时间尺度因子，v是等离子体碰撞频率，ωp是等离子体频率，ε0是相对介电常数。
Figure BDA0004088166790000171
步骤12更新计算整个计算区域的电磁场分量系数标准差的辅助变量，具体更新公式为：
Figure BDA0004088166790000172
Figure BDA0004088166790000173
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
步骤13更新计算观测点处的电磁场分量均值及标准差，具体更新公式为：
Figure BDA0004088166790000174
其中，U表示电磁场分量Ex,Ey,Hz，Uq表示q阶电磁场分量系数，
Figure BDA0004088166790000175
是q阶加权拉盖尔多项式，
Figure BDA0004088166790000176
是带有时间尺度因子s＞0的扩展时间，
Figure BDA0004088166790000177
是q阶拉盖尔多项式。
步骤14：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。
本发明的有益效果是：
①本发明一种高精度快速计算随机磁化等离子体电磁统计特性的方法，在直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯韦方程，使得在更新计算整个计算区域的电磁场分量系数时不涉及到时间步长，只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程中时间步长可以取得比柯西稳定性条件限制的时间步长更大，使得计算速度更快；
②本发明一种高精度快速计算随机磁化等离子体电磁统计特性的方法，在求解随机磁化等离子体媒质中电磁场分量系数时，对电磁场分量的均值和方差进行了求解，使得它在计算时比传统的蒙特卡洛方法更简单、计算速度更快，而且可以对大区域的电磁场问题进行求解；
③本发明一种高精度快速计算随机磁化等离子体电磁统计特性的方法，由于使用了蒙特卡洛相关系数法，计算精度得到了较大的提高。
附图说明
图1是本发明所用方法的流程图；
图2是本发明的方法与传统蒙特卡洛方法和使用不同相关系数的FDTD方法在观测点处电场标准差的时域波形对比图；
图3是本发明的方法与使用MC-CC方法的S-ADE-FDTD方法在观测点处电场标准差的误差系数对比图；
图4是本发明的方法和其他几种方法的计算时间对比；
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明是一种高精度快速计算随机磁化等离子体电磁统计特性的方法，所依据的原理为：首先导出在磁化等离子体中电磁场所满足的麦克斯韦方程；然后利用一种随机WLP-FDTD方法推导出电磁场分量均值和方差的更新方程；最后求解观测点处的电磁场分量。
在求电磁场分量系数更新方程时，首先需要推导出磁化等离子体中电磁场所满足的麦克斯韦方程；
在SC扩展坐标系下，波在磁化等离子体中传播的麦克斯韦方程可写为：
Figure BDA0004088166790000191
其中，E表示电场强度，H表示磁场强度，j为虚数单位，ωp为等离子体频率，v是等离子体碰撞频率，ωb是电子回旋频率，看作一个常数，μ0为真空磁导率，ε0为真空介电常数，
Figure BDA0004088166790000192
为修正后的微分算子，可以写成：
Figure BDA0004088166790000193
sx、sy和sz是坐标扩展变量，可以表示成：
sζ＝kζ+σζ/(αζ+jωε) (3)
其中ζ表示x、y、z，kζ、σζ和αζ为PML的有关的参数。
假设外部磁场为z轴上，于是(1)式中第三式为：
Figure BDA0004088166790000201
Figure BDA0004088166790000202
Figure BDA0004088166790000203
然后，求出电磁场分量系数的更新方程；
为了计算方便，引入下面几个辅助变量：
Figure BDA0004088166790000204
Figure BDA0004088166790000205
将(3)代入(5)，然后利用jω→t的变换，可以得到十二组方程，这里仅给出第一个方程：
Figure BDA0004088166790000206
由于电磁场分量和其对时间的一阶偏导可以展开成一系列的电磁场分量系数与加权拉盖尔多项式的函数之和，公式如下：
Figure BDA0004088166790000207
上式中U表示电磁场分量Eζ,Hζ，Uq表示q阶电磁场分量系数，
Figure BDA0004088166790000208
是q阶加权拉盖尔多项式，
Figure BDA0004088166790000209
是带有时间尺度因子s＞0的扩展时间，
Figure BDA00040881667900002010
是q阶拉盖尔多项式。将(7)代入(1)中，然后让方程两边同乘以
Figure BDA00040881667900002011
可以得到：
Figure BDA0004088166790000211
Figure BDA0004088166790000212
Figure BDA0004088166790000213
Figure BDA0004088166790000214
Figure BDA0004088166790000215
Figure BDA0004088166790000216
上面式中，q是加权拉盖尔多项式的阶数，Dx、Dy和Dz分别是沿x、y和z方向上的微分算子，
Figure BDA0004088166790000217
和
Figure BDA0004088166790000218
是q阶电场分量系数，
Figure BDA0004088166790000219
和
Figure BDA00040881667900002110
是q阶磁场分量系数，Ciζ,i＝1,2；ζ＝x,y,z是与坐标网格有关的PML系数，计算式为：
Figure BDA00040881667900002111
Figure BDA00040881667900002112
和
Figure BDA00040881667900002113
是拉盖尔域中电磁场分量和辅助变量的低阶和，公式如下：
Figure BDA0004088166790000221
上式中的辅助变量
Figure BDA0004088166790000222
的计算式为：
Figure BDA0004088166790000223
在磁化等离子体中，WLPs域的极化电流密度为：
Figure BDA0004088166790000224
Figure BDA0004088166790000225
Figure BDA0004088166790000226
其中
Figure BDA0004088166790000227
将式(19)代入到(10-12)式，可得
Figure BDA0004088166790000228
Figure BDA0004088166790000229
Figure BDA0004088166790000231
Figure BDA0004088166790000232
其中
Figure BDA0004088166790000233
Figure BDA0004088166790000234
Figure BDA0004088166790000235
Figure BDA0004088166790000236
然后是对电磁场方差分量进行求解，首先，对式(10)取方差，有：
Figure BDA0004088166790000237
其中相关系数近似为：
Figure BDA0004088166790000238
因此有：
Figure BDA0004088166790000241
其中
Figure BDA0004088166790000242
是
Figure BDA0004088166790000243
的标准差，有
Figure BDA0004088166790000244
辅助变量的计算公式为：
Figure BDA0004088166790000245
接下来我们对极化电流密度J的更新方程进行推导，假设
Figure BDA0004088166790000248
(29)
将式(29)代入到式(22)中，得到：
Figure BDA0004088166790000246
Figure BDA0004088166790000249
(30)
其中，相关系数应用以下近似值：
Figure BDA0004088166790000247
由此可得极化电流密度的标准差方程为：
Figure BDA0004088166790000251
Figure BDA0004088166790000252
其中σ{v},σ{ωp}分别是碰撞频率和等离子体频率的标准差。
将式(32)和式(33)代入到式(26)中，可得到电场的标准差方程为：
Figure BDA0004088166790000253
Figure BDA0004088166790000254
Figure BDA0004088166790000261
Figure BDA0004088166790000262
其中
Figure BDA0004088166790000263
上面式中，我们通过使用S-FDTD方法，对电磁场分量均值和标准差的方程进行推导，完成了对随机磁化等离子体中电磁分量的计算；同时使用了WLP-FDTD方法，既能消除柯西稳定性条件的限制，而且又能解决使用较大的时间步长时产生很大的色散误差这个难题。而且使用相关系数法(MC-CC)提高了算法的精度，最后通过(9)式解得观测点的电磁场分量。
本发明提供一种高精度快速计算随机磁化等离子体电磁统计特性的方法，具体按照以下步骤实施：
步骤1输入模型文件，具体为：
计算区域大小Nx×Ny×Nz，其中Nx为x方向的网格数，Ny为y方向的网格数；Nz为z方向的网格数；空间步长Δζ，ζ＝x,y,z，x为横坐标，y为纵坐标，z为竖坐标；时间步长Δt；真空中的电导率σ，磁导率μ0，介电常数ε0；吸收边界SC-PML与相关参数κζmax，σζmax，αζmax；其中，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax/σopt取值范围为(0，12]；仿真计算时长Tf；迭代次数k，k≥0且为整数；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[109，1013]；观测点；场源参数。
步骤2初始化参数和设置参数具体为：
初始化的参数包括：
将整个计算区域的电磁场均值系数
Figure BDA0004088166790000271
整个计算区域的电磁场标准差系数
Figure BDA0004088166790000272
整个计算区域的电磁场分量系数的和
Figure BDA0004088166790000273
整个计算区域的辅助变量
Figure BDA0004088166790000274
拉盖尔多项式
Figure BDA0004088166790000275
均初始化为零，其中Fη＝Eζ,Hζ，ζ＝x,y,z，
Figure BDA0004088166790000276
初始化PML系数(c1ζ,c2ζ)，具体为：
c1ζ＝1/(1+0.5ε0s)
c2ζ＝0
其中，ζ＝x,y,z，ε0是空气中的介电常数，s为时间尺度因子，s取值范围为[109，1013]；
设置的参数具体为：
设置CFS-PML吸收边界的参数，具体为：
σζ＝σζmax|ζ-ζ0|m/dm
κζ＝1+(κζmax-1)|ζ-ζ0|m/dm
αζ＝αζmaxζ0/d
式中ζ＝x,y,z，ζ0为PML层与非PML截面位置，d是PML吸收边界的厚度，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax根据σopt来设置，σζmax/σopt取值范围为(0，12]；σopt＝(m+1)/150πΔζ，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δζ取值范围
Figure BDA0004088166790000281
λ为源的波长；
设置PML系数，具体为：
Figure BDA0004088166790000282
步骤3中所添加的场源的表达式为：
Ey(t)＝exp(-(t-t0)2/τ2)
其中，t0,τ为场源参数。
步骤4更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure BDA0004088166790000291
具体更新公式为：
Figure BDA0004088166790000292
Figure BDA0004088166790000301
Figure BDA0004088166790000311
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000312
Figure BDA0004088166790000321
Figure BDA0004088166790000322
Figure BDA0004088166790000323
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤5更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure BDA0004088166790000324
ζ＝x,y,z，具体更新公式为：
Figure BDA0004088166790000325
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000326
Figure BDA0004088166790000331
Figure BDA0004088166790000332
Figure BDA0004088166790000333
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤6更新计算整个计算区域的x,y,z方向上的电流密度分量系数均值
Figure BDA0004088166790000334
ζ＝x,y,z，具体更新公式为：
Figure BDA0004088166790000335
Figure BDA0004088166790000336
Figure BDA0004088166790000337
其中，
Figure BDA0004088166790000341
c1ζ＝1/(κζαζ+σζ+0.5κζε0s),c2ζ＝(2αζ/ε0s+1)c1ζ,c3＝ε0/μ0，c4＝2/(ε0s)，
i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
步骤7更新计算整个计算区域的电磁场分量系数均值的辅助变量，具体更新公式为：
Figure BDA0004088166790000342
Figure BDA0004088166790000343
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
步骤8使用MC-CC方法计算所需的相关系数，具体公式如下：
Figure BDA0004088166790000344
步骤9更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure BDA0004088166790000345
ζ＝x,y,z；具体更新公式为：
Figure BDA0004088166790000351
Figure BDA0004088166790000361
Figure BDA0004088166790000371
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000381
Figure BDA0004088166790000382
Figure BDA0004088166790000383
Figure BDA0004088166790000384
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤10更新计算整个计算区域的x,y,z方向上磁场分量系数标准差
Figure BDA0004088166790000385
ζ＝x,y,z；具体更新公式为：
Figure BDA0004088166790000386
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure BDA0004088166790000391
Figure BDA0004088166790000392
Figure BDA0004088166790000393
Figure BDA0004088166790000394
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
步骤11更新计算整个计算区域的x,y,z方向上电流密度分量系数标准差
Figure BDA0004088166790000395
ζ＝x,y,z,具体更新公式为：
Figure BDA0004088166790000401
Figure BDA0004088166790000402
Figure BDA0004088166790000403
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
s是时间尺度因子，v是等离子体碰撞频率，ωp是等离子体频率，ε0是相对介电常数。
Figure BDA0004088166790000411
步骤12更新计算整个计算区域的电磁场分量系数标准差的辅助变量，具体更新公式为：
Figure BDA0004088166790000412
Figure BDA0004088166790000413
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
步骤13更新计算观测点处的电磁场分量均值及标准差，具体更新公式为：
Figure BDA0004088166790000414
其中，U表示电磁场分量Ex,Ey,Hz，Uq表示q阶电磁场分量系数，
Figure BDA0004088166790000415
是q阶加权拉盖尔多项式，
Figure BDA0004088166790000416
是带有时间尺度因子s＞0的扩展时间，
Figure BDA0004088166790000417
是q阶拉盖尔多项式。
步骤14：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。
实施例
对非均匀随机磁化等离子体中的电波仿真计算
按照本发明的方法步骤进行实施，如图2所示，实验中整个计算区域为600×600网格，网格是非均匀划分的。四个边界采用10层网格的PML吸收边界，每个网格的尺寸为150μm，空气在计算区域的x和y方向都在第11-120,481-590个网格里，空气的网格尺寸和PML吸收边界一样，而中间的等离子体分为三部分，左右两边等离子体的参数一样，分别在第121-240,361-480个网格里，厚度为3.75mm，其z轴空间步长Δz＝37.5μm，两边等离子体的具体参数为：
μv＝1.2×1010rad/s,μωp＝2.625×1011rad/s,σ{v}＝1.2×108rad/s,
σ{ωp}＝2.625×109rad/s,ωb＝3×1011rad/s
其中μv,μωp分别为碰撞频率和等离子体频率的平均值。
中间等离子体位于第241-360个网格里，厚度为1.5mm,其z轴空间步长Δz＝9.375μm，中间等离子体的具体参数为：
μv＝2×1010rad/s,μωp＝3.142×1011rad/s,σ{v}＝2×108rad/s,
σ{ωp}＝3.142×109rad/s,ωb＝3×1011rad/s
其中μv,μωp分别为碰撞频率和等离子体频率的平均值。
计算中所加的源位于第50个网格处，所加场源的表达式如下：
Ey(t)＝exp(-(t-t0)2/τ2)
其中，t0＝20ps,,τ＝5ps。时间步长Δt＝0.5ps,加权拉盖尔多项式的阶数q＝200，时间扩展因子s＝1.256×1012，整个仿真时间为Tf＝0.5ns。
采用本发明方法计算的观测点处的电场分量标准差σ{Ey}与采用传统MC方法及使用不同相关系数的FDTD方法计算的结果参见图2，从图2中可见本发明方法与传统MC方法计算结果一致，同时精度高于其他使用相关系数为0.2和0.5时的方法，验证了本发明方法的正确性。
图3是本发明的方法与使用MC-CC方法的S-ADE-FDTD方法在观测点处电场标准差的误差系数对比图，从图3中可见，本发明方法比使用相关系数法的的S-ADE-FDTD方法的计算精度更高；图4为几种方法所需的运行时间对比，可以明显看到使用了MC-CC方法的S-WLP-FDTD算法，在保持了精度的情况下，大大减少了所需的运行时间，验证了本发明方法的有效性。

Claims (14)

1.一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，具体按照以下步骤实施：
步骤1：输入模型文件；
步骤2：初始化参数和设置参数；
步骤3：添加场源到y方向上的电场分量系数中；
步骤4：更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure FDA0004088166780000011
步骤5：更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure FDA0004088166780000012
步骤6：更新计算整个计算区域的x,y,z方向上电流密度分量系数均值
Figure FDA0004088166780000013
步骤7：更新计算整个计算区域的电磁场分量系数均值的辅助变量；
步骤8：使用MC-CC方法计算所需的相关系数；
步骤9：更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure FDA0004088166780000014
步骤10：更新计算整个计算区域的x,y,z方向上磁场分量系数标准差
Figure FDA0004088166780000015
步骤11：更新计算整个计算区域的x,y,z方向上电流密度分量系数标准差
Figure FDA0004088166780000016
步骤12：更新计算整个计算区域的电磁场分量系数标准差的辅助变量；
步骤13：更新计算观测点处的电磁场分量均值及标准差；
步骤14：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。
2.根据权利要求1所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤1输入模型文件，具体为：计算区域大小Nx×Ny×Nz，其中Nx为x方向的网格数，Ny为y方向的网格数；Nz为z方向的网格数；空间步长Δζ，ζ＝x,y,z，x为横坐标，y为纵坐标，z为竖坐标；时间步长Δt；真空中的电导率σ，磁导率μ0，介电常数ε0；吸收边界SC-PML与相关参数κζmax，σζmax，αζmax；其中，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax/σopt取值范围为(0，12]；仿真计算时长Tf；迭代次数k，k≥0且为整数；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[109，1013]；观测点；场源参数。
3.根据权利要求1所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤2初始化参数和设置参数具体为：
初始化的参数包括：
将整个计算区域的电磁场均值系数
Figure FDA0004088166780000021
整个计算区域的电磁场标准差系数
Figure FDA0004088166780000022
整个计算区域的电磁场分量系数的和
Figure FDA0004088166780000023
整个计算区域的辅助变量
Figure FDA0004088166780000024
拉盖尔多项式
Figure FDA0004088166780000025
均初始化为零，其中Fη＝Eζ,Hζ，ζ＝x,y,z，
Figure FDA0004088166780000031
初始化PML系数(c1ζ,c2ζ)，具体为：
c1ζ＝1/(1+0.5ε0s)
c2ζ＝0
其中，ζ＝x,y,z，ε0是空气中的介电常数，s为时间尺度因子，s取值范围为[109，1013]；
设置的参数具体为：
设置CFS-PML吸收边界的参数，具体为：
σζ＝σζmax|ζ-ζ0|m/dm
κζ＝1+(κζmax-1)|ζ-ζ0|m/dm
αζ＝αζmaxζ0/d
式中ζ＝x,y,z，ζ0为PML层与非PML截面位置，d是PML吸收边界的厚度，κζmax取整数，κζmax取值范围为[1，60]；αζmax取值范围为[0，1)；σζmax根据σopt来设置，σζmax/σopt取值范围为(0，12]；σopt＝(m+1)/150πΔζ，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δζ取值范围
Figure FDA0004088166780000032
λ为源的波长；
设置PML系数，具体为：
Figure FDA0004088166780000033
4.根据权利要求3所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤3中所添加的场源的表达式为：
Ey(t)＝exp(-(t-t0)2/τ2)
其中，t0,τ为场源参数。
5.根据权利要求4所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤4更新计算整个计算区域的x,y,z方向上电场分量系数均值
Figure FDA0004088166780000041
具体更新公式为：
Figure FDA0004088166780000042
Figure FDA0004088166780000051
Figure FDA0004088166780000061
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure FDA0004088166780000062
Figure FDA0004088166780000071
Figure FDA0004088166780000072
Figure FDA0004088166780000073
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
6.根据权利要求5所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤5更新计算整个计算区域的x,y,z方向上磁场分量系数均值
Figure FDA0004088166780000074
具体更新公式为：
Figure FDA0004088166780000075
Figure FDA0004088166780000076
Figure FDA0004088166780000077
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure FDA0004088166780000081
Figure FDA0004088166780000082
Figure FDA0004088166780000083
Figure FDA0004088166780000084
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
7.根据权利要求6所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤6更新计算整个计算区域的x,y,z方向上的电流密度分量均值
Figure FDA0004088166780000085
具体更新公式为：
Figure FDA0004088166780000091
Figure FDA0004088166780000092
Figure FDA0004088166780000093
其中，
Figure FDA0004088166780000094
c1ζ＝1/(κζαζ+σζ+0.5κζε0s),c2ζ＝(2αζ/ε0s+1)c1ζ,c3＝ε0/μ0，c4＝/(ε0s)，
i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
8.根据权利要求7所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤7更新计算整个计算区域的电磁场分量系数的辅助变量，具体更新公式为：
Figure FDA0004088166780000101
Figure FDA0004088166780000102
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
9.根据权利要求8所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤8使用MC-CC方法计算所需的相关系数，具体公式如下：
Figure FDA0004088166780000103
10.根据权利要求9所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤9更新计算整个计算区域的x,y,z方向上电场分量系数标准差
Figure FDA0004088166780000104
具体更新公式为：
Figure FDA0004088166780000111
Figure FDA0004088166780000121
Figure FDA0004088166780000131
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure FDA0004088166780000141
Figure FDA0004088166780000142
Figure FDA0004088166780000143
Figure FDA0004088166780000144
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
11.根据权利要求10所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤10更新计算整个计算区域的x,y,z方向上磁场分量标准差
Figure FDA0004088166780000145
具体更新公式为：
Figure FDA0004088166780000146
Figure FDA0004088166780000147
Figure FDA0004088166780000148
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
κζ，σζ和αζ是吸收边界SC-PML的相关参数，s是时间尺度因子，ε0和μ0是真空中的磁导率和介电常数。
Figure FDA0004088166780000151
Figure FDA0004088166780000152
Figure FDA0004088166780000153
Figure FDA0004088166780000154
c1ζ＝1/(κζαζ+σζ+0.5κζε0s)，c2ζ＝(2αζ/ε0s+1)c1ζ
c3＝ε0/μ0，c4＝2/(ε0s)
12.根据权利要求11所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤11更新计算整个计算区域的x,y,z方向上电流密度分量标准差
Figure FDA0004088166780000155
具体更新公式为：
Figure FDA0004088166780000161
Figure FDA0004088166780000162
Figure FDA0004088166780000163
其中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；k表示竖坐标上的第k个计算网格。
s是时间尺度因子，v是等离子体碰撞频率，ωp是等离子体频率，ε0是相对介电常数。
p1＝2ωb/(s+2v),
Figure FDA0004088166780000171
13.根据权利要求12所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤12更新计算整个计算区域的电磁场分量系数的辅助变量，具体更新公式为：
Figure FDA0004088166780000172
Figure FDA0004088166780000173
其中，Fη＝Eζ,Hζ；ζ＝x,y,z。
14.根据权利要求13所述的一种高精度快速计算随机磁化等离子体电磁统计特性的方法，其特征在于，所述步骤13更新计算观测点处的电磁场分量均值及标准差，具体更新公式为：
Figure FDA0004088166780000174
其中，U表示电磁场分量Ex,Ey,Hz，Uq表示q阶电磁场分量系数，
Figure FDA0004088166780000175
是q阶加权拉盖尔多项式，
Figure FDA0004088166780000176
是带有时间尺度因子s＞0的扩展时间，
Figure FDA0004088166780000177
是q阶拉盖尔多项式。
CN202310143038.3A 2023-02-21 2023-02-21 一种高精度快速计算随机磁化等离子体电磁统计特性的方法 Pending CN116227185A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310143038.3A CN116227185A (zh) 2023-02-21 2023-02-21 一种高精度快速计算随机磁化等离子体电磁统计特性的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310143038.3A CN116227185A (zh) 2023-02-21 2023-02-21 一种高精度快速计算随机磁化等离子体电磁统计特性的方法

Publications (1)

Publication Number Publication Date
CN116227185A true CN116227185A (zh) 2023-06-06

Family

ID=86572617

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310143038.3A Pending CN116227185A (zh) 2023-02-21 2023-02-21 一种高精度快速计算随机磁化等离子体电磁统计特性的方法

Country Status (1)

Country Link
CN (1) CN116227185A (zh)

Similar Documents

Publication Publication Date Title
Dehghan et al. Application of semi‐analytic methods for the Fitzhugh–Nagumo equation, which models the transmission of nerve impulses
CN110276109B (zh) 一种高超声速飞行器等离子体鞘套电磁特性的仿真方法
CN113158527B (zh) 一种基于隐式fvfd计算频域电磁场的方法
CN111159637B (zh) 一种应用于磁化等离子体计算的电磁波时域精细积分方法
CN111783339B (zh) 电磁波在随机色散介质中传播的pce-fdtd方法
Graves et al. Influence of modeling and simulation on the maturation of plasma technology: Feature evolution and reactor design
Colella et al. A conservative finite difference method for the numerical solution of plasma fluid equations
CN107016184B (zh) 一种二维高精度迭代的非磁化等离子体中的实现方法
CN111800932B (zh) 一种等离子体放电过程模拟方法及系统
Zhang et al. Space and phase resolved ion energy and angular distributions in single-and dual-frequency capacitively coupled plasmas
Wang et al. Sheath thickness evaluation for collisionless or weakly collisional bounded plasmas
CN106777472B (zh) 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法
CN117332658B (zh) 一种各向异性时变等离子体的电磁特性确定方法及系统
CN117787056B (zh) 一种适用于Debye模型电磁场传播的高精度计算方法
CN112307639B (zh) 一种基于高品质算法的贝伦格完全匹配层仿真方法
CN108090296B (zh) 基于高阶辛紧致格式的波导全波分析方法
Wolf et al. A particle-in-cell method for the simulation of plasmas based on an unconditionally stable field solver
CN116227185A (zh) 一种高精度快速计算随机磁化等离子体电磁统计特性的方法
CN116401921B (zh) 一种各项异性磁化等离子体媒质处理方法及系统
CN104820660B (zh) 一种扩展柱坐标系下完全匹配吸收边界的实现方法
Chizhonkov et al. Numerical simulation of the breaking effect in nonlinear axially-symmetric plasma oscillations
Xiao et al. Application of High Dimensional B-Spline Interpolation in Solving the Gyro-Kinetic Vlasov Equation Based on Semi-Lagrangian Method
Wu et al. Theoretical and experimental studies of the bell-jar-top inductively coupled plasma
JP2008031515A (ja) マイクロ波プラズマ解析プログラム
CN116882185A (zh) 一种矩形感性耦合等离子体的三维流体模拟方法及系统

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