Analysis of Harmonic Propagation in Power Systems Using Standing Waves

—In many countries, the proportion of high-voltage cables is increasing to address the issue of public resistance to transmission system extensions. A major technical challenge in this is the propagation and ampliﬁcation of existing background harmonics caused by system resonances at lower frequencies. This paper presents a method to explain harmonic propagation in transmission systems using theory of standing waves. By analyzing the characteristics of standing waves, one can ﬁnd useful information on the nature and extent of harmonic propagation. It is shown that the presented method describes the standing harmonic voltage on power system lines with different topology and therefore, how a simple change in grid conﬁguration affects the harmonic propagation. The results are veriﬁed with good agreement against simulations with commercial power system software. Guides to advance this study to meshed systems is given


I. INTRODUCTION
Electricity systems are undergoing significant changes in many countries in recent years due to the integration of renewable energy sources.At transmission level, this typically leads to far-reaching structural changes to the grid, and in this context these changes and expansions are, among other things, limited by the public's opposition to both substations and overhead line installations.To counter the latter, underground cables are to a large extent used in combination with overhead lines or completely as replacement for these.
One of the greater challenges of using underground cables (UGC) instead of overhead lines (OHL) is a considerable increase in harmonic distortions in the 50-2500 Hz (or 60-3000 Hz) frequency range.A practical example from Denmark is given in [1], where significantly increased harmonics are measured in two geographically remote substations following commissioning of an 8 km long 400 kV cable circuit.Here, Energinet, the Danish Transmission System Operator (TSO), commissioned two parallel 400 kV UGC systems in the Vejle-Aadal area to partially replace a single overhead line system.Fig. 1 shows a segment of the existing transmission grid in Denmark as of 2020.The location of UGC is shown in the zoom in Fig. 1.Increased 11th and 13th harmonic voltages in two substations Trige (TRI) and Fraugde (FGD) were measured after commissioning, 90 km and 80 km electrical length, respectively, from the new UGC.Measurements of the 11th harmonic voltages as 10-minute mean values are shown in Fig. 2, as documented by [1].The increased harmonic distortion Figure 2. "Measurements of 11th harmonic as 10-minute averages values at the Fraugde and Trige 400 kV substations from week 10, 2017, until week 19, 2018.The date and time stamp indicates commissioning of the 8 km underground cable section, leading to increased 11th harmonic voltage" [1].
has been linked to harmonic emissions from Konti-Skan and Great Belt LCC HVDC converter stations located in Vester 978-1-6654-1639-9/22/ $31.00 ©2022 IEEE Hassing (VHA) and in FGD, respectively.It is evident that a relatively small and simple change in the existing transmission system, the replacement of 8 km of OHL with UGC, resulted in major changes in harmonic propagation.Energinet have developed a new simulation model for better estimation of the harmonic distortions [2], [3].This new method and the model have been validated against measurements of harmonic voltages at specific nodes in the Danish transmission system.However, the observations and the work in [1]- [3] still leave us with basic questions such as: • How does a meshed grid generally behave with regard to harmonic propagation and amplification?• How far can a harmonic distortion spread in a meshed system?• How does changes in the grid affect harmonic propagation and distortion?Answers to the above questions will greatly help to better understand harmonic propagation in a meshed system.There are several paths a harmonic current from a given source can flow to a given node and thereby cause a harmonic distortion at this node and also along the entire path.In addition, the presence of several sources will add up as a vector sum of harmonics.It is also a difficult challenge that the number of harmonic sources today and in the future is unknown and, consequently, it is difficult to compute in the design state the development in harmonic distortions.
As there are several paths from the sources, the topology and the individual line types will influence how harmonic currents is distributed.There is a need for a method for determining and, more important, understand the propagation of harmonics with focus on the impact of the individual transmission lines.Based on this knowledge, we will be able to forecast from which nodes a given new source is most likely to result in the greatest harmonic distortion, the impact of new or changed lines, or even the impact of an overall restructuring and strengthening of the transmission system on the propagation.Once this is known, the TSO will to a greater extent be able to understand the propagation and choose an advantageous location for mitigating equipment.
It has been indicated by [4], that the differences in harmonic magnitudes along a line are due to standing wave effects and that shifting of the resonant frequencies is caused by line terminations.This is a well known phenomenon within communication and microwave theory, but no research has been found in this field related to power system harmonics.
This paper presents observations of correlations between harmonic propagation in cases where a system is changed topologically, in the form of replacement of an OHL with an UGC or vice versa.The main contribution of this paper is the use of the theory of standing waves to explain the questions posed above.The complexity of this subject can however, not be described in this paper alone so the focus is on the effect of grid changes on the propagation and distortion.This type of explanation to harmonic propagation cannot be achieved using today's load flow based harmonic studies which only provides information in nodes and no information on the impact of harmonics along lines nor any information on the propagation of selected harmonics into specific nodes.Furthermore, the proposed method helps in understanding the phenomena of harmonic amplification.
In future work, the results will be used to develop a harmonic propagation model for meshed systems, which will be able to determine the distortions scattered throughout the system as seen from individual sources and topological changes in the system.This paper is organized as follows.Characteristics of harmonic impedance for OHL/UGC lines are discussed in section II.The concept of standings waves is introduced in section III.A developed model are presented in section IV.In section V, case studies are presented to demonstrate the unique information that can be provided by the theory of standing waves.The obtained results are discussed in section VI along with discussions on how to advance this study in regards to meshed systems.Lastly, in section VII, the conclusions are summarized.

