**The MatLabTM Software Application for Electrical Engineering Simulations and Power System**

Leonardo da S. Lessa, Afonso J. Prado, Rodrigo Cleber Silva, Sérgio Kurokawa, Luiz F. Bovolato and José Pissolato Filho

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/48555

## **1. Introduction**

20 Will-be-set-by-IN-TECH

*of Computational and Applied Mathematics*, Vol. 6, No. 1, pages 19-26, 1980.

*Transactions on Education*, Vol. 16, No. 3, pages 152-156, 1973.

*European Journal of Physics*, Vol. 8, No. 2, pages 88-92, 1987.

*Annalen der Physik*, Vol. 301, No. 6, pages 440-452, 1898.

2nd edition, Addison Wesley, Boston.

Press, Berkeley.

[6] Dormand J.R. & Prince P.J. (1980). A family of embedded Runge-Kutta formulae, *Journal*

[7] Feynman R.P. (2005). *The Feynman lectures on physics: the definitive and extended edition*,

[8] Hoefer W.J.R. (1973). On teaching relativity to undergraduate electrical engineers, *IEEE*

[9] Luneburg R.K. (1964). *Mathematical theory of optics*, 1st edition, University of California

[10] Petouris A., Humphries J., Russell P., Vourdas A. & Jones G.R (1994). Computer aided design for teaching modern optics to electronic engineering undergraduates, *IEE Proceedings Science, Measurement & Technology*, Vol. 141, No. 3, pages 171-176, 1994. [11] Sader U. & Jodl H.J. (1987). Computer-aided physics-particles in electromagnetic fields,

[12] Wien W. (1898). Untersuchungen über die electrische entladung in verdünnten gasen,

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.

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 435

the parameters of which longitudinal rows can be considered constant and independent of

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

Thus, in this chapter, it is intended, using the MatLabTM, 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-

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 MatLabTM. The cascade will also be implemented in software such as EMTP, used for simulations of electromagnetic transients in power systems. Then the results

Besides this mentioned modeling, for MatLabTM 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 MatLabTM 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

It is a brief equation that relates the numerical method for lumped parameters used by the

11 1 () () () () () ()

*dx t y t dt x t dx y t dt x t x y t dt kk k*

<sup>0</sup> <sup>0</sup> <sup>1</sup> ( ) ( ) *<sup>t</sup> <sup>t</sup> x t x y t dt <sup>k</sup>*

0 0 0

*x t t*

*x t t*

( ) ( ) *dx t yt k dt* (1)

0

(3)

(2)

obtained with MatlabTM and EMTP are compared, considering the single-phase line.

phase line, which considers the presence of soil and peel effect.

frequency.

to each method.

or finally:

**2. Trapezoidal integration** 

ATP. Is an arbitrary differential equation given by:

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

amplitudes than the real ones.

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 MatLabTM, 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 singlephase 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 MatLabTM. 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 MatlabTM and EMTP are compared, considering the single-phase line.

Besides this mentioned modeling, for MatLabTM 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 MatLabTM 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{d\mathbf{x}(t)}{dt} \tag{1}$$

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

$$d\mathbf{x}(t) = \frac{1}{k}y(t)dt \Longrightarrow \int\_{x\_0}^{x} \mathbf{x}(t)d\mathbf{x} = \frac{1}{k}\int\_{t\_0}^{t} y(t)dt \Rightarrow \mathbf{x}(t) - \mathbf{x}\_0 = \frac{1}{k}\int\_{t\_0}^{t} y(t)dt\tag{2}$$

or finally:

434 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

using formulas derived from the Bessel equations.

line models that are developed directly in the time domain.

Electromagnetic Transients Program, or EMTP simply.

be represented the same as a cascade of π circuits.

any computer language or mathematical software.

transmission lines with non-linear elements.

ground and skin effects.

and inductive elements pure.

domain.

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

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

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

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

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

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

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

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

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

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

$$\mathbf{x}(t) = \mathbf{x}\_0 + \frac{1}{k} \int\_{t\_0}^{t} y(t)dt\tag{3}$$

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

$$\int\_{\mathbf{x}(t-\Delta t)}^{\mathbf{x}(t)} \mathbf{x} d\mathbf{x} = \frac{1}{k} d\mathbf{x}(t) = \frac{1}{k} \int\_{t-\Delta t}^{t} y(t) dt \implies \mathbf{x}(t) = \mathbf{x}(t - \Delta t) + \frac{1}{k} \int\_{t-\Delta t}^{t} y(t) dt \tag{4}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 437

It may be noted that the term historical hTR(t-Δ t) depends on x (t) and y (t) after t seconds,

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

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

*R L*

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

*G* 2

*C* 2

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

'

'

*<sup>d</sup> L L*

*<sup>n</sup>* (9)

*<sup>n</sup>* (10)

*<sup>d</sup> R R*

approximate manner and following a series of restrictions, as a cascade of π circuits.

*C* 2

t, is already known during the solution.

**3. Modeling cascada of π circuits** 

