Bulk FDTD Simulation of Distributed Corona Effects and Overvoltage Profiles for HSIL Transmission Line Design

Abstract

Power system load growth and transmission corridor constraints are driving industry activity in the area of high surge impedance loading (HSIL). Examples include compact structure design and uprating existing transmission lines. Recent research relating electric field uniformity to transmission line capacity and critical flashover voltage underscored the need for better overvoltage data to quantify insulation margins for HSIL design. To that end, this work extends the finite difference time domain (FDTD) method with distributed corona losses to transmission lines with bundled conductors. The model was adapted for practical use in high-volume statistical transient simulation and applied to an example 500 kV line. Transients included line energization and trapped charge reclosing. Overvoltage profiles and statistical distributions were generated from 9500 simulations obtained by random breaker close timing and variation in line length and altitude. Distributed corona losses reduced 98th percentile line-to-ground switching overvoltages by 4%–14% of nominal. The estimated line-to-ground switching surge flashover probability was 54%–80% lower with corona loss. Corona had less impact on line-to-line overvoltages, but the effects were still notable. Results highlight the importance of considering detailed overvoltage profiles and accounting for corona loss attenuation when seeking to carefully quantify insulation design margins.

1. Introduction

Transmission owners frequently face the challenge of accommodating power system load growth despite transmission corridor constraints and lack of access to new corridors. Consequently, it is important to use existing corridors as efficiently as possible [1,2,3]. This has led to ongoing research and the development of methods for uprating existing transmission lines and designing compact transmission line structures [4,5]. Recent research has explicitly connected increased power capacity to reduced tolerance of voltage surges. This has highlighted the importance of carefully quantifying all factors affecting overvoltage distributions [6]. Such information will help engineers better understand design margins and maximize capacity while maintaining reliability. The objective of the present work was to implement a practical simulation approach capable of generating data not normally available for use in the design of overhead transmission line insulation. These data include detailed switching overvoltage profiles with high spatial resolution and quantified impacts of distributed corona losses. The approach was demonstrated through bulk statistical simulation of an example 500 kV transmission line. Results were used to draw preliminary conclusions about the possible importance of having such data for use in refining insulation design as part of capacity optimization. The remainder of the introduction provides additional detail regarding the background and context of the research. Future research will assess preliminary conclusions by applying the model to a broader set of design cases.

It has been shown that for relatively long transmission lines, the power capacity of the line can be reasonably approximated as its “surge impedance loading” [7,8]. More specifically, the approximate capacity can be written as shown in (1).

$$P_{SIL}={\left|V_{LL}\right|}^2/Z_0$$

(1)

where $V_{LL}$ is the line-to-line voltage (in kilovolts) and $Z_0$ is the positive sequence “surge impedance” of the transmission line. Historically, the need to increase power capacity of transmission lines has driven the increase in transmission voltage. If it is not possible to increase the voltage, a smaller capacity improvement can be achieved by reducing the transmission line’s surge impedance.

This idea has been investigated for some time [9,10,11,12,13,14,15] and is known in the literature as high surge impedance loading (HSIL). Experience has shown that capacity increases of approximately 30% can be obtained with designs not significantly different from traditional designs. One recent example of this is the AEP BOLD design [16]. A good review of the work in this area can be found in [17].

It was recently shown that the surge impedance of a three-phase transmission line can be directly related to the spatial distribution of the electric field between its phases [6]. The more spatially uniform the electric field between phase conductors, the smaller the surge impedance and, hence, the larger the surge impedance loading. However, it was also shown in [6] that electric field uniformity is inversely related to the line-to-line critical flashover voltage (CFO). Since CFO reduction is limited by the statistical switching overvoltage distribution of the transmission line, it follows that the switching surge distribution likewise imposes a limit on the capacity improvement available from HSIL designs. It is this fact that drives the need for the research presented in this paper.

This observation is not entirely new; however, it has not been explicitly related to the spatial distribution of the electric field between phase conductors nor to reduced surge impedance. It was, for example, recognized in [18] that for HSIL designs, more care has to be taken with respect to protection against surges. The authors of [18] specifically discussed the potential need for surge arresters, circuit breakers with pre-insertion resistors, or controlled closing of the breaker poles.

The specific finding in [6] was that the ratio of line-to-line CFO to operating voltage is smaller for designs with more spatially uniform electric fields. In other words, the transmission line becomes less tolerant of surges the more spatially uniform the electric field (i.e., the higher the surge impedance loading). It follows that a detailed understanding of the surge voltage profile is important, particularly when seeking to carefully quantify design margins to optimize capacity.

It is interesting to note that these issues are related to transmission line corona. A transmission line is designed so that, at its operating voltage, the effects of conductor corona (e.g., audible or radio noise) are just tolerable [6]. Corona effects are generally considered limiting constraints for compact transmission structures, yet, there is a possible benefit. Corona losses, which become significant for voltages higher than the operating voltage, generally act to attenuate surges that propagate along a transmission line [19]. Thus, corona losses can result in smaller probability of flashover.

