CN113887104A - 一种基于mhd的电磁波与磁化等离子体作用数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于MHD的电磁波与磁化等离子体作用数值模拟方法，包括如下步骤：步骤1，确定MHD模型：步骤2，初始化参数：步骤3，基于时域有限差分法更新磁场分量；步骤4，基于高斯消元法更新计算等离子体速度；步骤5，基于约束插值曲线法更新计算等离子体密度分布：步骤6，采用时域有限差分离散格式计算电场分量；步骤7，记录采样点的电磁场、等离子体速度、等离子密度数据；步骤8，将n+1赋值给n，并判断n是否达到预设时间步数，若未达到预设值，则返回步骤3，若达到预设值，则结束。本发明所公开基于MHD的电磁波与磁化等离子体作用数值模拟方法，可实现实时多物理参变量(电场、磁场、等离子体密度和等离子体速度等)的演变特征监测。

Description

一种基于MHD的电磁波与磁化等离子体作用数值模拟方法

技术领域

本发明属于计算电磁学技术领域，特别涉及该领域中的一种基于MHD的电磁波与磁化等离子体作用数值模拟方法。

背景技术

磁流体理论(Magnetohydrodynamics，MHD)是描述高频电磁波与磁化等离子体相互作用的重要理论之一，主要涉及电场、磁场、等离子体(电子和离子)速度、等离子体密度等多物理参数。建立多物理参量的时空离散方案、创建物理参量的更新模拟方法是模拟高频电磁波与磁化等离子体相互作用的难题。考虑到MHD理论涉及麦克斯韦方程，可把时域有限差分方法中的Yee结构作为首选，但仍需定义每个格点的物理性质，赋予不同位置的等离子体特征。同时还需要解决等离子体速度的收敛性问题，以及由于离子体密度更新中对流项的存在而导致求解过程中出现的不稳定现象。因此建立一套基于MHD的电磁波与磁化等离子体作用数值方法一直是等离子体数值模拟的热点及难点。

发明内容

本发明所要解决的技术问题就是提供一种基于MHD的电磁波与磁化等离子体作用数值模拟方法，建立多物理参量的时空离散方案，创建电场、磁场、等离子体速度、等离子密度的更新模拟方法，以解决等离子体速度的收敛性问题以及等离子体密度更新中的不稳定现象，具有稳定和低耗散性的明显优势。

本发明采用如下技术方案：

一种基于MHD的电磁波与磁化等离子体作用数值模拟方法，其改进之处在于，包括如下步骤：

步骤1，确定MHD模型：

步骤11，建立电磁波与磁化等离子体相互作用的控制方程组，如下式：

其中，下标α表示等离子体的种类：电子e或氧离子O⁺，

分别指x，y，z方向的单位向量，

表示磁场扰动，

表示电场扰动，

是随时间变化的等离子体速度，N_α为等离子体的数密度，式(1)中μ是磁导张量，σ^m为磁损耗张量，式(2)中ε为介电张量，σ^e为电损耗，q_e＝-e，

分别表示带电荷量，电流

由电子运动和离子运动形成，即

式(4)中m_α分别为m_e电子质量，

氧离子质量，B₀表示外加磁场，忽略入射电磁波产生的扰动磁场，ν_α为碰撞频率，k_B是玻尔兹曼常数，T_α是等离子体温度；

步骤12，选取二维平面xOz为外加磁场所在的平面，忽略上述物理参数在y方向梯度的变化

ψ为任意参量，只考虑外加磁场平面上梯度的变化，即垂直方向z方向和x方向的梯度变化，将式(1)至(4)展开，如下式：

步骤2，初始化参数，包括确定计算参量时域和空间域的离散化方案和输入模型文件：

步骤21，确定计算参量时域和空间域的离散化方案：

确定时域递推方案采用蛙跳交替的离散方式，更新循环为：

二维空间域是由Nx×Nz个Yee单元构成的网格，以单元的大小作为分辨率，每个场分量的实际位置与标志具有如下关系：

