用快速边界元法得到大型复杂飞行器电场分布的方法

摘要

本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法，技术特征在于：采用边界元法解决飞行器设计领域所存在的飞行器大规模静电场的快速计算分析问题。采用奈氏离散化方法，把介质界面划分成具有参数方程的光滑曲面，通过参数方程把积分点投影到界面上形成边界积分点；基于积分点的多尺度分组构造小波变换矩阵，对边界元系数矩阵进行矩阵压缩、计算和存储。克服了现有小波边界元法存在的系数矩阵计算复杂、计算量大、精度不高的缺点，有效降低了边界元分析的内存占用量和计算时间，使其可解决自由度上百万的大规模工程问题。本发明可应用于航空航天领域的大型复杂飞行器电场计算分析，以及其它工程领域的电场分析问题。

说明书

技术领域

本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法，飞行器设计领域，具体属于飞行器大规模电磁场的快速分析设计方法。

背景技术

在航空航天领域，飞行器结构的充、放电现象会干扰飞行器的通讯、导航，甚至造成飞行器结构破坏等严重后果。因此，对飞行器结构充、放电过程的计算分析是飞行器设计中必须考虑的问题，其核心是通过求解拉普拉斯方程，确定飞行器的静电场。这个问题具有求解区域复杂、求解规模大的特点。

现有飞行器静电场的计算方法大致有两类：

第一类是采用有限元法。由于飞行器结构的静电场本质上是一个无限域问题，采用有限元法时必须通过设置人工边界，得到包含飞行器的有限空间作为分析区域，计算出整个区域内的电学量(如电势、电荷密度等)，而实际用到的只有飞行器表面的电学量。因此，这种方法的计算效率和精度较低，而且较费时。

第二类是采用常规边界元法，主要是求解定义在飞行器结构中不同介质界面上的拉普拉斯边界积分方程。与有限元法相比，该方法的自由度数目小、精度高。但常规边界元法形成稠密的系数矩阵，求解过程的内存占用量和计算时间与自由度数目N的平方成比例，即为O(N²)量级，无法应用到飞行器大规模静电场的快速计算分析中。

现有的降低常规边界元法内存占用量和计算时间的方法也可归为两类：一类是基于低秩逼近的，如快速多极方法、预修正快速傅里叶变换方法等，这类方法的特点是不显式计算系数矩阵，只是在迭代求解过程中对系数矩阵进行分块低秩分解，并完成矩阵-向量乘积运算，但它迭代一次的时间较长，当系数矩阵性态不好时，需要的迭代次数较多，因此总的计算时间会很长。另一类是基于小波压缩的，这类方法在计算系数矩阵之前首先对其进行小波压缩，大大减少非零元素的数目，然后再显式计算并存储一个稀疏的系数矩阵，从而达到降低边界元法的内存占用量和计算时间的目的。与低秩逼近方法相比，小波压缩方法具有迭代时间短(大约小于前者的1％)，系数矩阵预处理简便等优势。

但是，现有基于小波压缩的边界元法(以下简称“小波边界元法”)中都采用配点法和伽辽金(Galerkin)法作为离散化方法，主要存在三个问题：①需要计算单重或二重曲面积分，奇异积分处理复杂，系数矩阵计算量较大；②受边界剖分和形函数逼近方法的限制，难于得到高精度的结果；③飞机、航天器等的表面通常由大量光滑曲面片组成，划分单元会人为破坏边界的连续性，所以配点法和伽辽金法的计算效率较低。现有基于配点法和伽辽金法的小波边界元法所存在的系数矩阵计算方法复杂、计算量大、精度受单元剖分所限制、对复杂飞行器电场计算效率较低的不足。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种用快速边界元法得到大型复杂飞行器电场分布的方法，是一种基于奈氏离散化方法的小波边界元法，并将其应用于大型复杂飞行器的电场计算分析中，以有效降低计算量、提高结果精度和求解效率。

技术方案

一种用快速边界元法得到大型复杂飞行器电场分布的方法，其特征在于步骤如下：

步骤1：按照飞行器的实际尺寸建立几何模型：x＝Υ_i(ξ)，i＝1，…，N_s，其中x为飞行器边界坐标，ξ为参数坐标，Υ_i为第i块曲面的参数方程，N_s飞行器边界曲面个数，所述的飞行器边界曲面为光滑的；

步骤2在每块光滑曲面上布置积分点：将参数的取值范围ξ变换到标准参考区域[0，1]×[0，1]上，以参数方程x＝Υ_i(ξ)将标准参考区域[0，1]×[0，1]上的Gauss积分点ξ_k(k＝1，…，n_i)投影到飞行器的边界曲面上，得到x_k＝Υ_i(ξ_k)和与x_k对应的Gauss权系数ω_k；所述的Gauss积分点数n_i大于10；

