Electromagnetic dispersion modeling and sensitivity analysis for HVDC power cables

(1)

Research Report

Stefan Gustafsson, Sven Nordebo, Börje Nilsson 2012-11-13

URI: urn:nbn:se:lnu:diva-22296

Electromagnetic dispersion modeling and

sensitivity analysis for HVDC power cables

(2)

Electromagnetic dispersion modeling and sensitivity analysis for HVDC power cables

Stefan Gustafsson, Sven Nordebo and B¨ orje Nilsson

Abstract

This paper addresses electromagnetic wave propagation in High Voltage Direct Current (HVDC) power cables. An electromagnetic model, based on long (10 km or more) cables with a frequency range of 0 to 100 kHz, is derived.

Relating the frequency to the propagation constant a dispersion relation is formulated using a recursive approach. The propagation constant is found numerically with normalized residue calculation. The paper is concluded with a sensitivity analysis of the propagation constant with respect to the electrical parameters _r (the real relative permittivity) and σ (the conductivity).

1 Introduction

The growing environmental awareness, the need of low maintenance costs, high availability and few faults [1] has led to an increasing demand for High Voltage Direct Current (HVDC) power cables and they are today an important part of the power transmission infrastructure. Delivering power from off-shore wind power plants, feeding oil/gas platforms with power and cost-effective transporting of power over long distances are all examples where power cables play an important role [2].

Once powered on and in service there are many areas of applications where a study in signal transmission and corresponding dispersion is of interest. Such areas of applications could be power line communication (PLC) techniques, simulation and analysis of lightning and switching overvoltages, transient based protection, fault localization and Partial Discharge (PD) surveillance.

High voltage cables are designed to transport very high currents and voltages, and in order to do so they need to be designed in a very specific manner [3], see figure 1. The innermost layer is a copper or aluminum conductor made out of compact stranded round wires. The next three layers is the actual insulation system, where the conductor- and insulation screen are made out of a semi-conductive polymer and the insulation is made out of a cross linked HVDC polymer. Outside this insulation system there is a metallic layer, in this case a lead alloy sheath. The lead sheath serves as a radial water sealing. Over the lead sheath there is a polyethylene sheath, the inner jacket, and it is designed to act as a mechanical protection for

(3)

the lead sheath. Another mechanical protection is the tensile armour which consists of round galvanized steel wires. The purpose of the armour is to protect the cable during installation and against impacts and abrasion if the cable is not burred to a safe depth in the seabed. The outermost layer is a polypropylene yarn serving as a semi-wet covering.

In addition to the layers mentioned above there are other thin layers in the cable, some of them made of swelling material acting as a longitudinal water sealing, but these layers are assumed to be of less importance when constructing a mathematical model of the cable.

Figure 1: A cutaway drawing of a submarine cable.

Each layer in the cable has its own specific electrical properties. To separate them from each other they are indexed from 0 to 8, see figure 2 and table 1 below, where r is the real relative permittivity, σ is the conductivity and ρ represents the distance from the center of the cable.

The electromagnetic model presented in this paper is focused on very long (10 km or more) HVDC power cables and the frequency range is about 0-100 kHz. The derivation of the model starts with Maxwell’s equations [4, 5]. The electromagnetic fields are then expanded in cylindrical vector waves, see Appendix A. The model is based on the first propagating mode, the TM01 mode (the quasi-TEM mode), but

(4)

with a brief comparison with the second propagating mode, the TM₀₂mode, see figure 5 below. The model includes the boundary conditions related to the boundaries between the different layers in the cable. The field components corresponding to these boundary conditions are assembled into a square matrix A(γ) and this matrix dictates the values of the propagation constant γ in the cable. For each frequency the value of γ is found by giving it an initial guess, based on a visual estimation of γ in the complex plane (see figure 3), and then computed numerically by using residue calculus.

The last part in this paper addresses the question to which of the electrical parameters have the most influence on the propagation constant, and hence which parameters, if not correctly chosen, could affect the electromagnetic model nega- tively. Consequently, a sensitivity analysis of the propagation constant is made with respect to the electrical parameters _r and σ.

_r0

_r1

_r2

_r3

_r4

_r5

_r6

_r7

_r8 σ₀ σ₁ σ₂ σ₃ σ₄ σ₅ σ₆ σ₇ σ₈ ρ0

ρ1

ρ2

ρ3

ρ4

ρ₅ ρ6

ρ7

Figure 2: Geometrical and material definitions of the multi-layered cable. Addi- tional information regarding the parameters can be found in table 1.

(5)

2 The electromagnetic model

2.1 The wave equations

In the frequency domain the source free Maxwell’s equations can be expressed as

∇ × E = −iωµ₀µH

∇ × H = iωε₀E, (1)

where E and H are the electric- and magnetic fields, ω = 2πf is the angular frequency and f is the frequency. The permeability of free space, the permittivity of free space and the real relative permeability are denoted µ₀, ε₀ and µ, respectively.

