FDTD Electromagnetic and Thermal Simulation of a Metal Oxide Varistor Element Considering the Temperature Dependence of Its Resistivity

: Electromagnetic ﬁelds and the heat of a metal oxide varistor (MOV), in which a lightning impulse current ﬂows, are calculated using the ﬁnite-difference time-domain (FDTD) method. The MOV is represented with small rectangular parallelepiped cells, each of which has a resistivity dependent on electric ﬁeld and temperature. For this purpose, the expression of resistivity as a function of the electric ﬁeld, proposed previously, is extended to include the dependence on temperature. The temperature dependence is based on voltages across an MOV for impulse currents of 0.5 to 10 kA at temperatures in a range from about 300 to 900 K, measured by Andoh et al. (2000). FDTD-calculated waveform of voltage across the MOV agrees well with the corresponding measured one for a short impulse current with a magnitude of about 4 kA and a duration of about 30 µ s. In addition, the temperature on the surface of the MOV agrees well with the corresponding measured one. Further, calculations are carried out for the MOV with a nonuniform resistivity distribution, which roughly simulates deterioration or degradation of the MOV, for a long duration current having a magnitude of 5 kA. The proposed expression of resistivity, given as a function of electric ﬁeld and temperature, is useful in studying electro-thermal calculations, which can provide insights into causes of MOV damages.


Introduction
Metal oxide varistors (MOVs) or zinc oxide (ZnO) varistors are widely installed in electrical systems to protect their equipment against overvoltages, such as lightning surge voltages (e.g., [1][2][3]). MOVs behave as insulators for normal-operation voltages, but as conductors for overvoltages owing to their highly nonlinear voltage-versus-current (V-I) properties. Therefore, overvoltages are not applied to the equipment connected in parallel with an MOV.
In [4], the distribution of voltage along the surface of a metal oxide (MO) arrester has been studied using the surface charge simulation method, since the nonuniform voltage distribution or locally intense electric field might damage the arrester or shorten the lifetime. On the basis of the calculated results, a countermeasure to improve the voltage distribution has been proposed. Similarly, in [5], the distribution of electric field around an MO arrester has been analyzed using the finite element method (FEM) and applying artificial neural networks. In [6], the distribution of voltage along an MO arrester has been studied using the FEM, and also related recent works have been finely reviewed.
Intensive heat generated in an MOV, in which a high current flows, might also damage the MOV, shorten the lifetime or cause partial melting. Therefore, electro-thermal calculations of MOVs or MO arresters have been carried out (e.g., [7][8][9]). Andoh et al. [10] have shown that the resistivity of MOV material starts decreasing significantly with increasing temperature at some critical high temperature. If the distribution of high current over the cross-section of an MOV is not uniform, the local concentration of high current might be more significant because of the resistivity decrease with increasing temperature and might cause partial melting [10].
To study the transient nonuniform temperature rise in an MOV, in which an impulse current flows, the finite-difference time-domain (FDTD) method [11] is of use. This is because the resistivity of each of the small cells, which compose the MOV, can be controlled as a function of the electric field and temperature there.
In this paper, the voltage generated across an MOV, in which a lightning impulse current flows, and the temperature distribution are calculated using the FDTD method. The MOV is represented with small rectangular parallelepiped cells, each of which has a resistivity dependent on electric field and temperature. The temperature dependence is based on voltages across an MOV for impulse currents of 0.5 to 10 kA at temperatures in a range from about 300 to 900 K, measured by Andoh et al. [10]. The heat is calculated with the heat transport equation from FDTD-calculated electric and magnetic fields. The aim of this work is to extend the expression of resistivity as a function of electric field, proposed previously [12], to include the dependence on temperature, and to show its usefulness in electro-thermal calculations of an MOV with an initially uniform or nonuniform (from the beginning) resistivity distribution, in which an impulse high current flows. The results will provide insights into causes of MOV damages. Figure 1 shows the V-I property of an MOV element, measured at Otowa Electric Co. Ltd., Osaka, Japan, and employed in this paper. The voltage per millimeter (V 1mA ), appearing across this MOV element when a direct current with a magnitude of 1 mA flows in, is 500 V.