is shown in Figure 2.

**Figure 2.** A π circuit.

connected in cascade.

**Figure 3.** Line represented by a cascade of π circuits.

*G* 2

transverse dispersion, respectively. These parameters are written as:

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:

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

or:

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

in other words:

$$\mathbf{x}(t) = \mathbf{C}\_{TR}\mathbf{y}(t) + h\_{TR}\left(t - \Delta t\right) \tag{7}$$

where:

$$\begin{cases} \mathbf{C}\_{TR} = \frac{\Delta t}{2k} \\\\ h\_{TR}(t - \Delta t) = \mathbf{x}(t - \Delta t) + \frac{\Delta t}{2k} \mathbf{y}(t - \Delta t) \end{cases} \tag{8}$$

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

#### **3. Modeling cascada of π circuits**

436 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

curve, limited by the instants t and t-Δt as shown in Figure 1.

( ) ( )

or:

in other words:

where:

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

1 1 <sup>1</sup> () () () ( ) () *x t <sup>t</sup> <sup>t</sup>*

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

(a) (b)

Solving equation (2.5) by the trapezoidal rule integration, which allows a linear variation

1 1 () ( ) () ( ) () ( ) <sup>2</sup>

*yt yt t x t x t t y t dt x t t <sup>t</sup> k k*

> () () ( ) ( ) 2 2 *t t xt yt xt t yt t k k*

> > ( )( ) ( ) <sup>2</sup>

*<sup>t</sup> h t t xt t yt t <sup>k</sup>*

(5)

*xt C yt h t t* () () *TR TR* (7)

(6)

(8)

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

between y (t-Δt) and y (t) of the integrated, as shown in figure 1 (b), as:

*t t t*

2

*TR*

*<sup>t</sup> <sup>C</sup> k*

*TR*

*xt t t t t t xdx dx t y t dt x t x t t y t dt k k <sup>k</sup>* (4)

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' \frac{d}{n} \tag{9}$$

$$L = L' \frac{d}{n} \tag{10}$$

$$G = G' \frac{d}{n} \tag{11}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 439

The models proposed by would become more complete if the effect of frequency on the

So, the parameters of a transmission line can be synthesized by means of a circuit of the type

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

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

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

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 vk(t), respectively. It is also considering that the currents ik0, ik1, ..., ikm are circulating through the inductors L0, L1, L2, ..., Lm, 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 vk(t) state variable. From the currents and voltages in the circuit of Figure 5 it can be

0 0 00 1 1

*di i R Ri ut v dt L L L L* 

1 11 ( )

(13)

*j j kj k*

longitudinal parameters of the line was inserted in them.

**Figure 5.** Cascade of π circuits considering the effect of frequency.

cascade with any quantity of these circuits.

0 0

*m m k k*

*j j L*

shown in Figure 4.

Figure 5.

of frequency.

determined:

$$C = C' \frac{d}{n} \tag{12}$$

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.

438 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

conductance per unit cross-sectional line in length.

can be considered constant and independent of frequency.

are found by solving the set of equations.

amplitudes than are actually.

the parameters of the longitudinal line.

'

'

*<sup>n</sup>* (11)

*<sup>n</sup>* (12)

*<sup>d</sup> G G*

*<sup>d</sup> C C*

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

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

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

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

The representation of transmission lines through cascades of π circuits, taking into account

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,

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

**4. Line representation with parameters of frequency dependent** 

the effect of frequency, is usually implemented in EMTP-type programs.

the quality of results obtained from simulations may be compromised.

**Figure 4.** Circuit that represents the longitudinal parameters of the line.

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 vk(t), respectively. It is also considering that the currents ik0, ik1, ..., ikm are circulating through the inductors L0, L1, L2, ..., Lm, 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 vk(t) state variable. From the currents and voltages in the circuit of Figure 5 it can be determined:

$$\frac{d\dot{t}\_{k0}}{dt} = \frac{\dot{t}\_{k0}}{L\_0} \left( -\sum\_{j=1}^{m} R\_j \right) + \frac{1}{L\_{1,0}} \left( \sum\_{j=1}^{m} R\_j \dot{t}\_{k\_j} \right) + \frac{1}{L\_0} u(t) - \frac{1}{L\_0} v\_k \tag{13}$$

$$\frac{d\dot{i}\_{k1}}{dt} = \frac{R\_1}{L\_1}\dot{i}\_{k0} - \frac{R\_1}{L\_1}\dot{i}\_{k1} \tag{14}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 441

*k k k k km k x iii i v* (21)

*T*

(22)

0

(23)

012

1 2 0

In the equations (4.8) and (4.9), BT e xkT correspond to transposed B and xk, respectively.

*k k k k k m k*

*dx di di di di dv*

*dt dt dt dt dt dt* 

The obtained results show that the vesctor xk have (m + 2) elements and the matrix A is a

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

*A A*

 

*n n*

The elements x1, x2, …, xn 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

> *Akk A*

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

00 0

 

1 0 0

0

*Aik*

1, 2

*C ik kn*

