一種考慮動態(tài)磁滯效應(yīng)的高效穩(wěn)定時域有限元計(jì)算方法

魏 鵬 陳 龍 賁 彤 井立兵 張 獻(xiàn)

（1.三峽大學(xué)電氣與新能源學(xué)院 宜昌 443002 2.湖北省微電網(wǎng)工程技術(shù)研究中心（三峽大學(xué)） 宜昌 443002 3.省部共建電工裝備可靠性與智能化國家重點(diǎn)實(shí)驗(yàn)室（河北工業(yè)大學(xué)） 天津 300130）

0 引言

時域有限元方法廣泛應(yīng)用于電工裝備的運(yùn)行狀態(tài)預(yù)測與損耗評估之中[1]，準(zhǔn)確預(yù)測電機(jī)、變壓器等電工裝備的時域暫態(tài)損耗分布特性不僅對降低設(shè)備能耗、改善局部過熱具有重要意義，同時考慮損耗對暫態(tài)磁場求解的反饋?zhàn)饔茫瑢M(jìn)一步分析電工裝備的瞬時功率平衡、暫態(tài)勵磁特征等電磁特性至關(guān)重要[2-4]。

目前，對于電工裝備鐵心損耗的有限元計(jì)算方法可大致分為兩類：一是基于場計(jì)算結(jié)果的后處理方法，該類方法利用Steinmetz 損耗模型或Bertotti統(tǒng)計(jì)學(xué)損耗分離理論對僅單值非線性磁場的計(jì)算結(jié)果進(jìn)行后處理，得到每一個單元的損耗，最終進(jìn)行積分得到電工裝備整體的損耗[5-7]。這一類方法主要用于對電機(jī)或變壓器的損耗進(jìn)行頻域分析，并通過傅里葉變換來分析含有諧波激勵時的鐵心損耗。Ansys 公司D.Lin 等基于等效橢圓環(huán)思想，將上述方法進(jìn)行了時域瞬時功率損耗計(jì)算擴(kuò)展[8]。該方法的優(yōu)勢在于在損耗計(jì)算過程中無需迭代，在磁路簡單的結(jié)構(gòu)中具有較好的精度，但由于沒有考慮鐵心磁滯、渦流損耗、剩余損耗等效動態(tài)磁場對鐵心中磁場分布的反饋?zhàn)饔茫陔姍C(jī)等復(fù)雜磁路電工裝備的損耗計(jì)算中將帶來較大誤差。另一類鐵心損耗計(jì)算方法則是在時域有限元分析中，考慮動態(tài)磁滯效應(yīng)對每一瞬時的功率損耗進(jìn)行精確求解[9-10]。這一類方法可最大限度地模擬鐵心的真實(shí)磁化過程，在處理含有磁通密度波形畸變工況下的損耗計(jì)算問題時具有較高精度。同時，在考慮鐵心動態(tài)磁滯效應(yīng)后，鐵心中動態(tài)磁滯回線波形的精確模擬可在時域中對勵磁電流的波形進(jìn)行實(shí)時反饋，進(jìn)而可對電機(jī)或變壓器中暫態(tài)運(yùn)行特性進(jìn)行有效預(yù)測。因此，電工裝備損耗的精確預(yù)測與運(yùn)行狀態(tài)的有效評估都離不開在時域有限元分析時對動態(tài)磁滯特性的影響的考慮。

