CN105808968A - 时域航空电磁数值模拟中c-pml边界条件加载方法 - Google Patents

Abstract

本发明涉及一种时域航空电磁数值模拟中C‑PML边界条件加载方法，通过引入散度方程作为磁场垂直分量的控制方程，推导出C‑PML层内的电场和磁场表达式，确定了C‑PML层分界面上的电磁场无反射条件。利用变限积分方法构建C‑PML层内的时域卷积变步长递归算法，避免了对每一时刻电磁场值的存储，并基于有限差法推导相应的电场和磁场的差分迭代格式，最终实现了时域航空三维晚期高精度、长时窗计算。本发明目的在于可以克服时域航空电磁数值模拟中计算区域的截断误差问题，可以更高效、准确地计算三维时域电磁响应。

Description

时域航空电磁数值模拟中C-PML边界条件加载方法

技术领域

本发明涉及一种地球物理勘探领域的时域电磁响应精准计算方法，尤其是对时域航空电磁响应的精准数值计算。

背景技术

时域航空瞬变电磁法(Time domainAirborne Transient electromagneticmethods)发射系统与接收系统均在空中，相比于时域地面瞬变电磁法(Timedomain Ground Transient electromagnetic methods)，时域航空电磁法能够更加适应复杂多变的地质环境，目前已被大量的用于金属矿、煤炭、油气的勘查，水文地质调查和环境监测等领域。

在时域有限差分(Finite-difference Time-domain)算法中，由于计算机容量限制存在截断边界误差，而吸收边界条件(Absorbing Boundary Condition)的引用大大降低了电磁法正演计算中截断边界处的反射，提高了计算精度。

中国专利CN 102207987A公开了一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，采用OpenCL简化FDTD方法对电磁场仿真的设计与实现，并通过UPML吸收边界条件来提高了计算精度。

中国专利CN 104794289 A公开了一种扩展直角坐标系下完全匹配吸收边界的实现方法，采用WLP-FDTD方法可以求解精细结构模型下的电磁场问题并引入CFS-PML提高计算精度，有效吸收低频与凋落的电磁波。

中国专利CN 104237944A公开了一种适用于交错网格有限差分的全吸收PML方法，通过参数优化，吸收C-PML和M-PML优势，实现了PML在求解TTI介质时的高效吸收，得到与复杂的大斜度井声波测井吻合的数值模拟结果。

中国专利CN 102722651A公开了一种二维柱坐标完全匹配吸收边界的实现方法，把CFS引入柱坐标PML吸收边界以改善PML对低频凋落波与掠射波的吸收效果，很方便的与CFS参数相结合，有效的吸收向外传播的电磁波。

日本专利JP2004239736公开了一种不均匀网格的电磁场仿真方法,应用UPML吸收边界条件来吸收入射波,插值方法求得各时刻虚拟磁场分量来进行电磁场分析。

以上所述方法公布了利用吸收边界条件提高电磁响应计算精度的方法，但国内外专利均未涉及到变步长的C-PML吸收边界条件加载方法及散度条件在C-PML层中的处理办法，对于时域航空三维瞬变电磁计算来说，边界条件的加载是非常重要的，如何在提高计算精度的同时减少计算时间是待本领域技术人员迫切解决的一个技术问题。

发明内容

本发明所要解决的技术问题在于提供一种时域航空电磁数值模拟中C-PML边界条件加载方法，通过引入散度方程作为磁场垂直分量的控制方程，推导出C-PML层内的电场和磁场表达式，确定C-PML层分界面上的电磁场无反射条件。利用变限积分方法构建C-PML层内的时域卷积变步长递推算法，避免了对每一时刻电磁场值的存储，并基于有限差法推导相应的电场和磁场的差分迭代格式。利用GPU并行技术对电磁场的迭代过程进行加速优化，最终实现了时域航空三维晚期高精度、长时窗计算。本发明可以克服三维电磁场数值模拟中计算区域的截断误差问题，适用于时域三维航空电磁响应有限差分的快速、精准计算。

本发明是这样实现的，一种时域航空电磁数值模拟中C-PML边界条件加载方法：

