CN103310069A - 面向时域有限差分电磁计算的载体网格划分方法 - Google Patents

Abstract

本发明涉及一种面向时域有限差分电磁计算的载体网格划分方法，用离散电尺寸精度和载体表面逼近精度共同控制网格大小，所划分的网格模型与载体模型几何逼近程度高（特别是低频即大波长的情况下），以此网格模型进行时域有限差分电磁计算，计算结果较为准确；为划分网格模型提供了表面逼近精度控制参数，避免了现有的直接剖分方法网格划分精度控制的盲目性，使网格划分计算量得到有效掌控；当需要对网格加密时，只需对最底层黑结点作八叉树递归剖分即可，网格加密容易实现。

Description

面向时域有限差分电磁计算的载体网格划分方法

技术领域

本发明属于电磁兼容技术领域，涉及电磁场数值分析所用的网格划分方法，具体地说是面向时域有限差分电磁计算的载体网格划分方法，可用于时域有限差分法计算电磁场问题的研究。

背景技术

随着现代科学技术的快速发展，机动式通信系统的载体如车辆、飞机、舰船上装载了各种功能不同的电子设备及各种接收和发射信号的天线，由于各类无线电子设备和通信设备的密度的增加，导致通信之间的电磁信号更加密集，频谱重叠率逐渐增大，使得电磁环境日益恶化。无线设备之间的相互干扰以及无线电频率之间的相互影响这样的问题就必须得到重视。在这种电磁环境存在的情况下，为了确保通信设备能够正常发挥自己的功能，那么通信系统的电磁兼容问题就不能忽视。

机动式通信系统归结为载体多天线系统，研究这类系统的电磁兼容性具有很大的工程意义和实用价值。电磁兼容性预测技术能够有效地完成系统或设备的电磁兼容设计并且达到一个很好的效果。电磁兼容性预测技术是通过理论方法分析和评估系统或设备的电磁兼容的程度。这项技术通常应用在系统或设备研制的方案设计阶段和工程研制阶段。电磁兼容的主要目的是分析不兼容的相对比较薄弱的环节、评价系统的兼容安全裕度、为防护设计以及方案的修改提供依据。当系统或设备的功能设计方案初次形成以后，就可根据电磁兼容性的相关要求以及指标，对目标设计方案进行电磁兼容性预测分析、干扰程度计算、电磁耦合仿真计算以及对可能存在的干扰源以及敏感设备的电磁敏感度进行分析。电磁兼容性预测分析的应用是伴随着电磁兼容性设计开展的，贯穿于系统研制的整个过程，它具有计算速度快、成本低、参数修改方便、可以多次反复的计算、预测成功率较高等优点。

电磁兼容性预测技术的核心是电磁数值计算以及为数值计算提供离散数据的网格划分。电磁数值计算的任务是基于Maxwell方程组，建立逼近实际问题的连续型数学模型，然后采用相应的数值计算方法，经离散化处理，将连续型数学模型转化为等价的离散型数学模型，由离散数值构成代数方程组，应用有效的代数方程组解法，求解出该数学模型的数值解。现阶段主要的电磁数值计算方法为有限元法（FEM）、时域有限差分法（FDTD）、矩量法（MOM）等。对于这三种数值方法来说，具体实施时矩量法最难以实现，因为矩量法的积分方程比较复杂且基函数的选取也需一定的技巧。而有限元法在具体实施时需要注意稀疏矩阵的存储方式以及基函数的选择。与这两种方法比较，时域有限差分法的实施最简单。时域有限差分法直接将Maxwell方程在时间和空间上采用中心差分方式将其离散，形成一组在空间上采样、时间上递推的离散方程组。在时间轴上逐步推进的求解，有很好的稳定性和收敛性，在空间上它以Yee元胞为空间电磁场离散单元，将Maxwell方程转化为差分方程，是对Maxwell方程在离散空间的直接描述，结合计算机计数能够处理十分复杂的电磁问题。由于其离散方程递推格式上的特点，使其很容易实现高性能的并行运算。而且从应用方面来说，经过多年的发展，时域有限差分法已比较成熟，在很多方面已得到应用，其在天线辐射特性计算、微波电路分析、散射体雷达散射截面等方面有广泛的应用。与有限元法、矩量法相比，时域有限差分法很适合求解电磁的时域瞬态响应和宽频带响应问题。