At present, traditional electromagnetic transients programs (EMTP) can account for transmission line conductor losses and earth losses. They also allow modeling of applicable lumped nonlinear devices such as surge arresters. However, they do not easily account for corona loss, which is a distributed nonlinear phenomenon [20]. Various corona models have been proposed for EMTP-type simulation [21,22,23], but the industry standard software has yet to incorporate a model available for the typical user. Consequently, industry practice does not presently include corona losses in bulk statistical switching overvoltage studies. A large portion of the flashover risk from an overvoltage distribution is attributable to a relatively small percentage of the highest magnitude surges [24]. Thus, any attenuation of the highest surges by corona losses could be a welcome side effect that is not typically accounted for in insulation design studies.

At the outset, attenuation of the voltage by corona may seem too small compared to that introduced by line arresters which are often installed at the line terminals. However, the effective electrical reach of arrester protection is limited. Transient overvoltage maxima often occur away from the ends of arrester terminated lines [25]. In addition, arresters do not have much impact until voltages exceed 2.0 per unit. A conventional study approach is to divide EMTP line models into a handful of segments with voltage probes at the junctions between segments [26]. Overvoltage profiles are interpolated from these few measurement points. Historically, this has been satisfactory, but profiles with higher spatial resolution could provide detailed information to help optimize line design for highly constrained situations. Few studies are found in the literature that discuss high-resolution overvoltage profiles in the context of statistical switching studies. There is also little information regarding the nature of phase-to-phase voltage surge distributions which could be the controlling case for HSIL lines, especially for tower designs with no grounded conductors between phases [16,17].

A finite difference time domain (FDTD) model was developed to estimate the impact of distributed corona losses on transmission line flashover probability. The model was also used to determine if detailed information from switching overvoltage profiles with high spatial resolution could benefit transmission optimization studies. The research highlights the differences between line-to-ground and line-to-line overvoltage profile characteristics. The model was demonstrated through analysis of a realistic 500 kV transmission line. Bulk simulations were performed to generate batches of switching surge data for statistical calculations. Corona losses reduced transient overvoltages by between 4% and 14% of nominal for 98th percentile line-to-ground exposures (Section 4.5). Results vary with overvoltage severity and corona onset conditions as affected by conductor and bundle geometry, altitude, atmospheric conditions, etc. The corresponding line-to-ground flashover probability was reduced by 54%–80% (Section 4.6). In general, line-to-line overvoltage profiles are less affected by corona and are flatter than those of line-to-ground exposures.

The research demonstrates that more detailed information about overvoltage profiles and distributed corona losses can benefit rigorous HSIL optimization. The information could also help determine, with greater certainty, whether costs must be incurred to mitigate transients through such means as pre-insertion resistors or controlled breaker closing schemes. Other factors, such as lightning and contamination performance, are also important for transmission line insulation design [24], but are outside the scope of this study. The next two sections summarize the example design scenarios and the FDTD model. The final two sections present results and conclusions.

2. Design Scenarios

A 500 kV transmission line of varying length was selected as an initial demonstration system for the simulation model.

Table 1. The 500 kV transmission line characteristics.

Parameter	Values	Notes
Nominal/maximum voltage	500 kV/550 kV	Typical max as per [28]
Structure type	H-Frame
Phase spacing	10 m
Average conductor height	16 m
Conductor	ACSR 1272 Bittern	21 outer strands
Bundling	Triple	45.72 cm spacing
Shield wires	None
Line length	50, 100, 250, 500, and 800 km
Soil resistivity	100 ohm-m	Typical as per [29]
Earth relative permittivity	10	Typical as per [30]
Altitudes considered	0 and 2000 m

summarizes the characteristics of the line. Phase spacing is relatively narrow, selected to be somewhat representative of a structure that borders on being compact. Conductor size and bundling are selected such that audible noise and radio interference are just within the recommended limits [27] for altitudes up to 2000 m. Studies were performed with line lengths varying from 50 to 800 km.

The switching scenarios completed for the research are summarized in

Table 2. Study scenarios.

	Transient Type	Altitude (meters)	Corona Enabled?	50 km	100 km	250 km	500 km	800 km
1	Energization	n/a 1	N	✓	✓	✓	✓	✓
2	Trapped Charge Reclosing	n/a 1	N	✓	✓	✓	-	-
3	Energization	0	Y	✓	✓	✓	✓	✓
4	Trapped Charge Reclosing	0	Y	✓	✓	✓	-	-
5	Energization	2000	Y	✓	✓	✓	-	-
6	Trapped Charge Reclosing	2000	Y	✓	✓	✓	-	-

¹ Altitude not an important parameter for these cases since corona losses were not enabled.

. A batch of 500 simulations was run for each case resulting in a total dataset of 9500 simulations. Simulations consisted of three-phase energization transients and three-phase trapped charge reclosing transients. These are common transients for transmission line insulation studies [31]. The latter is somewhat academic for 500 kV transmission lines as many utilities use single-pole reclosing to improve system stability [32]. Single-pole reclosing reduces the chance of a trapped charge situation. However, trapped charge cases represent a reasonable upper bound for transient overvoltage severity. The trapped charge condition assumed ±1.0 per unit voltage on each of phases A and C, prior to circuit breaker closing. Phase B was assumed to have been the faulted phase and, therefore, had its initial voltage set to zero. Simulations assumed that phase B fault had successfully cleared before reclosing.

Breaker closing, whether for energization or reclose, incorporated random timing. The first pole to operate was selected at random, with each phase having equal probability of being first. The point-on-wave at which the first breaker closes was selected with a uniform probability distribution from 0 to 360 degrees. The delays until close of each of the remaining breakers were based on a Gaussian distribution with a standard deviation ($\sigma )$ of 1.33 ms (3$\sigma$= 4 ms).

