趙冬一，李 凡，湯 霖

（1.深圳ABB銀星避雷器有限公司，廣東深圳 518110；2.國家高低壓電器質(zhì)量監(jiān)督檢驗中心，甘肅天水 741018；3.中國電力科學研究院，武漢 4300741）

0 引言

從1967年日本松下電器開創(chuàng)MOA（metal-oxide surge arresters，MOA）技術(shù)的新紀元以來[1]，歷經(jīng)40多年的光輝歷程。MOA已經(jīng)成為電力系統(tǒng)絕緣配合和過電壓保護的主要裝置。特別是隨著特高壓交、直流輸電工程的建設(shè)，MOA技術(shù)的發(fā)展已經(jīng)到了一個新的水平[2-5]。數(shù)值仿真計算技術(shù)在MOA的研究發(fā)展中起到了十分重要的作用[6]。

基于MOA及其核心元件MOV（metal oxide varistors，MOV）的研究及實際工程技術(shù)數(shù)值計算的要求，需要二類仿真模型：1）限制暫態(tài)（或暫時）電壓時的計算模型。2）正常工作時的計算模型。因為MOV具有高非線性電阻特性以及其電阻、電容、介電常數(shù)都隨溫度、電場變化而變化，是一個非常復(fù)雜的多物理耦合場。自上個世紀80年代，國際上專家、學者開始用數(shù)值仿真計算技術(shù)研究MOA[7-8]。經(jīng)過近40年的發(fā)展，取得了很多成果。

筆者從3個方面對此進行綜述：1）基于等效電路模型的數(shù)值仿真計算；2）基于電、機械與熱的兩或多物理量（場）耦合的MOV數(shù)值仿真計算；3）整只MOA的電場、溫度場或多物理場耦合的數(shù)值仿真計算。

1 基于等效電路模型的MOA數(shù)值仿真計算

基于電力系統(tǒng)設(shè)計、分析的需要，當系統(tǒng)有非線性電阻的MOA時，通常在暫態(tài)網(wǎng)絡(luò)分析系統(tǒng)中加入其等效電路模型。常見的等效電路模型見圖1[9-12]。

圖1 不同用途的MOA等效電路模型Fig.1 Equivalent circuit model for different purposes of MOA

一般情況下，在預(yù)擊穿區(qū)域采用直流U-I的特性來校驗等效電路模型，而在翻轉(zhuǎn)區(qū)域，采用沖擊電流（例如：8/20 μs等）特性來校驗等效電路模型。這一特性可以用函數(shù)I=k Uα來進行模擬逼近，并且需要考慮電容CP和電感LS的頻響特性[13]。如果要求對數(shù)函數(shù)用于極性變化的雙極性模擬，那么進行計算時數(shù)值收斂可能難以接近電流零點。可以使用簡化的PSpice元件，如二極管和電壓源（圖1（c））[14]。學者們還研究了使用較復(fù)雜的等效電路模型去更準確地模擬MOV的特性，例如，對于低電流區(qū)域，可能需要考慮電阻的溫度系數(shù)以及電容的頻率系數(shù)，見圖1（d）[15]。

采用具有磁滯回線的非線性電阻模型可用于提高模擬計算的精度[16]。在U-I特性曲線的大電流翻轉(zhuǎn)區(qū)域內(nèi)，非線性電阻特性取決于沖擊電流波前的陡度，IEEE工作組3.4.11給出了等效電路模型[17]（以下簡稱IEEE模型），見圖2?；贗EEE模型的各種簡化模型也相繼出現(xiàn)，Pinceti[18]和Fernande[19]模型是其中兩種比較準確度高而且參數(shù)確定簡單的MOV動態(tài)模型。

圖2 IEEE公布的MOV等效電路模型Fig.2 IEEE published MOV equivalent circuit model

國內(nèi)對這一方面的研究也比較早、比較深入。西安交通大學張桂紅[20]等研究了MOV在各種沖擊電流下的等值電路模型，特別是采用圖2（a）陡波電流沖擊下的電路模型，指出IEEE工作組3.4.11給出了等效電路模型考慮了MOV的動態(tài)伏安特性。中國電力科學研究院的李志兵[21]等在分析IEEE推薦模型和東芝提出模型的基礎(chǔ)上，對安裝在GIS中的MOA在特快速過電壓下的相應(yīng)特性進行了計算，表明模型參數(shù)的取值對計算結(jié)果有很大的影響。最近上海交通大學何雨薇[22]等采用非線性電阻模型和3種常用動態(tài)模型（IEEE模型、Pinceti模型和Fernandez模型），估算不同用途的MOV在8/20、30/60和1 μs的陡波沖擊電流下殘壓和吸收能量，并與試驗測量結(jié)果對比分析。這是非常有益的嘗試。

2 基于電、機械與熱的兩或多物理量（場）耦合的MOV數(shù)值仿真計算

2.1 基于熱、機械耦合的MOV數(shù)值仿真計算

1984年，日本學者Eda[23]第一次系統(tǒng)建立了一個計算模型（MOV尺寸：?130 mm×10 mm），利用導(dǎo)熱微分方程和MOV的U-I特性方程，進行了有限差分法計算。

