Three-Dimensional FDTD-Based Simulation of Lightning-Induced Surges in Secondary Circuits With Shielded Control Cables Over Grounding Grids in Substations

Nowadays, plenty of sensitive electronic devices are installed in secondary circuits in substations. When lightning strikes a substation and transmission line, lightning currents flow into the grounding structure of the substation through shield wires and lightning surge arresters in the primary circuit, and lightning overvoltages are induced on control cables in the secondary circuit owing to the effect of transient ground potential rises of the grounding structure and electromagnetic coupling. Shielded control cables are one of the countermeasures against lightning overvoltage induced on control cables; thus, it is useful to evaluate the effectiveness of shielded control cables for designing protection measures properly. Recently, the finite-difference time-domain (FDTD) method, which is one of the full-wave numerical approaches, has been commonly used for simulating electromagnetic transient phenomena in grounding structures, such as the grounding grids of substations. In this work, first, we studied the modeling of tape-type shielded control cables used in substations for FDTD-based electromagnetic transient analyses. Then, using the control cable model, we simulated a test platform of primary and secondary circuits in a substation and calculated the transient response of a grounding grid and voltages induced on the control cable when a lightning impulse current was injected into the grounding grid. We then compared the calculated results with the measured waveforms for validation.

In this work, first, we studied the modeling of tape-type shielded control cables used in substations for FDTD-based electromagnetic transient analyses. Then, using the control cable model, we simulated a test platform of primary and secondary circuits in a substation and calculated the transient response of a grounding grid and voltages induced on the control cable when a lightning impulse current was injected into the grounding grid. We then compared the calculated results with the measured waveforms for validation.

I. INTRODUCTION
E LECTROMAGNETIC disturbances are induced in secondary circuits in substations and power plants owing to lightning and switching transients arising during the operation of disconnectors and circuit breakers [1]. Nowadays, plenty of sensitive electronic devices, which are more vulnerable to electromagnetic transients than conventional analog-type devices, are installed in secondary circuits. To protect secondary circuits properly, electromagnetic disturbances and the effectiveness of countermeasures have been studied through measurements in actual substations or test platforms (e.g., [2], [3], [4], [5], [6], [7], [8], and [9]). On the other hand, numerical simulations are also practical tools for predicting electromagnetic transient phenomena and evaluating the effectiveness of countermeasures. Traditionally, circuit-theory-or transmission-line (TL)-theory-based simulation techniques have been studied for the electromagnetic transient analysis of secondary circuits [7], [8], [9], [10], [11], [12], [13], [14], [15], and simulation techniques for modeling high-frequency voltage and current transformers (VT and CT, respectively) have been proposed based on circuit elements and the use of the black-box technique [16], [17], [18], [19]. In Japan, a survey was carried out to investigate electromagnetic disturbances of panel-type equipment occurring in power plants, substations, and other electric power facilities for about ten years up to 1999, which revealed that 72% of disturbances were presumed to be caused by lightning [20], [21].
Recently, the application of full-wave numerical approaches, such as the finite-difference time-domain (FDTD) method [22], the method of moments (MoM) [23], the hybrid electromagnetic model [24], the finite-element method (FEM) [25], and so forth, has been useful for analyzing electromagnetic transient phenomena in 3-D structures and grounding systems since they do not require the assumption of transverse electromagnetic mode (TEM) propagation [26]. In [27], the responses of a grounding grid and coaxial shielded cables to AC currents and lightning strikes were studied using a MoM-based simulation technique coupled with circuit equations. In [28], magnetic fields around the metallic sheath of a braided cable and the shielding effect of a metallic channel against magnetic fields to protect circuits from electromagnetic interferences caused by disconnectors and circuit breakers were studied by FEM simulations.
The FDTD method has a great advantage compared with other methods in terms of its capabilities to model inhomogeneous electrical parameters in soil and nonflat ground surfaces in a straightforward manner and to incorporate nonlinear phenomena such as the characteristics of surge arresters [29], [30] and the occurrence of back-flashover at arcing horns of TL towers [31], [32]. The development of new simulation techniques to apply the FDTD method to electromagnetic transient analyses and the proliferation of high-performance computers have enhanced the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ efficacy of the FDTD method for lightning transient simulations, and the FDTD method has been applied to practical transient analyses of, for example, power plants, substations, TLs, distribution lines, and buildings [33], [34], [35], [36].
In previous studies regarding the application of the FDTD method to the lightning transient analysis of substations, the main focus was put on the analysis of lightning transients entering substations owing to lightning strikes and of the transient responses of grounding structures (e.g., [31], [37], [38], [39], and [40]). Such lightning transients propagate in the primary circuits and the grounding structure of substations, and then electromagnetic disturbances are induced in the secondary circuits through the effect of surge transitions from the primary to secondary circuits at instrument transformers such as VTs and CTs, the effect of transient ground potential rises (GPRs) of the grounding structure and electromagnetic coupling due to the lightning currents propagating through the grounding structure, the effect of direct coupling from the lightning transients propagating through the primary circuits, and so forth [20]. In [41], voltages induced in secondary circuits were calculated using the FDTD method when a transient voltage was applied to primary circuits in a test platform of primary and secondary circuits in a substation, whereas the frequency-dependent effect of a surge transition at VT and CT was taken into account on the basis of the black-box technique, and the FDTD simulations were validated by comparison with measurements.
In this study, we focused on voltages induced in secondary circuits owing to the effect of the GPRs of grounding structures. First, for FDTD-based lightning transient simulations, we studied the modeling of tape-type shielded control cables, which are one of the countermeasures against electromagnetic disturbances in secondary circuits and are commonly installed in extra-high-voltage and high-voltage substations in Japan. Second, the applicability of FDTD simulations of the transient analysis of voltages induced on shielded control cables using the test platform of primary and secondary circuits, which was mainly composed of a grounding grid, gas-insulated switchgear (GIS) model, control cable, and digital-type protection relay, when a lightning impulse current was injected into the grounding grid, was validated by comparing the simulated results with measured waveforms.

