A Review about the Modeling and Simulation of Electro-Quasistatic Fields in HVDC Cable Systems

: In comparison to high-voltage alternating current (HVAC) cable systems, high-voltage direct current (HVDC) systems have several advantages, e.g., the transmitted power or long-distance transmission. The insulating materials feature a non-linear dependency on the electric ﬁeld and the temperature. Applying a constant voltage, space charges accumulate in the insulation and yield a slowly time-varying electric ﬁeld. As a complement to measurements, numerical simulations are used to obtain the electric ﬁeld distribution inside the insulation. The simulation results can be used to design HVDC cable components such that possible failure can be avoided. This work is a review about the simulation of the time-varying electric ﬁeld in HVDC cable components, using conductivity-based cable models. The e ﬀ ective mechanisms and descriptions of charge movement result in di ﬀ erent conductivity models. The corresponding simulation results of the models are compared against measurements and analytic approximations. Di ﬀ erent numerical techniques show variations of the accuracy and the computation time that are compared. Coupled electro-thermal ﬁeld simulations are applied to consider the environment and its e ﬀ ect on the resulting electric ﬁeld distribution. A special case of an electro-quasistatic ﬁeld describes the drying process of soil, resulting from the temperature and electric ﬁeld. The e ﬀ ect of electro-osmosis at HVDC ground electrodes is considered within this model.


Introduction
For the transmission of high electric power over long distances, high-voltage direct current (HVDC) is used due to less losses in comparison to high-voltage alternating current (HVAC) [1]. Vanishing capacitive and inductive losses in HVDC power cables result in no theoretical limit for the cable length. A HVDC cable system consists of the cable itself, the cable joint and the cable termination. The cable joint is used to connect two cable segments. Cable terminations are installed at the end of a cable or connected e.g., to a "gas-insulated substation" (GIS) [1][2][3][4].
Worldwide HVDC projects show that the two main insulation materials are mass-impregnated paper (MI) and extruded materials. Extruded insulations are polymeric materials such as different types of polyethylene (PE), e.g., "cross-linked polyethylene" (XLPE) or "low-density polyethylene" (LDPE). The applications of extruded cables have increased in recent years due to significant advantages in manufacturing, installation, transportation and maintenance [5]. The world's first HVDC cable was put into service in 1954, with a voltage of 80 kV and an insulation of MI. Polymeric insulation materials appeared in the late 1960s and have been used for HVAC cables. In 1999, the first commercial polymeric HVDC cable was put into service, with a voltage of 80 kV [6,7]. Today, the voltage has where ε0 = 8.854 × 10 −12 As/(Vm) is the dielectric constant, εr is the relative permittivity and   ,| | TE  is the electric field E and temperature T dependent electric conductivity, a stationary electric field is obtained at approximately 10•τ. Due to the non-linear electric conductivity, a stationary field may be obtained after hours or weeks [1,3,13,14]. For example, Figure 1a shows the electric field within a MI insulation with conductor temperature T(ri) = Ti = 55 °C, sheath temperature T(ra) = Ta = 35 °C and applied voltage U = 450 kV. A stationary electric field is seen above the black dashed line. The cylindrical geometry is depicted in Figure 1b, where the electric field is computed along the red evaluation line due to its radial dependency only [11,13].
(a) (b) Figure 1. (a) Time-dependent electric field within a mass-impregnated paper (MI) insulation. The stationary electric field is obtained above the black dashed line; (b) geometry of the cylindrical insulation. The electric field is computed along the red evaluation line [11,13].
During the operation time, the electric field varies from a purely capacitive field (t = 0), determined by the permittivity, to a purely resistive field (t ≈ 4 h), determined by the electric conductivity [11]. With accumulated space charges, the electric field decreases in the vicinity of the conductor and increases near the sheath, called the "inversion of field" [8]. Due to the field-changing behavior of accumulated space charges, the electric field may exceed the breakdown strength of the dielectric that can cause reliability problems or failure of HVDC cable components. Furthermore, high electric field values and increased losses within the insulation may result in a thermal breakdown [1,15,16].
In addition to measurements, numerical simulations are used to determine the space charge and electric field distribution inside the insulation. The simulated electric field and charge distribution can be used to prevent HVDC cable components from possible failure.
In the last years, several research and review papers (e.g., [1,4,11,15]) have been published within the field of HVDC cables systems, where the modeling and the numerical simulation of the The stationary electric field is obtained above the black dashed line; (b) geometry of the cylindrical insulation. The electric field is computed along the red evaluation line [11,13].
During the operation time, the electric field varies from a purely capacitive field (t = 0), determined by the permittivity, to a purely resistive field (t ≈ 4 h), determined by the electric conductivity [11]. With accumulated space charges, the electric field decreases in the vicinity of the conductor and increases near the sheath, called the "inversion of field" [8]. Due to the field-changing behavior of accumulated space charges, the electric field may exceed the breakdown strength of the dielectric that can cause reliability problems or failure of HVDC cable components. Furthermore, high electric field values and increased losses within the insulation may result in a thermal breakdown [1,15,16].
In addition to measurements, numerical simulations are used to determine the space charge and electric field distribution inside the insulation. The simulated electric field and charge distribution can be used to prevent HVDC cable components from possible failure.
Energies 2020, 13, 5189 3 of 42 In the last years, several research and review papers (e.g., [1,4,11,15]) have been published within the field of HVDC cables systems, where the modeling and the numerical simulation of the electric field are mainly of minor interest. Different conductivity models and simulation methods are found throughout the literature. This work is a review about the simulation of the time-varying electric field in HVDC cable components, using conductivity-based cable models. Based on the effective mechanisms and descriptions of charge movement, different conductivity models are presented. With a comparison of the simulation results against measurements and analytic solutions, the main results of the electric field in HVDC cable systems are summarized.
The simulation of the electric field under constant stress covers various topics, which makes it necessary to divide this paper into several parts. After the introduction, the next section consists of the charge movement within the insulations and the resulting non-linear electric conductivity. Section 3 gives a short review about the accumulation of charges and their influence on the electric field. The necessary equations and possible pseudo codes to obtain the electric and thermal fields are described in Section 4. Simulation results of the time-varying and the stationary electric fields are presented in Section 5. The results in Section 5 neglect the effect of interfaces or moving charges, which are seen in Section 6. A comparison between different simulation techniques, focusing the simulation time and the accuracy, is presented in Section 7. A special case of a time-varying electric field in the vicinity of ground electrodes is described in Section 8, followed by the conclusions and an outlook into further research areas in Section 9.

