CN102332055A - 一种极低频电磁波的仿真计算方法 - Google Patents

Abstract

本发明公开了一种极低频电磁波的仿真计算方法，属于极低频电磁波计算方法技术领域。本发明所述方法包括亚网格划分步骤、电磁场的布置步骤、计算步骤；如下的内容：亚网格划分步骤针对多种地形学网格模型的不同的亚网格划分方案；电磁场的布置步骤是对于基于三角形地形学网格、基于六边形地形学网格的电磁场分布结合到一起，形成基于地形学网格划分的亚网格元胞；计算步骤基于不同亚网格划分方案引入的内插或近似算法。本发明不仅提高了全球极低频电磁波算法的计算精度和有效分辨率，同时也大大提高了算法的计算效率。

Description

一种极低频电磁波的仿真计算方法

技术领域

本发明涉及一种基于地形学时域有限差分算法的亚网格技术，尤其涉及一种提高极低频电磁波算法的计算效率、精度和局部分辨率的方法。

背景技术

众所周知，极低频电磁波(频率低于3000Hz)的仿真计算需要计算地球-电离层系统中的电磁波传播，建模所需的计算机资源庞大，而日益严格的仿真要求又需要较高的建模分辨率，所以针对当前主流的极低频时域算法提出一种既能保证其计算效率，又能有效的提高其局部分辨率的技术是非常必要的。

近年来，针对极低频传播的问题的应用研究逐渐受到重视[1]。对于这个频段的电磁波来说，传统的电磁波传播模型(解析方法)给出的ELF电磁波传播的粗略规律和真实的实验结果存在较大的误差，已经无法适应当前科学研究要求[2]。新型的数值方法(这里主要是指新型的基于时域有限差分方法[3]的数值解法)针对解析方法的局限性，引入较为复杂的地质和电磁环境，很大程度上解决了复杂的ELF电磁波传播等问题[1、3]。但是由于计算手段的限制，特别是当前个人计算机和工作站的计算能力的限制，局部高精度和高分辨要求的复杂电磁环境的ELF频段电磁波传播问题往往得不到有效的解决。这是因为对于目前的最有效的极低频电磁算法，由于其多数都是基于传统时域有限差分方法的变形[1]，所以整个模型的分辨率统一。如果要提高全球模型中的局部研究区域的精度，则整个电磁模型的分辨率都要提高，这大大影响了算法的计算效率。

比如在地震预测的研究过程中，由于地震电磁效应影响范围大，所以整个地球-电离层系统都需要进行仿真建模。与此同时，震区的电磁环境非常复杂，每一个地质因素都与即将发生的全球电磁现象有直接的关系，所以对震区的网格精度有效高的要求。这样，在整个地球电离层模型的仿真的基础上，由于网格的高精度要求，使得网格数量非常大，仿真时间和资源消耗巨大。此时，采用局部亚网格技术，在震区内采用高精度的亚网格建模方案，在震区外采用精度较低的粗网格，这样可以有效地大量的减少计算消耗，提高仿真效率。

参考文献：

[1]J.Simpson.Current and future applications of 3-D global Earth-ionosphere models basedon the full-vector Maxwell’s equations FDTD method.Surveys in Geophysics，2009，30：105-130.

[2]L.Sevgi，F.Akleman，and L.B.Felsen，“Groundwave propagation modeling：problem-matched analytical formulations and direct numerical techniques，”IEEE Antennas andPropagation Magazine，vol.44，no.1，pp.55-75，2002.

[3]A.Taflove and S.C.Hagness，Computational electrodynamics：the finite-differencetime-domain method，3rd ed.Artech House，2005.

发明内容

本发明所要解决的技术问题是提高当前主流的用于求解极低频电磁波传播问题的计算效率、精度和局部分辨率，为此而提出了基于地形学时域有限差分算法的亚网格技术。

本发明为解决上述技术问题采用以下技术方案：

一种极低频电磁波的仿真计算方法，包括亚网格划分步骤、电磁场的布置步骤、计算步骤；其中，

A、所述亚网格划分步骤如下：

A-1，基于三角形地形学网格划分：将原始主球面三角形网格的各个边的中点沿球面弧两两相连接，形成一种对主三角形网格的球面划分，划分出4个球面三角形作为亚网格元胞；