1)、采用准静态条件下，散度方程作为磁场垂直分量H_z的控制方程，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程；

2)、根据C-PML层中散度和旋度运算，推导电、磁场控制方程的表达式；

3)、推导具有散度方程的C-PML层和计算区域分界面电磁波无反射的条件，设置C-PML层参数；

4)、利用变限积分方法，构建C-PML层内的时域卷积项ψ变步长递归算法；

5)、基于有限差分法，结合卷积项递归算法，推导控制方程在整个计算区域的差分迭代格式；

6)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率等参数，实现含有C-PML边界条件的电、磁场迭代运算；

7)、迭代计算结束后，提取磁场各分量响应，并对计算结果进行显示。

进一步地，步骤2中，采用散度方程作为磁场垂直分量的控制方程，以保证解的稳定性与唯一性。其在C-PML层表达形式为：

$\frac{1}{s_{mx}}\frac{\∂H_x}{\∂x}+\frac{1}{s_{my}}\frac{\∂H_y}{\∂y}+\frac{1}{s_{mz}}\frac{\∂H_z}{\∂z}=0---\left(1\right)$

式(1)中S_mi(i＝x、y、z)为各个方向C-PML层中的磁场的延伸坐标矩阵式，其表达式为k_mi、α_mi(i＝x、y、z)为C-PML层引入的参数，σ_pmi为延伸坐标中磁导率。

经过时-频转换后，式(1)可以表达成：

$\frac{1}{k_{mx}}\frac{\∂H_x}{\∂x}+\frac{1}{k_{my}}\frac{\∂H_y}{\∂y}+\frac{1}{k_{mz}}\frac{\∂H_z}{\∂z}+{\mathrm{\ψ}}_{hzx}+{\mathrm{\ψ}}_{hzy}+{\mathrm{\ψ}}_{hzz}=0---\left(2\right)$

式(2)中ψ_hzi(i＝x、y、z)为求解H_z时与对各个方向偏导数有关的卷积项。

${\mathrm{\ψ}}_{hzi}=\left(\frac{\mathrm{\δ}\left(t\right)}{k_{mi}}-\frac{{\mathrm{\σ}}_{pmi}}{{{\mathrm{\ϵk}}_{mi}}^2}\exp \left(-\left(\frac{{\mathrm{\σ}}_{pmi}}{{\mathrm{\ϵk}}_{mi}}+\frac{{\mathrm{\α}}_{pmi}}{\mathrm{\ϵ}}\right)t\right)\mathrm{\μ}\left(t\right)\right)*\frac{\∂H_i}{\∂i}---\left(3\right)$

式(3)中ε为人工介电常数，δ(t)为单位脉冲函数，μ(t)为单位阶跃函数。

进一步地，步骤3中，电磁波在C-PML层无反射的原理是C-PML层中波阻抗与问题空间波阻抗相等，以任意方向扩散的二维TMx极化波为例：

平面电磁波解的形式为：

E_x＝E₀e^{jω(t-ay-βz)}(4)

代入Maxwell方程推导出C-PML层波阻抗Z₁的表达式为：

问题空间波阻抗表达式为：

$Z_2=\sqrt{\frac{\mathrm{\μ}j\mathrm{\ω}}{j\mathrm{\ω}\mathrm{\ϵ}+\mathrm{\σ}}}---\left(6\right)$

令在入射角为任意值情况下Z₁＝Z₂，得出平面波在C-PML层无反射的条件为：

s_ei＝s_mi(7)

而从而得出平面波在C-PML层无反射的条件为：

$k_{ei}=k_{mi},\frac{{\mathrm{\σ}}_{pei}}{{\mathrm{\σ}}_{pmi}}=\frac{{\mathrm{\α}}_{ei}}{{\mathrm{\α}}_{mi}}=\frac{\mathrm{\ϵ}}{\mathrm{\μ}}---\left(8\right)$