A. Modeling of Shielded Control Cables for FDTD-Based Lightning Transient Simulations
Several techniques have been developed for representing coaxial cables [42], [43], [44], pipe-type power cables [45], and coaxial cables with multiple conductive layers [46] for FDTD-based electromagnetic transient simulations, and these techniques have been applied to the electromagnetic transient analysis of power cables in a substation [37], power cables in a distribution line [45], [46], and a TL tower with a power cable [47]. In this study, the metallic sheath of a shielded control cable is simulated by a lossy thin wire using the technique proposed in [42] and the 3-D FDTD method, whereas electromagnetic transient phenomena inside the metallic sheath are solved by applying the TL theory.
To simulate the frequency-dependent conductor internal impedance of the metallic sheath, a series of lumped impedances are placed along a thin wire to represent the metallic sheath in the 3-D FDTD method. The expressions to update the electric fields on the lumped impedances to take into account the conductor internal impedance of a tubular conductor representing the metallic sheath are given in [42]. The radius of the thin wire is represented using the technique presented in [48], where the relative permittivity and permeability are modified in the updated four radial electric fields and four rotational magnetic fields adjacent to the thin wire.
In this study, we assume that an electrostatic coupling through shields is less important and that the surface transfer impedance is more dominant than the surface transfer admittance. Voltages are assumed to be induced on the N cores of the shielded control cable owing to the effect of the surface transfer impedance of the metallic sheath. By assuming that weak coupling and the propagation of transients inside the metallic sheath are the TEM mode, voltages and currents inside the metallic sheath are expressed on the basis of TL theory as follows [49]. [50], [51], [52]: where V(l, t) and I(l, t) are the vectors of voltages and currents along the cores of the control cable, respectively. L and C denote the per-unit-length (p.u.l.) inductance (N × N) and capacitance matrices (N × N), respectively. The second term in (1) represents the frequency-dependent loss due to the finite conductivities of the cores. Z t (N vector) and I s (l, t) are, respectively, the surface transfer impedance and sheath current. As will be presented in Section II-B, the surface transfer impedance of the shield of a control cable shows frequencydependent characteristics, and thus the right-hand side of (1) is calculated as follows. The relationship between input i(s) and output signals v(s) is expressed using a transfer function H(s) as follows: Assuming that v(s) tends to zero when frequency f → Ý, the transfer function can be approximated using the following expression: where p k denotes the kth pole and R k denotes residue of p k . Here, by replacing H(s) and i(s) with the surface transfer impedance Z t (s) and the sheath current I s (s), respectively, and by applying the trapezoidal rule for a time discretization Δt c , we can obtain the following expressions to calculate the right-hand side of (1), i.e., V s (l, t) = Z t * I s (l, t) [53]: where K r and K c are the numbers of the identified real poles and complex-conjugate pole pairs, respectively. p k and p k denote the real poles and complex poles with positive imaginary parts, respectively, whereas the residues of p k and p k are denoted by R k and R k , respectively. Here, p k , p k , R k , and R k are identified using the XFIT software based on a frequency partitioning fitting [53], [54] in XTAP (expandable transient analysis program), which is an electromagnetic transient simulation program developed by the Central Research Institute of Electric Power Industry (CRIEPI) based on circuit theory. In this study, as mentioned above, the metallic sheath is represented by a thin wire in the 3-D FDTD method, and the sheath current I s can be obtained from the magnetic fields around the thin wire on the basis of Ampere's law. Telegrapher's equations (1) and (2) are solved by applying a 1-D FDTD method with a time discretization Δt c and a cell size Δl. Note that it is not necessary to set Δt c and Δl to the same values as Δt and Δs, which, respectively, denote the time discretization and cell size used in the 3-D FDTD method. Δt c and Δt are required to satisfy the Courant stability condition [55]. In this study, we obtain I s (l, t) in (1) from the calculated result of I s in the 3-D FDTD method using the first-order basis function in space and the second-order basis function in time. As an example, Fig. 1 shows a procedure to update the electric and magnetic fields (E and H) in the 3-D FDTD method and the voltage and current (V and I) in the 1-D FDTD method, when Δt c is set to be a half value of Δt. The superscript n of the variables represents the time nΔt. In the updating procedure in Fig. 1, the sheath current I s calculated by the 3-D FDTD method is required to update the current I in the 1-D FDTD method. For example, the update of I n−1/2 needs I s n , and thus in this study, it is obtained using the following expression: Substituting ξ = 1/2, which corresponds to t = nΔt, into the above given equations gives I s n . Although actually, control cables used in substations are composed of twisted wires, the cores of the shielded control cable are treated approximately as parallel-wire pairs in distributed-parameter line models on the basis of the assumption that we can simulate electromagnetic transient phenomena inside the metallic sheaths of the shielded control cables using uniform distributed-parameter line models representing parallelwire pairs, where the inductance and capacitance matrices are approximated to be uniform, because the geometrical position of the cores and metallic sheaths of shielded control cables in a cross section is almost the same.
In this study, the FDTD simulations were carried out using an electromagnetic transient simulation code developed by CRIEPI on the basis of the 3-D FDTD method, which is called VSTL REV (Virtual Surge Test Lab. Restructured and Extended Version) [56], [57], and it can be run on a high-performance parallel computer with graphics processing units (GPUs) using compute unified device architecture (CUDA), [58].