A-2，基于六边形地形学网格划分：将原始主球面六边形网格中相邻三个主六边形互相交接的三个边去掉，以三个面积大小接近为原始六边形1/4的三个新的六边形，划分出3个球面六边形作为亚网格元胞；

A-3，将所述步骤A-1划分的三角形亚网格及步骤A-2划分的六边形亚网格进行叠加，在径向上将三角形亚网格和六边形亚网格相结合，形成三维的亚网格元胞，成为对整个地球-电离层系统的完整亚网格划分；

B、所述电磁场的布置步骤如下：

B-1，对于基于三角形地形学网格划分，将被亚网格化的主网格上的1个径向场和3个切向场去除，以新的亚网格上的3个径向场和9个切向场分别布置电磁场；

B-2，对于基于六边形地形学网格划分，将被亚网格化的3个主网格上的相邻3条边上的3个切向场去除，以新的3个亚网格上的3个径向场和9个切向场、6个虚拟切向场分别布置电磁场；

B-3，将步骤B-1、步骤B-2的电磁场分布结合到一起，形成基于地形学网格划分的亚网格元胞；

C、计算步骤如下：

C-1，对于基于三角形地形学网格划分，在原始的主网格和新的亚网格内的电磁场计算沿用传统的经典时域有限差分算法，在原始的主网格和新的亚网格边界处的电磁场计算采用如下方法：

C-1-1，对于边界处主网格的径向磁场有：

${\mathrm{Hr}}^{n+0.5}={\mathrm{Hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{norm}}}\×\{\mathrm{\Σ}({\mathrm{Et}}^n\×\mathrm{\Δ}{l}_{\mathrm{norm}})+\mathrm{\Σ}({\mathrm{et}}^n\×\mathrm{\Δ}{l}_{\mathrm{sub}})\}---\left(1\right)$

其中Hr为主网格径向磁场，Et为主网格切向电场，et为亚网格上的切向电场，Δt为时间步长，ΔS_norm、Δl_norm、Δl_sub，分别为主网格面积、主网格边长、亚网格边长，μ为磁导率，n为时间步数，取值由所需运行时间决定；

C-1-2，对于边界处的亚网格切向电场有：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{hr}}^{n+0.5}mH{r}^{n+0.5})---\left(2\right)$

其中，hr为亚网格上的径向磁场，δ_sub为主网格中心到亚网格中心的距离，ε为介电系数；

C-2，对于基于六边形地形学网格划分，在原始的主网格和新的亚网格内的电磁场的计算沿用传统的经典时域有限差分算法，在原始的主网格和新的亚网格边界处的电磁场采用以下方法：

C-2-1，对于边界处主网格的径向电场有：

${\mathrm{Er}}^{n+1}={\mathrm{Er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{norm}}}\×\left\{\mathrm{\Σ}{\mathrm{Ht}}^{n+0.5}\×l_{\mathrm{norm}}\±\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\×0.5\×l_{\mathrm{sub}}\right\}---\left(3\right)$

其中Er为边界上主网格径向电场，Ht为主网格切向磁场，ht为亚网格上的切向磁场；

C-2-2，对于边界处的亚网格切向电场有：

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{er}}^{n}mE{r}^n)---\left(4\right)$

其中er为亚网格上的径向电场；

C-2-3，对于亚网格和主网格边界上的虚拟切向磁场有两种计算方法：

①，同向场近似法：

htv0≈Ht₀                   (5)

②，周围场插值法：

$\mathrm{htv}0\≈-\frac{1}{2}\mathrm{ht}0+\frac{1}{6}({\mathrm{Ht}}_0+{\mathrm{Ht}}_1+\mathrm{htv}1)$ (6)

$\≈\frac{7}{35}{\mathrm{Ht}}_0+\frac{5}{35}{\mathrm{Ht}}_1+\frac{3}{35}\mathrm{ht}1-\frac{18}{35}\mathrm{ht}0$

其中，ht0，ht1为亚网格切向磁场，Ht₀，Ht₁为主网格切向磁场，htv0，htv1为亚网格虚拟磁场；

