潘 超 何雯婷 郝文波 徐冰亮

（1. 東北電力大學(xué)電氣工程學(xué)院 吉林 132012 2. 黑龍江省電力科學(xué)研究院 哈爾濱 150030）

1 引言

變壓器直流偏磁條件下的時域磁場計算以穩(wěn)態(tài)電磁場分析為基礎(chǔ)[1-3]，通常直接求解有限元微分方程，數(shù)據(jù)運(yùn)算規(guī)模龐大，精確性較低[4,5]。文獻(xiàn)[6,7]采用諧波平衡有限元法對單元內(nèi)磁矢量位的各次諧波同時求解再疊加，結(jié)果具有較高的精確性；這種方法原理復(fù)雜，占用資源較多[8]。文獻(xiàn)[9]通過計算分解的方法來優(yōu)化資源，雖然減少了系統(tǒng)資源的占用率，但是計算時間較傳統(tǒng)方法更長?；陔娐?磁路耦合的方法能夠反映出變壓器直流偏磁時的勵磁特性，滿足工程分析需要，但是沒有充分考慮漏磁情況，導(dǎo)致準(zhǔn)確性降低[10]。場路耦合方法分為直接耦合與間接耦合兩種方式，目前的大部分研究都采用直接方法實(shí)現(xiàn)場與路的耦合，將磁場方程與電路方程聯(lián)立計算，效率較低[11]。上述方法大都建立二維模型模擬變壓器鐵心磁場分布，具有較高的精度，但二維場或?qū)ΨQ場僅可作為局部場域的近似處理，只能用于特定運(yùn)行方式或模型特性的分析，無法模擬變壓器直流偏磁條件下的內(nèi)部磁場分布[12]。這些問題在變壓器等電磁設(shè)備的時域電磁計算中尤為突出，并亟需解決。

另一方面，傳統(tǒng)時域計算通常采用單步低階方法，其截斷誤差較大，且由于時間迭代造成誤差累積，易引起求解振蕩。若細(xì)化空間網(wǎng)格或縮短步長，雖然可以在一定程度上減小振蕩，但計算規(guī)模會大幅度增加[13,14]。研究表明，變壓器瞬態(tài)磁場和時域場路計算中均存在計算穩(wěn)定性、效率及精確性的問題[15,16]。文獻(xiàn)[17]討論了變壓器時域場路計算過程中非線性勵磁、步長選擇及時間離散在直流擾動下的變化情況。但是關(guān)于變壓器時域場路耦合計算的穩(wěn)定性、精度及效率等優(yōu)化問題，目前罕有文獻(xiàn)研究。

本文結(jié)合變壓器直流偏磁問題，推導(dǎo)計算收斂時的穩(wěn)定域邊界，確定時間步長與計算穩(wěn)定性的關(guān)系，研究變壓器三維時域場路耦合計算效率與精確性的優(yōu)化方法。根據(jù)穩(wěn)定性分析提出一種自適應(yīng)變步長算法，在計算效率和精確性方面與定步長算法對比，并通過實(shí)驗(yàn)測量對其計算結(jié)果進(jìn)行驗(yàn)證。

2 時域場路耦合

利用時域場路耦合方法計算變壓器磁場和等效電路，動態(tài)電感和時域電流為關(guān)鍵耦合參數(shù)，分別由三維棱邊元磁場模型和等效電路模型計算。

建立變壓器磁場模型，不考慮磁滯效應(yīng)，棱邊有限元法采用矢量磁位A，得到非線性磁場方程

變壓器磁鏈方程

式中，ψ 為磁鏈向量；i為繞組電流向量；LS為靜態(tài)電感矩陣，表示激勵電流與磁鏈的穩(wěn)態(tài)關(guān)系。

實(shí)際上，變壓器勵磁過程中電感與電流均體現(xiàn)出時變非線性的特點(diǎn)[18]。由u=dψ /dt，推導(dǎo)變壓器等效電路的時域響應(yīng)方程

式中，u為交流電壓向量；UDC為直流電壓；LD為考慮勵磁變化的動態(tài)電感矩陣，需要通過磁場模型計算獲得。

