Abstract
This paper presents a multiple time-step solver based on a time-domain explicit integration algorithm for improving the computational speed of high-voltage transmission line electromagnetic transient (EMT) simulation. For weakly rigid EMT models of high-voltage transmission lines, the previous precise Runge-Kutta integration method with a small time step is adopted; for strongly rigid nonlinear EMT models of high-voltage transmission lines, a an improved precise integration method using a large time step is used for the solver of EMT simulation. Practical simulations for overvoltages of high-voltage transmission line show that multiple time-step EMT solver can be applied to different cases of EMT simulation for transmission line.
1. Introduction
As is known to all, the EMT simulation program (EMTP) [1, 2] is widely used to calculate the overvoltage of power systems. Overvoltage calculation often is implemented by solving different ordinary differential or partial differential equations that are used to describe different electrical phenomena or processes. Therefore, the numerical algorithms of the differential equation become a key part of the EMT simulation. In EMT-type simulation tools including EMTP, PSCAD/EMTDC [1], NETOMAC, and RTDS, EMTP uses the implicit trapezoidal integration method with second-order accuracy and A-stability to discretize the differential equation of basic electrical components (inductors and capacitors), and the discrete time-domain solution of state variables is obtained by solving the algebraic equation [2, 3].
The implicit trapezoidal integration method is limited in integration accuracy and numerical stability. On the one hand, when the power grid topology changes abruptly in EMT simulation, voltage or current simulation waveform will produce a numerical oscillation phenomenon [2, 3]. A typical and simple simulation example is the calculation of overvoltage of overhead transmission line without load being switched suddenly. When the trapezoidal method is used to calculate the overvoltage at the end of the transmission line, serious numerical oscillation occurs in the simulation result. This simple case shows that the trapezoidal method cannot accurately simulate the EMT process for the change of topology in power grids. On the other hand, with the continuous access to distributed clean energy, the number of power grid nodes has increased dramatically, resulting in a very large calculation scale for EMT simulation. For such systems, the computation effect of implicit numerical integration algorithms does not perform as fast as explicit integration algorithms. This is mainly because the mathematical model of distributed circuit elements contains many power electronic devices, limiting links, and multiscale transient processes [4, 5].
To solve this problem, multirate simulation and parallel computing are used to perform EMTs of large power grids [6, 7]. A Parareal [8] parallel computation method is proposed to be applied to EMT simulation, and it is verified that the parallel algorithm can obtain an effective acceleration ratio, but the algorithm cannot suppress the numerical oscillation. Graphics processing units are used to accelerate EMTs of power systems and adopt thread-oriented conversion and automatic code generation technology [9]. Reference [10] proposed the multirate EMT simulation using latency technology. This method uses the extended implicit trapezoidal integration method to realize the synchronization process between subnets with different simulation steps, but its computational amount is large. Reference [11] proposed the concept of “double-layer grid separation” for the topology characteristics of AC-DC hybrid system. This method greatly improved the simulation efficiency of AC-DC hybrid power grids.
To replace the role of implicit trapezoidal integration in EMT simulation, researchers have begun to seek the application of higher-order and A-stable or L-stable numerical algorithms in EMT simulation. The two-stage diagonally implicit Runge-Kutta method (2S-DIRK) [12] is applied to EMT simulation. The 2S-DIRK is L-stable as same as the implicit Euler method, so it can effectively avoid numerical oscillation [2]. Two-stage three-order Radau IIA method [13] is used for a new parallel method for solving nonlinear EMT state-space equation at multiple time points, and then the Newton method is used to solve it iteratively. The algorithm not only has high computational efficiency but also can avoid numerical oscillation.
Both the precise integration method [14] and the differential quadrature (DQ) method [15] have high accuracy and strong stability, and both of them have been widely used in the field of computational mechanics [16, 17]. In order to improve the applicable scope of the precise integration method, the precise Runge-Kutta method is proposed, which is a new algorithm formed by the precise integration method and the explicit Runge-Kutta method with the fourth order [18]. The calculation efficiency and precision of the precise Runge-Kutta method need to be improved, so it cannot be directly used in the EMT simulation of the power systems. In order to improve the EMT simulation efficiency via an improved precise integration method, this paper makes full use of the advantages of the DQ and precise integration method to construct a new method to solve the problem. The improved precise integration method has high accuracy and good stability, and it can use a large time step to improve the simulation efficiency of EMTs. Compared with the implicit trapezoidal integration using a small time step, the calculation error of this method meets the needs of practical engineering calculation and the calculation efficiency is considerable.
Compared with the implicit numerical algorithm, the explicit numerical algorithm has higher execution efficiency. However, the existing explicit methods, such as the explicit Euler method, have too low numerical precision to be directly used in power system transient simulation. Fortunately, the new method proposed in this paper is essentially an explicit single-step high-precision integration method. The idea of applying the improved precise integration method to electromagnetic transient simulation proposed in this paper is the biggest advantage over the existing implicit electromagnetic transient simulation algorithm [19, 20] in calculation efficiency, which will be confirmed later in this paper.
The remainder of this paper is organized as follows. In Section 2, an improved precise integration method based on DQ with different grids is presented in detail. In Section 3, numerical simulations for overvoltages of high-voltage transmission lines with different time scale and complexity are given to verify the accuracy and effectiveness of the improved precise integration method. Finally, the conclusions derived from the study are presented in Section 4.
2. Improved Precise Integration Method
Generally, the state-space model of EMTs in power systems can be described as an initial value problem of the first-order ordinary differential equation:where represents time variable; is a function of time and state variable ; is the value of state variable at initial time .
Formula (1) can always be rewritten into ordinary differential equations consisting of homogeneous and nonhomogeneous parts by a certain transformation, andwhere is the coefficient matrix related to system parameters. To the right of the medium sign in equation (2), the first item is called homogeneous part, and the second item is called nonhomogeneous term. The exact solution of the homogeneous equation of equation (2) on the computer is as follows:
In the single-step integration, equation (3) has the following discrete recursive calculation format:
The precise integration method (PIM) [14] proposed by Zhong uses additive theorem and incremental storage to calculate exponential matrix :where (where ). For very small , we can use Taylor series to expand the exponential matrix in (5) into M terms (where M = 4); then, there are
Substituting formula (6) into formula (5) yields
According to the matrix multiplication, the following recursive operations are performed iteratively:where . At the end of the cycle,
In the single integration step , the solution of the nonhomogeneous dynamic equation (2) can be expressed as follows:
In formula (10), the exponential matrix in the first term on the right side of equation (10) has been obtained by the precise integration method, while the second integration is related to the characteristics of the power system, which is called the Duhamel integration term. Because the calculation of the first term can be achieved by PIM on computer, the numerical error mainly comes from the numerical calculation error of Duhamel integration term. In this paper, Duhamel integration terms are calculated by time-domain differential quadrature scheme based on different grids.
2.1. Uniform Grid
Taking the time-domain differential quadrature method [15] as an example, the integration format of Duhamel integration term in equation (10) is derived. For the initial value problem , the s-th formula of the DQ scheme with s-stage s-order can be expressed as follows:where is the time step; is the grid point; are the integration coefficients related to the grid points of DQM.
The DQM uses uniform grid, and the formula of Duhamel integration term in equation (10) is as follows:where .
The explicit numerical method can be used for estimation and then the second equation in formula (11) can be directly used for single-step numerical integration. Specifically, when , the four-order explicit Runge-Kutta method is used to calculate the estimated value as follows:where ;
The other values of , , and are calculated in turn according to formula (13), and the approximate value of Duhamel integration term in equation (12) is obtained as follows:where ; ; ; ; ; ; ; there is
It can be seen from formula (16) that only the first exponential matrix T3 can be calculated by the PIM in the calculation of this method, so that the amount of calculation will not be too large because of the increase in the number of integration nodes. After finding the approximate value of the Duhamel integration term, the approximate value of can be obtained by substituting it into equation (10).
2.2. Nonequidistant Grids
The commonly used nonequidistant grids are Legendre grid, Chebyshev grid, and Chebyshev-Gauss-Lobatto grid [21]. Consider the distribution of Chebyshev grid points on regularized interval [0, 1], where
Similarly, when , the four-order explicit Runge-Kutta method is used to calculate the estimated value as follows:
In equation (18), ;
Then the approximation of the Duhamel integration term is as follows:where ; ; ; , where ; ; ; .
The approximation of can be obtained by substituting into equation (10). Obviously, using nonuniform grids needs to calculate exponential matrices, which requires more extra calculation than using uniform grid.
3. Numerical Simulations
3.1. Accuracy Test for Improved Precise Integration Method
As shown in Figure 1, this circuit represents the equivalent schematic diagram of a 515 kV bus switching AC filter in a substation. is the equivalent inductance of the AC grid; represents the equivalent resistance of the AC grid; denotes the equivalent electromotive force of the AC system; is the initial phase of power supply es. The capacitor C represents a simplified equivalent circuit of an AC filter bank.
When the circuit breaker K is closed at t = 0 s, the AC filter is fully discharged and there is no current in the branch of inductor. In Figure 1, there are
In general, the system is usually an underdamped or attenuated oscillation system. At this time, the characteristic equation of differential equation (21) has a pair of conjugate complex roots with a negative real part. Under this condition, the voltage across the AC filter can be obtained aswhere ; ;
In Figure 1, the decay constant of resistance, inductance, and capacitor series circuit is ; resonance angular frequency of this circuit is ; and natural angular frequency of this circuit is [22, 23].
In this simulation, the lumped electrical parameters are , , and . The proposed methods were used to solve equation (22) and implemented with the MATLAB scripting language. Then, we calculated the voltage of the capacitor C by improved precise integration method (with uniform grid and s = 4), precise Runge-Kutta integration method (four-order explicit Runge-Kutta is used with precise integration method), the two-stage diagonally implicit Runge-Kutta (2S-DIRK with second order in accuracy), and the trapezoidal method (TR) or critical damping adjustment [19]. Figures 2–4 are the computational results of this simulation test calculated by precise differential integration method (PDIM, h = 0.1 ms), precise Runge-Kutta integration method [18, 24] (PRKM, h = 0.1 ms), the two-stage diagonally implicit Runge-Kutta (2S-DIRK, h = 0.01 ms), and critical damping adjustment (CDA, h = 0.01 ms) or TR (h = 0.01 ms). And the curve of difference values by these three methods relative to the true solution is also given in Figures 3 and 4. The simulation starts from the zero initial state. The simulation ends at t = 0.1 s.
As shown in Figures 2–4, PDIM has a better computational result compared with TR and 2S-DIRK in the bigger time step. It seems that the PRKM is as accurate as PDIM for this case. And even with large time step, the computational accuracy of PDIM is two orders of magnitude higher than that of the trapezoidal method and 2S-DIRK. In addition, in order to further verify the effectiveness of the proposed method in suppressing numerical oscillations, a nonlinear reactor switching circuit is used for simulation analysis. In Figure 5, the lumped electrical parameters are equivalent resistance, , and equivalent inductor, . The piecewise linear current-flux curve of is given in Figure 6. The peak voltage of the sinusoidal voltage source is and . The PDIM and the trapezoidal method (TR) give completely different results and their differences are quite large. The result obtained by the trapezoidal method is shown in Figure 7, where sustained numerical oscillations are observed while the results obtained by PDIM are not. The calculated results obtained using a time step of 0.01 ms are shown in Figure 7.
(a)
(b)
(c)
3.2. Single-Phase Transmission Line with Nonlinear Inductance Load
This case is a high-voltage transmission line with a nonlinear inductance load as shown in Figure 8. In Figure 8, is the excitation AC voltage source; and are the internal resistance and internal inductance of the voltage source at the sending end; the switch S is suddenly closed at time t = 0 s. The total length of the line is L = 100 km, and the distributed parameters of the transmission line are resistance , inductance , and capacitance . The nonlinear load is composed of load resistance and the nonlinear inductor in parallel connection, where the relationship between the magnetic link and current of the nonlinear inductor is ; the electrical parameters in this case are shown in Table 1.
The electrical model used to describe the EMT process of the high-voltage transmission line shown in Figure 8 is the telegraph equation. The telegraph equation is a hyperbolic partial differential equation, and it is necessary to transform it into ordinary differential equations in the form of equation (2) before different numerical methods are used for EMT simulation. We divide the whole transmission line into M sections uniformly as described in Figure 9. The Π-type cascade equivalent circuit model was carried out on the transmission lines in Figure 8. So, the resistance, inductance, and capacitance of each section are as follows:
As shown in Figure 9, it is easy to establish the following first-order linear ordinary differential equations (ODEs) for the transmission line using Kirchhoff’s voltage and current law [22, 23]:
Considering the boundary conditions at the head and end of the transmission line in Figure 9, equation (26) is arranged into the following matrix form:where is the constant sparse matrix; is sparse column vector, which is the excitation source of EMT simulation for the transmission line;
In Figure 8, the number of interval segments is after the transmission line is discretized in spatial domain. The CDA and the PDIM (with uniform grid and s = 4) are used to solve equation (27) to obtain the voltage variation curve of the nonlinear inductance load. The switch S is initially open and is closed at t = 0 s. The CDA method uses small time step and PDIM uses a larger time step to calculate transient voltage of the transmission line. And the sending and receiving end voltage are shown in Figures 10 and 11.
As shown in Figures 10 and 11, computational results of terminal voltage waveform by the CDA and PDIM are in good agreement. However, CDA is an implicit method, and the Newton-Raphson formula must be used to calculate the voltage waveform of propagative transmission line. During the calculation, the Newton-Raphson formula solves the nonlinear algebraic equations with two iterations. This process is time-consuming when using small step simulations. The simulation efficiency of the two methods is compared in Table 2. The simulation platform is MATLAB R2012a. The tablet PC processor is AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx 2.10 GHz. The tablet PC uses a 64-bit operating system, and the capacity of RAM is 8 GB. The simulation starts from the zero initial state except for , which is given a small initial value so that the calculation can be performed. The total simulation ends at t = 60 ms. The simulation acceleration ratio of this example is defined as the ratio of the small time-step simulation time of the trapezoidal method (TR) to the large time-step simulation time of PDIM.
As can be observed in Table 2, the PDIM is significantly more efficient at handling nonlinear EMT simulations than the trapezoidal method. Obviously, as the simulation time step of PDIM increases, the corresponding acceleration ratio also changes significantly.
In Figure 12, the calculation result of PDIM using large time step is almost consistent with the simulation waveform of trapezoid method with small time step, which shows that PDIM has good numerical stability and high precision for nonlinear EMT models.
(a)
(b)
3.3. Lightning Overvoltage Calculation of Substation Bus
This case is a simulation example of lightning tower overvoltage calculation. Figure 13 is a simplified equivalent circuit diagram of lightning overvoltage calculation of substation bus. During lightning stroke, the lightning channel is simulated by resistor parallel ideal current source, and the resistance is the resistance of lightning channel; when lightning stroke hits the top of tower, the ideal switch closes; represents the impact grounding resistance of tower; the tower is modeled by lossless transmission line, whose wave impedance and wave velocity are and , respectively, and the transmission line length is = 50 m. Lightning current is simulated by double exponential wave. Its wavefront time and half peak time are and the peak value of lightning current is 100 kA.
The expression of lightning current is as follows:where , , and .
In this case, the electrical model of lightning overvoltage simulation for tower is established by using telegraph equation. After spatial interpolation and discretization using the fourth- and second-order interpolation formulas [25, 26], the following ordinary differential equations are obtained by taking the number of space segments N = 30:where constant coefficient matrices and are the input excitation sources of overvoltage at the top of lightning tower.
Equation (31) is solved by PDIM (with Chebyshev grid and s = 4) and the trapezoidal method, respectively. The simulation step of the two methods is h = 0.01 μs. The simulation results are shown in Figure 14. As shown in Figure 14, when , the voltage value of at the end of the line is about 984.5 kV and 951.6 kV for the head of the line. In Figures 15 and 16, with the time prolonging, it can be seen that the voltage and current of sending end and receiving end at the end of the transient process are almost the same. And the steady-state current values are almost 99.86 and 98.43 kA, which are the current values at the beginning and end of the transmission line, respectively, which shows the correctness of the simulation results of this case that its real steady-state current value of the lossless transmission line is near 99 kA.
In the calculation of lightning overvoltage, it can be seen from Figure 16 that, because of the fast-changing rate of double exponential lightning current, the two algorithms can only accurately simulate the changing waveform of lightning overvoltage by using smaller simulation steps. As can be observed in difference value of two methods in Figure 15, the simulation results of PDIM and the trapezoidal method are almost similar. This example shows that, for systems with very fast change frequency, PDIM is also competent.
4. Conclusions
Aimed at the simulation efficiency of EMT simulation for overvoltages of the high-voltage transmission line, an improved precise integration method based on DQM is proposed in this paper. The improved precise integration method inherits the characteristics of high precision and strong stability of the PIM and DQM. And PDIM improves the approximate calculation method of Duhamel integration term in the calculation of nonhomogeneous differential equations by the traditional PIM. Compared with the numerical results of CDA method or the trapezoidal method with small time step, the advantage of PDIM with larger time step is verified in the simulation efficiency of EMT simulation for high-voltage transmission lines.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the support from the National Natural Science Foundation of China (NSFC) through its Grant no. 52007103, Natural Science Foundation of Hubei Province through its Grant no. 2019CFB331 and Science, and Technology Project of State Grid Corporation of China through its Grant no. 5200-201956111A-0-0-00.