The breaker close timing distribution was prepared independently of the simulations. The same distribution was “played back” for each simulation batch in order to improve comparison of study variables. Input timing parameters were recorded for each simulation such that results of interest from later data analysis could be recreated as needed. Circuit breaker pre-strike and re-strike were neglected. The total simulation time for each run was set to 4 ms plus at least 6 times the wave travel time for one length of the line. This ensured sufficient time to record the highest peak voltage for each transient. It is possible that spurious peaks might have occurred later than the maximum simulation time, but an examination of the time of each observed peak indicates that such an occurrence would be very rare and unlikely to have a significant impact on statistical results. The 500 and 800 km cases were included for academic purposes. Energizing 500 kV lines this long would not normally be done without surge mitigation measures.

3. Simulation Method and Model Validation

The FDTD approach was selected because of its inherent ability to handle broadband signals and frequency-dependent and nonlinear components [33]. Specifically, a spatially one-dimensional, constant-parameter model was used. Detailed coverage of this method is beyond the scope of the present paper. Here, the basics of the method are summarized and differences in the approach compared to other implementations are highlighted. The reader is referred to the works of Celozzi, Rachidi, Paul, Kunz, and others for details regarding the implementation of the method [33,34,35,36,37,38,39].

The basic idea of the approach is illustrated in Figure 1. Space and time are discretized, and electric and magnetic field points are offset from each other in both space and time. In the case of constant-parameter models, voltages and currents are surrogates for electric and magnetic fields. Each conductor is treated as a spatially one-dimensional problem. Space and time dependencies are established through discretization of the telegrapher equations with interaction between adjacent conductors modeled via mutual impedance terms. The resulting “update” equations for voltages and currents at each node are iteratively calculated with constraints imposed by boundary conditions at the line terminals.

Using the notation in [36], an example voltage update equation is shown in (2). Included below the equation are the dimensions of the matrices for the three-conductor case. A detailed current update equation is not shown here. The summary equation in (3) shows that current is a function of the adjacent voltage terms from the previous time step and all past currents for that node (via the convolution term with the transient impedance). These dependencies are illustrated in Figure 1 at the circled nodes.

$$\left[V_k^{n+1}\right]=\left[V_k^n\right]-{\left(\frac{\Delta z}{\Delta t}\left[c\right]\right)}^{-1}\left(\left[I_k^{n+1/2}\right]-\left[I_{k-1}^{n+1/2}\right]\right)$$

(2)

$$\left[3\times 1\right] \left[3\times 1\right] \left[3\times 3\right] \left[3\times 1\right] \left[3\times 1\right]$$

$$\left[I_k^{n+3/2}\right]=f\left(\left[V_{k+1}^{n+1}\right],\left[V_k^{n+1}\right],\left[Z\left(n\right)\right]*\left[I\left(n\right)\right]\right)$$

(3)

$$\left[3\times 1\right] \left[3\times 1\right] \left[3\times 3\right] \left[3\times 3\right] \left[3\times 1\right]$$

where $c$ is the capacitance matrix, $Z\left(n\right)$ is the time dependent impedance, and * is the convolution operator.

The next subsections list specific FDTD modeling challenges encountered in the research and describe the approaches used to mitigate them.

3.1. Computation Speed

The target dataset of 9500 simulations meant that the model had to be as efficient as possible in order to complete batches in a reasonable amount of time. This was particularly important for long line cases, where higher node counts and longer simulation times presented greater computational burden. This was accomplished through the selected computing platform and the simplifications discussed below. In the end, the total time to simulate all batches for the 50 through 250 km cases was about three hours. Simulation of the limited case set for the 500 and 800 km cases took a total of about 3 h. These times are for a business laptop with Intel^® Core™ i7 8th generation processor.

3.1.1. Computing Platform

The Julia programming language was ultimately selected as the platform for the model. Julia is a high-level language specifically developed for scientific computing and large-scale linear algebra operations [40]. Julia is a compiled language which means that the computation speed for carefully implemented models can approach that of system level languages such as C++. The FDTD approach is very iterative, so fast loop handling for large arrays containing 3 × 3 matrices was achieved through Julia’s StaticArrays.jl package [41].

3.1.2. Limited Conductor Count

As a tradeoff to maximize speed, the explicit conductor count was limited to three, with ground being implicitly accounted for in the impedance matrices. Shield wires were neglected. The experience of the authors is consistent with [42] in that shield wires have a minor impact on transient results. Shield wires should be included in actual design studies; however, doing so was not necessary for this research. The triple bundle conductors were represented by an equivalent single conductor for each phase. This is discussed further in a subsequent section.

3.1.3. Recursive Convolution

A critical factor for speed is computational burden. Conductor internal impedance and ground impedance are functions of frequency. Hence, in the time domain, they are modeled as transient impedances and are included in a numerical convolution term with the conductor current. As illustrated in Figure 1 and Equation (2), convolution at any given point in the FDTD model requires all past values for the current at each node. Therefore, as a simulation progresses, the number of floating-point operations to compute voltage or current at each spatial node also increases. In [34], Celozzi applies a recursive convolution technique to the FDTD method. In this technique, transient impedances are each approximated as a sum of exponential terms. Each term has unique coefficients. The recursive property of exponential functions allows computation at a given spatial node to be represented as a function of only the previous time step, minimizing storage requirements and floating-point operations. It is important to note that each element of the 3 × 3 earth return matrix is a time varying impedance. Asymmetry in phase conductor geometry means that each of the nine impedance terms in the matrix requires its own set of exponential coefficients. Since phase conductors are nearly always the same size and type, a single set of exponential terms can be used for the internal impedance terms for each conductor.

