# A Time-Domain Explicit Integration Algorithm for Fast Overvoltage Computation of High-Voltage Transmission Line

Authors: Lei Zhang; Zhongyi Xu; Jing Ye
Journal: Mathematical Problems in Engineering
 (2020)

---

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


---

## Body

## 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 Section2, 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:(1)x˙t=fxt,t∈t0,T,xt0=x0,where t represents time variable; ft,xt is a function of time t and state variable x; x0 is the value of state variable at initial time t0.Formula (1) can always be rewritten into ordinary differential equations consisting of homogeneous and nonhomogeneous parts by a certain transformation, and(2)x˙=Hx+gx,t,where H∈Rn×n 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:(3)x=expH⋅tx0.In the single-step integration, equation (3) has the following discrete recursive calculation format:(4)xk+1=T×xk,k=0,1,2,….The precise integration method (PIM) [14] proposed by Zhong uses additive theorem and incremental storage to calculate exponential matrix T=exph×H:(5)T=expH⋅h=expH⋅hmm,where m=2N (where N=15). For very small η=h/m, we can use Taylor series to expand the exponential matrix T in (5) into M terms (where M = 4); then, there are(6)expH⋅η=I+Ta,Ta=Hη+Hη22!+Hη33!+Hη44!.Substituting formula (6) into formula (5) yields(7)T=I+Tam=I+Ta2N−1⋅I+Ta2N−1.According to the matrix multiplication, the following recursive operations are performed iteratively:(8)Ta,i=2∗Ta,i−1+Ta,i−1∗Ta,i−1,i=1,2,…,N,where Ta,0=Ta. At the end of the cycle,(9)T=I+Ta,N.In the single integration steptk,tk+1, the solution of the nonhomogeneous dynamic equation (2) can be expressed as follows:(10)xk+1=Txk+∫tktk+1eHtk+1−τgx,τdτ.In formula (10), the exponential matrix T=exph×H 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 z˙=wt,z, the s-th formula of the DQ scheme with s-stage s-order can be expressed as follows:(11)zk+1=zk+h∑j=1sbjwtk+cjh,z˜j,where h is the time step; cj is the grid point; bj are the integration coefficients related to the grid points of DQM.The DQM uses uniform grid, and the formula of Duhamel integration termx¯k+1 in equation (10) is as follows:(12)x¯k+1=h∑i=1sbieHtk+1−t˜igx˜k+i/s,t˜i,i∈1,s,where t˜i=tk+i×h/s.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 s=4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x˜k+i/si=1,2,…,s as follows:(13)x˜k+1/4=xk+h4×6S1+2S2+2S3+S4,where S1=Hxk+gtk,xk;(14)S2=Hxk+h4×2S1+gtk+h4×2,xk+h4×2S1,S3=Hxk+h4×2S2+gtk+h4×2,xk+h4×2S2,S4=Hxk+h4S3+gtk+h4,xk+h4S3.The other values ofx˜k+1/2, x˜k+3/4, and x˜k+1 are calculated in turn according to formula (13), and the approximate value of Duhamel integration term x¯k+1 in equation (12) is obtained as follows:(15)x¯k+1=hb1K1+b2K2+b3K3+b4K4,where K1=T1gx˜k+1/4,t˜1; K2=T2gx˜k+1/2,t˜2; K3=T3gx˜k+3/4,t˜3; K4=gx˜k+1,t˜4; T1=eH×3h/4; T2=eH×h/2; T3=eH×h/4; there is(16)T2=T3×T3,T1=T2×T3,T=T1×T3.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 x¯k+1 of the Duhamel integration term, the approximate value of xk+1 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(17)ck=121−cos2k−12s−1π,k∈1,s−1,c0=0,cs=1.Similarly, whens=4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x˜k+cii=1,2,…,s as follows:(18)x˜k+ci=x˜k+ci−1+cih6S1+2S2+2S3+S4.In equation (18), x˜k+c0=xk;(19)S1=Hx˜k+ci−1+gtk+ci−1h,x˜k+ci−1,S2=Hx˜k+ci−1+cih2S1+gtk+ci−1h+cih2,x˜k+ci−1+cih2S1,S3=Hx˜k+ci−1+cih2S2+g˜tk+ci−1h+cih2,x˜k+ci−1+cih2S2,S4=Hx˜k+ci−1+cihS3+gtk+ci−1h+cih,x˜k+ci−1+cihS3.Then the approximationx¯k+1 of the Duhamel integration term is as follows:(20)x¯k+1=hb˜1K˜1+b˜2K˜2+b˜3K˜3+b˜4K˜4,where K˜1=T1gx˜k+c1,t^1; K˜2=T2gx˜k+c2,t˜2; K˜3=T3gx˜k+c3,t3; K˜4=gx˜k+1,t˜4, where T1=eH1−c1h; T2=eH×h1−c2; T3=eH×h1−c3; T=T1×T3.The approximation ofxk+1 can be obtained by substituting x¯k+1 into equation (10). Obviously, using nonuniform grids needs to calculate s−1 exponential matrices, which requires more extra calculation than using uniform grid.