B. Calculation of Inductance and Capacitance Matrices of Shielded Control Cables
The p.u.l. inductance matrix L in (1) is computed using a conductor-subdivision-based method for calculating the series impedance matrix of arbitrary cross-sectional conductors [59], [60], [61]. In the conductor-subdivision-based method, the cross sections of the cores and sheath of the shielded cable are divided into square subconductors in this study, and the inductance matrix can be obtained from the calculated self-inductance and mutual inductance of the subconductors on the basis of geometric mean distances on the assumption that the current density is uniform in each subconductor. On the other hand, the p.u.l. capacitance matrix C in (2) is calculated using an electrostatic field analysis code developed by CRIEPI on the basis of the 2-D surface charge method [62], [63], which is one of the boundary-dividing methods, such as the boundary element method (BEM), and is also called the indirect BEM. The analysis code can simulate conductors and dielectrics, the surfaces of which are divided into curved line elements. The shapes of curved line elements are represented by the third-order shape function [64], whereas the surface charge density on each element is represented by the fourth-order basis function [65].

C. Measurement of Surface Transfer Function
Using an experimental setup, shown in Fig. 2, we measured the surface transfer impedance of a shielded control cable. The shielded control cable was composed of a tape-type shield and two cores with a cross section of 3.5 mm 2 , which are hereafter referred to as cores #1 and #2, respectively. The cores had a twisted-wire structure. As shown in Fig. 2, the shielded control cable with a length of 5 m was placed horizontally at a height of 0.4 m over a large grounded copper plate. At one end of the control cable, to inject a current into the metallic sheath, a pulse generator was connected to the metallic sheath via an amplifier. At this end, both cores were connected to the metallic sheath. At the other end, the metallic sheath was connected to the copper plate via a resistor of 10 Ω. By injecting a current with a frequency of 10 kHz to 5 MHz with the pulse generator, we measured the sheath current I s (t) and the voltage V c (t) between core #1 or #2 and the metallic sheath.
From the measured results, we obtained the surface transfer impedance Z t (f) of the shielded control cable using the following expression: where F(·) and L denote the Fourier transform and the length of the cable, respectively. Fig. 3 shows the measured surface transfer impedances of the shielded control cable for cores #1 and #2. As shown in Fig. 3, the surface transfer impedances for cores #1 and #2 are very similar to each other in terms of amplitude and phase. For example, the differences of the amplitudes for cores #1 and #2 are less than 6.7%. The amplitudes of the surface transfer impedances are almost flat at frequencies of less than a few tens of kHz, but at higher frequencies, the amplitudes increase owing to the effect of the diffusion of the electromagnetic field through the metallic sheath.