0 00 0 0 0 0 0 0 0

43 3 2

000 0

 

*nn nn*

0 0 00

22 21

*AAA*

12 11 1

*n n n n nn*

1

*A A*

*n n nn*

(24)

(25)

*T*

*k*

dimension n(m +2). The A matrix can be written as:

*A A AAA A A A A <sup>A</sup>*

0

AKK at main diagonal of the matrix A is written as:

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

value <sup>1</sup> *<sup>C</sup>* .

The structure is:

11 12 21 22 23 32 33

*x*

(m + 2) order square matrix.

$$\frac{d\dot{i}\_{k2}}{dt} = \frac{R\_2}{L\_2}\dot{i}\_{k0} - \frac{R\_2}{L\_2}\dot{i}\_{k2} \tag{15}$$

$$\frac{d\dot{I}\_{k\,m}}{dt} = \frac{R\_m}{L\_m}\dot{i}\_{k0} - \frac{R\_m}{L\_m}\dot{i}\_{k\,m} \tag{16}$$

$$\frac{dv\_k(t)}{dt} = \frac{2}{C}\dot{i}\_{k0} - \frac{G}{C}v\_k(t) \tag{17}$$

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

$$\stackrel{\bullet}{\chi} = A\chi + Bu\tag{18}$$

For one circuit, the A matrix is substituted by:

$$A\_{\pi} = \begin{bmatrix} \sum\_{i=0}^{m} R\_{i} & & & & & & & & \underline{R}\_{m} & & & \underline{1} \\ & \underline{I}\_{0} & & \underline{I}\_{0} & & \underline{R}\_{2} & & & & & \underline{1} \\ & & & & & & & & & \\ \hline \\ & \underline{R}\_{1} & & & & & & & & & \\ & \underline{R}\_{1} & & & & & & & & \\ & \underline{R}\_{2} & & & & & & & & \\ & \underline{R}\_{2} & & & & & & & & \\ & \underline{R}\_{2} & & & & & & & & \\ & & \vdots & & & & & & & & \\ & & \vdots & & & & & & & & \\ & & \underline{R}\_{m} & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & & \\ & & &$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 441

$$\mathbf{x}\_k^T = \begin{bmatrix} \mathbf{i}\_{k0} & \mathbf{i}\_{k1} & \mathbf{i}\_{k2} & \cdots & \mathbf{i}\_{km} & \mathbf{v}\_k \end{bmatrix} \tag{21}$$

$$\mathbf{x}\_k = \frac{d\mathbf{x}\_k}{dt} = \left[\frac{d\mathbf{i}\_{k0}}{dt} \quad \frac{d\mathbf{i}\_{k1}}{dt} \quad \frac{d\mathbf{i}\_{k2}}{dt} \quad \cdots \quad \frac{d\mathbf{i}\_{km}}{dt} \quad \frac{d\mathbf{v}\_k}{dt}\right]^T \tag{22}$$

In the equations (4.8) and (4.9), BT e xkT correspond to transposed B and xk, respectively.

The obtained results show that the vesctor xk 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 = \begin{bmatrix} 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\_{(n-3)(n-2)} & 0 & 0\\ 0 & \ddots & 0 & \ddots & A\_{(n-2)(n-2)} & A\_{(n-2)(n-1)} & 0\\ 0 & \ddots & \ddots & \ddots & A\_{(n-1)(n-2)} & A\_{(n-1)(n-1)} & A\_{(n-1)n}\\ 0 & 0 & 0 & \cdots & 0 & A\_{n(n-1)} & A\_{nn} \end{bmatrix} \tag{23}$$

The elements x1, x2, …, xn 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 AKK at main diagonal of the matrix A is written as:

$$A\_{kk} = A\_{\pi} \tag{24}$$

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 <sup>1</sup> *<sup>C</sup>* .

The structure is:

440 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

*k*

*k*

1 1 1

2 2 2

0 *k m m m*

0 ( ) <sup>2</sup> ( ) *<sup>k</sup> k k*

*dv t <sup>G</sup> i vt*

*x Ax Bu*

00 0 0 0

*R R R LL L L L*

0 00

0 0 0

0 0 0

<sup>2</sup> 00 0

<sup>1</sup> 00 00 *<sup>T</sup> <sup>B</sup>*

0

*L*

*C C*

*m m m m*

*R R L L*

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

*j m*

0 1 2

2 2 2 2

*R R L L*

1 1 0 0

*R R L L*

For one circuit, the A matrix is substituted by:

*j m*

*A*

*j*

*R*

*di R R i i*

*di R R i i*

*di R R i i*

0 1 1 1

0 2 2 2

*k km m m*

*k k*

*dt L L* (14)

*dt L L* (15)

*dt L L* (16)

*dt C C* (17)

(18)

0 0

*G*

(20)

1

(19)

*k k*

$$A\_{ik} = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \cdots & \vdots \\ 0 & \cdots & \ddots & \vdots \\ -1 \bigvee\_{C} & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \end{bmatrix} \tag{25}$$
  $i = k - 1, \quad 2 \le k \le n$ 

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