## 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 z˙=wt,z, the s-th formula of the DQ scheme with s-stage s-order can be expressed as follows:(11)zk+1=zk+h∑j=1sbjwtk+cjh,z˜j,where h is the time step; cj is the grid point; bj are the integration coefficients related to the grid points of DQM.The DQM uses uniform grid, and the formula of Duhamel integration termx¯k+1 in equation (10) is as follows:(12)x¯k+1=h∑i=1sbieHtk+1−t˜igx˜k+i/s,t˜i,i∈1,s,where t˜i=tk+i×h/s.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 s=4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x˜k+i/si=1,2,…,s as follows:(13)x˜k+1/4=xk+h4×6S1+2S2+2S3+S4,where S1=Hxk+gtk,xk;(14)S2=Hxk+h4×2S1+gtk+h4×2,xk+h4×2S1,S3=Hxk+h4×2S2+gtk+h4×2,xk+h4×2S2,S4=Hxk+h4S3+gtk+h4,xk+h4S3.The other values ofx˜k+1/2, x˜k+3/4, and x˜k+1 are calculated in turn according to formula (13), and the approximate value of Duhamel integration term x¯k+1 in equation (12) is obtained as follows:(15)x¯k+1=hb1K1+b2K2+b3K3+b4K4,where K1=T1gx˜k+1/4,t˜1; K2=T2gx˜k+1/2,t˜2; K3=T3gx˜k+3/4,t˜3; K4=gx˜k+1,t˜4; T1=eH×3h/4; T2=eH×h/2; T3=eH×h/4; there is(16)T2=T3×T3,T1=T2×T3,T=T1×T3.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 x¯k+1 of the Duhamel integration term, the approximate value of xk+1 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(17)ck=121−cos2k−12s−1π,k∈1,s−1,c0=0,cs=1.Similarly, whens=4, the four-order explicit Runge-Kutta method is used to calculate the estimated value x˜k+cii=1,2,…,s as follows:(18)x˜k+ci=x˜k+ci−1+cih6S1+2S2+2S3+S4.In equation (18), x˜k+c0=xk;(19)S1=Hx˜k+ci−1+gtk+ci−1h,x˜k+ci−1,S2=Hx˜k+ci−1+cih2S1+gtk+ci−1h+cih2,x˜k+ci−1+cih2S1,S3=Hx˜k+ci−1+cih2S2+g˜tk+ci−1h+cih2,x˜k+ci−1+cih2S2,S4=Hx˜k+ci−1+cihS3+gtk+ci−1h+cih,x˜k+ci−1+cihS3.Then the approximationx¯k+1 of the Duhamel integration term is as follows:(20)x¯k+1=hb˜1K˜1+b˜2K˜2+b˜3K˜3+b˜4K˜4,where K˜1=T1gx˜k+c1,t^1; K˜2=T2gx˜k+c2,t˜2; K˜3=T3gx˜k+c3,t3; K˜4=gx˜k+1,t˜4, where T1=eH1−c1h; T2=eH×h1−c2; T3=eH×h1−c3; T=T1×T3.The approximation ofxk+1 can be obtained by substituting x¯k+1 into equation (10). Obviously, using nonuniform grids needs to calculate s−1 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 Figure1, this circuit represents the equivalent schematic diagram of a 515 kV bus switching AC filter in a substation. Ls is the equivalent inductance of the AC grid; Rs represents the equivalent resistance of the AC grid; es 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.Figure 1
The equivalent circuit of a 515 kV bus switching AC filter.When the circuit breakerK 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(21)LsCd2ucdt2+RsCducdtuc+uc=es,(22)dⅆtiCuC=−RsLs−1Ls1C0iCuC+esLs0.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 as(23)uCt=K1sinωαt+K2cosωαte−αt+Ucmsinωt+φc,where est=Umsinωt+φs; ω=100πrad/s;(24)Ucm=UmRs2+ωL−1/ωC21ωC,φc=ϕs−π2−arctanωLs−1/ωCRs,K1=−αUcmsinφc−ωUcmcosφcωα,K2=−Ucmsinφc.In Figure1, the decay constant of resistance, inductance, and capacitor series circuit is α=Rs/2Ls; resonance angular frequency of this circuit is ω0=1/LsC; and natural angular frequency of this circuit is ωα=ω02−α2 [22, 23].In this simulation, the lumped electrical parameters areRs=1.59Ω, Ls=35.37mH, and C=3.124μF. 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.Figure 2
Computational result of capacitor voltage by precise differential integration method.Figure 3
The curve of computational difference values by these three methods.Figure 4
The curve of computational difference values by these three methods.As shown in Figures2–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, rs=10Ω, and equivalent inductor, ls=20mH. The piecewise linear current-flux curve of lnl is given in Figure 6. The peak voltage of the sinusoidal voltage source is 44.52kV and est=44.52sin120πtkV. 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.Figure 5
The equivalent circuit diagram of a saturated reactor charged by an AC grid.Figure 6
A nonlinear current-ﬂux curve of the saturated reactor.Figure 7
Calculated results of the saturated reactor obtained by PDIM and TR. (a) Local snapshot results of the saturated reactor voltage calculated by PDIM and TR. (b) Calculated result of the saturated reactor voltage obtained by PDIM. (c) Calculated result of the saturated reactor current obtained by PDIM.
(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 Figure8. In Figure 8, et is the excitation AC voltage source; R˜s and L˜s 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 R0, inductance L0, and capacitance C0. The nonlinear load is composed of load resistance RL and the nonlinear inductor LL in parallel connection, where the relationship between the magnetic link ϕL and current iL of the nonlinear inductor is ϕL=atanhb×iL; the electrical parameters in this case are shown in Table 1.Figure 8
Schematic of equivalent circuit of single-phase transmission line with nonlinear inductance load.Table 1
System parameters of example 3.2.
Electrical parametersParameter valuesR˜s2ΩL˜s0.06 HR00.07Ω/kmL02.08×10−3H/kmC012×10−9F/kmRL96Ωa8.40×102V·sb5.95×10−3A−1The electrical model used to describe the EMT process of the high-voltage transmission line shown in Figure8 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:(25)r=R0LM,l=L0LM,c=C0LM.Figure 9
Equivalent model of Π-type lumped parameter circuits.As shown in Figure9, 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]:(26)LSdi0dt+RSi0+u1=et,ldimtdt=umt−um+1t−rimt,m∈1,M,cdumtdt=im−1t−imt,iM+1=uM+1RL+iL,m∈1,M+1,LLdiLdt=uM+1,LL=atanhbiLiL.Considering the boundary conditions at the head and end of the transmission line in Figure9, equation (26) is arranged into the following matrix form:(27)y˙t=Ayt+μt,where A∈R2M+3×2M+3 is the constant sparse matrix; μt is 2M+3 sparse column vector, which is the excitation source of EMT simulation for the transmission line;(28)yt=i0i1…iMu1…uM+1iLT,μt=etLs0…00…0uM+1LL−iLT, (29)In Figure8, the number of interval segments is M=30 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 h=1.0×10−6s and PDIM uses a larger time step h^=10h to calculate transient voltage of the transmission line. And the sending and receiving end voltage are shown in Figures 10 and 11.Figure 10
Computational results of sending terminal voltage by improved PIM with uniform grid.Figure 11
Computational results of receiving terminal voltage by improved PIM with uniform grid.As shown in Figures10 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 iL, 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.Table 2
Comparison of calculation time between the two algorithms.
AlgorithmTime step (h)Time consumed (s)Speedup ratioTRh=1.0×10−6s46.551PDIMh=1.0×10−5s0.7562PDIMh=2.0×10−5s0.36129PDIMh=3.0×10−5s0.26178As can be observed in Table2, 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 Figure12, 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.Figure 12
Computational results of receiving terminal voltage by improved PIM with uniform grid using bigger time steps. (a) PDIM withh=2.0×10−5s. (b) PDIM with h=3.0×10−5s.
(a)(b)