步骤3积分点分组：取包含飞行器边界曲面上所有积分点的正方体格子，对格子进行逐层细分，直至每个格子中包含的积分点数目小于给定常数n₀，形成2^d叉树结构；所述2^d叉树根节点正方体格子的层数为0，最底层为最小的积分点分组、层序号为L；所述给定常数n₀取15～30；所述d为问题的维数，取为3；

步骤4构造小波变换矩阵Q_v：求出2^d叉树结构中每个结点对应的积分点组v所对应的小波变换矩阵Q_v，步骤如下：

步骤a：求最底层上的矩量矩阵：[M_v]_α，i＝(x_i-x_v)^α，|α|≤p_l，其中v为任意第L层格子，α为三重指标，x_i为v上第i个积分点，格子中心为x_v；所述p_l为l层上的小波消失矩阶数；

步骤b：对矩量矩阵M_v做奇异值分解得到v上的小波变换矩阵尺度函数变换矩阵和尺度函数矩量矩阵其中，和由Q_v中分别与零奇异值和非零奇异值对应的列组成，包含∑_v中与所有非零奇异值对应的列；

步骤c：求高层上的矩量矩阵：其中，v为第l层格子，l＝L-1，…，1，μ1，…μs表示v的所有子格子，T_v，μ为矩转换矩阵，C_α^β为二项式系数；

步骤d：求M_v的奇异值分解可得v的小波变换矩阵尺度函数变换矩阵和尺度函数矩量矩阵$M_v^{\phi }=U_v{\tilde{\Sigma }}_v;$

步骤5计算边界元系数矩阵A：

采用小波边界元法的非标准型离散方法，系数矩阵A的元素分布为

其中，子矩阵和(l＝L，…，1)；

所述A₁的计算步骤如下：

步骤I：按关系Near(v)＝{v′:level(v)＝level(v′)，dist(v，v′)＜ηmax{diam(v)，diam(v′)}}找出每个格子v的所有临近格子集合Near(v)，其中，η为给定常数，dist(v，v′)为v和v′的距离，diam(v)为v的特征尺寸，level(v)为v所在的层数；所述的η小于0.5；

步骤II：对第L层上任意两个临近的格子v和v′∈Near(v)，计算矩阵

其中，K_ij＝K(x_v，i，x_v′，j)，K(x，y)为静电场积分方程的积分核函数，x_v，i为v中的序号为i的积分点，n_v为v中积分点总数；

步骤III：计算子矩阵和其中的元素按如下方法计算

${\hat{Q}}_v^T{A}_{v,v^{\prime }}{\hat{Q}}_{v^{\prime }}={\hat{I}}_{v,v^{\prime }},$

${\hat{Q}}_v^T{A}_{v,v^{\prime }}{\tilde{Q}}_{v^{\prime }}={\tilde{I}}_{v,v^{\prime }},$

${\tilde{Q}}_v^T{A}_{v,v^{\prime }}{\hat{Q}}_{v^{\prime }}={\bar{I}}_{v,v^{\prime }},$

${\tilde{Q}}_v^T{A}_{v,v^{\prime }}{\tilde{Q}}_{v^{\prime }}=I_{v,v^{\prime }}$

得到的和分别为矩阵和中与格子v和v′对应的子块，为

和的结构与相同；

步骤IV：对第l层(l＝L-1，…，1)上任意两个临近的格子v和v′，计算矩阵

其中，μ和μ′分别表示v和v′的子格子，所述的其中，矩阵A_μ，μ′采用递归方法计算：

当μ和μ′的层数为L时，计算方法采用步骤II；

当μ和μ′的层数为l＝L-1，…，1时，计算方法采用步骤IV；

步骤V：采用步骤III计算子矩阵和(l＝L-1，…，1)，以及矩阵A₁：

其中，v_i和v′_j为第1层上的相邻格子；

步骤6：设定线性方程组Ax＝b，其中b的元素为飞行器边界曲面上每个积分点上给定的电压或电荷密度的数值，A为边界元系数矩阵，x为待求向量；所述待求向量x为飞行器边界曲面上待求的电压或电荷密度的数值；

步骤7：采用迭代方法快速求解线性方程组Ax＝b，得到x。

有益效果

本发明提出的用快速边界元法得到大型复杂飞行器电场分布的方法，有益效果是：

1、通过采用奈氏离散化方法，避免了数值积分的计算，简化了边界元法的系数矩阵计算过程。

