本发明属于飞机隐身设计领域，具体提供一种基于曲面三角形网格的飞机雷达散射截面积计算方法。
背景技术：
1、在飞机隐身设计领域，常采用电磁仿真方法得到飞机的雷达散射截面积，面积分方程方法(或边界元方法)作为一种高效数值方法被广泛应用；首先，对飞机表面进行网格(通常为三角形网格单元)离散，构造针对待求解电流系数的矩阵方程并进行求解，得到离散的电流系数，再基于电流系数求得远区散射电场，进而求得飞机雷达散射截面积。然而，在构造针对待求解电流系数的矩阵方程的过程中，计算阻抗矩阵时通常会出现难以处理的奇异性积分，这些奇异性积分往往不能直接利用数值积分来进行求解，并且这些奇异性积分的形式多种多样，很难通过一种方法来统一求解；此外，为了更好地拟合飞机的几何外形，通常需要采用曲面网格(一般为曲面三角形网格)进行离散，此时，这些奇异性积分定义在曲面网格区域上，使得其计算更加困难，导致待求解电流系数计算结果不准确，使得飞机雷达散射截面积计算结果存在较大误差，无法实现设计目标。
技术实现思路
1、本发明的目的在于提供一种基于曲面三角形网格的飞机雷达散射截面积计算方法，用以提高计算效率和计算精度。对于飞机雷达散射截面积计算过程中奇异性积分，本发明提出一种纯数值方法进行计算，该方法通过积分变量的函数变换，利用变换后雅克比因子弱化原始积分中被积函数的奇异性，实现数值积分时积分截断误差到零的快速收敛，从而有效的提高计算效率和计算精度，从而提高飞机雷达散射截面积的计算效率和计算精度。
2、为实现上述目的，本发明采用的技术方案为：
3、一种基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，包括以下步骤：
4、步骤1.采用曲面三角形网格离散飞机模型，得到所有曲面三角形网格单元的节点序号与节点坐标信息；
5、步骤2.根据网格信息建立针对待求解电流的矩阵方程，表示为：
6、[k]·{i}＝{b}
7、其中，k为阻抗矩阵，i为待求解电流系数向量，b为右端激励向量，[·]代表方阵，{·}代表向量；
8、以电场积分方程为例，k中第n行第m列的元素为knm，表示为：
9、
10、其中，r表示场点，rs表示源点，▽表示哈密顿算子，ω为角频率，k为波数，μ为磁导率，fn表示公共边n对应的测试函数，fm表示公共边m对应的基函数，表示格林函数，其梯度表示为：r＝|r-rs|。
11、步骤3.求解矩阵方程得到电流系数，再基于电流系数计算得到远区散射电场，进而计算得到飞机雷达散射截面积。
12、进一步的，步骤2中，阻抗矩阵k的计算过程为：
13、步骤2.1、计算场点与源点之间的距离，当距离小于等于预设阈值时，采用步骤2.2计算knm；否则，采用步骤2.3计算knm；
14、步骤2.2、首先，利用牛顿迭代法确定场点在源点所在曲面三角形网格单元上的投影点；然后，将曲面三角形网格单元投影到该投影点处的切平面上，得到投影三角形；再然后，将曲面三角形网格单元上的积分转换到投影三角形上，并对积分变量进行函数变换；最后，对变换后的积分变量采用高斯-勒让德积分方法进行计算，得到knm；
15、步骤2.3、采用高斯-勒让德积分方法计算得到knm。
16、更进一步的，步骤2.2中，利用牛顿迭代法确定场点在源点所在曲面三角形网格单元上的投影点，其具体过程为：
17、确定曲面三角形六个控制点的空间坐标：r1,r2,r3,r4,r5,r6，将该曲面三角形网格单元上的任意一点r′写成面积坐标，表示为：
18、r′(ξ1,ξ2,ξ3)＝ξ1(2ξ1-1)r1+ξ2(2ξ2-1)r2+ξ3(2ξ3-1)r3+4ξ1ξ2r4+4ξ2ξ3r5+4ξ1ξ3r6
19、其中，ξ1、ξ2与ξ3表示面积坐标，且ξ1+ξ2+ξ3＝1；r1,r2,r3,r4,r5,r6依次表示第1至第6个控制点的空间坐标；
20、利用六个控制点中三个顶点r1、r2与r3构成平面三角形，计算牛顿迭代方法中的初始点r0的面积坐标，计算方程为：
21、
22、其中，r表示场点的空间坐标；
23、求解上述方程得到初始点r0的空间坐标，计算得到初始点r0的面积坐标
24、
25、
26、
27、其中，表示该平面三角形的外法向单位矢量，area表示该平面三角形的面积；
28、设置牛顿迭代方法的收敛门限，通过牛顿迭代方法求解以下方程：
29、
30、得到投影点rp的面积坐标
31、再进一步的，步骤2.2中，将曲面三角形网格单元投影到该投影点处的切平面上，具体过程为：
32、将投影点rp的面积坐标转换为空间坐标，带入下式：
33、
34、
35、
36、计算得到r1,r2,r3在切平面上的投影点的空间坐标，利用所得到的三个投影点构造出所需的投影三角形。
37、再进一步的，步骤2.2中，将曲面三角形网格单元上的积分转换到投影三角形上，具体过程为：
38、在面积分方程方法(或边界元方法)中出现的奇异性积分可以统一表示成：
39、
40、其中，i表示积分结果；s′表示源点所在的积分区域，在本发明中规定该区域为一个曲面三角形网格单元；r表示场点，rs表示源点，r＝|r-rs|表示场点到源点的距离；
41、f(r,rs)表示在积分区域上没有奇异性的函数，以步骤2中积分∫s′fm(rs)g(r，rs)ds′为例，表示为：则对于积分表示为：则
42、根据场源点的位置将奇异性积分划分为真实奇异性积分(场点和源点位于同一网格单元内)和近奇异性积分(场点和源点十分靠近，但不在同一网格单元内)，根据积分核的表达式来划分的话，一般划分为弱奇异性积分(t＝1)和强奇异性积分(t＝2)，t＝1,2；
43、将曲面三角形网格单元上的强近奇异性积分转换到参数三角形中进行积分求解，表示为：
44、
45、其中，表示转换到参数三角形当中的雅克比因子；
46、将上述积分转换到该投影三角形区域当中，表示为：
47、
48、其中，a0表示投影三角形的面积，s″表示投影三角形区域；这样一来便将原来在曲面三角形区域s′上的奇异性积分转换到平面投影三角形区域s″上的奇异性积分。
49、再进一步的，步骤2.2中，对积分变量进行函数变换，其具体过程为：
50、以投影点rp为中心，将投影三角形区域s″划分成三个子三角形，表示为：
51、
52、其中，s″i,i＝1,2,3表示三个子三角形区域；
53、然后对于每个子三角形，构造以rp为原点、rp对立边为x轴、rp到对立边的高为y轴的局部坐标系，则对于一个子三角形区域上的积分表示为：
54、
55、其中，hi表示第i个子三角形中rp到对立边的高，xmax与xmin表示在x方向上积分的上限与下限；
56、对上述积分做一次duffy变换，即令u＝x/y，得到：
57、
58、其中，umax与umin表示在变量u上积分的上限与下限；
59、最后对积分变量u和y分别作以下函数变换：
60、
61、
62、其中，d表示投影点rp到场点r的距离；
63、则ii转换为：
64、
65、
66、其中，pmax与pmin表示在变量p上积分的上限与下限，qmax与qmin表示在变量q上积分的上限与下限。
67、再进一步的，步骤2.2中，对变换后的积分变量采用高斯-勒让德积分方法进行计算，其具体过程为：
68、将子三角形区域上的积分ii的积分区间映射到[0,1]×[0,1]上：
69、
70、
71、将上面积分写成高斯-勒让德积分加权和的形式，得到最终的结果：
72、
73、其中，n表示积分的采样点数，wα表示加权系数，xα为积分采样点，α＝l,k；为了表达简洁，采用f(xl,xk)表示积分ii中的被积函数。
74、基于上述技术方案，本发明的有益效果在于：
75、本发明提供一种基于曲面三角形网格的飞机雷达散射截面积计算方法，针对飞机雷达散射截面积计算过程中涉及的在曲面三角形区域上的强近奇异性积分，提出一种高效纯数值方法，首先将积分区域从曲面三角形区域转移到平面投影三角形区域，然后在该平面三角形区域上进行函数变化，使得最终积分中被积函数的奇异性被变换后的雅克比因子所弱化，达到直接使用高斯-勒让德积分求解的形式，最终实现强近奇异性积分的高效率、高精度计算；该方法作为纯数值方法，不受被积函数形式的影响，适用于弱奇异性积分或强奇异性积分，尤其能够处理曲面三角形区域上的奇异性积分，对于极端恶化的三角形(如包含极大钝角的三角形或极小锐角的三角形)区域，同样能够实现精确计算；并且，该方法具有更快的误差收敛速度，能够在较少积分采样点的情况下实现高计算精度；由此，本发明能够有效提高飞机雷达散射截面积的计算效率和计算精度。
1.一种基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，包括以下步骤：
2.根据权利要求1所述基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，步骤2.2中，利用牛顿迭代法确定场点在源点所在曲面三角形网格单元上的投影点，其具体过程为：
3.根据权利要求2所述基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，步骤2.2中，将曲面三角形网格单元投影到该投影点处的切平面上，其具体过程为：
4.根据权利要求3所述基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，步骤2.2中，将曲面三角形网格单元上的积分转换到投影三角形上，表示为：
5.根据权利要求4所述基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，步骤2.2中，对积分变量进行函数变换，其具体过程为：
6.根据权利要求5所述基于曲面三角形网格的飞机雷达散射截面积计算方法，其特征在于，步骤2.2中，对变换后的积分变量采用高斯-勒让德积分方法进行计算，其具体过程为：
