本發(fā)明涉及電子電氣領域和計算機輔助工程的數(shù)值仿真方法領域，具體的說是一種計算靜磁場的正則化有限元數(shù)值方法，該技術可用于解決電子電氣設備中靜磁場計算的問題。

背景技術：

靜磁場數(shù)值仿真技術在電子電氣領域有重要的工業(yè)應用。作為靜磁場數(shù)值計算的一種常用方法，因為有效地克服了節(jié)點元在材料交界面出精度差問題，詳見文獻【jinjm.thefiniteelementmethodinelectromagnetics[m].johnwiley&sons,2014】，因此磁矢勢的棱邊單元法被廣泛使用，這種單元的一個特點是，未知自由度沿棱邊方向，大小為所求場量沿棱邊方向的積分。放松了節(jié)點元中未知量的法向連續(xù)性，只要求其切向分量連續(xù)，這被認為更符合物理定律，并可以得到更高的計算精度。在靜磁場磁矢勢控制方程中，常常需要施加額外的規(guī)范來保證磁矢勢的唯一性。節(jié)點元可以通過罰函數(shù)法在有限元列式中施加coulomb規(guī)范，以保證較好的矩陣性態(tài)，使迭代求解器較快地收斂。然而，棱邊元由于其本身的構造特點，不能如節(jié)點元那樣方便地施加庫侖規(guī)范即coulomb規(guī)范詳見文獻【bíróo.edgeelementformulationsofeddycurrentproblems[j].computermethodsinappliedmechanicsandengineering,1999,169(3):391-405】。

在現(xiàn)有的在棱邊單元法中，一種常見的做法是施加樹-余樹規(guī)范詳見文獻【magnetostaticswithedgeelements:anumericalinvestigationinthechoiceofthetree[j].magnetics,ieeetransactionson,1994,30(5):2877-2880】，這種方法在有限元網(wǎng)格的生成樹對應的棱邊上施加磁矢勢約束條件；由于生成樹對應的棱邊數(shù)正好等于矩陣的零特征值個數(shù)，因此通過在施加生成樹上的棱邊約束，就可以消除矩陣的奇異性；但是樹-余樹規(guī)范常常會使得所得到的有限元矩陣的條件數(shù)變大，導致迭代器收斂緩慢，甚至會帶來偽解。

另一種常見的做法是，無需施加規(guī)范，直接對離散方程采用迭代法求解，如共軛梯度法(cg)或不完全喬列斯基共軛梯度法(iccg)，進行求解詳見文獻【improvementofconvergencecharacteristicoficcgmethodforthea-methodusingedgeelements[j].magnetics,ieeetransactionson,1996,32(3):804-807】，但是此種方法的缺陷是為了保證迭代的收斂性，一般還要額外對激勵電流進行協(xié)調(diào)化處理，否則也會導致收斂緩慢或求解失??；此外，若激勵電流的形狀比較復雜，將難以對電流進行協(xié)調(diào)化處理。

最常用的方法是混合單元法，混合單元法作為處理含約束問題的一種有效方法，在電磁場求解中也得到應用，它可以消除棱邊元離散時雙旋度算子的零空間問題。它最早應用在力學問題中，例如不可壓固體和不可壓流體問題中。對于靜磁場問題，kikuchi均提出了相應的混合列式詳見文獻【mixedformulationsforfiniteelementanalysisofmagnetostaticandelectrostaticproblems[j].japanjournalofappliedmathematics,1989,6(2):209-221】。其中，散度約束條件通過引入lagrange乘子被施加到原問題中；并對磁矢勢或電場采用棱邊單元基函數(shù)插值，對標量lagrange乘子采用節(jié)點元基函數(shù)插值；然而，混合單元法離散得到的方程常是一個病態(tài)的鞍點問題。在這里所討論的靜磁場問題中，塊(1，1)處的k矩陣甚至還是高度奇異的。這無論對直接法還是迭代法的求解都帶來較大的挑戰(zhàn)詳見文獻【benzim,golubgh,liesenj.numericalsolutionofsaddlepointproblems[j].actanumerica,2005,14:1-137】。

技術實現(xiàn)要素：