式(7)中S_ei(i＝x、y、z)为各个方向C-PML层中的电场的延伸坐标矩阵，S_mi(i＝x、y、z)为各个方向C-PML层中的磁场的延伸坐标矩阵式，式(8)中ε为人工介电常数，μ为问题空间磁导率，σ_pei与σ_pmi为C-PML层中的电导率与磁导率，k_ei、k_mi、α_ei、α_mi(i＝x、y、z)为C-PML层中引入的参数。

进一步地，步骤4中，在迭代时间步长Δt递增情况下，从卷积项的时频转换关系入手，利用变限积分积分的方法，推导出C-PML层中卷积项ψ的递推关系式为(以ψ_exy为例)：

$\begin{array}{c}{\mathrm{\ψ}}_{exy}\left(n+1\right)={\mathrm{\ξ}}_{ey}\left(t\right)*\frac{\∂H_z}{\∂y}={\∫}_0^{t_{n+1}}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)\frac{\∂H_z\left(t_{n+1}-\mathrm{\τ}\right)}{\∂y}d\mathrm{\τ}\\ =b_{ey}\left({\mathrm{\Δt}}_n\right){\mathrm{\ψ}}_{exy}\left(n\right)+b_{ey}\left({\mathrm{\Δt}}_n\right){\∫}_{t_n}^{t_{n+1}}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)\frac{\∂H_z\left(t_n-\mathrm{\τ}\right)}{\∂y}d\mathrm{\τ}+{\∫}_{t_n}^{t_{n+1}}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)\frac{\∂H_z\left(t_{n+1}-\mathrm{\τ}\right)}{\∂y}d\mathrm{\τ}\\ -b_{ey}\left({\mathrm{\Δt}}_n\right)\frac{1}{\mathrm{\Δ}y}{\∫}_0^{{\mathrm{\Δt}}_n}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)d\mathrm{\τ}\left(H_z^{n-1+\frac{1}{2}}\left(i,j,k\right)-H_z^{n-1+\frac{1}{2}}\left(i,j-1,k\right)\right)\end{array}---\left(9\right)$

用τ-Δt_n取代τ代入等式(9)右边第二项，并将第二项与第三项整合得出：

$\begin{array}{c}{\mathrm{\ψ}}_{exy}(n+1)=b_{ey}\left({\mathrm{\Δt}}_n\right){\mathrm{\ψ}}_{exy}\left(n\right)+{\∫}_{t_n-{\mathrm{\Δt}}_n}^{t_{n+1}}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)\frac{\∂H_z\left(t_{n+1}-\mathrm{\τ}\right)}{\∂y}d\mathrm{\τ}\\ -b_{ey}\left({\mathrm{\Δt}}_n\right)\frac{1}{\mathrm{\Δ}y}{\∫}_0^{{\mathrm{\Δt}}_n}{\mathrm{\ξ}}_{ey}\left(\mathrm{\τ}\right)d\mathrm{\τ}\left(H_z^{n-1+\frac{1}{2}}\left(i,j,k\right)-H_z^{n-1+\frac{1}{2}}\left(i,j-1,k\right)\right)\end{array}---\left(10\right)$

用取代τ代入等式(10)右边第二项，并整理得出：

$\begin{array}{c}{\mathrm{\ψ}}_{exy}\left(n+1\right)=b_{ey}\left({\mathrm{\Δt}}_n\right){\mathrm{\ψ}}_{\mathrm{ex}y}\left(n\right)-b_{ey}\left({\mathrm{\Δt}}_n\right)a_{ey}\left({\mathrm{\Δt}}_n\right)\left(H_z^{n-1+\frac{1}{2}}\left(i,j,k\right)-H_z^{n-1+\frac{1}{2}}\left(i,j-1,k\right)\right)\\ +\left(a_{ey}\left(t_{n+1}\right)-a_{ey}\left(t_n\right)\right)\left(H_z^{n-1+\frac{1}{2}}\left(i,j,k\right)-H_z^{n-1+\frac{1}{2}}\left(i,j-1,k\right)\right)\\ +\left(a_{ey}\left(t_n\right)-a_{ey}\left(t_n-{\mathrm{\Δt}}_n\right)\right)\left(H_z^{n+\frac{1}{2}}\left(i,j,k\right)-H_z^{n+\frac{1}{2}}\left(i,j-1,k\right)\right)\end{array}---\left(11\right)$