目前关于时域有限差分的网格划分方法分为非均匀网格划分和均匀网格划分。在非均匀网格划分方面，论文“一种改进的FDTD网格剖分算法”（电讯技术，第50卷第5期，2010年5月，72-75，作者：陈晓羽等）通过读取模型获得轴线上的不连续分界点，将整个空间沿X轴、Y轴和Z轴方向各自分成多个区间；然后通过各区间的长度及其在模型中所处的位置，来确定该不同区间剖分时所采用的具体网格生成算法，由此得到一种新型的分区递变非均匀网格；论文“准三维时域有限差分中非均匀网格划分的高效新算法”（北京广播学院学报,第1卷第2期，1994年6月，30-37，作者：卿安永等）通过对计算区域按一定规则进行准均匀划分，在保证计算精度的前提下，削减了所占用的内存和运行时间。在均匀网格划分方面，论文“车载多天线系统的FDTD方法研究”（西安电子科技大学硕士学位论文，2009年9月，作者：周姗）通过构造车体及天线的包围长方体，基于空间位置枚举方法进行网格划分，采用辅助线段与载体面元平面的交点并判定交点是否在面元内的奇偶判别法来判定网格点的去留；论文“一种基于面元模型的FDTD自动网格产生技术”（计算机应用，第26卷第z2期，2006年12月，119-120，作者：普鑫）采用一个完全包含带剖分目标的长方体为计算空间，选定三个坐标轴方向的空间步长△x，△y，△z，计算空间离散为一系列网格，取网格的中心点为网格点，并假定当且仅当网格点位于目标内部时认为该网格被目标媒质完全填充；论文“基于CAD的电磁场数值计算中的建模方法”（西安电子科技大学硕士学位论文，2008年8月，作者：张友强）仅提出用八叉树对载体进行递归细分，缺乏实施过程描述。可以看出，已有研究成果存在以下问题：（1）对载体表面形状逼近精度低（特别是在低频即大波长的情况下），而载体形状对天线间的耦合度、天线方向图等电参数有较大的影响；（2）直接剖分导致在网格尺寸小时的网格划分计算量较大；（3）当需要对网格进行加密时，不能利用已有网格进行划分处理，而必须按新的网格尺寸或空间步长重新进行网格划分。

发明内容

本发明的目的在于克服上述现有网格划分中存在的问题，提供一种面向时域有限差分电磁计算的载体网格划分方法，以便能划分出既符合时域有限差分法离散电尺寸精度要求又满足载体表面逼近精度要求的均匀立方体网格。

本发明的目的是这样实现的：

面向时域有限差分电磁计算的载体网格划分方法，其特征是：包括如下步骤：

步骤101：将数据模型中的天线数据、三角形面元数据和四边形面元数据分别读入天线数据链表Ant_list，三角形面元数据链表Tria_list,四边形面元数据链表Quad_list中，并将三角形面元和四边形面元的所有顶点和边的数据信息分别存放在Vertex和Side链表中；

步骤102：计算天线数gw_num,三角形面元数tria_num和四边形面元数quad_num，将Vertex链表中重复的顶点去除，将Side链表中重复的边去除，去除重复的顶点数据存放在链表Vertexno中，去除重复的边数据存放在链表Sideno中，并将不重复顶点个数记为vertexno_num，不重复边条数记为sideno_num；

步骤103：通过遍历链表Ant_list、Tria_list、Quad_list计算载体天线模型在三维坐标系下的x、y、z坐标方向的各自最大值和最小值，以构造载体的逼近参考体，其中，z坐标方向的最大值取包含天线底端点个数最多的平面的z坐标值进行计算。逼近参考体的边界坐标记为che_xmax、che_xmin、che_ymax、che_ymin、che_zmax、che_zmin；

步骤104：遍历链表Ant_list、Tria_list、Quad_list，计算出载体天线模型在三维坐标系下的x、y、z坐标方向的最大值和最小值,以构造整个模型的包围盒，包围盒的边界坐标记为x_max,x_min,y_max,y_min,z_max,z_min。进一步计算出载体天线模型在x、y、z坐标方向的跨度值，即x_range、y_range、z_range，其中x_range=x_max-x_min，y_range=y_max-y_min,z_range=z_max-z_min；

步骤105：以包围盒的体中心点为中心，以x_range、y_range、z_range中最大值为边长d，建立第一个立方体结点，即八叉树的根结点，建立黑结点链表node_black和灰结点链表node_gray，保存根结点至灰结点链表node_gray；

步骤106：根据电磁波在媒质中的波长λ和包围盒边长d计算八叉树的理论递归深度值max_deep=int((1+lgd-lgλ)/lg2)，其中，max_deep的计算满足网格单元边长不超过0.1λ的约束关系；

步骤107：剖分链表node_gray和node_black中存储的网格结点（被剖分的结点称为父结点，剖分得到的八个结点称为子结点），将结点按三个方向的中垂面分割成八个体积相等的子结点，其中node_gray链表中的结点划分后，转至步骤108，node_black链表中的结点划分后，转至步骤113；

步骤108：判断子结点与载体面元是否相交，如果相交，转至步骤112，如果不相交，转至步骤109；

步骤109：判断子结点体中心点及八个顶点是否在载体之内，如果全部在载体内，转至步骤110，如果全部不在载体内，转至步骤111；

步骤110：结点为黑结点，转至步骤113；

步骤111：结点为白结点，转至步骤114；

步骤112：结点为灰结点，转至步骤115；

步骤113：子黑结点保存在node_black_re，转至步骤116；

步骤114：白结点直接去除；

步骤115：子灰结点保存在node_gray_re，转至步骤117；

步骤116：清空链表node_black，将node_black_re链表中的结点数据赋给node_black链表，清空链表node_black_re，转至步骤118；