鑒于已有技術存在的缺陷，本發(fā)明的目的是要提供一種計算靜磁場的正則化有限元數(shù)值方法，該方法通過采用正則化方法來求解靜磁場問題，從而具有較好的精度和效率。

為了實現(xiàn)上述目的，本發(fā)明的技術方案：

一種計算靜磁場的正則化有限元數(shù)值方法，其特征在于，包括如下步驟：

步驟1、建立與待分析的靜磁場問題相對應的磁矢勢方程，即采用磁矢勢方程來描述靜磁場，所述磁矢勢方程在ω中為：

curl(νcurla)＝js(1)

其中，a表示磁矢勢，js表示線圈中的激勵電流密度，ν表示磁阻率，ω表示求解區(qū)域，且在邊界γb上，對應的邊界條件：

n×a＝0(2a)

在邊界γh上，對應的邊界條件：

n×curla＝0(2b)

其中，n表示邊界的單位外法向量，同時為保證磁矢勢a的唯一性，在上述求解區(qū)域ω引入磁矢勢a的散度約束條件，所述散度約束條件選用coulomb規(guī)范：即在ω中滿足

diva＝0(3)；

步驟2、建立所述磁矢勢方程即公式(1)、邊界條件即公式(2a)及(2b)及約束條件即公式(3)，所對應的約束變分公式；

步驟3、對步驟2得到的所述約束變分公式，采用有限元方法進行數(shù)值離散，以得到所對應的離散的鞍點系統(tǒng)；

步驟4、將所對應的鞍點系統(tǒng)轉(zhuǎn)為一個無約束的最小化問題，并縮減掉所對應的鞍點系統(tǒng)的lagrange乘子自由度，得到一個對稱正定的系統(tǒng)，并采用預處理的共軛梯度法對所對應的對稱正定系統(tǒng)進行求解；

步驟5、對所獲得的計算結(jié)果進行后處理。

進一步的所述步驟2中具體包括如下步驟：

對公式(1)、(2a)以及(2b)，建立對應的自然變分公式，即使得式(1)、(2)的解對應泛函π取駐值：

δπ(a)＝∫ωνcurla·curl(δa)dω-∫ωδa·jsdω＝0(4)

式(4)中，δ為變分號；

對所述磁矢勢場函數(shù)施加約束條件，即引入lagrange乘子λ對所述磁矢勢場函數(shù)施加式(3)所示的約束條件，進而構造出所對應的修正泛函：

π^*＝π-∫ωλ·divadω(5)

且式(5)中所引入的lagrange乘子λ在邊界γb上滿足λ＝0的條件；并使得修正泛函π^*取駐值的條件為：

δπ+∫ωδa·gradλdω＝0(6a)

∫ωa·gradδλdω＝0(6b)。

進一步的所述步驟3中具體包括如下步驟：

對求解域ω剖分網(wǎng)格，并對磁矢勢a的磁矢勢場函數(shù)采用矢量基函數(shù)插值，對lagrange乘子λ采用標量基函數(shù)插值：

a＝w^tξ(7a)

λ＝n^tη(7b)

其中，w是磁矢勢a的矢量基函數(shù)，n是lagrange標量乘子λ的標量基函數(shù)，ξ表示磁矢勢a在求解域ω內(nèi)的棱邊自由度，η表示lagrange乘子λ在求解域ω內(nèi)的節(jié)點自由度；

將(7)式代入到(6a)、(6b)中，并令δa為矢量基函數(shù)w，δλ為標量基函數(shù)n，得到對應的有限元方程即所述離散的鞍點系統(tǒng)：

式中，k為雙旋度算子對應的離散矩陣，其表示為

k＝∫ωνcurlw·(curlw)^tdω；

g為散度算子對應的離散矩陣，其表示為

g＝∫ωw·(gradn)^tdω；

j為載荷向量，其表示為

j＝∫ωw·jsdω。

進一步的所述步驟4中具體包括如下步驟：

首先將求解(8)式所述的鞍點系統(tǒng)等效為下式(9)所示的含約束的最小化問題：

且滿足gξ＝0(9)