II. SCREENING STUDY OF HARMONIC IMPEDANCE
Electrical signals and the characteristics of network components can be represented either in the time domain or the frequency domain.Harmonic studies can be carried out in one of these domains or their hybrid combination, depending on study requirements [5] [6].For harmonic studies related to a meshed transmission system, harmonic load flow calculations and frequency sweep calculations are typically used.Both types of studies are limited to presenting the results in the form of voltage levels and impedance characteristics at specific nodes (or substations).
A research project DANPAC 2020 (DANish Power system with underground Ac Cables) has been initiated by Energinet in order to gain more knowledge about the increased use of UGC in the meshed transmission system, as the overall goal is to be able to use a significant increased amount of UGC in the Danish transmission system.One of the first studies in this research project naturally involved studying impedance characteristics for different line configurations, in this case nine study cases has been investigated.A simple model consisting of two lines, one line from node "A" to "B" and another line from node "B" to "C" has been used as shown in Fig. From this, we observed that the location of parallel resonances in the frequency spectrum on Fig. 4 highly depends on the line type and their termination.If we observe the differences between study case 3 (OHL followed by UGC) and study case 4 (UGC followed by OHL), it can be seen that the parallel resonance frequencies and impedance is highly dependent on the topology.For study case 3, the impedance peaks are approximately 5 times higher than for study case 4, indicating that UGC has a damping effect on the impedance.It is also evident, that the first resonance frequency appears at lower frequency (≈ 250 Hz) when the first line section is an UGC compared to an OHL section (≈ 500 Hz).This phenomenon is of course not a new discovery, but merely the recognition that this physical behaviour is not fully known today.The focus will be on the presented study cases 1, 3, 4 and 7 in Fig. 4.