步骤117：清空链表node_gray，将node_gray_re链表中的结点数据赋给node_gray链表，清空链表node_gray_re，转至步骤118；

步骤118：查看结点剖分是否已达到理论递归深度，如果结点剖分已达到理论递归深度，剖分停止，转至步骤119，否则，转至步骤107，继续剖分；

步骤119：从已达到理论递归深度的灰结点链表node_gray和黑结点链表node_black中，通过遍历链表中的结点数据，计算全部结点在x、y、z方向的最大值与最小值，即xmax、xmin、ymax、ymin、zmax、zmin，由此构成一个结点体包围盒；

步骤120：计算结点体包围盒与逼近参考体对应的六个面的逼近距离，即xmax-che_xmax、che_xmin-xmin、ymax-che_ymax、che_ymin-ymin、zmax-che_zmax、che_zmin-zmin；

步骤121：计算平均逼近距离aver=[(xmax-che_xmax)+(che_xmin-xmin)+(ymax-che_ymax)+(che_ymin-ymin)+(zmax-che_zmax)+(che_zmin-zmin)]/6；

步骤122：比较平均逼近距离aver与表面逼近精度ε（其大小可根据需要设置）的大小，如果aver≤ε，结点剖分最终结束，转至步骤123；否则，转至步骤107，继续剖分结点；

步骤123：将链表node_gray里的结点视为黑结点，然后合并链表node_black、node_gray中的结点信息，并保存到node_last_list链表中，该链表中的结点即为最终的网格划分模型。

所示的步骤109中空间点是否在载体内的判断方法，包括以下步骤：

步骤201：读取检测点坐标数据；

步骤202：将载体顶点与检测点比较，找出顶点数最少的一侧，以检测点为端点且沿该侧坐标轴方向做射线；

步骤203：确定过检测点且垂直于射线的平面，称此平面为截切面；

步骤204：判断载体所有顶点是否在平面射线方向的反方向一侧，若在，转步骤210，否则转步骤205；

步骤205：平面将载体分为两侧部分，保留射线方向的一侧；

步骤206：将平面与载体的边相交的交点按其连接关系形成若干个截切多边形，这些截切多边形与载体相应保留部分构成若干个多面体；

步骤207：计算过检测点的射线穿过多面体的个数；

步骤208：判断个数是否为奇数，个数为奇数时，转至步骤209，否则转至步骤210；

步骤209：检测点在载体内部；

步骤210：检测点不在载体内。

所述的步骤202中射线方向的确定，包含如下步骤：

步骤301：读取检测点的坐标数据；

步骤302：将载体n个顶点的x坐标值分别与检测点的x坐标值作比较，计算出小于检测点x坐标值的顶点数，记为nx，则大于或等于检测点x坐标的顶点数为n-nx；利用同样的方法计算出ny，n-ny，nz，n-nz；

步骤303：计算k=min（nx，n-nx，ny，n-ny，nz，n-nz）；

步骤304：判断k是否等于0，如果k=0，转至步骤305，如果k不等于0，转至步骤306；

步骤305：检测点不在载体内，转至步骤317；

步骤306：判断k是否等于nx，如果k=nx，转至步骤307，如果k不等于nx，转至步骤308；

步骤307：射线方向平行于x轴且指向x轴负方向，转至步骤317；

步骤308：判断k是否等于n-nx，如果k=n-nx，转至步骤309，如果k不等于n-nx，转至步骤310；

步骤309：射线方向平行于x轴且指向x轴正方向，转至步骤317；

步骤310：判断k是否等于ny，如果k=ny，转至步骤311，如果k不等于ny，转至步骤312；

步骤311：射线方向平行于y轴且指向y轴负方向，转至步骤317；

步骤312：判断k是否等于n-ny，如果k=n-ny，转至步骤313，如果k不等于n-ny，转至步骤314；

步骤313：射线方向平行于y轴且指向y轴正方向，转至步骤317；

步骤314：判断k是否等于nz，如果k=nz，转至步骤315，如果k不等于nz，转至步骤316；

步骤315：射线方向平行于z轴且指向z轴负方向，转至步骤317；

步骤316：射线方向平行于z轴且指向z轴正方向；

步骤317：将载体中包含k个顶点的所有边保存到链表SB1中，k个顶点数据保存到VB1中，并保存k值。

所述的步骤206中将载体保留部分与截切面构成若干个多面体，包含如下步骤：

步骤401：将步骤317中的链表SB1分成两部分，链表SBF1存放完全在截切面保留侧的边数据，链表SBF2存放与截切面相交的边数据；

步骤402：计算截切面与链表SBF2中的边相交的交点，并保存交点所在载体的边的数据信息；

步骤403：判断交点彼此之间的相邻关系；

步骤404：将所有交点中彼此之间具有相邻关系的交点归为一组，每一组交点对应一个截切多边形，若所有交点彼此之间都具有相邻关系，则截切面切割载体为一个截切多边形，否则为多个截切多边形，截切多边形保存在多边形链表polygon中；