式(11)中，ψ_exy(n)表示在t_n时刻，求解E_x时与对y方向偏导有关的卷积项。表示在时刻时刻，在网格节点i、j、k处的z方向磁场值,Δt_n＝t_n+1-t_n。b_ey(x)、a_ey(x)为涉及到的中间变量其表达式为:

$a_{ey}\left(x\right)=\frac{{\mathrm{\σ}}_{pey}}{\mathrm{\Δ}y({\mathrm{\σ}}_{pey}{k}_{ey}+{\mathrm{\α}}_{pey}{k}_{ey}^2)}\left(b_{ey}\right(x)-1)---\left(12\right)$

$b_{ey}\left(x\right)=\exp (-(\frac{{\mathrm{\σ}}_{pey}}{{\mathrm{\ϵk}}_{ey}}+\frac{{\mathrm{\α}}_{pey}}{\mathrm{\ϵ}}\left)x\right)---\left(13\right)$

进一步地，步骤5中，有限差分方法离散磁场垂直分量H_z的迭代方程时，由于其控制方程中的卷积项ψ_hzz与当前时刻下一网格中的H_z有关，需将ψ_hzz的递推式代入控制方程重新整理H_z的迭代格式，使其只与前时刻H_z值有关:

$\begin{array}{c}{H}_z^{n+\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)=H_z^{n+\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2},k+1\right)+B(H_z^{n-1+\frac{1}{2}}\left(i+\frac{1}{2},j+\frac{1}{2},k+1\right)-H_z^{n-1+\frac{1}{2}}(i+\frac{1}{2},j+\\ \frac{1}{2},k\left)\right)+\frac{A}{k_{my}\mathrm{\Δ}y}\left(H_y^{n+\frac{1}{2}}\left(i+\frac{1}{2},j+1,k+\frac{1}{2}\right)-H_y^{n+\frac{1}{2}}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)\right)+\frac{A}{k_{mx}\mathrm{\Δ}x}\left(H_x^{n+\frac{1}{2}}\right(i+1,j+\frac{1}{2},k+\frac{1}{2})-\\ H_x^{n+\frac{1}{2}}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right))+{\mathrm{A\ψ}}_{hzx}^{n+\frac{1}{2}}+{\mathrm{A\ψ}}_{hzy}^{n+\frac{1}{2}}+{\mathrm{Ab}}_{mz}({\mathrm{\Δt}}_n){\mathrm{\ψ}}_{hzz}^{n-\frac{1}{2}}\end{array}---\left(14\right)$

其中：

$A={\left(a_{mz}\left(t_{n-1+\frac{1}{2}}\right)-a_{mz}\left(t_{n-1+\frac{1}{2}}-{\mathrm{\Δt}}_{n-1+\frac{1}{2}}\right)+\frac{1}{k_{mz}\mathrm{\Δ}z}\right)}^{-1}---\left(15\right)$

$B=A\left(a_{mz}\right(t_{n+\frac{1}{2}})-a_{mz}(t_{n-1+\frac{1}{2}})-b_{mz}({\mathrm{\Δt}}_{n-1+\frac{1}{2}}\left)a_{mz}\right({\mathrm{\Δt}}_{n-1+\frac{1}{2}}\left)\right)---\left(16\right)$

${\mathrm{\Δt}}_{n-1+\frac{1}{2}}=t_{n+\frac{1}{2}}-t_{n-1+\frac{1}{2}}---\left(17\right)$

$a_{mz}\left(x\right)=\frac{{\mathrm{\σ}}_{pmz}}{\mathrm{\Δ}y({\mathrm{\σ}}_{pmz}{k}_{mz}+{\mathrm{\α}}_{pmz}{k}_{mz}^2)}\left(b_{mz}\right(x)-1)---\left(18\right)$

