Fast calculation of steady‐state charge distribution in high voltage power cables

Modern high voltage direct current (HVDC) power cables consist of nonlinear electric field and temperature depending insulation materials, eg, polyethylene (PE) or mass‐impregnated paper (MI). Within the insulation material, slowly time varying electric fields are the result of accumulated space charges. Due to the long operation time of power cables, the stationary charge distribution and the resulting additional field stresses inside the insulation material need to be taken into account, to ensure a stable operation of the cable system. The long accumulation time of charges in comparison to fast charge dynamics usually requires extended computation times within transient electro‐quasistatic simulation models. An alternative approach to calculate the final steady‐state charge distribution arises from the simultaneous solution of both, the stationary current and the electrostatic potential problem, which are nonlinearly coupled by Ohm's law. Each of the two potential equations are discretized using the finite difference method. Utilizing a nonlinear fixed point iteration for the electric potential, the steady‐state solution of the charge distribution is computed until convergence is obtained. During this pseudo time stepping process, oscillations may occur that yield incorrect results of the charge distribution. These oscillations are eliminated utilizing underrelaxation strategies. Numerical tests on a cylindrical and a plane parallel problem, using both approaches for the calculation of the steady‐state charge distribution, show a faster computation with the nonlinear fixed point iteration, reducing the computation time by at least one order of magnitude.