然而，目前在有限元分析中常用的磁滯模型，如J-A 模型[11-13]、Preisach 模型[14-17]等均是與磁化速率無關(guān)的磁滯模型，并不能直接考慮電機(jī)或變壓器疊片宏觀渦流損耗與剩余損耗等動態(tài)因素對鐵心磁場分布的影響。為此，需要進(jìn)一步將二者的等效動態(tài)磁場考慮到有限元的計(jì)算方程之中。傳統(tǒng)方法是將動態(tài)效果等效為電導(dǎo)率乘以磁矢位A對時間的一階導(dǎo)數(shù)，但該種等效僅適合均勻?qū)嵭膮^(qū)域渦流路徑的求解，并不適合疊片鐵心結(jié)構(gòu)[18-19]。對于鐵心疊片結(jié)構(gòu)，進(jìn)行實(shí)際三維建模將會帶來無法承受的計(jì)算成本。為此，E.Dlala 等提出了一種耦合一維疊片有限元方法的二維時域有限元計(jì)算方法[20-21]。在疊片維度，忽略邊緣效應(yīng)，認(rèn)為磁通密度僅是疊片厚度方向上的函數(shù)，進(jìn)而僅需要求解一維有限元方程，即可得到疊片的渦流損耗等效磁場。然而，該方法實(shí)質(zhì)上是在二維有限元計(jì)算中嵌套了一維有限元程序，仍具有較高的求解難度。為了對上述求解方法進(jìn)行簡化，目前常采用Bertotti 模型的變形形式，即在耦合靜態(tài)磁滯模型的過程中，通過耦合等效微分磁阻率來考慮渦流損耗磁場與剩余損耗磁場的效果，該方法避免了一維有限元的求解，在一定程度上大大簡化了計(jì)算，但在每次迭代過程中仍需更新磁阻率和剛度矩陣，導(dǎo)致計(jì)算效率低下[22-23]。因此，在進(jìn)一步求解含有等效磁阻率的方程時，需要進(jìn)一步耦合非線性磁場計(jì)算中的固定點(diǎn)迭代技術(shù)。在僅考慮靜態(tài)磁滯模型的有限元計(jì)算中，固定點(diǎn)迭代方法已經(jīng)取得了較好的收斂效果與計(jì)算速度[24-26]，并已在目前的商業(yè)軟件中取得了較好的計(jì)算效果。然而，在引入渦流損耗磁場與剩余損耗磁場后，收斂將難以得到保證，特別是在計(jì)算至磁通密度零點(diǎn)附近時間步時算法極易發(fā)散，分析其原因主要是渦流損耗、剩余損耗等效動態(tài)磁場的介入改變了固定點(diǎn)法中磁通密度與磁場強(qiáng)度的收斂過程，使得傳統(tǒng)微分磁阻率的定義無法準(zhǔn)確描述收斂路徑的斜率，計(jì)算得到微分磁阻率過低導(dǎo)致考慮動態(tài)磁滯效應(yīng)的有限元方程求解不收斂。

為解決在考慮動態(tài)磁滯特性的時域有限元計(jì)算難以收斂的問題，本文提出一種基于等效磁阻率的固定點(diǎn)迭代求解策略，進(jìn)而提出一種耦合動態(tài)磁滯特性的時步有限元分析方法。首先，為了獲得電工鋼片的動態(tài)、靜態(tài)磁滯特性，基于愛潑斯坦方圈法搭建了一維磁特性測量系統(tǒng)；其次，考慮到逆Preisach模型可以考慮磁化歷史具有較高的計(jì)算精度，同時在耦合有限元計(jì)算時易于數(shù)值實(shí)現(xiàn)的特點(diǎn)，本文構(gòu)建了基于參數(shù)化解析一階回轉(zhuǎn)曲線的逆 Preisach模型，并進(jìn)行了參數(shù)辨識；再次，分析了渦流損耗和剩余損耗對算法收斂特性的影響，提出基于等效磁阻率的固定點(diǎn)迭代策略，并在有限元迭代過程中加入松弛因子保證穩(wěn)定性；最后，以愛潑斯坦方圈為例，進(jìn)行了考慮動態(tài)磁滯效應(yīng)的時域有限元分析，求解了瞬時功率損耗特性，并進(jìn)行了實(shí)驗(yàn)對比驗(yàn)證。結(jié)果驗(yàn)證了本文所提方法的準(zhǔn)確性、高效性與穩(wěn)定性。

