Charging and Discharging Current Measurements and Impact of Polarization Dynamics on Electric Field Modeling in HVDC Paper-Insulated Cables

Accurate modeling and simulation of electric field transients in HVDC cables is an important support to optimize insulation system design and to evaluate the influence of voltage transients and steady-state conditions on accelerated ageing mechanisms and insulation reliability. Traditionally, field models considering time-independent permittivity and conductivity are used, but this approach neglects polarization mechanisms and charge trapping-detrapping phenomena. This article includes polarization dynamics in the field model and shows that its impact on transient electric field simulations in HVDC paper-insulated cables can be significant. A method is presented to infer the model parameters from experimental polarization and depolarization current measurements.


I. INTRODUCTION
High Voltage Direct Current (HVDC) systems constitute already an important pillar of electric power grids, and their presence is expected to grow massively in the future. A fundamental element of HVDC systems is insulation, which is the weakest component of any electric system and it must be properly managed to prevent fast ageing and premature failures. This paper focuses on HVDC cables, generally employed for long distance power transmission.
There are two major factors causing extrinsic, accelerated electrical insulation ageing under DC supply [1], that is, space charge (SC) and partial discharges (PD). The former is responsible for the distortion of the Laplacian electric field that may lead to considerable local field magnification, increasing unpredictably electrical stress and, thus, reducing life [2]- [4]. PD take place in defects distributed inside insulation and containing low-density matter. They occur when the nominal voltage in the defect is higher than the inception voltage [5]- [7], which is the minimum voltage required to trigger a PD in a defect. Their impact on ageing comes from The associate editor coordinating the review of this manuscript and approving it for publication was Ali Raza . their magnitude and repetition rate, that is, the number of discharges in the unit of time.
In order to infer intrinsic and extrinsic accelerated aging mechanism [8], accurate simulation of electric field inside insulation, both in steady state and transient conditions, is necessary. Indeed, stress magnitude and distribution may change considerably from transient conditions, when it is mostly driven by permittivity, to steady-stated DC, when it is driven by conductivity [9]. This may impact both the PD inception voltage and the cumulative (intrinsic) aging process through the change of the voltage endurance coefficient, n, of the life-line model [3], [8] being E the electrical stress, L life, ES 0 the reference electric stress (generally close to the electric strength) and t 0 the relevant failure time at applied field E = ES 0 . Equation (1) provides a straight life-line in a log − log coordinate system (log(E) vs log(L)). It is noteworthy that n is higher in DC than in AC, and much lower in the presence of PD [9], thus the life reduction during transient, especially if in the presence of PD, may VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ be not negligible compared to the design DC life conceived for pure steady-state conditions. It must be also noted that voltage (and load) transients will occur even often in DC cables, due e.g., to voltage polarity inversions that are needed to change power flow direction, cable energization, ripple and harmonics superimposed to DC. Traditional electric field simulation approaches [2], [10], [11] are based on the following quasi-electrostatic model: where ρ is the charge density, ε is the constant dielectric permittivity and σ is the temperature-and field-dependent electric conductivity responsible for SC accumulation under DC excitation. This model has been used successfully to compute SC and electric field profiles in HVDC cables, but it is based on static (i.e. non time-dependent) permittivity and conductivity. Actually, they both vary with time, due to polarization mechanisms (ε) [12], and charge trappingdetrapping phenomena (σ ).
The focus of the work presented in the following is, indeed, to include the polarization dynamic, that is, a time dependence of permittivity, in the electric field modelling, and show how it can impact on the simulation accuracy of electric field transients. The model parameters are derived from experimental measurements of polarization and depolarization currents performed at different field and temperatures on specimens of dielectric material.
Even though the HVDC cable insulation technology has been moving recently towards solid and homogeneous dielectrics [13], [14], oil-paper and laminated insulation has been and is still used, especially at the largest values of supply voltage. Therefore, this paper deals with the mature impregnated oil-paper insulation technology. This choice is not only supported by the existence of a large number of oil-paper insulated HVDC cables installed all around the world that are experiencing operating conditions they were not designed for [15] (such as high voltage and load transient rate), but also by the availability of previous studies on the modelling of dielectric polarization mechanisms in paper-oil insulation which can be used as a reference for this investigation [16].