D. Validation of Technique for Representing Shielded Control Cables
First, as shown in Fig. 4, using a calculation arrangement to simulate the experimental setup used for measuring the surface transfer impedance of the shielded control cable, we calculated the surface transfer impedance of the shielded control cable by the FDTD method with the technique for representing shielded control cables. The dimensions of the analysis space to be solved by the FDTD method were 15 m × 10 m × 5 m. The bottom surface of the analysis space was a perfectly conducting ground plane, whereas the other five surfaces were treated as absorbing boundaries on the basis of Liao's formulation of the second order [66]. The shielded control cable was modeled using the technique presented in the previous sections: The poles and residues in (12) and (13) in Section II-A to simulate the surface transfer impedances for cores #1 and #2 were obtained by applying the XFIT to the measured results shown in Fig. 3. A thin wire to represent the metallic sheath was placed at a height of 0.4 m. At one end of the thin wire, a current source to produce a Gaussian pulse was connected, where the conductors representing the cores in Telegrapher's equations were grounded. At the other end, a lumped-parameter resistance of 10 Ω was connected [67]. The conductivities of the thin wire representing the metallic sheath and the conductors representing the cores were set to be 5.4 × 10 7 S/m. Fig. 5(a) and (b) shows calculation models used for calculating the p.u.l. inductance and capacitance matrices in (1) and (2), respectively. In the calculation of the inductance matrix, as shown in Fig. 5(a), the cross sections of the cores and metallic sheath of the control cable were divided uniformly into 0.02-mm square subconductors, and the total number of the subconductors was 25 832. In the calculation of the capacitance matrix, as shown in Fig. 5(b), the surfaces of the conductors representing each of the core and metallic sheath were divided into 40 and 800 curved line elements, respectively. The vinyl  sheath of each core was simulated by a dielectric with a relative permittivity of 4.7, and the surface of each dielectric was divided into 800 curved line elements. The total number of variables used in the electrostatic field analysis was 9918.
Injecting a Gaussian-pulse current into the thin wire representing the metallic sheath, we calculated the surface transfer impedances for cores #1 and #2 using (16) from the calculated results of the sheath current I s and the voltages V c of each core. Fig. 6 shows the calculated and measured surface transfer impedances for core #1. For the sake of brevity, the calculated surface transfer impedance for core #2 is not shown here, but we confirmed that the calculated surface transfer impedances agree well with the measured results for both cores #1 and #2.
Second, we measured voltages induced on the cores when a current was injected into the metallic sheath of the shielded control cable. Fig. 7 shows an experimental setup used for measuring the voltages induced on the cores. Although this experimental setup was similar to the one shown in Fig. 2, the setup was changed as follows.
1) A resistor of 50 Ω was inserted between the two cores at each end. 2) At the end close to the pulse generator, the cores were disconnected from the metallic sheath. 3) At the other end, one of the cores was connected to the metallic sheath to simulate the condition that one of the cores is grounded in VT and CT circuits in secondary circuits in substations. By injecting a current with a frequency of 10 kHz to 5 MHz with the pulse generator, we measured the sheath current I s (t) and the induced voltage V r (t) across a resistor of 50 Ω. In addition, in the FDTD method, we updated the calculation arrangement in Fig. 4 to simulate the above-mentioned experimental setup used to measure the induced voltages.
Using the updated calculation arrangement, we injected a Gaussian-pulse current into the thin wire representing the metallic sheath of the control cable, and we calculated the sheath  current I s (t) and the induced voltage V r (t). From the measured and calculated results, we obtained the normalized induced voltages V i (f) using the following expression: Fig. 8 shows the calculated and measured normalized induced voltages. As shown in Fig. 8, although the calculated and measured induced voltages are in good agreement at frequencies higher than around 100 kHz in terms of amplitude and phase, the differences between the two results are larger at lower frequencies. The differences between the measured and calculated induced voltages at lower frequencies arise from the modeling error of the effect of the surface transfer impedances. For lower frequencies, such as below 1 kHz, the transient phenomena inside the metallic sheath of the shielded control cable can be approximately represented using a simple equivalent circuit shown in Fig. 9, where, for simplicity, the cores and sheath are treated as lossless conductors and the phases of the surface transfer impedances are assumed to be zero (see Fig. 3). As shown in Fig. 9(a), if the sources on cores #1 and #2 are assumed to be almost same due to the effect of the sheath current through the surface transfer impedances , then in such an almost balanced condition, the induced voltage across the resistance is quite small. This is the main reason why the amplitude of the induced voltage in the FDTD simulation gets smaller for lower frequencies in Fig. 8. As shown in Fig. 8, at lower frequencies less than some hundred kHz, the amplitude of the measured induced voltage is almost constant, and thus the induced voltage at these frequencies is presumed to be caused by a slightly unbalanced condition between the sources on cores #1 and #2 due to the effect of the surface transfer impedances, which are difficult to correctly obtain in the measurements generally, as shown in Fig. 9(b).
To incorporate such an unbalanced condition of the effect of the surface transfer impedances in a simple manner, here, we set the amplitude of the surface transfer impedance for core #2 (Z t2 ) using the following expression based on the measured surface transfer impedance for core #1 (Z t1 ), assuming that the phase of Z t2 is the same as that of Z t1 : Using the simple equivalent circuit shown in Fig. 9(b) and the measured induced voltage at a frequency of 1 kHz, we estimated an approximate value of k = 1.023, and we finally set k to be 1.03 at frequencies of 1-50 kHz and 1.0 at higher frequencies to obtain better results. Using the measured surface transfer impedance for core #1 and the surface transfer impedance for core #2 estimated using (18), we calculated the normalized induced voltage. Fig. 10 shows the comparison between the amplitude and phase of the measured and calculated normalized induced voltages when |Z t2 | was estimated using (18). As shown in Fig. 10, the differences between the measured and calculated induced voltages at lower frequencies are significantly improved by taking into account the effect of unbalances of the surface transfer impedances for cores #1 and #2 using (18), and the calculated induced voltage agrees well with the measured result over wide frequencies. The above-mentioned results validate the technique to represent tape-type shielded control cables commonly used in substations for FDTD-based electromagnetic transient simulations.