将步骤C-1及步骤C-2所述的两种平面二维算法叠加，得到三维亚网格算法，计算切向电场分量、径向电场分量、虚拟切向磁场分量，其中径向电场分量和虚拟切向磁场的计算方法与C-1，C-2所示步骤类同。

进一步的，本发明的一种极低频电磁波的仿真计算方法，步骤C-3所述三维亚网络算法的具体计算步骤为：

(1)计算主六边形网格中心的径向电场Er，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Er}}^{n+1}={\mathrm{Er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Ht}}^{n+0.5}\×l_{\mathrm{norm}})---\left(7\right)$

(2)计算亚六边形网格中心的径向电场er，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{er}}^{n+1}={\mathrm{er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{ht}}^{n+0.5}\×l_{\mathrm{sub}})---\left(8\right)$

(3)采用公式(3)计算与亚网格相邻的主网格上的径向电场；

(4)计算主三角形网格上的切向电场Et，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Et}}^{n+1}={\mathrm{Et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Ht}}^{n+0.5}\right)---\left(9\right)$

其中δ_norm、Δr分别为主网格中心距离、径向空间步长。

(5)计算亚网格上的切向电场et：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(10\right)$

(6)利用公式(2)的三维形式来计算边界处的亚网格切向电场：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{hr}}^{n+0.5}\mathrm{mH}{r}^{n+0.5})+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(11\right)$

至此，电场计算结束，下面计算磁场部分：

(7)计算主三角形网格中心(不包括与亚网格相邻)的径向磁场Hr，有

${\mathrm{Hr}}^{n+0.5}={\mathrm{Hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Et}}^n\×\mathrm{\Δ}{L}_{\mathrm{norm}})---\left(12\right)$

(8)计算亚三角形网格中心的径向磁场hr，有

${\mathrm{hr}}^{n+0.5}={\mathrm{hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{et}}^n\×\mathrm{\Δ}{l}_{\mathrm{sub}})---\left(13\right)$

(9)计算主六边形网格边上的切向磁场Ht，有

${\mathrm{Ht}}^{n+0.5}={\mathrm{Ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Et}}^n\right)---\left(14\right)$

(10)计算亚六边形网格边上的切向磁场ht，有

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(15\right)$

(11)利用公式(4)的三维形式来计算边界处的亚网格切向磁场：

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{er}}^{n}m{\mathrm{Er}}^n)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(16\right)$

(12)利用公式(5)或者公式(6)计算虚拟切向磁场htv；

至此，全部的迭代过程完毕。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明借鉴传统时域有限差分方法中的亚网格思想，提出基于地形学的球面非均匀亚网格技术。采用这种技术，可以在不影响全球模型分辨率的前提条件下提高局部分辨率，做到“粗中有细”，从而显著提高算法在实际应用研究中的执行效率，使局部高分辨率高精度要求得到满足。同时，本技术在提高局部分辨率的同时有效的节省了计算消耗，减少了计算时间。本发明不仅提高了全球极低频电磁波算法的计算精度和有效分辨率，同时也大大提高了算法的计算效率。

附图说明

图1是亚网格划分方案示意图；

其中，图1(a)：三角形亚网格划分方案；图1(b)：六边形亚网格划分方案。

图2是六边形亚网格划分方案中水平(切向)虚拟场计算方案。

图3是地形学网格划分示意图(二维)；

其中图3(a)：三角形亚网格划分方案；图3(b)：六边形亚网格划分方案。

图4是地形学网格划分示意图(三维)。

图5是亚网格划分方案在整个地形学划分方案中的示意图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

1)准备工作。在基于地形学的球面已经完成的基础上，对于具有高分辨率和高精度要求的区域采用上面给出的划分方案进行亚网格再划分。对于图3(a)所示的三角形球面划分，如图1(a)所示，采用将一个三角形网格划分为4个小的球面三角形网格(亚网格)的方法，然后记录每个亚三角形的边长、面积。对于图3(b)所示的六边形球面划分，采用将原始三个将交接的主六边形的三个边去掉，取而代之以三个面积大小接近为原始六边形1/4的三个新的六边形(亚网格)的方法(如图1(b)所示)，然后记录每个亚六边形的边长、面积、与相邻主六边形和亚六边形的中心距离。对于图4所示的三维情况，亚三角形与亚六边形需要处于相同的位置(如图5所示)，并同时记录所有亚网格的信息。针对地形学时域有限差分算法模型的特点，在地球径向，不需要进行亚网格划分。