基于能量擾動的思想，根據(jù)系統(tǒng)能量變化計算動態(tài)電感參數(shù)。當(dāng)線圈電流增加δip（0≤δ≤1，單相變壓器中 1≤p≤2）時，將電源總能量與動態(tài)電感和電流關(guān)聯(lián)，得到

若由電流增量δi引起的場量變化為δH、δB，變壓器內(nèi)部系統(tǒng)的磁場能量增量為

由能量平衡原理，方程（5）和（6）的能量相等，則可計算動態(tài)電感LD。

方程（4）的求解可采用四階龍格庫塔法，由tk時刻的線圈電流ik計算tk+1時刻ik+1

式中，h為步長；s1～s4為步長內(nèi)的分段計算斜率。

利用繞組電流i與動態(tài)電感LD即可實(shí)現(xiàn)時域場路耦合，其迭代計算原理如圖1所示。

式中，μ 為磁導(dǎo)率；J為電流密度，需要通過電路模型計算得到。

應(yīng)用格林定理，得伽遼金加權(quán)余量方程

式中，Mm為權(quán)函數(shù)序列，且權(quán)函數(shù)與基函數(shù)相同。

圖1 時域場路耦合Fig.1 Time-domain magnetic field and electrical circuit coupling

3 自適應(yīng)變步長算法

時域電路模型（4）的離散通解為式中，τ 為時間常數(shù)，取決于狀態(tài)矩陣的特征值λmax，且τ =1/λmax。

另一方面，在tk時刻式（4）中的泛函可化為常微分方程，得到式中，v為泰勒級數(shù)項(xiàng)的組合函數(shù)。

i(tk+1) 與 ik+1一致收斂到 h4階，聯(lián)立式（8）和式（10）有

式中，G為傳遞函數(shù)；χ =hλmax。

時域場路耦合計算收斂所對應(yīng)的絕對穩(wěn)定域判據(jù)為│G│＜1。確定時域場路耦合計算的穩(wěn)定邊界為

式中，h0為穩(wěn)定邊界對應(yīng)的計算步長。

對時域場路耦合迭代計算進(jìn)行改進(jìn)，提出一種自適應(yīng)變步長算法，其基本原理為

（1）時域場路耦合計算受初值、步長等因素的影響，在起始階段存在過渡過程，為了減少計算耗時，根據(jù)狀態(tài)矩陣特征值，在計算穩(wěn)定的前提下增大步長，從而提高效率。

（2）經(jīng)過某段計算時間直至計算穩(wěn)定時，在每次磁場求解動態(tài)電感之后，由變壓器狀態(tài)矩陣計算特征值，根據(jù)穩(wěn)定域邊界確定對應(yīng)的h0，引入修正因子α(0＜α＜1)，確定下一時刻的計算步長，h=αh0，進(jìn)行迭代計算。

（3）合理增大磁場計算的時間間隔，即降低磁場求解頻率，增加電路迭代計算次數(shù)，w=βh0（80≥w≥40，β取整數(shù)），從而保證計算效率與精確性。

定義tk時刻的動態(tài)電感變化函數(shù)g(LD)，表示在該時間段內(nèi)動態(tài)電感參數(shù)的變化速率，用來評價磁場計算的精確性。

式中，max(LD) 為動態(tài)電感的最大值；ε 為變化率限值。

4 算例分析

采用時域場路耦合模型計算變壓器直流偏磁問題，編譯四階龍格庫塔法程序求解時域電路模型，利用Ansys軟件建立變壓器的八分之一模型（模型尺寸與實(shí)際比例為1∶1），對絕對穩(wěn)定域邊界和自適應(yīng)算法進(jìn)行研究分析。實(shí)驗(yàn)變壓器型號為 BK300，如圖 2所示，具體參數(shù)見表 1；鐵心硅鋼片型號為DW360—50，磁化曲線如圖3所示。

圖2 實(shí)驗(yàn)變壓器Fig.2 Transformer for experiment

表1 變壓器參數(shù)Tab.1 Parameters of transformer

圖3 鐵心磁化曲線Fig.3 Magnetization curve of the core

4.1 時域場路耦合定步長算法