步骤405：每一个截切多边形与载体相应保留部分构成一个多面体。

所述的步骤207中计算过检测点的射线穿过多面体的个数，包含如下步骤：

步骤601：遍历链表polygon内的所有多边形，判断检测点是否在多边形内，计算包含检测点的多边形个数；

步骤602：将包含检测点的多边形个数作为过检测点的射线穿过多面体的个数。

所述的步骤403中交点相邻关系的判断方法，包含如下步骤：

步骤501：读取两个交点及所在载体面元的边的数据信息；

步骤502：判断交点所在面元的边是否在载体同一个面元上，若交点所在面元的边在载体同一个面元上，转至步骤504，否则转至步骤503；

步骤503：两个交点不相邻；

步骤504：两个交点相邻。

本发明有如下优点：

（1）用离散电尺寸精度和载体表面逼近精度共同控制网格大小，所划分的网格模型与载体模型几何逼近程度高（特别是低频即大波长的情况下），以此网格模型进行时域有限差分电磁计算，计算结果较为准确；

（2）为划分网格模型提供了表面逼近精度控制参数，避免了现有的直接剖分方法网格划分精度控制的盲目性，使网格划分计算量得到有效掌控；

（3）当需要对网格加密时，只需对最底层黑结点作八叉树递归剖分即可，网格加密容易实现。

附图说明

图1是本发明的总流程图；

图2是空间点在载体内的判断流程图；

图3是确定以检测点为端点的射线方向流程图；

图4是构成截切多面体的流程图；

图5是交点相邻关系判断流程图；

图6是载体天线模型示例图；

图7是载体天线模型的立方体根结点示意图；

图8是第二次划分后的结点效果图；

图9是第三次划分后的结点效果图；

图10是第四次划分后的结点效果图；

图11是第五次划分后的结点效果图；

图12是第六次划分后的结点效果图，其中图12a是从载体前部观察到的效果图，图12b是从载体后部观察到的效果图。

具体实施方式

本发明的网格划分是在所建立的几何模型基础上进行的，几何模型包括载体模型和天线模型。载体模型是由三角形面元、四边形面元围成的封闭几何体，且三角形、四边形面元的顶点顺序由载体外观察时应为逆时针方向；天线模型为若干个与载体顶面相连接的直线段，每个直线段代表一根鞭天线并由底端和顶端两个顶点描述。

本发明研究载体的网格划分。按照时域有限差分法均匀网格划分的要求，所划分的立方体网格单元的边长不超过0.1λ（λ为电磁波在媒质中传播的波长）；同时，划分的网格体主要表面与载体主要表面的平均距离（即表面逼近精度）不超过ε（其大小可根据需要设置）。

本发明重点研究对载体天线模型的载体网格划分。本发明对载体模型的要求是：载体模型是由三角形面元、四边形面元围成的封闭几何体，且三角形、四边形面元的顶点顺序由载体外观察时应为逆时针方向。

参照图1，本发明的载体网格划分包括如下步骤：

步骤101：将数据模型中的天线数据、三角形面元数据和四边形面元数据分别读入天线数据链表Ant_list，三角形面元数据链表Tria_list,四边形面元数据链表Quad_list中，并将三角形面元和四边形面元的所有顶点和边的数据信息分别存放在Vertex和Side链表中；

步骤102：计算天线数gw_num,三角形面元数tria_num和四边形面元数quad_num，将Vertex链表中重复的顶点去除，将Side链表中重复的边去除，去除重复的顶点数据存放在链表Vertexno中，去除重复的边数据存放在链表Sideno中，并将不重复顶点个数记为vertexno_num，不重复边条数记为sideno_num；

步骤103：通过遍历链表Ant_list、Tria_list、Quad_list计算载体天线模型在三维坐标系下的x、y、z坐标方向的各自最大值和最小值，以构造载体的逼近参考体，其中，z坐标方向的最大值取包含天线底端点个数最多的平面的z坐标值进行计算。逼近参考体的边界坐标记为che_xmax、che_xmin、che_ymax、che_ymin、che_zmax、che_zmin；

步骤104：遍历链表Ant_list、Tria_list、Quad_list，计算出载体天线模型在三维坐标系下的x、y、z坐标方向的最大值和最小值,以构造整个模型的包围盒，包围盒的边界坐标记为x_max,x_min,y_max,y_min,z_max,z_min。进一步计算出载体天线模型在x、y、z坐标方向的跨度值，即x_range、y_range、z_range，其中x_range=x_max-x_min，y_range=y_max-y_min,z_range=z_max-z_min；

步骤105：以包围盒的体中心点为中心，以x_range、y_range、z_range中最大值为边长d，建立第一个立方体结点，即八叉树的根结点，建立黑结点链表node_black和灰结点链表node_gray，保存根结点至灰结点链表node_gray；

步骤106：根据电磁波在媒质中的波长λ和包围盒边长d计算八叉树的理论递归深度值max_deep=int((1+lgd-lgλ)/lg2)，其中，max_deep的计算满足网格单元边长不超过0.1λ的约束关系；