1 準(zhǔn)靜態(tài)與動態(tài)磁滯特性測量

為建立磁滯模型及驗(yàn)證有限元計(jì)算準(zhǔn)確性，基于 IEC 60404—2 標(biāo)準(zhǔn)建立了一維磁性能測量平臺，對牌號為B30P105 的取向硅鋼片進(jìn)行準(zhǔn)靜態(tài)與動態(tài)磁滯特性測量。電工鋼片一維磁特性測量平臺如圖1 所示，該平臺包括寬頻功率放大器、電壓前置放大器、愛潑斯坦方圈、多功能數(shù)據(jù)采集卡和 LabVIEW 控制系統(tǒng)。本實(shí)驗(yàn)采用的硅鋼片規(guī)格及測試條件見表1。分別設(shè)定正弦交流電壓激勵的頻率為5 Hz 與50 Hz，測量得到準(zhǔn)靜態(tài)磁滯回線和50 Hz 動態(tài)磁滯回線，結(jié)果如圖2 和圖3 所示。

表1 硅鋼片規(guī)格及測試條件Tab.1 Specifications and test conditions of silicon steel sheets

圖1 電工鋼片一維磁特性測量平臺Fig.1 Measuring platform for 1-D magnetic propertiesof electrical steel sheets

圖2 5 Hz 準(zhǔn)靜態(tài)磁滯回線測量結(jié)果Fig.2 Measured quasi-static hysteresis loop (5 Hz)

圖3 50 Hz 動態(tài)磁滯回線測量結(jié)果Fig.3 Measured dynamic hysteresis loop (50 Hz)

2 靜態(tài)逆Preisach 磁滯模型的構(gòu)建

本文采用離散形式的靜態(tài)逆Preisach 磁滯模型描述硅鋼片的磁滯特性，該模型采用一種一階回轉(zhuǎn)曲線（First Order Reversal Curves, FORCs）的解析計(jì)算方法，并基于少量準(zhǔn)靜態(tài)磁滯回線測量數(shù)據(jù)實(shí)現(xiàn)解析式的參數(shù)辨識，從而構(gòu)造出辨識逆Everett 函數(shù)所需的一階回轉(zhuǎn)曲線，并進(jìn)一步建立靜態(tài)逆Preisach 磁滯模型。解析一階回轉(zhuǎn)曲線構(gòu)造示意圖如圖4 所示：一條起始于準(zhǔn)靜態(tài)極限磁滯回線下降支的一階回轉(zhuǎn)曲線R-P-N-T，該曲線起點(diǎn)為回轉(zhuǎn)點(diǎn)R，極限磁滯回線上升支Ha和下降支Hd交匯于頂點(diǎn)T。該一階回轉(zhuǎn)曲線構(gòu)造時需首先確定極限磁滯回線在回轉(zhuǎn)點(diǎn)R 處寬度ΔHrev和R 點(diǎn)與T 點(diǎn)的垂直高度ΔBrev為

圖4 解析一階回轉(zhuǎn)曲線構(gòu)造示意圖Fig.4 Construction ofanalytical first order reversal curve

圖4 中一階回轉(zhuǎn)曲線上點(diǎn)P 所在磁通密度為BP，磁場強(qiáng)度HP的計(jì)算式為

式中，ΔHout為極限磁滯回線在高度為BP處的寬度；x表征不同的P 點(diǎn)在一階回轉(zhuǎn)曲線上的相對位置，定義為

a、b、c為待擬合的模型參數(shù)。當(dāng)滿足a＞0，0＜b＜1，c＞0 等條件時，可使得一階回轉(zhuǎn)曲線斜率始終為正且不會超出最大極限磁滯回線。同時，以上模型參數(shù)與回轉(zhuǎn)點(diǎn)R 在下降支上的位置有關(guān)，此處定義位置系數(shù)β為

式中，ΔBout為最大環(huán)磁滯回線的高度。