定义E_x、H_x、U_αx的大小为Nx×Nzp1的零矩阵、E_z、H_z、U_αz的大小为Nxp1×Nz的零矩阵，E_y、H_y、U_αy的大小为Nx×Nz的零矩阵，此处Nzp1＝Nz+1，Nxp1＝Nx+1；

步骤22，输入模型文件：

具体输入的参数包括：计算区域大小x方向网格数Nx，z方向网格数Nz；空间步长Δx，Δz；时间步长Δt；真空磁导率μ₀；真空介电系数ε₀；玻尔兹曼常数k_B；带电荷量e；电子质量m_e，离子质量

碰撞频率ν_α；外加磁场B_0x、B_0z；电子密度分布N_e0，离子密度分布N_i0；电磁回旋频率ω_ce；入射电磁波的位置Source_x，Source_z；频率f₀；入射波波长λ；确定电磁波的波形，包括余弦函数调制的高斯波形

τ为以参数决定高斯脉冲在时域和频域的宽度，

为时间移动，余弦波形Wave(t)＝cos(2πf(t-t₀))；确定入射电磁波的极化E_{x_wave}＝Wave(t)tan(φ)，E_{y_wave}＝Wave(t)cot(φ)；电磁波CPML吸收边界，其相关参数PML层的厚度d，σ_max，κ_max，a_max取值范围为[0,0.05]；采样点的坐标为sample_x，sample_z)；

步骤3，基于时域有限差分法更新磁场分量

首次循环中，利用步骤2所得电场分量，非首次循环，电场分量系数为步骤6所得，更新计算整个计算区域的磁场分量

具体为：

步骤4，基于高斯消元法更新计算等离子体速度

首次循环，采用步骤2所得的电场分量，等离子体速度分量，等离子体密度参数，非首次循环，采用步骤5所得电子和离子密度分布，步骤6所得电场分量，基于高斯消元法更新整个计算区域等离子体速度分量

具体为：

将速度更新公式简化为：

化简式(13)左边为：

其中，

同理，化简式(13)右边为：

其中，

由此更新方程简化为：

即：

式中

采用高斯消元法将式(19)的左手边系数矩阵最终化简为三角矩阵，即：

式中，

等离子体速度的更新式为：

步骤5，基于约束插值曲线法更新计算等离子体密度分布

利用步骤4所得等离子体速度分量参数，更新计算整个计算区域的离子体密度分布，在二维空间中逐步在各维度上应用约束插值曲线法，

具体为：首先计算中间变量

为对应的数密度

在x、z方向的偏导，C(x,Δt)为x方向上应用CIP方法，S(x,Δt)为在x方向上采用中心差分所得结果，根据中间变量得到下一时间点等离子体密度的结果：

具体x方向的输运方程为：

化为一般化的运输方程：

f＝N_αx，u＝U_αx，

对上式作偏导，有：

式中，

引入三阶插值：

F_i(x)＝A1_iX³+A2_iX²+g_iX+f_i (30)

X＝x-x_i，假设下边界x_i＝0，因此三阶插值为：

F_i(x)＝A1_ix³+A2_ix²+g_ix+f_i (31)

其中，三阶、二阶系数分别为：

D≡x_iup＝-Δx·sgn(u_i)，下标iup＝i-sgn(u_i)，u_i表示i位置的各方向速度大小，判断函数sgn(u_i)为：

将x＝0,x＝D代入三阶插值函数F_i(x)，得边界条件：

根据对流特征，可得：

F(x,t+Δt)＝F(x-uΔt,t) (36)

因此通过上式可得：

其中，ξ＝-μ_iΔt；

步骤6，采用时域有限差分离散格式计算电场分量

采用步骤3所得等离子体磁场分量参数，步骤4所得等离子体速度参数，以及步骤5所得等离子体密度参数，计算更新电场分量

具体为：

步骤7，记录采样点的电磁场、等离子体速度、等离子密度数据；

步骤8，将n+1赋值给n，并判断n是否达到预设时间步数，若未达到预设值，则返回步骤3，若达到预设值，则结束。

本发明的有益效果是：

