赵冬一，李 凡，汤 霖

（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]。如果要求对数函数用于极性变化的双极性模拟，那么进行计算时数值收敛可能难以接近电流零点。可以使用简化的PSpice元件，如二极管和电压源（图1（c））[14]。学者们还研究了使用较复杂的等效电路模型去更准确地模拟MOV的特性，例如，对于低电流区域，可能需要考虑电阻的温度系数以及电容的频率系数，见图1（d）[15]。

采用具有磁滞回线的非线性电阻模型可用于提高模拟计算的精度[16]。在U-I特性曲线的大电流翻转区域内，非线性电阻特性取决于冲击电流波前的陡度，IEEE工作组3.4.11给出了等效电路模型[17]（以下简称IEEE模型），见图2。基于IEEE模型的各种简化模型也相继出现，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）。此时，单个晶粒可以承载50%的总泄漏电流。在击穿区域内，电流变得更集中，见图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时间尺度的电行为和分钟、甚至小时时间尺度的热行为的相互耦合进行建模。总模拟时间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的研究上，与国际上还存在较大差距。