變壓器一次側(cè)空載電流峰值為Im，接入直流源UDC時的直流電流為 IDC。對空載運(yùn)行 IDC=0、IDC=25%Im、IDC=50%Im、IDC=100%Im時的直流偏磁情況進(jìn)行實(shí)驗(yàn)和仿真計算。無直流時，穩(wěn)定邊界對應(yīng)步長 h0=6.0×10-4s，步長h取 5.0×10-4s時計算收斂；隨著直流電流的升高，鐵心勵磁飽和程度加深，動態(tài)電感變化加劇，狀態(tài)矩陣特征值增大，穩(wěn)定性降低。當(dāng) IDC=50%Im時，h0=2.57×10-4s，h若取5.0×10-4s不滿足穩(wěn)定性要求；選取步長 h=2.50×10-4s，計算收斂。

分別采用基于時域場路耦合模型的定步長算法和自適應(yīng)變步長算法進(jìn)行計算，結(jié)果如圖4所示。

圖4 電流仿真計算和實(shí)驗(yàn)測量結(jié)果Fig.4 Results of currents of simulation and experiment

空載運(yùn)行時一次電流 i1即為勵磁電流 ie，由于磁場模型未考慮磁滯效應(yīng)，根據(jù) u→dψ/dt→dΦ/dt→B-H→i1的電磁耦合關(guān)系，i1為對稱波，如圖4a、圖4b所示，兩種方法計算的電流誤差較小。當(dāng)存在直流時，變壓器勵磁受直流水平影響，ie波形在正負(fù)半周不再對稱，隨著直流電流增大，變壓器勵磁飽和程度加深，ie波形畸變嚴(yán)重，正負(fù)半周的不對稱程度加劇。由圖4c可知，空載運(yùn)行直流偏磁時電流的計算結(jié)果與實(shí)驗(yàn)測量基本相同；兩者存在的誤差可能由磁滯所導(dǎo)致，在直流偏磁情況下兩者誤差更小，從而驗(yàn)證了本文所提方法是正確可行的。

設(shè)計低通濾波和時域差分模塊測量變壓器瞬時電感波形，并與計算的動態(tài)電感對比，如圖5所示。

圖5 動態(tài)電感Fig.5 Results of dynamic inductances of experiment and calculation

由圖5可以確定動態(tài)電感變化與勵磁非線性的對應(yīng)關(guān)系。無直流時，L1在勵磁的正負(fù)半周為對稱波形，L1的波峰、波谷分別對應(yīng)于勵磁的非飽和狀態(tài)與飽和狀態(tài)。ie接近零值時，鐵心勵磁處于不飽和區(qū)，L1數(shù)值趨于最大；當(dāng)ie趨于各半周內(nèi)的極值時，鐵心勵磁飽和程度逐漸加深，L1數(shù)值趨于最小。當(dāng)存在直流時，L1受直流水平影響，波形在正負(fù)半周不對稱，隨著直流電流增大，變壓器勵磁飽和程度加深，i1畸變嚴(yán)重，L1波形在正負(fù)半周的不對稱程度加劇。為了便于對比，計算時采用單調(diào)雙曲函數(shù)對磁化曲線進(jìn)行擬合，圖 5c的動態(tài)電感波形為固定步長計算結(jié)果，兩種算法計算的電感在整體波形上并沒有明顯差異，因此這里不再給出。

變壓器八分之一鐵心模型磁場計算結(jié)果如圖6所示。IDC=0時，鐵心最大磁感應(yīng)強(qiáng)度Bmax為1.24T，IDC=25%Im時，Bmax升至1.25T。通過計算可知，隨著直流水平升高，變壓器勵磁飽和程度加深，鐵心磁通密度增加。

圖6 鐵心磁場計算結(jié)果Fig.6 Magnetic distribution in iron core

采用定步長算法計算變壓器正常運(yùn)行，研究h、w對電感變化率和計算速度的影響。分析不同方案的計算穩(wěn)定性、精確性及效率，結(jié)果見表2。

表2 定步長不同方案的計算結(jié)果Tab.2 Results of cases with fixed-step algorithm