The complex relative permittivity is denoted and is represented as

= _r− iσ

kη₀, (2)

where r is the real relative permittivity, σ is the conductivity; the wave number of free space is k = ω/c₀, c₀ = 1/√

µ₀ε₀ is the speed of light in free space, and η₀ = pµ0/ε₀ is the wave impedance of free space. The time-convention used here is e^iωt.

Assuming a source free environment, and by applying the ∇×-operator on the equations in (1), the vector Helmholtz wave equations can be described as

(∇²+ k²µ)E = 0

(∇²+ k²µ)H = 0. (3)

Since the cable, and thus the region of interest, is cylindrical in its structure, it is appropriate to represent the study in cylindrical coordinates, here denoted by (ρ, φ, z) and with the transversal coordinate vector ρ = ρˆρ. If we assume that the electromagnetic wave in the cable is propagating in the z-direction, then we let the space-convention of that propagating wave be given by e^−γz, where γ is the propagation constant representing a particular mode [4, 5]. The electric- and magnetic fields E and H can then be written as E = E(ρ, γ)e^−γz and H = H(ρ, γ)e^−γz. Equation (3) will then become

(∇²_t + k²µ + γ²)E(ρ, γ) = 0

(∇²_t + k²µ + γ²)H(ρ, γ) = 0, (4) where ∇²_t is the transversal part of the Laplace operator.

2.2 The boundary conditions and the dispersion relation

The electromagnetic field components, at the boundary between two layers in the cable, not including ρ = 0 or ρ = ∞, are subjected to the appropriate boundary conditions. Each circle in figure 2 is a boundary between two different layers. The

(6)

electromagnetic fields on each side of the boundary are related to each other by the continuity of the tangential components, as [4, 5]