1999年，美國研究凝聚態(tài)物理的學者Bartkowiak、Mahan與Hubbell公司的Comber一起[24]，從陶瓷材料的熱應(yīng)力角度出發(fā)，對建立的MOV熱機模型進行計算，得出了不同波形和幅值下沖擊電流應(yīng)力下MOV的應(yīng)力變化情況，揭示了MOV能量耐受能力的熱機原理。見圖3[24]。

1997 年，美國固體材料學者Vojta、Clarke[25]等分析了一種新的MOV能量吸收下的失效模式。在極快速能量注入（高達：107k/s）的情況下，焦耳熱產(chǎn)生的慣性應(yīng)力會在預(yù)存微結(jié)構(gòu)下傳播和振蕩，應(yīng)力振蕩幅值是溫度對時間的二階導(dǎo)數(shù)，并形成正反饋，將導(dǎo)致電阻片折斷，見圖4[25]。

圖3 電站型MOA用MOV平均能量耐受能力與試驗電流峰值之間關(guān)系的仿真與試驗Fig.3 The simulated and mean measured energy handling capabilities of station-class-arrester disks as a function of the peak test current

圖4 大電流沖擊應(yīng)力下典型的“中平面”斷裂現(xiàn)象Fig.4 Typical“midplane”crack under high current impulse stress

2000年，加拿大學者Boggs與日本東芝公司的Andoh、Nishiwaki合作，在中國電磁暫態(tài)計算專家Kuang的幫助下，深入研究了電極對MOV能量耐受的影響，提出了優(yōu)化電極設(shè)計的技術(shù)思路，見圖5[26]。

2000年，澳大利亞Lengauer等[27]采用有限元分析法進行了3D仿真計算，見圖6[28]，并與文獻[28]比較，提出了MOV的設(shè)計要求，。

將來在這一方面的研究，應(yīng)主要針對熱應(yīng)力產(chǎn)生的彈性沖擊波進行仿真模擬計算。需要研究實際MOA設(shè)計中存在的電、熱、機械應(yīng)力的不均勻性、不同材料的熱特性、內(nèi)部芯體的固定壓緊接觸力以及MOA外套熱特性等，可能對運行中MOA壽命的影響。

2.2 二維和三維MOV導(dǎo)電機理的微觀模型計算

圖5 耦合熱電模擬得到的MOV的等溫等高線圖Fig.5 Coupled thermoelectric simulation of isothermal contour plots of MOV

圖6 電阻片的應(yīng)力波形圖（1D和3D）Fig.6 Stress oscillation in centre of the disc（1D and 3D）

利用模型仿真計算進行MOV導(dǎo)電機理的研究一直是專家、學者們關(guān)注的重點。借助數(shù)學工具，利用MOV的實驗數(shù)據(jù)建立的仿真模型可以讓人們看到實驗中不易想象的機理，對弄清楚MOV的電子水平導(dǎo)電機理起到了很大的作用。經(jīng)過反復(fù)、多次、長時間的試驗研究，人們總結(jié)出：對于多晶壓敏陶瓷材料的電輸運機理的實際模型，必須考慮材料的以下固有特征，1）MOV是由粒徑6～15 μm、雜亂無序的3D網(wǎng)絡(luò)多晶體結(jié)構(gòu)。見圖7[29]。2）每個晶界因為配方或其他工藝不一樣而具有局部的非線性U-I特征，并且是有方向的。3）每個晶界處的熱量產(chǎn)生以及由此產(chǎn)生的局部溫度分布都會影響局部晶粒和整體MOV溫度相關(guān)的非線性U-I特性。4）在工頻或沖擊電壓應(yīng)力的作用下，晶界空間電荷的介電響應(yīng)現(xiàn)象及其對微觀局部電壓分布有一定的影響。5）晶粒微結(jié)構(gòu)內(nèi)的不同晶相和晶界勢壘矢量方向產(chǎn)生的晶粒水平的微觀機械應(yīng)力[30]以及壓電效應(yīng)[31]。高非線性U-I特性、典型預(yù)擊穿區(qū)電壓負溫度系數(shù)（negative temperature coefficient，NTC）現(xiàn)象以及多晶體微觀結(jié)構(gòu)無序分布結(jié)構(gòu)結(jié)合，致使MOV在傳輸電流時呈現(xiàn)強烈的“細絲”流道，這可以從電致發(fā)光[32]或先進的紅外成像技術(shù)[33]觀察中直接觀察到，見圖 8[32-33]。

1996年，美國橡樹嶺國家實驗室學者Bartkowiak等[34]利用Voronoi網(wǎng)絡(luò)方法，對MOV的2D隨機分布非線性網(wǎng)絡(luò)進行了數(shù)值模擬，以便理解MOV多晶體陶瓷中的微結(jié)構(gòu)不規(guī)則性對整體非線性U-I特性的基本影響。由隨機分布的種子點產(chǎn)生的Voronoi細胞結(jié)構(gòu)見圖9[34]。Voronoi網(wǎng)絡(luò)方法計算的不同荷電率下電流導(dǎo)通路徑見圖10[34]。1997年，美國加州大學巴巴拉分校工程學院的Agnes Vojta等[35]也公布了類此網(wǎng)絡(luò)方法的研究成果。