III. VOLTAGE INDUCED ON A CONTROL CABLE OVER A GROUNDING GRID OWING TO A LIGHTNING IMPULSE CURRENT
In this section, using the technique for representing shielded control cables described in Section II, we calculated the voltages induced on a control cable in a secondary circuit by a lightning impulse current injected into a grounding grid. We then compared the calculated results with the measured waveforms obtained using a test platform of primary and secondary circuits in a substation. Fig. 11 shows a test platform of primary and secondary circuits in a substation, which was mainly composed of two grounding grids, a GIS model, a control cable, and a digital-type protection relay. This test platform was similar to the experimental setup used in a previous study [41]. The two grounding grids are referred to as the main grounding grid and subgrounding grid, respectively, hereafter. The size of the main grounding grid was 30 m × 60 m, whereas the size of the subgrounding grid was 6 m × 4 m. Both grounding grids were composed of bare copper wires with a cross section of 60 mm 2 , which were buried at a depth of 0.5 m. The two grounding grids were connected via the same type of grounding wire at a depth of 0.5 m.

A. Test Platform of Primary and Secondary Circuits in a Substation
The GIS model comprised a VT, a CT, and a gas-insulated bus model. The VT and CT were the actual equipment for a 66-kV nominal voltage. Similar to a real GIS, in the gas-insulated bus model, a center conductor with a radius of 0.05 m was supported by some post-type spacers. The internal and external radii of the external sheath were 0.155 and 0.16 m, respectively. The external sheath was connected to the main grounding grid at 13 positions. The relative permittivity of SF 6 is almost one, and the presence of SF 6 does not practically influence surge propagation. Thus, to simplify the measurements, the internal space of the gas-insulated bus model was filled with air, not SF 6 .
In a control building above the subgrounding grid, the digitaltype protection relay, which was also the actual equipment for a nominal voltage of 66 kV, was installed. The secondary circuit of the CT was connected to the signal-input terminal of the protection relay, where two capacitors of 0.5 μF in series were inserted as countermeasures against transients (see Fig. 12), using an unshielded or shielded control cable. In the case of the shielded control cable, its metallic sheath was connected to the grounding grids at both ends. The type of shielded control cable was the same as the one used in Section II, whereas the unshielded control cable was also composed of two cores with a cross section of 3.5 mm 2 . One of the cores of each control cable was grounded at the protection relay, as shown in Fig. 12. The cross section of the unshielded control cable is shown in Fig. 15. The unshielded and shielded control cables were placed on the surface of the soil.
To inject a lightning impulse current into the grounding grid, a pulse generator was placed near the grounding grid and connected to the grounding grid via an insulated vinyl (IV) wire. A current injection wire (IV wire with a cross section of 2 mm 2 ) at a height of 0.5 m was connected to the other terminal of the pulse generator. The far end of the current injection wire was positioned more than 150 m away from the grounding grid, and its far end was grounded with a matched resistance.
Using the above-mentioned experimental setup, we injected a lightning impulse current into the main grounding grid and measured the following items: 1) the current injected into the grounding grid at point A (see Fig. 13), 2) the GPRs of the main grounding grid at points B and C, 3) the voltage induced on the control cable and the current flowing through the grounding wire of the shielded control cable at the protection relay (point D). To obtain the GPRs of the main grounding grid, as shown in Fig. 13, a voltage reference wire with a cross section of 2 mm 2 was placed at a height of 1.5 m. At its far end, an earth electrode was buried, and the voltage reference wire was connected through a matched resistance to the earth electrode. As represented by dotted lines in Fig. 13, the voltage reference wire was extended to point B or C for measuring the GPR of the grounding grid, that is, the voltage difference between the grounding grid and the near end of the voltage reference wire.