本发明所公开基于MHD的电磁波与磁化等离子体作用数值模拟方法，通过对电磁波与磁化等离子体相互作用效应的仿真，可实现实时多物理参变量(电场、磁场、等离子体密度和等离子体速度等)的演变特征监测，使得理论模拟和实验结果比对成为可能，这对深入理解等离子体中波-波相互作用和波-粒相互作用具有重要的理论科学价值。

附图说明

图1是本发明方法的流程示意图；

图2是二维元胞示意图；

图3是等离子体速度更新格点示意图；

图4a是在更新步为10000时，电场的模拟结果示意图；

图4b是在更新步为10000时，电子速度的模拟结果示意图；

图4c是在更新步为10000时，等离子体速度的模拟结果示意图；

图4d是在更新步为10000时，电子密度的模拟结果示意图；

图4e是在更新步为10000时，等离子体密度的模拟结果示意图；

图5是本发明实施例1中x方向的电场谱分析示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明公开了一种基于MHD的电磁波与磁化等离子体作用数值模拟方法，如图1所示，包括如下步骤：

步骤1，确定MHD模型：

步骤11，建立电磁波与磁化等离子体相互作用的控制方程组，如下式：

其中，下标α表示等离子体的种类：电子e或氧离子O⁺，

分别指x，y，z方向的单位向量，

表示磁场扰动，

表示电场扰动，

是随时间变化的等离子体速度，N_α为等离子体的数密度，式(1)中μ是磁导张量，σ^m为磁损耗张量，式(2)中ε为介电张量，σ^e为电损耗，q_e＝-e，

分别表示带电荷量，电流

由电子运动和离子运动形成，即

式(4)中m_α分别为m_e电子质量，

氧离子质量，B₀表示外加磁场，忽略入射电磁波产生的扰动磁场，ν_α为碰撞频率，k_B是玻尔兹曼常数，T_α是等离子体温度；在本实施例中，可将外加磁场设置为0，亦或根据不同物理场景将热压力

或阻力项N_αm_αν_αU_α设置为0。

步骤12，选取二维平面xOz为外加磁场所在的平面，忽略上述物理参数在y方向梯度的变化

ψ为任意参量，只考虑外加磁场平面上梯度的变化，即垂直方向z方向和x方向的梯度变化，将式(1)至(4)展开，如下式：

步骤2，初始化参数，包括确定计算参量时域和空间域的离散化方案和输入模型文件：

步骤21，确定计算参量时域和空间域的离散化方案：

确定时域递推方案采用蛙跳交替的离散方式，更新循环为：

二维空间域是由Nx×Nz个Yee单元构成的网格，图2给出了二维元胞示意图，以单元的大小作为分辨率，每个场分量的实际位置与标志具有如下关系：

定义E_x、H_x、U_αx的大小为Nx×Nzp1的零矩阵、E_z、H_z、U_αz的大小为Nxp1×Nz的零矩阵，E_y、H_y、U_αy的大小为Nx×Nz的零矩阵，此处Nzp1＝Nz+1，Nxp1＝Nx+1；

步骤22，输入模型文件：

具体输入的参数包括：计算区域大小x方向网格数Nx，z方向网格数Nz；空间步长Δx，Δz；时间步长Δt；真空磁导率μ₀；真空介电系数ε₀；玻尔兹曼常数k_B；带电荷量e；电子质量m_e，离子质量m_O+；碰撞频率ν_α；外加磁场B_0x、B_0z；电子密度分布N_e0，离子密度分布N_i0；电磁回旋频率ω_ce；入射电磁波的位置Source_x，Source_z；频率f₀；入射波波长λ；确定电磁波的波形，包括余弦函数调制的高斯波形

τ为以参数决定高斯脉冲在时域和频域的宽度，

为时间移动，余弦波形Wave(t)＝cos(2πf(t-t₀))；确定入射电磁波的极化E_{x_wave}＝Wave(t)tan(φ)，E_{y_wave}＝Wave(t)cot(φ)；电磁波CPML吸收边界，其相关参数PML层的厚度d，σ_max，κ_max，a_max取值范围为[0,0.05]；采样点的坐标为sample_x，sample_z)；