參數(shù)a、b、c可表示為β的多項(xiàng)式，即

式中，y為多項(xiàng)式系數(shù)向量。

以上多項(xiàng)式系數(shù)通過準(zhǔn)靜態(tài)磁滯回線測量數(shù)據(jù)進(jìn)行辨識，由于對稱性，辨識過程僅需要磁滯回線的上升支。逆Everett 函數(shù)值可利用上升支一階回轉(zhuǎn)曲線Hforc通過式（7）計(jì)算。

磁通密度幅值為Bm的同心磁滯回線上升支可通過逆Everett 函數(shù)進(jìn)行插值運(yùn)算得到。

辨識程序中一共使用9 條磁滯回線，將每條磁滯回線上升支根據(jù)磁通密度均勻等分為50 個點(diǎn)?；谑剑?）、式（8）可得到由一階回轉(zhuǎn)曲線計(jì)算磁滯回線的數(shù)值關(guān)系，進(jìn)而建立適應(yīng)度函數(shù)F為

式中，Hmeas、Hcal分別為磁滯回線磁場強(qiáng)度的測量值和計(jì)算值；Bki為辨識程序中第k條磁滯回線上第i個點(diǎn)的磁通密度值。采用遺傳算法基于適應(yīng)度函數(shù)進(jìn)行參數(shù)尋優(yōu)，結(jié)果見表2。獲得全部參數(shù)后，采用該解析式模擬所需一階回轉(zhuǎn)曲線數(shù)據(jù)，進(jìn)一步可計(jì)算逆Everett 函數(shù)E(bu,bv)，通過解析一階回轉(zhuǎn)曲線辨識的逆Everett 函數(shù)如圖5 所示。

表2 B30P105 型電工鋼片多項(xiàng)式參數(shù)辨識結(jié)果Tab.2 The identified polynomial parameters of B30P105 elctrical steel sheet

圖5 通過解析一階回轉(zhuǎn)曲線辨識的逆Everett 函數(shù)Fig.5 Inverted Everett functionidentified by analytical first order reversal curves

通過上述逆Everett 函數(shù)辨識結(jié)果，可計(jì)算出一個周期內(nèi)不同磁通密度下電工鋼片的靜態(tài)磁滯回線。B30P105 型取向硅鋼的準(zhǔn)靜態(tài)磁滯回線實(shí)驗(yàn)結(jié)果和模型計(jì)算結(jié)果對比如圖6 所示，可以看出，本文所構(gòu)建的逆Preisach 磁滯模型具有一定的全局意義，可準(zhǔn)確模擬不同磁通密度幅值下的準(zhǔn)靜態(tài)磁滯回線，為進(jìn)一步考慮動態(tài)磁滯特性的有限元計(jì)算提供了可靠的模型基礎(chǔ)。

圖6 準(zhǔn)靜態(tài)磁滯回線測量值與計(jì)算值對比Fig.6 Comparison between measured and calculated quasi-static hysteresis loops

3 考慮動態(tài)磁滯特性的時域有限元方法

3.1 動態(tài)磁場的表征方法

Bertotti 基于磁疇理論和統(tǒng)計(jì)學(xué)原理提出損耗分離理論，將鐵磁材料單位體積的總損耗Ptot分為磁滯損耗Phys、渦流損耗Peddy和異常損耗Pexc三項(xiàng)。即

式中，σ為材料的電導(dǎo)率；d為硅鋼片厚度；S為橫截面積；G為無量綱耦合常數(shù)，G=0.135 6；V0為與材料內(nèi)部磁性單元有關(guān)的統(tǒng)計(jì)參數(shù)，同磁通密度幅值Bm相關(guān)，可通過不同頻率下的磁滯回線實(shí)驗(yàn)數(shù)據(jù)擬合得到。

磁場損耗W與B、H之間的關(guān)系為

將式（11）代入式（10）可得動態(tài)磁場表達(dá)式為