B. Calculation Model of the Test Platform
Fig. 14 shows an analysis space used to simulate the experimental setup in the FDTD method. The dimensions of the analysis space were 330 m × 380 m × 300 m, where the bottom space with a thickness of 150 m was treated as soil, and the relative permittivity of the soil was set to 30 [40]. Absorbing boundaries based on Liao's formulation of the second order were applied to all the external surfaces to assume an open space [66].
The two grounding grids, voltage difference wire, and current injection wire were simulated using the thin-wire representation techniques [48], [68]. A current source was placed near the main grounding grid to inject a current and connected to the grounding grid using a thin wire. The thin wire representing the current injection wire was connected to the other terminal of the current source, whereas the far end of this thin wire was attached directly to the absorbing boundary of the side surface simply to simulate a matched condition. On the other hand, the earth electrode at the far end of the voltage reference wire was modeled by a thin wire buried in the soil, because the distance between the grounding grids and the end of the voltage reference wire was important for reproducing the measured GPRs. The thin wire representing the voltage reference wire was connected through a lumped-parameter resistance corresponding to a matched resistance to the thin wire representing the earth electrode. The external sheath of the GIS model was simulated by combining some conductor plates, whereas the center conductor was simulated by a thin wire. The circular cross section of the external sheath was modeled using a staircase approximation.
The shielded control cable was represented using the technique presented in Section II. As explained in Section II-A, in this technique, voltage and current traveling on the cores inside the shielded control cable were solved by applying 1-D TL theory on the assumption of the TEM mode, whereas the electromagnetic transient phenomena outside the metallic sheath of the shielded control cable were solved by the 3-D FDTD method, where the metallic sheath and the grounding wires of the metallic sheath at both ends were straightforwardly modeled by thin wires. The thin wire used to model the metallic sheath was placed at a height of approximately 0.2 m. On the other hand, the unshielded control cable, which had a twisted-wire structure, was modeled using the technique for representing twisted-wire pairs on the basis of the hybrid technique combining the FDTD method and TL theory [69]. In this technique, voltages and currents induced on a twisted-wire pair illuminated by a nonuniform incident electromagnetic field are calculated by solving Agrawal et al.'s model [70], which is one of the field-to-transmission coupling techniques [71]. In the Agrawal et al. model, exciting electric fields are incorporated into TLtheory-based equations. The electric fields along the conductors in the twisted-wire pair in the absence of the conductors are obtained by interpolating the electric fields calculated by the 3-D FDTD method, where the twisted-wire pair is treated as an ideal bifilar helix and the position coordinate on each wire is given by mathematical expressions. In [69], the induced voltages on a twisted-wire pair and a parallel-wire pair were calculated using the above-mentioned hybrid technique, and the calculated results were validated through comparison with results obtained by the conventional TL-theory-based simulation technique. It was also confirmed that the effectiveness of twisting wires to reduce the electromagnetic coupling was well reproduced using the hybrid technique of the 3-D FDTD method and TL theory.
Similar to the shielded control cable, the p.u.l. inductance and capacitance matrices of the unshielded control cable were calculated using the conductor-subdivision-based method and the electrostatic field analysis code. Fig. 15(a) and (b) shows models used for calculating the inductance and capacitance matrices, respectively. In the calculation of the inductance matrix, the cross sections of each core (conductivity: 5.4 × 10 7 S/m) and the ground plane (conductivity: 10 8 S/m) were divided into 0.04-mm square subconductors, and the total number of the subconductors was 10 644. On the other hand, in the calculation of the capacitance matrix, the surface of each core was divided into 40 curved line elements, whereas the surface of each dielectric (relative permittivity: 4.7) was divided into 800 curved line elements. The presence of the ground plane was taken into account using image theory. The total number of variables used in the electrostatic field analysis was 13 118. Unlike in the case of the shielded control cable, to take into account the effect of the twisted-wire structure of the unshielded control cable, the above-mentioned calculation models were rotated around the center axis of the wire pair at an interval of 20°, and the inductance and capacitance matrices were obtained by averaging the calculated results.
In the case of the unshielded control cable, the ground impedance in the TL-theory-based equations was approximately given by the following expressions: where h c is the height of a conductor, and σ g , and ε rg are the conductivity and relative permittivity of soil, respectively. ε 0 and μ 0 are, respectively, the permittivity and permeability in vacuum. Here, h c was approximated to be the height of the center axis of the wire pair.
In the TL-theory-based equations, the frequency-dependent lumped admittances simulating the output impedance of the CT and the input impedance of the protection relay were connected to the ends of the distributed-parameter line model of the control cables. The frequency characteristics of the output impedance of the CT and the input impedance of the protection relay are presented in [41].
The analysis space was divided into 630 × 1001 × 307 nonuniform cells (0.05-2 m) for calculating the GPRs and 566 × 912 × 307 nonuniform cells (0.05-2 m) for calculating the voltages induced on the control cables, whereas the time discretization was set to 48.1 ps. On the other hand, the TL-theorybased equations were solved by applying the 1-D FDTD method, where the cell size was 1.38 m and the time discretization was 12 ps. In this study, we carried out the FDTD simulations using a GPU-based parallel computer with four GPUs (NVIDIA, P100), and the calculation times were 6.4 and 6.6 h for calculating the induced voltages on the unshielded and shielded control cables, respectively.