The number of terms in the exponential approximation depends on the desired accuracy and the method used to identify coefficients for each exponential term. In [36], Paul compares three methods for calculating exponential coefficients. The Matrix Pencil method is shown to have the best accuracy with fewest terms. Using the Matrix Pencil algorithm presented in [43], the transient impedances for the present study were suitably represented by a sum of 5 terms for the conductor internal impedance and 2 terms for each of the 9 earth return impedances.

The conductor internal transient impedance terms are calculated using equations in [36]. The ground impedance terms are based on equations in [44]. Within a given time step, the transient impedances are assumed constant, and have the average value for that time step. It is noted here that Tossani derives the full Sunde expression for ground impedance [44]. The Carson approximation would normally be adequate for frequencies in the switching surge range, but the model was also developed for other purposes. Since the earth return impedance drops so rapidly, the zero-time point becomes an important contributor to the convolution integral, so it is important to be as accurate as possible when establishing ground impedance for the first time step.

3.2. Distributed Dynamic Corona Model for Bundled Conductors

The next modeling challenge relates to implementing distributed corona losses. The approach is an adaptation of a distributed dynamic corona model discussed in [35,45,46]. The original process in the literature focuses on the case of one conductor per phase, and is summarized in the next subsection. A discussion of adaptations necessary for the present research follows.

3.2.1. Summary of Dynamic Corona Capacitance from Literature

The dynamic corona capacitance calculation consists of 4 primary steps:

• First, calculate the corona onset gradient for positive and negative polarity using Peek’s formula (4) with atmospheric correction via (5) [46,47]. These onset gradients are the conductor surface electric fields corresponding to corona onset.

  $$E_c=m ·E_0·\delta ·f·\left(1+\frac{K}{\sqrt{\delta ·r_0}}\right)$$

  (4)

  $$\delta =\frac{P}{760}·\frac{293}{273+T}=e^{\left(\raisebox{1ex}{$-Alt$}\!\left/ \!\raisebox{-1ex}{$8150$}\right.\right)}·\frac{293}{273+T}$$

  (5)

  where $E_c$ is the corona onset gradient, $r_0$ is the conductor radius, $m$ is a surface irregularity factor, $E_0$ is a reference electric field, K is an empirical constant,$\delta$ is the atmospheric correction factor, $f$ is a constant accounting for polarity, $P$ is the pressure in Torr, $T$ is the temperature in °C, and $Alt$ is the altitude above sea level in meters.
• Second, calculate the corona onset voltages for positive and negative polarity using (6). These are the conductor voltages corresponding to the positive and negative corona onset gradients from the previous step. This equation can be derived from first principles considering an isolated conductor above a perfect conducting ground plane.

  $$V_c=E_c·r_0\frac{2·h·r_0}{2·h}ln\left(\frac{2·h-r_0}{r_0}\right)$$

  (6)

  where $V_c$ is the corona onset voltage, $E_c$ is the corona onset gradient from the previous step, $r_0$ is the conductor radius, and $h$ is the height of the conductor above ground. Units for $r_0$ and $h$ must be consistent with those of $E_c$.
• Third, during simulation, monitor the transient voltage on each differential segment of each conductor. If the voltage rises above the corona onset voltage, use Equation (7) to calculate an equivalent conductor radius representing a cylinder that encloses the conductor and a region of free charge produced by corona. This equation is derived from first principles assuming an isolated conductor above a perfect conducting ground plane with the assumption of constant electric field ($\alpha ·E_c)$ between the conductor surface out to radius $r_c$ which defines the corona boundary in air (see Figure 2).

  $$r_c=\frac{\left(r_0+\frac{V}{\alpha ·E_c}\right)}{1+\frac{\left(2·h-r_c\right)}{2·h}·ln\left(\frac{2·h-r_c}{r_c}\right)}$$

  (7)

  where $h$ is the height of the conductor above ground, $V$ is the simulated voltage of the conductor segment from the most recent time step, $r_c$ is the equivalent radius of the corona cylinder, $r_0$ is the geometric radius of the actual conductor, $E_c$ is the corona onset gradient, and $\alpha$ is a multiplier (typically about 0.9) which accounts for the fact that after corona onset, and the electric field at the surface of the conductor drops slightly [46]. Note that (7) requires an iterative solution since $r_c=f\left(r_c\right)$. A simple Gaussian iteration exhibited good convergence with less than 10 iterations.
• Fourth, calculate the total charge on the conductor and in the corona cylinder. Then, calculate the effective capacitance as $c_{eff}=dq/dv\approx \Delta q/\Delta v$. The change in charge and voltage are found by comparing results of the most recent time step with that of the previous time step. This dynamically updated capacitance is calculated for each discrete line segment and each time step as long as the voltage is above the critical voltage and increasing in magnitude. If voltage decreases (even if still above the critical voltage), the capacitance is approximated as the geometric capacitance [19,48].

3.2.2. Process Adaptations for Present Research