圖7 MOV的SEM照片F(xiàn)ig.7 SEM photograph of MOV

圖8 MOV的電致發(fā)光和微觀紅外熱像Fig.8 SEM photograph of MOV

圖9 由隨機分布的種子點產(chǎn)生的Voronoi細胞結(jié)構(gòu)Fig.9 Voronoi cell structure generated from randomly distributed seed points

2002年，韓國Hanyang University的Se-Won HAN等[36]也利用Voronoi網(wǎng)絡(luò)方法計算MOV的電流密度和能量吸收能力，來模擬多晶體中的微觀結(jié)構(gòu)無序參數(shù)（例如晶粒尺寸和晶界）對擊穿失效的影響。2007年，美國學者Guogang Zhao等[37]建立了與時間相關(guān)的電熱隨機Voronoi網(wǎng)絡(luò)，來研究MOV響應(yīng)高幅值沖擊電流時的內(nèi)部加熱效應(yīng)和可能出現(xiàn)的故障機理，認為勢壘擊穿電壓的均勻分布和正態(tài)分布二者之間的差異很小。

圖10 Voronoi網(wǎng)絡(luò)方法計算的不同荷電率下電流導(dǎo)通路徑Fig.10 Voronoi network method to calculate the current conduction path under different charge rates

中國學者在這一方面也進行了大量的研究，得出了一些有益的研究成果。清華大學的何金良教授等[38-40]采用Voronoi網(wǎng)絡(luò)方法數(shù)值模擬，基于MOV的微觀電路和傳熱模型，分析了施加電流后因MOV的微觀結(jié)構(gòu)的不均勻性而引起的電流、熱場及應(yīng)力的不均勻性，解釋了能量吸收能力的不均勻性現(xiàn)象，見圖11[39]和圖12[39]。2010年，他們利用該方法研究了基于Voronoi網(wǎng)絡(luò)的、包含晶界的實際傳導(dǎo)機制以及晶界的電容的數(shù)值模型，仿真MOV的沖擊電流響應(yīng)與實驗結(jié)果波形非常吻合，見圖13[41]。

圖11 模擬電流定位現(xiàn)象（從白色到黑色的灰度級譜表示通過顆粒的電流的相對值，顆粒中的顏色越深，通過顆粒的電流越大）Fig.11 Simulated current localization phenomena（the graylevel spectrum from white to black represents the relative value of the current passing through a grain.The darker the color in a grain is，the MOV current passes through the grain）

圖12 MOV內(nèi)部的模擬溫度、熱應(yīng)力分布情況Fig.12 MOV inside the simulation temperature，thermal stress distribution

圖13 MOV顯微結(jié)構(gòu)內(nèi)的電流分布Fig.13 The current distribution inside the microstructure of the ZnO varistor

從圖中便可看出，盡管是二維數(shù)值仿真模擬計算，但是它使我們初步了解到影響MOV特性參數(shù)的主要因素，這些因素是很難從直覺預(yù)測到的。最近，奧地利學者Michael Hofstatter[42]、德國學者K.Bavelis[43]等將這種模擬已經(jīng)擴展到3D微結(jié)構(gòu)模型的研究中。

2014年，K.Bavelis通過3D微觀結(jié)構(gòu)模型和改進的等效電路方法研究MOV內(nèi)的電流傳輸機理。他們引入非線性等效電路模型表示微觀結(jié)構(gòu)中的電流通道的晶界和一定數(shù)值的電阻，見圖14[43]。建立的非均勻粒度分布的晶粒模型見圖15[43]，其圓柱形MOV模型的晶粒結(jié)構(gòu)無效晶界分布分別為5%、10%和20%，伏安曲線有較大的不同。他們進一步研究了在施加不同荷電率時，包含10%無效晶界分布的MOV內(nèi)的電流分布狀態(tài)，見圖16[43]。在預(yù)擊穿區(qū)域內(nèi)，觀察到不均勻的電流分布，見圖16（a）。此時，單個晶粒可以承載50%的總泄漏電流。在擊穿區(qū)域內(nèi)，電流變得更集中，見圖16（b），幾乎所有的電流都沿一條寬度不大于晶粒尺寸兩倍的狹窄路徑流動。隨著電壓增加到翻轉(zhuǎn)區(qū)域內(nèi)，電流分布又再次變得更均勻，見圖16（c）。該區(qū)域的MOV響應(yīng)本質(zhì)上是線性電阻，因此，正如預(yù)期的那樣，電流分布中非均勻粒度分布的影響不太明顯。

圖14 六角形晶粒的晶界電阻和電流分布和分布式晶粒電導(dǎo)表示方法的等效電路表示Fig.14 Grain boundary resistors for a hexagonal grain and current distribution and equivalent circuit representation in the distributed grain conductance approach

圖15 圓柱形MOV模型的晶粒結(jié)構(gòu)表示和分別為5%，10%和20%無效晶界分布計算的U-I特性曲線Fig.15 The grain structure representation of the cylindrical MOVmodelandthecalculatedU-I characteristiccurves for 5%，10%and 20%of the ineffective grain boundary distribution，respectively