變壓器時域場路耦合計算過程中的電感參數(shù)具有時變性，導(dǎo)致勵磁飽和時對應(yīng)的穩(wěn)定域邊界較小。直流偏磁時隨著直流分量增大，變壓器勵磁飽和程度加劇，對計算穩(wěn)定性的要求也更為嚴(yán)格。當(dāng)?shù)嬎愕臅r間步長滿足穩(wěn)定域條件時，可以保證計算收斂，但是按該步長計算將占用較多資源和時間，效率較低?？偨Y(jié)其規(guī)律，計算的穩(wěn)定性由時間步長確定，滿足絕對穩(wěn)定域時，步長對計算效率和精確性的影響很小，此時效率和精確性主要取決于磁場的求解頻率。不難看出，定步長算法在計算變壓器非線性磁場時存在效率和精確性的優(yōu)化問題。

4.2 自適應(yīng)變步長算法

采用變步長算法進(jìn)行計算，并與定步長算法的效率和精確性對比，結(jié)果見表3。

表3 不同方法的計算結(jié)果Tab.3 Results of different algorithms

結(jié)果表明，兩種方法都能模擬變壓器直流偏磁時的勵磁過程。當(dāng)時間步長選定為0.25ms時，雖然能夠在一定程度上保證計算的穩(wěn)定性和精確性，但是整個過程將會耗費(fèi)大量時間。定步長算法雖然能夠模擬變壓器的勵磁過程，但不能充分考慮其隨時間變化的勵磁特性，因此在計算時可能存在不收斂的問題；同時由于該算法的h和w固定，導(dǎo)致其計算效率和精確性無法改善。自適應(yīng)變步長方法能夠模擬變壓器直流偏磁時勵磁的時變非線性，時域計算每次迭代的h都能自動滿足絕對穩(wěn)定域邊界，因此具有良好的計算穩(wěn)定性；在此基礎(chǔ)上合理選取 w能達(dá)到優(yōu)化計算效率和精確性的目的。

對于實(shí)際工程中普遍關(guān)注的大型電力變壓器直流偏磁問題，定步長算法可能受直流水平影響而不再適用，因而更能凸顯出自適應(yīng)優(yōu)化算法的計算效率和精確性。目前作者還未能獲取足夠的實(shí)際數(shù)據(jù)以進(jìn)行大型變壓器的直流偏磁計算，因此該內(nèi)容將是下一步研究的重點(diǎn)。

5 結(jié)論

結(jié)合變壓器直流偏磁問題，基于時域場路耦合計算的穩(wěn)定性分析提出一種自適應(yīng)變步長算法，通過研究得出以下結(jié)論：

（1）變壓器非線性勵磁具有時變性，直流偏磁時勵磁飽和程度加深，狀態(tài)矩陣特征值增大，對應(yīng)的穩(wěn)定域邊界變小。時域場路耦合計算時必須首先考慮穩(wěn)定性問題，其計算步長應(yīng)滿足收斂條件，在此基礎(chǔ)上合理選取磁場求解頻率能夠有效提高計算精確性和效率。

（2）與定步長時域場路耦合算法相比，自適應(yīng)變步長算法計算變壓器直流偏磁時能夠有效模擬內(nèi)部勵磁的時變非線性，具有良好的計算穩(wěn)定性，并能實(shí)現(xiàn)計算效率與精確性的合理優(yōu)化。空載直流偏磁的電流計算結(jié)果與實(shí)驗(yàn)基本相同，驗(yàn)證了該方法的正確性。

[1] C M Liu, L G Liu, R Pirjola. Geomagnetically induced currents in the high-voltage power grid in China[J]. IEEE Transactions on Power Delivery, 2010,24(4): 2368-2374.

[2] 郭滿生, 梅佳華, 張喜樂, 等. 直流偏磁條件下單相三柱電力變壓器的損耗計算[J]. 電工技術(shù)學(xué)報,2010, 25(7): 67-71.Guo Mansheng, Mei Jiahua, Zhang Xile, et al.Calculation of losses in three phase limb transformer under DC-biasing[J]. Transactions on China Electrotechnical Society, 2010, 25(7): 67-71.

[3] 趙志剛, 劉福貴, 張俊杰, 等. 直流偏磁條件下變壓器勵磁電流的實(shí)驗(yàn)與分析[J]. 電工技術(shù)學(xué)報,2010, 25(4): 71-76.Zhao Zhigang, Liu Fugui, Zhang Junjie, et al.Measurement and analysis of magnetizing current in DC-biased transformers[J]. Transactions on China Electrotechnical Society, 2010, 25(4): 71-76.