0 1 *L* . The structure is:

$$A\_{ik} = \begin{bmatrix} 0 & \cdots & 0 & 1 \\ \vdots & \ddots & \cdots & 0 \\ \vdots & \cdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix} \\ \tag{26}$$
  $i = k + 1, \quad 1 \le k \le n - 1$ 

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 443

Equation (5.2) is applied Simpson's rule for solving the state equation of state of the type

<sup>1</sup> 6

 <sup>2</sup> <sup>6</sup> *<sup>k</sup> <sup>t</sup> a I A xt*

> <sup>3</sup> 2

 4 1 4 6 *k Mk <sup>t</sup> a But ut ut*

> 1 2 *k k*

**Figure 6.** Approximation of the derivative of the function to be integrated (solid curve) by a function of

 <sup>1</sup> 1 2 34 <sup>1</sup> 2 2

<sup>6</sup> *k k xt xt k k k k* (35)

*t t*

*M*

3 *<sup>M</sup>*

shown in (5.1). In this case, the matrices *a1*, *a2*, *a3* and *a4* are written as:

*<sup>k</sup>* 1 12 3 4 *xt a a a a* (28)

*<sup>t</sup> aI A* (29)

(30)

*<sup>t</sup> a Ay t* (31)

*k k* <sup>1</sup> *tt t* (33)

*<sup>t</sup>* (34)

(32)

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

It is defined the following terms:

the second degree (dashed curve).

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

**5.2. Runge-Kutta's rule** 

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 0 1 *L* .

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 MatLabTM 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'\left(t\_{k+1}\right) + 4y'\left(t\_M\right) + y'\left(t\_k\right) \right] \Delta t \tag{27}$$

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

442 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

*ik*

*A*

element, which is the first array element and it has the value

simulation of transient phenomena in power systems.

0 1 *L* 

. The structure is:

described in the next item

**5.1. Simpson's rule** 

function as shown in Figure 6. From Figure 1, it is obtained:

**5. Other numeric methods** 

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

1 0 0

 

00 0 1, 1 1

*ik kn*

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

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

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 MatLabTM 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

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

> 1 1 <sup>1</sup> ' 4' '

<sup>3</sup> *k k Mk <sup>y</sup> t yt yt yt t* (27)

0

(26)

0 1 *L* .

*L*

0

$$\propto \left( t\_{k+1} \right) = a\_1 \left( a\_2 + a\_3 + a\_4 \right) \tag{28}$$

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 *a1*, *a2*, *a3* and *a4* are written as:

$$a\_1 = \left(I - \frac{\Delta t}{6}A\right) \tag{29}$$

$$a\_2 = \left(I + \frac{\Delta t}{6}A\right) \mathbf{x}\left(t\_k\right) \tag{30}$$

$$a\_3 = \frac{2\Lambda t}{3} A \, y\left(t\_M\right) \tag{31}$$

$$
\Delta a\_4 = \frac{\Delta t}{6} B \left[ \mu \left( t\_k \right) + 4 \mu \left( t\_M \right) + \mu \left( t\_{k+1} \right) \right] \tag{32}
$$

It is defined the following terms:

$$
\Delta t = t\_{k+1} - t\_k \tag{33}
$$

$$t\_M = \frac{t\_{k+1} + t\_k}{2} \tag{34}$$

**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:

$$\chi\left(t\_{k+1}\right) = \chi\left(t\_k\right) + \frac{1}{6}\left(k\_1 + 2k\_2 + 2k\_3 + k\_4\right) \tag{35}$$

The terms in the last equation are described by:

$$k\_1 = \Delta t \left[ A \ge \left( t\_k \right) + B \, u \left( t\_k \right) \right] \tag{36}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 445

Computational effort (simulation time) Heun Simpson Runge-Kutta

important for analysis in very short simulation times. The next table shows the simulation

Table 1 shows the time relationship between the applied numerical methods. It is possible to carry out these comparisons because the MatlabTM 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 MatLabTM is an adequate tool for developing new models and numeric routines used for analysis and simulations of transient

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

The model that represents the single-phase line was implemented in a microcomputer using

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

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

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

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

time comparisons obtained with the MatLabTM software.

**Table 1.** Simulation time required by numerical methods

**6. Model implementation of single phase line** 

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

parameters are fitted by means of rational functions.

Step of calculating

phenomena in transmission lines.

with results obtained with EMTP.

the software MatLabTM.

frequency.

point of the line.

(µs)

$$k\_2 = \Delta t \left\{ A \left[ \mathbf{x} \left( t\_k \right) + \frac{k\_1}{2} \right] + B \boldsymbol{u} \left( t\_k \right) \right\} \tag{37}$$

$$k\_3 = \Delta t \left\{ A \left[ x \left( t\_k \right) + \frac{k\_2}{2} \right] + B u \left( t\_k \right) \right\} \tag{38}$$