( n × (Eˆ 1 − E2) = 0 ˆ

n × (H₁− H₂) = 0, (5)

where ˆn is the unit normal, E₁ and H₁ are the fields on one side of the boundary, and E2 and H2 are the fields on the other side of the boundary.

Before employing the boundary conditions we note that the electric- and magnetic fields in equation (4) originate from

E = 1 2πi

Z i∞

−i∞

" _∞ X

m=−∞

ˆ

ρE_ρm+ ˆφE_φm+ ˆzE_zm

#

e^imφe^−γzdγ (6) and

H = 1

2πωµ₀µ Z i∞

−i∞

" _∞ X

m=−∞

ˆ

ρH_ρm+ ˆφH_φm+ ˆzH_zm

#

e^imφe^−γzdγ. (7) The field components in equation (6) and (7) can be derived from a general discussion regarding cylindrical vector wave expansion, see Appendix A. The tangential fields are











E_φm = _2πi¹

−a_1mH⁽¹⁾⁰m (κρ) − b_1mH⁽²⁾⁰m (κρ) − a_2m^imγ_kκρH⁽¹⁾m (κρ) − b_2m^imγ_kκρH⁽²⁾m (κρ) E_zm = _2πi¹

a_2m^κ_kH⁽¹⁾m (κρ) + b_2m^κ_kH⁽²⁾m (κρ) H_φm= _2πη¹

0µ

−a_1m^imγ_kκρH⁽¹⁾m (κρ) − b_1m^imγ_kκρH⁽²⁾m (κρ)

−a_2mµH⁽¹⁾⁰m (κρ) − b_2mµH⁽²⁾⁰m (κρ) H_zm= _2πη¹

0µ

a_1m^κ_kH⁽¹⁾m (κρ) + b_1m^κ_kH⁽²⁾m (κρ) ,

(8) The Hankel function H^(j)0m (κρ), of order m, is of first kind (j = 1) or second kind (j = 2), and the prime sign represents a differentiation with respect to the argument.

The focus of this paper is concerned with the Transversal Magnetic modes of order m = 0, here denoted TM_0n. The most important mode is the first mode, TM₀₁, which is also called the quasi-TEM mode. All the other modes of order m 6= 0, including the Transversal Electric modes, TE_0n, will essentially be cut-off when the frequency is in the lower domains. The only field components that will remain non-zero after this treatment is E_z0, H_φ0 and E_ρ0, where the first two can be found in equation (8) if m is set to 0. Now, let us formulate the boundary conditions in the cable, starting with the boundary indicated ρ₀ in figure 2. The materials on each side of the boundary differ in characteristics, and to separate them from each other they will be indexed with an i, where i = 0, 1, ..., N + 1. That is, the complex valued permittivity and the radial wave number is defined by

_i = _{r i}− iσi

k η₀, i = 0, 1, ..., N + 1 (9)

(7)

and

κ_i =p

k²_iµ_i+ γ², i = 0, 1, ..., N + 1. (10) Hence, the boundary conditions at ρ₀, using equation (5), (8), (9) and (10), are given by

( a₁κ₁H⁽¹⁾₀ (κ₁ρ₀) + b₁κ₁H⁽²⁾₀ (κ₁ρ₀) − a₀κ₀J₀(κ₀ρ₀) = 0 a₁₁H⁽¹⁾₁ (κ₁ρ₀) + b₁₁H⁽²⁾₁ (κ₁ρ₀) − a₀₀J₁(κ₀ρ₀) = 0.

(11) The coefficients a_2m and b_2m in equation (8) are now expressed as a₁, b₁ and a₀ where the index tells us in which layer the coefficient belongs. In equation (11) the derivative of the Hankel- and Bessel functions are written as −H⁽¹⁾⁰₀ = H⁽¹⁾₁ ,

−H⁽²⁾⁰₀ = H⁽²⁾₁ and −J⁰₀ = J₁, respectively. Note also that ₀ is the complex valued permittivity in layer 0, not to be mixed up with the permittivity constant denoted ε₀ in equation (1). Derived in the same way as above the intermediate boundaries have boundary conditions given by











a_iκ_iH⁽¹⁾₀ (κ_iρ_i−1) + b_iκ_iH⁽²⁾₀ (κ_iρ_i−1)

−a_i−1κ_i−1H⁽¹⁾₀ (κ_i−1ρ_i−1) − b_i−1κ_i−1H⁽²⁾₀ (κ_i−1ρ_i−1) = 0 a_i_iH⁽¹⁾₁ (κ_iρ_i−1) + b_i_iH⁽²⁾₁ (κ_iρ_i−1)

−a_i−1_i−1H⁽¹⁾₁ (κ_i−1ρ_i−1) − b_i−1_i−1H⁽²⁾₁ (κ_i−1ρ_i−1) = 0

(12)

where i = 2, 3, ..., N . The last boundary conditions related to the boundary at radius ρ_N are given by

( b_{N +1}κ_{N +1}H⁽²⁾₀ (κ_{N +1}ρ_N) − a_Nκ_NH⁽¹⁾₀ (κ_Nρ_N) − b_Nκ_NH⁽²⁾₀ (κ_Nρ_N) = 0 bN +1N +1H⁽²⁾₁ (κN +1ρN) − aNNH⁽¹⁾₁ (κNρN) − bNNH⁽²⁾₁ (κNρN) = 0.

(13)

The field components of the boundary conditions in equation (11), (12) and (13) can all be assembled into a square matrix A(γ) of size (2N + 2) × (2N + 2), and the corresponding coefficients, a_i and b_i where i relates to the layer, can be expressed as a vector x. The boundary conditions can now be written as A(γ)x = 0, to which we want to find a non-trivial solution. One way to do that is to solve

detA(γ) = 0, (14)

which is referred to as the dispersion relation since it relates the propagation constant γ to the frequency.

2.3 The recursive form

In this subsection we will derive an explicit expression for the dispersion relation mentioned in equation (14), and one way to do so is to represent the dispersion relation in a recursive way. Functions that we will use more frequently are, for the

(8)

sake of convenience, defined as











a_i(κ_i) = H⁽¹⁾₁ (κ_iρ_i−1)H⁽²⁾₀ (κ_iρ_i) − H⁽²⁾₁ (κ_iρ_i−1)H⁽¹⁾₀ (κ_iρ_i) b_i(κ_i) = H⁽²⁾₀ (κ_iρ_i−1)H⁽¹⁾₀ (κ_iρ_i) − H⁽¹⁾₀ (κ_iρ_i−1)H⁽²⁾₀ (κ_iρ_i) c_i(κ_i) = H⁽¹⁾₁ (κ_iρ_i−1)H⁽²⁾₁ (κ_iρ_i) − H⁽²⁾₁ (κ_iρ_i−1)H⁽¹⁾₁ (κ_iρ_i) d_i(κ_i) = H⁽²⁾₀ (κ_iρ_i−1)H⁽¹⁾₁ (κ_iρ_i) − H⁽¹⁾₀ (κ_iρ_i−1)H⁽²⁾₁ (κ_iρ_i),

(15)

where i ≥ 1.

Further, let us define two determinants of order zero, f₀(γ, k) and g₀(γ, k), as







f0(γ, k) = −J0(κ0ρ0) g₀(γ, k) = −k₀

κ₀ J₁(κ₀ρ₀), (16)

and in the same way, the determinants of first order, as

f₁(γ, k) =

−J₀(κρ) κ₁H⁽¹⁾₀ (κ₁ρ₀) κ₁H⁽²⁾₀ (κ₁ρ₀)

−k₀

κ₀ J₁(κ₀ρ₀) k₁H⁽¹⁾₁ (κ₁ρ₀) k₁H⁽²⁾₁ (κ₁ρ₀) 0 −κ₁H⁽¹⁾₀ (κ₁ρ₁) −κ₁H⁽²⁾₀ (κ₁ρ₁)

, (17)

and

g₁(γ, k) =

−J₀(κρ) κ₁H⁽¹⁾₀ (κ₁ρ₀) κ₁H⁽²⁾₀ (κ₁ρ₀)

−k₀ κ0

J₁(κ₀ρ₀) k₁H⁽¹⁾₁ (κ₁ρ₀) k₁H⁽²⁾₁ (κ₁ρ₀) 0 −k₁H⁽¹⁾₁ (κ₁ρ₁) −k₁H⁽²⁾₁ (κ₁ρ₁)

. (18)

A recursive representation of the first order determinants, using equation (15) and (16), can be written as

( f₁(γ, k) = −k₁κ₁a₁(κ₁)f₀(γ, k) − κ²₁b₁(κ₁)g₀(γ, k)

g₁(γ, k) = −k²²₁c₁(κ₁)f₀(γ, k) − k₁κ₁d₁(κ₁)g₀(γ, k), (19) where the determinants (17) and (18) have been evaluated along their last row.

Higher order determinants can be derived in the same way, thus

f_i(γ, k) =

−J0(κ0ρ0) κ1H⁽¹⁾₀ (κ1ρ0) κ1H⁽²⁾₀ (κ1ρ0)

−k₀

κ₀ J₁(κ₀ρ₀) k₁H⁽¹⁾₁ (κ₁ρ₀) k₁H⁽²⁾₁ (κ₁ρ₀)

. .. . ..

−x⁽¹⁾_i−1,i−1 −x⁽²⁾_i−1,i−1 x⁽¹⁾_i,i−1 x⁽²⁾_i,i−1

−y_i−1,i−1⁽¹⁾ −y_i−1,i−1⁽²⁾ y_i,i−1⁽¹⁾ y_i,i−1⁽²⁾

−x⁽¹⁾_i,i −x⁽²⁾_i,i

(20)

(9)

and

gi(γ, k) =

−J₀(κ₀ρ₀) κ₁H⁽¹⁾₀ (κ₁ρ₀) κ₁H⁽²⁾₀ (κ₁ρ₀)

−k₀

κ₀ J₁(κ₀ρ₀) k₁H⁽¹⁾₁ (κ₁ρ₀) k₁H⁽²⁾₁ (κ₁ρ₀)

. .. . ..

−x⁽¹⁾_i−1,i−1 −x⁽²⁾_i−1,i−1 x⁽¹⁾_i,i−1 x⁽²⁾_i,i−1

−y⁽¹⁾_i−1,i−1 −y⁽²⁾_i−1,i−1 y_i,i−1⁽¹⁾ y_i,i−1⁽²⁾

−y_i,i⁽¹⁾ −y_i,i⁽²⁾

(21)

are the ith order determinants where ( x⁽¹⁾_i,j = κ_iH⁽¹⁾₀ (κ_iρ_j)

y_i,j⁽¹⁾ = k_iH⁽¹⁾₁ (κ_iρ_j) ,

( x⁽²⁾_i,j = κ_iH⁽²⁾₀ (κ_iρ_j) y⁽²⁾_i,j = k_iH⁽²⁾₁ (κ_iρ_j)

(22)

and i ≥ 2 and j = i − 1. Following the recursive relation in (19) the ith order determinants can now be derived as

( f_i(γ, k) = −k_iκ_ia_i(κ_i)f_i−1(γ, k) − κ²_ib_i(κ_i)g_i−1(γ, k)

g_i(γ, k) = −k²²_ic_i(κ_i)f_i−1(γ, k) − k_iκ_id_i(κ_i)g_i−1(γ, k). (23) The recursive relation (23) is valid for 1 ≤ i ≤ N , but we also have to take the last boundary condition, the boundary between the outer layer and the external region, into consideration. The complete determinant, detA(γ) from equation (14), considering all the boundaries and all the layers, including the external region, can then be written as

detA(γ) = k_{N +1}H⁽²⁾₁ (κ_{N +1}ρ_N)f_N(γ, k) − κ_{N +1}H⁽²⁾₀ (κ_{N +1}ρ_N)g_N(γ, k). (24)

2.4 The algorithm, finding the zeros

The first step towards finding a zero to equation (14) is to plot the dispersion relation in a complex plane. The zeros can then be found visually and given an approximate complex value. A typical example of such a plot can be seen in figure 3, where the upper graph shows a plot of log |detA(γ)| and the lower graph shows a plot of its arguments. In both the graphs two zeros can be found. The first zero, closer to the origin, is the quasi-TEM mode which we will call TM₀₁, while the other zero represents the next mode and will be called TM₀₂.

Suppose that an approximate value ˆγ1 of the true propagation constant γ1 is given. The exact value of the propagation constant γ₁ can then be found by normalized residue [6],

γ₁ = I

C

γ detAdγ I

C

1 detAdγ

. (25)

(10)

5 10 15 20 0

5 10 15 20

20 20.5 21 21.5 22 22.5 23

5 10 15 20

0 5 10 15 20

−2

−1 0 1 2

Im [γ] /k

Re [γ] · 2 · 10⁶· log e (dB/100km)

Im [γ] /k

Re [γ] · 2 · 10⁶ · log e (dB/100km)

log |detA(γ)| arg detA(γ)

Figure 3: The left graph shows a plot of log |detA(γ)| and the right graph shows a plot of its arguments. On the horizontal axis: Re[γ] · 2 · 10⁶log e, i.e. the attenuation of the mode expressed in dB/100km. On the vertical axis: Im[γ]/k. The Quasi-TEM mode, TM01, is located close to the origin, and the TM02 mode more to the right.

The electric properties of the cable used in this plot can be found in table 1. The frequency is 50 Hz.

A contour C is set up around ˆγ₁, where the closed loop C must circumscribe the true value γ₁, and no other zeros of A(γ) is allowed inside the loop. Equation (25) is applied when computing numerical examples in MATLAB. Computing the exact value of the propagation constant using normalized residue is a stable but slow method. A faster procedure is the secant method, but this method on the other hand is not that stable. So, a combination of the two methods could be a way to speed up the calculations.

The zeros in figure 3 will change position when the frequency is changed. So, for each new frequency follows a new complex value of the propagation constant γ₁, and in order to find it we need a new initial guess of ˆγ₁. But, instead of once again doing a new plot and visually find the initial guess of the propagation constant, it is now possible to make use of the earlier calculated value of the propagation constant.

This process will only be helpful if the change of the frequency is relatively small.

A numerical algorithm to find and compute the different zeros when changing the frequency can look like this:

1. Plot the dispersion function for a specific frequency, k₀. Choose which zero to compute, and give that zero an initial guess of its propagation constant, ˆγ₁. 2. Circumscribe ˆγ₁ with a contour C small enough to avoid other zeros and

possible branch-cuts.

3. Compute the true value of the propagation constant, γ₁, by using equation (25). End the algorithm here or move on to step 4 if a new γ₁ with a new frequency, close to the previous one, is requested.

4. Change the frequency a relatively small step, ∆k. Determine an initial guess ˆ

γ1(k + ∆k)¹ based on the previous frequency, k, in step 3. Return to step 2 above, or end the algorithm.

1The real part of the initial guess ˆγi(ki−1+ ∆k) has the same real part as γi−1(ki−1). The

(11)

Using the algorithm above on the first propagating mode, TM₀₁, in a frequency range of 0-100 kHz, yields the result shown in figure 4. The upper plot relates the frequency to the real part of the propagation constant and the lower plot relates the frequency to the normalized wave velocity, where the speed of light is set to 1. The exact electrical properties in all of the layers in the cable are difficult to determine, instead focus has been added on properties that we believe are of most importance.

The modeling parameters used in the calculations can be seen in table 1.

It is also interesting to compare the first and the second mode with each other, i.e. the two modes TM₀₁ and TM₀₂. Let the frequency range be 0-1000 Hz, using the algorithm above the propagation constants for the two different modes can be computed and plotted in the same graph, see figure 5. The choice to focus only on the first mode, TM₀₁, in this paper is supported by the big difference in the attenuation between the two modes at higher frequencies².

0 2 4 6 8 10

x 10⁴ 0

10 20 30 40

0 2 4 6 8 10

x 10⁴ 0.2

0.3 0.4 0.5 0.6 0.7

Attenuation, dB/100km

f [ Hz]

Normalized wave velocity

f [ Hz]

Figure 4: Upper plot: on the vertical axis: Re[γ] · 2 · 10⁶log e, i.e. the attenuation of the TM01 mode expressed in dB/100km. On the horizontal axis: the frequency, 0-100 kHz. Lower plot: the normalized wave velocity, c/c0= k/Im[γ(k)] of the TM01

mode on the vertical axis, and the frequency 0-100 kHz on the horizontal axis.

imaginary part of the initial guess ˆγ_i(k_i−1+ ∆k) is computed as Im[γ_i]·((k_i−1+ ∆k)/k_i−1), and this is based on the almost linear relation between k and Im[γ].

2At lower frequencies the TM₀₂ mode may have more influence over the propagating electromagnetic wave.

(12)

Layer R[mm] Permit. Conduct. Permeab.

Inner conductor ρ₀= 24.3 _r0 = 1 σ₀= σ_Cu= 5.8 · 10⁷ µ₀= 1 Conductor screen ρ₁= 26.0 _r1 = 2.3 σ₁= σ_s = 1 µ₁= 1 Insulation ρ₂= 42.0 _r2 = 2.3 σ₂= 0 µ₂= 1 Insulation screen ρ₃= 43.9 _r3 = 2.3 σ₃= σ_s = 1 µ₃= 1 Lead sheath ρ₄= 46.9 _r4 = 1 σ₄= σ_{P b}= 4.6 · 10⁶ µ₄= 1 Inner jacket sheath ρ₅= 49.5 _r5 = 2.3 σ₅= 0 µ₅= 1 Armour ρ₆= 53.5 _r6 = 1 σ₆= σ_{F e} = 1.1 · 10⁶ µ₆= 1*

Outer serving ρ₇= 58.5 _r7 = 2.3 σ₇= 0 µ₇= 1 Exterior region (air) ρ₈= ∞ _r8 = 1 σ₈= 0 µ₈= 1 Table 1: Electric properties of the layers in the cable. The subscript indicates in which layer the parameter is located, e.g. ρ₀, _r0, σ₀and µ₀all belong to layer 0, i.e.

the conductor, the innermost layer. The R[mm] is the radius, Permit. is the real relative permittivity, Conduct. is the conductivity and Permeab. is the relative permeability.

* In this paper a value of 1 is used. A more accurate value is approximately 100.

0 200 400 600 800 1000

0 20 40 60 80

0 200 400 600 800 1000

0 0.2 0.4 0.6 0.8

f [ Hz]

Figure 5: A solid line represent the TM01 mode and a dashed line represent the TM02mode. Upper plot: on the vertical axis: Re[γ] · 2 · 10⁶log e, i.e. the attenuation expressed in dB/100km. On the horizontal axis: the frequency, 0-1000 Hz. Lower plot: the normalized wave velocity c/c0 = k/Im[γ(k)] on the vertical axis, and the frequency 0-1000 Hz on the horizontal axis.

(13)

3 Sensitivity analysis

As was mentioned in the previous chapter it is difficult to determine the exact electrical properties of the layers in the cable. But, in order to know if a small change in one of the electric properties will give a significant change in the propagation constant, it is appropriate to do a parameter sensitivity study. The result of such a study will at least give us guidance in which layers and materials to focus on.

3.1 Computation of sensitivity

Variables of interest when performing the parameter sensitivity study are the real relative permittivity, _r, and the conductivity, σ. The relative permittivity and conductivity of each layer is indexed as in equation (9). Gather all the variables of interest in a parameter vector θ =r0, ..., _{r(N +1)}, σ₀, ...., σ_{N +1} and let γ(θ) represent the propagation constant related to the TM₀₁ mode. The dispersion relation in equation (14) can now be expressed as

detA(γ(θ), θ) = 0, (26)

clarifying the explicit dependence of the complex propagation constant and the parameter vector. If we differentiate equation (26) with respect to one of the variables in the parameter vector θ, the following expression is obtained

∂

∂θ_idetA(γ(θ), θ) = ∂

∂γdetA(γ, θ)∂γ(θ)

∂θ_i + ∂

∂θ_idetA(γ, θ) = 0, (27) and rearranging the equation yields

∂γ(θ)

∂θ_i =

− ∂

∂θ_idetA(γ, θ)

∂

∂γdetA(γ, θ)

. (28)

For each frequency ω, the TM₀₁ mode has its specific propagation constant γ(θ), and these two variables are used when computing equation (28) in MATLAB. In addition to this, one helpful trick when computing equation (28) is to express it using the relation

∂

∂θdetA = tr

M∂A^T

∂θ

, (29)

where M is the cofactor matrix of A. Equation (29) is derived in appendix B.

Equation (28) will give us information about the propagation constant γ(θ) and how it will be perturbed when the material parameter θ_i is changed. If no, or very little, change occurs it can be assumed that θ_i, and its default value, only has a slight impact on the propagation constant γ(θ). On the other hand, if the impact on γ(θ) is large we can assume that θ_i and its associated layer and material plays a larger role in the model of the cable.

(14)

3.2 The conductivity

A sensitivity analysis of the propagation constant γ(θ) with respect to the conductivity σ is of interest. Some of the materials in the cable have none or close to zero conductivity and for that reason they are of less interest. Therefore, focus will essentially be on the materials with conductivity, i.e. σ₀, σ₁, σ₃, σ₄ and σ₆. But, even among these conductivities there are big differences, see table 1. It is then appropriate to look at their relative change which can be expressed as

dγ = ∂γ

∂σdσ = ∂γ

∂σσdσ

σ . (30)

Equation (30) yields that if σ is changed e.g. 10 percent (see figure 15) from its original value, i.e. ∆_relσ = dσ/σ = 0.1, the relative change in dγ can be expressed as

dγ = ∂γ

∂σσ∆_relσ. (31)

Splitting γ in one real and one imaginary part yields dRe{γ} · 2 · 10⁶· log e = ∂Re{γ}

∂σ σ∆_relσ · 2 · 10⁶· log e dIm{γ}1

k = ∂Im{γ}

∂σ σ∆_relσ1 k,

(32)

where the equations have been scaled with 2 · 10⁶· log e (dB/100km) and 1/k respectively to match the scale factors in figure 3, and σ represents either σ_Cu, σ_s, σ_{P b} or σ_{F e} depending on which layer is studied. Hence, the sensitivity functions of interest

are ∂γ

∂σ_Cuσ_Cu, ∂γ

∂σ_sσ_s, ∂γ

∂σ_{P b}σ_{P b}, ∂γ

∂σ_{F e}σ_{F e}, (33)

and when scaled and illustrated in a graph, see figure 6, 7 and 8 below.

(15)

0 1 2 3 4 5 6 7 8 9 10 x 10⁴

−12

−10

−8

−6

−4

−2 0 2

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−2

−1.5

−1

−0.5 0 0.5

2 · 10⁶· log e · ^∂Re[γ]_∂σ · σ

f [ Hz]

1

k · ^∂Im[γ]_∂σ · σ

f [ Hz]

Figure 6: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the conductivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the conductivity. Each colour corresponds to a sensitivity function, σ0= σCu (inner conductor, yellow), σ1= σs(conductor screen, magenta), σ3 = σs (insulation screen, dashed magenta which is hidden behind the solid magenta line), σ4 = σP b (lead sheath, green), σ6 = σF e (armour, black). The frequency range is 0-100 kHz.

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

−3

−2.5

−2

−1.5

−1

−0.5 0

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

−2

−1.5

−1

−0.5 0

2 · 10⁶· log e · ^∂Re[γ]_∂σ · σ

f [ Hz]

1

k · ^∂Im[γ]_∂σ · σ

f [ Hz]

Figure 7: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the conductivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the conductivity. Each colour corresponds to a sensitivity function, σ0= σCu (inner conductor, yellow), σ1= σs(conductor screen, magenta), σ3 = σs (insulation screen, dashed magenta which is hidden behind the solid magenta line), σ4 = σP b (lead sheath, green), σ6 = σF e (armour, black). The frequency range is 0-10 kHz.

(16)

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

−3

−2

−1 0

x 10⁻⁵

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

−2.5

−2

−1.5

−1

−0.5 0 0.5

x 10⁻⁷

2 · 10⁶· log e · ^∂Re[γ]_∂σ · σ

f [ Hz]

1

k · ^∂Im[γ]_∂σ · σ

f [ Hz]

Figure 8: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the conductivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the conductivity. Each colour corresponds to a sensitivity function, σ1= σs(conductor screen, magenta), σ3= σs(insulation screen, dashed magenta). The frequency range is 0-10 kHz.

3.3 The real relative permittivity

A parameter sensitivity study of the propagation constant with respect to the real relative permittivity _r is also of interest. Although the real relative permittivities in the different layers almost have the same values (1 or 2.3), see table 1, we will hold on to the same method, as in the case of the conductivity, and express the sensitivity functions in terms of their relative change. The relative change can then be expressed as

dγ = ∂γ

∂r

_rd_r

r

. (34)

Similarly, as above, γ is split in one real and one imaginary part and the sensitivity function is scaled as was done in equation (32).

Figure 9 shows the situation when all the sensitivity functions with their associated _ri parameters are plotted in the same graph, i = 0, 1, ..., 8. The cyan coloured curve represents the most sensitive real relative permittivity parameter r2 whether we are talking about low or high frequencies. The blue curve, representing _r5, also shows signs of increasing sensitivity at a frequency around 10 kHz. The rest of the sensitivity functions hide behind each other close to the horizontal axis. One con- clusion we can draw from this is that the sensitivity functions differ a lot. In figure 10, 11, 12, 13 and 14 some of the sensitivity functions have been gradually removed for easier interpretation.

(17)

0 1 2 3 4 5 6 7 8 9 10 x 10⁴ 0

5 10 15

0 1 2 3 4 5 6 7 8 9 10

x 10⁴ 0

1 2 3

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 9: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity. Each colour corresponds to a sensitivity function, _r0 (inner conductor, yellow), _r1 (conductor screen, magenta), _r2(insulation, cyan), _r3(insulation screen, dashed magenta), _r4 (lead sheath, green), _r5 (inner jacket sheath, blue), _r6 (armour, black), _r7 (outer serving, red), r8 (exterior region, air, dashed red). Most of the functions can not be seen because they are hidden behind each other close to the horizontal axis, see figure 10−14 for more information. The frequency range is 0-100 kHz.

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−0.5 0 0.5 1

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−0.01 0 0.01 0.02 0.03

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 10: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity.

Each colour corresponds to a sensitivity function, _r5(inner jacket sheath, blue), _r7 (outer serving, red), _r8 (exterior region, air, dashed red). The frequency range is 0-100 kHz.

(18)

0 1 2 3 4 5 6 7 8 9 10 x 10⁴

−0.01 0 0.01 0.02 0.03 0.04

0 1 2 3 4 5 6 7 8 9 10

x 10⁴ 0

0.01 0.02 0.03

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 11: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity. Each colour corresponds to a sensitivity function, _r7 (outer serving, red), _r8 (exterior region, air, dashed red). The frequency range is 0-100 kHz.

0 1 2 3 4 5 6 7 8 9 10

x 10⁴ 0

5 10

x 10⁻⁵

0 1 2 3 4 5 6 7 8 9 10

x 10⁴ 0

2 4 6

x 10⁻⁵

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 12: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity. Each colour corresponds to a sensitivity function, r1 (conductor screen, magenta), r3

(insulation screen, dashed magenta), r7 (outer serving, red). The frequency range is 0-100 kHz.

(19)

0 1 2 3 4 5 6 7 8 9 10 x 10⁴ 0

2 4 6 8

x 10⁻¹⁰

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−20

−15

−10

−5 0

x 10⁻¹²

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 13: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity. Each colour corresponds to a sensitivity function, r1 (conductor screen, magenta), r3

(insulation screen, dashed magenta), r4 (lead sheath, green). The frequency range is 0-100 kHz.

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−3

−2

−1 0 1 2 3

x 10⁻¹²

0 1 2 3 4 5 6 7 8 9 10

x 10⁴

−5 0 5

x 10⁻¹⁵

2 · 10⁶· log e · ^∂Re[γ]_∂

r · _r

f [ Hz]

1

k · ^∂Im[γ]_∂

r · _r

f [ Hz]

Figure 14: Upper plot: the real part of the scaled (2 · 10⁶· log e) sensitivity function with respect to the real relative permittivity. Lower plot: the imaginary part of the scaled (1/k) sensitivity function with respect to the real relative permittivity. Each colour corresponds to a sensitivity function, r0 (inner conductor, yellow), r1 (conductor screen, magenta), _r3 (insulation screen, dashed magenta), _r4 (lead sheath, green), _r6(armour, black). The frequency range is 0-100 kHz.

(20)

3.4 Summary, conclusions and future work

A sensitivity analysis of the propagation constant has been presented in this paper. The propagation constant was derived from the first propagating mode, the quasi-TEM mode, in a HVDC power cable. The analysis was done with respect to some of the electrical parameters in the different layers in the cable. The electrical parameters of interest were the real relative permittivity _r and the conductivity σ, and these had different values and importance depending on which layer they belonged to.

When the conductivities in the layers where compared with each other, results showed that σ4, the conductivity in the lead alloy sheath, was the parameter dis- playing the most sensitivity, i.e. a change in σ₄’s default value had more impact on the propagation constant than if a change was made in any of the other conductivities and their default values. In the same way, the real relative permittivities in the different layers where compared with each other. The parameter most sensitive to a change was _r2, i.e. the real relative permittivity in the insulation layer. Thus, the two most sensitive parameters were σ4 and r2, and this was true for all frequencies in a frequency range of 0-100 kHz.

The fact that the propagation constant, to some extent, does react to a relatively small change, see figure 15, in some of the electrical parameters shows that it is of great interest to find more accurate parameter values.

The relative permeability µ₆ found in the armour was in this paper set to 1. A more likely value is approximately 100, a comparison can be seen in figure 16. Thus, future work and improvements should focus on finding more accurate values of the electrical parameters in the cable.

(21)

0 1 2 3 4 5 6 7 8 9 10 x 10⁴ 0

10 20 30 40

0 1 2 3 4 5 6 7 8 9 10

x 10⁴ 0

0.2 0.4 0.6

f [ Hz]

Figure 15: Upper plot: on the vertical axis: Re[γ] · 2 · 10⁶log e, i.e. the attenuation expressed in dB/100km. On the horizontal axis: the frequency, 0-100 kHz. Lower plot: on the vertical axis: the normalized wave velocity c/c0 = k/Im[γ(k)]. On the horizontal axis: the frequency 0-100 kHz. Blue colour is the original (default) value of the propagation constant. Red colour represents a 10 percent value increase of σ4. Green colour represents a 10 percent value increase of r2. Black colour represents a 10 percent value increase in both the σ4 and the r2 parameter. In the lower graph the blue and red lines are very close to each other, and the same applies to the green and the black lines.

0 0.5 1 1.5 2 2.5 3 3.5 4

x 10⁴ 0

5 10 15 20

0 0.5 1 1.5 2 2.5 3 3.5 4

x 10⁴ 0

0.1 0.2 0.3 0.4 0.5 0.6

f [ Hz]

Figure 16: Upper plot: on the vertical axis: Re[γ] · 2 · 10⁶log e, i.e. the attenuation expressed in dB/100km. On the horizontal axis: the frequency, 0-40 kHz. Lower plot: on the vertical axis: the normalized wave velocity c/c0 = k/Im[γ(k)]. On the horizontal axis: the frequency 0-40 kHz. Blue colour indicates the original (default) value (µ6 = 1) of the propagation constant. The cyan coloured line indicates the propagation constant when µ6 = 100. At higher frequencies, i.e. 40-100 kHz, the lines only differ insignificantly.