II. MODEL
The dielectric displacement field D is related to the electric field E via where ε 0 is the vacuum electric permittivity and P is the polarization vector. Within the theory of dielectrics [11], P can be expressed as the superposition of a certain number of polarization mechanisms that respond to the application of an electric field with different characteristic times: P ∞ accounts for polarization mechanisms at very high frequencies (as atomic and electronic polarization) and it is related to E by where χ ∞ is the high-frequency dielectric susceptibility. Due to technical reasons, from DC current measurements it is not possible to retrieve information on polarization mechanisms with characteristic times lower than few seconds, so that they can be incorporated in P ∞ . In this way the electric displacement field becomes For the remaining polarization processes we assume a predominant Debye-type relaxation model, where χ k is the process susceptibility and τ k is its time constant.
The dependence on field and temperature of the electric conductivity in (2) is modelled by an empirical law derived from Arrhenius' model of chemical reaction rate [2], [16], [17]: The electro-quasistatic model can then be written as where D is given by (6). Note that if the ''slow'' polarization processes are negligible, i.e.
then D ε ∞ E, the second equation in (9) becomes irrelevant and the model reduces to (2).

III. MODEL PARAMETERS DERIVATION
The model parameters are permittivity ε ∞ , conductivity σ , susceptibilities χ k and time constants τ k . They can be inferred from conduction current measurements performed on flat specimens of the material under investigation: they are energized by a step DC voltage, which is kept until the measured polarization current reaches steady state; then they are shortcircuited, and the depolarization current is measured. The steady-state polarization current is used to retrieve the conductivity parameters: measurements at different temperatures and electric fields are necessary to estimate σ 0 , α and β in (4). The depolarization current helps to get information on the nature of the polarization processes and of the trapping dynamic [18]. Permittivity ε ∞ can be obtained from high frequency AC measurements [12], e.g., 50 Hz for the purpose of this study. If the specimen is sufficiently thin, e.g., 1 mm or less, and it can be assumed that no space charge is accumulated within it (the electric field is below the threshold for space charge accumulation [19]), the specimen can be modelled with the equivalent circuit [16] depicted in Fig. 1. Assuming that the tested specimen has thickness d and area A, its equivalent circuit components are related to the model parameters as: These relationships constitute the link between the parameters of the circuit in Fig. 1 and the ones of the differential model (9). Once the circuital parameters are identified, the corresponding differential ones are retrieved and used to perform the transient electric field simulations with (9). The flowchart in Fig. 2 helps to visualize this process. If a step DC voltage of amplitude V DC is applied between nodes a and b in the above circuit, the resulting polarization current is, for t > 0, if after a time t pol the voltage is removed and the sample is short-circuited, the depolarization current in absolute value is given by where the expression is written assuming that t = 0 is the instant of voltage removal, and It should be noted that the step DC voltage and the short-circuit are actually exponential transients with the duration of few tens of seconds. As we already stated in the previous section however, we are looking for long-range polarization mechanisms that manifest themselves when the source voltage is already at DC, thus the assumption of instantaneous commutation seems to be acceptable.
Equations (12) and (13) indicate that C k and τ k must be inferred from the depolarization current measurement. The method is illustrated with the support of Fig. 3. Taking the logarithm of current i k (t), defined in (12), we have: with ξ = ln t. Let i m (t) be the measured depolarization current; since it is a discharge current, we can expect it to be monotonically decreasing as a function of time, therefore the derivative of I m (t) = ln i m (t) in a given time instantt k can be written as with m > 0. A simple way to derive the time constant τ k comes from the condition which becomes, from (14) and (15): The time instantt k is arbitrary. Once the N time constants τ k have been estimated, capacitances C k can be determined by imposing, for eacht k , the following condition: Using (12) and (13), current i dep t k in (18) can be written as which yields the linear system where The implementation of this method is iterative. One starts to choose a timet N which is close to the maximum measured depolarization time and computes τ N and C N . Thent N −1 <t N is picked in such a way that i dep (t) approximates the measured current; note that C N and C N −1 are computed by solving (19), so the value of C N changes from the initial step. The process is iterated until time instantst 1 , . . . ,t N are got and the total depolarization current i dep (t) is a good approximation of the measured one. Some remarks have to be made: • the method is empirical, but the advantage is that only time instantst k are empirically selected; • this method can fit the measured current with relatively few time constants, which means that the main polarization processes can be successfully identified: this is much better than what is usually achieved with simpler interpolation methods [20]- [22], which require a great number of time constants that have not any physical meaning; • as demonstrated in the following very good fittings can be obtained with the proposed method, but the outcome may also be used to perform a non-linear least squares optimization that helps improving the quality of the result [23]. An experimental depolarization current, measured in [24], from a paper-oil sample is fitted in Fig. 4, to show the effectiveness of the method. The experiment was performed at V DC = 24 kV, on a specimen with d = 0.8 mm and A = 120 × 60 mm 2 . In Table 1 the fitting parameters are summarized. It is worth noting that the highest susceptivity is χ 3 , which is the one associated to the longest time constant: indeed, greater dipoles (hence electric permittivities) are associated to the slower polarization mechanisms, such as Maxwell-Wagner-Sillars polarization, or electrode polarization [12].
Fitting parameter values for the experimental depolarization current measured in [24].  To verify the validity of the model presented in Fig. 1 we should be able to describe the experimental polarization current using the parameters in Table 1. Resistance R is obtained from the experimental steady-state polarization current I pol,∞ , The fitting of the proposed model to the measured polarization current from [24] is shown in Fig. 5. The initial transient cannot be modelled properly because experimental data in the first 10 seconds are not available [24], but for longer times the fitting is very good.

