The MatLab Software Application for Electrical Engineering Simulations and Power System

1. Introduction

Transmission lines are the elements of the electric power system that connect the load to the generation joining the production facilities of energy over large geographic areas. You could say that the transmission of electricity is one of the most important contributions that engineering has offered to the modern civilization.

The distribution of current potential differences and the transfer of energy along a transmission line can be analyzed by various methods, and it is expected that all lead to the same result. In engineering problems, in general, it can not indiscriminately apply a single formula for solving a specific problem, without the full knowledge of the limitations and simplifications accepted in its derivation. It is worth mentioning that this circumstance would lead to its misuse. The so-called mathematical solutions of physical phenomena typically require simplifications and idealizations.

Therefore, there are several models that represent the transmission lines and can be classified as to the nature of their parameters in the model parameters constant and variable parameters the models with frequency.

The parameters in the models in terms of frequency are easy to use, but can not adequately represent the line in the entire range of frequencies at which these phenomena are transient in nature. In most cases, these designs increase the amplitude of the high order harmonics, distorting the waveform peaks and producing exaggerated errors.

For adequate representation of the transmission line should be considered that the longitudinal line parameters are strongly frequency dependent, including the variable parameters in the models with frequency, the sum of the effect of soil, developed by Carson and Pollaczek, with skin effect, whose behavior as a function of frequency can be calculated using formulas derived from the Bessel equations.

Models with variable parameters in terms of frequency are considered more accurate when compared to models that consider the constant parameters. The variation depends upon the frequency. This variation may be represented by series and parallel combination of resistive and inductive elements pure.

As transmission lines are inserted in an electrical system that has many nonlinear elements and, thus, it is difficult to represent them in the frequency domain, there is a preference for line models that are developed directly in the time domain.

Another factor, that makes the model lines developed directly in the time domain are most commonly used, is that the majority of programs that perform simulations electromagnetic transients in electrical systems require that the system components are represented in time domain.

Probably, the first model to represent the transmission line directly in the time domain is designed for H. W. Dommel, which was based on the method of the characteristics or Bergeron’s method. The model has combined the method of characteristics with the trapezoidal method numerical integration, resulting in an algorithm able to simulate transient electromagnetic networks whose parameters are distributed or discrete. This algorithm has undergone successive developments and is now known as the Electromagnetic Transients Program, or EMTP simply.

In situations where one wants to simulate the propagation of electromagnetic waves resulting from operations carried out maneuvers and switching the transmission lines, it can be represented the same as a cascade of π circuits.

In this model, each segment consists of a combination of series and parallel circuits composed by resistors and inductors. This results in an equivalent resistance and an inductance that depend on the frequency. It is considered in hese equivalent elements the ground and skin effects.

Due to the fact that EMTP-type programs are not easy to use, several authors suggest describing the currents and voltages in the cascade of π circuits by means of state variables. The state equations are then transformed into differential equations and can be solved using any computer language or mathematical software.

The representation of the line by means of state variables can be used in teaching basic concepts of wave propagation in transmission lines, the analysis of the distribution of currents and voltages along the line, and the simulation of electromagnetic transients in transmission lines with non-linear elements.

Although the technique of state variables is widely used in the representation of transmission lines, it can be seen in recent publications, that it was only used to represent the parameters of which longitudinal rows can be considered constant and independent of frequency.

However, it is recognized today that the use of constant parameters to represent the line in the whole frequency range, present in the signals during the occurrence of disturbances in it, may result in responses in which the high frequency harmonic components have larger amplitudes than the real ones.

Thus, in this chapter, it is intended, using the MatLab^TM, to enter the effect of frequency on a line represented by π circuits connected in cascade and obtain the currents and voltages on the line from the use of the technique of state variables. The method is applied in a single-phase line, which considers the presence of soil and peel effect.

This line is approximated by a cascade of π circuits which will then be represented by state equations. The state equations, which are the voltages and currents along the line, will then be simulated in MatLab^TM. The cascade will also be implemented in software such as EMTP, used for simulations of electromagnetic transients in power systems. Then the results obtained with Matlab^TM and EMTP are compared, considering the single-phase line.

Besides this mentioned modeling, for MatLab^TM software users, the main contribution of this chapter is to demonstrate the application of this program for analyzing and simulating transient phenomena in power systems. It is important for undergraduate students for understanding the basic concepts of wave propagation and transmission lines. On the other hand, specific programs, like EMTP type programs, use a specific numeric method for solving the linear systems through numeric integration. With MatLab^TM basic tools, the graduate students can apply different numeric methods for solving the same problem, comparing the results and analyzing the accuracy, precision and computational time related to each method.

2. Trapezoidal integration

It is a brief equation that relates the numerical method for lumped parameters used by the ATP. Is an arbitrary differential equation given by:

$$y(t)=k\frac{dx(t)}{dt}$$E1

The solution between the instants t₀ and t is obtained by the following:

$$dx(t)=\frac{1}{k}y(t)dt\Rightarrow {\displaystyle \underset{x_0}{\overset{x}{\int }}x(t)}dx=\frac{1}{k}{\displaystyle \underset{t_0}{\overset{t}{\int }}y(t)}dt\Rightarrow x(t)-x_0=\frac{1}{k}{\displaystyle \underset{t_0}{\overset{t}{\int }}y(t)}dt$$E2