步骤3，基于时域有限差分法更新各磁场分量

首次循环中，利用步骤2所得电场分量，非首次循环，电场分量系数为步骤6所得，更新计算整个计算区域的各磁场分量

具体为：

步骤4，基于高斯消元法更新计算各等离子体速度

首次循环，采用步骤2所得的电场分量，等离子体速度分量，等离子体密度参数，非首次循环，采用步骤5所得电子和离子密度分布，步骤6所得电场分量，基于高斯消元法更新整个计算区域等离子体速度分量

具体为：

将速度更新公式简化为：

等离子体速度更新在空间节点上离散需注意两方面，一方面是高斯消元法更新时，不同方向的速度更新在实际位置中要保持一致；另一方面是热压力项的微分需要与求解的速度位置保持一致。

如图3所示，N_α(i,j)所处的实际位置为①，该位置所在的物理量有U_αy(i-1,j-1)、E_y(i-1,j-1)。以位置①为参考，位置②、③分别为代入高斯消元法计算更新方程的U_αz(i,j)、E_z(i,j)以及U_αx(i,j)、E_x(i,j)。整个计算过程中，这一点必须满足才能达到稳定。式(13)等式右边的热压力项

需与计算的等离子体速度位置一致。更新U_αx(i,j)、U_αz(i,j)时，涉及x和z方向的偏导，分别需要在②、③号位置离散，即：

化简式(13)左边为：

其中，

同理，化简式(13)右边为：

其中，

由此更新方程可简化为：

即：

式中

与式(14)一致，

与式(15)一致；

采用高斯消元法将式(19)的左手边系数矩阵最终化简为三角矩阵，即：

式中，

等离子体速度的更新式为：

步骤5，基于约束插值曲线法更新计算等离子体密度分布

利用步骤4所得等离子体速度各分量参数，更新计算整个计算区域的离子体密度分布，即在二维空间中逐步在各维度上应用约束插值曲线法，即

具体为：首先计算中间变量

为对应的数密度

在x、z方向的偏导，C(x,Δt)为x方向上应用CIP方法，S(x,Δt)为在x方向上采用中心差分所得结果，根据中间变量得到下一时间点等离子体密度的结果，即：

具体x方向的输运方程为：

化为一般化的运输方程，即：

即f＝N_αx，u＝U_αx，

对上式作偏导，有：

式中，

引入三阶插值：

F_i(x)＝A1_iX³+A2_iX²+g_iX+f_i (30)

X＝x-x_i，本项目中假设下边界x_i＝0，因此三阶插值为：

F_i(x)＝A1_ix³+A2_ix²+g_ix+f_i (31)

其中，三阶、二阶系数分别为：

D≡x_iup＝-Δx·sgn(u_i)，下标iup＝i-sgn(u_i)，u_i表示i位置的各方向速度大小，判断函数sgn(u_i)为：

将x＝0,x＝D代入三阶插值函数F_i(x)，得边界条件：

根据对流特征，可得：

F(x,t+Δt)＝F(x-uΔt,t) (36)

因此通过上式可得：

其中，ξ＝-μ_iΔt；

步骤6，采用时域有限差分离散格式计算各电场分量

采用步骤3所得等离子体各磁场分量参数，步骤4所得等离子体速度参数，以及步骤5所得等离子体密度参数，计算更新各电场分量

具体为：

步骤7，记录采样点的电磁场、等离子体速度、等离子密度数据；

步骤8，将n+1赋值给n，并判断n是否达到预设时间步数，若未达到预设值，则返回步骤3，若达到预设值，则结束。

实施例1，为了验证本发明方法的正确性和有效性，本实施例采用典型的EISCAT电离层等离子体实验条件，予以证明本发明方法的可行性，获得了二维空间下毫秒时间内电磁场与等离子体相互作用过程中，电场、磁场、等离子体密度以及等离子体速度的数值模拟结果。