再將式(9)轉(zhuǎn)為相應的無約束最小化問題

并且式(10)中不含有l(wèi)agrange乘子自由度η，即lagrange乘子自由度被減縮掉；式(10)中，γ為正則化系數(shù)；使得式(10)的變分為0，得到以下的方程

[k+γg^tg]{ξ}＝{j}(11)

則當γ>0時，式(11)中的矩陣k+γg^tg是對稱正定的，即獲得所對應的對稱正定系統(tǒng)，并采用預處理共軛梯度法對所對應的對稱正定系統(tǒng)進行求解。

進一步的所述參數(shù)γ值確定方法包括：

式(12)中，hmin為所述剖分網(wǎng)格中的最小網(wǎng)格尺寸。

與現(xiàn)有技術相比，本發(fā)明的有益效果：

本發(fā)明具有較好的精度和效率，具體的其與傳統(tǒng)的樹-余樹規(guī)范棱邊單元法相比，其無需在有限元網(wǎng)格中生成樹即可保證解的唯一性，且所得到的矩陣有更好的性態(tài)，因此采用迭代求解器求解時，可以在較少的迭代步內(nèi)完成計算收斂，并且可以避免樹-余樹規(guī)范可能導致的偽解；其與無規(guī)范棱邊單元法相比，考慮到了非協(xié)調(diào)電流的影響，當存在非協(xié)調(diào)電流時，不會影響到迭代器收斂，因此本發(fā)明有更好的健壯性，適合處理任意的電流激勵。

附圖說明

圖1為所述塊狀腔體的有限元模型圖；

圖2(a)為所述矩陣k的譜分布圖；

圖2(b)為所述矩陣kr的譜分布圖；

圖2(c)為所述矩陣kt的譜分布圖；

圖3為所述不協(xié)調(diào)激勵電流源示意圖；

圖4為所述樹-余樹規(guī)范棱邊元方法得到的磁感應強度矢量圖的俯視圖；

圖5為所述正則化方法得到的磁感應強度矢量圖的俯視圖；

圖6為本發(fā)明所述方法對應的步驟流程圖。

具體實施方式

為使本發(fā)明的目的、技術方案和優(yōu)點更加清楚，下面將結(jié)合本發(fā)明實施例中的附圖，對本發(fā)明的技術方案進行清楚、完整地描述，顯然，所描述的實施例是本發(fā)明一部分實施例，而不是全部的實施例?；诒景l(fā)明中的實施例，本領域普通技術人員在沒有做出創(chuàng)造性勞動前提下所獲得的所有其他實施例，都屬于本發(fā)明保護的范圍。

本發(fā)明提出了一種采用正則化方法來求解靜磁場問題的方法，如圖6，其包括如下步驟：

步驟1、建立與待分析的靜磁場問題相對應的磁矢勢方程，即采用磁矢勢方程來描述靜磁場，所述磁矢勢方程在ω中為：

curl(νcurla)＝js(1)

其中，a表示磁矢勢，js表示線圈中的激勵電流密度，ν表示磁阻率，ω表示求解區(qū)域，相應的邊界條件為即在邊界γb上，對應的邊界條件：

n×a＝0(2a)

在邊界γh上，對應的邊界條件：

n×curla＝0(2b)

其中，n表示邊界的單位外法向量，同時為保證磁矢勢a的唯一性，在上述求解區(qū)域ω引入磁矢勢a的散度約束條件，所述散度約束條件選用coulomb規(guī)范：即在ω中滿足

diva＝0(3)；

步驟2、建立所述磁矢勢方程即公式(1)、邊界條件即公式(2a)及(2b)及約束條件即公式(3)，所對應的約束變分公式；進一步的，所述步驟2中具體包括如下步驟：對公式(1)、(2a)以及(2b)，基于自然變分原理建立對應的自然變分公式，即使得式(1)、(2a)以及(2b)的解對應泛函π取駐值：

δπ(a)＝∫ωνcurla·curl(δa)dω-∫ωδa·jsdω＝0(4)

式(4)中，δ為變分號；