$b_{mz}\left(x\right)=\exp (-(\frac{{\mathrm{\σ}}_{pmz}}{{\mathrm{\ϵk}}_{mz}}+\frac{{\mathrm{\α}}_{pmz}}{\mathrm{\ϵ}}\left)x\right)---\left(19\right)$

本发明与现有技术相比，有益效果在于：针对时域航空电磁数值模拟计算中引入散度条件作为控制方程，以及迭代时间步长非均匀无法引用常规卷积递推公式等问题，本发明推导出C-PML层内基于散度方程的磁场表达式及其差分迭代格式，确定C-PML层分界面上的电磁场无反射条件，利用变限积分方法构建C-PML层内的时域卷积变步长递推算法，避免了对每一时刻电磁场值的存储，利用GPU并行技术对电磁场的迭代过程进行加速优化，达到了对时域航空三维电磁响应进行快速、精准计算的目的。

附图说明

图1是时域航空电磁数值模拟中C-PML边界条件加载方法示意图；

图2是时域航空电磁响应计算模型示意图；

图3是时域航空感应电动势的不同边界条件吸收效果对比图；

图4是加载C-PML边界条件后感应电动势的计算相对误差；

图5是均匀半空间模型加载C-PML边界条件吸收效果切片图；

图6是均匀半空间模型加载狄利克雷边界条件吸收效果切片图。

具体实施方式

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

实施例

参见图1，一种时域航空电磁数值模拟中C-PML边界条件加载方法，包括：

1)、采用散度方程作为磁场垂直分量的控制方程，准静态条件下、无源Maxwell旋度方程作为其他电磁场分量的控制方程结合非均匀Yee式网格和DuFort-Frankel方法构建显式差分格式，网格数目为107×107×53，其中x、y方向上的网格数目均为107个，z方向上网格数53个，最小和最大网格步长分别为10m和120m，除地面外计算区域其余5面最外层8个网格加载C-PML层。

2)在整个计算区域内设置电导率、磁导率、C-PML系数等参数，算例中模型为均匀半空间模型，电导率设置为0.01西门子/米，磁导率设置为真空磁导率。

3)、将计算区域的电性参数代入到发射在空中的地下电场响应表达式中，计算初始时刻的x、y方向的电场响应。根据Maxwell旋度方程计算下半个时刻磁场响应。

4)、将参与迭代的CPU序列转化为GPU序列。

5)、计算当前时刻迭代时间步长。

${\mathrm{\Δt}}_n=\mathrm{\α}\sqrt{\frac{{\mathrm{\μ\σt}}_n}{6}}{\mathrm{\Δ}}_{min}---\left(1\right)$

Δ_min为最小网格步长，α取值0.1。

6)、将电场初始场值代入迭代方程中在整个计算区域更新电场值E_x、E_y、E_z，保存前一时刻与当前时刻电场值。计算各分量电场卷积项ψ_e，将卷积项加到对应区域的电场值中。

7)、将磁场初始场值代入迭代方程中在整个计算区域更新磁场值H_x、H_y、H_z，保存前一时刻与当前时刻磁场值。计算各分量磁场卷积项ψ_m，将卷积项加到对应区域的磁场值中。

8)、向上延拓方法将地面电磁场延拓到空中，求得空中接收线圈处电磁场值。

9)、是否完成全部迭代，如未完成，重复5-8步骤，如果完成全部迭代，将计算结果从GPU释放出来，并对计算结果进行显示，完成全部计算。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种时域航空电磁数值模拟中C-PML边界条件加载方法其特征在于，包括如下步骤：

1)、采用准静态条件下，散度方程作为磁场垂直分量H_z的控制方程，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程；

2)、根据C-PML层中散度和旋度运算，推导电、磁场控制方程的表达式；

3)、推导具有散度方程的C-PML层和计算区域分界面电磁波无反射的条件，设置C-PML层参数；

4)、利用变限积分方法，构建C-PML层内的时域卷积项ψ变步长递归算法；

5)、基于有限差分法，结合卷积项递归算法，推导控制方程在整个计算区域的差分迭代格式；

6)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率等参数，实现含有C-PML边界条件的电、磁场迭代运算；