The first adaptation of the above process is calculation of the corona onset gradient. Peek’s formula is suited for a single cylindrical conductor but is less accurate for the bundled case. In this research Equations (8)–(10) from [49] were used to obtain a better approximation of corona onset for bundled conductors.

$$E_c=m·31.53·K·\left(1+\frac{A}{K^{\alpha }{n}_s^b{r}_0^c}\right)\left(\frac{1}{\eta }\right)$$

(8)

$$n_s=\frac{1}{1+\sqrt{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$\left(1-\cos \left(\frac{2\pi }{n}\right)\right)$}\right.}}$$

(9)

$$K={\delta }^{1.01}\left(1+0.08\left({\left(\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$11$}\right.\right)}^{0.72}-1\right)\right)$$

(10)

where $E_c$ is the corona onset gradient in kV_peak/cm (8) ; $m$ is a surface irregularity factor (set to 0.6 for a weathered conductor) ; K is an atmospheric correction factor (10); $A$,$\alpha$,$b$, and$c$ are polarity dependent empirical coefficients where, for positive polarity, $A$ = 0.425,$\alpha$ = 0.40,$b$ = 0.26, and $c$ = 0.43, and for negative polarity,$A$ = 0.375,$\alpha$ = 0.49,$b$ = 0.30, and$c$ = 0.45; $\eta$ (Eta) is a field enhancement coefficient accounting for stranding ($\eta \approx 1.38 for ACSR 1272 Bittern)$; $r_0^{ }$is the subconductor radius in cm; $n_s$ is the stranding ratio (9); $n$ is the number of strands in the outer layer of the subconductor, $\delta$ is the relative air density as calculated in (5); and $H$ is the absolute humidity in g/m³ (set equal to 10 for the example 500 kV line in this study).

The next step is to determine the radius of the equivalent conductor that approximates the bundle. Ultimately, since the corona cylinder radius calculations are used to determine a dynamic capacitance, the equivalent conductor should have the same total charge as the bundled conductors for the same voltage. This condition is met by the geometric mean radius of the bundle which can be calculated with Equations (11) and (12) from [27].

$$r_{eq}={\left(N·r_0^{’}·{\left(R_b\right)}^{N-1}\right)}^{\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N$}\right.\right)}$$

(11)

$$R_b=\frac{s}{\left(2·sin\left(\pi /N\right)\right)} \mathrm{for} N>1$$

(12)

where $r_{eq}$ is the equivalent radius giving the same Q–V characteristics as the bundle, $N$ is the number of subconductors, $r_0^{’}$ is the subconductor geometric mean radius ($\approx 0.7788·r_0)$, $R_b$ is the radius of the bundle, and $s$ is the bundle spacing.

While the equivalent conductor has the same charge as the bundle for the same voltage, the surface electric field of the equivalent is lower than that of the bundle. This is illustrated in the finite element electric field simulation in Figure 3. Here, the ACSR 1272 Bittern triple bundle and the equivalent conductor ($r_{eq}$ = 0.1406 m) are both energized to 449 kV_peak (corresponding to 550 kV_rms line-to-line). The resulting total surface charge is 4.52 µC/m for the bundled conductor and 4.51 µC/m for the equivalent conductor (the slight difference is due to rounding error in model inputs). The maximum surface electric field of the bundled conductors is nearly three times that of the equivalent conductor.

This difference must be accounted for when calculating the equivalent corona cylinder radius in Step 3 of the above process. This is achieved by adapting the following equations from [47].

$$\frac{{\rho }_t}{V_{pklg}}=\frac{2·\pi ·{ϵ}_0}{\left(\frac{2·h}{r_{eq}}\right)}$$

(13)

$$E_{avg}=\frac{{\rho }_t}{2·\pi ·{ϵ}_0·r_0}$$

(14)

$$E_{ma\times }=E_{avg}·\left[1+\left(N-1\right)·\frac{r_0}{R_b}\right]$$

(15)

where ${\rho }_t$ is the total charge density of the bundle, $V_{pklg}$ is the peak line-to-ground voltage, ${ϵ}_0$ is the permittivity of free space, $E_{avg}$ is the average subconductor surface electric field for the bundle, $E_{ma\times }$ is the average of the bundle subconductor surface electric field maxima, and all other variables are as defined previously. Let $E_{ma\times }=E_c$ and $V_{pklg}=V_c$ and use (13)–(15) to derive the following expression for voltage in terms of the corona onset gradient calculated in (8) and the equivalent conductor radius calculated in (11).

$$V_c=\frac{E_c·n·r_0·ln\left(\frac{2·h}{r_{eq}}\right)}{\left[1+\left(N-1\right)·\frac{r_0}{R_b}\right]}$$

(16)

Equation (16) replaces (6) in Step 2 of the process. Next, find the equivalent corona onset gradient for the equivalent conductor radius, $r_{eq}$. This can be done by substituting $r_{eq}$ for $r_0$ and $E_{ceq}$ for $E_c$ in (6) and then solving for $E_{ceq}$, giving (17).

$$E_{ceq}=\frac{V_c}{r_{eq}\frac{2·h·r_{eq}}{2·h}ln\left(\frac{2·h-r_{eq}}{r_{eq}}\right)}$$

(17)

Finally, $r_{eq}$ and $E_{ceq}$ are substituted to revise (7) in Step 3, as shown in (18). This equation is solved iteratively in order to find the effective corona cylinder radius as discussed in Step 3.