III. STANDING WAVES
A mathematical model of a uniform transmission line could be a model with series impedance and shunt admittance uniformly distributed along the line, such that the observed line has the same specific line constants per unit of length.A theoretically correct treatment of the transmission conditions for such a line can only be carried out by solving the differential equations of voltage and current.The solution for determining currents and voltages everywhere on a line with arbitrary constants, the line being applied to a sinusoidal voltage with a fixed frequency, can be obtained using the theory of travelling waves [7].This theory is well described in [7]- [9] or similar.
Any current and voltage distribution on an electrical line can be interpreted as the resultant of two travelling waves propagating in opposite directions on the line and each of them satisfies simple laws, that for uniform lines the propagation UGC UGC occurs at a constant velocity and that the relationship between current and voltage for the individual waves is constant, designated as the characteristic impedance Z 0 (at a given frequency).The incident and reflected waveforms will superimpose and it is said that stationary waveforms or standing waves are present [7]- [9].Lets look at an illustrative example in Fig. 5. Here, four standing voltage waves on a 90 km OHL (study case 1) are shown at the frequencies 700 Hz, 2100 Hz, 3500 Hz and 4800 Hz, respectively.As it can be seen, more wave peaks and valleys appear on the line as the frequency increases.As it is evident that harmonics are simply standing waves on a line, then the wave theory can describe how a given harmonic voltage (and current) propagates along one or more interconnected lines.If the wave theory can describe the propagation, it must also necessarily be able to explain max and min values at selected points along a line or at the termination point of the line.

IV. MODEL DESCRIPTION
We define the line topology used for the four cases as given in Table I, which are dealt with in this paper.The line primary constants are given in the Appendix and the characteristic impedance calculated at 1000 Hz for OHL and UGC, respectively, is as given in Table II.A schematic diagram of the simple model used is shown in Fig. 6.The points marked as "A", "B" and "C" are nodes (i.e. the points where the characteristic impedance may change).Line 1 is between node "A" and "B" and line 2 is between node "B" and "C".Γ coefficients, Z 0x is the complex characteristic impedance of the line, γ x is the complex propagation constant of the line and d x is the line length.For all four cases, d 1 = d 2 = 45 km.Z G is the source impedance, and Z L is the load.A sinusoidal voltage V G at a frequency of 1000 Hz is applied, leading to an incident voltage V i of 1∠0 pu.The objective of this model is to determine the standing voltage waves on lines "A-B" and "B-C" when an incident voltage wave of 1 pu is sent in at node "A".In all cases, Z G will be effectively matched with Z 01 , which leads to no reflection at "A".

A. Transmission and Reflection Coefficients
The transmission of the voltage wave at some interface is quantified trough the transmission coefficient for voltage, defined as the ratio between the transmitted voltage amplitude V t→ and the incident voltage amplitude V i→ , and is given as [8]: As a definition in this paper, the arrow indicates the direction of propagation of the wave.In this context, Z 02 is the impedance in which the incident wave is transmitted, and Z 01 is the impedance from which the incident wave comes.As regards reflection coefficients, these can be determined as [8]: In ( 2), the designations for the reflection of the incident wave from node "A" to node "B" are shown.In the event that Z 01 = Z 02 , Γ B will have a value of 0, i.e. there is no reflection of the incident wave, but this is absorbed completely by the subsequent medium.By using (1) and ( 2), the transmission and reflection coefficients for the four cases can be determined as shown in Table III.Neither the coefficients Γ A← nor Γ A→ were calculated, since it is not expected that waves are reflected from node "A", as the source Z G is effectively matched with Z 01 in all cases.
It can be seen that for the open-end termination, the numerical value of the reflection coefficient at "C" is 1 and the transmission coefficient is 2. For the short circuited termination, the numerical value of the reflection coefficient at "C" is -1 and the transmission coefficient is 0. For case 3 and 4, reflections occur at "B", which is not the case for 1 and 7. When the coefficients have been calculated, they can be used to determine the magnitude of the reflected voltage wave V r← and the transmitted voltage wave V t→ at a node where an impedance mismatch is present.For transient state, this should thus mean that when a continuous wave front with amplitude V i→ arrives at node "B" from the left, there will be a partial reflection and the value of the reflected voltage is V i→ Γ B← .At the same time, there will also be a partial transmission of the voltage wave, with a similar amplitude V i→ T B→ This wave, which arrives at node "C", will be partially reflected back to node "B", partially transmitted towards "A" and also reflected back towards "C".The transmitted part is V i→ T B→ Γ C← T B← e 2γx•dx , where the exponential expression describes the phase shift of the voltage wave, which goes from "B" to "C" and then returns to "B".The propagation constant γ x is in this case related to "B-C".This is illustrated in Fig. 7.For standing waves (or stationary waves) dealt with in this paper, only Γ B← and T B← at node "B" are used for determining the standing wave on line "A-B".