### 3.3. Lightning Overvoltage Calculation of Substation Bus

This case is a simulation example of lightning tower overvoltage calculation. Figure13 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 re is the resistance of lightning channel; when lightning stroke hits the top of tower, the ideal switch K¯ closes; Rch=10Ω represents the impact grounding resistance of tower; the tower is modeled by lossless transmission line, whose wave impedance and wave velocity are 100Ω and 2.7×108m/s, respectively, and the transmission line length is L1 = 50 m. Lightning current is simulated by double exponential wave. Its wavefront time and half peak time are 2.6/50μs and the peak value of lightning current iLt is 100 kA.Figure 13
A simplified model for lightning overvoltage calculation.The expression of lightning current is as follows:(30)iLt=aILe−αt−e−βt,where a=1.058, α=1.5×10−2μs−1, and β=1.86μs−1.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:(31)x.=Hx+δt,where constant coefficient matrices H∈R61×61 and δt 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 t=4μs, the voltage value of ut 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.Figure 14
Computational results of voltage waveform by PDIM with Chebyshev nodes.Figure 15
Receiving end voltage difference waveform of PDIM with Chebyshev nodes and TR.Figure 16
Computational results of current waveform by PDIM with Chebyshev nodes.In the calculation of lightning overvoltage, it can be seen from Figure16 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.