步骤107：剖分链表node_gray和node_black中存储的网格结点（被剖分的结点称为父结点，剖分得到的八个结点称为子结点），将结点按三个方向的中垂面分割成八个体积相等的子结点，其中node_gray链表中的结点划分后，转至步骤108，node_black链表中的结点划分后，转至步骤113；

步骤108：判断子结点与载体面元是否相交，如果相交，转至步骤112，如果不相交，转至步骤109；

步骤109：判断子结点体中心点及八个顶点是否在载体之内，如果全部在载体内，转至步骤110，如果全部不在载体内，转至步骤111；

步骤109中空间点是否在载体内的判断方法，参照图2，包括以下步骤：

步骤201：读取检测点坐标数据；

步骤202：将载体顶点与检测点比较，找出顶点数最少的一侧，以检测点为端点且沿该侧坐标轴方向做射线；

步骤203：确定过检测点且垂直于射线的平面（称此平面为截切面）；

步骤204：判断载体所有顶点是否在平面射线方向的反方向一侧，若在，转步骤210，否则转步骤205；

步骤205：平面将载体分为两侧部分，保留射线方向的一侧；

步骤206：将平面与载体的边相交的交点按其连接关系形成若干个截切多边形，这些截切多边形与载体相应保留部分构成若干个多面体；

步骤207：计算过检测点的射线穿过多面体的个数；

步骤208：判断个数是否为奇数，个数为奇数时，转至步骤209，否则转至步骤210；

步骤209：检测点在载体内部；

步骤210：检测点不在载体内。

步骤202中射线方向的确定，参照图3，包含如下步骤：

步骤301：读取检测点的坐标数据；

步骤302：将载体n个顶点的x坐标值分别与检测点的x坐标值作比较，计算出小于检测点x坐标值的顶点数，记为nx，则大于或等于检测点x坐标的顶点数为n-nx；利用同样的方法计算出ny，n-ny，nz，n-nz；

步骤303：计算k=min（nx，n-nx，ny，n-ny，nz，n-nz）；

步骤304：判断k是否等于0，如果k=0，转至步骤305，如果k不等于0，转至步骤306；

步骤305：检测点不在载体内，转至步骤317；

步骤306：判断k是否等于nx，如果k=nx，转至步骤307，如果k不等于nx，转至步骤308；

步骤307：射线方向平行于x轴且指向x轴负方向，转至步骤317；

步骤308：判断k是否等于n-nx，如果k=n-nx，转至步骤309，如果k不等于n-nx，转至步骤310；

步骤309：射线方向平行于x轴且指向x轴正方向，转至步骤317；

步骤310：判断k是否等于ny，如果k=ny，转至步骤311，如果k不等于ny，转至步骤312；

步骤311：射线方向平行于y轴且指向y轴负方向，转至步骤317；

步骤312：判断k是否等于n-ny，如果k=n-ny，转至步骤313，如果k不等于n-ny，转至步骤314；

步骤313：射线方向平行于y轴且指向y轴正方向，转至步骤317；

步骤314：判断k是否等于nz，如果k=nz，转至步骤315，如果k不等于nz，转至步骤316；

步骤315：射线方向平行于z轴且指向z轴负方向，转至步骤317；

步骤316：射线方向平行于z轴且指向z轴正方向；

步骤317：将载体中包含k个顶点的所有边保存到链表SB1中，k个顶点数据保存到VB1中，并保存k值。

步骤206中将载体保留部分与截切多边形构成若干个多面体，参照图4，包含如下步骤：

步骤401：将步骤317中的链表SB1分成两部分，链表SBF1存放完全在截切面保留侧的边数据，链表SBF2存放与截切面相交的边数据；

步骤402：计算截切面与链表SBF2中的边相交的交点，并保存交点所在载体的边的数据信息；

步骤403：判断交点彼此之间的相邻关系；

步骤404：将所有交点中彼此之间具有相邻关系的交点归为一组，每一组交点对应一个截切多边形，若所有交点彼此之间都具有相邻关系，则截切面切割载体为一个截切多边形，否则为多个截切多边形，截切多边形保存在多边形链表polygon中；

步骤405：每一个截切多边形与载体相应保留部分构成一个多面体。

步骤403中交点相邻关系的判断方法，参照图5，包含如下步骤：

步骤501：读取两个交点及所在载体面元的边的数据信息；

步骤502：判断交点所在面元的边是否在载体同一个面元上，若交点所在面元的边在载体同一个面元上，转至步骤504，否则转至步骤503；

步骤503：两个交点不相邻；

步骤504：两个交点相邻。

步骤207中计算过检测点的射线穿过多面体的个数，包含如下步骤：

步骤601：遍历链表polygon内的所有多边形，判断检测点是否在多边形内，计算包含检测点的多边形个数；

步骤602：将包含检测点的多边形个数作为过检测点的射线穿过多面体的个数。

步骤110：结点为黑结点，转至步骤113；

步骤111：结点为白结点，转至步骤114；

步骤112：结点为灰结点，转至步骤115；