B. Standing Wave Model
The previous section showed that when an incident wave on a transmission line reaches the end of the line, the wave is either fully absorbed by a matched load, reflected all the way back to the sender's end, or it is partially reflected, while some of the incoming wave is transmitted to the subsequent medium.The latter is the case at any time where an impedance mismatch is present between the transmission line on one side of a node and the transmission line or load at the other side [8], [9].According to [8], the standing voltage wave V at a distance x along a line can generally be determined as: where θ is the angle of the reflection coefficient Γ, V i is the incident voltage and γ is the propagation constant.Similarly, the current can be determined as [8]: In matrix form, ( 3) and ( 4) may be expressed as: We define [Γ E ] as: for x is an equivalent reflection coefficient determined as described in section IV-A.Only the voltages are considered in this paper.In relation to the study cases discussed, the voltage along lines "A-B" and "B-C" has been calculated by using (5) with two different loads:

C. Verification of the Model
It has not been possible to verify the method against real life measurements.The model, on the other hand, has been implemented in a commercial power system simulation software with a view to comparing the results.There is currently a general perception that frequency domain simulations in power system analysis software typically produce valid results, and this method has therefore been chosen to verify the presented model.In this simulation software, the two lines "A-B" and "B-C" from Fig. 6 with a total length of 90 km are modelled using the line parameters in Appendix.The impedance and admittance are calculated simply as Z = R + jωL and Y = G + jωC, respectively, which is also the case for the wave model and hence, the results are comparable.The two 45 km lines are modelled as a series of lumped parameter line sections, each of a length of 3 km.It is known that distributed models are more accurate than a single lumped parameter model.However, according to [5], cascading sections improves the accuracy of the lumped parameter model up to a certain frequency, depending on the circuit length and the number of sections.In this case, 15 cascaded sections are used per line and it is shown in [5], that 10 cascaded sections of 100 km UGC gives sufficient results for at least 1000 Hz.For the simple verification done here, that is, verification of voltage amplitude along the line and the corresponding wave shape, this is assumed to be a sufficient approach.A 1 pu 1000 Hz voltage source is applied and harmonic load flow were carried out.The results are compared with results obtained with (5).

A. Open-end Termination
The voltage along the line at open-end (OC) termination for the four cases are shown in Fig. 8.For case 1 and case 7, Z 01 = Z 02 and the effective reflection coefficient at node "B" is 0. The standing voltage wave is in this case determined by the load.For cases 3 and 4, where Z 01 is different from Z 02 , the standing voltage wave on line 1 is determined by the effective reflection coefficient at node "B" and the standing voltage wave on line 2 is determined by the effective reflection coefficient at node "C".It can be seen that the voltage at the endpoint (i.e.x=90 km) for all four cases shows a maximum voltage amplitude, which is as expected when terminated with an infinite load [8], [9].For case 1 and case 7, where the reflection coefficient at node "B" is 0, it can be seen that the voltage curve is smooth in the transition between line 1 and line 2 (at x=45 km), whereas cases 3 and 4 show a transition between the standing wave on line 1 and the standing wave on line 2 and that these two waves look different.Specific for case 3, the wave voltage and the simulation voltage deviates at the first few kilometres.It can also be seen that the voltage along the lines at specific distances moves towards zero.The distance from the end point (x=90 km or x=45 km) to the point on the lines where the voltage moves towards zero coincides with 1/4 wave length for the line type in question as shown in Table II.That is approx 21 km for UGC and approx 71 km for OHL at 1000 Hz.