or finally:

$$x(t)=x_0+\frac{1}{k}{\displaystyle {\int }_{t_0}^{t}y(t)dt}$$E3

The numerical solution of the equation is obtained at discrete instants of time, given:

$${\displaystyle {\int }_{x(t-\Delta t)}^{x(t)}xdx=\frac{1}{k}}dx(t)=\frac{1}{k}{\displaystyle {\int }_{t-\Delta t}^{t}y(t)dt}\Rightarrow x(t)=x(t-\Delta t)+\frac{1}{k}{\displaystyle {\int }_{t-\Delta t}^{t}y(t)dt}$$E4

(4)

The solution of equation (2.4) consists of numerical integration of a discrete function y(t) which must be performed in discrete steps size Δt, which solution is the area under the curve, limited by the instants t and t-Δt as shown in Figure 1.

Figure 1.

(a) function y(t) slight to be integrated; (b) trapezoidal rule integration.

Solving equation (2.5) by the trapezoidal rule integration, which allows a linear variation between y (t-Δt) and y (t) of the integrated, as shown in figure 1 (b), as:

$$x(t)=x(t-\Delta t)+\frac{1}{k}{\displaystyle {\int }_{t+\Delta t}^{t}y(t)dt=x(t-\Delta t)+\frac{1}{k}}\frac{y(t)+y(t-\Delta t)}{2}\Delta t$$E5

or:

$$x(t)=\frac{\Delta t}{2k}y(t)+\left[x(t-\Delta t)+\frac{\Delta t}{2k}y(t-\Delta t)\right]$$E6

in other words:

$$x(t)=C_{TR}y(t)+h_{TR}\left(t-\Delta t\right)$$E7

where:

$$\{\begin{array}{l}{C}_{TR}=\frac{\Delta t}{2k}\\ h_{TR}(t-\Delta t)=x(t-\Delta t)+\frac{\Delta t}{2k}y(t-\Delta t)\end{array}$$E8

It may be noted that the term historical h_TR(t-Δ t) depends on x (t) and y (t) after t seconds, t, is already known during the solution.

3. Modeling cascada of π circuits

Discussion will now represent a mathematical model for a transmission line using an electrical circuit. With this model, you can make a study of the behavior of a transmission line for energizing and switching transient simulations. A transmission line, whose parameters can be considered independent of frequency, can be represented in an approximate manner and following a series of restrictions, as a cascade of π circuits.

Each segment of π circuit consists of one resistor and one series inductance in series. Completing the circuit, there are components in parallel: capacitance and conductance. It is shown in Figure 2.

Figure 2.

A π circuit.

It’s represented a transmission line through this model, connects to n π circuits in cascade. So Figure 3 shows a single-phase transmission line length d represented by n π circuits connected in cascade.

Figure 3.

Line represented by a cascade of π circuits.

In Figure 3, the parameters R and L are the resistance and inductance of the longitudinal line, respectively. The parameters G and C are the capacitance and conductance of transverse dispersion, respectively. These parameters are written as:

$$R=R\text{'}\frac{d}{n}$$E9

$$L=L\text{'}\frac{d}{n}$$E10

$$G=G\text{'}\frac{d}{n}$$E11

$$C=C\text{'}\frac{d}{n}$$E12

In the equations (3.1) to (3.4), R 'and L' are the resistance and inductance of the longitudinal row per unit length, respectively while the terms G 'and C' are the capacitance and conductance per unit cross-sectional line in length.

Using this representation of the line, a state model is formulated for the energy system that uses the voltages on the capacitors and the currents through the inductors as state variables. The system that describes the state equations is transformed into a set of differential equations whose solution is given by the use of trapezoidal integration. The state variables are found by solving the set of equations.

Although the technique of state variables is widely used in the representation of transmission lines, it is applied only to depictions of longitudinal lines whose parameters can be considered constant and independent of frequency.

However, it is recognized today that the use of constant parameters to represent the line in the whole frequency range, present in the signals during the occurrence of disturbances in it, may result in responses in which the high frequency harmonic components have larger amplitudes than are actually.

4. Line representation with parameters of frequency dependent

The representation of transmission lines through cascades of π circuits, taking into account the effect of frequency, is usually implemented in EMTP-type programs.

Another drawback of the type EMTP programs is that they limit the amount of π circuit that can be used to represent the line. Thus, depending on the length of the line to be depicted, the quality of results obtained from simulations may be compromised.

To circumvent the difficulties mentioned, it is suggested to describe the cascade of π circuits by means of state equations. However, some authors disregarded the effect of frequency on the parameters of the longitudinal line.

Figure 4.

Circuit that represents the longitudinal parameters of the line.

The models proposed by would become more complete if the effect of frequency on the longitudinal parameters of the line was inserted in them.

So, the parameters of a transmission line can be synthesized by means of a circuit of the type shown in Figure 4.

It’s used a cascade of π circuits to represent a transmission line taking into account the effect of frequency on the longitudinal parameters. In this case, every π circuit will be shown in Figure 5.

Figure 5.

Cascade of π circuits considering the effect of frequency.

In Figure 5, the RL parallel associations are as many as necessary to represent the variation of parameters in each decade of frequency that will be considered. First, arrays are displayed in a state for a line represented by a single π circuit, whereas the effect of frequency is synthesized by n RL associations. Then, the results will be extended to a line represented by a cascade of n π circuit, considering n RL associations to synthesize the effect of frequency.

Before being certain state equations for a line represented by a cascade of n π circuits, it will be shown in detail the development of equations of state considering only one π circuit.

Then, the development done for a single π circuit element can be extended to a generic cascade with any quantity of these circuits.

Considering, as shown in the Figure 5, a transmission line represented by a single π circuit, the effect of frequency on the longitudinal parameters is represented by n RL associations.

In line shown in Figure 5, the voltages at terminals A and B are u(t) and v_k(t), respectively. It is also considering that the currents i_k0, i_k1,..., i_km are circulating through the inductors L₀, L₁, L₂,..., L_m, respectively. These currents are dependent time functions and the notations do not include the time dependence for more simplicity. This simplification is also applied to the v_k(t) state variable. From the currents and voltages in the circuit of Figure 5 it can be determined:

$$\frac{d{i}_{k0}}{dt}=\frac{i_{k0}}{L_0}\left(-{\displaystyle \sum _{j=1}^m{R}_j}\right)+\frac{1}{L_{L0}}\left({\displaystyle \sum _{j=1}^m{R}_j{i}_{kj}}\right)+\frac{1}{L_0}u(t)-\frac{1}{L_0}{v}_k$$E13

$$\frac{d{i}_{k1}}{dt}=\frac{R_1}{L_1}{i}_{k0}-\frac{R_1}{L_1}{i}_{k1}$$E14

$$\frac{d{i}_{k2}}{dt}=\frac{R_2}{L_2}{i}_{k0}-\frac{R_2}{L_2}{i}_{k2}$$E15

$$\frac{d{i}_{km}}{dt}=\frac{R_m}{L_m}{i}_{k0}-\frac{R_m}{L_m}{i}_{km}$$E16

$$\frac{d{v}_k(t)}{dt}=\frac{2}{C}{i}_{k0}-\frac{G}{C}{v}_k(t)$$E17

The equations (4.1) to (4.5), describing the circuit shown in figure 4.2, can be written as:

$$\stackrel{•}{x}=Ax+Bu$$E18

For one circuit, the A matrix is substituted by:

$$A_{\pi }=\left[\begin{array}{cccccc}-\frac{{\displaystyle \sum _{j=0}^{j=m}{R}_j}}{L_0}& \frac{R_1}{L_0}& \frac{R_2}{L_0}& \cdots & \frac{R_m}{L_0}& -\frac{1}{L_0}\\ \frac{R_1}{L_0}& -\frac{R_1}{L_0}& 0& \cdots & 0& 0\\ \frac{R_2}{L_2}& 0& -\frac{R_2}{L_2}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & 0& 0\\ \frac{R_m}{L_m}& 0& 0& \cdots & -\frac{R_m}{L_m}& 0\\ \frac{2}{C}& 0& 0& \cdots & 0& -\frac{G}{C}\end{array}\right]$$E19

$$B^T=\left[\begin{array}{cccc}\frac{1}{L_0}& 0& 0& \begin{array}{cc}\dots & \begin{array}{cc}0& 0\end{array}\end{array}\end{array}\right]$$E20

$$x_k^T=\left[\begin{array}{cccccc}{i}_{k0}& i_{k1}& i_{k2}& \cdots & i_{km}& v_k\end{array}\right]$$E21

$$\stackrel{•}{x_k}=\frac{d{x}_k}{dt}={\left[\begin{array}{cccccc}\frac{d{i}_{k0}}{dt}& \frac{d{i}_{k1}}{dt}& \frac{d{i}_{k2}}{dt}& \cdots & \frac{d{i}_{km}}{dt}& \frac{d{v}_k}{dt}\end{array}\right]}^T$$E22

In the equations (4.8) and (4.9), B^T e x_k^T correspond to transposed B and x_k, respectively.

The obtained results show that the vesctor x_k have (m + 2) elements and the matrix A is a (m + 2) order square matrix.

Based on the equations and the results for one π circuit, it can be extended the analysis to a cascade of π circuits. Thus, the matrix A will have an order of n(m +2) and the vector x has dimension n(m +2). The A matrix can be written as:

$$A=\left[\begin{array}{ccccccc}{A}_{11}& A_{12}& 0& \cdots & 0& 0& 0\\ A_{21}& A_{22}& A_{23}& \ddots & \ddots & \ddots & 0\\ 0& A_{32}& A_{33}& \ddots & 0& \ddots & 0\\ 0& 0& A_{43}& \ddots & A_{\left(n-3\right)\left(n-2\right)}& 0& 0\\ 0& \ddots & 0& \ddots & A_{\left(n-2\right)\left(n-2\right)}& A_{\left(n-2\right)\left(n-1\right)}& 0\\ 0& \ddots & \ddots & \ddots & A_{\left(n-1\right)\left(n-2\right)}& A_{\left(n-1\right)\left(n-1\right)}& A_{\left(n-1\right)n}\\ 0& 0& 0& \cdots & 0& A_{n\left(n-1\right)}& A_{nn}\end{array}\right]$$E23

The elements x₁, x₂, …, x_n are describe by equation (4.9). In equation (4.11), A is a tridiagonal matrix which elements are square matrices of order (m +2). In this case, a generic element A_KK at main diagonal of the matrix A is written as:

$$A_{kk}=A_{\pi }$$E24

The A _π matrix is defined in equation (4.7).

An element of any upper subdiagonal in equation (4.11) is a square matrix of order (m +2) which the only one nonzero element is located in the first column of last row and has the value$$\left(\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$C$}\right.\right)$$.

The structure is:

$$\begin{array}{l}{A}_{ik}=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ ⋮& \ddots & \cdots & ⋮\\ 0& \cdots & \ddots & ⋮\\ -\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$C$}\right.& 0& \cdots & 0\end{array}\right]\\ i=k-1,2\le k\le n\end{array}$$E25

The subdiagonal elements in the equation (4.11) are square matrices of order (m +2). These arrays have a single nonzero element which is in the last column of first row. It is a value of$$\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$L_0$}\right.\right)$$. The structure is:

$$\begin{array}{l}{A}_{ik}=\left[\begin{array}{cccc}0& \cdots & 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$L_0$}\right.\\ ⋮& \ddots & \cdots & 0\\ ⋮& \cdots & \ddots & ⋮\\ 0& 0& \cdots & 0\end{array}\right]\\ i=k+1,1\le k\le n-1\end{array}$$E26

Considering a cascade of π circuits, the vector B has the dimension of n(m +2) and if it is connected a u(t) source at the beginning of the line, the vector B has a single nonzero element, which is the first array element and it has the value$$\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$L_0$}\right.\right)$$.

The state equation, that describes a line representation by a cascade of π circuits cant be solved by numerical methods, like Euler and Heun methods. Other numeric methods are described in the next item

5. Other numeric methods

Besides the trapezoidal integration, there are also other numeric methods that can be used for solving state equations related to the modeling of transient phenomena in transmission lines. Among them, it may be mentioned, for example, the Simpson’s and the Runge-Kutta’s method. The following items will describe these methods. In this case, the MatLab^TM is important because it facilitates the comparisons among the mentioned methods. With this software, undergraduate students can analyze the application of different numeric methods for solving an engineering problem, besides the analysis about wave propagation and the transmission line modeling. Graduated students can analyze what it is the best option for specific transmission line characteristics and develop the numeric method related to the simulation of transient phenomena in power systems.

5.1. Simpson’s rule