對所述磁矢勢場函數(shù)施加約束條件，即引入lagrange乘子λ對所述磁矢勢場函數(shù)施加式(3)所示的約束條件，進而構造出所對應的修正泛函：

π^*＝π-∫ωλ·divadω(5)

且式(5)中所引入的lagrange乘子λ在邊界γb上滿足λ＝0的條件；于是，原泛函π的含約束的駐值問題就轉(zhuǎn)化為修正泛函π^*的無約束的駐值問題。π^*的駐值條件是它的一次變分為零，即使得修正泛函π^*取駐值的條件為：

δπ+∫ωδa·gradλdω＝0(6a)

∫ωa·gradδλdω＝0(6b)。

步驟3、對步驟2得到的所述約束變分公式，采用有限元方法進行數(shù)值離散，以得到所對應的離散的鞍點系統(tǒng)；進一步的，所述步驟3中具體包括如下步驟：對求解域ω剖分網(wǎng)格(如將求解域ω劃分成為四面體或六面體網(wǎng)格)，并對磁矢勢a的磁矢勢場函數(shù)采用矢量基函數(shù)插值，對lagrange乘子λ采用標量基函數(shù)插值：

a＝w^tξ(7a)

λ＝n^tη(7b)

其中，其中，w是磁矢勢a的矢量基函數(shù)，n是lagrange標量乘子λ的標量基函數(shù)，ξ表示磁矢勢a在求解域ω內(nèi)的棱邊自由度，η表示lagrange乘子λ在求解域ω內(nèi)的節(jié)點自由度；它們的表達形式參見文獻【jinjm.thefiniteelementmethodinelectromagnetics[m].johnwiley&sons,2014】；

將(7)式代入到(6a)、(6b)中，并令δa為矢量基函數(shù)w，δλ為標量基函數(shù)n，得到對應的有限元方程即所述離散的鞍點系統(tǒng)：

式中，k為雙旋度算子對應的離散矩陣，其表示為

k＝∫ωνcurlw·(curlw)^tdω；

g為散度算子對應的離散矩陣，且所述離散矩陣為滿秩矩陣，且ker(k)∩ker(g)＝{0}，

g＝∫ωw·(gradn)^tdω；

j為載荷向量，其表示為

j＝∫ωw·jsdω。

步驟4、步驟4、將所對應的鞍點系統(tǒng)轉(zhuǎn)為一個無約束的最小化問題，并縮減掉所對應的鞍點系統(tǒng)的lagrange乘子自由度，得到一個對稱正定的系統(tǒng)，并采用預處理的共軛梯度法對所對應的對稱正定系統(tǒng)進行求解；進一步的，所述步驟4中具體包括如下步驟：首先求解(8)式所述的鞍點系統(tǒng)等效為下式(9)所示的含約束的最小化問題即lagrange乘子自由度被減縮掉：

且滿足gξ＝0(9)

再將式(9)轉(zhuǎn)為相應的無約束最小化問題，且γ需要足夠大，確定參數(shù)γ值滿足

并且式(10)中不含有l(wèi)agrange乘子自由度η，即lagrange乘子自由度被減縮掉；||·||²表示該矩陣的2-范數(shù)，γ為正則化系數(shù)；使得式(10)的梯度為0，得到以下的方程

[k+γg^tg]{ξ}＝{j}(11)

則當γ>0時，式(11)中的矩陣k+γg^tg是對稱正定的，即獲得所對應的對稱正定系統(tǒng)，并采用預處理共軛梯度法對所對應的對稱正定系統(tǒng)進行求解。

更進一步的，原則上，我們可以采用預處理共軛梯度法來迭代求解線性方程(11)。但是，為了得到逼近原問題，我們需要設定γ→∞。但是，這會導致矩陣的條件數(shù)過大，出現(xiàn)病態(tài)。通過比較(9)和(1)-(3)，可以識別出約束方程gξ＝0與離散的coulomb規(guī)范對應，其中g可以看作離散的散度算子。在正則化式(11)中γ→∞表明coulomb規(guī)范條件是嚴格施加的，這使得(11)式等價于原來的鞍點問題。但是，顯然γ→∞導致矩陣病態(tài)，不能獲得預期的解。實際上，coulomb規(guī)范條件(9)并不是唯一的選擇。為了保證磁矢勢的唯一性，我們只需要保證磁矢勢的散度被定義即可。因此，從計算層面來看，離散的coulomb規(guī)范條件要求(13)式中的γ→∞，是沒有必要的，我們可以放松這個約束條件，選取一個有限大小的γ。這里，我們采用一種簡單確定γ大小的方法：