$$k\_4 = \Delta t \left\{ A \left[ \mathbf{x} \left( t\_k \right) + k\_3 \right] + B \boldsymbol{u} \left( t\_k \right) \right\} \tag{39}$$

#### **5.3. Comparisons among the numeric methods**

The MatlabTM 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 MatLabTM software.

**Table 1.** Simulation time required by numerical methods

444 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

<sup>1</sup> *k k k t Ax t Bu t* (36)

*k k* (39)

**Heun Simpson Runge-Kutta**

(37)

(38)

<sup>1</sup>

<sup>2</sup>

<sup>2</sup> 2 *k k <sup>k</sup> k t A x t Bu t*

<sup>3</sup> 2 *k k <sup>k</sup> k t A x t Bu t*

*k t A x t k Bu t* 4 3

The MatlabTM 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

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

**0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6**

**Time [ms]**

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

The terms in the last equation are described by:

**5.3. Comparisons among the numeric methods** 

**Figure 7.** Comparisons among the three numerical methods.

methods: trapezoidal rule, Simpson's rule and Range-Kutta's one.

a single-phase transmission line.

**-10**

**0**

**10**

**20**

**Volt [kV]**

**30**

**40**

**50**

Table 1 shows the time relationship between the applied numerical methods. It is possible to carry out these comparisons because the MatlabTM 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 MatLabTM 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 MatLabTM.

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 *x Ax Bu* [ ] [] , it's used Heun formula (trapezoidal integration method). This method is widely used in simulations of electromagnetic transients in power systems

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 447

In the equivalent circuit, the values were considered in the frequency range from 101 to 106, 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

> 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

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

longitudinal parameters of single-phase transmission line.

**Figure 9.** Circuit that represents the longitudinal parameters of the line.

**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** 

recombination, also releases excess energy.

Figure 8 is shown a block diagram of the algorithm of the program developed for singlephase 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\_{\rm VIT}(oo) = R\_0 + joI\_0 + \frac{joR\_1}{jo + \frac{R\_1}{L\_1}} + \frac{joR\_2}{jo + \frac{R\_2}{L\_2}} + \frac{joR\_3}{jo + \frac{R\_3}{L\_3}} + \frac{joR\_4}{jo + \frac{R\_4}{L\_4}}\tag{40}$$

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 101 to 106 for it is possible to calculate the longitudinal impedance synthesized.

In the equivalent circuit, the values were considered in the frequency range from 101 to 106, 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.


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

446 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

simulations of electromagnetic transients in power systems

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

0 0

 

( ) *VFIT*

Table 2.

impedance synthesized.

**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:

*jR jR j R j R Z R jL RRRR jj j <sup>j</sup> LL L <sup>L</sup>* 

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

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 101 to 106 for it is possible to calculate the longitudinal

1 2 3 4

1234 1 2 3 4

 

(40)

 

phase 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 *x Ax Bu* [ ] []

it's used Heun formula (trapezoidal integration method). This method is widely used in

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

## **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.

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 449

*C*

0,22 1,2 *r* (42)

(43)

(41)

*C*

*se v V*

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,

1

In equation (7.1), CC 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, CV is the minimum voltage required to the corona effect accurrance and η is a coefficient

Thus, given the presence of the corona effect and the effect of frequency, a differential

<sup>2</sup> <sup>1</sup> *<sup>C</sup>*

*v* 

Gary's model considers that the corona effect is manifested only if the voltage VC 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 Vp at this point must satisfy the

*C C <sup>V</sup> G k*

*<sup>v</sup> <sup>C</sup> se v V*

*C*

*V*

where currents and voltages along line is described by state variables.

*C*

*C*

where: r the radius of the conductor in centimeters.

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:

defined as:

following conditions:

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

0

  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 nonlinear 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:

$$\mathbf{C}\_{\mathcal{C}} = \begin{cases} \mathbf{C}\eta \left(\frac{\upsilon}{V\_{\mathcal{C}}}\right)^{\eta - 1} & \text{se} & \upsilon \ge V\_{\mathcal{C}} \\\\ 0 & \text{se} & \upsilon < V\_{\mathcal{C}} \\\\ \end{cases} \tag{41}$$

In equation (7.1), CC 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, CV is the minimum voltage required to the corona effect accurrance and η is a coefficient defined as:

$$
\eta = 0, \mathcal{D}\mathcal{D}r + 1, \mathcal{D} \tag{42}
$$

where: r the radius of the conductor in centimeters.

448 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

contributing to the understanding of the basic mechanism of the corona.

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,

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

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

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

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

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

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 nonlinear 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

made of in the laboratory of a conductor under corona (project called Tidd 500 kV).

**7.2. Corona effect representation** 

and uses abacuses.

capacitance and dynamic tension.

ionization.

this nature yet.

**7.3. Gary's model** 

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\_{\mathbb{C}} = k\_{\mathbb{C}} \left( \frac{1 - V\_{\mathbb{C}}}{v} \right)^2 \tag{43}$$

Gary's model considers that the corona effect is manifested only if the voltage VC 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 Vp at this point must satisfy the following conditions:

$$dV\_p > V\_C \quad \text{and} \quad \frac{dV\_p}{dt} > 0 \tag{44}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 451

simulation period is 600 µs. The voltage in the initial of the line is 1 kV. It can also be

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 MatLabTM program, using the trapezoidal integration method. Considering the frequency independent line parameters, the circuit of Figure 4 is reduced to the R0 and L0

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

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 R4 value per

> 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

> > 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

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

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

elements. In this case, the values are: R0 = 0.05 Ω/km and L0 = 1 mH/km.

[H], respectively. In the Figs 12-21, it is used the values of Table 4.

**Table 3.** Values of line parameter used in all simulations without frequency.

**Table 4.** Values of line parameter used in simulations with 200 circuits.

increase of the number of RL parallel blocks.

considered as 1 pu.

length unit.

point.

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.

450 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

*p C*

voltage.

**8. Testing the model** 

C′ =11,11nF/km.

**representation** 

and frequency dependent ones.

u(t)

**Swith S**

**Figure 11.** Single phase line representation with opened terminal.

structure is completed using a RL series block as shown in Figure 9.

**9. Simulation analysis of transmission lines with single phase** 

*V V and*

0 *<sup>p</sup>*

(44)

*dV*

*dt*

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

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

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

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

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 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 MatLabTM program, using the trapezoidal integration method. Considering the frequency independent line parameters, the circuit of Figure 4 is reduced to the R0 and L0 elements. In this case, the values are: R0 = 0.05 Ω/km and L0 = 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 R4 value per length unit.


**Table 3.** Values of line parameter used in all simulations without frequency.


**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.

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 453

**Figure 15.** Synthesis of inductance per unit length using only R'4L'4 parallel block.

**103**

**0.05 <sup>6</sup>**

**Frequency**

**104**

**104**

**105**

**105**

**10**

**10**

**102**

**102**

**101**

**101**

**2.21 2.22 2.23 2.24 2.25 2.26 2.27**

**2.25**

**2.3**

**2.35**

**H/km**

**2.4**

**2.45**

**H/km**

**0.05**

**0.05**

**0.05**

**H/km**

**0.05**

**0.05**

**0.05**

**Figure 16.** Synthesis of inductance per unit length using R'1L'1 and R'4L'4 parallel blocks.

**103**

**2. 2 <sup>6</sup>**

**Frequency**

**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

**10<sup>1</sup> 102 103 10<sup>4</sup> 105 10 2.2 <sup>6</sup>**

**Frequency**

**Figure 12.** Synthesis of resistance per unit length using R'4L'4 parallel block.

**Figure 13.** Synthesis of resistance per unit length using R'1L'1 and R'4L'4 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'4L'4 parallel block.

**Figure 12.** Synthesis of resistance per unit length using R'4L'4 parallel block.

**110.5**

**101**

**101**

**40**

**60**

**80**

**Ohm /km** **100**

**120**

**140**

**20**

**40**

**60**

**Ohm /km** **80**

**100**

**120**

**111**

**111.5**

**Ohm/km**

**112**

**112.5**

**<sup>10</sup><sup>1</sup> <sup>10</sup><sup>2</sup> <sup>10</sup><sup>3</sup> <sup>10</sup><sup>4</sup> <sup>105</sup> <sup>10</sup><sup>6</sup> <sup>110</sup>**

**Frequency**

**Figure 13.** Synthesis of resistance per unit length using R'1L'1 and R'4L'4 parallel blocks.

**103**

**104**

**104**

**0 <sup>6</sup>**

**Frequency**

**105**

**105**

**10**

**10**

**102**

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

**103**

**20 <sup>6</sup>**

**Frequency**

**102**

**Figure 16.** Synthesis of inductance per unit length using R'1L'1 and R'4L'4 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.

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 455

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

**Figure 20.** Current at the end of the transmission line without the effect of frequency.

**0 1 2 3 4 5 6**

**Time [s]**

**x 10-4**

**x 10-4**

**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

**0 1 2 3 4 5 6**

**Time [s]**

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

consider the frequency influence.

**-1**

**<sup>6</sup> x 10-4**

**-0.5**

**0**

**0.5**

**Current [A]**

**1**

**1.5 x 10-3**

initial line terminal.

and 200 π circuits, respectively.

**Current [A]**

**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.

454 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

through the line.

circuits.

circuits.

make up the output voltage shown.

**0**

**0.5**

**1**

**Voltage [pu]**

**1.5**

**2**

**-1 -0.5 0 0.5 1 1.5 2 2.5 3**

**Voltage [pu]**

because the inclusion of the frequency influence.

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

**Figure 18.** Energization of the transmission line without the effect of frequency considering 200 π

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

**0 1 2 3 4 5 6**

**Time [s]**

**x 10-4**

**x 10-4**

**Figure 19.** Voltage at the end of the transmission line with the effect of frequency considering 200 π

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

**0 1 2 3 4 5 6**

**Time [s]**

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

**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.

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 457

**x 10 -4**

