专利名称：用快速边界元法得到大型复杂飞行器电场分布的方法
技术领域：
本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法，飞行器设 计领域，具体属于飞行器大规模电磁场的快速分析设计方法。
背景技术：
在航空航天领域，飞行器结构的充、放电现象会干扰飞行器的通讯、导航，甚至造成飞行器结构破坏等严重后果。因此，对飞行器结构充、放电过程的计算分析是飞行器设计 中必须考虑的问题，其核心是通过求解拉普拉斯方程，确定飞行器的静电场。这个问题具有 求解区域复杂、求解规模大的特点。现有飞行器静电场的计算方法大致有两类第一类是采用有限元法。由于飞行器结构的静电场本质上是一个无限域问题，采 用有限元法时必须通过设置人工边界，得到包含飞行器的有限空间作为分析区域，计算出 整个区域内的电学量(如电势、电荷密度等)，而实际用到的只有飞行器表面的电学量。因 此，这种方法的计算效率和精度较低，而且较费时。第二类是采用常规边界元法，主要是求解定义在飞行器结构中不同介质界面上的 拉普拉斯边界积分方程。与有限元法相比，该方法的自由度数目小、精度高。但常规边界 元法形成稠密的系数矩阵，求解过程的内存占用量和计算时间与自由度数目N的平方成比 例，即为0(N2)量级，无法应用到飞行器大规模静电场的快速计算分析中。现有的降低常规边界元法内存占用量和计算时间的方法也可归为两类一类是基 于低秩逼近的，如快速多极方法、预修正快速傅里叶变换方法等，这类方法的特点是不显式 计算系数矩阵，只是在迭代求解过程中对系数矩阵进行分块低秩分解，并完成矩阵-向量 乘积运算，但它迭代一次的时间较长，当系数矩阵性态不好时，需要的迭代次数较多，因此 总的计算时间会很长。另一类是基于小波压缩的，这类方法在计算系数矩阵之前首先对其 进行小波压缩，大大减少非零元素的数目，然后再显式计算并存储一个稀疏的系数矩阵，从 而达到降低边界元法的内存占用量和计算时间的目的。与低秩逼近方法相比，小波压缩方 法具有迭代时间短(大约小于前者的)，系数矩阵预处理简便等优势。但是，现有基于小波压缩的边界元法(以下简称“小波边界元法”)中都采用配 点法和伽辽金(Galerkin)法作为离散化方法，主要存在三个问题①需要计算单重或二重 曲面积分，奇异积分处理复杂，系数矩阵计算量较大；②受边界剖分和形函数逼近方法的限 制，难于得到高精度的结果；③飞机、航天器等的表面通常由大量光滑曲面片组成，划分单 元会人为破坏边界的连续性，所以配点法和伽辽金法的计算效率较低。现有基于配点法和 伽辽金法的小波边界元法所存在的系数矩阵计算方法复杂、计算量大、精度受单元剖分所 限制、对复杂飞行器电场计算效率较低的不足。
发明内容
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种用快速边界元法得到大型复杂飞 行器电场分布的方法，是一种基于奈氏(Nystr0m)离散化方法的小波边界元法，并将其应 用于大型复杂飞行器的电场计算分析中，以有效降低计算量、提高结果精度和求解效率。技术方案一种用快速边界元法得到大型复杂飞行器电场分布的方法，其特征在于步骤如 下步骤1 按照飞行器的实际尺寸建立几何模型χ = !,(ξ),! = 1,…，Ns，其中χ 为飞行器边界坐标，I为参数坐标，Ti为第i块曲面的参数方程，Ns飞行器边界曲面个数， 所述的飞行器边界曲面为光滑的；步骤2在每块光滑曲面上布置积分点将参数的取值范围ξ变换到标准参考区域
X W，1]上，以参数方程Χ=Κξ)将标准参考区域
X
上的Gauss积分 点€k(k= 1，…，ni)投影到飞行器的边界曲面上，得到Xk= Ti(Ik)和与Xk对应的Gauss 权系数；所述的Gauss积分点数Iii大于10 ；步骤3积分点分组取包含飞行器边界曲面上所有积分点的正方体格子，对格子 进行逐层细分，直至每个格子中包含的积分点数目小于给定常数IV形成2d叉树结构；所述 2d叉树根节点正方体格子的层数为0，最底层为最小的积分点分组、层序号为L ；所述给定 常数Iitl取15 30 ；所述d为问题的维数，取为3 ；步骤4构造小波变换矩阵Qv 求出树结构中每个结点对应的积分点组ν所对 应的小波变换矩阵Qv，步骤如下步骤a:求最底层上的矩量矩阵[Mv] a,i= (Xi-Xv) a, a | Sp1，其中ν为任意第 L层格子，α为三重指标，XiSv上第i个积分点，格子中心为\；所述？1为1层上的小波 消失矩阶数；步骤b 对矩量矩阵Mv做奇异值分解Mv = /νΣνρν7'，得到ν上的小波变换矩阵 、尺 度函数变换矩阵众和尺度函数矩量矩阵M丨= /Λ，其中，4和众由Qv中分别与零奇异值和 非零奇异值对应的列组成，a =[QV,QV]，良包含Σ ν中与所有非零奇异值对应的列；步骤c 求高层上的矩量矩阵MV=[7^,…,UdiagiM^···,^^}，其中，ν为 第1层格子，1 = L-1，…，1，μ1，表示V的所有子格子，Tv, μ为矩转换矩阵， [^U 二 W- 广 Ca 0为二项式系数；步骤d 求虬的奇异值分解Mv = UνΣνρ丨，可得ν的小波变换矩阵众、尺度函数变换
矩阵α和尺度函数矩量矩阵似丨；步骤5计算边界元系数矩阵A 采用小波边界元法的非标准型离散方法，系数矩阵A的元素分布为<formula>formula see original document page 7</formula>其中，子矩阵i,、4和4(1 = L，...，1)；所述A1的计算步骤如下步骤I 按关系 Near(v) = {ν ‘ level (v) = level (ν ‘ ), dist(v, ν ‘) < nmax{diam(v), diam(v' )}}找出每个格子ν的所有临近格子集合Near (ν)，其中，η 为给定常数，dist(v，v')为ν和ν'的距离，diam(ν)为ν的特征尺寸，IeveKv)为ν所 在的层数；所述的η小于0.5;步骤II 对第L层上任意两个临近的格子ν和ν' e Near (ν)，计算矩阵
<formula>formula see original document page 7</formula>其中，Kij = K(xv,i;xv, ,j),K(x,y)为静电场积分方程的积分核函数，为ν中的 序号为i的积分点，nv为ν中积分点总数；步骤III 计算子矩阵足、為和不，其中的元素按如下方法计算Ql AvyQv, =Ivy,Ql AvyQv, =Ivy,QlAvyQvl=Tvy,QlAvyQv = ^v,V得到的I,、7V V,和J1^分别为矩阵入、足和劣中与格子ν和ν'对应的子块，入为
<formula>formula see original document page 7</formula>為和劣的结构与么相同；步骤IV 对第1层(1 = L-I，…，1)上任意两个临近的格子ν和ν'，计算矩阵
<formula>formula see original document page 8</formula>
其中，μ和μ ‘分别表示ν和ν'的子格子，所述的二迄^^么，其中，矩阵
Αμ,μ,采用递归方法计算当μ和μ ‘的层数为L时，计算方法采用步骤II ；当μ和μ ‘的层数为1 = L_l，…，1时，计算方法采用步骤IV ；步骤V 采用步骤III计算子矩阵為、為和孑(1 = L-I,…，1)，以及矩阵A1
<formula>formula see original document page 8</formula>其中，Vi和ν'」为第1层上的相邻格子；步骤6 设定线性方程组Ax = b，其中b的元素为飞行器边界曲面上每个积分点上 给定的电压或电荷密度的数值，A为边界元系数矩阵，χ为待求向量；所述待求向量χ为飞 行器边界曲面上待求的电压或电荷密度的数值；步骤7 采用迭代方法快速求解线性方程组Ax = b，得到χ。有益效果本发明提出的用快速边界元法得到大型复杂飞行器电场分布的方法，有益效果 是1、通过采用奈氏(Nystr0m)离散化方法，避免了数值积分的计算，简化了边界元
法的系数矩阵计算过程。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π ε 0|x-y|)，ε Q = 8. 8541878 为真空电容率，f (χ)为已知 函数，ο (y)为待求函数。应用飞行器电场分析的快速数值方法的主要步骤如下1.建立几何模型。可利用CAD软件(如Pro-E等)，按照实际尺寸建立飞行器边 界曲面模型，并将其保存为IGES格式，其中包含了每块曲面的参数方程，图2(a)显示了飞 行器边界各参数曲面的轮廓；2.布置积分点。从IGES模型文件中读取边界曲面参数方程，将每块曲面剖分为若 干曲面四边形子块(称之为“单元”)，保证每个单元有独立的参数方程，如图2(b)所示，总 单元数目为122。将单元参数方程转化到标准参考区域
X
上。根据计算精度要 求选取合适的高斯积分点数，并将定义在
X W，l]上的高斯积分点通过参数方程投 影到单元上，形成单元积分点。图3给出了飞行器翼面上一个单元、与7点高斯积分公式对 应的所有积分点的分布情况；3.积分点分组。本实施例是三维问题(d = 3)，可按以下步骤建立积分点的多尺 度分组(1)选取包含飞行器上所有积分点的正方体格子，如图4(a)所示，将其层数设为 0，显然此格子非空(包含积分点的格子是非空的)，将其作为八叉树的根结点；(2)层数1从0开始递增，如果某个第1层格子所包含的积分点数目大于某个给定 的数目Iitl = 30，则将其由各边中点等分为8个子格子，略去不包含积分点的格子，就得到了 第1+1层格子；(3)重复步骤(2)中的细分过程，直至每个最小格子中包含的积分点数目都不大于IV设最大层数为L。以上分组过程可由图4说明。图4(a)显示了第0层格子及其一次细分的情况，图 4(b)为第0层格子右-前-下方位的子格子的一次细分情况，图4(c)是由上述多尺度细分 所建立的八叉树结构，层数L = 2，它具有两个特征一是每个结点中包含一定数目的积分点，二是所有子格子中积分点之和就是父格子中的积分点。4.构造小波变换矩阵Qv。由于L层至第1层，逐层构造。取小波消失矩阶数与层 数的关系为P1 = L-1+1。如下两步构造(1)计算L = 2层()上的Qv。该层格子中包含的积分点数目均不大于η。= 30， 消失矩阶数为P2= 1。设该层格子ν中的积分点为Xl，x2，…，xn，如图5 (a)所示。计算矩 量矩阵Mv:
<formula>formula see original document page 10</formula>作奇异值分解从,=UvYvQl，得到ν的小波变换矩阵&、尺度函数变换矩阵众，同时 还得到ν的尺度函数的矩量矩阵M= = UvIv，存储Qv和Μνφ。(2)计算1层上的Qv。该层消失矩阶数为2。对每个格子V，设其子格子为μ ^·· μ s，如图5(b)所示。以步骤(1)为初始，每个子格子中已计算出尺度函数的矩量矩阵，即 Μη1φ,…，于是可由式Mv=[7；^,-,^JdiagiMj1,-,MjJ计算当前格子ν的矩量矩阵Mv。上式中转换矩阵！；，吣只与ν和μ的中心坐标有 关[H= Cf W-Xv 广"同样，计算Mv[，即可得ν的小波变换矩阵Qv和尺度函数的矩量矩阵 Mt=H通过以上两步，可对八叉树结构中除根结点外的每个结点构造出一个小波变换矩 阵Qv和尺度函数矩量矩阵Μνφ，它们将在下面的边界元系数矩阵计算中用到。 5.计算边界元系数矩阵Α，步骤如下(1)对每个非空格子V，计算临近格子集合NeaHv)，其中参数η = 0. 4 ；(2)对第2层上任意两个临近的格子ν和ν' e Near (ν)，计算矩阵Av, v,，其元素 为[Av, v, Jij = 1/(4 π ε J xvji-xv,,」|)其中，Xva为ν中的序号为i的积分点；(3)计算子矩阵Λ、為和孑，其中的元素按如下方法计算QlAyQv=Iy,Ql AvyQv, =Ivy,Q Av yQ, =Iy^QlAvyQv, = Ivy得到的、7"v v,和Zv v,分别为矩阵&、為和孑中与格子ν和ν'对应的子块。(4)对第1层上任意两个临近的格子ν和ν'，计算矩阵<formula>formula see original document page 11</formula>其中，μ和μ'分别表示ν和ν'的子格子，它们的层数为L = 2，I μ , μ ,的计算 方法为=Q7mA^Q,其中，矩阵Αμ,μ,按照第(2)步的方法计算，即[Αμ,μ, Jij = 1/(4 π ε0|χμ,「χμ, J)(5)按与第(3)步相同的方法计算子矩阵4 ,A1、^和A1。按上述步骤计算出的边界元系数矩阵具有如图6所示的稀疏结构，其中黑色区域 为非零元素，其余均为零元素。6.设定线性方程组Ax = b。计算飞行器电容时，在边界上施加IV电压，即积分方 程中f (X) = 1，于是，向量b的元素均为1。7.求解线性方程组Ax = b得到X，其中包含飞行器边界所有积分点上的电荷密度 σ (Xi)数值。飞行器电容为C = / s σ (χ) dx = Σ σ (Xi)式中，Wi为已知的积分权系数。
权利要求
一种用快速边界元法得到大型复杂飞行器电场分布的方法，其特征在于步骤如下步骤1按照飞行器的实际尺寸建立几何模型x＝γi(ξ)，i＝1，…，Ns，其中x为飞行器边界坐标，ξ为参数坐标，γi为第i块曲面的参数方程，Ns飞行器边界曲面个数，所述的飞行器边界曲面为光滑的；步骤2在每块光滑曲面上布置积分点将参数的取值范围ξ变换到标准参考区域
×
上，以参数方程x＝γi(ξ)将标准参考区域
×
上的Gauss积分点ξk(k＝1，…，ni)投影到飞行器的边界曲面上，得到xk＝γi(ξk)和与xk对应的Gauss权系数ωk；所述的Gauss积分点数ni大于10；步骤3积分点分组取包含飞行器边界曲面上所有积分点的正方体格子，对格子进行逐层细分，直至每个格子中包含的积分点数目小于给定常数n0，形成2d叉树结构；所述2d叉树根节点正方体格子的层数为0，最底层为最小的积分点分组、层序号为L；所述给定常数n0取15～30；所述d为问题的维数，取为3；步骤4构造小波变换矩阵Qv求出2d叉树结构中每个结点对应的积分点组v所对应的小波变换矩阵Qv，步骤如下步骤a求最底层上的矩量矩阵[Mv]α，i＝(xi-xv)α，|α|≤pl，其中v为任意第L层格子，α为三重指标，xi为v上第i个积分点，格子中心为xv；所述pl为l层上的小波消失矩阶数；步骤b对矩量矩阵Mv做奇异值分解得到v上的小波变换矩阵尺度函数变换矩阵和尺度函数矩量矩阵其中，和由Qv中分别与零奇异值和非零奇异值对应的列组成，包含∑v中与所有非零奇异值对应的列；步骤c求高层上的矩量矩阵其中，v为第l层格子，l＝L-1，…，1，μl，…μs表示v的所有子格子，Tμ为矩转换矩阵，Cαβ为二项式系数；步骤d求Mv的奇异值分解可得v的小波变换矩阵尺度函数变换矩阵和尺度函数矩量矩阵步骤5计算边界元系数矩阵A采用小波边界元法的非标准型离散方法，系数矩阵A的元素分布为其中，子矩阵和(l＝L，…，1)；所述A1的计算步骤如下步骤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)，计算矩阵其中，Kij＝K(xv，i，xv′，j)，K(x，y)为静电场积分方程的积分核函数，xv，i为v中的序号为i的积分点，nv为v中积分点总数；步骤III计算子矩阵和其中的元素按如下方法计算 <mrow><msubsup> <mover><mi>Q</mi><mo>^</mo> </mover> <mi>v</mi> <mi>T</mi></msubsup><msub> <mi>A</mi> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><msub> <mover><mi>Q</mi><mo>^</mo> </mover> <msup><mi>v</mi><mo>′</mo> </msup></msub><mo>=</mo><msub> <mover><mi>I</mi><mo>^</mo> </mover> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><mo>,</mo> </mrow> <mrow><msubsup> <mover><mi>Q</mi><mo>^</mo> </mover> <mi>v</mi> <mi>T</mi></msubsup><msub> <mi>A</mi> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><msub> <mover><mi>Q</mi><mo>~</mo> </mover> <msup><mi>v</mi><mo>′</mo> </msup></msub><mo>=</mo><msub> <mover><mi>I</mi><mo>~</mo> </mover> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><mo>,</mo> </mrow> <mrow><msubsup> <mover><mi>Q</mi><mo>~</mo> </mover> <mi>v</mi> <mi>T</mi></msubsup><msub> <mi>A</mi> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><msub> <mover><mi>Q</mi><mo>^</mo> </mover> <msup><mi>v</mi><mo>′</mo> </msup></msub><mo>=</mo><msub> <mover><mi>I</mi><mo>‾</mo> </mover> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><mo>,</mo> </mrow> <mrow><msubsup> <mover><mi>Q</mi><mo>~</mo> </mover> <mi>v</mi> <mi>T</mi></msubsup><msub> <mi>A</mi> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub><msub> <mover><mi>Q</mi><mo>~</mo> </mover> <msup><mi>v</mi><mo>′</mo> </msup></msub><mo>=</mo><msub> <mi>I</mi> <mrow><mi>v</mi><mo>,</mo><msup> <mi>v</mi> <mo>′</mo></msup> </mrow></msub> </mrow>得到的和分别为矩阵和中与格子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)，以及矩阵A1其中，vi和v′i为第1层上的相邻格子；步骤6设定线性方程组Ax＝b，其中b的元素为飞行器边界曲面上每个积分点上给定的电压或电荷密度的数值，A为边界元系数矩阵，x为待求向量；所述待求向量x为飞行器边界曲面上待求的电压或电荷密度的数值；步骤7采用迭代方法快速求解线性方程组Ax＝b，得到x。FSA00000075206200011.tif,FSA00000075206200012.tif,FSA00000075206200013.tif,FSA00000075206200014.tif,FSA00000075206200015.tif,FSA00000075206200016.tif,FSA00000075206200017.tif,FSA00000075206200018.tif,FSA00000075206200019.tif,FSA00000075206200021.tif,FSA00000075206200022.tif,FSA00000075206200023.tif,FSA00000075206200024.tif,FSA00000075206200025.tif,FSA00000075206200026.tif,FSA00000075206200027.tif,FSA00000075206200028.tif,FSA00000075206200029.tif,FSA00000075206200031.tif,FSA00000075206200032.tif,FSA00000075206200037.tif,FSA00000075206200038.tif,FSA00000075206200039.tif,FSA000000752062000310.tif,FSA000000752062000311.tif,FSA000000752062000312.tif,FSA000000752062000313.tif,FSA000000752062000314.tif,FSA000000752062000315.tif,FSA000000752062000316.tif,FSA000000752062000317.tif,FSA000000752062000318.tif,FSA000000752062000319.tif,FSA000000752062000320.tif
全文摘要
本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法，技术特征在于采用边界元法解决飞行器设计领域所存在的飞行器大规模静电场的快速计算分析问题。采用奈氏离散化方法，把介质界面划分成具有参数方程的光滑曲面，通过参数方程把积分点投影到界面上形成边界积分点；基于积分点的多尺度分组构造小波变换矩阵，对边界元系数矩阵进行矩阵压缩、计算和存储。克服了现有小波边界元法存在的系数矩阵计算复杂、计算量大、精度不高的缺点，有效降低了边界元分析的内存占用量和计算时间，使其可解决自由度上百万的大规模工程问题。本发明可应用于航空航天领域的大型复杂飞行器电场计算分析，以及其它工程领域的电场分析问题。
文档编号G06F17/50GK101833597SQ20101014225
公开日2010年9月15日 申请日期2010年4月8日 优先权日2010年4月8日
发明者文立华, 校金友 申请人:西北工业大学