## 3.1. Accuracy Test for Improved Precise Integration Method

As shown in Figure1, this circuit represents the equivalent schematic diagram of a 515 kV bus switching AC filter in a substation. Ls is the equivalent inductance of the AC grid; Rs represents the equivalent resistance of the AC grid; es 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.Figure 1
The equivalent circuit of a 515 kV bus switching AC filter.When the circuit breakerK 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(21)LsCd2ucdt2+RsCducdtuc+uc=es,(22)dⅆtiCuC=−RsLs−1Ls1C0iCuC+esLs0.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 as(23)uCt=K1sinωαt+K2cosωαte−αt+Ucmsinωt+φc,where est=Umsinωt+φs; ω=100πrad/s;(24)Ucm=UmRs2+ωL−1/ωC21ωC,φc=ϕs−π2−arctanωLs−1/ωCRs,K1=−αUcmsinφc−ωUcmcosφcωα,K2=−Ucmsinφc.In Figure1, the decay constant of resistance, inductance, and capacitor series circuit is α=Rs/2Ls; resonance angular frequency of this circuit is ω0=1/LsC; and natural angular frequency of this circuit is ωα=ω02−α2 [22, 23].In this simulation, the lumped electrical parameters areRs=1.59Ω, Ls=35.37mH, and C=3.124μF. 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.Figure 2
Computational result of capacitor voltage by precise differential integration method.Figure 3
The curve of computational difference values by these three methods.Figure 4
The curve of computational difference values by these three methods.As shown in Figures2–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, rs=10Ω, and equivalent inductor, ls=20mH. The piecewise linear current-flux curve of lnl is given in Figure 6. The peak voltage of the sinusoidal voltage source is 44.52kV and est=44.52sin120πtkV. 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.Figure 5
The equivalent circuit diagram of a saturated reactor charged by an AC grid.Figure 6
A nonlinear current-ﬂux curve of the saturated reactor.Figure 7
Calculated results of the saturated reactor obtained by PDIM and TR. (a) Local snapshot results of the saturated reactor voltage calculated by PDIM and TR. (b) Calculated result of the saturated reactor voltage obtained by PDIM. (c) Calculated result of the saturated reactor current obtained by PDIM.
(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 Figure8. In Figure 8, et is the excitation AC voltage source; R˜s and L˜s 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 R0, inductance L0, and capacitance C0. The nonlinear load is composed of load resistance RL and the nonlinear inductor LL in parallel connection, where the relationship between the magnetic link ϕL and current iL of the nonlinear inductor is ϕL=atanhb×iL; the electrical parameters in this case are shown in Table 1.Figure 8
Schematic of equivalent circuit of single-phase transmission line with nonlinear inductance load.Table 1
System parameters of example 3.2.
Electrical parametersParameter valuesR˜s2ΩL˜s0.06 HR00.07Ω/kmL02.08×10−3H/kmC012×10−9F/kmRL96Ωa8.40×102V·sb5.95×10−3A−1The electrical model used to describe the EMT process of the high-voltage transmission line shown in Figure8 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:(25)r=R0LM,l=L0LM,c=C0LM.Figure 9
Equivalent model of Π-type lumped parameter circuits.As shown in Figure9, 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]:(26)LSdi0dt+RSi0+u1=et,ldimtdt=umt−um+1t−rimt,m∈1,M,cdumtdt=im−1t−imt,iM+1=uM+1RL+iL,m∈1,M+1,LLdiLdt=uM+1,LL=atanhbiLiL.Considering the boundary conditions at the head and end of the transmission line in Figure9, equation (26) is arranged into the following matrix form:(27)y˙t=Ayt+μt,where A∈R2M+3×2M+3 is the constant sparse matrix; μt is 2M+3 sparse column vector, which is the excitation source of EMT simulation for the transmission line;(28)yt=i0i1…iMu1…uM+1iLT,μt=etLs0…00…0uM+1LL−iLT, (29)In Figure8, the number of interval segments is M=30 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 h=1.0×10−6s and PDIM uses a larger time step h^=10h to calculate transient voltage of the transmission line. And the sending and receiving end voltage are shown in Figures 10 and 11.Figure 10
Computational results of sending terminal voltage by improved PIM with uniform grid.Figure 11
Computational results of receiving terminal voltage by improved PIM with uniform grid.As shown in Figures10 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 iL, 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.Table 2
Comparison of calculation time between the two algorithms.
AlgorithmTime step (h)Time consumed (s)Speedup ratioTRh=1.0×10−6s46.551PDIMh=1.0×10−5s0.7562PDIMh=2.0×10−5s0.36129PDIMh=3.0×10−5s0.26178As can be observed in Table2, 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 Figure12, 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.Figure 12
Computational results of receiving terminal voltage by improved PIM with uniform grid using bigger time steps. (a) PDIM withh=2.0×10−5s. (b) PDIM with h=3.0×10−5s.
(a)(b)