式中，Htot、Hhys、Heddy、Hexc分別為動態(tài)磁場總磁場強(qiáng)度、靜態(tài)磁滯磁場強(qiáng)度、渦流磁場強(qiáng)度、異常磁場強(qiáng)度；δB為符號函數(shù)，δB=sign(dB/dt)=±1。然而，由于較強(qiáng)的非線性特征，上述磁場難以直接在有限元分析中進(jìn)行迭代耦合計(jì)算，需進(jìn)一步研究其有限元迭代形式。

3.2 動態(tài)磁滯特性在有限元中的直接耦合形式

基于后向差分法，總磁場強(qiáng)度的三個分量的時域離散形式為

式中，t為時間步；Δt為步長間隔；vh為靜態(tài)磁滯微分磁阻率，定義為

相關(guān)麥克斯韋磁場方程為

結(jié)合方程式（13）～式（19）可得基于矢量磁位A的控制方程為

式中，ve為等效磁阻率，表達(dá)式為

3.3 基于等效磁阻率的固定點(diǎn)迭代策略

固定點(diǎn)法通過將非線性磁滯關(guān)系分為線性和非線性兩部分，通過在迭代過程內(nèi)不斷修改非線性部分以達(dá)到收斂，在解決磁滯非線性問題中具有較好的收斂性。為此，本文考慮將固定點(diǎn)迭代技術(shù)引入考慮動態(tài)磁滯特性的有限元迭代過程之中，在考慮靜態(tài)磁滯特性的基礎(chǔ)上，引入考慮動態(tài)磁場分量影響的等效磁阻率，采用式（21）中等效磁阻率代替?zhèn)鹘y(tǒng)固定點(diǎn)法的微分磁阻率vFP。應(yīng)用固定點(diǎn)技術(shù)，總磁場強(qiáng)度Htot和磁通密度B的非線性關(guān)系可表示為

式中，R為類磁化余量。

結(jié)合麥克斯韋方程，可得控制方程為

在二維求解域Ω中，將方程式（23）展開為離散形式為

式中，N為離散單元的基函數(shù)。將方程式（24）改寫為矩陣形式，即

式中，S為剛度矩陣；A為磁矢位向量；Y為電流區(qū)域系數(shù)向量；I為繞組電流向量；G為余量R的旋度矩陣。該方法中，渦流損耗等效磁場和剩余損耗等效磁場表達(dá)在類磁化余量R的計(jì)算中。通過這種方式，不需要大幅改動有限元方程，使得算法流程更簡潔，同時，將動態(tài)特性考慮在余量中，便于在以后的工作中對渦流和剩余損耗計(jì)算模型的改進(jìn)與修正。而且，在每一時刻的迭代過程，僅需計(jì)算一次等效磁阻率ve和剛度矩陣S，顯著提升了計(jì)算效率。

為保證固定點(diǎn)算法收斂，磁阻率v應(yīng)滿足式（26）的必要條件[25,27]。

在傳統(tǒng)固定點(diǎn)法中，所有網(wǎng)格單元采用統(tǒng)一的磁阻率，即

式中，vmin、vmax分別為磁滯回線上斜率最小值與最大值。該方法稱為GCM（global-coefficient method）。顯然，采用式（27）計(jì)算的磁阻率總能滿足算法收斂的必要條件。然而，GCM 采用統(tǒng)一的磁阻率導(dǎo)致計(jì)算時間過長，需采用更具效率的微分磁阻率進(jìn)行迭代，微分磁阻率vFP定義為

式中，dHtot、dB分別為磁場強(qiáng)度和磁通密度的差分。