$$r_c=\frac{\left(r_{eq}+\frac{V}{\alpha ·E_{ceq}}\right)}{1+\frac{\left(2·h-r_c\right)}{2·h}·ln\left(\frac{2·h-r_c}{r_c}\right)}$$

(18)

3.3. Numerical Stability

Improving numerical stability was another challenge experienced during development of the model. Simulations were run in automated batches that could take up to a couple of hours for the longest line cases. It was important that the model be stable and avoid numerical oscillations or other instability that would interrupt simulation flow. Many of the FDTD examples in the literature are excited with controlled waveshapes similar to the standard lighting and switching impulse curves. In this research, the model had to maintain computation through a wide range of switching transients with their attendant wave reflections and coupling between phases and ground. This was particularly challenging considering the dynamic corona capacitance with its nonlinear behavior and on/off thresholds. Three key items helped achieve good numerical stability.

3.3.1. Selection of Spatial Step ($\Delta \mathrm{z}$) and Time Step ($\Delta \mathrm{t}$)

FDTD simulations require that the Courant stability limit be satisfied [38]. This is accomplished by observing the inequality $c\Delta t\le \Delta z$, where $c$ is the speed of light. This is somewhat at odds with the speed requirements and is one disadvantage of the FDTD method. EMTP-type switching surge studies can often use a longer time step. The authors found the following time steps gave stable performance in simulation: 1.67 µs ($\Delta z$ = 500 m) for lines 100 km or less and 2.5 µs ($\Delta z$ = 750 m) for the 250 km and longer lines. From the standpoint of slow-front switching surges, these time and space discretizations allowed the model to effectively approximate a fully distributed approach.

3.3.2. Alternate Dynamic Capacitance Calculation

Recall that the capacitance calculation in Step 4 of the distributed dynamic corona process called for calculation of $c_{eff}=\Delta q/\Delta v$. The $\Delta v$ term is prone to rapid change from small numerical oscillations in the voltage signal. This introduces chatter and greater risk of numerical instability. Therefore, instead of a $\Delta q/\Delta v$ calculation, the equivalent corona cylinder radius from Step 3 of the process is used to directly update diagonal terms of the potential coefficient matrix of each spatial segment. These are then inverted to find the respective capacitance matrix.

3.3.3. Digital Filtering

Even with the alternate capacitance calculation described above, the dynamic capacitance was still prone to chatter caused by rapid changes of the corona onset logic input signals. Assertion and deassertion of corona state resulted in sudden capacitance changes over 50%. The logic inputs consist of a voltage magnitude measurement and the voltage trend (increasing or decreasing). A simple low-pass digital filter was implemented for each input. Rather than using only the voltage from the last time step as the voltage magnitude indication, the voltages of the last two time steps are averaged. Likewise, the voltage trend input looks at the voltage difference over two time steps rather than just one.

Another filter was placed in series with the dynamic capacitance signal output. A single-pole recursive low-pass filter was used [50]. Filter form and parameters are shown in (19).

$$y\left[n\right]=a_0·\times \left[n\right]+b_1·y\left[n-1\right]$$

(19)

where $a_0$ and $b_1$ are filter parameters, $y\left[n\right]$ is the calculated output, $\times \left[n\right]$ is the filter input, and $y\left[n-1\right]$ is the output from the previous time step. Figure 4 shows the response of the filter to a noisy step input illustrative of the possible non-ideal changes in the dynamic capacitance signal.

.

The digital filters described above introduce a small delay. This would be problematic for steep-front waveforms but is tolerable for switching surge type transients.

3.3.4. Arrester Approximation

Typical high-voltage transmission lines are terminated with arresters connected line-to-ground to protect substation equipment from incoming surges. These have a significant impact on the terminal voltage and reflected wave characteristics. The nonlinear volt–current curve for arresters requires an iterative solution which was found to be a source of instability when imposed as a constraint to the FDTD line terminal boundary conditions. Since a detailed model of the arrester itself was unnecessary for the research, an approximation was implemented using an exponential function. The Matrix Pencil method [43] was used to find a best fit with respect to the realistic volt–current curves. Slight manual adjustments were then made based on visual inspection of the curves. Figure 5 shows the volt–current curve of the approximate model.

3.4. Model Validation

The performance of the model was validated by comparing results to EMTP simulations and test data from the literature. EMTP simulations of line energization transients and reclosing transients compared well to FDTD model results with corona losses disabled. A reasonable comparison including corona losses was obtained with the Big Eddy to Chemawa 230 kV EMTP simulations and test data in [22]. Here, the authors compared measured data at line terminals with simulations involving a lumped Siliciu corona model applied at six locations along the line.

Charge–voltage curves from the FDTD model for nonbundled conductors matched those of [46] which were based on test data in [19]. Limited charge–voltage curve data for bundled conductors were found in the literature. However, general trends in curve characteristics between single and bundled conductors from the FDTD model were similar to those in [48]. Finally, simulated results from the FDTD model were consistent with test observations in [51], indicating measured switching overvoltages on a 204 km test line were consistently about 0.1 per unit less than simulated values due to corona loss.

4. Results