步骤113：子黑结点保存在node_black_re，转至步骤116；

步骤114：白结点直接去除；

步骤115：子灰结点保存在node_gray_re，转至步骤117；

步骤116：清空链表node_black，将node_black_re链表中的结点数据赋给node_black链表，清空链表node_black_re，转至步骤118；

步骤117：清空链表node_gray，将node_gray_re链表中的结点数据赋给node_gray链表，清空链表node_gray_re，转至步骤118；

步骤118：查看结点剖分是否已达到理论递归深度，如果结点剖分已达到理论递归深度，剖分停止，转至步骤119，否则，转至步骤107，继续剖分；

步骤119：从已达到理论递归深度的灰结点链表node_gray和黑结点链表node_black中，通过遍历链表中的结点数据，计算全部结点在x、y、z方向的最大值与最小值，即xmax、xmin、ymax、ymin、zmax、zmin，由此构成一个结点体包围盒；

步骤120：计算结点体包围盒与逼近参考体对应的六个面的逼近距离，即xmax-che_xmax、che_xmin-xmin、ymax-che_ymax、che_ymin-ymin、zmax-che_zmax、che_zmin-zmin；

步骤121：计算平均逼近距离aver=[(xmax-che_xmax)+(che_xmin-xmin)+(ymax-che_ymax)+(che_ymin-ymin)+(zmax-che_zmax)+(che_zmin-zmin)]/6；

步骤122：比较平均逼近距离aver与表面逼近精度ε（其大小可根据需要设置）的大小，如果aver≤ε，结点剖分最终结束，转至步骤123；否则，转至步骤107，继续剖分结点；

步骤123：将链表node_gray里的结点视为黑结点，然后合并链表node_black、node_gray中的结点信息，并保存到node_last_list链表中，该链表中的结点即为最终的网格划分模型。

本发明所涉及的载体上的天线为鞭天线，鞭天线的截面尺寸远小于载体大小，为了减少计算量，提高计算效率，本发明对鞭天线采用细导线模型。通过半解析数值方法模拟细导线周围场的变化。细导线附近的环向磁场和径向电场假定均按1/r的规律变化，其中r为距导线中心的距离。对天线的模拟，需对天线附近包含鞭天线在内四个元胞的八个磁场分量进行处理，解决这个问题，需要通过导线内部及含有细导线元胞表面的环路场的表示以及法拉第定律求解。天线馈电方式选用同轴线馈电，该方法是将激励设置划分出来使它成为激励空间，在这个空间中有与天线馈线相同的传输线，而所要研究的天线结构处于另一个网格空间之内(天线结构空间)，通过在激励空间设置激励源，将能量耦合到天线网格中进行迭代计算，从而将激励网格空间与天线网格空间关联起来。而天线底端与载体表面相连，对载体与天线进行模拟，就必须将天线网格空间与载体网格空间关联起来，所以，网格划分时也就必须将天线包含在内。

本发明面向时域有限差分法实现对载体的网格划分，采用八叉树原理，划分的网格结点分为黑结点、灰结点和白结点三类。八叉树（octree）是一种用于描述三维空间的树状数据结构。八叉树的每个结点表示一个正方体的体积元素，每个结点有八个子结点，这八个子结点所表示的体积元素加在一起就等于父结点的体积。一般中心点作为结点的分叉中心。八叉树的定义是：若不是空树的话，那么树中的任一结点的子结点恰好只有八个或零个。在具体的剖分过程中，我们首先将载体天线结构所在空间用一个立方体来表示，即八叉树是从一个立方体开始的，首先将它划分为八个子立方体，然后再按照一定的规则继续地分解，将某些子立方体再划分为八个子立方体，依此递归划分，直到满足一定的划分终止条件。整个划分过程可以用一个树状的结构来表示，每个结点的子结点有8个或者0个。

仿真实例

本发明对图6所示的载体天线模型进行网格划分，该模型包含4根天线，编号为1、2、3、4，其中1号天线的长度为3米，2号、3号、4号天线长度为1.8米，载体面元为18个一般四边形面片，4根天线底端点同在一个四边形面片内。

遍历载体天线模型数据，计算模型在x，y，z坐标轴的跨度值，取其最大跨度值为边长，建立根结点，根结点的边长为4.79米（如图7所示）。

设表面逼近精度ε=0.02米，当波长λ=11.5米时，计算的理论递归深度为2。剖分结点，保留黑结点和灰结点，当递归深度为2时，黑结点和灰结点总数为16个（如图8所示），计算表面平均逼近距离aver为0.21米，未达到预定的精度0.02米；继续剖分，递归深度为3时，删除白结点，黑结点和灰结点总数为92个（如图9），计算表面平均逼近距离为0.111米，未达到预定的精度0.02米；继续剖分，递归深度为4时，黑结点和灰结点总数为656个（如图10所示），表面逼近距离为0.039米，但仍未满足要求；继续剖分，递归深度为5时，黑结点和灰结点总数为3930个（如图11所示），表面逼近距离仍为0.039米，此时，划分的网格逼近主要体现在载体斜面的逼近，而载体六个面（前侧面、后侧面、顶面、底面、左侧面、右侧面）的逼近缓慢，继续剖分结点后，当递归深度为6时，黑结点和灰结点总数为29994个（如图12所示，其中图12a是从载体前部观察到的效果图，图12b是从载体后部观察到的效果图），表面逼近距离为0.0165米，达到预定精度要求，划分结束。合并链表node_black、node_gray，此时，链表node_gray里的结点视为黑结点，输出黑结点得到最终的网格划分模型。可以看出，在设定合适的表面逼近精度时，所划分的网格模型达到了电尺寸离散与几何保形的统一，为时域有限差分计算提供了较为准确的模型。

