國際上基于數值仿真計算金屬氧化物避雷器的研究進展

趙冬一，李 凡，湯 霖

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

0 引言

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

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

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

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

基于電力系統設計、分析的需要，當系統有非線性電阻的MOA時，通常在暫態網絡分析系統中加入其等效電路模型。常見的等效電路模型見圖1[9-12]。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 二維和三維MOV導電機理的微觀模型計算

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

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

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

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

圖7 MOV的SEM照片Fig.7 SEM photograph of MOV

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

圖9 由隨機分布的種子點產生的Voronoi細胞結構Fig.9 Voronoi cell structure generated from randomly distributed seed points

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

圖10 Voronoi網絡方法計算的不同荷電率下電流導通路徑Fig.10 Voronoi network method to calculate the current conduction path under different charge rates

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

圖11 模擬電流定位現象（從白色到黑色的灰度級譜表示通過顆粒的電流的相對值，顆粒中的顏色越深，通過顆粒的電流越大）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內部的模擬溫度、熱應力分布情況Fig.12 MOV inside the simulation temperature，thermal stress distribution

圖13 MOV顯微結構內的電流分布Fig.13 The current distribution inside the microstructure of the ZnO varistor

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

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

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

圖15 圓柱形MOV模型的晶粒結構表示和分別為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特性的簡化，導致2D與3D的電流路徑有很大的預期差異，特別是尺寸較小的樣品和考慮了比例較高的晶界勢壘不均勻性的情況下。他們認為，3D材料模型可以更詳細和更真實地分析MOV的電流分布狀態。遺憾的是國內缺乏這方面的研究報道。

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

3.1 整只MOA在持續運行電壓下的數值仿真計算

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

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

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

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

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

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

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

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

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

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

2011年，ABB的Oliver Fritz等[49]基于COMSOL多物理場耦合模型，模擬研究了高壓GISMOA的熱行為。通過試驗數據擬合出MOV與溫度有關的非線性電導率，MOV的介電常數也以相同的方式確定，見圖19[48]。為了考慮從MOA柱到周圍氣體環境的熱傳遞，定義了MOA外表面的等效傳熱系數，并通過實驗數據來確定。他們指出了計算所需要處理的難點：對在ms時間尺度的電行為和分鐘、甚至小時時間尺度的熱行為的相互耦合進行建模?？偰M時間5 h，10 min通過重新計算電場并因此產生的焦耳熱來進行下一步的傳熱計算。因此，他們首先進行了瞬態熱量計算，再通過基于電場評估而重新計算的MOV的電阻率及功率損耗等。如此循環，最后分析評估整個MOA的溫度分布。Fritz等人表示：結果與試驗結果基本一致。與測量的結果相比，最大模擬溫度有些增加。這是由于保守的邊界條件導致的。用功率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]發表了他們在提高計算非線性電場效率方面研究的成果。引入“envelope”擬設準靜電場解算器，在網格內確定非線性材料特性的方式，自適應時間步進確保收斂，避免由時間步長之間的過度電荷轉移引起的大誤差以及通過場的適當空間和/或時間平均，在一階網格上，可以在合理的時間內獲得準確的解決方案。鄭重等[51]依此研究了空腔結構MOA的芯體調整金屬墊塊對散熱的影響。從芯體到外套壁為固定均勻的傳熱系數，而忽略熱引起的氣體對流局部傳熱的影響。按照IEC標準要求的動作負載試驗程序進行加載。改變傳熱系數和金屬墊塊的位置。鄭重等假設實際持續運行電壓后MOA很快達到熱平衡。他們重點關注了金屬墊塊對MOA的散熱效能。金屬墊塊在降低功率耗散和改善MOV的熱傳輸方面起著重要作用。2010年，已經到華北電力大學工作的鄭重教授與Steven A.Boggs、東芝公司的Toshiya Imai和Susumu Nishiwaki等[52]繼續進行了金屬墊塊對MOA散熱影響的研究，重點關注復合外套MOA的散熱機理。金屬墊塊配置對MOA散熱影響計算見圖21[52]。因此，為了表征樣本的傳熱參數，擬合了單獨熱冷卻過程的測量數據。他們得出結論，熱參數以及散熱器元件的位置有助于熱穩定性。這兩種效應可能無法在保守的縮小比例模型中正確表示。遺憾的是，鄭重等沒有對仿真計算結果進行試驗驗證。

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

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

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

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

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

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

3.2 整只MOA在限制暫態過電壓的數值仿真計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中，Gr、Pr分別為 Grashof數和 Prandtl數。

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

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

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

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

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

3.3.3 建立電、熱耦合數學計算模型

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

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

圖27 MOA電場-熱場耦合系統關系圖Fig.27 MOA electric field-thermal field coupling system diagram

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

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

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

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

圖29 在持續運行電壓下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]。計算結果見圖30[57]，初始態和注入能量后的電位分布計算結果見圖31[57]?？梢钥吹?，在沖擊電流注入能量后，最上部元件電場強烈的不均勻性。這是因為沖擊后的大量能量導致的明顯溫度差異，并且因為電導率σ=σ（E，T）的變化，該部位的電場值也下降了。

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

圖30 模擬IEC規定的動作負載試驗中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 總結及展望

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

1）基于MOV的溫度函數的電導率特性的仿真數學計算是當前該領域研究重點和熱門。

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

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

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

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