圖16 包含10%無效晶界的圓柱形模型中的電流分布狀況Fig.16 Current distribution within a cylindrical sample

K.Bavelis比較了2D和3D模型。與2D方法相比，在3D模擬中，降低了晶界勢壘擊穿電壓和取消了U-I特性的簡化，導(dǎo)致2D與3D的電流路徑有很大的預(yù)期差異，特別是尺寸較小的樣品和考慮了比例較高的晶界勢壘不均勻性的情況下。他們認為，3D材料模型可以更詳細和更真實地分析MOV的電流分布狀態(tài)。遺憾的是國內(nèi)缺乏這方面的研究報道。

3 整只MOA的電場、熱場或多物理場耦合的數(shù)值仿真計算

3.1 整只MOA在持續(xù)運行電壓下的數(shù)值仿真計算

整只MOA在持續(xù)運行電壓工作條件下，由于其高介電常數(shù)和非常低的電導(dǎo)率而主要表現(xiàn)為電容性。類似線路絕緣子和多端口高壓斷路器，其電位分布受到了對地雜散電容（以及其他帶電部分）的影響。對于特高壓MOA、GIS用MOA及全屏蔽MOA尤其明顯。但是MOA與線路絕緣子不同，不均勻的電位分布會導(dǎo)致MOA熱的分布不均勻，特別是荷電率較高或多柱內(nèi)并聯(lián)的MOA更需要給予重視。至今為止，有關(guān)MOA標準（IEC、IEEE、JEC、GB）只考慮電位分布不均勻是不全面的。建議將來修訂標準時，應(yīng)規(guī)定MOA的溫度分布不均勻系數(shù)的上限值。這也應(yīng)是整只MOA在持續(xù)運行電壓下的數(shù)值計算模型努力發(fā)展的方向。

為了模擬電位或熱的不均勻分布，數(shù)值計算模型類型如下：

1）僅采用電性能方面的特征參數(shù)表示：a）純電容表示：最不均勻的電位分布情況，特別是MOA高壓端，如IEC60099-4:2014的附錄F；b）電容+非線性電阻（準靜態(tài)電場）表示：接近MOA實際電位分布。

2）采用電、熱耦合的特征參數(shù)表示：a）采用具有熱回路的等效集總參數(shù)電路來計算MOA熱或溫度分布；b）采用準靜電場與熱場相互耦合的模型類型：包括隨溫度、電場函數(shù)的電容/電阻效應(yīng)，時間維度上的功耗函數(shù)的電場和溫度效應(yīng)。

經(jīng)常用的MOA的U-I特性曲線，其電流值是電壓峰值瞬間的瞬時電流值，也即阻性電流的峰值。但是，在數(shù)值模擬計算中使用U-I特性曲線時，要求的是每時間步長電壓變化下的相應(yīng)電流值。應(yīng)該指出，到目前為止，建立在一個50 Hz半周期內(nèi)的隨場強變化而變化的MOA電導(dǎo)率數(shù)學計算模型和軟件程序還是困難重重。2010年，Chrsiten等[44]嘗試用受電場影響的電容來替代固定電容的解決方案；2015年，德國的Sébastien Blatt等[45]引入滯回特性來解決此問題，見圖17[45]。但此問題仍然懸而未決。

圖17 MOV的電位移密度-電場強度曲線Fig.17 MOV D-E curve

1989年，德國V.Hinriche等[46]通過基于集總元件電路和熱路的仿真研究了瓷外套MOA在工頻持續(xù)運行電壓下的電和熱行為，數(shù)學計算模型見圖18[46]?；趯嶒炑芯亢屠碚摻?，引入電和熱模型。將計算結(jié)果與試驗結(jié)果進行了比較。但是，詳細模擬計算過程沒有公布。

2005年，Zeller等[47]開發(fā)了一個FEMLAB數(shù)值模擬程序。他們建立了10 kV MOA的電加熱速率和靜止的傳熱耦合模型，用以計算其熱穩(wěn)定極限?？紤]到自由對流和輻射傳熱的模擬將不會收斂，簡化了模型，并將傳熱近似為線性或?qū)?shù)傳熱系數(shù)。由于MOV的非線性特性，他們設(shè)定了起始溫度并用MOA計算點處的近似行為來計算。該模擬方法能夠基于已建立的近似值預(yù)測給定MOA的熱穩(wěn)定極限，可以評估出熱穩(wěn)定極限數(shù)值大致范圍，以便進一步通過實驗來精確該值定位。但是，他們沒有考慮由于MOA的非線性U-I特性引起的函數(shù)收斂引起的動態(tài)過程，不過他們建議在計算點采用線性模型來改進這種模擬方法。

圖18 等效計算模型Fig.18 Equivalent calculation model

2010年，Sj?stedt等[48]采用基于準靜電場有限元法（FEM）方法3D模型，討論了UHV MOA的電壓分布設(shè)計。模型由純電容-電阻元件組成。計算過程為單步步長，沒有考慮耦合電熱行為。