电磁波频率采用f₀＝4.2MHz，垂直入射进入xoz平面，线极化波E_x(t)＝1.5sin(2πf₀t)。计算区域大小x方向网格数Nx＝50，z方向网格数Nz＝970；空间步长Δx＝5.948m，Δz＝5.948m；时间步长Δt＝7.014638×10^-9s；真空磁导率μ₀＝4π×10^-7H/m；真空介电系数ε₀＝8.85×10^-12F/m；玻尔兹曼常数k_B＝1.38×10^-23J/K；带电荷量e＝1.602×10^-19C；电子质量m_e＝9.1094×10^-31kg，离子质量

碰撞频率ν_e＝500、ν_i＝6；外加磁场B_0x＝sin(12°)×4.8×10^-5T、B_0z＝-cos(12°)×4.8×10^-5T；电子、离子初始密度分布N_e0(z)＝N_crit(1+(z-z_crit)/L_z)，N_crit＝2.1778×10¹¹m^-3，z_crit＝4327m，L_z＝20km；电磁回旋频率ω_ce＝8.4424×10⁶rad/s；入射电磁波的位置Source_z＝51；频率f₀＝4.2MHz；入射波波长λ＝71.3792m；入射电磁波的波形线极化波E_x(t)＝1.5sin(2πf₀t)；电磁波CPML吸收边界，其相关参数PML层的厚度d＝50，σ_max＝0.0149，κ_max＝1，a_max＝0；采样点的坐标(26，486)。

采用本发明方法，循环步数为10000步时，计算结果如图4a—4e所示，由于电场与磁场具有对偶性，仅显示了电场、电子速度、离子速度、电子密度以及离子密度的变化。可以明显看到在高度z_O＝4627m处产生了明显的O波反射现象，在O波反射高度下，电磁波是圆偏振态，包含泵波激发的线偏振态E_x和法拉第磁旋感应导致的E_y，到达O波反射高度后电波由原来的圆偏振态变为沿地磁场方向的线偏振态，Ez幅度的增长是反射回波不断叠加而形成的驻波结构，即电磁波的肿胀效应。同时，可以看到O波反射点处垂向离子速度U_iz和电子速度U_ez对比示意图，离子-电子速度比约为29460，符合离子-电子质量比29376，这也是本发明方法正确的证据之一。对电场Ex进行傅里叶变化后结果如图5所示，可以看到主频发射频率为4.2MHz，该处峰值为电子在泵波驱动下的振荡运动符合理论预期。综上所述，从O波的反射高度、电磁波模式转换(Faraday磁旋效应)、电场频谱分析以及电子-离子速度比，四个方面验证了本发明方法的正确性和有效性。

Claims (1)

1.一种基于MHD的电磁波与磁化等离子体作用数值模拟方法，其特征在于，包括如下步骤：

步骤1，确定MHD模型：

步骤11，建立电磁波与磁化等离子体相互作用的控制方程组，如下式：

其中，下标α表示等离子体的种类：电子e或氧离子O⁺，

分别指x，y，z方向的单位向量，

表示磁场扰动，

表示电场扰动，

是随时间变化的等离子体速度，N_α为等离子体的数密度，式(1)中μ是磁导张量，σ^m为磁损耗张量，式(2)中ε为介电张量，σ^e为电损耗，q_e＝-e，q_O+＝e分别表示带电荷量，电流

由电子运动和离子运动形成，即

式(4)中m_α分别为m_e电子质量，

氧离子质量，B₀表示外加磁场，忽略入射电磁波产生的扰动磁场，ν_α为碰撞频率，k_B是玻尔兹曼常数，T_α是等离子体温度；

步骤12，选取二维平面xOz为外加磁场所在的平面，忽略上述物理参数在y方向梯度的变化

ψ为任意参量，只考虑外加磁场平面上梯度的变化，即垂直方向z方向和x方向的梯度变化，将式(1)至(4)展开，如下式：

步骤2，初始化参数，包括确定计算参量时域和空间域的离散化方案和输入模型文件：

步骤21，确定计算参量时域和空间域的离散化方案：

确定时域递推方案采用蛙跳交替的离散方式，更新循环为：

二维空间域是由Nx×Nz个Yee单元构成的网格，以单元的大小作为分辨率，每个场分量的实际位置与标志具有如下关系：