Claims (6)

1.面向时域有限差分电磁计算的载体网格划分方法，其特征是：

包括如下步骤：

步骤101：将数据模型中的天线数据、三角形面元数据和四边形面元数据分别读入天线数据链表Ant_list，三角形面元数据链表Tria_list,四边形面元数据链表Quad_list中，并将三角形面元和四边形面元的所有顶点和边的数据信息分别存放在Vertex和Side链表中；

步骤102：计算天线数gw_num,三角形面元数tria_num和四边形面元数quad_num，将Vertex链表中重复的顶点去除，将Side链表中重复的边去除，去除重复的顶点数据存放在链表Vertexno中，去除重复的边数据存放在链表Sideno中，并将不重复顶点个数记为vertexno_num，不重复边条数记为sideno_num；

步骤103：通过遍历链表Ant_list、Tria_list、Quad_list计算载体天线模型在三维坐标系下的x、y、z坐标方向的各自最大值和最小值，以构造载体的逼近参考体，其中，z坐标方向的最大值取包含天线底端点个数最多的平面的z坐标值进行计算。逼近参考体的边界坐标记为che_xmax、che_xmin、che_ymax、che_ymin、che_zmax、che_zmin；

步骤104：遍历链表Ant_list、Tria_list、Quad_list，计算出载体天线模型在三维坐标系下的x、y、z坐标方向的最大值和最小值,以构造整个模型的包围盒，包围盒的边界坐标记为x_max,x_min,y_max,y_min,z_max,z_min。进一步计算出载体天线模型在x、y、z坐标方向的跨度值，即x_range、y_range、z_range，其中x_range=x_max-x_min，y_range=y_max-y_min,z_range=z_max-z_min；

步骤105：以包围盒的体中心点为中心，以x_range、y_range、z_range中最大值为边长d，建立第一个立方体结点，即八叉树的根结点，建立黑结点链表node_black和灰结点链表node_gray，保存根结点至灰结点链表node_gray；

步骤106：根据电磁波在媒质中的波长λ和包围盒边长d计算八叉树的理论递归深度值max_deep=int((1+lgd-lgλ)/lg2)，其中，max_deep的计算满足网格单元边长不超过0.1λ的约束关系；

步骤107：剖分链表node_gray和node_black中存储的网格结点，将结点按三个方向的中垂面分割成八个体积相等的子结点，其中node_gray链表中的结点划分后，转至步骤108，node_black链表中的结点划分后，转至步骤113；

步骤108：判断子结点与载体面元是否相交，如果相交，转至步骤112，如果不相交，转至步骤109；

步骤109：判断子结点体中心点及八个顶点是否在载体之内，如果全部在载体内，转至步骤110，如果全部不在载体内，转至步骤111；

步骤110：结点为黑结点，转至步骤113；

步骤111：结点为白结点，转至步骤114；

步骤112：结点为灰结点，转至步骤115；

步骤113：子黑结点保存在node_black_re，转至步骤116；

步骤114：白结点直接去除；

步骤115：子灰结点保存在node_gray_re，转至步骤117；

步骤116：清空链表node_black，将node_black_re链表中的结点数据赋给node_black链表，清空链表node_black_re，转至步骤118；

步骤117：清空链表node_gray，将node_gray_re链表中的结点数据赋给node_gray链表，清空链表node_gray_re，转至步骤118；

步骤118：查看结点剖分是否已达到理论递归深度，如果结点剖分已达到理论递归深度，剖分停止，转至步骤119，否则，转至步骤107，继续剖分；

步骤119：从已达到理论递归深度的灰结点链表node_gray和黑结点链表node_black中，通过遍历链表中的结点数据，计算全部结点在x、y、z方向的最大值与最小值，即xmax、xmin、ymax、ymin、zmax、zmin，由此构成一个结点体包围盒；

步骤120：计算结点体包围盒与逼近参考体对应的六个面的逼近距离，即xmax-che_xmax、che_xmin-xmin、ymax-che_ymax、che_ymin-ymin、zmax-che_zmax、che_zmin-zmin；

步骤121：计算平均逼近距离aver=[(xmax-che_xmax)+(che_xmin-xmin)+(ymax-che_ymax)+(che_ymin-ymin)+(zmax-che_zmax)+(che_zmin-zmin)]/6；