Thermal Calculation
The temperature T (K) of an MOV element is obtained by solving the heat equation [13,14] given below: where t is time (s), Pd is the absorbed power per unit volume (W/m 3 ), ρm is the mass density (kg/m 3 ), Cm is the specific heat (J/(kg K)), κ is the thermal conductivity (W/(Mk),α = κ/(ρmCm) is the thermal diffusivity (m 2 /s). The discretized equation of (3) is given as follows:  A simple mathematical expression, which has only three adjustable constants to represent nonlinear properties of MOV material, has been proposed for FDTD-based electromagnetic and surge calculations in [12]. It has been shown that the resistivity versus electric-field (ð-E) properties, which are based on the measured voltage-versus-current properties, are approximated reasonably well by the following expression: where c 0 , c 1 , and c 2 are adjustable constants, and are determined by applying the leastsquare method to the entire measured property of log 10 ρ versus log 10 E.
The following expression of resistivity as a function of electric field, based on (1), is used in the time update equations for the electric fields in the FDTD calculation: This expression is useful in FDTD calculations because it can represent measured ρ-E properties approximately and does not require iterative calculations to find for any value of E.

Thermal Calculation
The temperature T (K) of an MOV element is obtained by solving the heat equation [13,14] given below: where t is time (s), P d is the absorbed power per unit volume (W/m 3 ), ρ m is the mass density (kg/m 3 ), C m is the specific heat (J/(kg K)), κ is the thermal conductivity (W/(Mk), α = κ/(ρ m C m ) is the thermal diffusivity (m 2 /s). The discretized equation of (3) is given as follows: where ∆t' is the thermal-calculation time step. The absorbed power per unit volume at the time step n + 1 is expressed as follows: The power P d absorbed at the location (i, j, k) at the time step n + 1 is calculated by where ∆t is the time step of electromagnetic-field calculation, P x m (i,j,k), P y m (i,j,k), and P z m (i,j,k) are instantaneous spatial average values of neighboring four components at the time step m of the absorbed power in x, y, and z directions, n s is the first time step number in the thermal-calculation time step ∆t', n e is the last time step number in ∆t'.
To the surface of the MOV element, the convective boundary condition [14] is applied. First, (4) is solved using the adiabatic boundary condition. Then, the amount of heat flowing out to the surrounding air during ∆t' is calculated as follows: where ∆Q n+1 (i, j, k) is the heat transfer rate (W) at the time step of n + 1, ∆T n+1 (i, j, k) is the surface temperature of the material (K) at the time step of n + 1, T air (i, j, k) is the ambient temperature, which is set to 293 K, h is the convection heat transfer rate, which is set to 10 W/(m 2 , and ∆S is the surface area of the boundary of the FDTD cell. The temperature drop ∆T ((K)) on the material surface at the time step n + 1 is calculated by where ∆V (m 3 ) is the volume of the FDTD cell. The temperature of the surface boundary ∆T' n+1 (i, j, k) (K), which is updated considering the temperature drop due to the heat transfer to the surrounding air, is provided below:

Temperature Dependence of MOV Material
Andoh et al. [10] have found that, for 8/20 µs impulse currents of 0.5 to 10 kA, the resistivity of an MOV element increases slightly with increasing temperature up to about 550 or 600 K, but it starts decreasing significantly when the temperature exceeds about 700 K (Figure 7 of [10]). This temperature dependence can be reproduced by making two constants c 0 and c 1 of Equation (1) be parameters c 0 (T) and c 1 (T) dependent on temperature T (K) as follows: where T m is the temperature at which the MOV starts melting (K), and T 0 is the ambient temperature (K). In this work, T m and T 0 are set to 900 K and 293 K, respectively. Figure 2 shows ρ-E properties measured at an initial temperature of 293 K, and the corresponding curve drawn using Equations (10)- (12). Values of the three adjustable constants are c 0 = −1.78 (negative value), c 1 = 8.44 × 10 51 , and c 2 = −68.9 (negative value), respectively, which are determined by applying the least-square method to the entire measured ρ-E property. ρ-E properties are also shown at T = 493, 693 and 893 K drawn using the same expressions with the same constants as the above.
where ΔV (m 3 ) is the volume of the FDTD cell. The temperature of the surface boundary ΔT' n+1 (i, j, k) (K), which is updated considering the temperature drop due to the heat transfer to the surrounding air, is provided below:

Temperature Dependence of MOV Material
Andoh et al. [10] have found that, for 8/20 μs impulse currents of 0.5 to 10 kA, the resistivity of an MOV element increases slightly with increasing temperature up to about 550 or 600 K, but it starts decreasing significantly when the temperature exceeds about 700 K (Figure 7 of [10]). This temperature dependence can be reproduced by making two constants c0′ and c1′ of Equation (1) be parameters c0 (T) and c1 (T) dependent on temperature T (K) as follows: where Tm is the temperature at which the MOV starts melting (K), and T0 is the ambient temperature (K). In this work, Tm and T0 are set to 900 K and 293 K, respectively. Figure 2 shows ρ-E properties measured at an initial temperature of 293 K, and the corresponding curve drawn using Equations (10)- (12). Values of the three adjustable constants are c0′ = −1.78 (negative value), c1′ = 8.44 × 10 51 , and c2′ = −68.9 (negative value), respectively, which are determined by applying the least-square method to the entire measured ρ-E property. ρ-E properties are also shown at T = 493, 693 and 893 K drawn using the same expressions with the same constants as the above.  Figure 3 shows temperature dependences of voltage per millimeter across the MOV element for impulse currents of 2, 3, and 5 kA. In the calculations, the resistivity of each cell at each time step is determined with Equations (10)- (12). It appears from Figure  3 that the resistivity increases slightly with increasing temperature up to about 550 or 600 K, and then at about 700 K it starts decreasing significantly. This reproduces well the measured temperature dependence described in the beginning of this section, which shows the validity of the proposed model. It is assumed that the thermal conductivity is constant. Since the dependence of thermal conductivity of MOV is smaller than that of electrical conductivity or resistivity [15], influences of the variation of thermal conductivity with temperature on the calculated voltages and temperatures are expected to be less significant. the resistivity increases slightly with increasing temperature up to about 550 or 600 K, and then at about 700 K it starts decreasing significantly. This reproduces well the measured temperature dependence described in the beginning of this section, which shows the validity of the proposed model. It is assumed that the thermal conductivity is constant. Since the dependence of thermal conductivity of MOV is smaller than that of electrical conductivity or resistivity [15], influences of the variation of thermal conductivity with temperature on the calculated voltages and temperatures are expected to be less significant.  Figure 4a shows the side view of a test circuit comprising an MOV element, a lumped current source, and perfectly conducting lead wires. Figure 4b shows plan and side views of the MOV element. The element thickness is d = 1.5 mm, and the cross-sectional area is S = 144 mm 2 . This type of MOV is usually used for systems of about 300 to 400 V, and its maximum discharge capacity of lightning impulse current is about 4 to 5 kA. The MOV element is represented with 1 mm × 1 mm × 0.3 mm cells. Each of the MOV cells has a resistivity in each of the x, y, and z directions, depending on the electric field and temperature there. The relative permittivity of the MOV is set to 800 because MOV material has a relative permittivity value of several hundreds [7]. In addition, the thermal conductivity is set to 29.2 W/(m K), the specific heat is set to 480 J/(kg K), and the mass density is set to 5.394 × 10 3 kg/m 3 [16]. The thickness of each of upper and lower silver resin electrodes is 0.3 mm, and the cross-sectional area is S = 100 mm 2 . Its conductivity, relative permittivity, specific heat, thermal conductivity, and mass density are set to 6.1 × 10 7 S/m, 1, 500 J/(kg K), 25 W/(m K), and 6.1 × 10 3 kg/m 3 , respectively. Note that the edge of the silver resin electrode is covered with resin. The resin resistivity, relative permittivity, specific heat, thermal conductivity, and mass density are set to 1.8 × 10 8 Ωm, 9.2920 J/(kg K), 0.81 W/(m K), and 1.8 × 10 3 kg/m 3 , respectively. This model is accommodated in a working volume of 40 mm × 60 mm × 18 mm. Liao's 2nd-order absorbing boundary condition [17] is applied to each surface of the working volume.

Modeling
The time step of electromagnetic calculation is set to t = 0.456 ps on the basis of the Courant-Friedrichs-Lewy condition to keep the numerical stability, and the time step of thermal calculation is set to ∆t' = 0.182 μs (=400,000 t). The maximum observation time is set to 32 μs. 3. Modeling Figure 4a shows the side view of a test circuit comprising an MOV element, a lumped current source, and perfectly conducting lead wires. Figure 4b shows plan and side views of the MOV element. The element thickness is d = 1.5 mm, and the cross-sectional area is S = 144 mm 2 . This type of MOV is usually used for systems of about 300 to 400 V, and its maximum discharge capacity of lightning impulse current is about 4 to 5 kA. The MOV element is represented with 1 mm × 1 mm × 0.3 mm cells. Each of the MOV cells has a resistivity in each of the x, y, and z directions, depending on the electric field and temperature there. The relative permittivity of the MOV is set to 800 because MOV material has a relative permittivity value of several hundreds [7]. In addition, the thermal conductivity is set to 29.2 W/(m K), the specific heat is set to 480 J/(kg K), and the mass density is set to 5.394 × 10 3 kg/m 3 [16]. The thickness of each of upper and lower silver resin electrodes is 0.3 mm, and the cross-sectional area is S = 100 mm 2 . Its conductivity, relative permittivity, specific heat, thermal conductivity, and mass density are set to 6.1 × 10 7 S/m, 1, 500 J/(kg K), 25 W/(m K), and 6.1 × 10 3 kg/m 3 , respectively. Note that the edge of the silver resin electrode is covered with resin. The resin resistivity, relative permittivity, specific heat, thermal conductivity, and mass density are set to 1.8 × 10 8 Ωm, 9.2920 J/(kg K), 0.81 W/(m K), and 1.8 × 10 3 kg/m 3 , respectively. This model is accommodated in a working volume of 40 mm × 60 mm × 18 mm. Liao's 2nd-order absorbing boundary condition [17] is applied to each surface of the working volume.
The time step of electromagnetic calculation is set to t = 0.456 ps on the basis of the Courant-Friedrichs-Lewy condition to keep the numerical stability, and the time step of thermal calculation is set to ∆t' = 0.182 µs (=400,000 t). The maximum observation time is set to 32 µs.   Figure 5 shows the waveform of current injected in the MOV element in the experiment. The approximated waveform, which is used as the output waveform of the current source in the FDTD calculation, is also shown in Figure 5. Note that initial high-frequency oscillations superposed on the measured waveform are electromagnetic noises probably   Figure 5 shows the waveform of current injected in the MOV element in the experiment. The approximated waveform, which is used as the output waveform of the current source in the FDTD calculation, is also shown in Figure 5. Note that initial high-frequency oscillations superposed on the measured waveform are electromagnetic noises probably generated from the impulse current generator, which are neglected in the calculation.   Figure 5 shows the waveform of current injected in the MOV element in the experiment. The approximated waveform, which is used as the output waveform of the current source in the FDTD calculation, is also shown in Figure 5. Note that initial high-frequency oscillations superposed on the measured waveform are electromagnetic noises probably generated from the impulse current generator, which are neglected in the calculation.  Figure 6 shows measured and calculated waveforms of voltage generated across the MOV. Both the waveforms calculated with and without the temperature dependence agree well with the measured one. The temperature on the side surface of the MOV was about 293 K just before the measurement, and therefore in the calculation it was set to 293 K.

Time [μs]
Computed with temperature dependence Computed without temperature dependence Measured  Figure 6 shows measured and calculated waveforms of voltage generated across the MOV. Both the waveforms calculated with and without the temperature dependence agree well with the measured one. The temperature on the side surface of the MOV was about 293 K just before the measurement, and therefore in the calculation it was set to 293 K.   Figure 5 shows the waveform of current injected in the MOV element in the experiment. The approximated waveform, which is used as the output waveform of the current source in the FDTD calculation, is also shown in Figure 5. Note that initial high-frequency oscillations superposed on the measured waveform are electromagnetic noises probably generated from the impulse current generator, which are neglected in the calculation.  Figure 6 shows measured and calculated waveforms of voltage generated across the MOV. Both the waveforms calculated with and without the temperature dependence agree well with the measured one. The temperature on the side surface of the MOV was about 293 K just before the measurement, and therefore in the calculation it was set to 293 K.

Time [μs]
Computed with temperature dependence Computed without temperature dependence Measured    Figure 7 shows the time variations of surface and internal temperatures of the MOV calculated with (solid lines) and without (dashed lines) temperature dependence. It follows from Figure 7 that the temperature attains its maximum value at a time of 32 μs, when the injected current decreases to zero.    Table 1 gives the measured and calculated maximum temperature values on the surface of the MOV element. The measurement was carried out using a thermographic camera with a framerate of 30 fps (Avionics TVS-700). Therefore, the measured temperature is probably the maximum within several tens of milliseconds after the impulse current was injected. Since the injected current has a short duration (about 30 μs) and no additional heat is generated after about 30 μs, the comparison of measured and calculated temperatures in such different time scales is not much unreasonable. The maximum temperatures, calculated with and without the temperature dependence, on the MOV side surface are both 356 K. These are close to the corresponding measured temperature, 359 K. Since the significant decrease in resistivity does not occur up to about 550 or 600 K, these two calculated temperatures are almost the same. In addition, the electrical energy converted into thermal energy, calculated integrating the FDTD-calculated MOV voltage multiplied by current over time, is about 60 J, which is almost the same as the corresponding one (about 59 J) calculated from the measured voltage and current waveforms.    Figure 7 shows the time variations of surface and internal temperatures of the MOV calculated with (solid lines) and without (dashed lines) temperature dependence. It follows from Figure 7 that the temperature attains its maximum value at a time of 32 μs, when the injected current decreases to zero.    Table 1 gives the measured and calculated maximum temperature values on the surface of the MOV element. The measurement was carried out using a thermographic camera with a framerate of 30 fps (Avionics TVS-700). Therefore, the measured temperature is probably the maximum within several tens of milliseconds after the impulse current was injected. Since the injected current has a short duration (about 30 μs) and no additional heat is generated after about 30 μs, the comparison of measured and calculated temperatures in such different time scales is not much unreasonable. The maximum temperatures, calculated with and without the temperature dependence, on the MOV side surface are both 356 K. These are close to the corresponding measured temperature, 359 K. Since the significant decrease in resistivity does not occur up to about 550 or 600 K, these two calculated temperatures are almost the same. In addition, the electrical energy converted into thermal energy, calculated integrating the FDTD-calculated MOV voltage multiplied by current over time, is about 60 J, which is almost the same as the corresponding one (about 59 J) calculated from the measured voltage and current waveforms.  Figure 5. The initial temperature is set to 293 K. Table 1 gives the measured and calculated maximum temperature values on the surface of the MOV element. The measurement was carried out using a thermographic camera with a framerate of 30 fps (Avionics TVS-700). Therefore, the measured temperature is probably the maximum within several tens of milliseconds after the impulse current was injected. Since the injected current has a short duration (about 30 µs) and no additional heat is generated after about 30 µs, the comparison of measured and calculated temperatures in such different time scales is not much unreasonable. The maximum temperatures, calculated with and without the temperature dependence, on the MOV side surface are both 356 K. These are close to the corresponding measured temperature, 359 K. Since the significant decrease in resistivity does not occur up to about 550 or 600 K, these two calculated temperatures are almost the same. In addition, the electrical energy converted into thermal energy, calculated integrating the FDTD-calculated MOV voltage multiplied by current over time, is about 60 J, which is almost the same as the corresponding one (about 59 J) calculated from the measured voltage and current waveforms. It appears from above that the temperature-dependent resistive property of the MOV element within the considered range can be expressed or reproduced appropriately using Equations (10)-(12).

Discussion
Calculations are carried out for the MOV having a nonuniform resistivity distribution from the beginning, which roughly simulates deterioration or degradation of the MOV. The initial resistivity in each of the x, y, and z directions of each cell is set randomly to a value in a range from 0.99 ρ to 1.01 ρ (±1% variation), from 0.95 ρ to 1.05 ρ, or from 0.9 ρ to 1.1 ρ using a pseudorandom number generator of FORTRAN. Note that the use of multivariate random resistivity distribution considering some correlation among three directions can be one of the options, although a uniform random resistivity distribution is employed in this work. The initial temperature is set to 293 K. The injected ramp current has a risetime of 1 µs and a magnitude of 5 kA (constant after 1 µs). Figure 9 shows calculated waveforms of voltage across the MOV. In Figure 9, the waveform of voltage is also shown across the MOV, having uniform resistivity initially (±0% variation). The voltages drop abruptly to zero in a time range of about 45 to 50 µs, when the temperature exceeds about 1000 K, and the resistivity is significantly low owing to the temperature dependence specified by Equations (10)- (12). If no temperature dependence of resistivity is considered, the voltage does not drop to zero, but keeps about 1.2 kV.
It appears from above that the temperature-dependent resistive property of the MOV element within the considered range can be expressed or reproduced appropriately using Equations (10)-(12).

Discussion
Calculations are carried out for the MOV having a nonuniform resistivity distribution from the beginning, which roughly simulates deterioration or degradation of the MOV. The initial resistivity in each of the x, y, and z directions of each cell is set randomly to a value in a range from 0.99 ρ to 1.01 ρ (±1% variation), from 0.95 ρ to 1.05 ρ , or from 0.9 ρ to 1.1 ρ using a pseudorandom number generator of FORTRAN. Note that the use of multivariate random resistivity distribution considering some correlation among three directions can be one of the options, although a uniform random resistivity distribution is employed in this work. The initial temperature is set to 293 K. The injected ramp current has a risetime of 1 μs and a magnitude of 5 kA (constant after 1 μs). Figure 9 shows calculated waveforms of voltage across the MOV. In Figure 9, the waveform of voltage is also shown across the MOV, having uniform resistivity initially (±0% variation). The voltages drop abruptly to zero in a time range of about 45 to 50 μs, when the temperature exceeds about 1000 K, and the resistivity is significantly low owing to the temperature dependence specified by Equations (10)- (12). If no temperature dependence of resistivity is considered, the voltage does not drop to zero, but keeps about 1.2 kV. Figure 9. Waveforms of voltage across the MOV element with uniform or nonuniform resistivity distribution from the beginning, calculated using the FDTD method. The variation of resistivity is in a range of ±0% (uniform), ±1%, ±5%, or ±10%. Table 2 gives FDTD-calculated time interval and generated thermal energy until the voltage across the MOV dropped to zero. The resistivity variation of the MOV element little influences the generated thermal energy.  ± 0% ± 1% ± 5% ± 10% Figure 9. Waveforms of voltage across the MOV element with uniform or nonuniform resistivity distribution from the beginning, calculated using the FDTD method. The variation of resistivity is in a range of ±0% (uniform), ±1%, ±5%, or ±10%. Table 2 gives FDTD-calculated time interval and generated thermal energy until the voltage across the MOV dropped to zero. The resistivity variation of the MOV element little influences the generated thermal energy.  Figure 10 shows temperature distributions on the upper surface of the MOV at times of 15, 30, 43.2, and 45 µs for the resistivity variation of ±10%. The temperature distribution is almost uniform until it goes up to about 600 K. Beyond it, the current density becomes more significant where the resistivity is lower because of the locally higher temperature rise. In addition, it shows that the temperature dependence of MOV resistivity, which decreases significantly above 700 K [10], is well simulated.

±1%
49.4 266 ±5% 47.0 257 ±10% 45.2 249 Figure 10 shows temperature distributions on the upper surface of the MOV at times of 15, 30, 43.2, and 45 μs for the resistivity variation of ±10%. The temperature distribution is almost uniform until it goes up to about 600 K. Beyond it, the current density becomes more significant where the resistivity is lower because of the locally higher temperature rise. In addition, it shows that the temperature dependence of MOV resistivity, which decreases significantly above 700 K [10], is well simulated.

Conclusions
The voltage generated across an MOV, in which an impulse current was injected, and the temperature distribution were calculated with the FDTD method. The MOV was represented with small cells, each of which had a resistivity dependent on electric field and temperature. For this purpose, an expression of resistivity as a function of electric field and temperature was proposed. Calculated waveform of voltage across the MOV agreed well with the corresponding measured one for an impulse current with a magnitude of about 4 kA and a duration of about 30 μs. In addition, the temperature on the surface of the MOV agreed well with the corresponding measured one. Further, calculations were carried out for the MOV with a nonuniform resistivity distribution, which roughly simulated deterioration or degradation of the MOV, for a long duration current having a magnitude of 5 kA. The initial resistivity of each cell was set randomly to a value in a range of ±1%, ±5% or ±10% variation. The calculated voltages dropped abruptly to zero in a time range of about 45 to 50 μs, when the temperature exceeded about 1000 K and the resistivity was significantly low owing to the temperature dependence specified by Equations (10)- (12). The temperature distribution had been almost uniform until it went up to about 600 K. Beyond it, the current density became more significant, where the resistivity was lower because of the locally higher temperature rise.
The proposed expression of temperature-dependent resistivity can be applied to impulse currents having any magnitude and waveshape, and it can work properly even beyond 1000 K. The MOV model, represented with small cells, can be extended to other MOVs or arresters having different geometries. It appears that the proposed expression of resistivity and the FDTD model of MOV are useful in studying electro-thermal calculations, which can provide insights into causes of MOV damages.

Conclusions
The voltage generated across an MOV, in which an impulse current was injected, and the temperature distribution were calculated with the FDTD method. The MOV was represented with small cells, each of which had a resistivity dependent on electric field and temperature. For this purpose, an expression of resistivity as a function of electric field and temperature was proposed. Calculated waveform of voltage across the MOV agreed well with the corresponding measured one for an impulse current with a magnitude of about 4 kA and a duration of about 30 µs. In addition, the temperature on the surface of the MOV agreed well with the corresponding measured one. Further, calculations were carried out for the MOV with a nonuniform resistivity distribution, which roughly simulated deterioration or degradation of the MOV, for a long duration current having a magnitude of 5 kA. The initial resistivity of each cell was set randomly to a value in a range of ±1%, ±5% or ±10% variation. The calculated voltages dropped abruptly to zero in a time range of about 45 to 50 µs, when the temperature exceeded about 1000 K and the resistivity was significantly low owing to the temperature dependence specified by Equations (10)- (12). The temperature distribution had been almost uniform until it went up to about 600 K. Beyond it, the current density became more significant, where the resistivity was lower because of the locally higher temperature rise.
The proposed expression of temperature-dependent resistivity can be applied to impulse currents having any magnitude and waveshape, and it can work properly even beyond 1000 K. The MOV model, represented with small cells, can be extended to other MOVs or arresters having different geometries. It appears that the proposed expression of resistivity and the FDTD model of MOV are useful in studying electro-thermal calculations, which can provide insights into causes of MOV damages.