Raw results consist of an array of voltages for each of the 500 simulations in a batch. The length of each array corresponds to the number of spatial steps. The voltages in the array are the simulation maxima for a given spatial node. These values are obtained by comparing voltages of all phases over all time steps at each spatial node. Normally, only about 100 or 200 simulations are required in a batch [31], but a larger sample size was used to reduce statistical margin of error. Results presented here include:

• Example plots of the raw output voltage profiles for line energization and trapped charge reclosing cases (Section 4.1). This section includes histograms for a typical cross section of overvoltage data at a given spatial node on the line.
• An example three-phase transient plot showing the impact of corona on a switching transient. The data include the dynamic capacitance response and corresponding charge–voltage curve (Section 4.2).
• Example profiles comparing results with and without distributed corona losses enabled. These plots also show the 98th percentile data calculated from the batch distributions (Section 4.3).
• A brief summary of unexpected results and other cases of interest (Section 4.4).
• Tabular data summarizing the corona impact for all cases analyzed (Section 4.5).
• Tabular data comparing flashover estimates with and without distributed corona losses (Section 4.6).

4.1. Raw Output Data

Figure 6 and Figure 7 show example plots of raw output data from a batch of 500 simulations completed for the 100 km line length case. The compression of the voltage profile toward the open end of the line (right side) in Figure 6 is due to the terminal arresters. Since arresters are connected line-to-ground, the impact on line-to-line voltages is largely negligible. Figure 8 shows the distribution of voltages for both datasets at a point 80 km from the closing end (left side, 0 km position). The line-to-line voltage distributions tend to have a negative (right modal) skew.

4.2. Example Transient Plot of Voltage Attenuation by Distributed Corona Losses

Figure 9 is an example transient case comparing waveforms with and without distributed corona losses. The largest peak on phase C is clearly reduced by corona. It is also of interest to note that the added capacitance from corona introduces a slight delay or shift in the transient. This has implications as discussed in Section 4.4. Figure 10 is a plot of the dynamic capacitance response of phase C from Figure 8. In this example, the dynamic capacitance peaks at nearly 50% above the baseline geometric capacitance. Figure 11 shows the charge–voltage curve for the phase C peak. Energy loss is proportional to the area enclosed by the curve [19]. The curve does not intersect the origin due to the initial trapped charge conditions of the simulation.

4.3. Example Overvoltage Profiles Comparing Corona Impacts (Maximum and 98th Percentile Data)

Figure 12 and Figure 13 are examples of the overvoltage profile data produced in the simulations. Line-to-ground and line-to-line voltage profiles such as these were generated for each case listed in Table 2. The plots show overall maxima (100th percentile) and 98th percentile data. The latter is useful for filtering out spurious peaks of low statistical significance. Results show substantial variation in the shape of the maximum line-to-ground voltage profiles. As seen in comparing Figure 12 and Figure 13, results of line-to-ground cases tend to be less variable. This is due, in part, to capacitive coupling between phases. When the voltage changes sharply on one phase, coupling causes a similar change in adjacent phases. This reduces the differential change between conductors.

It is clear in Figure 12 that corona losses have a notable impact on transient overvoltage severity for line-to-ground voltages. The fact that the corona impact is also prominent in the 98th percentile profile over most of the length of the line indicates the impact is not limited to just the highest peaks. Line-to-line voltage profiles show less reduction due to corona. This is reasonable as separate phases will largely be in different states of corona severity, especially considering the underlying power frequency voltage on which the transients are superimposed.

Results from several cases highlight the importance of resolving voltage profiles with high spatial resolution, with Figure 14 being one example. The solid vertical lines represent a possible approach for selection of terminal and intermediate measurement points in an EMTP study to create overvoltage distributions for insulation design. The dashed line shows the expected voltage profile for such an approach. It is evident that the result would miss the highest voltages and underreport flashover risk. This would still be true (though to a lesser extent) if additional probes are placed at the midpoints between locations already specified in the plot. Targeted placement of a small number of probes could provide an adequate approximation, but without a detailed profile, the optimal locations are unknown. Most profiles in this study indicated that a higher probe density should be used for the last quarter of the line if line-to-ground voltages are being measured. Further analysis of such data could provide better rule-of-thumb guidance for probe placement in EMTP-type studies.

4.4. Cases of Interest

One unexpected result came from the study. In Figure 15, the maximum line-to-line overvoltage for a portion of the profile was worse with corona losses than it was without. Further investigation showed that the result was legitimate from a simulation standpoint, and not the consequence of numerical instability. In these cases, the dynamic corona capacitance induced a delay such that a negative-going transient on one phase overlapped with a positive-going impulse on another phase. The overlap did not occur in the case with corona disabled. It is clear from Figure 15 that such an occurrence is very rare, since there is no accompanying rise in the 98th percentile data. A similar result was observed in the line-to-line voltage profile of the higher altitude 100 km case.

Other results of interest in Figure 16 include overvoltage profiles for the 500 and 800 km cases. Transient voltages generally increase with line length. The results illustrate why long extra-high-voltage (EHV) and ultra-high-voltage (UHV) lines are rarely energized without pre-insertion resistors or controlled closing schemes.

4.5. Tabular Summary of Voltage Attenuation by Distributed Corona Losses

Table 3. Percent reduction of maximum line-to-ground voltages due to corona (percent of nominal).

	Line Energization		Trapped Charge Reclosing