7)、迭代计算结束后，提取磁场各分量响应，并对计算结果进行显示。

2.按照权利要求1所述的一种时域航空电磁数值模拟中C-PML边界条件加载方法，其特征在于：

步骤2中，采用散度方程作为磁场垂直分量的控制方程，以保证解的稳定性与唯一性，散度方程在C-PML层表达形式为：

$\frac{1}{S_{mx}}\frac{\∂H_x}{\∂x}+\frac{1}{S_{my}}\frac{\∂H_y}{\∂y}+\frac{1}{S_{mz}}\frac{\∂H_z}{\∂z}=0---\left(1\right)$

式(1)中S_mi(i＝x、y、z)为各个方向C-PML层中的磁场的延伸坐标矩阵式，其表达式为k_mi、α_mi(i＝x、y、z)为C-PML层引入的参数，σ_pmi为C-PML层中的磁导率；

经过时频转换后，式(1)可以表达成：

$\frac{1}{k_{mx}}\frac{\∂H_x}{\∂x}+\frac{1}{k_{my}}\frac{\∂H_y}{\∂y}+\frac{1}{k_{mz}}\frac{\∂H_z}{\∂z}+{\mathrm{\ψ}}_{hzx}+{\mathrm{\ψ}}_{hzy}+{\mathrm{\ψ}}_{hzz}=0---\left(2\right)$

式(2)中ψ_hzi(i＝x、y、z)为求解H_z时，对应各个方向偏导数有关的卷积项。

3.按照权利要求1所述的一种时域航空电磁数值模拟中C-PML边界条件加载方法，其特征在于：

步骤3中，根据电磁波在C-PML层中无反射条件，为C-PML层中的波阻抗与计算区域的波阻抗相等，以任意方向扩散的二维TMx极化波为例，将平面电磁波的解形式代入Maxwell方程中，推导出C-PML层波阻抗Z₁的表达式，并令其在入射角为任意值情况下与计算区域波阻抗相等，最终推出平面波在C-PML层无反射的条件为：

s_ei＝s_mi (3)

即:

式(3)中S_ei(i＝x、y、z)为各个方向C-PML层中的电场延伸坐标矩阵，其表达式为ε为人工介电常数，σ_pei与σ_pmi为C-PML层中的电导率与磁导率，μ为计算区域磁导率，k_ei、k_mi、α_ei、α_mi(i＝x、y、z)为C-PML层中引入的参数。

4.按照权利要求1所述的一种时域航空电磁数值模拟中C-PML边界条件加载方法，其特征在于：

步骤4中，在迭代时间步长Δt递增情况下，从卷积项的积分表达式入手，利用变限积分的方法，推导出C-PML层中卷积项ψ的递推关系式为(以ψ_exy为例)：

$\begin{array}{c}{\mathrm{\ψ}}_{exy}(n+1)=b_{ey}\left({\mathrm{\Δt}}_n\right){\mathrm{\ψ}}_{exy}\left(n\right)-b_{ey}\left({\mathrm{\Δt}}_n\right)a_{ey}\left({\mathrm{\Δt}}_n\right)\left(H_z^{n-1+\frac{1}{2}}\right(i,j,k)-H_z^{n-1+\frac{1}{2}}(i,j-1,k\left)\right)\\ +\left(a_{ey}\right(t_{n+1})-a_{ey}(t_n\left)\right)\left(H_z^{n-1+\frac{1}{2}}\right(i,j,k)-H_z^{n-1+\frac{1}{2}}(i,j-1,k\left)\right)\\ +\left(a_{ey}\right(t_n)-a_{ey}(t_n-{\mathrm{\Δt}}_n\left)\right)\left(H_z^{n+\frac{1}{2}}\right(i,j,k)-H_z^{n+\frac{1}{2}}(i,j-1,k\left)\right)\end{array}---\left(5\right)$

式(5)中，ψ_exy(n)表示在t_n时刻，求解E_x时与对y方向偏导有关的卷积项，表示在时刻时刻，在网格节点i、j、k处的z方向磁场值，Δt_n＝t_n+1-t_n，b_ey(x)、a_ey(x)为涉及到的中间变量其表达式为:

$\begin{array}{c}{a}_{ey}\left(x\right)=\frac{{\mathrm{\σ}}_{pey}}{\mathrm{\Δ}y({\mathrm{\σ}}_{pey}{k}_{ey}+{\mathrm{\α}}_{pey}{k}_{ey}^2)}\left(b_{ey}\right(x)-1)---\left(6\right)\\ b_{ey}\left(x\right)=\exp (-(\frac{{\mathrm{\σ}}_{pey}}{{\mathrm{\ϵk}}_{ey}}+\frac{{\mathrm{\α}}_{pey}}{\mathrm{\ϵ}}\left)x\right)---\left(7\right)\end{array}.$

5.按照权利要求1所述的一种时域航空电磁数值模拟中C-PML边界条件加载方法，其特征在于：

步骤5中，有限差分方法离散磁场垂直分量H_z的迭代方程时，由于其控制方程(2)中的卷积项ψ_hzz与当前时刻下一网格中的H_z有关，需将ψ_hzz的递推式代入控制方程重新整理H_z的迭代格式，使其只与前时刻H_z值有关:

$\begin{array}{c}{H}_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2},k)=H_z^{n+\frac{1}{2}}(i+\frac{1}{2},j+\frac{1}{2},k+1)+B\left(H_z^{n-1+\frac{1}{2}}\right(i+\frac{1}{2},j+\frac{1}{2},k+1)-H_z^{n-1+\frac{1}{2}}(i+\frac{1}{2},j+\\ \frac{1}{2},k\left)\right)+\frac{A}{k_{my}\mathrm{\Δ}y}\left(H_y^{n+\frac{1}{2}}\right(i+\frac{1}{2},j+1,k+\frac{1}{2})-H_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k+\frac{1}{2}\left)\right)+\frac{A}{k_{mx}\mathrm{\Δ}x}\left(H_x^{n+\frac{1}{2}}\right(i+1,j+\frac{1}{2},k+\frac{1}{2})-\\ H_x^{n+\frac{1}{2}}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right))+{\mathrm{A\ψ}}_{hzx}^{n+\frac{1}{2}}+{\mathrm{A\ψ}}_{hzy}^{n+\frac{1}{2}}+{\mathrm{Ab}}_{mz}\left({\mathrm{\Δt}}_n\right){\mathrm{\ψ}}_{hzz}^{n-\frac{1}{2}}\end{array}---\left(8\right)$

其中:

$A={\left(a_{mz}\right(t_{n-1+\frac{1}{2}})-a_{mz}(t_{n-1+\frac{1}{2}}-{\mathrm{\Δt}}_{n-1+\frac{1}{2}})+\frac{1}{k_{mz}\mathrm{\Δ}z})}^{-1}---\left(9\right)$

$B=A\left(a_{mz}\right(t_{n+\frac{1}{2}})-a_{mz}(t_{n-1+\frac{1}{2}})-b_{mz}({\mathrm{\Δt}}_{n-1+\frac{1}{2}}\left)a_{mz}\right({\mathrm{\Δt}}_{n-1+\frac{1}{2}}\left)\right)---\left(10\right)$

${\mathrm{\Δt}}_{n-1+\frac{1}{2}}=t_{n+\frac{1}{2}}-t_{n-1+\frac{1}{2}}---\left(11\right)$

$a_{mz}\left(x\right)=\frac{{\mathrm{\σ}}_{pmz}}{\mathrm{\Δ}y({\mathrm{\σ}}_{pmz}{k}_{mz}+{\mathrm{\α}}_{pmz}{k}_{mz}^2)}\left(b_{mz}\right(x)-1)---\left(12\right)$

$b_{mz}\left(x\right)=\exp (-(\frac{{\mathrm{\σ}}_{pmz}}{{\mathrm{\ϵk}}_{mz}}+\frac{{\mathrm{\α}}_{pmz}}{\mathrm{\ϵ}}\left)x\right)---\left(13\right)$

公式(8)到(10)涉及到的变量前面均有提及。