對于標(biāo)量磁滯模型，二維電磁場量間的關(guān)系可通過取模值降階為一維問題。如圖7 所示，藍(lán)色實(shí)線表示靜態(tài)磁滯回線，藍(lán)色虛線表示動態(tài)磁滯曲線，曲線L 為固定點(diǎn)迭代過程中動態(tài)磁滯磁場Htot與磁通密度B的軌跡。由于磁場分量Heddy、Hexc的計(jì)算是基于磁通密度B的后向差分對時間的偏導(dǎo)數(shù)，因此，該曲線L 起始于t時刻靜態(tài)磁滯回線上的點(diǎn)P1，止于t+1 時刻動態(tài)磁滯回線上點(diǎn)P3。顯然，由于渦流損耗和剩余損耗的影響，曲線L 的斜率明顯不同于動態(tài)磁滯回線的斜率，在磁通密度的零點(diǎn)附近，兩者相差很大。傳統(tǒng)的固定點(diǎn)迭代法采用式（28）來計(jì)算微分磁阻率vFP，由于只考慮到動態(tài)磁滯磁場的軌跡，其計(jì)算結(jié)果偏小，往往不滿足收斂的必要條件，在迭代計(jì)算中算法極易發(fā)散。相反，在本文所提式（21）中等效磁阻率ve的計(jì)算考慮了渦流損耗和剩余損耗的等效磁場，反映了鐵心動態(tài)磁滯損耗對迭代過程的影響大小，可準(zhǔn)確估計(jì)動態(tài)B-Htot軌跡L 的斜率，從而保證迭代過程穩(wěn)定收斂。在有限元程序中ve根據(jù)前兩個時間點(diǎn)的結(jié)果進(jìn)行估算，即

圖7 動態(tài)磁滯收斂路徑示意圖Fig.7 Schematic diagramof dynamic hysteresis convergence path

3.4 考慮動態(tài)磁滯特性的場-路耦合有限元分析

鑒于電機(jī)、變壓器等大多數(shù)電氣設(shè)備工作于50 Hz正弦交流電壓激勵下，為反映動態(tài)磁場和電壓源的相互作用，需在時域有限元分析中進(jìn)行“場-路”耦合。由于愛潑斯坦方圈實(shí)質(zhì)可看作一臺等比的空載變壓器，可以滿足對變壓器鐵心損耗分析的需要。為驗(yàn)證本文所提方法的有效性，本文以愛潑斯坦方圈為例構(gòu)建了二維有限元計(jì)算模型。

在模型的計(jì)算過程中，以50 Hz 正弦交流電壓作為激勵，從電路角度分析，不考慮勵磁線圈漏磁和匝間電容，則端口電路方程為

式中，U為激勵電壓；Uin為繞組上的感應(yīng)電動勢；r為線路電阻；I為勵磁電流。對于一階三角形單元，Uin的計(jì)算式為

式中，Φ為繞組磁通；Nc為線圈匝數(shù)；Sc為繞組截面積；Ωc為電流區(qū)域；se為單元面積；lz為二維求解系統(tǒng)縱向深度；Ae為單元節(jié)點(diǎn)磁矢位。將式（31）代入式（30）得到

式中，U為激勵電壓向量；r為支路電阻矩陣。采用Crank-Nicolson 差分格式，將式（25）、式（32）聯(lián)立求解。

通過求解方程式（33）得到磁矢勢向量，然后根據(jù)B=?×A計(jì)算出每個單元的磁通密度大小。如果計(jì)算結(jié)果未達(dá)到收斂條件eB=∣Bi+1-Bi∣＜Bε（相對殘差閾值一般取εB= 1× 1 0-5），則通過磁滯模型計(jì)算靜態(tài)磁滯磁場Hhys、渦流磁場Heddy、剩余損耗磁場Hexc，進(jìn)而計(jì)算總磁場強(qiáng)度Htot，并通過松弛因子λ對Htot進(jìn)行修正。

式中，i為迭代次數(shù)。若磁通密度B相鄰迭代誤差滿足收斂條件，則判斷是否達(dá)到最大時間步，若t=tmax，則停止迭代，輸出全部計(jì)算結(jié)果，否則跳轉(zhuǎn)至下一時間步。整個有限元迭代過程如圖8 所示。