2011年，ABB的Oliver Fritz等[49]基于COMSOL多物理場耦合模型，模擬研究了高壓GISMOA的熱行為。通過試驗數(shù)據(jù)擬合出MOV與溫度有關(guān)的非線性電導(dǎo)率，MOV的介電常數(shù)也以相同的方式確定，見圖19[48]。為了考慮從MOA柱到周圍氣體環(huán)境的熱傳遞，定義了MOA外表面的等效傳熱系數(shù)，并通過實驗數(shù)據(jù)來確定。他們指出了計算所需要處理的難點：對在ms時間尺度的電行為和分鐘、甚至小時時間尺度的熱行為的相互耦合進行建模。總模擬時間5 h，10 min通過重新計算電場并因此產(chǎn)生的焦耳熱來進行下一步的傳熱計算。因此，他們首先進行了瞬態(tài)熱量計算，再通過基于電場評估而重新計算的MOV的電阻率及功率損耗等。如此循環(huán)，最后分析評估整個MOA的溫度分布。Fritz等人表示：結(jié)果與試驗結(jié)果基本一致。與測量的結(jié)果相比，最大模擬溫度有些增加。這是由于保守的邊界條件導(dǎo)致的。用功率1 kW加熱MOA 5 h后的溫度分布見圖20[49]。

圖19 在6個荷電率下，測量和擬合出的單個MOV的功率密度Fig.19 Measured and fitted power loss density of a single varistor block for six different values of applied electric field

圖20 用功率1 kW加熱MOA 5 h后的溫度分布Fig.20 Temperature distribution after five hours with a constant heating power of 1 kW

2002年，在加拿大多倫多大學學習的中國學者鄭重與Steven A.Boggs[50]發(fā)表了他們在提高計算非線性電場效率方面研究的成果。引入“envelope”擬設(shè)準靜電場解算器，在網(wǎng)格內(nèi)確定非線性材料特性的方式，自適應(yīng)時間步進確保收斂，避免由時間步長之間的過度電荷轉(zhuǎn)移引起的大誤差以及通過場的適當空間和/或時間平均，在一階網(wǎng)格上，可以在合理的時間內(nèi)獲得準確的解決方案。鄭重等[51]依此研究了空腔結(jié)構(gòu)MOA的芯體調(diào)整金屬墊塊對散熱的影響。從芯體到外套壁為固定均勻的傳熱系數(shù)，而忽略熱引起的氣體對流局部傳熱的影響。按照IEC標準要求的動作負載試驗程序進行加載。改變傳熱系數(shù)和金屬墊塊的位置。鄭重等假設(shè)實際持續(xù)運行電壓后MOA很快達到熱平衡。他們重點關(guān)注了金屬墊塊對MOA的散熱效能。金屬墊塊在降低功率耗散和改善MOV的熱傳輸方面起著重要作用。2010年，已經(jīng)到華北電力大學工作的鄭重教授與Steven A.Boggs、東芝公司的Toshiya Imai和Susumu Nishiwaki等[52]繼續(xù)進行了金屬墊塊對MOA散熱影響的研究，重點關(guān)注復(fù)合外套MOA的散熱機理。金屬墊塊配置對MOA散熱影響計算見圖21[52]。因此，為了表征樣本的傳熱參數(shù)，擬合了單獨熱冷卻過程的測量數(shù)據(jù)。他們得出結(jié)論，熱參數(shù)以及散熱器元件的位置有助于熱穩(wěn)定性。這兩種效應(yīng)可能無法在保守的縮小比例模型中正確表示。遺憾的是，鄭重等沒有對仿真計算結(jié)果進行試驗驗證。

圖21 金屬墊塊配置對MOA散熱影響計算圖Fig.21 Metal pad configuration on the MOA cooling effect calculation diagram

2006年，Clemens等[53]首先嘗試一種單步計算方法。2008年，Hinrichsen等[54]進一步改進了單步方法，在MEQSICO的軟件包進行實施，用以計算MOA的電壓和溫度分布。首先，用Galerkin FEM公式，基于準靜電場和熱傳導(dǎo)方程，建立了電、熱耦合計算模型。模型分別以純電容C和可變電阻R表示，成功計算了IE60099-4:2014規(guī)定的動作負載試驗程序。隨后，對其溫度的分布進行研究。圖22中[54]，（a）為空外套值；（b）為電容為60 pF；（c）為 電容為10 nF。遺憾的是沒有進行實際測試并進行比較。這種方法與標準IEC 60099-4—2014中提出的兩步法相比，在計算時間、實用性和熱耦合效應(yīng)等方面，具有一定的優(yōu)點。該程序需要進一步發(fā)展和驗證。

圖22 MOA的溫度、電壓分布圖Fig.22 MOA temperature and voltage map

2011年，山東大學溫惠等[55]在計算串聯(lián)補償電容器組保護用MOV時，建立了3D溫度場有限元模型，利用ANASYSY軟件計算了不同散熱結(jié)構(gòu)效果，為MOV的本體結(jié)構(gòu)設(shè)計提供了理論依據(jù)。2014年，德國Darmstadt大學的 Sp?ck Leigsnering發(fā)表了碩士論文《Modeling and Simulation of a Station Class Surge Arrester》[56]，這是一篇很重要的文獻！