Length (km)	Alt. = 0 m	Alt. = 2000 m	Alt. = 0 m	Alt. = 2000 m
50	5.44%	8.08%	7.58%	10.01%
100	6.71%	9.89%	8.90%	11.41%
250	7.49%	10.32%	12.55%	16.28%
500	9.82%	-	-	-
800	17.10%	-	-	-

,

Table 4. Percent reduction of maximum line-to-line voltages due to corona (percent of nominal).

	Line Energization		Trapped Charge Reclosing
Length (km)	Alt. = 0 m	Alt. = 2000 m	Alt. = 0 m	Alt. = 2000 m
50	2.35%	3.93%	6.60%	6.94%
100	3.20%	5.39%	5.44%	−5.43% 1
250	4.20%	6.64%	7.97%	−5.08% 1
500	9.02%	-	-	-
800	15.16%	-	-	-

¹ Negative due to the first case of interest discussed in Section 4.4.

,

Table 5. Percent reduction of 98th percentile line-to-ground voltages due to corona (percent of nominal).

	Line Energization		Trapped Charge Reclosing
Length (km)	Alt. = 0 m	Alt. = 2000 m	Alt. = 0 m	Alt. = 2000 m
50	4.32%	6.51%	6.72%	8.96%
100	4.63%	7.06%	9.59%	12.07%
250	5.62%	8.37%	10.58%	14.07%
500	8.20%	-	-	-
800	9.27%	-	-	-

and

Table 6. Percent reduction of 98th percentile line-to-line voltages due to corona (percent of nominal).

	Line Energization		Trapped Charge Reclosing
Length (km)	Alt. = 0 m	Alt. = 2000 m	Alt. = 0 m	Alt. = 2000 m
50	2.33%	4.04%	1.13%	2.78%
100	2.82%	4.76%	5.92%	7.10%
250	3.70%	6.08%	6.37%	8.07%
500	5.99%	-	-	-
800	11.42%	-	-	-

summarize the overall impact of corona losses on transient overvoltages. The values in these tables are in percent of the nominal line-to-ground voltage of 408.248 kV_peak. They are obtained by averaging the difference between voltage profiles generated with and without corona losses then dividing by the nominal voltage and converting to percent. Certain profile segments were excluded from the calculation. For example, in many cases, the steep rise in voltage near the closing end of the line (at 0 miles) would tend to exaggerate the influence of the corona losses. In general, impacts are high enough to be worth considering as part of the design process for HSIL applications.

4.6. Estimated Impact of Distributed Corona Losses on Switching Surge Flashover Rate

Table 7. Estimated flashover comparison with and without distributed corona losses (switching surge flashover rate (SSFOR)—normalized to cases without corona losses).

	Line	Operation	Flashover Path 1	Altitude	CFO 2 (kV)	Estimated SSFOR 3
						No Corona Losses	With Corona Losses
1	50 km	Energization	LG	0	946	1.0	0.46
2	50 km	Reclose	LG	2000	1206	1.0	0.34
3	50 km	Reclose	LL	2000	1756	1.0	0.83
4	100 km	Energization	LG	0	1007	1.0	0.42
5	100 km	Reclose	LG	2000	1277	1.0	0.22
6	100 km	Reclose	LL	2000	1802	1.0	0.68
7	250 km	Energization	LG	0	1076	1.0	0.35
8	250 km	Reclose	LG	2000	1438	1.0	0.20
9	250 km	Reclose	LL	2000	1837	1.0	0.53

¹ Line-to-line or line-to-ground flashover path. ² Estimated CFO to give an SSFOR of 1.0 for the case without distributed corona losses. ³ Flashovers per 100 operations based on the specified CFO.

shows results of calculations to estimate the impact of corona losses on flashover rates for the example 500 kV line cases. The result for each row was obtained from analysis of 500 voltage points selected at a location near an overvoltage profile maximum. Switching surge flashover rates were calculated using the methods described in [24]. The first step was to calculate the critical flashover voltage (CFO) that would give a baseline flashover rate of 1.0 flashovers per 100 switching operations for the case without corona. The resulting CFO was then applied to the overvoltage distribution with corona losses. The reductions are noteworthy. However, any reduction in strike distance or structure clearance will depend on the specific situation and may be negligible for some design cases.

5. Conclusions

The multiconductor FDTD method with distributed dynamic capacitance was successfully extended to the case of bundled phase conductors. The efficiency and practicality of the model was demonstrated through high-volume statistical simulation of switching surges on an example 500 kV transmission line of varying length.

Detailed switching overvoltage profiles obtained through distributed techniques, such as the FDTD method, provide information that could be important for refining the design of high-voltage insulation. The cases in this work focused on energization and trapped charge reclosing, but the method could also be used to gain additional insight for switching surges limited by pre-insertion resistors or controlled closing schemes. Results illustrate how simplified overvoltage profiles could underreport flashover probability.

Corona losses have a notable impact on transient overvoltages, particularly line-to-ground exposures. This concept has been known for many years, but it has historically been difficult to quantify the impact because of difficulty in modeling this nonlinear phenomenon in a way that is practical for high-volume simulation. The research demonstrated that such an approach is practical with modern computing capability, reasonable simplifications, and the application of techniques such as digital filtering to improve simulation stability.

Detailed overvoltage profiles that account for voltage attenuation by distributed corona losses provide valuable information for projects seeking to carefully quantify insulation design margins for optimization of transmission capacity in HSIL applications.