C. Comparison Between Calculated and Measured Results
By injecting a lightning impulse current into the main grounding grid, first, we calculated and measured the current injected into the main grounding grid and the GPRs of the main grounding grid using the voltage reference wire, that is, the voltage differences between the grounding grid and the voltage reference wire, in the absence of the control cables. Fig. 16 shows the calculated and measured injected currents at point A. Fig. 17(a) and (b) shows the calculated and measured GPRs at points B and C obtained using the voltage reference wire, respectively. As shown in Fig. 16, the peak values of the calculated and  measured injected currents are normalized to 1 A. As shown in Fig. 17, the calculated waveforms of the GPRs of the main grounding grid are in good agreement with the measured results, where the differences between the calculated and measured peak values are less than 2.2%. These results confirmed that the FDTD simulations well reproduced the transient responses of the grounding grids, which induce voltages on control cables when lightning currents flow into grounding grids.
Second, we connected the control cable to the CT and protection relay, and we calculated and measured the voltage induced on the control cable in the cases of the unshielded and shielded control cables. Fig. 18 shows the calculated and measured voltages induced on the unshielded control cable. In the case of the unshielded control cable, the voltage was induced on the control cable mainly by the electrostatic coupling due to the effect of the transient GPR, because one of the cores of the control cable was grounded at the protection relay, as shown in Fig. 12, and the effect of the electromagnetic coupling was suppressed by the twisted-wire structure of the control cable. Fig. 19(a) and (b) shows the calculated and measured waveforms of the voltages induced on the shielded control cable and the current flowing through the grounding wire of the metallic sheath of the control cable at the protection relay, respectively. As shown in Fig. 18, the calculated voltage induced on the unshielded control cable is in good agreement with the measured result, where the difference between the peak induced voltages is 9.8%. As shown in Fig. 19, in the case of the shielded control cable, the calculated current flowing through the metallic sheath agrees well with the measured result. Although the difference between the measured and calculated peak voltages induced on the shielded control cable is 31.8%, which is larger than that in the case of the unshielded control cable, presumably owing to the modeling error of the surface transfer impedance arising from deformation of the tape shield when it was connected in the test platform, the calculated waveform of the induced voltage agrees well with the measured result, which works for evaluating the effectiveness of the use of the shielded control cable for suppressing induced voltages.
The above-mentioned results validate the applicability of FDTD-based electromagnetic transient analysis to the calculation of voltages induced on unshielded and shielded control cables by lightning impulse currents flowing into grounding structures, which are useful for evaluating the efficacy of the installation of shielded control cables and for designing effective countermeasures against lightning threats to protect secondary circuits in substations.