Sp?ck Leigsnering采用FEM數(shù)學計算方法，利用準靜電場和瞬態(tài)熱傳導(dǎo)方程的耦合計算模型，采用多速率時間步長方案的策略，分析了典型的電站型550 kV瓷外套MOA，見圖23[56]。模型考慮了非線性等效材料的輻射和自然對流，使MOA內(nèi)部氣腔的傳熱模型得到了進一步發(fā)展。根據(jù)IEC操作動作負載試驗的要求，他們計算了試樣在連續(xù)注入沖擊能量及后續(xù)的暫時過電壓下的應(yīng)變。其中瓷外套表面的熱傳遞采用了一種假設(shè)邊界條件來處理。計算結(jié)果得到了實測數(shù)的驗證。

圖23 被計算的MOA情況Fig.23 Calculated MOA situation

3.2 整只MOA在限制暫態(tài)過電壓的數(shù)值仿真計算

對于在計算限制暫態(tài)過電壓時的完整MOA整只數(shù)學模型，常見的方法是將MOA看作集總元件電路。IEC、IEEE提供許多不同的模型。這些模型主要由電容C，電感L，可變電阻R，恒壓源及可能有二極管組成，如圖1、2所示。當在電力系統(tǒng)絕緣配合的研究中，MOA對緩波前沖擊過電壓（slow-front overvoltages，SFO），快波前過電壓（fast-front overvoltages，F(xiàn)FO）和特快速波前過電壓（very-fast-front overvoltages，VFFO）的響應(yīng)是不同的。對于不同的標準沖擊電流波形，如 30/60 μs，8/20 μs，4/10 μs，1/10 μs或工頻正弦半波，需要用不同的U-I特性來表示，見圖24[43]。MOA的保護水平受行波特性的影響很大，特別是在FFO和VFFO下。在一些狀態(tài)下特別是低電壓系統(tǒng)中，MOA的高壓引線和接地引線的電感效應(yīng)是主要因素，保護水平主要取決于安裝方式而非MOA保護水平。

圖24 不同沖擊電流波形下的MOA的U-I特性曲線Fig.24 Current distribution within a cylindrical sample

當采用數(shù)值模擬計算MOA的能量吸收能力時，單或者重復(fù)脈沖能量吸收能力和熱吸收能力是有區(qū)別的。關(guān)于這一方面國際上的研究進展，筆者在文獻[57]中進行了詳細的綜述，不再贅述。我們需要重點關(guān)注，IEC TC37和Cigre Working Group A3.17在德國達姆施塔特理工大學，由Hinrichsen教授主持的一系列研究成果[58-61]。這些成果得到了國際上的認可，并在IEC60099-4:2014[62]中得到了體現(xiàn)，即，重復(fù)電荷轉(zhuǎn)移額定值Qrs和熱電荷轉(zhuǎn)移額定值Qth（在配電MOA的情況下）或熱能額定值Wth（在站MOA的情況下）等定義和術(shù)語。

MOA受到脈沖后的溫度升高數(shù)值，可以通過單位體積吸收的能量（一般是TNA研究的輸出結(jié)果）和其比熱容Cth進行計算。通常，Cth≈2.59 J/（cm3·k），具有0.0044 J/（cm3·K2）的溫度系數(shù)[63]。

2011年，中國電力科學研究院賀子鳴等[64]利用有限元法建立了多柱并聯(lián)結(jié)構(gòu)MOA散熱3D計算模型，采用2 ms方波電流試驗方法在較短時間內(nèi)給芯體注入較大能量以模擬系統(tǒng)MOA的暫態(tài)溫升過程。計算與試驗結(jié)果對比表明，誤差＜5%。

Sp?ck Leigsnering在文獻[56]也給出了這一方面的研究成果，是在1臺普通四核PC（2.4 GHz）上花了一周的時間。我們在3.3中給予綜述。

國內(nèi)專家、學者在這一方面也進行了大量的研究。但是，還處于探索階段[65-67]。

3.3 整只MOA電場、熱場及其耦合的數(shù)值仿真計算梳理

從上述文獻的研究思路來看，對于MOA的數(shù)學仿真計算的建模，要從3個方面分別進行考慮：1）利用靜電場或準靜電場，簡化麥克斯韋（Maxwell）電磁方程組[68]；2）采用非線性等效材料在MOA氣腔內(nèi)傳熱近似建立MOA的熱模型，包括輻射和自然對流；3）建立電、熱耦合數(shù)學計算模型。它在考慮電、熱不同時間尺度的情況下，可以通過多步長時間積分技術(shù)有效解決。

3.3.1 利用靜電場或準靜電場，簡化麥克斯韋（Maxwell）電磁方程組

MOA在正常持續(xù)運行電壓下，表現(xiàn)為阻容特性。MOA的材料特征參數(shù)為介電常數(shù)ε（r）和電導(dǎo)率σ=σ（r，E，T）。一般情況下，1 m長的MOV芯體的在MOV的預(yù)擊穿區(qū)域，σ=0.2 μs/m；在翻轉(zhuǎn)區(qū)域，σ=0.01 μs/m。則