**x 10-4**

**x 10-5**

**Figure 24.** Voltage at the end of the line with the effect of frequency considering 200 π circuits and

**0 1 2 3 4 5 6**

**Time [s]**

**Figure 25.** Voltage at the end of the line with the effect of frequency considering 200 π circuit and using

**0 1 2 3 4 5 6**

**Time [s]**

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

**All RL parallel blocks One RL parallel block Three RL parallel blocks**

**4 5 6 7 8 9 10**

**Time [s]**

using R4L4 parallel block.

**0**

**0.5**

**1**

**Voltage [kV]**

**1.5**

**2**

R0L0, R2L2 e R4L4 parallel blocks.

**1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 2.05**

**Voltage [kV]**

**0**

**0.5**

**1**

**Voltage [kV]**

**1.5**

**2**

frequency.

**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 22.** Specific portion of the simulation at the end of the line with 100 π circuits and without the

**3 4 5 6 7 8 9 10**

**Time [s]**

**x 10-5**

**x 10-5**

Based on the mentioned results, it concludes that increasing the number of π circuits leads

**Figure 23.** Specific portion of the simulation at the end of the line with 200 π circuits and without the

**3 4 5 6 7 8 9 10**

**Time [s]**

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

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

values of Tables 3 and 1, as well as, 3 and 4 RL parallel blocks, respectively.

the voltage in the line end has a smoothly curve in the time domain.

effect of frequency.

effect of frequency.

to a decrease in the numerical oscillations.

**1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5**

**Voltage [kV]**

**1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5**

**Voltage [kV]**

**Figure 24.** Voltage at the end of the line with the effect of frequency considering 200 π circuits and using R4L4 parallel block.

**Figure 25.** Voltage at the end of the line with the effect of frequency considering 200 π circuit and using R0L0, R2L2 e R4L4 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{aligned} C\_{\odot} &= \mathbf{C} \cdot \eta \cdot \left(\frac{V}{V\_{\odot}}\right)^{\eta - 1} \\ \eta &= 0.22 \cdot R\_{\rm COD} + 1.2 \end{aligned} \tag{45}$$

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 459

**Figure 28.** The time domain simulations with the corona effect for VCORONA = 0,5 V.

**Figure 29.** The time domain simulations with the corona effect for VCORONA = 0,7 V.

**Figure 30.** The time domain simulations with the corona effect for VCORONA = 0,9 V.

Based on Figure 10, it's simulated the corona effect with effect of frequency for some

relations between V and VC. The V value is 1 pu to the following shown simulation.

In this case, RCOND is the conductor phase radius in [cm], C is the transversal line capacitance and the CC is the new value of the transversal line capacitance because the corona effect. In this chapter, the RCOND value is 2.54 cm. The simulations are carried out for some relations between V and VC. 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 VCORONA = 0,35 V.

**Figure 28.** The time domain simulations with the corona effect for VCORONA = 0,5 V.

*C*

between V and VC. The V value is 1 pu for all following shown simulations.

**Figure 27.** The time domain simulations with the corona effect for VCORONA = 0,35 V.

*<sup>V</sup> C C*

 

In this case, RCOND is the conductor phase radius in [cm], C is the transversal line capacitance and the CC is the new value of the transversal line capacitance because the corona effect. In this chapter, the RCOND value is 2.54 cm. The simulations are carried out for some relations

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

in mathematical software.

introduced by

transmission lines.

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

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

0.22 1.2

*V R*

*C COND* 1

(45)

**Figure 29.** The time domain simulations with the corona effect for VCORONA = 0,7 V.

**Figure 30.** The time domain simulations with the corona effect for VCORONA = 0,9 V.

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

The MatLabTM Software Application for Electrical Engineering Simulations and Power System 461

parameters. These parameters have their characteristics distributed along the line and this is

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 MatLabTM software is an excellent

tool for basic and profound studies of transient phenomena in transmission lines.

*Transient Electromagnetic Studies Laboratory, Department of Electrical Engineering,* 

*Department of Electrical Engineering, UNICAMP – The State University of Campinas,* 

Chipman, R. A. Teoria e problemas de linhas de transmissão. São Paulo: Mc Graw Hill do

Dommel H.W., Electromagnetic Transients Program. Reference Manual (EMTP Theory

Dommel, H.W. Digital computer solution of electromagnetic transients in single and multiphase networks. IEEE Trans. On Power App. And Systems, v. PAS-88, n.4, p. 388- 399, 1969. Faria, A.B.; Washington, L.A.; Antônio, C.S. Modelos de linhas de transmissão no domínio das fases: estado da arte. In: CONGRESSO BRASILEIRO DE AUTOMÁTICA, 14, 2002,

Fuchs, R. D. Transmissão de energia elétrica: linha aéreas; teoria das linhas em regime

Greenwood, A. Electrical transients in power systems. New York: John Wiley&Sons, 1971. 558 p. Hedman, D. E. Teorias das linhas de transmissão-II. 2.ed. Santa Maria: Edições UFSM, 1983.

J. R. Marti, "Accurate modelem of frequency-dependent transmission lines in 105 electromagnetic transient simulations", IEEE Trans. on Power Apparatus and Systems,

Kurokawa, S.; Yamanaka, F. N. R.; Prado, A. J.; Pissolato Filho, J.. Using state-space techniques to represent frequency dependent single-phase lines directly in time domain. In: THE 2008 IEEE/PES Transmission and Distribution Conference and Exposition: Latin America, 2008, Bogotá. Proceedings, Bogotá:[s.n.], 2008. p. 312-316,.

permanente, 2.ed. Rio de Janeiro: Livros Técnicos e Científicos, 1979. 582 p.

Leonardo da S. Lessa, Afonso J. Prado, Rodrigo Cleber Silva,

Book), Bonneville Power Administration, Portland, 1986.