IV. SIMULATIONS AND DISCUSSION
In this section we show that polarization processes in the differential model can impact the transient simulation of a typical HVDC paper-insulated cable, whose parameters are given in Table 2. The details of the numerical implementation will be given in a separate article [25].
The insulating material is that discussed in the previous section. From the experimental data in [24] it comes out that ε ∞ = 4ε 0 . Paper measurements in [24] are performed only at room temperature, T ref = 20 • C, thus for the temperature coefficient in (8) we use a typical value for paper from [16], α = 1.05 × 10 4 K. The average insulation electric field, from 45158 VOLUME 9, 2021  the cable parameters in Table 2, is V n /d ins = 12.5 kV/mm. This field is sufficiently small, and oil-paper conductivity is large enough, to be able to neglect the dependence of conductivity on field [2]. In conclusion, the expression of the conductivity adopted in the simulations is: 1.0194 × 10 −14 S/m. In general capacitances C k and time constants τ k may be temperaturedependent too, but without measurements at temperatures other than T ref we must consider them constant. The external temperature is fixed at 20 • C, while different internal temperatures are evaluated (30 • C, 40 • C and 50 • C), corresponding to different loading conditions. The cable is supplied by the voltage waveform depicted in Fig. 6. Field behavior during energization, inversion and discharge transients is simulated. To simulate realistic switching, the voltage jumps in Fig. 5 are implemented by an exponential law, i.e.
where τ g = 0.5 s and V = +V n , −2V n , +V n in the three switching times t 0 = 0, 5·10 4 , 10 5 [s]. The insulation electric field profiles at different instants are shown in Figs. 7 to 10, to highlight the differences between the ''improved'' model (9) and the ''standard'' one (2). Fig. 7 shows the resistive steady-state field at the end of the energization transient. The two models predict the same profile because the steady-state field is determined only by the electric conductivity σ and the cable geometry. Figs. 8, 9 and 10 show the electric field profiles in the three simulated transients. In all cases, one can expect that the transient predicted by the improved model is slower than the one predicted by the standard model: this is due to the polarization processes that are considered by (9).
To summarize, let be the difference between the electric field obtained with the improved (E imp ) and standard (E std ) models as a function of radial position and time, divided by the average electric field, E avg = 12.5 kV/mm. The maximum amplitude of E(x, t)   is depicted in Fig. 11 for each time instant, The figure shows δ E (t) at various temperature drops T = T int − T ext , with T ext = 20 • C. As can be seen, its values increase with temperature. Moreover, during the inversion transient δ E (t) is approximately twice than in the energization and discharge ones: this is due to the voltage step of the commutation, which is 2V n in case of polarity inversion and V n in the other two stages. The space and time dependence of E(x, t) is summarized in Figs. 12 and 13, where the cases with T = 10 • C and T = 30 • C are shown, respectively.  In the blue regions E < 0, thus E std > E imp , while in the violet ones E > 0, which means E imp > E std . The figures show that the difference between E imp and E std , regardless of the sign, intensifies immediately after a commutation and it is spread all over the insulating region, not just confined in narrow areas. The colored regions in Fig. 13 are more intense than the ones in Fig. 12, meaning that E becomes larger as temperature drop T increases.
In the white regions E = 0, therefore E std = E imp . This situation is reached when both models are at steady-state and the electric field profile for both E std and E imp is determined by the condition which holds for homogeneous materials.