即所述參數(shù)γ值確定方法包括：

式(12)中，hmin為所述剖分網(wǎng)格中的最小網(wǎng)格尺寸，最后，我們得到一個對稱、正定的系統(tǒng)，從而可以采用預處理共軛梯度法進行迭代求解。

步驟5、對所獲得的計算結(jié)果進行后處理。

下面通過兩個案例來說明采用本發(fā)明方法的實施方式，或者通過與其它方法進行對比，說明采用本發(fā)明方法的優(yōu)勢。

案例1：案例描述：如圖1，考慮一個充滿均勻空氣的塊狀腔體，腔體的幾何尺寸是[0,1]×[0,0.5]×[0,0.75]。；采用有限元網(wǎng)格剖分，得到3543個四面體；其腔體的外表面定義γb邊界條件，施加磁矢勢為0；最后，在這個有限元模型中包含388個自由節(jié)點和3488個自由棱邊。

相應的優(yōu)勢說明：我們知道離散矩陣的譜性質(zhì)和krylov子空間方法的迭代收斂性之間有很重要的關系。一般來說，離散矩陣的特征值聚攏在單一的特征值附近是比較理想的，這時我們可以獲得較快的收斂效果；具體的，圖2是采用不同的規(guī)范方法得到的譜分布，其中，無規(guī)范棱邊單元法對應系數(shù)矩陣k，本發(fā)明所述正則化方法對應系數(shù)矩陣kr，以及樹-余樹規(guī)范棱邊單元法對應系數(shù)矩陣kt；如表1所示，我們可以看到矩陣k是半正定的，0特征值的個數(shù)等于自由節(jié)點的個數(shù)；當電流激勵的右端項不協(xié)調(diào)時，矩陣k的奇異性限制了無規(guī)范棱邊單元法的使用；矩陣kr和矩陣kt都是正定矩陣，但是矩陣kr的條件數(shù)要小于矩陣kt的條件數(shù)。

表1塊狀腔體的不同離散矩陣的譜性質(zhì)

案例2：案例描述：考慮正方體空氣盒[0,1]×[0,1]×[0,1]中放置一個長方體狀[0.4,0.6]×[0.4,0.6]×[0.2,0.8]的激勵電流的靜磁場問題；其中激勵電流密度的大小為1a/m²，沿著激勵電流區(qū)域的長邊方向，如圖3所示；正方體空氣盒的外邊界全部施加γb類型的邊界。由于該激勵電流在長方體區(qū)域的頂端和末端不能滿足散度為零的條件，因此該激勵電流在物理上是不協(xié)調(diào)性的。相應的實施步驟為：首先，建立上述靜磁場問題相對應的幾何模型，基于上述內(nèi)容，對幾何模型剖分網(wǎng)格，并設定磁導率、單元類型、邊界條件以及電流激勵等求解信息；基于本發(fā)明所述步驟2-5，獲得該模型的磁感應強度分布。

相應的優(yōu)勢說明：我們分別采用樹-余樹規(guī)范棱邊單元法和本發(fā)明所述方法，計算該問題。具體計算結(jié)果如圖4和圖5所示。我們可以發(fā)現(xiàn)，樹-余樹規(guī)范棱邊單元法計算出現(xiàn)了偽解，明顯不符合靜磁學規(guī)律，而本發(fā)明方法則得到了良好的結(jié)果。

以上所述，僅為本發(fā)明較佳的具體實施方式，但本發(fā)明的保護范圍并不局限于此，任何熟悉本技術領域的技術人員在本發(fā)明揭露的技術范圍內(nèi)，根據(jù)本發(fā)明的技術方案及其發(fā)明構思加以等同替換或改變，都應涵蓋在本發(fā)明的保護范圍之內(nèi)。