2、通过应用电学结构中不同介质界面的参数方程，实现了边界元法中边界曲面的精确表示，避免了常规边界元法中的单元剖分，提高了边界元分析的精度。

3、通过对边界元系数矩阵的小波变换和矩阵压缩，将其转化为稀疏矩阵，非零元素数目降低到O(N)量级，大大降低了边界元分析的内存占用量和计算时间，因此本发明可广泛用于大型复杂飞行器的大规模电场快速分析中，提高设计分析效率数十倍以上。

利用本发明的基本原理和计算方法，可以突破现有小波边界元法所存在的系数矩阵计算复杂、计算量大、精度不高的限制，并可以有效降低边界元分析的内存占用量和计算时间，扩大边界元分析的规模，对航空航天领域的大型复杂飞行器的电场计算分析，以及其它工程领域的电场分析问题具有重要意义。

附图说明

图1是本发明的流程图；

图2是飞行器边界的参数曲面轮廓图；

图3是飞行器翼面的一个单元上、与7点高斯积分公式对应的所有积分点分布图；

图4(a)是包含飞行器上所有积分点的第0层格子及其一次细分的示意图；图4(b)是第0层格子右-前-下方位子格子一次细分的示意图；图4(c)是两次细分对应的八叉树结构示意图；

图5(a)是第L层格子中积分点的分布示意图；图5(b)是高层格子与其子格子的关系示意图；

图6是飞行器电场数值方法所得到的系数矩阵的稀疏结构图，其中黑色区域为非零元素，其余均为零元素。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

以某飞行器模型的电容计算问题为例，按照本发明技术方案进行实施，给出了详细的实施过程。

计算一飞行器模型的电容，核心是在飞行器边界S上求解边界积分方程

∫_SK(x，y)σ(y)dy＝f(x)

其中，K(x，y)＝1/(4πε₀|x-y|)，ε₀＝8.8541878为真空电容率，f(x)为已知函数，σ(y)为待求函数。应用飞行器电场分析的快速数值方法的主要步骤如下：

1.建立几何模型。可利用CAD软件(如Pro-E等)，按照实际尺寸建立飞行器边界曲面模型，并将其保存为IGES格式，其中包含了每块曲面的参数方程，图2(a)显示了飞行器边界各参数曲面的轮廓；

2.布置积分点。从IGES模型文件中读取边界曲面参数方程，将每块曲面剖分为若干曲面四边形子块(称之为“单元”)，保证每个单元有独立的参数方程，如图2(b)所示，总单元数目为122。将单元参数方程转化到标准参考区域[0，1]×[0，1]上。根据计算精度要求选取合适的高斯积分点数，并将定义在[0，1]×[0，1]上的高斯积分点通过参数方程投影到单元上，形成单元积分点。图3给出了飞行器翼面上一个单元、与7点高斯积分公式对应的所有积分点的分布情况；

3.积分点分组。本实施例是三维问题(d＝3)，可按以下步骤建立积分点的多尺度分组：

(1)选取包含飞行器上所有积分点的正方体格子，如图4(a)所示，将其层数设为0，显然此格子非空(包含积分点的格子是非空的)，将其作为八叉树的根结点；

(2)层数l从0开始递增，如果某个第l层格子所包含的积分点数目大于某个给定的数目n₀＝30，则将其由各边中点等分为8个子格子，略去不包含积分点的格子，就得到了第l+1层格子；

(3)重复步骤(2)中的细分过程，直至每个最小格子中包含的积分点数目都不大于n₀，设最大层数为L。

以上分组过程可由图4说明。图4(a)显示了第0层格子及其一次细分的情况，图4(b)为第0层格子右-前-下方位的子格子的一次细分情况，图4(c)是由上述多尺度细分所建立的八叉树结构，层数L＝2，它具有两个特征：一是每个结点中包含一定数目的积分点，二是所有子格子中积分点之和就是父格子中的积分点。

4.构造小波变换矩阵Q_v。由于L层至第1层，逐层构造。取小波消失矩阶数与层数的关系为p_l＝L-l+1。如下两步构造：

(1)计算L＝2层()上的Q_v。该层格子中包含的积分点数目均不大于n₀＝30，消失矩阶数为p₂＝1。设该层格子v中的积分点为x₁，x₂，…，x_n，如图5(a)所示。计算矩量矩阵M_v：