IV. CONCLUSION
In this article, we applied the FDTD method, which is advantageous for analyzing the transient responses of grounding structures, to the analysis of lightning-induced transients in secondary circuits in substations. First, we studied the technique for modeling tape-type shielded control cables, which are commonly used in substations as countermeasures against transients, on the basis of the hybrid technique of the 3-D FDTD method and TL theory. We confirmed that the use of surface transfer impedances estimated by taking into account the effect of unbalances between surface transfer impedances for multicores improved the accuracy of simulations and the calculated voltages induced on a shielded control cable with multicores by sheath currents agree well with measured results over a wide range of frequencies from 1 kHz to 5 MHz, which is useful for analyzing lightning transient phenomena.
Second, we modeled a test platform consisting of a GIS, a CT, a digital protection relay, and a control cable set up over grounding grids by the 3-D FDTD method using the above-mentioned technique. We calculated the transient response of the grounding grid and voltages induced on unshielded and shielded control cables with multicores, when a lightning impulse current was injected into the grounding grid. We confirmed the applicability of FDTD-based transient analyses to the computation of lightning transient phenomena in the grounding structure and secondary circuits and to the analysis of the effectiveness of the shielded control cable for suppressing induced voltages by comparison with measured waveforms obtained using the experimental setup of the above-mentioned test platform.