High-Voltage Direct Current (HVDC) Cable Insulation Materials and the Non-Linear Electric Conductivity
MI as insulation material consists of high-density paper that is lapped and then impregnated with a high-viscosity mineral oil [5]. XLPE and LDPE are based on PE, where e.g., XLPE is manufactured within an additional cross-linking process. PE is a semicrystalline polymer that consists of amorphous and crystalline regions. Within the crystalline regions, polymer chains contain an orderly configuration and align themselves in crystalline lamellae. During the crystallization, many lamellae may grow from a single nucleus and form a spherulite within which the lamellae are radially aligned. The amorphous regions between the lamellae and spherulites, consists e.g., of chain ends or branches [3,[17][18][19]. The final insulation material consists of the dielectric (PE) and two semiconducting layers (semicon), where all three layers are usually extruded within the same production step [1].
Using the band diagram to describe both materials, single states (traps) form the conduction band and the valence band, resulting in a mean band gap energy W. These traps are the result of oil and cellulose molecules in MI, the semicrystalline structure of PE and different impurities that are seen in both insulations, coming from the manufacturing process. Charges are caught within the traps and facing a potential barrier, when trying to escape. For example, the band diagram of XLPE is seen in Figure 2a, where charges are trapped and released.
The traps are called "physical traps" or "shallow traps" and "chemical traps" or "deep traps", depending on the trap depth. For XLPE, shallow traps have a depth of 10 −2 -1 eV, while the depth of deep traps is >1 eV [17,20,21]. The density of the traps is high near the conduction and valence band and decreases with increasing distance, as illustrated in Figure 2b.
Depending on the trap depth E Trap , charges remain in traps that may range from a few seconds to years. The time t Trap is given by: where k B = 1.38 × 10 −23 J/K is the Boltzmann constant and 1/t Trap,0 ≈ (k B ·T)/h P is the "escape frequency", with h P = 6.626 × 10 −34 Js as the Plank constant [17,20,[22][23][24].
morphology, the local electric field and temperature or the humidity. In literature, the conductivity is commonly given with a dependency on the local electric field and temperature only (see e.g., [14,17,20,23,25]), which are the most important factors [26]. With (2), the electric conductivity is described by:   1 A, 1 1 and similar to (3) and (4), whereas the dependence on the electric field differs [17,20]. In (5), the constants are E A,3 , K 3 and γ 3 that have to be determined by measurements as well. The conductivity models in (3)-(5) are described by analogies to the semiconductor technology. Based on worldwide experience and a narrow range of the operation temperature, the electric conductivity is commonly given in a double exponential form with: where the constants are κ 0 , α and β [22,[30][31][32]. For MI, α ≈ 0.1 • C −1 and β ≈ 0.03 mm/kV [30]. In the case of XLPE, α ≈ 0.1 • C −1 and β ≈ 0.1 mm/kV [31]. Thus, the temperature dependency of both materials is approximately equal, but the electric field dependency is 3.5 times higher in XLPE [22,30,31]. A different formulation of (6) is: where the constants κ 0 and α are equal to (6), but the electric field dependency is determined by the constants E Ref and v. Considering a cable geometry, as depicted in Figure 1b, the approximations E Ref = U·exp(−1)/(r a − r i ) and v = U·β/(r a − r i ), with r i as the radius of the conductor and r a as the radius of the sheath, hold true [25,26,30]. There is still a gap of knowledge to describe the charge transport and the electric conductivity of the insulating material in detail, due to the complex morphology. This results in many different possibilities to describe non-linear variations over temperature and electric field. In [33], further conductivity models are shown. The conductivity models (3)-(7) are valid for high field conduction ( → E > 10 kV/mm), which comes e.g., from injection processes at the electrode-insulation-interface [8,20,34]. Below the electric field threshold (10 kV/mm), the dependency on the electric field is negligible [23]. Using the non-linear electric conductivity, space charge accumulation is described. As mentioned in the next section, the charges are characterized by their position within the insulation, having different effects on the electric field.

Space Charges, Surface Charges and Charge Packets
Space charges accumulate within the material due to injection processes, the ionization of impurities, the polarization of the material, and a varying ratio of permittivity over conductivity [8,35,36]. Homocharges are located near the electrode with the same polarity and heterocharges are located near the electrode with the opposite polarity. The effect of both charge distributions on the resulting electric field, between planar electrodes with distance D and applied voltage U, is depicted in Figure 3. Homocharges decrease the electric field at the electrodes and increase the field within the bulk. Heterocharges increase the field at the electrodes, while it is decreased in the bulk (see Figure 3b,c) [2].
Homocharges are injected charges that are trapped close to the electrode, if charge injection exceeds the drift [2,37]. The presence of some impurities, e.g., cumyl alcohol, a byproduct of the XLPE crosslinking agent dicumyl peroxide (DCP), increase the trap density near the electrodes and supports the accumulation of homocharges [38].
Heterocharges are formed with and without a temperature gradient. In the absence of a temperature gradient, heterocharge accumulation results from the ionization of impurities. A spatial varying distribution of permanent dipole molecules form an inhomogeneous polarization and, thus, heterocharges [36,39]. Charges that have dissociated from molecules move to the electrodes and may also form heterocharges. With a temperature gradient, charges are injected at the high temperature electrode and move to the low temperature electrode, where they might not be extracted [2,40,41].
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies [8,35,36]. Homocharges are located near the electrode with the same polarity and heterocharges are located near the electrode with the opposite polarity. The effect of both charge distributions on the resulting electric field, between planar electrodes with distance D and applied voltage U, is depicted in Figure 3. Homocharges decrease the electric field at the electrodes and increase the field within the bulk. Heterocharges increase the field at the electrodes, while it is decreased in the bulk (see Figure 3b,c [2]. Homocharges are injected charges that are trapped close to the electrode, if charge injection exceeds the drift [2,37]. The presence of some impurities, e.g., cumyl alcohol, a byproduct of the XLPE crosslinking agent dicumyl peroxide (DCP), increase the trap density near the electrodes and supports the accumulation of homocharges [38]. (b) schematic homocharge and heterocharge distribution; (c) effect of the charge distribution on the resulting electric field [2,8].
Surface charges at both electrodes are the result of the polarization of the material under an applied electric field and the presence of bulk space charges [40].
Injected charges may also form charge packets (charge pulses) that move to the opposite electrode. Measured charge packets are mainly found in polymeric insulation materials [42,43]. A description of the physics of charge pulses is reviewed in [43]. According to [44], charge packets are observable at electric fields of 10-50 kV/mm. At electric fields of >40 kV/mm a repetitive injection of charges is seen and several pulses move through the insulation. In [42], measured charge packets are divided into "fast" and "slow" packets. Fast charge packets have a mobility of up to 10 −10 m 2 /(Vs) and a density of 0.07-0.1 C/m 3 . By contrast, slow charge packets have a mobility of 10 −16 -10 −14 m 2 /(Vs) with a density >0.1 C/m 3 [45]. The semicon interface between the polymeric insulation and the electrode shows a blocking behavior for charge packets, resulting in the accumulation of heterocharges and in a decrease of the lifetime [42,46].

Numerical Simulation of Charge Transport and the Electric Field in Direct Current (DC) Cable Insulations
Numerical simulations are a powerful tool to determine the charge transport and the corresponding electric field distribution within the insulation material and furthermore a complement to measurements. Two models are commonly described in literature to compute the charge transport and the electric field in DC insulations, where the corresponding equations are presented in Section 4.1 and the discretization with possible pseudo codes is shown in Section 4.2.

Equations to Simulate the Time-Varying Charge and Electric Field Distribution
The macroscopic description ("conductivity model") of charge accumulation is based on the electric conductivity of the insulation and simulates an average value of positive and negative charge carriers (see e.g., [11,14]). The conduction is described at different temperatures and electric fields, but it does not describe the involved mechanisms (e.g., trapping or detrapping of charges). These mechanisms are modeled within the microscopic description ("bipolar-charge-transport-model") [22,47].
The governing equations for the simulation, using the conductivity model, are the continuity equation, Poisson's equation and Ohm's law, i.e., div(ε 0 ε r grad ϕ) = −ρ, where m is the discrete time index and ∆t is the time step length that has to be lower than ∆t CFL , determined by the Courant-Friedrich-Levy (CFL) criterion, for stability reasons [52]. Similarly as with (11), the time discretization of the heat conduction Equation (16), using an explicit Euler method, yields the update scheme: where M δ is the density matrix, M c is the matrix of the specific heat capacity, M λ is the thermal conductivity matrix, q T = M κ (u T ,Φ)GΦGΦ is the vector of heat sources and b T contains the boundary conditions for the thermal problem [48,49,60,61]. A possible calculation scheme to obtain the electric field and the space charge density, using (8)- (10), is depicted in Figure 4 and using (11) in Figure 5. The computation starts with an assumed vanishing space charge distribution ρ = 0 C/m 3 (q = 0) and a constant temperature equal to the environment temperature T ∞ . The vanishing space charge density results in a purely capacitive field that is determined by the permittivity and the geometry only. The time integration stops if a predefined time t = t END or Φ m+1 − Φ m / Φ m < η is obtained, where η 1 is the stop threshold and Φ m is the absolute value of the vector Φ at the discrete time step t m [11,58,59].
The simulation of an HVDC cable system includes the cable, the cable joint and the cable termination. From a numerical perspective, the cable joint and the cable termination are equally difficult to simulate. Both include different field-grading techniques and they consists of multiple subcomponents, resulting in interfaces between different dielectrics. Thus, the following simulations include cables and cable joints only.
where Mδ is the density matrix, Mc is the matrix of the specific heat capacity, Mλ is the thermal conductivity matrix, qT = Mκ(uT,Φ)GΦGΦ is the vector of heat sources and bT contains the boundary conditions for the thermal problem [48,49,60,61]. A possible calculation scheme to obtain the electric field and the space charge density, using (8)- (10), is depicted in Figure 4 and using (11) in Figure 5.  . Numerical computation scheme, to obtain the electric field and the space charge density, using the space charge-oriented field formulation (8)-(10) [11,58,59].
Energies 2019, 12, x FOR PEER REVIEW 9 of 42 Figure 5. Numerical computation scheme, to obtain the electric field and the space charge density, using the scalar potential field formulation (11) [11,58,59].
The computation starts with an assumed vanishing space charge distribution ρ = 0 C/m 3 (q = 0) and a constant temperature equal to the environment temperature T∞. The vanishing space charge density results in a purely capacitive field that is determined by the permittivity and the geometry only. The time integration stops if a predefined time t = tEND or where η 1 is the stop threshold and m Φ is the absolute value of the vector Φ at the discrete time step t m [11,58,59]. The simulation of an HVDC cable system includes the cable, the cable joint and the cable termination. From a numerical perspective, the cable joint and the cable termination are equally difficult to simulate. Both include different field-grading techniques and they consists of multiple subcomponents, resulting in interfaces between different dielectrics. Thus, the following simulations include cables and cable joints only.

Simulation of the Electric Field and the Space Charge Density within a Cable Insulation
Simulation results of a cable insulation are presented in this section, showing the behavior of the electric field and the space charge density under different thermal stress. In Section 5.1, the time-varying electric field is simulated and approximated, using the analytic solution of the stationary electric field. Furthermore, the stationary electric field and the corresponding charges are Figure 5. Numerical computation scheme, to obtain the electric field and the space charge density, using the scalar potential field formulation (11) [11,58,59].

Simulation of the Electric Field and the Space Charge Density within a Cable Insulation
Simulation results of a cable insulation are presented in this section, showing the behavior of the electric field and the space charge density under different thermal stress. In Section 5.1, the time-varying electric field is simulated and approximated, using the analytic solution of the stationary electric field. Furthermore, the stationary electric field and the corresponding charges are analyzed under different temperature gradients. The computation of the stationary electric field may result in instabilities during the iteration process, as shown in Section 5.2. Considering a varying environment of the cable, varying temperature and electric field distributions are computed and presented in Section 5.3.

Transient and Stationary Electric Field under Various Thermal Stresses
Several authors have studied space charge phenomena in DC power cables, after one of the first works of Lau in 1970. The space charge evolution in time has been described for the first time, using a non-linear electric conductivity, depending on the temperature only [62]. Within Eoll's work in 1975, the stationary electric field is described, considering the dependency on the temperature and on the electric field and the losses within the insulation [30]. In [62] and [30] the electric conductivity is described by (6), where the expression exp β · → E = 1 is used in [62], due to a temperature dependency only. Equation (6) or the approximation (7) are also used in the following years (see e.g., [11,14,25,26,63,64]).
Considering a MI cable insulation, depicted in Figure 1b, the transient electric field stress and space charge distribution during the load of the cable is shown in Figure 6. The electric conductivity is given by (6) and the constants for the computation are summarized in Table 1. A convective heat flux → q is assumed at the outer sheath and described by where α out is the heat transmission coefficient [65][66][67][68]. Using Figure 6, the maximum space charge density ( Figure 6c) and thus, the maximum field inversion (Figure 6b) are obtained for a stationary temperature distribution, showing the maximum temperature gradient (Figure 6a).  Figure 1b, the transient electric field stress and space charge distribution during the load of the cable is shown in Figure 6. The electric conductivity is given by (6) and the constants for the computation are summarized in where αout is the heat transmission coefficient [65][66][67][68]. Using Figure 6, the maximum space charge density ( Figure 6c) and thus, the maximum field inversion (Figure 6b) are obtained for a stationary temperature distribution, showing the maximum temperature gradient (Figure 6a).  In [62], it is assumed that the thermal time constant is low in comparison to the electric time constant (1) and a stationary temperature is already obtained at t = 0 [64]. Thus, in several studies only the stationary temperature distribution is considered [11,14,25,26]. In [62], it is assumed that the thermal time constant is low in comparison to the electric time constant (1) and a stationary temperature is already obtained at t = 0 [64]. Thus, in several studies only the stationary temperature distribution is considered [11,14,25,26].  [11,30,65,66].

Constant
Value (6)) 0.03 mm/kV δ·c p 2.5 × 10 6 J/(m 3 K) r a 42.4 mm v (for (7)) 0.7031 α out 5 W/(K·m 2 ) ε r 3.5 E Ref (for (7)) 8.622 kV/mm --With (7), the stationary electric field is given in a closed analytic form. An analytic derivation of the stationary electric field and charge density starts from Gauss' law (9), where Ohm's law (10) is used to replace the electric field with the electric conductivity and the current density. Applying the chain rule and using the stationary current div → J = 0, one obtains: With (24), a varying electric conductivity, either by a varying electric field and/or temperature or other processes will automatically result in accumulated charges [14,26]. For a radial symmetric electric field E(r), space charge ρ(r) and temperature profile T(r), (24) is reduced to: With the Gauss law (9) to replace the charge density ρ(r) in (25) and an assumption of a constant permittivity (ε r = const.), the stationary electric field is the solution of: With the electric conductivity (7), the electric field is obtained by solving: Considering negligible insulation losses and with constant parameters δ·c p and λ, the stationary temperature is the solution of (16) for ∂/∂t = 0 [69] and is given by: With (28), the solution of (27) is: where C is a constant. Using the definition a b = exp(b·ln(a)), (29) is reduced to: Energies 2020, 13, 5189

of 42
With the applied voltage U, the constant C is: In [3,62], the temperature is described by the thermal conductivity λ of the material and the heat losses of the conductor P V and is given by: With (32) the electric field is: Applying Gauss' law, the corresponding stationary space charge distribution is: Utilizing the stationary electric field E(r) and the initial electric field at t = 0 (E 0 (r) = U·r −1 ·ln(r a /r i ) −1 ), a rough estimation of the transient electric field is obtained. With (11) and Figure 6b, a reasonable assumption for the time dependence of the transient electric field is an exponential function. The stationary electric field E(r) is equivalently given by (30) or (33), respectively. With an exponential approach and the time constant τ in (1), the transient electric field is: With an electric field depending time constant, an approximation for τ is determined with the mean electric field U/(r a − r i ). Using (7), the constants in Table 1 and the temperature distribution (28), with T i = 50 • C and T a = 35 • C, the numerically computed solution and the relative error of the approximation (35) are depicted in Figure 7. The relative error of (35) has the highest values in the vicinity of the conductor and the sheath, during the inversion of the field. The maximum error is about +13.58%. If the electric field approaches the stationary solution, the error decreases and is close to zero.
The stationary electric fields (30) within a MI insulation for a constant conductor temperature of T i = 50 • C and a varying sheath temperature T a are depicted in Figure 8 [11,22]. As depicted in Figure 8a, the inversion of the electric field is seen at temperature gradients higher than 5 • C. Temperature gradients higher than 20 • C result in a sheath electric field that is higher than the maximum electric field at t = 0, which is about 32 kV/mm. The space charge distribution in Figure 8b for a temperature gradient of ∆T = 0 • C results from the dependency of the conductivity on the electric field. At temperature gradients higher than 15 • C the charge density increases towards the sheath, resulting from the constant δ E that is <−1 and the negative exponent of the radius r in (34).
Within the insulation, a homogeneous electric field is desirable. To obtain a homogeneous electric field, δ E needs to vanish, resulting in E(r) = U/(r a − r i ). With (30), δ E vanishes at a temperature gradient T i − T a = ln(r a /r i )/α ≈ 6 • C. Utilizing (33) this temperature gradient corresponds to losses of P V = 2πλ/α ≈ 10.5 W/m, with λ MI = 0.167 W/(Km) [30]. With losses higher than 10.5 W/m, the field inversion increases and with lower losses the electric field is higher at the conductor in comparison to the field at the sheath [70]. The temperature coefficient α is approximately equal for XLPE and MI, but the thermal conductivity of XLPE is higher (λ XLPE = 0.27-0.32 W/(Km)) in comparison to MI. Using an equal XLPE cable insulation, the losses are P V ≈ 17-20 W/m to obtain a homogeneous electric field. Thus, a better heat dissipation of XLPE results in higher transmitted power until the risk of high Energies 2020, 13, 5189 13 of 42 field values at the outer sheath occurs [5]. According to [25,26] the charge density is separated into two parts. One corresponds to the temperature dependency of the electric conductivity (ρ T ) and one corresponds to the electric field dependency (ρ E ). With a given temperature gradient both charge parts are of the opposite sign. With increasing field dependency (v) the effect of field inversion decreases, but the net charge density (|ρ T | + |ρ E |) increases. In comparison to MI, the electric field dependency of XLPE is 3.5 times higher [30,31]. Thus, the net charge density (|ρ T | + |ρ E |) is higher in polymeric cable insulations, whereby XLPE usually operate under lower applied electric fields [2,3].
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies Utilizing the stationary electric field E(r) and the initial electric field at t = 0 (E0(r) = U•r −1 •ln(ra/ri) −1 ), a rough estimation of the transient electric field is obtained. With (11) and Figure 6b, a reasonable assumption for the time dependence of the transient electric field is an exponential function. The stationary electric field E(r) is equivalently given by (30) or (33), respectively. With an exponential approach and the time constant τ in (1), the transient electric field is: With an electric field depending time constant, an approximation for τ is determined with the mean electric field U/(ra − ri). Using (7), the constants in Table 1 and the temperature distribution (28), with Ti = 50 °C and Ta = 35 °C, the numerically computed solution and the relative error of the approximation (35) are depicted in Figure 7. The relative error of (35) has the highest values in the vicinity of the conductor and the sheath, during the inversion of the field. The maximum error is about +13.58%. If the electric field approaches the stationary solution, the error decreases and is close to zero. The stationary electric fields (30) within a MI insulation for a constant conductor temperature of Ti = 50 °C and a varying sheath temperature Ta are depicted in Figure 8 [11,22]. As depicted in Figure  8a, the inversion of the electric field is seen at temperature gradients higher than 5 °C. Temperature gradients higher than 20 °C result in a sheath electric field that is higher than the maximum electric field at t = 0, which is about 32 kV/mm. The space charge distribution in Figure 8b for a temperature gradient of ΔT = 0 °C results from the dependency of the conductivity on the electric field. At temperature gradients higher than 15 °C the charge density increases towards the sheath, resulting from the constant δE that is <−1 and the negative exponent of the radius r in (34).
Within the insulation, a homogeneous electric field is desirable. To obtain a homogeneous electric field, δE needs to vanish, resulting in E(r) = U/(ra − ri). With (30), δE vanishes at a temperature gradient Ti − Ta = ln(ra/ri)/α ≈ 6 °C. Utilizing (33) this temperature gradient corresponds to losses of PV = 2πλ/α ≈ 10.5 W/m, with λMI = 0.167 W/(Km) [30]. With losses higher than 10.5 W/m, the field inversion increases and with lower losses the electric field is higher at the conductor in comparison to the field at the sheath [70]. The temperature coefficient α is approximately equal for XLPE and MI, but the thermal conductivity of XLPE is higher (λXLPE = 0.27-0.32 W/(Km)) in comparison to MI. Using an equal XLPE cable insulation, the losses are PV ≈ 17-20 W/m to obtain a homogeneous electric field. Thus, a better heat dissipation of XLPE results in higher transmitted power until the risk of high field values at the outer sheath occurs [5]. According to [25,26] the charge density is separated into two parts. One corresponds to the temperature dependency of the electric conductivity (ρT) and one corresponds to the electric field dependency (ρE). With a given temperature gradient both charge parts are of the opposite sign. With increasing field dependency (v) the effect of field inversion decreases, but the net charge density (|ρT| + |ρE|) increases. In comparison to MI, the electric field dependency of XLPE is 3.5 times higher [30,31]. Thus, the net charge density (|ρT| + |ρE|) is higher in polymeric cable insulations, whereby XLPE usually operate under lower applied electric fields [2,3].

Fast Calculation of the Steady State Charge Distribution
In the case of overvoltage (e.g., impulse voltage) or other effects, a transient simulation is mandatory. Under the assumption of no disturbances and a pure stationary electric field, the transient computation the field stress is often not necessary, because of the long operation time of power cables [1,11,25,26]. Thus, for the reliability of the insulation material it is sufficient to analyze the stationary electric field that is above the black dotted line in Figure 1a (∂/∂t ≈ 0). Excluding the conductivity model (7), the stationary electric field is not given in a closed analytic form. In addition to the explicit time integration of (8)- (10) or the explicit/implicit time integration of (11), the stationary electric potential is also obtained by a pseudo time stepping of the stationary current: The corresponding charge distribution is obtained from Poisson's Equation (9).
Using an explicit Euler method, the time step length Δt has to satisfy the CFL-criterion. With an implicit Euler method, there is no limit for Δt and the stationary electric field is theoretically

Fast Calculation of the Steady State Charge Distribution
In the case of overvoltage (e.g., impulse voltage) or other effects, a transient simulation is mandatory. Under the assumption of no disturbances and a pure stationary electric field, the transient computation the field stress is often not necessary, because of the long operation time of power cables [1,11,25,26]. Thus, for the reliability of the insulation material it is sufficient to analyze the stationary electric field that is above the black dotted line in Figure 1a (∂/∂t ≈ 0). Excluding the conductivity model (7), the stationary electric field is not given in a closed analytic form. In addition to the explicit time integration of (8)- (10) or the explicit/implicit time integration of (11), the stationary electric potential is also obtained by a pseudo time stepping of the stationary current: The corresponding charge distribution is obtained from Poisson's Equation (9).
Using an explicit Euler method, the time step length ∆t has to satisfy the CFL-criterion. With an implicit Euler method, there is no limit for ∆t and the stationary electric field is theoretically obtained after the first iteration. Resulting from the rather low number of degrees of freedom, the explicit Euler method allows a fast computation of many time steps, while the benefit of the implicit method is reduced, because of solving a non-linear system. By using different acceleration techniques, the computation time, using an explicit Euler method, is further reduced, [60,61].
During the pseudo time stepping of (36), oscillations may occur, resulting from high electric field values and the non-linear electric conductivity. A diverging electric field is seen in Figure 9a, resulting from these oscillations [71]. A planar MI insulation with thickness D = 20 mm is considered. The electric conductivity is described by (6) and the conductivity constants are seen in Table 1. The conductor temperature T i = 50 • C and the sheath temperature T a = 35 • C. Using an applied voltage of 470 kV, oscillations at the sheath start to increase after 5 iterations. High electric field values result in increased conductivity values, having an effect on the computation of the electric potential ϕ in (36). Decreasing the voltage from 470 kV to 400 kV yields lower and converging field values ( Figure 9b). Further possibilities to reduce the charge density and thus, the field strength are a lower temperature gradient and lower constants α and β [71]. Introducing an under-relaxation parameter ω, with ω  [0,1], the oscillations are minimized and the electric potential at the iteration step m + 1 is: Using ω < 1, more iterations are required to obtain a stationary solution, but the result does not diverge. As depicted in Figure 10b, with U = 470 kV and ω = 0.95, the electric field after the first iteration is slightly reduced to 44.5 kV/mm and convergence is obtained after 30 iterations [71]. Introducing an under-relaxation parameter ω, with ω ∈ [0, 1], the oscillations are minimized and the electric potential at the iteration step m + 1 is: Using ω < 1, more iterations are required to obtain a stationary solution, but the result does not diverge. As depicted in Figure 10b, with U = 470 kV and ω = 0.95, the electric field after the first iteration is slightly reduced to 44.5 kV/mm and convergence is obtained after 30 iterations [71].
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies and the electric potential at the iteration step m + 1 is: Using ω < 1, more iterations are required to obtain a stationary solution, but the result does not diverge. As depicted in Figure 10b, with U = 470 kV and ω = 0.95, the electric field after the first iteration is slightly reduced to 44.5 kV/mm and convergence is obtained after 30 iterations [71]. A comparison between the explicit Euler time integration and the fixed point iteration, to obtain the measured stationary charge distribution in [72,73], shows a faster computation using the fixed point iteration. Depending on the electric field and spatial resolution, the fixed point iteration yields a faster computation of about one order of magnitude. High electric fields result in fast charge movements that need to be considered with a small time step, using the explicit time integration. This results in more time steps and a higher computation time until a stationary solution is obtained. Fast charge dynamics are neglected using the fixed point iteration, which yields a lower computation time [71]. A comparison between the explicit Euler time integration and the fixed point iteration, to obtain the measured stationary charge distribution in [72,73], shows a faster computation using the fixed point iteration. Depending on the electric field and spatial resolution, the fixed point iteration yields a faster computation of about one order of magnitude. High electric fields result in fast charge movements that need to be considered with a small time step, using the explicit time integration. This results in more time steps and a higher computation time until a stationary solution is obtained. Fast charge dynamics are neglected using the fixed point iteration, which yields a lower computation time [71].

Electric Fields in Power Cables, Considering the Environment
HVDC power cables are typically buried within the earth or the sea bed. The thermal conductivity of the insulation material and the environment have an effect on the resulting temperature gradient within the insulation and thus, on the operation current and voltage. The calculations of the thermal rating of power cables go back to the late 19th and early 20th centuries, due to the advent of electricity. Simplified formulas are derived to compute the temperature of the cable sheath and in the vicinity of the cable, by using analogous descriptions to electrostatic fields [74][75][76]. In [77,78], temperature gradients within the cable insulation and the sheath are also considered to obtain the conductor and sheath temperature. Furthermore, a humidity-dependent electric conductivity of the soil is used.
To obtain the stationary temperature distribution, (16) is solved with ∂/∂t = 0. Typically the density δ and the specific heat capacity c p of the insulation material or the soil are constant, but the thermal conductivity λ of soil or air depends on the humidity and on the temperature. A coupling between the thermal and the electric field results from the thermal heat source κ T, As insulation losses have only a negligible influence on the resulting temperature distribution, the non-linear coupling can be reduced, resulting in a simplified temperature calculation [51,58,69].
A MI cable, buried 1.5 m, is considered and depicted in Figure 11 [77][78][79]. The electric conductivity is given by (6) and the necessary parameters are taken from Table 1, with r out = 52.4 mm and T i = 55 • C.
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies between the thermal and the electric field results from the thermal heat source   T ,| E | E  . As insulation losses have only a negligible influence on the resulting temperature distribution, the non-linear coupling can be reduced, resulting in a simplified temperature calculation [51,58,69]. A MI cable, buried 1.5 m, is considered and depicted in Figure 11 [77][78][79]. The electric conductivity is given by (6) and the necessary parameters are taken from Table 1, with rout = 52.4 mm and Ti = 55 °C.
where T is the absolute temperature in Kelvin [51,67]. The thermal conductivity of soil is more complex to describe, but a good assumption is a dependency on the humidity only [77,78]. For moderate humidity values, the thermal conductivity attains values between 0.47 W/(K•m) and 2.1 W(K•m) and increases with increasing humidity [51,67]. Depending on the soil thermal conductivity, the temperature gradient within the insulation and the soil itself varies. This results in different Temperature-independent thermal conductivities of MI and the outer PE sheath are assumed with λ MI = 0.167 W/(K·m) and λ PE = 0.3-0.4 W/(K·m) [30,67,80]. Within the range of −30 • C and +90 • C, the thermal conductivity of air λ Air depends on the temperature and is given by: where T is the absolute temperature in Kelvin [51,67]. The thermal conductivity of soil is more complex to describe, but a good assumption is a dependency on the humidity only [77,78]. For moderate humidity values, the thermal conductivity attains values between 0.47 W/(K·m) and 2.1 W(K·m) and increases with increasing humidity [51,67]. Depending on the soil thermal conductivity, the temperature gradient within the insulation and the soil itself varies. This results in different temperatures at the earth-air interface (x = 5.3 m, y = 2r out + 1.5 m). For example, with a soil thermal conductivity of 0.47 W/(K·m), the temperature at the earth-air interface is T = 32.4 • C and with 2.1 W/(K·m), the temperature is reduced to T = 26.7 • C. Within the insulation, an increased thermal conductivity results in an increased temperature gradient and, thus, in an increased field stress at the sheath ( Figure 12). Due to the large distance of 1.5 m between the cable and the earth-air interface, the temperature distribution within the insulation is approximately radial symmetric. The temperature within the insulation and the outer sheath is depicted in Figure 12a. The vertical black dotted line shows the interface between the insulation and the outer sheath. The corresponding electric field is seen in Figure 12b.
The temperature gradient within the insulation increases from 8 • C (0.47 W/(K·m)) to 18.3 • C (2.1 W/(K·m)), resulting in a decreasing electric field at the conductor of 27.5% and an increasing field at the sheath of 24.2%. With a varying environment temperature (T ∞ ), the effect on the insulation temperature also varies. With increasing T ∞ , the temperature gradient and the field stress at the outer sheath decreases. By contrast, the soil around the cable heats up and the temperature at the earth-air interface assumes its maximum value for dry soil (λ = 0.47 W/(K·m)). With decreasing T ∞ , the effect on the environment is reduced, but the insulation is facing higher stress levels [51]. Considering a buried cable pair, a higher influence on the environment is seen, due to the additional heat losses of the second cable (see Figure 13). Within the insulation, an increased thermal conductivity results in an increased temperature gradient and, thus, in an increased field stress at the sheath ( Figure 12). Due to the large distance of 1.5 m between the cable and the earth-air interface, the temperature distribution within the insulation is approximately radial symmetric. The temperature within the insulation and the outer sheath is depicted in Figure 12a. The vertical black dotted line shows the interface between the insulation and the outer sheath. The corresponding electric field is seen in Figure 12b. The temperature gradient within the insulation increases from 8 °C (0.47 W/(K•m)) to 18.3 °C (2.1 W/(K•m)), resulting in a decreasing electric field at the conductor of 27.5% and an increasing field at the sheath of 24.2%. With a varying environment temperature (T∞), the effect on the insulation temperature also varies. With increasing T∞, the temperature gradient and the field stress at the outer sheath decreases. By contrast, the soil around the cable heats up and the temperature at the earth-air interface assumes its maximum value for dry soil (λ = 0.47 W/(K•m)). With decreasing T∞, the effect on the environment is reduced, but the insulation is facing higher stress levels [51]. Considering a buried cable pair, a higher influence on the environment is seen, due to the additional heat losses of the second cable (see Figure 13). The temperature gradient within the insulation increases from 8 °C (0.47 W/(K•m)) to 18.3 °C (2.1 W/(K•m)), resulting in a decreasing electric field at the conductor of 27.5% and an increasing field at the sheath of 24.2%. With a varying environment temperature (T∞), the effect on the insulation temperature also varies. With increasing T∞, the temperature gradient and the field stress at the outer sheath decreases. By contrast, the soil around the cable heats up and the temperature at the earth-air interface assumes its maximum value for dry soil (λ = 0.47 W/(K•m)). With decreasing T∞, the effect on the environment is reduced, but the insulation is facing higher stress levels [51]. Considering a buried cable pair, a higher influence on the environment is seen, due to the additional heat losses of the second cable (see Figure 13). The temperature of the earth-air interface increases of about 4 • C, compared to a single cable. Furthermore, the temperature distribution within the insulation is no longer radially symmetric, as seen in Figure 13b. Between both conductors ("Area 1"), the temperature is higher compared to "Area 2", but the temperature gradient is higher in "Area 2". Utilizing a metallic sheath between the insulation and the outer sheath, the temperature outside of the cable shows negligible changes, but results in an approximately radially symmetric temperature distribution within the insulation. Due to the high thermal conductivity of the metallic sheath in comparison to the insulation or the outer sheath, the non-symmetry of the temperature, seen in Figure 13b, is reduced. The temperature has approximately equal values at r = r a . The metallic sheath has a negligible influence on the environment and the temperature at the earth-air-interface remains unaffected. Thus, one-dimensional electric field simulations are also applicable for cable pairs, if a metallic sheath is considered.
The simulation results in Section 5 do not include the effects of interfaces, where charges accumulate due to a high trap density. These effects are modeled with an additional spatial variation, presented in the following section.

Simulation of Space Charge Effects at Interfaces and Surfaces and of Moving Charge Packets
The blocking behavior of interfaces is not considered in the commonly used conductivity models (3)- (7). The modeling of charges at interfaces is given in Section 6.1, where analytic and simulation results of the electric field in a cable and a dual-dielectric interface are shown. Using space charge measurements, an empirical conductivity equation, considering heterocharge accumulation in polymeric cable insulation, is developed in Section 6.2. In Section 6.3, the transient simulation of charge packets and a comparison to measurements are seen. In [81,82], heterocharge build up at the electrodes is described with an extended conductivity equation, that has additional spatial variations towards the electrodes. A closed analytic description of the stationary electric field due to homo-or heterocharges at the electrodes within a cable insulation is shown in [70]. In [81], the additional spatial variations towards the electrodes are described with: where f con (r) yields spatial variations at the conductor ("con") and f sh (r) at the sheath ("sh"). The conductivity increases or decreases by a factor of n towards the electrodes and ζ is a constant that defines the gradient of f (r). With (7), the total conductivity is: Equal to (27), the electric field is the solution of: describing the effects at both electrodes by the additional term f (r). Using the assumption f con (r = r a ) = 1 and f sh (r = r i ) = 1, results in f (r) = f con (r) at the conductor and f (r) = f sh (r) at the sheath. With the temperature distribution (28), the solution of (41) is: Equal to (31), the constant K is computed with the applied voltage U and given by: where γ con = a −1/(1+v) (a−1)(ζ con +r i ) [a·(a 1/(1+v) −1)]ζ con (1+v) , A con = a −1/(1+v) −1 (ζ con +r i ) γcon , a = 1 + (n con − 1) exp(−1) , Energies 2020, 13, 5189

of 42
A comparison of (42) and (43) with (30) and (31) shows the additional term f (r) that is influenced by the electric field dependency only [70]. For example, the stationary electric field and charge distribution for n con = n sh = 1.5 (n con = n sh = 0.75), T i = 50 • C, T a = 35 • C, and ζ con = ζ sh = (r a − r i )/10 are depicted in Figure 14 [11,22]. The conductivity of MI is described by (7) and the necessary constants are given in Table 1. Resulting from the negative exponent of f (r) in (42), n > 1 corresponds to a decreasing electric field and, thus, in additional accumulated homocharges and n < 1 describes additional heterocharges. With n = 1 additional homo-or heterocharge distributions vanish and bulk effects are only considered (see Figure 14b).
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies con con con 1/ (1 ) con i con con ,, A comparison of (42) and (43) with (30) and (31) shows the additional term f(r) that is influenced by the electric field dependency only [70]. For example, the stationary electric field and charge distribution for ncon = nsh = 1.5 (ncon = nsh = 0.75), Ti = 50 °C, Ta = 35 °C, and ζcon = ζsh = (ra − ri)/10 are depicted in Figure 14 [11,22]. The conductivity of MI is described by (7) and the necessary constants are given in Table 1. Resulting from the negative exponent of f(r) in (42), n > 1 corresponds to a decreasing electric field and, thus, in additional accumulated homocharges and n < 1 describes additional heterocharges. With n = 1 additional homo-or heterocharge distributions vanish and bulk effects are only considered (see Figure 14b). A heterocharge distribution affects the electric field at the sheath, while a homocharge distribution needs to be considered for the electric field at the conductor. Using a classical "line-commutated current source converter" (CSC) HVDC system, the direction of the current remains constant and the polarity of the applied voltage changes, to invert the direction of the power flow. This polarity reversal results in high electric field stress at the conductor, where the negative electric fields due to the space charges and the applied voltage are superimposed [2,4,83]. To avoid high field values, XLPE cables are used in modern "self-commutated voltage source converters" (VSC) HVDC systems [2,83]. A heterocharge distribution affects the electric field at the sheath, while a homocharge distribution needs to be considered for the electric field at the conductor. Using a classical "line-commutated current source converter" (CSC) HVDC system, the direction of the current remains constant and the polarity of the applied voltage changes, to invert the direction of the power flow. This polarity reversal results in high electric field stress at the conductor, where the negative electric fields due to the space charges and the applied voltage are superimposed [2,4,83]. To avoid high field values, XLPE cables are used in modern "self-commutated voltage source converters" (VSC) HVDC systems [2,83].

Modeling of Charges Close to Interfaces of Different Dielectrics
Interfaces between two different dielectrics are found in cable joints and cable terminations. Accumulated interfacial charges are described by the Maxwell-Wagner-Sillars polarization, if constant permittivity and conductivity values are used within both dielectrics [14,[84][85][86]. Due to the electric field depending on conductivity, the Maxwell-Wagner-Sillars model is not applicable [14,85]. Furthermore, space charges accumulate in the vicinity of interfaces, resulting from physical and chemical disorders and thus, a high trap density [14].
Space charge measurements within two dielectrics of XLPE (ε r,XLPE = 2.3) and ethylene propylene rubber (EPR, ε r,EPR = 2.9) are found in [84]. The applied voltage is U = 30 kV, the conductor radius is r i = 1.8 mm, with a temperature of T i = 64 • C and the sheath radius is r a = 3.9 mm, with T a = 42 • C. The interface radius is r Int = 3.25 mm. For the stationary temperature distribution, the thermal conductivity is λ XLPE = 0.27 W/(K·m) and λ EPR = 0.3 W/(K·m). The space charge distribution is obtained after t = 20,000 s [84,87]. For the electric conductivity, (7) is used, with the constants given in [14]. A sketch of the dual-dielectric-interface is seen in Figure 15a.
with EEPR(rInt) being the electric field within the EPR material at the interface (r = rInt) and EXLPE(rInt) is the electric field in the XLPE material at the interface [84,86]. Measurements and simulation results with and without accumulated charges close to the interface are seen in Figure 15b, where the corresponding electric field is depicted in Figure 15c. The simulation results are obtained with (39) and (45)  The simulation results, considering charges close to the interface show fewer differences to the measurements, compared to the conductivity model (7). Furthermore, the accumulated charges in the vicinity of the interface (rInt = 3.25 mm) increase the electric field of 20%, highlighting the importance of accurate modelling of HVDC components with multiple subcomponents. To consider accumulated charges in the vicinity of interfaces, an additional conductivity gradient is applied [70]. Equal to the modelling of homo-and heterocharge in cables (39), accumulated charges close to the interface are described with: where n and ζ have the equal meaning, compared to (39), but f XLPE is only valid for r i ≤ r ≤ r Int and f EPR is only valid for r Int ≤ r ≤ r a . The interfacial charges δ Int are computed with: with E EPR (r Int ) being the electric field within the EPR material at the interface (r = r Int ) and E XLPE (r Int ) is the electric field in the XLPE material at the interface [84,86]. Measurements and simulation results with and without accumulated charges close to the interface are seen in Figure 15b, where the corresponding electric field is depicted in Figure 15c. The simulation results are obtained with (39) and (45) [70]: The simulation results, considering charges close to the interface show fewer differences to the measurements, compared to the conductivity model (7). Furthermore, the accumulated charges in the vicinity of the interface (r Int = 3.25 mm) increase the electric field of 20%, highlighting the importance of accurate modelling of HVDC components with multiple subcomponents.

Empirical Conductivity Equation for the Simulation of Heterocharges in Polymeric Cable Insulations
A different approach to model the electric conductivity with additional spatial variations in the vicinity of the electrodes is described in [88][89][90]. The conductivity is given by: where K con and K sh yield conductivity variations at the conductor and the sheath. The constant r x defines the distance between the conductor (r i ) and the position of the highest gradient (K con − K sh = 0.5 at r = r x ) and χ defines the gradient of K con and K sh in the vicinity of both electrodes, which has an effect on the magnitude and the shape of the resulting heterocharge distribution [88,89]. Considering a planar insulation, r i = 0 and r a = D. To define the constants r x and χ, different measurements are taken from literature. The empirical conductivity Equation (47) is valid for heterocharges only, because of considered low applied electric fields (<20 kV/mm) only [88,89]. Furthermore, measured heterocharges in XLPE and LDPE insulations, at an applied average electric field below 20 kV/mm, are found in [73,[91][92][93]. For simplicity, χ con = χ sh = χ is assumed. To obtain the constants r x and χ, a "conductivity gradient region" ∆ is defined and yields r x = ∆/2 and 10χ ≈ ∆ [88,89].
A common measurement technique is the pulsed electro acoustic (PEA) method, where surface charges and bulk space charges are subjected to a filtering process during the measurement, giving inaccurate results at the conductor and the sheath [15]. Thus, a Gaussian curve instead of a delta function is obtained, if a space charge free material under voltage load is measured [40,[93][94][95]. An example of a space charge measurement within an XLPE cable insulation is seen in Figure 16, also indicating heterocharges and "filtered" surface charges at both electrodes [73]. where Kcon and Ksh yield conductivity variations at the conductor and the sheath. The constant rx defines the distance between the conductor (ri) and the position of the highest gradient (Kcon − Ksh = 0.5 at r = rx) and χ defines the gradient of Kcon and Ksh in the vicinity of both electrodes, which has an effect on the magnitude and the shape of the resulting heterocharge distribution [88,89]. Considering a planar insulation, ri = 0 and ra = D. To define the constants rx and χ, different measurements are taken from literature.
A common measurement technique is the pulsed electro acoustic (PEA) method, where surface charges and bulk space charges are subjected to a filtering process during the measurement, giving inaccurate results at the conductor and the sheath [15]. Thus, a Gaussian curve instead of a delta function is obtained, if a space charge free material under voltage load is measured [40,[93][94][95]. An example of a space charge measurement within an XLPE cable insulation is seen in Figure 16, also indicating heterocharges and "filtered" surface charges at both electrodes [73]. For an accurate comparison between simulated space charges and a measured charge distribution, surface charges at both electrodes need to be considered. Surface charges at the conductor δ+ and the sheath δ-are derived from [40] and approximately described by: Figure 16. Example of a space charge measurement within an XLPE cable insulation. "Filtered" surface charges appear as a Gaussian curve [73].
For an accurate comparison between simulated space charges and a measured charge distribution, surface charges at both electrodes need to be considered. Surface charges at the conductor δ + and the sheath δ − are derived from [40] and approximately described by: for a cylindrical geometry. Induced charges, resulting from accumulated space charges, represent the first term on the right-hand side in (48). The second term represents the capacitance charges induced by the polarization of the material under an applied voltage. To combine the surface charges and the bulk space charge, δ + and δ − are divided by the spatial discretization ∆h [96]. Finally, to compare simulations and measurements, the simulation results need to be converted to an equivalent signal, using a Gaussian filter [40].
For the calibration of the filter, the space charge density vanishes within the bulk. The voltage is applied for a short period of time, to compare the space charge measurements and simulation results at t ≈ 0 s, for example seen in Figure 17 [95,97]. The black dashed line are the measurements, the black line is the unfiltered simulation result and the red line is the filtered simulation result. With no bulk space charges, (48) is reduced to δ + = ε 0 ε r E(r i ) and δ − = ε 0 ε r E(r a ) [97]. In Figure 17, the applied voltage is U = 90 kV, r i = 4.5 mm, r a = 9 mm and ε r = 2.3 [14]. For the calibration of the filter, the space charge density vanishes within the bulk. The voltage is applied for a short period of time, to compare the space charge measurements and simulation results at t ≈ 0 s, for example seen in Figure 17 [95,97]. The black dashed line are the measurements, the black line is the unfiltered simulation result and the red line is the filtered simulation result. With no bulk space charges, (48) is reduced to δ+ = ε0εrE(ri) and δ− = ε0εrE(ra) [97]. In Figure 17, the applied voltage is U = 90 kV, ri = 4.5 mm, ra = 9 mm and εr = 2.3 [14]. Several space charge measurements within a planar and a cylindrical insulation material of XLPE and LDPE are seen in [72,73,91,93]. By minimizing the difference between the simulated filtered space charge distribution and the measured charge distribution, the constants for rx = Δ/2, 10χ ≈ Δ and Δ are summarized in Table 2. The measurements in [91] and the corresponding simulations are for example seen in Figure 18. To describe the electric conductivity of XLPE and LDPE, (4) is used, where the constants are given in [14,56,89,98,99].
The constants χ and rx in (47) are approximated with χ = (0.25•D)/10 and rx = (0.25•D)/2 for a planar insulation and with χ = (0.25•(ra − ri))/10 and rx = (0.25•(ra − ri))/2 for a cylindrical insulation [88,89]. Utilizing a constant applied voltage with positive and negative polarity, the polarity of the charge density also changes, but the magnitude and the shape remains unaffected. This is called the "mirror image effect" and seen in Figure 18a-f. The effect is discussed in [97], where possible explanations are a spatially varying polarization of the dielectric and the injection of charges at the electrodes [100,101]. Several space charge measurements within a planar and a cylindrical insulation material of XLPE and LDPE are seen in [72,73,91,93]. By minimizing the difference between the simulated filtered space charge distribution and the measured charge distribution, the constants for r x = ∆/2, 10χ ≈ ∆ and ∆ are summarized in Table 2. The measurements in [91] and the corresponding simulations are for example seen in Figure 18. To describe the electric conductivity of XLPE and LDPE, (4) is used, where the constants are given in [14,56,89,98,99]. Table 2. Distance constant χ and position r x for the simulation of space charge measurements in [72,73,91,93]. The absolute value of the applied voltage is |U| = 20 kV.

Ref.
χ The constants χ and r x in (47) are approximated with χ = (0.25·D)/10 and r x = (0.25·D)/2 for a planar insulation and with χ = (0.25·(r a − r i ))/10 and r x = (0.25·(r a − r i ))/2 for a cylindrical insulation [88,89]. Utilizing a constant applied voltage with positive and negative polarity, the polarity of the charge density also changes, but the magnitude and the shape remains unaffected. This is called the "mirror image effect" and seen in Figure 18a-f. The effect is discussed in [97], where possible explanations are a spatially varying polarization of the dielectric and the injection of charges at the electrodes [100,101].
Energies 2020, 13, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies Accumulated heterocharges are approximately seen in one quarter of the insulation thickness, independent of the insulation material or the thickness. One possible assumption for accumulated heterocharges within this region might be the resolution of the measurement technique. On the other hand, the PEA method has a resolution of 1.6 μm for a one dimensional planar insulation with a thickness of 25-27,000 μm and 0.1-1 mm for a cable insulation with a thickness of 3.5-20 mm [94]. Thus, the resolution is accurate enough to separate between heterocharges and bulk space charges. Differences between the simulations and the measurements in Figure 18 may come from the constants to describe (4) or the approximate description of the surface charges in (48). The electric conductivity of polymeric insulation materials differs with the preparation, which affects the constants of the conductivity models (3)- (7). Conductivity measurements of XLPE, taken from several references, are depicted e.g., in [95]. The measurements differ of about one order of magnitude and consequently different constants are obtained, when using the same conductivity model (4). The distribution of a measured space charge density shows a dependency on many factors, e.g., the electric conductivity, the local electric field or the electrode material. Consequently, it is very difficult to simulate the charge distribution of different references with a unique conductivity model, even if it is the same material [95,102,103]. Differences between a simulated and Accumulated heterocharges are approximately seen in one quarter of the insulation thickness, independent of the insulation material or the thickness. One possible assumption for accumulated heterocharges within this region might be the resolution of the measurement technique. On the other hand, the PEA method has a resolution of 1.6 µm for a one dimensional planar insulation with a Energies 2020, 13, 5189 24 of 42 thickness of 25-27,000 µm and 0.1-1 mm for a cable insulation with a thickness of 3.5-20 mm [94]. Thus, the resolution is accurate enough to separate between heterocharges and bulk space charges.
Differences between the simulations and the measurements in Figure 18 may come from the constants to describe (4) or the approximate description of the surface charges in (48). The electric conductivity of polymeric insulation materials differs with the preparation, which affects the constants of the conductivity models (3)- (7). Conductivity measurements of XLPE, taken from several references, are depicted e.g., in [95]. The measurements differ of about one order of magnitude and consequently different constants are obtained, when using the same conductivity model (4). The distribution of a measured space charge density shows a dependency on many factors, e.g., the electric conductivity, the local electric field or the electrode material. Consequently, it is very difficult to simulate the charge distribution of different references with a unique conductivity model, even if it is the same material [95,102,103]. Differences between a simulated and a measured space charge distribution are minimized, if conductivity and space charge measurements of the material are provided. Then, an error that might occur due to the use of conductivity measurements of different references, is reduced.

Transient Simulation of Charge Packets within a Cross-Linked Polyethylene (XLPE) Cable Insulation
Space charge and conductivity measurements of a medium voltage (MV) XLPE cable are given in [14]. These measurements show space charges, moving from the conductor to the sheath and a buildup of heterocharges (see Figure 19l). In addition to a varying conductivity, space charges accumulate due to a varying permittivity as well (see (24)). Analogously to [81] the permittivity is described by: where n con,ε and n sh,ε are the increasing factors at the conductor ("con") and the sheath ("sh") and ζ con,ε and ζ sh,ε are constants defining the gradient of ε r (r) at both electrodes. The constant bulk permittivity is ε r,Bulk . A reasonable modelling approach for moving charges is an electric conductivity, having the shape of a Gaussian pulse and moving through the insulation [45,95]. With the non-linear bulk electric conductivity κ T, → E the total electric conductivity is: where n con,κ and n sh,κ are the increasing factors, ζ con,κ and ζ sh,κ constants for the shape of the pulse and v κ is the velocity of the Gauss pulse (velocity of the moving charge pulse). The constants for the simulation results are summarized in Table 3. The cylindrical geometry has a conductor radius of r i = 4.5 mm, a sheath radius of r a = 9 mm and a bulk permittivity of ε r,Bulk = 2.25. A time-independent temperature distribution is used and described by (28). For the bulk electric conductivity, the hopping model (4) is used, with the constants E A,2 = 1.48 eV, K 2 = 1·10 14 A/m 2 and γ 2 = 2 × 10 −7 m/V [14,81]. A comparison between the measurements and the simulation results, at the time t = 0 s, t = 10.000 s and t = 20.000 s, is seen in Figure 19. The results are obtained along the red evaluation line in Figure 1b. The used constants in Table 3 yield the best fit between the measurements in [14] and the simulation results. Increasing factors n sh,ε = n sh,κ = 1 indicate no increase/decrease of the permittivity and conductivity, which makes the corresponding constants ζ sh,ε and ζ sh,κ unnecessary. The conductivity factor n con,κ increases with voltage and temperature, where possible reasons are detrapped charges and thermally activated and electric field assisted injection [20].  (48), together with the filtering process, yields the highest differences at both electrodes, while within the bulk, the differences remain low (see Figures  18 and 19).    (49) and (50), for each applied voltage and temperature gradient [14,95].
moving regions, III. a shape that is approximately maintained during motion, IV. a periodic process (repetitive injection).
The simulation results in Figure 19 show space charge regions moving through the insulation that are especially seen in Figure 19k,l (I), (II). During the motion, the charge pulse slightly changes shape (III), but a periodic repetition (IV) is not seen. In [42], a repetitive injection of charge packets is seen at electric fields >40 kV/mm, whereas the electric fields in the simulations are between 5 kV/mm and 30 kV/mm.
A comparison between the simulated charge density and the measurements is depicted in Figure 19 and shows a good agreement with the ocular norm. For the reliability of the cable insulation, the electric field needs to be determined. Thus, the mean difference between the simulated and the measured electric field is shown in Table 4. The mean difference lies between 1% and 22%. For example a mean difference between the filtered simulation results and the measurements of 12.5% (Figure 19l) is depicted in Figure 20. In comparison to the conductivity model in [14], the mean differences are in the same range in (a) and (b) and below the values of [14] in (c) and (d). With (49) and (50), one possibility of high values in Table 4 are the filtered surface charges at both electrodes. The approximation in (48), together with the filtering process, yields the highest differences at both electrodes, while within the bulk, the differences remain low (see Figures 18 and 19).
Compared to the unfiltered simulation results in Figure 20, the measurements show lower values at both electrodes, due to the filtering process of the measurement technique. The surface charges and the bulk space charges are "mixed" yielding a Gaussian-like shape at both electrodes (see Figure 16). During the simulation, these mixed charges result in a reduction to the electric field and can be misinterpreted as only being homocharges (see Figures 19 and 20). Compared to the maximum initial electric field of 28.85 kV/mm, the unfiltered electric field at t = 20,000 s is about 32 kV/mm (increase of 11%). The electric field at the sheath is further increasing due to the charge pulse moving closer to the sheath, resulting in the accumulation of heterocharges (see Figure 14) [42]. A comparison between the different conductivity models in Sections 5 and 6 is depicted in Table 5. Table 4. Mean difference between the simulated and the measured electric field, using (49) and (50) and the results in [14,95]. Equations (49) and (50) 5.3% 6.1% 10.2% [14] 5.3% 6.8% 8.6% Equations (49) and (50) Table 4. Mean difference between the simulated and the measured electric field, using (49) and (50) and the results in [14,95]. t = 0 s t = 10.000 s t = 20.000 s (a) Equations (49) and (50) 5.3% 6.1% 10.2% [14] 5.3% 6.8% 8.6%
Compared to the unfiltered simulation results in Figure 20, the measurements show lower values at both electrodes, due to the filtering process of the measurement technique. The surface charges and the bulk space charges are "mixed" yielding a Gaussian-like shape at both electrodes (see Figure 16). During the simulation, these mixed charges result in a reduction to the electric field and can be misinterpreted as only being homocharges (see Figures 19 and 20). Compared to the maximum initial electric field of 28.85 kV/mm, the unfiltered electric field at t = 20,000 s is about 32 kV/mm (increase of 11%). The electric field at the sheath is further increasing due to the charge pulse moving closer to the sheath, resulting in the accumulation of heterocharges (see Figure 14) [42]. A comparison between the different conductivity models in Sections 5 and 6 is depicted in Table 5.
Commonly used in literature are the conductivity models, represented in the Equations (3)-(7). Space charge accumulation and electric field variations within the insulation, resulting from temperature and electric field changes, are described. The stationary electric field is given in a closed analytic form. The electric field is connected to the applied temperature gradient, but only a mean charge density of one sign is computed. Effects at interfaces and surfaces are neglected. Such effects are described with additional spatial variations, e.g., (39) or (45) and (40), but the necessary constants are difficult to determine. Using (47), the additional constants to describe effects in the vicinity of the electrodes are evaluated with space charge measurements. By contrast, (47) is limited to Figure 20. Corresponding electric field of the filtered and unfiltered simulation results in Figure 19l, compared with measurements at the moment t ≈ 20,000 s. [95].
Commonly used in literature are the conductivity models, represented in the Equations (3)-(7). Space charge accumulation and electric field variations within the insulation, resulting from temperature and electric field changes, are described. The stationary electric field is given in a closed analytic form. The electric field is connected to the applied temperature gradient, but only a mean charge density of one sign is computed. Effects at interfaces and surfaces are neglected. Such effects are described with additional spatial variations, e.g., (39) or (45) and (40), but the necessary constants are difficult to determine. Using (47), the additional constants to describe effects in the vicinity of the electrodes are evaluated with space charge measurements. By contrast, (47) is limited to heterocharges only and, equal to (39) or (45) and (40), only valid for the stationary case. Transient processes are simulated with (49) and (50), where additional spatial variations move within the insulation. The models (49) and (50) reduce the difference between simulation and measurement, but the additional constants vary over time and are determined by measurements. Table 5. Comparison of different conductivity models and their description of space charge accumulation within the insulation.

Model Description of Charge Dynamics Limitations
Equations (3)- (7) Charge accumulation due to a temperature gradient.
Computation of an average charge density of one sign only. No effects at interfaces and surfaces.
Equations (39) or (45) and (40) Charge distribution within the insulation and at interfaces and surfaces.
Limited to the stationary case and the additional constants are difficult to determine Equation (47) Description of the stationary bulk and heterocharge distribution in polymeric insulations.
Limited to the stationary case and heterocharges only.
Equations (49) and (50) Simulation of transient processes, homo and heterocharges at interfaces and charges within the bulk.
Many additional constants are used that change over time and are determined by space charge measurements.
The insulation materials in Sections 5 and 6 are considered without defects. It is often difficult to determine the defects and their electrical and thermal material characteristics. Commonly, an air-filled void with a vanishing electric conductivity (κ Void = 0), in the absence of partial discharges (PD), is assumed [104][105][106].
Different computation approaches are presented in Figures 4 and 5, to obtain the results in Sections 5 and 6. These approaches result in different accuracies of the stationary charge and electric field distribution that are compared in the following section.

Accuracy of the Electric Field Computation within Cables and Cable Joints
Compared to the scalar potential field formulation ( Figure 5), the space charge-oriented field formulation in Figure 4 uses an additional averaging process to compute the charge density. Utilizing the cable geometry in Table 1 and the conductivity model (7) to describe a MI insulation, the analytic solution of the stationary electric field is given by (30), where T i = 50 • C and T a = 35 • C [11,59]. The analytic solution, together with the simulation results, using Figures 4 and 5 are depicted in Figure 21a. Both formulations show a good accuracy in the ocular norm. Utilizing the space charge-oriented field formulation (8)-(10), inaccurate results are especially seen at the conductor, while using (11) the error has its maximum value at the sheath (see Figure 21b). With (8)-(10), the maximum error is 0.723%, with an average error of 0.3%. With (11), the maximum error (0.677%) and the average error (0.25%) are slightly lower. Equal to the results in Section 5.2, high electric field values are obtained during the time integration, resulting in inaccurate field values. Increasing the stationary electric field by decreasing the field dependency v (see (30)), the maximum error at the conductor is >1.5% using (8)-(10), but utilizing (11), the maximum error at the sheath has only slightly increased (0.74%) [59].
Cable joints and terminations are geometrically more complex than cable insulations, due to different dielectric materials and the resulting interfaces. Resulting from these interfaces and other factors it is widely believed that cable joints and terminations are the least reliable components of a HVDC cable system [4,84]. Furthermore, the use of field-grading materials (FGM), results in difficulties while simulating such components. Such materials have a strongly non-linear dependency on the electric field and a lower dependency on the temperature [107,108]. In [107], the electric conductivity is given by:   In Figure 22b,c, the initial electric field (t = 0 s) and the stationary electric field are seen. Especially in the vicinity of the triple point, the field-grading material reduces the electric field stress. The stationary field distribution in Figure 22c is computed using the scalar potential field formulating (11), while the space charge formulation (8)-(10) yields unstable results at t m=13 = 13Δt, as seen in Figure 22d,e. The geometry of the simulated cable joint is depicted in Figure 22a [59,109]. The simulation time t END = 2 µs, with a time step of ∆t = 0.0133 µs, yields a relative difference of Φ m+1 − Φ m / Φ m · 100 % < 1 % between the discrete time steps m and m + 1. For simplicity and due to the short time t END , the electric conductivity of XLPE and insulating silicone rubber (LSR) are assumed to be constant and given by κ XLPE = 10 −15 S/m and κ LSR = 5 × 10 −13 S/m. The relative permittivity of both materials is ε r,XLPE = 2.3 and ε r,LSR = 3.5 [59,109]. to Δt = 0.0066 μs. With (11), an averaging process is needed to compute the electric field E = -grad(φ). With (8)-(10), an additional averaging process is needed to compute the charge density (see (8)). Due to the additional averaging process and the non-linear electric conductivity, unstable electric field values are delivered [59]. In Figure 22b,c, the initial electric field (t = 0 s) and the stationary electric field are seen. Especially in the vicinity of the triple point, the field-grading material reduces the electric field stress. The stationary field distribution in Figure 22c is computed using the scalar potential field formulating (11), while the space charge formulation (8)- (10) yields unstable results at t m=13 = 13∆t, as seen in Figure 22d,e.
High electric field values occur in the vicinity of the triple point. The sharp edge in the cable joint geometry results in a compression of the potential lines and thus, in a high potential gradient. Using (8)-(10), electric field values >5 kV/mm are simulated at t m=12 = 12∆t and the electric conductivity has values of 10 −2 S/m (Figure 22e). Consequently, the chosen time step size yields unstable field values. Stable field values are computed with (8)-(10), if the time step size is reduced to ∆t = 0.0066 µs. With (11), an averaging process is needed to compute the electric field → E = −grad(ϕ). With (8)- (10), an additional averaging process is needed to compute the charge density (see (8)). Due to the additional averaging process and the non-linear electric conductivity, unstable electric field values are delivered [59].

The Electric Field of Ground Electrodes, Considering the Effect of Electro-Osmosis
A special case of an electro-quasistatic field is seen in the vicinity of HVDC ground electrodes, where the time varying electric field is the result of a drying process of the soil and its humidity dependency. A short overview of the computation is given in this section.
Common system configurations of HVDC transmission are monopole or bipole [83,110,111]. To use the earth or the sea for the return current, ground electrodes within the earth or on the sea bed are installed. Such electrodes are placed in different depths to ensure high conductivity values and are surrounded by a coke bed to reduce electrolytic corrosion and thus, the loss of electrode material [110,112].
Under the influence of the resulting electric field of ground electrodes buried within the earth, the water around the anode drift away and the soil around the electrode dries out due to the effect of electro-osmosis [110,112]. Electro-osmosis needs to be avoided due to a decrease of the electric conductivity and an increase of the electric field and the surface step voltage. To prevent soil from drying, some HVDC projects have arrangements to add water, if the soil moisture concentration falls below a certain value [110].
Resulting from the electric conductivity of earth that depends on the humidity and the temperature distribution, the electric field varies both in space and time, respectively. This special case of an electro-quasistatic field, shows no accumulated space charge and results from an electric conductivity that depends on a varying humidity concentration.
Water-saturated soil behaves as an electrolyte, where charges are injected into the soil if electrodes are connected. If the charges move under an applied electric field, they carry water molecules by exerting a viscous drag on the water around them. Due to more mobile positive charges than negative ones, there is a net water flow from the anode to the cathode [113][114][115]. The continuity equation of water flow describes the time-varying humidity concentration (volumetric water content) θ with: where → q W is the water flow [116,117]. Darcy's law describes the water flow under the influence of an applied electric field with: where k is the electro-osmotic hydraulic conductivity (also called "electro-osmotic conductivity" or "electro-osmotic permeability"). Due to the constant current of an HVDC system, the electric potential ϕ and the electric field → E = −grad ϕ are the solution of the stationary current problem (36). With losses within the conductor and the soil, the temperature is evaluated from (16).
The electric conductivity of soil is difficult to determine due to the dependency on several factors like the humidity, the temperature distribution, the composition and the size of silt, clay and sand particles. From [118,119], the electric conductivity is approximately correlated to the humidity and Energies 2020, 13, 5189 31 of 42 has a low temperature dependence. Using measurements from [120][121][122][123][124], the electric conductivity is described with the model: where the constants are κ 0,θ = 2.5 × 10 −1 S/m, b θ = 1.5 and κ s = 2 × 10 −5 S/m [117]. To consider the temperature dependency, (54) is extended by a correction factor g(T) and finally given by: with the absolute temperature in kelvin [119]. Measurements in [113,114,125] show the electro-osmotic conductivity with a humidity dependency only. Due to no water flow for a vanishing electro-osmotic conductivity, k is described by: The constants are k 0,θ = 8.5 × 10 −9 m 2 /(Vs) and a θ = 1.4. Heat losses within the conductor of the ground electrode and within the soil result in a temperature increase where the thermal conductivity decreases with decreasing humidity [89]. Additionally, the thermal conductivity shows a dependency on the temperature over a temperature range of 2 • C-92 • C in [126]. Measurements of the thermal conductivity λ in [77] and [126] are approximated by: with the temperature dependent variables α θ = −0.024·(T-273.15) + 3.049, β θ = 0.067·(T-273.15) + 13.867 and γ θ = 6.67 × 10 −4 ·(T-273.15) + 0.299. The measurements of the soil electric conductivity, electroosmotic conductivity and thermal conductivity, together with their corresponding models (55)-(57) are seen in Figure 23 [127]. The time-varying humidity θ, the electric scalar potential φ and consequently the electric field are computed within a coupled electro-thermal-electro-osmotic simulation of (16), (36), (52) and (53) [129]. Thus, humidity variations immediately affect the thermal conductivity and the temperature, resulting in a weak coupling between the quantities θ and T in Figure 24 [127]. Simulation results in [127] indicate electro-osmosis, with time constants up to years. Thus, it is only considered in the close vicinity of the electrode [129]. High electric field values also increase the surface step voltage right at the electrode. For short operation times of a few days, electro-osmosis might be omitted within the model, used for the computation of the surface step voltage, as its effects are negligible.  The electro-osmotic conductivity is in the range of 10 −9 m 2 /(Vs). The effect of thermo-osmosis, due to a temperature gradient, is equally described as (53). By contrast, the thermo-osmotic conductivity is about 3 orders of magnitude lower than k and can be neglected in HVDC applications [128].
The time-varying humidity θ, the electric scalar potential ϕ and consequently the electric field are computed within a coupled electro-thermal-electro-osmotic simulation of (16), (36), (52) and (53). The coupled problem can only be solved with an explicit time integration, where Figure 24 shows a possible pseudo code. Equal to Figures 4 and 5, the computation stops if a predefined time t = t END or θ m+1 − θ m / θ m < η, where η 1 is the stop threshold, are obtained. Compared to the thermal problem, the time constant of the transient electro-osmosis process is higher [129]. Thus, humidity variations immediately affect the thermal conductivity and the temperature, resulting in a weak coupling between the quantities θ and T in Figure 24 [127].

Conclusions and Outlook
This review showed the advances in the calculation and simulation of HVDC cables in the past few years. With an overview of the charge transport in commonly used insulation materials, a description of the conductivity models in use and the tools necessary to simulate HVDC cable components, with or without consideration of the environment, were presented. Finally, results of the charge and electric field behavior in HVDC cable insulations were shown and compared against analytic solutions and measurements.
Utilizing a conductivity-based cable model, the space charge-oriented field formulation and the scalar potential oriented field formulation were able to compute the time-varying electro-quasistatic field and space charge distribution within the insulation. Different models were found in literature to describe the non-linear electric conductivity and considered bulk effects only. These models were able to describe e.g., the "mirror image effect" or the effect of "field inversion". The conductivity and the permittivity models were extended to consider effects at interfaces and surfaces and to model moving charge packets. Utilizing the extended models, differences between simulations and measurements were reduced, compared to commonly used models. Furthermore, a stationary heterocharge distribution was found to accumulate in 25% of the insulation thickness, independent of the insulation material.
Fast charge movements were considered with a low time step length during the time integration, using an explicit Euler method. To reduce the computation time, pseudo time stepping of the stationary current was applied, where oscillations were reduced with under-relaxation techniques.
Compared to the scalar potential formulation, the space charge formulation showed a lower accuracy, when computing complex geometries or materials with a strongly non-linear electric conductivity. The space charge formulation yielded unstable field values, due to an additional averaging process, when computing the charge density. Simulation results in [127] indicate electro-osmosis, with time constants up to years. Thus, it is only considered in the close vicinity of the electrode [129]. High electric field values also increase the surface step voltage right at the electrode. For short operation times of a few days, electro-osmosis might be omitted within the model, used for the computation of the surface step voltage, as its effects are negligible.

Conclusions and Outlook
This review showed the advances in the calculation and simulation of HVDC cables in the past few years. With an overview of the charge transport in commonly used insulation materials, a description of the conductivity models in use and the tools necessary to simulate HVDC cable components, with or without consideration of the environment, were presented. Finally, results of the charge and electric field behavior in HVDC cable insulations were shown and compared against analytic solutions and measurements.
Utilizing a conductivity-based cable model, the space charge-oriented field formulation and the scalar potential oriented field formulation were able to compute the time-varying electro-quasistatic field and space charge distribution within the insulation. Different models were found in literature to describe the non-linear electric conductivity and considered bulk effects only. These models were able to describe e.g., the "mirror image effect" or the effect of "field inversion". The conductivity and the permittivity models were extended to consider effects at interfaces and surfaces and to model moving charge packets. Utilizing the extended models, differences between simulations and measurements were reduced, compared to commonly used models. Furthermore, a stationary heterocharge distribution was found to accumulate in 25% of the insulation thickness, independent of the insulation material.
Fast charge movements were considered with a low time step length during the time integration, using an explicit Euler method. To reduce the computation time, pseudo time stepping of the stationary current was applied, where oscillations were reduced with under-relaxation techniques.
Compared to the scalar potential formulation, the space charge formulation showed a lower accuracy, when computing complex geometries or materials with a strongly non-linear electric conductivity. The space charge formulation yielded unstable field values, due to an additional averaging process, when computing the charge density.
The electric field of HVDC ground electrodes is a special case of an electro-quasistatic field. Due to the effect of electro-osmosis, the electric conductivity and the electric field within soil varies in time, resulting from the constant current. Electro-osmosis is a slow process and needs to be considered only in the close vicinity of the electrode.
A description of charge dynamics in unmodified, and modified insulation materials are still part of ongoing research. Modified insulations consist of polymeric materials with added nano-dielectrics that show an influence on the temperature and electric field dependency and consequently on the charge distribution within the bulk and at interfaces. Another area of research includes liquid insulation materials, e.g., transformer oil, where the drift and the diffusion of charges need to be considered. Conductivity models are not applicable to describe both. The bipolar charge-transport model may describe drift and diffusion, but needs to be improved or simplified, making it usable for complex geometries.

Conflicts of Interest:
The authors declare no conflict of interest.