$M_v=\left(\begin{array}{cccc}{(x_1-x_v)}^{\left(\mathrm{0,0,0}\right)}& {(x_2-x_v)}^{\left(\mathrm{0,0,0}\right)}& ...& {(x_n-x_v)}^{\left(\mathrm{0,0,0}\right)}\\ {(x_1-x_v)}^{\left(\mathrm{1,0,0}\right)}& {(x_2-x_v)}^{\left(\mathrm{1,0,0}\right)}& ...& {(x_n-x_v)}^{\left(\mathrm{1,0,0}\right)}\\ {(x_1-x_v)}^{\left(\mathrm{0,1,0}\right)}& {(x_2-x_v)}^{\left(\mathrm{0,1,0}\right)}& ...& {(x_n-x_v)}^{\left(\mathrm{0,1,0}\right)}\\ {(x_1-x_v)}^{\left(\mathrm{0,0,1}\right)}& {(x_2-x_v)}^{\left(\mathrm{0,0,1}\right)}& ...& {(x_n-x_v)}^{\left(\mathrm{0,0,1}\right)}\end{array}\right)$

作奇异值分解得到v的小波变换矩阵尺度函数变换矩阵同时还得到v的尺度函数的矩量矩阵存储Q_v和M_v^φ。

(2)计算1层上的Q_v。该层消失矩阶数为2。对每个格子v，设其子格子为μ1，…μs，如图5(b)所示。以步骤(1)为初始，每个子格子中已计算出尺度函数的矩量矩阵，即M_η1^φ，…，M_ηs^φ。于是可由式

$M_v=[T_{v,\mu 1},...,T_{v,\mathrm{\mu s}}]\mathrm{diag}\{M_{\mu 1}^{\phi },...,M_{\mathrm{\mu s}}^{\phi }\}$

计算当前格子v的矩量矩阵M_v。上式中转换矩阵T_v，μi只与v和μ的中心坐标有关：

${\left[T_{v,\mu }\right]}_{\alpha ,\beta }=C_{\alpha }^{\beta }{(x_{\mu }-x_v)}^{\alpha -\beta }$

同样，计算即可得v的小波变换矩阵Q_v和尺度函数的矩量矩阵$M_v^{\phi }=U_v{\tilde{\Sigma }}_v.$

通过以上两步，可对八叉树结构中除根结点外的每个结点构造出一个小波变换矩阵Q_v和尺度函数矩量矩阵M_v^φ，它们将在下面的边界元系数矩阵计算中用到。

5.计算边界元系数矩阵A，步骤如下：

(1)对每个非空格子v，计算临近格子集合Near(v)，其中参数η＝0.4；

(2)对第2层上任意两个临近的格子v和v′∈Near(v)，计算矩阵A_v，v′，其元素为

[A_v，v′]_ij＝1/(4πε₀|x_v，i-x_v′，j|)

其中，x_v，i为v中的序号为i的积分点；

(3)计算子矩阵和其中的元素按如下方法计算

${\hat{Q}}_v^T{A}_{v,v^{\prime }}{\hat{Q}}_{v^{\prime }}={\hat{I}}_{v,v^{\prime }},$

${\hat{Q}}_v^T{A}_{v,v^{\prime }}{\tilde{Q}}_{v^{\prime }}={\tilde{I}}_{v,v^{\prime }},$

${\tilde{Q}}_v^T{A}_{v,v^{\prime }}{\hat{Q}}_{v^{\prime }}={\bar{I}}_{v,v^{\prime }},$

${\tilde{Q}}_v^T{A}_{v,v^{\prime }}{\tilde{Q}}_{v^{\prime }}=I_{v,v^{\prime }}$

得到的和分别为矩阵和中与格子v和v′对应的子块。

(4)对第1层上任意两个临近的格子v和v′，计算矩阵

其中，μ和μ′分别表示v和v′的子格子，它们的层数为L＝2，I_μ，μ′的计算方法为

$I_{\mu ,{\mu }^{\prime }}={\tilde{Q}}_{\mu }^T{A}_{\mu ,{\mu }^{\prime }}{\tilde{Q}}_{{\mu }^{\prime }}$

其中，矩阵A_μ，μ′按照第(2)步的方法计算，即

[A_μ，μ′]_ij＝1/(4πε₀|x_μ，i-x_μ′，j|)

(5)按与第(3)步相同的方法计算子矩阵和A₁。

按上述步骤计算出的边界元系数矩阵具有如图6所示的稀疏结构，其中黑色区域为非零元素，其余均为零元素。

6.设定线性方程组Ax＝b。计算飞行器电容时，在边界上施加1V电压，即积分方程中f(x)＝1，于是，向量b的元素均为1。

7.求解线性方程组Ax＝b得到x，其中包含飞行器边界所有积分点上的电荷密度σ(x_i)数值。飞行器电容为

C＝∫_Sσ(x)dx＝∑_iw_iσ(x_i)

式中，w_i为已知的积分权系数。