Simpson’s rule is based on the assumption that, given short intervals of time the derivative of the function to be integrated, the function (y') can be approximated by a second degree function as shown in Figure 6.

From Figure 1, it is obtained:

$$y\left(t_{k+1}\right)=\frac{1}{3}\left[y\text{'}\left(t_{k+1}\right)+4y\text{'}\left(t_M\right)+y\text{'}\left(t_k\right)\right]\Delta t$$E27

From (5.1) and (4.6), it is obtained:

$$x\left(t_{k+1}\right)=a_1\left(a_2+a_3+a_4\right)$$E28

Equation (5.2) is applied Simpson's rule for solving the state equation of state of the type shown in (5.1). In this case, the matrices a₁, a₂, a₃ and a₄ are written as:

$$a_1=\left(I-\frac{\Delta t}{6}A\right)$$E29

$$a_2=\left(I+\frac{\Delta t}{6}A\right)x\left(t_k\right)$$E30

$$a_3=\frac{2\Delta t}{3}Ay\left(t_M\right)$$E31

$$a_4=\frac{\Delta t}{6}B\left[u\left(t_k\right)+4u\left(t_M\right)+u\left(t_{k+1}\right)\right]$$E32

It is defined the following terms:

$$\Delta t=t_{k+1}-t_k$$E33

$$t_M=\frac{t_{k+1}+t_k}{2}$$E34

Figure 6.

Approximation of the derivative of the function to be integrated (solid curve) by a function of the second degree (dashed curve).

5.2. Runge-Kutta’s rule

Runge-Kutta’s rule is defined by the following:

$$x\left(t_{k+1}\right)=x\left(t_k\right)+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)$$E35

The terms in the last equation are described by:

$$k_1=\Delta t\left[Ax\left(t_k\right)+Bu\left(t_k\right)\right]$$E36

$$k_2=\Delta t\left\{A\left[x\left(t_k\right)+\frac{k_1}{2}\right]+Bu\left(t_k\right)\right\}$$E37

$$k_3=\Delta t\left\{A\left[x\left(t_k\right)+\frac{k_2}{2}\right]+Bu\left(t_k\right)\right\}$$E38

$$k_4=\Delta t\left\{A\left[x\left(t_k\right)+k_3\right]+Bu\left(t_k\right)\right\}$$E39

5.3. Comparisons among the numeric methods

The Matlab^TM software is easily applied for the development of numeric routines. Using this characteristic, it is possible to compare different numerical methods. Figure 7 shows the results for the three methods investigated for the simulation of electromagnetic transients in a single-phase transmission line.

Figure 7.

Comparisons among the three numerical methods.

In the last figure, it is shown the result obtained with the application of a step voltage source on the initial transmission line. This is a 20 kV step voltage source. The transmission line is modeled as a mono-phase circuit and the linear system is solved using three numeric methods: trapezoidal rule, Simpson’s rule and Range-Kutta’s one.

Based on the results shown in Figure 7, the three numeric methods lead to similar results considering the method accuracy. Because of these results, it is necessary to analyze the other characteristics. An important characteristic for the application of the numeric methods for transient analysis in transmission lines is the simulation time. For example, this is important for analysis in very short simulation times. The next table shows the simulation time comparisons obtained with the MatLab^TM software.

Step of calculating (µs)	Computational effort (simulation time)
Heun	Simpson	Runge-Kutta
0,1	t1	22,52 t1	0,63 t1
0,5	t2	8,94 t2	0,23 t2
1,0	t3	8,40 t3	0,26 t3

Table 1.

Simulation time required by numerical methods

Table 1 shows the time relationship between the applied numerical methods. It is possible to carry out these comparisons because the Matlab^TM provides a simple way to replace the codes related to the numeric routine without changing the structure of the main code, maintaining the same computational flowchart. Using this characteristic, it is possible to determine the simulation time of each numeric method, because the computational time for introducing the data and other characteristics of the simulated transmission line is equal for all compared numeric methods. Because of this, the MatLab^TM is an adequate tool for developing new models and numeric routines used for analysis and simulations of transient phenomena in transmission lines.

6. Model implementation of single phase line

In this section, it is presented the implementation of a single-phase transmission line through a cascade of π circuits, considering the effect of frequency on their longitudinal parameters and using the concept of state variables. Then, the results obtained are compared with results obtained with EMTP.

6.1. Block diagram of the program for single-phase line

The model that represents the single-phase line was implemented in a microcomputer using the software MatLab^TM.

Data from single-phase transmission line are read in the first part of the program, then the parameters are calculated from the transmission line considering the influence of frequency.

After observing the behavior of single-phase line parameters as a function of frequency, it is necessary to represent this influence on the transmission line model proposed in item 4. This was done using the method called vector fitting and the longitudinal single phase line parameters are fitted by means of rational functions.

With these parameters synthesized and distributed in the proposed model of a single-phase line, it is possible to calculate the voltages and currents at the terminals of this line or at any point of the line.

The model represented by a cascade π circuits can be represented by a linear system that is represented by state variables. For the solution of the system represented by$$\stackrel{•}{x}=[A]x+[B]u$$, it’s used Heun formula (trapezoidal integration method). This method is widely used in simulations of electromagnetic transients in power systems

Figure 8 is shown a block diagram of the algorithm of the program developed for single-phase line.

Figure 8.

Block diagram of the routine developed for a single-phase.

6.2. Calculation of the parameters of the transmission line single phase

It’s possible to bring the longitudinal parameters of single-phase line by means of rational functions. Using the vector fitting method and allowing the fitting the frequency dependence of these parameters, this effect is entered in the discrete parameter model.

The equation that summarizes the parameters of longitudinal single-phase line is given by:

$$Z_{VFIT}(\omega )=R_0+j\omega L_0+\frac{j\omega R_1}{j\omega +\frac{R_1}{L_1}}+\frac{j\omega R_2}{j\omega +\frac{R_2}{L_2}}+\frac{j\omega R_3}{j\omega +\frac{R_3}{L_3}}+\frac{j\omega R_4}{j\omega +\frac{R_4}{L_4}}$$E40

The values of R and L in the equation (6.1) found by the method vector fitting are shown in Table 2.

From the values of Table 1, it can view a summary of the parameters of longitudinal line as follows: replacing the values of Table 1 in the expression (6.1) and attribute values are included in the frequency range 10¹ to 10⁶ for it is possible to calculate the longitudinal impedance synthesized.

In the equivalent circuit, the values were considered in the frequency range from 10¹ to 10⁶, because transients that occur in the transmission line are within this frequency range. So the four blocks RL in parallel, shown in Figure 9, are representing the influence of frequency on longitudinal parameters of single-phase transmission line.

Figure 9.

Circuit that represents the longitudinal parameters of the line.

Resistors (Ω/km)			Inductors (mH/km)
R’0	0,026	L’0		2,209
R’1	1,470	L’1		0,74
R’2	2,354	L’2		0,12
R’3	20,149	L’3		0,10
R’4	111,111	L’4		0,05

Table 2.

Values of elements R e L for a single-phase.

7. Transmission line with corona effect

7.1. General aspects about the corona

The corona discharge mechanism is an electrostatic phenomenon due to ionization in an insulating material, usually a gas, subject to electric field intensities above a critical level.

Electrical discharges in gases are usually triggered by an electric field that accelerates free electrons therein. When these electrons acquire enough energy from the electric field, they can produce new electrons from the collision with other atoms. It is the process of impact ionization. During its acceleration in the electric field, each free electron collides with atoms of oxygen, nitrogen and other present gases, missing, that collision, part of its kinetic energy. Occasionally, it can achieve an electron atom with sufficient force so as to excite it. Under these conditions, the atom is achieved to a higher energy state. The orbital state of one or more electrons and the electron moves colliding with the atom loses some of its energy to create this state. Subsequently, the atom can hit revert to its initial state, releasing the excess energy as heat, light, electromagnetic radiation and acoustic energy. An electron can also collide with a positive ion, converting it into neutral atom. This process, called recombination, also releases excess energy.

7.2. Corona effect representation

After the pioneering work of Peek (1915), several measurements have made on lines and experimental laboratories have examined the nature of the corona and its influence on wave propagation in transmission lines. These works have had fundamental importance, contributing to the understanding of the basic mechanism of the corona.

In 1954, Wagner et al. (1954) and Wagner and Lloyd et al. (1955) published two articles that would be a reference for future work on corona. Voltage experimental measurements were made of in the laboratory of a conductor under corona (project called Tidd 500 kV).

Through these experimental results, it has developed empirical formulas and procedures for considering the effects of attenuation and distortion in the spread of outbreaks. These formulas and procedures are based on the voltage gradient, and voltage curves are obtained from attenuation measurements and the power dissipation due to the corona effect. The usefulness of these methods is limited however, because it requires different approaches and uses abacuses.

The corona model can be divided into three classes: analog models, mathematical models and physical models. The electrical circuit analog models are designed to reproduce the geometric increase in the capacitance of the conductor to the voltage reaches critical ionization.

Like the analog models, mathematical models roughly reproduce the characteristics of drivers under the corona, but by means of mathematical equations. Several authors have empirical formulations for the variation of capacitance of the line based on functions and constants derived from measurements. Most models show a linear relationship between capacitance and dynamic tension.

Due to the complexity in describing mathematically the physical phenomena involved and the insufficient amount of data propagation in corona, it is not available generic models of this nature yet.

7.3. Gary’s model

The equations describing the corona effect is not easily implemented in the differential equations of the transmission line in order to obtain a solution formulation easyly. Thus, to obtain responses directly in the time domain, numerical models are used such as the finite difference method and the method of the characteristics. This last category of models has been developed to be implemented in EMTP-type programs. Some of these models use non-linear resistors and capacitors that are dependent on the voltage applied on them. However, most existing models of corona present satisfactory results only for a specific situation.

The mechanism of the corona effect can also be represented by the model of Gary, using a capacitance and a non-linear conductance to represent the accumulation and the pressure drops in the line. The capacitance and conductance mentioned above are variables with the applied voltage on them and are called corona capacitance (Cc) and corona conductance (Gc). Cc and Gc elements are obtained by known analytical function, and, therefore, this representation for the corona effect is called analytical model of the corona effect. This model for the corona effect can be inserted in lines represented by a cascade of circuits, where currents and voltages along line is described by state variables.

If the capacitance is represented by Gary’s corona model, it is defined as:

$$C_C=\{\begin{array}{l}C\eta {\left(\frac{v}{V_C}\right)}^{\eta -1}sev\ge V_C\\ \\ 0sev<V_C\end{array}$$E41

In equation (7.1), C_C is the corona capacitance, C is the geometric capacitance of the line segment represented by a π circuit, v is the voltage being applied to the line capacitance, C_V is the minimum voltage required to the corona effect accurrance and η is a coefficient defined as:

$$\eta =0,22r+1,2$$E42

where: r the radius of the conductor in centimeters.

Thus, given the presence of the corona effect and the effect of frequency, a differential element row can be represented as shown in Figure 10.

Figure 10.

A differential element of line considering the corona effects

The corona conductance to Gary’s model is defined as:

$$G_C=k_C{\left(\frac{1-V_C}{v}\right)}^2$$E43

Gary’s model considers that the corona effect is manifested only if the voltage V_C is greater than v and the rate of variation over time v is positive. Thus, if the corona effect is manifested at a given point P of the line, the voltage V_p at this point must satisfy the following conditions:

$$V_p>V_{C}and\frac{d{V}_p}{dt}>0$$E44

Thus, for the corona effect is present at a generic point of the line represented by a cascade of circuits, as shown in figure 10, it is necessary that the voltage at this point satisfying the two conditions shown in equation (7.4). If a condition is not met, this point will not have increased the capacitance and conductance representing corona. Therefore, the A matrix shown in equation (4.6) should be changed for each iteration as a function of cross-line voltage.

8. Testing the model

Checking the effectiveness of the developed model, it is simulated the energizing of a transmission line shown in figure 11, considering frequency independent line parameters and frequency dependent ones.

Figure 11.

Single phase line representation with opened terminal.

In Figure 11, S is a switch that be closed at time t = 0, energizing the line through a voltage source u(t). In the current procedure, the terminal B is open and the other terminal is powered by a constant voltage source. For frequency dependent line parameters, it is considered that the longitudinal parameters of the line per unit length can be perfectly summed up by a circuit consisting of four RL parallel blocks connected in series. The structure is completed using a RL series block as shown in Figure 9.

The values of R and L used to synthesize the effect of frequency on the longitudinal parameters of the line were obtained using the method proposed by [10] and are shown in Table 1. The parameters of the unit transverse line shown in figure 3 are G′=0,556 μS/km e C′ =11,11nF/km.

9. Simulation analysis of transmission lines with single phase representation

In all simulations, it is used the following values: the transmission line has 10 kilometers. It is represented through 200 π circuits. The time step used in the simulations is 50 ns and the simulation period is 600 μs. The voltage in the initial of the line is 1 kV. It can also be considered as 1 pu.

Since the values of R and L elements of the cascade of π circuits that describe the line are known, it can be obtained the state equations that describe the behavior of currents and voltages along the line. The simulations using the model proposed in this chapter were performed in MatLab^TM program, using the trapezoidal integration method. Considering the frequency independent line parameters, the circuit of Figure 4 is reduced to the R₀ and L₀ elements. In this case, the values are: R₀ = 0.05 Ω/km and L₀ = 1 mH/km.

The values of Table 9.1 are written per kilometer for the generic circuit. For the simulation, it used the values represented in the Table 9.2 that are resistance and inductance [] and [H], respectively. In the Figs 12-21, it is used the values of Table 4.

Figures 12-17 show the relationship between the number of RL parallel blocks and the inclusion of the frequency influence. In Figure 12, the resistance values are obtained using only one RL parallel block related to the high frequencies. Because of this, the resistance values for all frequency values are equal. In this case, the value is equal to the R₄ value per length unit.

Parameter	100 circuits	200 circuits
R	5 mΩ	2.5 mΩ
L	100 mH	50 mH
G	556 µS	278 µS
C	1.111 nF	555.5 nF

Table 3.

Values of line parameter used in all simulations without frequency.

Resistors(Ω)		Inductors(µH)
R0’	0.0013	L0’	110.45
R1’	0.0735	L1’	37
R2’	0.1177	L2’	6
R3’	1.0075	L3’	5
R4’	5.5556	L4’	2.5

Table 4.

Values of line parameter used in simulations with 200 circuits.

From the results of Figures 13 and 14, it is observed that the synthesis of the effect of frequency on the resistance can only be considered when inserted in the cascade of π circuits, at least, two RL parallel blocks, because each block is related to a frequency set point.

So, if more blocks are used, more frequency set points are obtained. It is confirmed through the results shown in Figures 15, 16 and 17 where the inductance values are presented. So, it is concluded that the synthesis of the longitudinal line parameters is improved with the increase of the number of RL parallel blocks.

Figure 12.

Synthesis of resistance per unit length using R’₄L’₄ parallel block.

Figure 13.

Synthesis of resistance per unit length using R’₁L’₁ and R’₄L’₄ parallel blocks.

Figure 14.

Synthesis of resistance per unit length using four RL parallel blocks.

Figure 15.

Synthesis of inductance per unit length using only R’₄L’₄ parallel block.

Figure 16.

Synthesis of inductance per unit length using R’₁L’₁ and R’₄L’₄ parallel blocks.

Figure 17.

Synthesis of inductance per unit length using four RL parallel blocks.

The result of the simulation made for the voltage input signal u(t) can be seen in Figure 18. It is shown the output voltage at the receiving end terminal of the line without the frequency influence. Using the routine without the influence of frequency from the results of Figure. 18, it is observed that there is a period of time related to the time of signal propagation through the line.

Figure 18.

Energization of the transmission line without the effect of frequency considering 200 π circuits.

Thus, it represents a time delay between input signal and output signal. After the delay, there are oscillations associated with wave reflections on the transmission line terminals that make up the output voltage shown.

Figure 19.

Voltage at the end of the transmission line with the effect of frequency considering 200 π circuits.

Using the routine with frequency influence in longitudinal parameters, it is obtained the results of Figure 20. In this figure, comparing to Figure 19, the voltage signal is attenuated because the inclusion of the frequency influence.

Figures 21 and 22 show the influence of frequency on the current results. Without the frequency influence, the obtained signal current is not attenuated and is highly modified by numeric oscillations (Figure 14). On the other hand, when the routine considers the influence of frequency (Figure 22), it is clear that the current signal could not contain those oscillations shown in Figure 21. So, those oscillations are numeric oscillations and they can be associated to the representation of the longitudinal line parameters which does not consider the frequency influence.

Figure 20.

Current at the end of the transmission line without the effect of frequency.

Figure 21.

Current at the end of the transmission considering the effect of frequency.

So, for the sequence of this work, it should be investigated what is the saturation point for the number of the π circuits and the number of the RL parallel blocks. It is carried out using the simulation results from several voltage signals that will be used as voltage sources at the initial line terminal.

For Figures 22 and 23, using the routine without the effect of frequency, it is analyzed a particular stretch of the simulation at the end of the line. It is shown the voltage values without frequency influence. In Figures. 22 and 23 it is used the values of Table 2, using 100 and 200 π circuits, respectively.

Figure 22.

Specific portion of the simulation at the end of the line with 100 π circuits and without the effect of frequency.

Based on the mentioned results, it concludes that increasing the number of π circuits leads to a decrease in the numerical oscillations.

Figure 23.

Specific portion of the simulation at the end of the line with 200 π circuits and without the effect of frequency.

In Figure 23, it is clearly a greater condensation of numerical oscillations in both x and y axes in contrast to Figure 23, where the period of the oscillations and the peak value are higher than those obtained in Figure 23. In Figures 24, 25 and 26, it is used the routine with the effect of frequency, analyzing the number of the RL parallel blocks inserted in the cascade of π circuits, searching for reducing of numerical oscillations. In Figures 24, 25 and 26, it was used the values of Tables 3 and 1, as well as, 3 and 4 RL parallel blocks, respectively.

From the results of Figures 24, 25 and 19, it is observed that the used routine including the effect of frequency in the cascade of π circuits obtains considerably fewer numerical oscillations when compared to Figures 22 and 23. From these results, it is evident that if it is increased the number of the RL parallel blocks, it is reduced the numerical oscillations and the voltage in the line end has a smoothly curve in the time domain.

Figure 24.

Voltage at the end of the line with the effect of frequency considering 200 π circuits and using R₄L₄ parallel block.

Figure 25.

Voltage at the end of the line with the effect of frequency considering 200 π circuit and using R₀L₀, R₂L₂ e R₄L₄ parallel blocks.

Figure 26.

Comparison of figures 24, 25 and 19 of the tension at the end of the line with the effect of frequency.

Concluding the analysis of the results, Figure 26 shows a comparison among the results of Figures 24, 25 and 19, leaving a clear decrease in the numerical oscillations. In this figure, it is shown a time range of the time simulation used at the last three figures. In this case, it is clearly shown that, increasing the number of RL parallel blocks, the numerical oscillations are decreased and the frequency influence is better reproduced through the state equations in mathematical software.

Considering the corona effect, Gary’s model is applied with the routine described in the previous items. It is used the line representation without frequency influence, because the representation with frequency influence has not hardly analyzed. The corona effect is introduced by

$$\begin{array}{l}{C}_C=C\cdot \eta \cdot {\left(\frac{V}{V_C}\right)}^{\eta -1}\\ \eta =0.22\cdot R_{COND}+1.2\end{array}$$E45

In this case, R_COND is the conductor phase radius in [cm], C is the transversal line capacitance and the C_C is the new value of the transversal line capacitance because the corona effect. In this chapter, the R_COND value is 2.54 cm. The simulations are carried out for some relations between V and V_C. The V value is 1 pu for all following shown simulations.

Considering Figures 27-31, the corona effect is related to the 10 π circuits in the middle line. It corresponds to the 500 m of the represented line. Figures 27-30 show results for different relative values of the corona voltage when compared to the line nominal voltage. In Figure 31, it is shown the comparisons among the results for the voltage values in the receiving end terminal. So, using a simple routine based on π circuits for transmission line representation, undergraduate students can analyze and simulate traveling wave phenomena in transmission lines.

Figure 27.

The time domain simulations with the corona effect for V_CORONA = 0,35 V.

Figure 28.

The time domain simulations with the corona effect for V_CORONA = 0,5 V.

Figure 29.

The time domain simulations with the corona effect for V_CORONA = 0,7 V.

Figure 30.

The time domain simulations with the corona effect for V_CORONA = 0,9 V.

Based on Figure 10, it’s simulated the corona effect with effect of frequency for some relations between V and V_C. The V value is 1 pu to the following shown simulation.

Figure 31.

Comparisons of the time domain simulations with the corona effect for some V_CORONA values.

Figure 32.

Comparisons of the time domain simulations with the corona effect with effect of frequency for some V_CORONA values.

10. Conclusions

It is described the application the MatLab^TM software in analysis and simulations of transient phenomena in transmission lines. Using the characteristics of this software, transmission lines are easily modeled as a mono-phase circuit. Transient simulations are also easily carried out. For these applications, it is used basic and simple tools of the MatLab^TM software. So, this software improves the analysis of the proposed problem, because it is possible to obtain several types of the graphic results that are not available in the specific programs for transient analysis like the EMTP programs. So, it is possible to analyze the resistance and inductance values that depend on the frequency when it is considered detailed transmission line models. It is possible to analyze the application of different numeric methods for solving the differential state equations by numeric integration routines. On the other hand, with a simple model of the transmission lines and the MatLab^TM software, it is possible to develop a routine that is used by undergraduate students, making easy the learning about important concepts as wave propagation, transient phenomena and transmission lines. This routine can be modified, introducing elements that are able to consider the frequency influence in the transmission line parameters. These parameters have their characteristics distributed along the line and this is considered in the mentioned routine.

The shown analysis and results can be used by undergraduate students for learning about the important concepts of power systems, transmission lines and wave propagation, for example. Related to the graduated students, it can be used for analyzing transient phenomena, developing transmission line models, improving numeric routines and comparing different numeric integration methods. So, the MatLab^TM software is an excellent tool for basic and profound studies of transient phenomena in transmission lines.