2)电磁场的布置步骤如下：

B-1，对于基于三角形地形学网格划分，将被亚网格化的主网格上的1个径向场和3个切向场去除，以新的亚网格上的3个径向场和9个切向场分别布置电磁场；

B-2，对于基于六边形地形学网格划分，将被亚网格化的3个主网格上的相邻3条边上的3个切向场去除，以新的3个亚网格上的3个径向场和9个切向场、6个虚拟切向场分别布置电磁场；

B-3，将步骤B-1、步骤B-2的电磁场分布结合到一起，形成基于地形学网格划分的亚网格元胞；

3)基于三维地形学时域有限差分法的亚网格具体计算步骤。在每一个时间步内，按照如下的计算方案执行：

(1)计算主六边形网格(不包括与亚网格相邻的主网格)中心的径向电场Er，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Er}}^{n+1}={\mathrm{Er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Ht}}^{n+0.5}\×l_{\mathrm{norm}})---\left(7\right)$

其中Er为主网格径向电场，Ht为主网格切向磁场，Δt为时间步长，ΔS_norm、ΔL_norm分别为主网格面积、主网格边长，ε为介电系数。下同，不再复述。

(2)计算亚六边形网格中心的径向电场er(如图1右所示)，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{er}}^{n+1}={\mathrm{er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{ht}}^{n+0.5}\×l_{\mathrm{sub}})---\left(8\right)$

其中er为亚网格径向电场，ht为亚网格切向磁场，ΔS_sub、Δl_sub分别为亚网格面积、亚网格边长。

(3)采用公式(3)计算与亚网格相邻的主网格上的径向电场。

(4)计算主三角形网格上的切向电场Et。采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Et}}^{n+1}={\mathrm{Et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Ht}}^{n+0.5}\right)---\left(9\right)$

其中Et为主网格切向电场，Hr为主网格径向磁场，Ht为主网格上的切向磁场，δ_norm、Δr为主网格中心距离、径向空间步长。

(5)计算亚网格上(不包括与主网格交界处)的切向电场et(如图1左所示)，有：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(10\right)$

其中et为亚网格切向电场，hr为亚网格径向磁场，ht为亚网格上的切向磁场。δ_sub为亚网格中心距离。

(6)利用公式(2)的三维形式(补充径向的切向磁场)来计算边界处的亚网格切向电场：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{hr}}^{n+0.5}\mathrm{mH}{r}^{n+0.5})+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(11\right)$

至此，电场计算结束，下面计算磁场部分。

(7)计算主三角形网格中心(不包括与亚网格相邻)的径向磁场Hr，有

${\mathrm{Hr}}^{n+0.5}={\mathrm{Hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Et}}^n\×\mathrm{\Δ}{L}_{\mathrm{norm}})---\left(12\right)$

符号说明同上，不再复述。

(8)计算亚三角形网格中心的径向磁场hr，有

${\mathrm{hr}}^{n+0.5}={\mathrm{hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{et}}^n\×\mathrm{\Δ}{l}_{\mathrm{sub}})---\left(13\right)$

符号说明同上，不再复述。

(9)计算主六边形网格边上的切向磁场Ht，有

${\mathrm{Ht}}^{n+0.5}={\mathrm{Ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Et}}^n\right)---\left(14\right)$

(10)计算亚六边形网格边上的切向磁场ht，有

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(15\right)$

(11)利用公式(4)的三维形式(补充径向的切向电场)来计算边界处的亚网格切向磁场：

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{er}}^{n}m{\mathrm{Er}}^n)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(16\right)$

(12)利用公式(5)或者公式(6)计算虚拟切向磁场htv，如图2所示。

至此，全部的迭代过程完毕。

3)补充说明

(1)以上给出的计算步骤是基于三维地形学时域有限差分法的亚网格迭代过程，当计算模型为二维三角形网格时或二维六边形网格时，只需将不存在的场分量和迭代方程中相应的纵向切向场分量去掉即可。