圖8 考慮動態(tài)磁滯特性的有限元分析流程Fig.8 Finite element iterative calculation process considering dynamic hysteresis characteristics

需要說明的是：不同于靜態(tài)磁滯磁場，動態(tài)磁場強(qiáng)度和磁通密度間的非線性關(guān)系不是一開始就確定的，由于渦流磁場和剩余損耗磁場不是直接通過有限元磁場方程表達(dá)，而是在固定點(diǎn)迭代過程加入，這意味著在算法迭代過程中，非線性磁場關(guān)系不斷改變，由動態(tài)磁場分量帶來的數(shù)值擾動使得磁矢位A和磁通密度B計(jì)算結(jié)果變化劇烈，造成收斂困難，尤其在磁通密度的零點(diǎn)附近，dB/dHtot較大，磁場強(qiáng)度的計(jì)算誤差引起的數(shù)值振蕩更大。因此松弛因子λ的引入是十分必要的。

4 結(jié)果與討論

4.1 愛潑斯坦方圈動態(tài)磁滯回線有限元模擬

為說明本文所提方法的有效性，在計(jì)算模型中，選取了愛潑斯坦方圈臂上的三角形單元作為磁滯回線與損耗的觀測點(diǎn)，其二維網(wǎng)格剖分情況如圖9 所示，三角網(wǎng)格M 為進(jìn)行實(shí)驗(yàn)數(shù)據(jù)與模型計(jì)算數(shù)據(jù)的對比的參考單元。在50 Hz 正弦條件下，圖10a、圖11a 分別給出了激勵電壓幅值為Umax=9 V 和Umax=13 V 時M 單元處磁滯回線測量值與計(jì)算值的對比?？梢钥闯?，再考慮動態(tài)磁滯特性后，本文所提方法具有較高的磁滯回線模擬精度。進(jìn)一步地，圖10b、圖11b 給出了相應(yīng)電壓激勵條件下愛潑斯坦方圈勵磁電流波形的模擬結(jié)果。由于考慮了動態(tài)磁場與磁滯效應(yīng)對電流波形的反饋?zhàn)饔?，本文所提方法可以較好地模擬勵磁電流時域波形。

圖9 愛潑斯坦方圈有限元剖分Fig.9 Finite element mesh of Epstein frame

圖10 Umax=9 V 電壓激勵下測量值與計(jì)算值對比Fig.10 Comparison between measured and calculated values under voltage excitation Umax=9 V

圖11 Umax=13 V 電壓激勵下測量值與計(jì)算值對比Fig.11 Comparison between measured and calculated values under voltage excitation Umax=13 V

為考察松弛因子λ對程序收斂性的影響，控制電壓激勵分別為Umax=9 V，調(diào)節(jié)松弛因子大小，記錄程序是否收斂，以及程序收斂情況下一周期內(nèi)的平均迭代次數(shù)和程序總計(jì)算時間。結(jié)果見表3，當(dāng)λ=0.8時，有限元迭代過程不收斂，在確保收斂的情況下，采用較大的松弛因子可減少每時步的平均迭代次數(shù)和總計(jì)算時間，提升計(jì)算效率，最短總計(jì)算時間為228 s。通常，為同時滿足穩(wěn)定性和計(jì)算速度的要求，設(shè)置λ=0.5～0.7 即可。

表3 9 V 電壓激勵時不同松弛因子對迭代次數(shù)和總迭代時間影響Tab.3 Effects of different relaxations on the iterations and computation time with voltage excitation amplitude of 9 V

4.2 磁通密度分布與瞬時損耗特性

為進(jìn)一步檢驗(yàn)算法的有效性，本文對愛潑斯坦方圈的磁通密度分布特性與瞬時功率損耗特性進(jìn)行了分析與計(jì)算。在t=0.005 s 時的磁通密度分布計(jì)算結(jié)果如圖12 所示?？梢钥闯?，主磁路上磁通密度分布較為均勻，側(cè)壁上磁通密度幅值計(jì)算結(jié)果與實(shí)驗(yàn)測量值基本一致；而在疊片的內(nèi)角處，磁通密度較大，在外角部分磁通密度較小。