Natal. Anais... Natal: [s.n.], 2002. p. 801-806.

vol. PAS-101, nº 1, pp. 147-155, January, 1982.

considered in the mentioned routine.

Sérgio Kurokawa and Luiz F. Bovolato

*UNESP – The University State of São Paulo, Brazil* 

**Author details** 

José Pissolato Filho

*Campinas, Brazil* 

**11. References** 

v. 2 e 3.

Brasil, 1976. 276p.

**Figure 31.** Comparisons of the time domain simulations with the corona effect for some VCORONA values.

**Figure 32.** Comparisons of the time domain simulations with the corona effect with effect of frequency for some VCORONA values.

#### **10. Conclusions**

It is described the application the MatLabTM 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 MatLabTM 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 MatLabTM 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 MatLabTM software is an excellent tool for basic and profound studies of transient phenomena in transmission lines.

## **Author details**

460 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

**Figure 31.** Comparisons of the time domain simulations with the corona effect for some VCORONA values.

**0,35V 0,5V 0,7V 0,8V 0,9V**

**Figure 32.** Comparisons of the time domain simulations with the corona effect with effect of frequency

**0 1 2 3 4 5 6**

**Tempo [s]**

**x 10-4**

It is described the application the MatLabTM 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 MatLabTM 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 MatLabTM 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

for some VCORONA values.

**-0.5**

**0**

**0.5**

**Voltage [pu]**

**1**

**1.5**

**2**

**10. Conclusions** 

Leonardo da S. Lessa, Afonso J. Prado, Rodrigo Cleber Silva, Sérgio Kurokawa and Luiz F. Bovolato *Transient Electromagnetic Studies Laboratory, Department of Electrical Engineering, UNESP – The University State of São Paulo, Brazil* 

José Pissolato Filho *Department of Electrical Engineering, UNICAMP – The State University of Campinas, Campinas, Brazil* 

## **11. References**


Kurokawa, S.; Yamanaka, F. N. R.; Prado, A. J. Representação de linhas de transmissão por meio de variáveis de estado levando em consideração o efeito da freqüência sobre os parâmetros longitudinais. Sba Controle& Automação, Campinas, v.18, n.3, p.337- 346, 2007.

**Chapter 19** 

© 2012 Monzani et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Monzani et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Using a Low Complexity Numeric Routine for** 

Rafael Cuerda Monzani, Afonso José do Prado, Sérgio Kurokawa,

Luiz Fernando Bovolato and José Pissolato Filho

Additional information is available at the end of the chapter

electromagnetic transients, at any point of the circuit.

basis of electromagnetic transients simulations in power systems.

http://dx.doi.org/10.5772/48507

programs of EMTP type are:

development since 1988.

**1. Introduction** 

**Solving Electromagnetic Transient Simulations** 

Different methods can be applied to accomplish the analysis of transmission lines. Many mathematical tools can be used, the main tools used are: circuits analysis with the use of Laplace or Fourier Transform, State Variables and Differential Equations. These tools can be included in a numeric routine in order to obtain voltage and current values in simulation of

The EMTP (ElectroMagnetic Transient Program) [1] is the main kind of this software. The prototype was developed in 60's by professionals of power system area led by Dr. Hermann Dommel (University of British Columbia, in Vancouver, B.C., Canada), and Dr. Scott Meyer (Bonneville Power Administration in Portland, Oregon, U.S.A.). Currently, EMTP is the

With EMTP type programs, the following analysis can be done: simulation of switching and lightning surges, transient and temporary overvoltage, electrical machines, resonance phenomena, harmonics, power quality and power electronics applications. The most known

 MicroTran Power Systems Analysis of the University of British Columbia, Vancouver, Canada, founded in 1987 by: Hermann W. Dommel, Jose R. Marti (University of British Columbia), and Luis Marti (University of Western Ontario, Hydro One Networks Inc.). PSCAD®, also known as PSCAD®/EMTDC™ of Manitoba HVDC Research Centre. Commercially available since 1993, PSCAD® is the result of continuous research and

 ATP has been continuously developed through international contributions by Drs. W. Scott Meyer and Tsu-huei Liu, the co-Chairmen of the Canadian/American EMTP User