(2)以上给出的步骤中默认电场位于六边形中心和三角形边上、磁场位于六边形边上和三角形中心，根据对偶性原则，可以将两种场分量互换，采用类似的迭代过程来实现亚网格计算方案。

(3)采用公式(5)或(6)得到的亚网格计算方案具有较好的稳定性特点，可以保证在4阶情况下的100000步稳定。

Claims (2)

1.一种极低频电磁波的仿真计算方法，其特征在于：包括亚网格划分步骤、电磁场的布置步骤、计算步骤；其中，

A、所述亚网格划分步骤如下：

A-1，基于三角形地形学网格划分：将原始主球面三角形网格的各个边的中点沿球面弧两两相连接，形成一种对主三角形网格的球面划分，划分出4个球面三角形作为亚网格元胞；

A-2，基于六边形地形学网格划分：将原始主球面六边形网格中相邻三个主六边形互相交接的三个边去掉，以三个面积大小接近为原始六边形1/4的三个新的六边形，划分出3个球面六边形作为亚网格元胞；

A-3，将所述步骤A-1划分的三角形亚网格及步骤A-2划分的六边形亚网格进行叠加，在径向上将三角形亚网格和六边形亚网格相结合，形成三维的亚网格元胞，成为对整个地球-电离层系统的完整亚网格划分；

B、所述电磁场的布置步骤如下：

B-1，对于基于三角形地形学网格划分，将被亚网格化的主网格上的1个径向场和3个切向场去除，以新的亚网格上的3个径向场和9个切向场分别布置电磁场；

B-2，对于基于六边形地形学网格划分，将被亚网格化的3个主网格上的相邻3条边上的3个切向场去除，以新的3个亚网格上的3个径向场和9个切向场、6个虚拟切向场分别布置电磁场；

B-3，将步骤B-1、步骤B-2的电磁场分布结合到一起，形成基于地形学网格划分的亚网格元胞；

C、计算步骤如下：

C-1，对于基于三角形地形学网格划分，在原始的主网格和新的亚网格内的电磁场计算沿用传统的经典时域有限差分算法，在原始的主网格和新的亚网格边界处的电磁场计算采用如下方法：

C-1-1，对于边界处主网格的径向磁场有：

${\mathrm{Hr}}^{n+0.5}={\mathrm{Hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{norm}}}\×\{\mathrm{\Σ}({\mathrm{Et}}^n\×\mathrm{\Δ}{l}_{\mathrm{norm}})+\mathrm{\Σ}({\mathrm{et}}^n\×\mathrm{\Δ}{l}_{\mathrm{sub}})\}---\left(1\right)$

其中Hr为主网格径向磁场，Et为主网格切向电场，et为亚网格上的切向电场，Δt为时间步长，ΔS_norm、Δl_norm、Δl_sub，分别为主网格面积、主网格边长、亚网格边长，μ为磁导率，n为时间步数，取值由所需运行时间决定；

C-1-2，对于边界处的亚网格切向电场有：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{hr}}^{n+0.5}mH{r}^{n+0.5})---\left(2\right)$

其中，hr为亚网格上的径向磁场，δ_sub为主网格中心到亚网格中心的距离，ε为介电系数；

C-2，对于基于六边形地形学网格划分，在原始的主网格和新的亚网格内的电磁场的计算沿用传统的经典时域有限差分算法，在原始的主网格和新的亚网格边界处的电磁场采用以下方法：

C-2-1，对于边界处主网格的径向电场有：

${\mathrm{Er}}^{n+1}={\mathrm{Er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{norm}}}\×\left\{\mathrm{\Σ}{\mathrm{Ht}}^{n+0.5}\×l_{\mathrm{norm}}\±\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\×0.5\×l_{\mathrm{sub}}\right\}---\left(3\right)$

其中Er为边界上主网格径向电场，Ht为主网格切向磁场，ht为亚网格上的切向磁场；

C-2-2，对于边界处的亚网格切向电场有：

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{er}}^{n}mE{r}^n)---\left(4\right)$

其中er为亚网格上的径向电场；

C-2-3，对于亚网格和主网格边界上的虚拟切向磁场有两种计算方法：

①，同向场近似法：

htv0≈Ht₀            (5)

②，周围场插值法：