因為τe＞τem＞τm，所以可采用準靜電場進行計算。要注意到，假設(shè)沖擊頻率在幾MHz范圍內(nèi)時，這種數(shù)學關(guān)系完全適用于MOV的非線性U-I特性。因此，準靜電場的方法可用于預(yù)擊穿區(qū)、擊穿區(qū)和翻轉(zhuǎn)區(qū)的MOV所有工作點。麥克斯韋（Maxwell）電磁方程組中的電荷守恒方程描述了準靜電場。從安培定律開始，div運算符應(yīng)用于方程的兩邊，用本構(gòu)方程代替D和J，并建立標量方程，然后獲得一個簡便的方程。

式（1）受到指定邊界處的電勢?或電流密度J邊界條件的影響。

3.3.2 采用非線性等效材料在MOA氣腔內(nèi)傳熱近似建立MOA的熱模型

MOA的熱傳遞有3種方式：傳導(dǎo)、對流和輻射。一般典型的瓷外套MOA是空腔結(jié)構(gòu)，芯體和外套之間的氣腔熱傳遞特性是非常關(guān)鍵的。熱輻射和自然對流傳熱決定著MOA在持續(xù)運行電壓下的電阻特性和暫態(tài)電壓下的熱穩(wěn)定性。對流傳熱的流體動態(tài)模型需要一個適合的建模方法和強大的計算資源。精確地計算MOA外表面的對流傳熱是比較困難的，因為影響外表對流的因素如傘形、風速等是不確定的。外部的輻射和對流熱傳遞可以在邊界條件中予以考慮。在MOA氣腔內(nèi)部，不規(guī)則的表面熱輻射和自然對流也需要考慮特定的求解算法。MOA氣腔內(nèi)部氣體的對流傳熱模型可以用一個具有一定導(dǎo)熱系數(shù)的介質(zhì)來等效：

式中：λeq為等效導(dǎo)熱系數(shù)；λair為空氣導(dǎo)熱系數(shù)；λconv為自然對流系數(shù)；λrad為熱輻射系數(shù)。包含熱容Cp和密度ρ的暫態(tài)熱傳遞方程：

式（3）描述了包含熱損耗率密度qj的物體溫度隨時間變化過程。對流傳熱由MOA的每節(jié)元件努賽爾數(shù)（Nueeslt number）Nu決定，表征一個封閉的垂直環(huán)形圓柱體。Nu數(shù)是無量綱的，由環(huán)形圓柱體的內(nèi)、外壁溫差ΔT、幾何尺寸（環(huán)形圓柱體的高度，內(nèi)、外徑）和溫度決定：

式中，Gr、Pr分別為 Grashof數(shù)和 Prandtl數(shù)。

等效熱輻射熱導(dǎo)率可以從斯特潘-玻爾茲曼定律（Stefan-Boltzmann law）中推導(dǎo)出來。假設(shè)熱導(dǎo)率在環(huán)形圓柱體中是恒定的，則

對于MOA來講，隨溫度變化的等效熱輻射熱導(dǎo)率可表示為

這樣計算MOA氣腔氣體的模型可以建立，見圖25[57]。

圖25 代表MOA氣腔的垂直環(huán)形圓柱體2D模型的幾何形狀Fig.25 The geometry of the 2D model of a vertical circular cylinder representing the MOA air chamber

在熱模型中，主要不確定性來自于所有內(nèi)部元件（例如：MOV芯體，瓷外套內(nèi)壁）和外部（例如瓷外套外壁，遮蔭情況的影響）表面的傳熱系數(shù)。對于無氣腔的MOA，情況得到了簡化，但在外套外表面上，對于不同的大氣邊界條件和遮蔭系數(shù)，實際的傳熱計算還是比較困難的。在大多數(shù)情況下，我們不得不只是暫時假設(shè)了外部接口處的傳熱系數(shù)值。

3.3.3 建立電、熱耦合數(shù)學計算模型

通過式（1）建立準靜電場（EQST）方程組來描述MOV工作情況，并用式（3）考慮其熱過程。通過MOV與溫度相關(guān)的電導(dǎo)率σ= σ（r，E，T）、焦耳熱耗散率q=q（r，E，T）E2=σ（r，E，T）[grad?]2，將電場與熱場進行耦合，見圖26[57]。這兩個方程可以采用弱耦合方法，如式（1）和式（3）是基于MOA幾何形狀的標準2維有限元離散化而連續(xù)求解的。前面我們知道分別表征MOA的電氣和熱過程具有差別很大的時間尺度。在MOA的正常工作區(qū)域內(nèi)，其電導(dǎo)率在一個AC電壓周波內(nèi)就有近6個數(shù)量級的變化，見圖26（d）。解決這些變化的必要時間步長通常為幾十μs，另一方面，MOA中的熱瞬態(tài)會在幾分鐘到幾小時內(nèi)才發(fā)生變化，見圖27（a）。Sp?ck-Leigsnering采用電場-熱場多步時間尺度耦合系統(tǒng)成功地解決此問題，見圖 27（b）[57]。

圖26 不同溫度下MOV的電場強度-電導(dǎo)率曲線Fig.26 Electric field-conductivity curves of MOV at different temperatures