定义E_x、H_x、U_αx的大小为Nx×Nzp1的零矩阵、E_z、H_z、U_αz的大小为Nxp1×Nz的零矩阵，E_y、H_y、U_αy的大小为Nx×Nz的零矩阵，此处Nzp1＝Nz+1，Nxp1＝Nx+1；

步骤22，输入模型文件：

具体输入的参数包括：计算区域大小x方向网格数Nx，z方向网格数Nz；空间步长Δx，Δz；时间步长Δt；真空磁导率μ₀；真空介电系数ε₀；玻尔兹曼常数k_B；带电荷量e；电子质量m_e，离子质量m_O+；碰撞频率ν_α；外加磁场B_0x、B_0z；电子密度分布N_e0，离子密度分布N_i0；电磁回旋频率ω_ce；入射电磁波的位置Source_x，Source_z；频率f₀；入射波波长λ；确定电磁波的波形，包括余弦函数调制的高斯波形

τ为以参数决定高斯脉冲在时域和频域的宽度，

为时间移动，余弦波形Wave(t)＝cos(2πf(t-t₀))；确定入射电磁波的极化E_{x_wave}＝Wave(t)tan(φ)，E_{y_wave}＝Wave(t)cot(φ)；电磁波CPML吸收边界，其相关参数PML层的厚度d，σ_max，κ_max，a_max取值范围为[0,0.05]；采样点的坐标为sample_x，sample_z)；

步骤3，基于时域有限差分法更新磁场分量

首次循环中，利用步骤2所得电场分量，非首次循环，电场分量系数为步骤6所得，更新计算整个计算区域的磁场分量

具体为：

步骤4，基于高斯消元法更新计算等离子体速度

首次循环，采用步骤2所得的电场分量，等离子体速度分量，等离子体密度参数，非首次循环，采用步骤5所得电子和离子密度分布，步骤6所得电场分量，基于高斯消元法更新整个计算区域等离子体速度分量

具体为：

将速度更新公式简化为：

化简式(13)左边为：

其中，

同理，化简式(13)右边为：

其中，

由此更新方程简化为：

即：

式中

采用高斯消元法将式(19)的左手边系数矩阵最终化简为三角矩阵，即：

式中，

等离子体速度的更新式为：

步骤5，基于约束插值曲线法更新计算等离子体密度分布

利用步骤4所得等离子体速度分量参数，更新计算整个计算区域的离子体密度分布，在二维空间中逐步在各维度上应用约束插值曲线法，

具体为：首先计算中间变量

为对应的数密度

在x、z方向的偏导，C(x,Δt)为x方向上应用CIP方法，S(x,Δt)为在x方向上采用中心差分所得结果，根据中间变量得到下一时间点等离子体密度的结果：

具体x方向的输运方程为：

化为一般化的运输方程：

f＝N_αx，u＝U_αx，

对上式作偏导，有：

式中，

引入三阶插值：

F_i(x)＝A1_iX³+A2_iX²+g_iX+f_i (30)

X＝x-x_i，假设下边界x_i＝0，因此三阶插值为：

F_i(x)＝A1_ix³+A2_ix²+g_ix+f_i (31)

其中，三阶、二阶系数分别为：

D≡x_iup＝-Δx·sgn(u_i)，下标iup＝i-sgn(u_i)，u_i表示i位置的各方向速度大小，判断函数sgn(u_i)为：

将x＝0,x＝D代入三阶插值函数F_i(x)，得边界条件：

根据对流特征，可得：

F(x,t+Δt)＝F(x-uΔt,t) (36)

因此通过上式可得：

其中，ξ＝-μ_iΔt；

步骤6，采用时域有限差分离散格式计算电场分量

采用步骤3所得等离子体磁场分量参数，步骤4所得等离子体速度参数，以及步骤5所得等离子体密度参数，计算更新电场分量

具体为：

步骤7，记录采样点的电磁场、等离子体速度、等离子密度数据；

步骤8，将n+1赋值给n，并判断n是否达到预设时间步数，若未达到预设值，则返回步骤3，若达到预设值，则结束。

Images