$\mathrm{htv}0\≈-\frac{1}{2}\mathrm{ht}0+\frac{1}{6}({\mathrm{Ht}}_0+{\mathrm{Ht}}_1+\mathrm{htv}1)$ (6)

$\≈\frac{7}{35}{\mathrm{Ht}}_0+\frac{5}{35}{\mathrm{Ht}}_1+\frac{3}{35}\mathrm{ht}1-\frac{18}{35}\mathrm{ht}0$

其中，ht0，ht1为亚网格切向磁场，Ht₀，Ht₁为主网格切向磁场，htv0，htv1为亚网格虚拟磁场；

C-3，将步骤C-1及步骤C-2所述的两种平面二维算法叠加，得到三维亚网格算法，计算切向电场分量、径向电场分量、虚拟切向磁场分量，其中径向电场分量和虚拟切向磁场的计算方法与C-1，C-2所示步骤类同。

2.根据权利要求1所述的一种极低频电磁波的仿真计算方法，其特征在于：步骤C-3所述三维亚网络算法的具体计算步骤为：

(1)计算主六边形网格中心的径向电场Er，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Er}}^{n+1}={\mathrm{Er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Ht}}^{n+0.5}\×l_{\mathrm{norm}})---\left(7\right)$

(2)计算亚六边形网格中心的径向电场er，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{er}}^{n+1}={\mathrm{er}}^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{ht}}^{n+0.5}\×l_{\mathrm{sub}})---\left(8\right)$

(3)采用公式(3)计算与亚网格相邻的主网格上的径向电场；

(4)计算主三角形网格上的切向电场Et，采用积分麦克斯韦方程和有限差分技术，得到：

${\mathrm{Et}}^{n+1}={\mathrm{Et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Ht}}^{n+0.5}\right)---\left(9\right)$

其中δ_norm、Δr分别为主网格中心距离、径向空间步长。

(5)计算亚网格上的切向电场et：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{hr}}^{n+0.5}\right)+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(10\right)$

(6)利用公式(2)的三维形式来计算边界处的亚网格切向电场：

${\mathrm{et}}^{n+1}={\mathrm{et}}^n+\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{hr}}^{n+0.5}\mathrm{mH}{r}^{n+0.5})+\frac{\mathrm{\Δt}}{\mathrm{\ϵ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{ht}}^{n+0.5}\right)---\left(11\right)$

至此，电场计算结束，下面计算磁场部分：

(7)计算主三角形网格中心的径向磁场Hr，有

${\mathrm{Hr}}^{n+0.5}={\mathrm{Hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{norm}}}\×\mathrm{\Σ}({\mathrm{Et}}^n\×\mathrm{\Δ}{L}_{\mathrm{norm}})---\left(12\right)$

(8)计算亚三角形网格中心的径向磁场hr，有

${\mathrm{hr}}^{n+0.5}={\mathrm{hr}}^{n+0.5}-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δ}{S}_{\mathrm{sub}}}\×\mathrm{\Σ}({\mathrm{et}}^n\×\mathrm{\Δ}{l}_{\mathrm{sub}})---\left(13\right)$

(9)计算主六边形网格边上的切向磁场Ht，有

${\mathrm{Ht}}^{n+0.5}={\mathrm{Ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{norm}}}\×\left(\mathrm{\Σ}{\mathrm{Er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{Et}}^n\right)---\left(14\right)$

(10)计算亚六边形网格边上的切向磁场ht，有

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×\left(\mathrm{\Σ}{\mathrm{er}}^n\right)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(15\right)$

(11)利用公式(4)的三维形式来计算边界处的亚网格切向磁场：

${\mathrm{ht}}^{n+0.5}={\mathrm{ht}}^{n-0.5}-\frac{\mathrm{\Δt}}{{\mathrm{\μ\δ}}_{\mathrm{sub}}}\×(\±{\mathrm{er}}^{n}m{\mathrm{Er}}^n)-\frac{\mathrm{\Δt}}{\mathrm{\μ\Δr}}\×\left(\mathrm{\Σ}{\mathrm{et}}^n\right)---\left(16\right)$

(12)利用公式(5)或者公式(6)计算虚拟切向磁场htv；

至此，全部的迭代过程完毕。

Images