圖27 MOA電場-熱場耦合系統(tǒng)關(guān)系圖Fig.27 MOA electric field-thermal field coupling system diagram

準靜電場的時間步長為Δtel，熱場的時間步長為Δtth。在AC電壓的一個短的時間周期內(nèi)，考慮到Δtel≤Δtth，電場可以達到一個穩(wěn)定態(tài)。假設(shè)溫度分布在這個熱時間步長內(nèi)（Δtth）持續(xù)保持一定，可以得到熱時間步長內(nèi)總的熱耗散q（T）Δtth。接著，可以在熱時間步長內(nèi)的計算式（3）。其計算結(jié)果得到的新的溫度分布結(jié)果，可以供下一個步驟進行計算電場分布。其中，考慮MOV的非線性電導(dǎo)率是非常關(guān)鍵的！

Sp?ck Leigsnering對圖26所示的樣機（但是，沒有均壓環(huán)），采用上述模型500 kV瓷外套MOA進行了仿真計算。計算顯示，在施加持續(xù)運行電壓Uc=345 kVr.m.s、5 h 后，溫度分布見圖28[57]。從圖中可看到，每節(jié)元件的溫度分布差異很大，上節(jié)元件與下節(jié)元件溫度相差近80℃，表明不同位置的元件工作在不同的U-I區(qū)域內(nèi)。圖28中虛線為實驗室MOA溫度分布實際測量值，與計算結(jié)果十分吻合。

圖28 在持續(xù)運行電壓下500 kV瓷外套MOA溫度分布圖Fig.28500 kV porcelain housing MOA temperature profile under continuous operating voltage

計算樣機的電位分布見圖29[57]，其中虛線為初始狀態(tài)下的電位分布曲線，實線為穩(wěn)定狀態(tài)下的電位分布曲線。我們可以觀察到一種“自均勻”現(xiàn)象。顯而易見，準靜電場熱模擬提供了更豐富的MOA信息。

圖29 在持續(xù)運行電壓下500 kV瓷外套MOA電位分布圖Fig.29500 kV porcelain housing MOA potential distribution diagram under continuous operating voltage

Sp?ck-Leigsnering對整只500 kV MOA，采用上述模型計算了MOA動作負載試驗的過程。選擇的時間步長見表1[57]。計算結(jié)果見圖30[57]，初始態(tài)和注入能量后的電位分布計算結(jié)果見圖31[57]?？梢钥吹剑跊_擊電流注入能量后，最上部元件電場強烈的不均勻性。這是因為沖擊后的大量能量導(dǎo)致的明顯溫度差異，并且因為電導(dǎo)率σ=σ（E，T）的變化，該部位的電場值也下降了。

表1 模擬動作負載試驗的時間步長參數(shù)Table 1 Time step parameter of simulation operation load test

圖30 模擬IEC規(guī)定的動作負載試驗中MOA溫度隨時間變化的分布圖Fig.30 Simulation of the distribution of MOA temperature over time in operating load test

圖31 500 kV瓷外套MOA電位分布圖Fig.31500 kV porcelain housing MOA potential distribution diagram

4 總結(jié)及展望

MOA是電力系統(tǒng)重要的過電壓保護裝置，目前已經(jīng)發(fā)展到了交流電力系統(tǒng)1000 kV、直流輸電系統(tǒng)±1100 kV電壓等級。在柔性交流輸電領(lǐng)域也得到了廣泛的應(yīng)用，例如串聯(lián)補償裝置TSC，靜止無功補償裝置（SVC）、潮流統(tǒng)一控制器（UPFS）等。但是，受到實驗條件的限制，我們對MOA的研究和試驗還局限于比例單元（特別是熱方面）。當今，隨著計算機技術(shù)的迅猛發(fā)展，給我們提供了關(guān)于整只MOA的電場和熱場全過程的研究條件，以便更經(jīng)濟、更高效地開展MOA研究和開發(fā)。

1）基于MOV的溫度函數(shù)的電導(dǎo)率特性的仿真數(shù)學計算是當前該領(lǐng)域研究重點和熱門。

2）整只MOA的數(shù)值模擬計算以及縮尺模型的實驗研究正在進行中，并可能成為標準化程序，從而實現(xiàn)安全、可靠的MOA的定制、設(shè)計。

3）在MOA實際結(jié)構(gòu)的熱-電場耦合的定量計算中，氣體計算流體動力學（CFD）熱部分與基于MOA內(nèi)部和外部固-氣界面處的熱傳遞的物理模型仍然存在問題。當計算機的能力得到進一步的提高時，我們必將解決此問題。

4）對于MOV元件的建模和仿真計算，將來主要的努力方向是改進其數(shù)學計算模型，以便它可以描述ZnO晶界與時間和溫度相關(guān)的特性。這樣我們就可以真實、定量地去描述，在給定的時間函數(shù)的電壓信號激勵下，MOV的電流和功率損耗響應(yīng)。我們在研究電場和溫度場的“自洽”（self-consistent）分布時，這些基本輸入是必須的！

5）必須看到中國在數(shù)值仿真計算MOA的研究上，與國際上還存在較大差距。