圖12 磁通密度分布（Umax=13 V、t=0.005 s）Fig.12 Distribution of magnetic flux density (Umax=13 V、t=0.005 s)

根據(jù)4.1 節(jié)中磁通密度B和動態(tài)磁場Htot計(jì)算結(jié)果，本文進(jìn)一步計(jì)算了Umax=9 V 和Umax=13 V 兩種激勵大小在參考單元M 處一周期內(nèi)的瞬時損耗功率tp并與測量值進(jìn)行對比，結(jié)果如圖13 所示。在t=0.00 2 s 和t=0.012 s 附近，測量與計(jì)算的瞬時損耗達(dá)到極大值；在t=0.005 s 和t=0.015 s 附近，瞬時損耗達(dá)到極小值。損耗負(fù)值來源于可逆磁化的貢獻(xiàn)?？梢钥闯觯疚乃岱椒ㄔ谀M瞬時功率方面具有較高的精度，瞬時損耗的最大誤差不超過6%。

圖13 動態(tài)磁滯瞬時損耗測量值與計(jì)算值對比Fig.13 Comparison between measured and calculated instantaneous loss when considering dynamic hysteresis

完成一個周期有限元計(jì)算后，積分得到每一個單元磁滯回線的面積并求和即可得到鐵心區(qū)域整體鐵心總損耗PFe為

式中，f為激勵的頻率；ne為鐵心離散單元總數(shù)。圖14 所示為不同磁通密度幅值時的鐵心總損耗計(jì)算值和測量值大小?？梢钥闯觯ㄟ^提出的有限元算法計(jì)算得到的一個周期內(nèi)的鐵心損耗與實(shí)驗(yàn)測量結(jié)果基本一致。

圖14 鐵心損耗測量值與計(jì)算值對比Fig.14 Comparison between measured and calculated power losses

5 結(jié)論

為解決時域有限元分析中考慮鐵心材料的動態(tài)磁滯特性出現(xiàn)難以收斂的問題，本文提出了一種基于等效磁阻率的高效穩(wěn)定改進(jìn)固定點(diǎn)迭代策略，并編寫了相應(yīng)的有限元程序，實(shí)現(xiàn)了硅鋼片的動態(tài)磁滯特性進(jìn)行模擬計(jì)算，并通過實(shí)驗(yàn)測量結(jié)果對比，驗(yàn)證了該有限元算法的有效性與準(zhǔn)確性，具體貢獻(xiàn)如下：

1）采用解析方法生成一階回轉(zhuǎn)曲線，辨識逆Everett 函數(shù)，建立逆Preisach 磁滯模型，可準(zhǔn)確模擬電工鋼片靜態(tài)磁滯特性。

2）采用固定點(diǎn)迭代技術(shù)耦合動態(tài)磁滯特性，通過分析渦流損耗和剩余損耗等效磁場對有限元迭代過程的影響，提出基于等效磁阻率的固定點(diǎn)策略，并在有限元迭代過程引入松弛因子，通過實(shí)驗(yàn)數(shù)據(jù)與計(jì)算結(jié)果的對比分析，驗(yàn)證了本文提出的方法可實(shí)現(xiàn)低頻動態(tài)磁滯磁場的準(zhǔn)確計(jì)算，并具有良好的計(jì)算速度和穩(wěn)定性。

3）通過場-路耦合有限元算法與動態(tài)磁滯模型結(jié)合，可充分考慮鐵心總損耗對電路和磁場的反饋效應(yīng)，實(shí)現(xiàn)瞬時損耗的準(zhǔn)確計(jì)算，最大計(jì)算誤差不超過6%，可滿足工程計(jì)算要求。