[4] Lu Shu, Liu Yilu. FEM analysis of DC saturation to assess transformer susceptibility to geomagnetically induced currents[J]. IEEE Transactions on Power Delivery, 1993, 8(3): 1367-1376.

[5] Mauricio Valencia Ferreira da Luz, Jean Vianei Leite,et al. Three-phase transformer modeling using a vector hysteresis model and including the eddy current and the anomalous losses[J]. IEEE Transactions on Magnetics, 2010, 46(8): 3201-3204.

[6] Zhao X J, Lu J W, et al. Analysis of DC bias phenomenon by the harmonic Balance Finite-Element method[J]. IEEE Transactions on Power Delivery,2011, 26(1): 475-485.

[7] Zhao X J, Li L, et al. Characteristic analysis of the square laminated core under dc-biased magnetization by the fixed-point harmonic-balanced FEM[J]. IEEE Transactions on Magnetics, 2012, 48(2): 747-750.

[8] Yamada R, Biringer P P, et al. Calculation of nonlinear eddy-current problems by the harmonic balance finite element method[J]. Transactions on Magnetics, 1991, 26(5): 4122-4125.

[9] 趙小軍. 基于諧波平衡有限元的變壓器直流偏磁特性研究[D]. 保定: 華北電力大學(xué), 2010.

[10] 李泓志, 崔翔, 盧鐵兵, 等. 變壓器直流偏磁的電路-磁路模型[J]. 中國電機(jī)工程學(xué)報, 2009, 29(27):119-125.Li Hongzhi, Cui Xiang, Lu Tiebing, et al. Electric circuit and magnetic circuit combined model of DC biased power transformer[J]. Proceedings of the CSEE, 2009, 29(27): 119-125.

[11] 郭健, 林鶴云, 徐子宏, 等. 單螺旋繞組變壓器支路電流的場路耦合計算及分析[J]. 電工技術(shù)學(xué)報,2010, 25(4): 65-70.Guo Jian, Lin Heyun, Xu Zihong, et al. Calculation and analysis of branch currents of single spiral winding transformer based on field circuit coupled method[J].Transactions on China Electrotechnical Society, 2010,25(4): 65-70.

[12] 潘超. 基于時域場路耦合模型的變壓器直流偏磁電磁特性研究[D]. 保定: 華北電力大學(xué), 2013.

[13] Shi Z W, Rajanathan C B. A method of approach to transient eddy current problems coupled with voltage sources[J]. IEEE Transactions on Magnetics, 1996,32(3): 1082-1085.

[14] D Rodger, N Allen. Calculation of transient 3D eddy currents in nonlinear media verification using a rotational test rig[J]. IEEE Transactions on Magnetics,1994, 30(5): 2988-2991.

[15] S L Ho, W N Fu. An incremental method for studying the steady state performance of induction motors using time stepping finite element model[J]. IEEE Transactions on Magnetics, 1997, 33(2): 1374-1377.

[16] S L Ho, S X Niu, et al. A power-balanced timestepping finite element method for transient magnetic field computation[J], IEEE Transactions on Magnetics,2012, 48(2): 291-294.

[17] 潘超, 王澤忠, 楊敬瑀, 等. 變壓器直流偏磁瞬態(tài)場路耦合計算穩(wěn)定性分析[J]. 電工技術(shù)學(xué)報, 2013,28(5): 82-87.Pan Chao, Wang Zezhong, Yang Jingyu, et al. Stability analysis based on transient magnetic-circuit coupled method for DC-biased transformer[J]. Transactions on China Electrotechnical Society, 2013, 28(5): 82-87.

[18] 王澤忠, 潘超, 周盛, 等. 基于棱邊有限元的變壓器場路耦合瞬態(tài)模型[J]. 電工技術(shù)學(xué)報, 2012,27(9): 146-152.Wang Zezhong, Pan Chao, Zhou Sheng, et al. Transient magnetic-circuit coupled model of transformer based on edge finite element method[J]. Transactions on China Electrotechnical Society, 2012, 27(9): 146-152.