用快速邊界元法得到大型複雜飛行器電場分布的方法

2023-05-07 20:19:51

專利名稱：用快速邊界元法得到大型複雜飛行器電場分布的方法
技術領域：
本發明涉及一種用快速邊界元法得到大型複雜飛行器電場分布的方法，飛行器設 計領域，具體屬於飛行器大規模電磁場的快速分析設計方法。
背景技術：
在航空航天領域，飛行器結構的充、放電現象會干擾飛行器的通訊、導航，甚至造成飛行器結構破壞等嚴重後果。因此，對飛行器結構充、放電過程的計算分析是飛行器設計 中必須考慮的問題，其核心是通過求解拉普拉斯方程，確定飛行器的靜電場。這個問題具有 求解區域複雜、求解規模大的特點。現有飛行器靜電場的計算方法大致有兩類第一類是採用有限元法。由於飛行器結構的靜電場本質上是一個無限域問題，採 用有限元法時必須通過設置人工邊界，得到包含飛行器的有限空間作為分析區域，計算出 整個區域內的電學量(如電勢、電荷密度等)，而實際用到的只有飛行器表面的電學量。因 此，這種方法的計算效率和精度較低，而且較費時。第二類是採用常規邊界元法，主要是求解定義在飛行器結構中不同介質界面上的 拉普拉斯邊界積分方程。與有限元法相比，該方法的自由度數目小、精度高。但常規邊界 元法形成稠密的係數矩陣，求解過程的內存佔用量和計算時間與自由度數目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 see original document page 7其中，子矩陣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 see original document page 7其中，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 see original document page 7為和劣的結構與麼相同；步驟IV 對第1層(1 = L-I，…，1)上任意兩個臨近的格子ν和ν'，計算矩陣
formula see original document page 8
其中，μ和μ 『分別表示ν和ν'的子格子，所述的二迄^^麼，其中，矩陣
Αμ,μ,採用遞歸方法計算當μ和μ 『的層數為L時，計算方法採用步驟II ；當μ和μ 『的層數為1 = L_l，…，1時，計算方法採用步驟IV ；步驟V 採用步驟III計算子矩陣為、為和孑(1 = L-I,…，1)，以及矩陣A1
formula see original document page 8其中，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 see original document page 10作奇異值分解從,=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 see original document page 11其中，μ和μ'分別表示ν和ν'的子格子，它們的層數為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計算子矩陣和其中的元素按如下方法計算 Q^ v T A v, v Q^ v = I^ v, v , Q^ v T A v, v Q~ v = I~ v, v , Q~ v T A v, v Q^ v = I v, v , Q~ v T A v, v Q~ v = I v, v 得到的和分別為矩陣和中與格子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日
發明者文立華, 校金友 申請人:西北工業大學