V. CONCLUSION
This article proposes an improved approach, compared to the conventional one, to model the electric field distribution in HVDC oil-paper cables during voltage transients. The novelty here is to incorporate the dynamics of dielectric polarization mechanisms into the fundamental equations and relevant simulations. A simple empirical method is proposed, which is based on a circuital model, to derive the differential model parameters from polarization and depolarization current measurements performed on appropriate insulation specimens. It is shown that the fitting of the model to experimental data can be accurate and that the most significant polarization processes under DC excitation can be successfully identified. The simulations performed on a typical paper insulated cable show that the impact of polarization dynamics on the electric field transient can be significant, in terms of both the transient duration and of its magnitude and profile across cable insulation. It is observed also that differences between the conventional model and the proposed one increase with voltage amplitude and temperature gradient. Future studies will focus on solid polymeric dielectrics, where differences with respect to the conventional model are expected to be larger. PAOLO SERI (Member, IEEE) was born in Macerata, Italy, June 1986. He received the M.Sc. degree in energy engineering, in 2012, and the Ph.D. degree in electrical engineering, in 2016, both from the University of Bologna. Since 2017, he has been part of the Laboratory of Innovative Materials for Electrical Systems (LIMES), University of Bologna, as a Research Fellow. He is currently an Assistant Professor with the Department of Electrical, Electronic and Information Engineering, DEI, Bologna University. He is also working on the topics of HVDC cables design, partial discharge detection and modelling, and characterization of dielectric materials.
GIAN CARLO MONTANARI is currently a Research Faculty III with the Center for Advanced Power Systems (CAPS), Florida State University, USA, and a Alma Mater Professor with Bologna University, Italy. He has been a Full Professor of electrical technology with the Department of Electrical, Electronic and Information Engineering, University of Bologna, where he is teaching courses on technology, reliability and asset management. Since 1979, he has been working in the field of aging and endurance of insulating materials and systems, diagnostics of electrical systems, asset management, and innovative electrical materials (magnetics, electrets, super-conductors, nano-materials). He has been engaged also in the fields of power quality and energy market, power electronics, reliability and statistics of electrical systems, and smart grid. He is also an author or a coauthor of more than 800 scientific articles. He was the recipient of the several awards, including the IEEE Ziu-Yeda, the Thomas W. Dakin, the Whitehead, the Eric Forster, and the IEC 1906 awards. He was the Founder and the President of spin-off Techimp, established in 1999.