步骤122：比较平均逼近距离aver与表面逼近精度ε的大小，如果aver≤ε，结点剖分最终结束，转至步骤123；否则，转至步骤107，继续剖分结点；

步骤123：将链表node_gray里的结点视为黑结点，然后合并链表node_black、node_gray中的结点信息，并保存到node_last_list链表中，该链表中的结点即为最终的网格划分模型。

2.根据权利要求1所述的面向时域有限差分电磁计算的载体网格划分方法，其特征是：所示的步骤109中空间点是否在载体内的判断方法，包括以下步骤：

步骤201：读取检测点坐标数据；

步骤202：将载体顶点与检测点比较，找出顶点数最少的一侧，以检测点为端点且沿该侧坐标轴方向做射线；

步骤203：确定过检测点且垂直于射线的平面；

步骤204：判断载体所有顶点是否在平面射线方向的反方向一侧，若在，转步骤210，否则转步骤205；

步骤205：平面将载体分为两侧部分，保留射线方向的一侧；

步骤206：将平面与载体的边相交的交点按其连接关系形成若干个截切多边形，这些截切多边形与载体相应保留部分构成若干个多面体；

步骤207：计算过检测点的射线穿过多面体的个数；

步骤208：判断个数是否为奇数，个数为奇数时，转至步骤209，否则转至步骤210；

步骤209：检测点在载体内部；

步骤210：检测点不在载体内。

3.根据权利要求2所述的面向时域有限差分电磁计算的载体网格划分方法，其特征是：所述的步骤202中射线方向的确定，包含如下步骤：

步骤301：读取检测点的坐标数据；

步骤302：将载体n个顶点的x坐标值分别与检测点的x坐标值作比较，计算出小于检测点x坐标值的顶点数，记为nx，则大于或等于检测点x坐标的顶点数为n-nx；利用同样的方法计算出ny，n-ny，nz，n-nz；

步骤303：计算k=min（nx，n-nx，ny，n-ny，nz，n-nz）；

步骤304：判断k是否等于0，如果k=0，转至步骤305，如果k不等于0，转至步骤306；

步骤305：检测点不在载体内，转至步骤317；

步骤306：判断k是否等于nx，如果k=nx，转至步骤307，如果k不等于nx，转至步骤308；

步骤307：射线方向平行于x轴且指向x轴负方向，转至步骤317；

步骤308：判断k是否等于n-nx，如果k=n-nx，转至步骤309，如果k不等于n-nx，转至步骤310；

步骤309：射线方向平行于x轴且指向x轴正方向，转至步骤317；

步骤310：判断k是否等于ny，如果k=ny，转至步骤311，如果k不等于ny，转至步骤312；

步骤311：射线方向平行于y轴且指向y轴负方向，转至步骤317；

步骤312：判断k是否等于n-ny，如果k=n-ny，转至步骤313，如果k不等于n-ny，转至步骤314；

步骤313：射线方向平行于y轴且指向y轴正方向，转至步骤317；

步骤314：判断k是否等于nz，如果k=nz，转至步骤315，如果k不等于nz，转至步骤316；

步骤315：射线方向平行于z轴且指向z轴负方向，转至步骤317；

步骤316：射线方向平行于z轴且指向z轴正方向；

步骤317：将载体中包含k个顶点的所有边保存到链表SB1中，k个顶点数据保存到VB1中，并保存k值。

4.根据权利要求2所述的面向时域有限差分电磁计算的载体网格划分方法，其特征是：所述的步骤206中将载体保留部分与截切面构成若干个多面体，包含如下步骤：

步骤401：将步骤317中的链表SB1分成两部分，链表SBF1存放完全在截切面保留侧的边数据，链表SBF2存放与截切面相交的边数据；

步骤402：计算截切面与链表SBF2中的边相交的交点，并保存交点所在载体的边的数据信息；

步骤403：判断交点彼此之间的相邻关系；

步骤404：将所有交点中彼此之间具有相邻关系的交点归为一组，每一组交点对应一个截切多边形，若所有交点彼此之间都具有相邻关系，则截切面切割载体为一个截切多边形，否则为多个截切多边形，截切多边形保存在多边形链表polygon中；

步骤405：每一个截切多边形与载体相应保留部分构成一个多面体。

5.根据权利要求2所述的面向时域有限差分电磁计算的载体网格划分方法，其特征是：所述的步骤207中计算过检测点的射线穿过多面体的个数，包含如下步骤：

步骤601：遍历链表polygon内的所有多边形，判断检测点是否在多边形内，计算包含检测点的多边形个数；

步骤602：将包含检测点的多边形个数作为过检测点的射线穿过多面体的个数。

6.根据权利要求4所述的面向时域有限差分电磁计算的载体网格划分方法，其特征是：所述的步骤403中交点相邻关系的判断方法，包含如下步骤：

步骤501：读取两个交点及所在载体面元的边的数据信息；

步骤502：判断交点所在面元的边是否在载体同一个面元上，若交点所在面元的边在载体同一个面元上，转至步骤504，否则转至步骤503；

步骤503：两个交点不相邻；

步骤504：两个交点相邻。

Images