Accumulated space charges result in a slowly time varying electric field, such that electric stresses under these DC conditions differ significantly from those of AC conditions. 3 Depending on the temperature in the insulation, it may take a few hours up to days until a stationary charge distribution is formed inside the insulating material after energizing the cable. 2 This is obtained by the time constant where κ is the electric conductivity, typically in the order of 10 −15 S/m, and ε = ε 0 ε r , with the dielectric constant ε 0 = 8.854 Á 10 −12 As/(Vm) and ε r is the relative permittivity. 4,5 The time to obtain a stationary charge and electric field distribution is derived from the time constant τ and approximately given by the time t END ≈ 10τ. On the other hand, power cables remain in service for several decades of years. 6 Thus, for the reliability assessment of the cable insulation only the steady-state charge and electric field distribution are needed. On the right side in Figure 1, a simplified twodimensional cable geometry is depicted. The cable consists of conductor, insulation, and sheath. Additional layers at the conductor and the sheath that are thin in comparison to the insulation thickness have negligible influence on the charge density. 2 The conductor has the radius r i and the temperature T i , while the sheath has the radius r a with the temperature T a . Because of the symmetry of the cable insulation, the space charge density and the electric field distribution only depend on the radius. On the left side in Figure 1, an example of the transient electric field distribution in space and time is depicted. At t ≈ 0 s, the electric field inside the insulation is determined by the geometry, the permittivity, and the applied external voltage. With increasing time, the electric field changes until a steady-state distribution (∂E ! =∂t ffi 0) is reached. The electric field above the black dotted line at ∂/∂t ffi 0 is called stationary electric field. The long accumulation time of charges in comparison to fast charge dynamics within the insulation (see Delpino et al 8 ), together with a fine time step resolution, usually requires a long computation time within transient electroquasistatic simulation models, using explicit type time integration schemes, until a stationary solution is obtained. 9 An alternative approach to calculate the steady-state charge distribution arises from the solution of the stationary current and the electrostatic potential problem, which are nonlinearly coupled by Ohm's law within the insulator material. 10 The electro-quasistatic field model for the transient computation of space charge and electric field distributions are found in Jeroense et al, 4 whereas the charge distribution calculated by using a fixed point iteration is found in Morshuis et al 11 and Hjerrild et al. 12 The field formulation in Jeroense et al 4 is mathematically equivalent to the electroquasistatic scalar potential field formulation derived directly from the continuity equation (see Clemens et al 13 and Steinmetz et al 14 ). Using the formulation in Jeroense et al, 4 the needed electric field and space charge distribution are computed separately in every time step, whereas only the electric potential is computed, utilizing the electro-quasistatic approximation. The space charge density is then derived from Gauss law.
In this paper, a comparison of the transient and the fixed point computation, to obtain the steady-state charge and electric field distribution, is discussed. Both the computation time and the accuracy are compared. In Section 2, the computation of the steady-state charge density and the electric field by a transient and a fixed point iteration is shown. In Section 3, oscillations that occur during the fixed point iteration are minimized, by introducing an underrelaxation parameter ω. Numerical tests on a cylindrical and a plane parallel insulation of polyethylene (PE), to obtain the computation time, are presented in Section 4, followed by a conclusion in Section 5.
where J is the leakage current density and T is the temperature. 15 To solve (2), the boundary condition at the conductor is the constant applied voltage φ(r i ) = U and at the sheath φ(r a ) = 0. At moderate temperatures, heat are neglected, and the temperature is obtained by the solution of div λgradT where λ is the thermal conductivity. 16 With the assumption of time independent and known conductor and sheath temperatures (T i and T a ), the boundary conditions are T(r i ) = T i at the conductor and T(r a ) = T a at the sheath. The electric conductivity of mass-impregnated paper (MI) and polyethylene (PE) is given by where κ 0 , α, and β are constants and determined by measurements. Furthermore, we assume that κ is not a tensor. 4,17 The spatially discretized formulations of (2) to (4) are and where q is the vector of discrete charges, K, whose entries depend on the vector φ e of electric potentials and on the vector φ T of temperatures, and K λ are the stiffness matrices, M ε is the mass matrix, and b e and b T contain the boundary conditions or sources. The discrete divergence matrix is G T and the discrete gradient matrix is G (see Clemens et al 13 and Steinmetz et al 14 ). Equation (7) is discretized in time using either an implicit or explicit Euler method. The explicit Euler method allows a faster computation, in comparison to implicit time integration schemes, using different acceleration techniques (see Richter et al 18,19 ). Furthermore, the electric field and space charge distributions only depend on the radial direction. Due to the one-dimensional problem, an explicit Euler method allows the fast computation of many time steps and the benefit of an implicit computation is reduced, because of solving a nonlinear system. The vector of electric potentials at the time step n and space charges at the time step n + 1 is obtained by and where n is the time stepping index. The time independent temperature is computed by In (10), Δt is the fixed time step length that must satisfy the Courant-Friedrichs-Lewy (CFL) condition Δt ≤ Δt CFL ≔ Δr/v, where v is the drift velocity of mobile carriers and Δr is the step size of the spatial discretization. 20 The drift velocity v = μÁ j E ! j depends on the magnitude of the electric field and the mobility of charge carriers μ. 21 A high mobility value results in a high velocity value and the charges need less time for the distance of one mesh cell Δr.
To resolute a process within Δr, the time step is chosen to be lower than Δt CFL . The computation stops if either t = t END or kq n + 1 − q n k/kq n + 1 k < ϵ min , where ϵ min ( 1 is the stop threshold. In Algorithm , the computation of (9) to (11) is depicted. The computation starts with the assumption of a vanishing space charge density at t = 0 s. Thus, q n = 0 = 0 C/m 3 and the electric potential φ n = 0 e is the solution of (9) and the temperature is obtained by (11). Using the electric potential and the temperature, the stiffness matrix K φ n = 0 e , φ T À Á is computed, to update the vector of space charges in the next time step n + 1 with (10).
Algorithm 1: Transient computation of space charge and electric field distribution, until steady-state.
1: initialize iteration counter: n = 0. 2: set initial values at time t n = 0 s : where the steady-state charge distribution is obtained by a pseudo time-stepping of (12) and κ E ! jT . With the discrete formulation, the electric potential at n + 1 is obtained by and the space charge density at n + 1 is computed by (6). Algorithm stops if kq n + 1 − q n k/kq n + 1 k < ϵ min or kq n + 1 − q n k/kq n + 1 k > ϵ max is true, where ϵ min, max are stop thresholds. The criteria ϵ min ( 1 consider the convergence and ϵ max ) 1 the divergence of the computed steady-state space charge density during the pseudo time stepping.To speed up convergence, (13) is commonly solved by Newton method. Due to the time-dependent sources (q in (6)), the simulations are compared against an analytic solution or measurements, to ensure accurate results.

| OSCILLATIONS AT THE FIXED POINT ITERATION
The nonlinear electric conductivity is very sensitive to electric field and temperature values. During the nonlinear fixed point iteration process, high electric field distributions might occur especially at the first iterations, resulting in oscillations. These oscillations result in a diverging electric field. For a more detailed analysis of these oscillations, the analytic solutions of the electric field and the space charge density of a one dimensional plane parallel insulation, with a thickness of L, is utilized (see Figure 2). Combining Gauss' law with Ohm's law and using the material equation of the electric flux density D ! = εÁE ! , the stationary space charge density is given by ρ , assuming a constant permittivity ε. Applying the chain rule and using (12), the charge distribution is reduced to where the electric conductivity is given by (5) and the temperature inside a plane parallel insulation is given by With the initial charge density ρ n = 0 = 0 C/m 3 , the electric field E = U/L and β dE x ð Þ dx = 0. Utilizing (14), the space charge distribution after the first iteration (n = 1) is constant and given by The charge density increases with the temperature gradient and with the applied voltage U. Using (2), the electric field after the first iteration (n = 1) is with At the conductor (x = 0), the electric field is below the initial electric field E = U/L, but at the sheath (x = L), the electric field is increased by 1 2 KL, which will also increase the electric conductivity. High electric field values result in oscillations that are seen on the left side in Figure 3, where a 300 μm thick mass-impregnated paper insulation, with κ 0 = 1 Á 10 −16 S/m, α = 0.1 C -1 , β = 0.03 Á 10 −6 m/V, ε r = 3.5 and an applied voltage of 9 kV, is simulated (see Jeroense et al 4 ). The temperature of the sheath is T a = 35 C and the temperature of the conductor is T i = 50 C. The oscillations start to increase after 15 iterations at the sheath and the electric field diverges. The high electric field values increase the conductivity and the space charge density (14) that increases again the electric field. A reduction of the applied voltage results in a decreased electric field and decreased oscillations. Decreasing the voltage from 9 to 8 kV, the electric field converges (see Figure 3 on the right side). With (14), further possibilities to reduce the oscillations are a reduced temperature gradient and reduced constants α and β. After the first iteration (n = 1), the electric field E(x) depends on the coordinate x. Thus, β dE x ð Þ dx is considered in (14) to compute the space charge density after the second iteration (n = 2) and further iterations.
, equal results of (14) to (18) are also obtained for a cylindrical geometry.
To minimize the oscillations, an underrelaxation parameter ω, with ω ∈ [0,1], is introduced. For example, the electric potential at the new iteration step n + 1 is obtained by The parameter ω is used to reduce the computed electric potential, and thus, the electric field in every iteration, during the pseudo time stepping. To obtain a stationary solution, more iterations are required depending on the choice of the underrelaxation parameter ω, but the result does not diverge. As seen in Figure 4 on the left side, after the first iteration, the maximum electric field at the sheath is approximately 53 kV/mm ( Figure 4 on the right side). The high electric field increases the electric conductivity at the sheath and finally results in an unstable result. Using ω = 0.9, the electric field at the sheath is slightly reduced (50.5 kV/mm) after the first iteration and convergence is obtained (see Figure 4 on the right side).
Due to the nonlinear dependency of the electric field on the electric conductivity, the optimal choice of ω depends on the constants α and β and thus, on the insulation material. For PE, the constants are κ 0 = 0.297 Á 10 −19 S/m, α = 0.094 C -1 , β = 0.121 Á 10 −6 m/V and ε r = 2.3, respectively (see Riechert et al 17 and Riechert 22 ). Comparing a F I G U R E 3 . Electric field in an 300 μm thick mass-impreganted insulation. With U = 9 kV, the oscillations increase after 15 iterations and the electric field diverges (left). With U = 8 kV, the electric field converges (right) polyethylene and mass-impregnated paper insulation, the constant α is approximately equal, but the constant β is approximately four times higher for polyethylene in comparison to mass-impregnated paper. Thus, the oscillations start at lower electric field values in comparison to mass-impregnated paper. In general, ω depends on the applied voltage, the temperature gradient, the constants α and β, the spatial resolution Δr, the insulation thickness and the geometry of the insulation (plane parallel or cylindrical). The dependancy of ω on the applied voltage and the temperature gradient is seen in Figure 5, where an 300 μm insulation of MI and PE, with a spatial resolution of Δr = 1 μm, is simulated. Due to the higher constant β of PE in comparison to MI, for example at U = 10 kV and T i − T a = 5 C, ω = 0.05 (PE) and ω = 0.9 (MI). Increasing the temperature gradient T i − T a = 15 C, ω = 0.04 (PE) and ω = 0.45 (MI).

| NUMERICAL VALIDATION
To validate both approaches, depicted in Algorithms and , a PE cable model, where the geometry is shown on the right side in Figure 1, is simulated. The symmetry of the cable results in a radial dependency of the electric field, the temperature, and the space charge density. An analytic expression is developed by rewriting the dependency of the electric field in (5).
where the average field strength E = U r a −r i , E 0 = E=exp 1 ð Þ and γ = β E (see Eoll 23 ). Using the approximation in (20), the one-dimensional stationary electric field (see Mazzanti et al 2 ) inside a cable insulation is given by F I G U R E 4 . Comparison of the electric field in an 300 μm thick massimpreganted insulation. With U = 9 kV and ω = 1, the electric field diverges (left). With U = 9 kV and ω = 0.9 the electric field converges F I G U R E 5 . Dependency of ω on the applied voltage and the temperature gradient for mass-impregnated paper and PE. The temperature gradient is with a one-dimensional time-independent temperature distribution The sheath temperature is T a = 35 C and the conductor temperature is T i = 50 C. The electric field and the space charge density are obtained by applying the finite difference method. The time step Δt that satisfies the CFL condition is obtained by Δt = Δ r/v. The velocity is approximated with v = Eμ. The mobility of charge carriers μ = 1 Á 10 −15 m 2 /Vs, the applied voltage U = 450 kV, the inner radius is r i = 23.2 mm, and the outer radius is r a = 42.4 mm. 4 C, the chosen value for t END is a reasonable assumption to obtain the stationary electric field (see Mazzanti et al 2 ). The stationary electric field, using the transient computation and the pseudo time stepping, is seen in Figure 6 on the left side, together with the analytic expression (21). The error between each numerical computation and the analytic solution is seen on the right side in the same figure. Figure 6 shows practically negligible differences between both numerical computations and the analytic solution. On the right side in Figure 6, the error indicates approximately equal values for both numerical computations of about −0.11% in the insulation. Furthermore, using the pseudo time stepping, the error is higher at the conductor and the sheath but lower within the bulk, in comparison to the transient computation. To compute the stationary solution, 30 iterations and ω = 0.1 are used. Depending on the resolution of the problem, a small time step results from a high mobility of charge carriers μ and a high electric field. Together with the low conductivity of the insulation material and thus, a high time t END until the stationary solution is obtained, the computation time is also high. For example, two more test problems are simulated to determine the increasing computation time, with increasing electric field. The corresponding measurements of two one-dimensional space charge distributions are taken from Bodega et al 25 and Wu et al. 26 The constants for the simulations are summarized in Table 1A. In Table 1B, the number of iterations and the computation time to obtain the stationary charge and electric field distribution are shown. The fixed point iteration shows a faster computation, in comparison to the transient computation for both problems. The fixed point iteration is faster by one order of magnitude, depending on the electric field and the spatial resolution. Algorithms and are equal due to their number of computations during one iteration process. High electric fields and fast charge dynamics are considered with a small time step in the explicit time integration scheme. This leads to more time steps, until F I G U R E 6 Comparison of the stationary electric field, computed by the transient and the fixed point iteration, with the analytic solution (left). Error between both computations and the analytic solution (right) steady-state is reached or a higher calculation time. The fast charge dynamic phase is neglected with the fixed point iteration for the steady-state solution, which leads to a lower calculation time.

| CONCLUSIONS
Two methods for the simulation of space charge distributions in high voltage insulations were compared. The transient simulation was shown to be slower, but no oscillations were observable during the iteration. The nonlinear fixed point iteration approach was faster, but at high electric field intensities, oscillations may occur within the iteration process. Introducing an underrelaxation parameter ω was used to reduce the calculated electric field and space charge density in every iteration, which decreased the oscillations but increased the number of iterations. Comparisons of the fixed point iteration scheme with the transient solution scheme showed a faster calculation of the steady-state solution with the nonlinear iteration scheme by one order of magnitude, depending on the applied electric field and the spatial resolution and thus, on the time step (CFL condition).