## 3.3. Lightning Overvoltage Calculation of Substation Bus

This case is a simulation example of lightning tower overvoltage calculation. Figure13 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 re is the resistance of lightning channel; when lightning stroke hits the top of tower, the ideal switch K¯ closes; Rch=10Ω represents the impact grounding resistance of tower; the tower is modeled by lossless transmission line, whose wave impedance and wave velocity are 100Ω and 2.7×108m/s, respectively, and the transmission line length is L1 = 50 m. Lightning current is simulated by double exponential wave. Its wavefront time and half peak time are 2.6/50μs and the peak value of lightning current iLt is 100 kA.Figure 13
A simplified model for lightning overvoltage calculation.The expression of lightning current is as follows:(30)iLt=aILe−αt−e−βt,where a=1.058, α=1.5×10−2μs−1, and β=1.86μs−1.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:(31)x.=Hx+δt,where constant coefficient matrices H∈R61×61 and δt 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 t=4μs, the voltage value of ut 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.Figure 14
Computational results of voltage waveform by PDIM with Chebyshev nodes.Figure 15
Receiving end voltage difference waveform of PDIM with Chebyshev nodes and TR.Figure 16
Computational results of current waveform by PDIM with Chebyshev nodes.In the calculation of lightning overvoltage, it can be seen from Figure16 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.



---

*Source: 8437617-2020-10-12.xml*