B. Short Circuited Termination
The results for termination with a short circuit (SC) appear from Fig. 9.It can be seen that the voltage at the endpoint (i.e.x=90 km) for all cases shows a minimum voltage amplitude, which is as expected for termination with a short circuit [8], [9].The transition behaviour between line 1 and line 2 for all cases are similar to those observed in Fig. 8 for an infinite load.It can be seen that the voltage along the lines at specific distances moves towards zero.In case of short circuit termination, this equals two times 1/4 wave length as can be seen from case 3 and case 7 (approximately 42 km).
For case 1, the minimum value can not be observed, as this would equal approximately 142 km.For case 4, the distance is approximately 21 km from node B, because of non-zero reflection at this node.

VI. DISCUSSION
The results in section V show that harmonic voltage along a line can be described with wave theory and that the results of using equations for standing waves are largely identical to harmonic load flow calculations with commercial power system software.This suggests that harmonic propagation can be considered as standing waves, as was indicated by [4].However, minor deviations in amplitude between the wave voltages and simulated voltages are seen and especially for case 3 with open-end termination, the voltages deviates at the first few kilometres, where the wave voltage increases whereas the simulation voltage decreases.
The results also show that the shape of the standing waves depends on the medium carrying the wave.It is thus possible to recognize the type of media from the waveform.This explains why the impedance characteristics differ for e.g.case 3 and case 4, as was illustrated in Fig. 4, i.e. a line with different wave length (at a specific frequency) results in another voltage at a specific point at the same distance x from the load.
With this method, it is possible to determine the propagation of the standing wave along lines with impedance mismatch at nodes, using the reflection and transmission coefficients as well as the propagation constant.It is also possible with this method to derive the expected location of impedance resonance at a given frequency based on the wavelength and the maximum and minimum values for voltage and current standing waves.This can be easily seen on Fig. 5, where the distance to a voltage peak or valley on a line can be observed.If a voltage peaks at a specific node, the current will be at a minimum and hence, a parallel resonance are present.
The method presented here may be advanced to meshed systems.Using the reflection and transmission coefficients, a "reaction matrix" can be constructed and by the use of ( 5), a better understanding of how system changes impact harmonic propagation can be achieved and used as a screening method to identify critical spots in the transmission system.

VII. CONCLUSION
In this paper a standing wave model for determining standing harmonic voltages on power system lines with impedance mismatch at a node was constructed.Methods for calculating equivalent reflection coefficients for use in the standing wave model were presented.The model was verified against simulations from commercial power system software and good agreement was found.It was shown that the presented method describes the standing harmonic voltage on lines in systems with different line types and therefore, how a simple change in grid configuration affects the harmonic propagation.

Figure 1 .
Figure 1.Segment of the West Danish transmission grid.

3 .Figure 3 .Figure 4 .
Figure 3. Simple model for screening study.theline is short circuited.Distributed OHL/UGC models with data from the Danish 400 kV system are used.The positivesequence impedance and especially the impedance peaks at resonance frequencies for each study case are shown in Fig.4.Study case 1 and 2 are frequency spectrums for OHL, study

Figure 5 .
Figure 5.The per-unit harmonic voltage as a function of frequency for study case 1.Four distinct voltage waves at the frequencies 700 Hz, 2100 Hz, 3500 Hz and 4800 Hz are shown.

Figure 6 .
Figure 6.The transmission line model.

2 Figure 8 .
Figure 8. Voltage along lines for the four cases with Z L = ∞.

2 Figure 9 .
Figure 9. Voltage along lines for the four cases with Z L = 0.

TABLE I .
LINE TOPOLOGY USED FOR CASES

TABLE II .
CHARACTERISTIC IMPEDANCE FOR OHL AND UGC MODELS

TABLE III .
REFLECTION/TRANSMISSION COEFFICIENTS AT "B" AND "C"