CN101833597B - 用快速边界元法得到大型复杂飞行器电场分布的方法 - Google Patents
用快速边界元法得到大型复杂飞行器电场分布的方法 Download PDFInfo
- Publication number
- CN101833597B CN101833597B CN 201010142251 CN201010142251A CN101833597B CN 101833597 B CN101833597 B CN 101833597B CN 201010142251 CN201010142251 CN 201010142251 CN 201010142251 A CN201010142251 A CN 201010142251A CN 101833597 B CN101833597 B CN 101833597B
- Authority
- CN
- China
- Prior art keywords
- matrix
- grid
- aircraft
- point
- boundary element
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Complex Calculations (AREA)
Abstract
本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法,技术特征在于:采用边界元法解决飞行器设计领域所存在的飞行器大规模静电场的快速计算分析问题。采用奈氏离散化方法,把介质界面划分成具有参数方程的光滑曲面,通过参数方程把积分点投影到界面上形成边界积分点;基于积分点的多尺度分组构造小波变换矩阵,对边界元系数矩阵进行矩阵压缩、计算和存储。克服了现有小波边界元法存在的系数矩阵计算复杂、计算量大、精度不高的缺点,有效降低了边界元分析的内存占用量和计算时间,使其可解决自由度上百万的大规模工程问题。本发明可应用于航空航天领域的大型复杂飞行器电场计算分析,以及其它工程领域的电场分析问题。
Description
技术领域
本发明涉及一种用快速边界元法得到大型复杂飞行器电场分布的方法,飞行器设计领域,具体属于飞行器大规模电磁场的快速分析设计方法。
背景技术
在航空航天领域,飞行器结构的充、放电现象会干扰飞行器的通讯、导航,甚至造成飞行器结构破坏等严重后果。因此,对飞行器结构充、放电过程的计算分析是飞行器设计中必须考虑的问题,其核心是通过求解拉普拉斯方程,确定飞行器的静电场。这个问题具有求解区域复杂、求解规模大的特点。
现有飞行器静电场的计算方法大致有两类:
第一类是采用有限元法。由于飞行器结构的静电场本质上是一个无限域问题,采用有限元法时必须通过设置人工边界,得到包含飞行器的有限空间作为分析区域,计算出整个区域内的电学量(如电势、电荷密度等),而实际用到的只有飞行器表面的电学量。因此,这种方法的计算效率和精度较低,而且较费时。
第二类是采用常规边界元法,主要是求解定义在飞行器结构中不同介质界面上的拉普拉斯边界积分方程。与有限元法相比,该方法的自由度数目小、精度高。但常规边界元法形成稠密的系数矩阵,求解过程的内存占用量和计算时间与自由度数目N的平方成比例,即为O(N2)量级,无法应用到飞行器大规模静电场的快速计算分析中。
现有的降低常规边界元法内存占用量和计算时间的方法也可归为两类:一类是基于低秩逼近的,如快速多极方法、预修正快速傅里叶变换方法等,这类方法的特点是不显式计算系数矩阵,只是在迭代求解过程中对系数矩阵进行分块低秩分解,并完成矩阵-向量乘积运算,但它迭代一次的时间较长,当系数矩阵性态不好时,需要的迭代 次数较多,因此总的计算时间会很长。另一类是基于小波压缩的,这类方法在计算系数矩阵之前首先对其进行小波压缩,大大减少非零元素的数目,然后再显式计算并存储一个稀疏的系数矩阵,从而达到降低边界元法的内存占用量和计算时间的目的。与低秩逼近方法相比,小波压缩方法具有迭代时间短(大约小于前者的1%),系数矩阵预处理简便等优势。
但是,现有基于小波压缩的边界元法(以下简称“小波边界元法”)中都采用配点法和伽辽金(Galerkin)法作为离散化方法,主要存在三个问题:①需要计算单重或二重曲面积分,奇异积分处理复杂,系数矩阵计算量较大;②受边界剖分和形函数逼近方法的限制,难于得到高精度的结果;③飞机、航天器等的表面通常由大量光滑曲面片组成,划分单元会人为破坏边界的连续性,所以配点法和伽辽金法的计算效率较低。现有基于配点法和伽辽金法的小波边界元法所存在的系数矩阵计算方法复杂、计算量大、精度受单元剖分所限制、对复杂飞行器电场计算效率较低的不足。
发明内容
要解决的技术问题
为了避免现有技术的不足之处,本发明提出一种用快速边界元法得到大型复杂飞行器电场分布的方法,是一种基于奈氏 离散化方法的小波边界元法,并将其应用于大型复杂飞行器的电场计算分析中,以有效降低计算量、提高结果精度和求解效率。
技术方案
一种用快速边界元法得到大型复杂飞行器电场分布的方法,其特征在于步骤如下:
步骤1:按照飞行器的实际尺寸建立几何模型:x=Υi(ξ),i=1,…,Ns,其中x为飞行器边界坐标,ξ为参数坐标,Υi为第i块曲面的参数方程,Ns飞行器边界曲面个数, 所述的飞行器边界曲面为光滑的;
步骤2在每块光滑曲面上布置积分点:将参数的取值范围ξ变换到标准参考区域[0,1]×[0,1]上,以参数方程x=Υi(ξ)将标准参考区域[0,1]×[0,1]上的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层上的小波消失矩阶数;
步骤5计算边界元系数矩阵A:
采用小波边界元法的非标准型离散方法,系数矩阵A的元素分布为
所述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中积分点总数;
步骤IV:对第l层(l=L-1,…,1)上任意两个临近的格子v和v′,计算矩阵
其中,μ和μ′分别表示v和v′的子格子,所述的 其中,矩阵Aμ,μ′采用递归方法计算:
当μ和μ′的层数为L时,计算方法采用步骤II;
当μ和μ′的层数为l=L-1,…,1时,计算方法采用步骤IV;
其中,vi和v′j为第1层上的相邻格子;
步骤6:设定线性方程组Ax=b,其中b的元素为飞行器边界曲面上每个积分点上 给定的电压或电荷密度的数值,A为边界元系数矩阵,x为待求向量;所述待求向量x为飞行器边界曲面上待求的电压或电荷密度的数值;
步骤7:采用迭代方法快速求解线性方程组Ax=b,得到x。
有益效果
本发明提出的用快速边界元法得到大型复杂飞行器电场分布的方法,有益效果是:
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|),ε0=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层格子所包含的积分点数目大于某个给定的数目n0=30,则将其由各边中点等分为8个子格子,略去不包含积分点的格子,就得到了第l+1层格子;
(3)重复步骤(2)中的细分过程,直至每个最小格子中包含的积分点数目都不大于n0,设最大层数为L。
以上分组过程可由图4说明。图4(a)显示了第0层格子及其一次细分的情况,图4(b)为第0层格子右-前-下方位的子格子的一次细分情况,图4(c)是由上述多尺度细分所建立的八叉树结构,层数L=2,它具有两个特征:一是每个结点中包含一定数目的积分点,二是所有子格子中积分点之和就是父格子中的积分点。
4.构造小波变换矩阵Qv。由于L层至第1层,逐层构造。取小波消失矩阶数与层数的关系为pl=L-l+1。如下两步构造:
(1)计算L=2层()上的Qv。该层格子中包含的积分点数目均不大于n0=30,消失矩阶数为p2=1。设该层格子v中的积分点为x1,x2,…,xn,如图5(a)所示。计算矩量矩阵Mv:
(2)计算1层上的Qv。该层消失矩阶数为2。对每个格子v,设其子格子为μ1,…μs,如图5(b)所示。以步骤(1)为初始,每个子格子中已计算出尺度函数的矩量矩阵,即Mη1 φ,…,Mηs φ。于是可由式
计算当前格子v的矩量矩阵Mv。上式中转换矩阵Tv,μi只与v和μ的中心坐标有关:
通过以上两步,可对八叉树结构中除根结点外的每个结点构造出一个小波变换矩阵Qv和尺度函数矩量矩阵Mv φ,它们将在下面的边界元系数矩阵计算中用到。
5.计算边界元系数矩阵A,步骤如下:
(1)对每个非空格子v,计算临近格子集合Near(v),其中参数η=0.4;
(2)对第2层上任意两个临近的格子v和v′∈Near(v),计算矩阵Av,v′,其元素为
[Av,v′]ij=1/(4πε0|xv,i-xv′,j|)
其中,xv,i为v中的序号为i的积分点;
(4)对第1层上任意两个临近的格子v和v′,计算矩阵
其中,μ和μ′分别表示v和v′的子格子,它们的层数为L=2,Iμ,μ′的计算方法为
其中,矩阵Aμ,μ′按照第(2)步的方法计算,即
[Aμ,μ′]ij=1/(4πε0|xμ,i-xμ′,j|)
按上述步骤计算出的边界元系数矩阵具有如图6所示的稀疏结构,其中黑色区域为非零元素,其余均为零元素。
6.设定线性方程组Ax=b。计算飞行器电容时,在边界上施加1V电压,即积分方程中f(x)=1,于是,向量b的元素均为1。
7.求解线性方程组Ax=b得到x,其中包含飞行器边界所有积分点上的电荷密度σ(xi)数值。飞行器电容为
C=∫Sσ(x)dx=∑iwiσ(xi)
式中,wi为已知的积分权系数。
Claims (1)
1.一种用快速边界元法得到大型复杂飞行器电场分布的方法,其特征在于步骤如下:
步骤1:按照飞行器的实际尺寸建立几何模型:x=γi(ξ),i=1,…,Ns,其中x为飞行器边界坐标,ξ为参数坐标,γi为第i块曲面的参数方程,Ns飞行器边界曲面个数,所述的飞行器边界曲面为光滑的;
步骤2在每块光滑曲面上布置积分点:将参数ξ的取值范围变换到标准参考区域[0,1]×[0,1]上,以参数方程x=Υi(ξ)将标准参考区域[0,1]×[0,1]上的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层上的小波消失矩阶数;
步骤5计算边界元系数矩阵A:
采用小波边界元法的非标准型离散方法,系数矩阵A的元素分布为
所述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中积分点总数;
步骤IV:对第l层(l=L-1,…,1)上任意两个临近的格子v和v′,计算矩阵
当μ和μ′的层数为L时,计算方法采用步骤II;
当μ和μ′的层数为l=L-1,…,1时,计算方法采用步骤IV;
其中,vi和v′i为第1层上的相邻格子;
步骤6:设定线性方程组Ax=b,其中b的元素为飞行器边界曲面上每个积分点上给定的电压或电荷密度的数值,A为边界元系数矩阵,x为待求向量;所述待求向量x为飞行器边界曲面上待求的电压或电荷密度的数值;
步骤7:采用迭代方法快速求解线性方程组Ax=b,得到x。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201010142251 CN101833597B (zh) | 2010-04-08 | 2010-04-08 | 用快速边界元法得到大型复杂飞行器电场分布的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201010142251 CN101833597B (zh) | 2010-04-08 | 2010-04-08 | 用快速边界元法得到大型复杂飞行器电场分布的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN101833597A CN101833597A (zh) | 2010-09-15 |
CN101833597B true CN101833597B (zh) | 2013-02-06 |
Family
ID=42717665
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN 201010142251 Expired - Fee Related CN101833597B (zh) | 2010-04-08 | 2010-04-08 | 用快速边界元法得到大型复杂飞行器电场分布的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN101833597B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103258119A (zh) * | 2013-04-19 | 2013-08-21 | 东北电力大学 | 一种高压变电站工频电场的评估方法 |
CN105184010A (zh) * | 2015-09-25 | 2015-12-23 | 天津城建大学 | 基于快速多极间接边界元法的高频地震波散射模拟方法 |
CN108932392B (zh) * | 2018-07-13 | 2023-05-26 | 湖南科技大学 | 基于改进三重互易边界元法的瞬态温度计算方法 |
CN112632825B (zh) * | 2020-12-22 | 2023-03-10 | 重庆大学 | 一种基于有限元超收敛性的静电场光滑有限元数值算法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2009115457A (ja) * | 2007-11-01 | 2009-05-28 | Ntt Docomo Inc | 電界強度分布または電力密度分布の評価方法、そのプログラム |
CN101539599A (zh) * | 2009-04-14 | 2009-09-23 | 国网电力科学研究院 | 数字式雷电探测方法及其装置 |
-
2010
- 2010-04-08 CN CN 201010142251 patent/CN101833597B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2009115457A (ja) * | 2007-11-01 | 2009-05-28 | Ntt Docomo Inc | 電界強度分布または電力密度分布の評価方法、そのプログラム |
CN101539599A (zh) * | 2009-04-14 | 2009-09-23 | 国网电力科学研究院 | 数字式雷电探测方法及其装置 |
Non-Patent Citations (3)
Title |
---|
张智诠等.一种新的计算三维电位分布的表面电荷法.《兵工学报》.2009,第30卷(第10期),第1334-1338页. * |
文立华等.基于表面细分技术的边界谱方法.《西北工业大学学报》.2005,第23卷(第06期),第720-724页. * |
校金友等.准消失矩变阶小波Galerkin边界元法.《西北工业大学学报》.2009,第27卷(第6期),第786-790页. * |
Also Published As
Publication number | Publication date |
---|---|
CN101833597A (zh) | 2010-09-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102129523B (zh) | 基于mda和mlssm的分析复杂目标电磁散射的方法 | |
Nakahashi et al. | Building-cube method for large-scale, high resolution flow computations | |
CN102033985A (zh) | 基于*-矩阵算法的高效时域电磁仿真方法 | |
CN102156764B (zh) | 一种分析天线辐射和电磁散射的多分辨预条件方法 | |
CN102081690B (zh) | 复杂电路的矩阵分解结合新奇异值分解方法 | |
Lallemand et al. | A multigrid finite element method for solving the two-dimensional Euler equations | |
CN103617367A (zh) | 电磁场-流场-温度场耦合计算中的异型网格映射方法 | |
CN101833597B (zh) | 用快速边界元法得到大型复杂飞行器电场分布的方法 | |
CN109190169B (zh) | 一种三维时域电磁学杂交时域间断伽辽金数值方法 | |
CN106530397B (zh) | 一种基于稀疏剖面地质轮廓线的地质面三维重建方法 | |
Owen et al. | Parallel hex meshing from volume fractions | |
CN103440683A (zh) | 一种基于三维散乱稠密点云的三角网格重构方法 | |
Yu et al. | Enhanced QMM-BEM solver for three-dimensional multiple-dielectric capacitance extraction within the finite domain | |
CN102708229A (zh) | 复杂分层媒质结构的矩阵分解结合新奇异值分解方法 | |
CN111507026A (zh) | 一种模拟节点达西渗透流速的双重网格多尺度有限单元法 | |
CN104198820A (zh) | 一种含块状介质的双层土壤接地电阻的计算方法 | |
Luo et al. | Adaptive edge-based finite element schemes for the Euler and Navier-Stokes equations on unstructured grids | |
Karahan et al. | Time‐dependent groundwater modeling using spreadsheet | |
CN104778293B (zh) | 非均匀介质目标电磁散射的体积分Nystrom分析方法 | |
CN103914431A (zh) | 一种计算各向异性结构雷达横截面的无网格法 | |
Luo | A finite volume method based on weno reconstruction for compressible flows on hybrid grids | |
Kamenetskiy et al. | Numerical evidence of multiple solutions for the Reynolds-averaged Navier-Stokes equations for high-lift configurations | |
Kim et al. | High-density mesh flow computations by building-cube method | |
CN106294283A (zh) | 基于泰勒级数展开的时域积分方程快速算法 | |
Makino et al. | Aerodynamic analysis of NASA common research model by block-structured cartesian mesh |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20130206 Termination date: 20160408 |