[b]1. Basic Method of Characteristic Difference[/b] The finite difference method is a major method for solving partial differential equations. Its key point is to discretize the problem in both spatial and temporal directions, approximate the derivative with the difference quotient, transform the differential equation into a difference equation, and then, starting from the initial conditions, proceed step-by-step over time to obtain the solution. When considering the influence of frequency characteristics and corona on transmission line parameters, the three-phase transmission line voltage matrix u = [uA(x,t), uB(x,t), uC(x,t)]T and the current matrix i = [iA(x,t), iB(x,t), iC(x,t)]T satisfy the following hyperbolic equation system: [img=247,76]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/31-2.jpg[/img] In the formula, x is a one-dimensional spatial coordinate, pointing from the beginning to the end of the line; t is the time coordinate; the reference direction of i is the same as x; uf is the additional voltage drop matrix simulating the influence of the frequency characteristics of line parameters per unit length, in the characteristic difference method uf = [ufA(x,t), ufB(x,t), ufC(x,t)]T; icor is the corona current matrix per unit length, icor = [icA(x,t), icB(x,t), icC(x,t)]T; C0 is the line capacitance matrix: [img=247,76]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/31-2.jpg[/img] L0 should satisfy the condition [img=286,76]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/31-3.jpg[/img] where c is the speed of light. The finite difference method requires that the line be uniformly segmented, and the uf and icor of each segment be moved to the nodes Z and Y formed by each segment, so that the actual line is discretized into a chain loop composed of several lossless line segments and nodes Z and Y. The voltage matrix and current matrix of adjacent nodes satisfy a certain relationship in time, that is, the forward traveling wave and the reverse traveling wave equation, as shown in Figure 1. In the figure, points p, d, and q are spatially adjacent nodes, h is the spatial step size, that is, the length of each segment after segmentation, Δt is the time step size, Δt＝h／c. Thus, the hyperbolic equations can be discretized into forward and reverse traveling wave equations: (ud-up) + z0(id-ip) = -h(uf+z0icor) (4) (ud-uq) - z0(id-iq) = h(uf-z0icor) (5) where z0 = [L0C0] - 1/2. Thus, if point d is the end point of the line, the forward or reverse traveling wave equation is solved by combining it with the end conditions; if point d is an internal node of the line, the above two equations are solved by combining them. In the calculation, the corona current can be calculated using the method provided in Reference 1. If the influence of corona is not considered, the matrix icor only needs to be set to zero. uf is the voltage drop matrix of the frequency characteristics of the line parameters, as shown in Figure 2. The parameters of each component in the figure can be determined using the method in Reference 2. [img=320,232]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/32-1.jpg[/img] [img=254,206]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/32-2.jpg[/img] [b]2 Correction for Instability in Numerical Calculation using the Characteristic Difference Method[/b] During the calculation, it was found that even small changes in the network data of the line parameters frequency characteristics could cause instability in the entire calculation, thus preventing the calculation from proceeding normally. Since uf in the forward and reverse traveling wave equations is the additional voltage drop matrix per unit length caused by the influence of the frequency characteristics of the line parameters, see Figure 2, [img=238,150]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/32-3.jpg[/img] After discretization using the three-point numerical differential formula, we have [img=353,150]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/32-4.jpg[/img] where t, t－Δt, and t－2Δt represent the values ​​of the variable at the current time, Δt ago, and 2Δt ago, respectively. In the finite difference method, explicit equations are generally used to simplify calculations. i1, i2, ..., in are taken as known currents i′1, i′2, ..., i′n before point Δt, or known currents i″1, i″2, ..., i″n at the same point in space before point Δt. This simplification is feasible in general calculations, but this method is highly sensitive to changes in the network parameters of the line frequency characteristics, easily causing numerical oscillations and leading to computational instability. From a physical perspective, since there are energy storage elements L1, L2, ..., Ln in the network, and currents i1, i2, ..., i″n flowing through them... The values ​​of current i1, i2, ..., in cannot change abruptly. Assigning them fixed values ​​forces these abrupt changes, leading to a maximum or even near-infinity in the derivative of the current, di/dt. This is the fundamental reason for its numerical instability. To address this instability, the following measures can be taken: treat the currents i1, i2, ..., in as variables, where i = i1 + i2 + ... + in is the value at the current calculation point at the current moment. This transforms the solution method into an implicit solution. Although the calculation process becomes more complex, the widespread use of high-performance computers makes the results more accurate and eliminates the problem of numerical oscillations. [b]3 Application of the Characteristic Difference Method in the Calculation of Overvoltage During Closing Operations in EHV Power Systems[/b] Due to limitations, extensive field measurements of overvoltage during closing operations in EHV power systems are not yet feasible. Compared to research methods utilizing physical simulation (such as TNA), numerical calculation methods are comparable in accuracy, comprehensiveness of considerations, and flexibility. Therefore, using numerical simulation to study the overvoltage during closing operations in EHV power systems, especially after eliminating the closing resistor, is a reasonable approach. Because the distribution of overvoltage during closing operations along the line differs significantly after eliminating the closing resistor, exhibiting a shape that is low at both ends and high in the middle, the characteristic difference method is a good choice for overvoltage calculation. Its greatest advantage is its ability to conveniently calculate the overvoltage distribution along the line. 3.1 Mathematical Models and Program Implementation of Each Component The main inductive components in a power system are the power supply equivalent inductance and shunt reactors. Theoretically, shunt reactors should be nonlinear inductors. However, since the frequency of switching overvoltages is higher than the power frequency, and the reactor saturation is not deep, in order to rigorously evaluate the surge arrester and simplify the calculation, the external characteristics of the shunt reactor can be assumed to be linear. Thus, when the reference direction of u [img=352,97]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/32-5.jpg[/img] in a three-phase system, a matrix is ​​formed, i.e. [img=321,40]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-1.jpg[/img] The volt-ampere characteristics of a three-phase switch are simulated using a time-controlled matrix: u＝RS i （10） Where R is the voltage across the switch; i is the current flowing through the switch; RS is [img=163,72]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-2.jpg[/img] If at time t＝t1, the A-phase switch is connected to the closing resistor Rc, and at time t＝t2, Rc is short-circuited, then [img=265,70]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-3.jpg[/img] Similarly, the switching simulation of the other two phases can be obtained. When t1＝t2, the case without the closing resistor can be simulated. Metal oxide surge arresters (MOAs) have no series or parallel gaps and are typical nonlinear components. MOAs can be simulated using piecewise linear simulation, which offers good numerical stability, as shown in Figure 3. [img=302,180]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-4.jpg[/img] The basic method for piecewise linear simulation is: U = Ek + RkI. The matrix form is U = Ek + RkI. When the number of segments reaches a certain level, the volt-ampere characteristics of the MOA can be simulated relatively accurately. Because accurate simulation of MOA is the foundation of the entire calculation, the piecewise linear simulation method used can be compared with the IEEE model [3] (see Figure 4). For MOA with the same volt-ampere characteristics, the accuracy of the piecewise linear model is verified. The current waveforms used are oblique wave and sine wave. The comparison results are shown in Table 1 (the energy dissipated by the IEEE model is the sum of all consumed active energy in Figure 4). The accuracy of the piecewise linear model can be seen from the table. [img=238,156]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-5.jpg[/img] 3.2 Digital simulation of the equivalent network calculation model and the connection of each component. Regarding the question of whether to use a complex network or a simple network, and whether to use a double circuit or a single circuit for calculation, EMTP was used for comparison calculation. The calculation results show that the closing overvoltage is the most serious when using a simple network with a single circuit. Therefore, the equivalent network calculation model adopts the structure of a simple network with a single circuit. [img=348,203]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-6.jpg[/img] As can be seen from the previous section, the mathematical model of each component of the power system has the following form: [img=333,121]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/33-7.jpg[/img] Suppose there are n components in the system: u1＝A1i1＋B1（15a） u2＝A2i2＋B2 （15b） … un＝Anin＋Bn（15n） When n components are connected in parallel, we have u1＝u2＝…＝un＝u（16） i1＋i2＋…＋in＝i（17） Therefore, we have A = [A1-1 + A-12 + ... + An-1]-1 (18) B = A1-1B1 + A-12 B2 + ... + An-1Bn (19) When n components are connected in series, we have i1 = i2 = ... = in = i (20) u1 + u2 + ... + un = u (21) Therefore, we have A = A1 + A2 + ... + An (22) B = B1 + B2 + ... + Bn (23) In this way, the entire iso-network can be transformed into a matrix equation, which can be solved using the characteristic difference method. 3.3 Determining the Initial Value of Single-Phase Reclosing: When the system trips due to a fault and recloses after the fault is cleared, the initial value of the system is not zero. This is because the faulty phase conductor to ground is equivalent to a capacitor, while the shunt reactor is an inductive element. The voltage and current in these energy storage elements cannot change abruptly. Therefore, the entire closing process is shown in Figure 5. The capacitor C and the inductor L both have initial values ​​u0 and i0 at t = 0. [img=319,160]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-1.jpg[/img] Strictly speaking, when calculating the overvoltage of single-phase reclosing, the entire transient process must be simulated, i.e., the steady state before the fault, the fault, the tripping, the fault clearing, and the reclosing. However, considering that the reclosing time is generally 0.3 s to 1.0 s, the initial value of the system is not zero. Since the transient process of the fault can be considered to have ended upon reclosing, calculations starting from the moment of reclosing will not introduce significant errors. The initial value for single-phase reclosing can be obtained by solving the steady-state value of the wiring when the system trips. Let phase A be the faulty phase, u01 be the recovery voltage at the beginning of the faulty phase, and u02 be the recovery voltage at the end of the faulty phase. From the fault boundary conditions: the current in phase A is zero, and the voltages at the breaks in phases B and C are zero, we have [img=323,76]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-2.jpg[/img] In the formula, the subscripts A, B, and C represent the three phases A, B, and C, respectively; the superscripts (0), (1), and (2) represent the zero sequence, positive sequence, and negative sequence, respectively; S = 1 represents the network structure shown in Figure 6. [img=331,117]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-4.jpg[/img][img=3 60,227]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-5.jpg[/img][align=left] [img=327,260]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-6.jpg[/img] [Image 1] [Image 2] [Image 3] [Image 4] [Image 5] [Image 6] [Image 7] A one-step equivalent, see Figure 7, then we have [Image 7] [Image 8] [Image 9] [Image 1] [Image 1] [Image 2] [Image 3] [Image 1] [Image 2] [Image 3] [Image 4] [Image 5] [Image 6] [Image 7] By combining the equations in this section and using the symmetrical component method, the recovery voltage u01 at the beginning of the circuit can be derived. [img=302,122]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/34-10.jpg[/img] 3.4 Implementation of the Calculation Program The characteristic difference method for eliminating numerical oscillations is used to study the technical conditions of overvoltage during closing in the EHV power system and the removal of the closing resistance of the line-side circuit breaker. The developed software must have the following functions: Based on the initial conditions of different power systems after calculation (power supply parameters, line parameters, whether closing resistors are installed, the arrangement of surge arresters, the arrangement of parallel reactors, etc.), the closing (single-phase reclosing) process is automatically simulated (closing time can be treated as uniform distribution and normal distribution respectively). The calculations include: phase-to-ground overvoltage, phase-to-phase overvoltage, maximum current flowing through each phase surge arrester, maximum half-wave energy, and cumulative energy absorbed in one closing operation; maximum phase-to-ground overvoltage, maximum phase-to-phase overvoltage, 2% phase-to-ground overvoltage, 2% phase-to-phase overvoltage, 50% phase-to-ground overvoltage and standard deviation, 50% phase-to-phase overvoltage and standard deviation, maximum energy dissipated by surge arresters, average energy dissipated by surge arresters and standard deviation, and statistical failure rate of line insulation (including phase-to-ground and phase-to-phase) in 120 operations. Figure 8 shows the PAD (Problem Analysis Diagram) of the software introduced in this paper. The software was tested and verified against Network 3 provided by CIGRE and compared with the results of TNA actual measurement (provided by China Electric Power Research Institute). The comparison results show that the calculation results of the software are reliable. [img=582,574]http://zszl.cepee.com/cepee_kjlw_pic/files/wx/dwjs/2000-10/35-1.jpg[/img] 4 Conclusion The characteristic difference method after adopting the implicit solution method can eliminate numerical oscillation. The software developed using this method can calculate the overvoltage of closing (single-phase reclosing), MOA dissipation energy and the statistical failure rate of line insulation. After comparison and verification, the calculation results of the developed software are reliable and can meet the needs of engineering calculation and engineering design. [b]References[/b] [1] Min Baomin. Corona model of transmission line and its application in calculation of switching overvoltage [C]. Xi'an: Xi'an Jiaotong University, 1985. [2] Wang Qinghai. Using zinc oxide surge arresters to limit phase-to-phase overvoltage during closing operation [C]. Xi'an: Xi'an Jiaotong University, 1987. [3] Shi Wei. Calculation of overvoltage in power systems [M]. Xi'an: Xi'an Jiaotong University Press, 1988.