[发明专利]一种基于MHD的电磁波与磁化等离子体作用数值模拟方法

权利要求书

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，若达到预设值，则结束。