Infinitesimal Method Based Calculation of Metro Stray Current in Multiple Power Supply Sections

The calculation of rail potential and stray current in the DC traction power system can guide the protection against the rail overvoltage problem and mitigate the electrolytic corrosion due to stray current. To obtain the accurate stray current distribution, the infinitesimal method based calculation of metro stray current in multiple power supply sections (PPSs) is studied in this paper. The through-type structure characteristics of power supply system and that of return current system which is constituted of rails - stray current collection mats - tunnel reinforcements - ground (RSTG) are considered in the proposed calculation method. In the traction power supply system, the currents injected into the rail by train, and absorbed by traction power substation (TPS) or regeneratively braking train, are solved by the power flow calculation which considers stray current. In the return current system, these currents are used as the boundary conditions for stray current calculation. Then the rail potential and stray current in multiple PPSs are solved by the infinitesimal method. The verification result shows that the proposed calculation method can effectively calculate the rail potential and stray current in multiple metro PPSs. At last, the novel calculation method is used to gain the dynamic stray current along the metro line and evaluate the influence factors of stray current.

INDEX TERMS Metro, stray current, rail potential, infinitesimal based calculation method, infinitesimal method. G d Conductance between rail and tunnel reinforcements G e Conductance between stray current collection mat and ground G f Equivalent conductance between rail and ground I s Equivalent current sources of train and TPS I T Train current S T Location of train I n Mesh current of traction power supply system V T Equivalent voltage of train R T Equivalent resistance of train ε Iteration accuracy U 1 Potential difference (PD) between rail and stray current collection mat U 2 PD between stray current collection mat and tunnel reinforcements U 3  In the metro system, the running rails are usually used as return paths for traction current. However, the rail can't be completely insulated from the ground, so a part of current flows into the ground, forming the stray current. Stray current accelerates the corrosion of the rails, the tunnel reinforcements, and other metal structures located around metro lines such as the underground pipeline. Therefore, in the DC metro system, a series of protection strategies for stray current is adopted in China [1]. Calculating the stray current and analyzing the distribution characteristics of the stray current are the bases of the protection strategies. To obtain the accurate stray current distribution, it is necessary to study the calculation method of stray current.

B. REVIEW AND MOTIVATION
At present, the calculation of stray current is mainly done by the simulation method [2]- [10] and the analytical method [10]- [22]. The simulation method is mainly implemented by simulation software, such as MATLAB, CDEGS, ANSYS. The effects of tunnel structure [2]- [5], crosstrack regeneration supply [6], the soil topologies [7]- [10] on stray current were studied based on the simulation software. However, when the train position changes, the model needs to be re-established in the simulation software. Therefore, the simulation method is more cumbersome for the analysis of dynamic stray current when train running. The analytical method makes up for the shortcoming of the simulation method. The analytical method is to calculate the current and potential by establishing the equivalent resistance model. According to the metal conductor topology of the DC metro system, the multi-layer resistance network model has been established [10]- [13]. Based on this, the current and potential of resistance and stray current are obtained by mathematical calculation method. Among these models, the soil topologies [10], elevated structures [11]- [12], train running mode [13] have been considered for stray current calculation and analysis. In existing resistance network models, the train current is provided and absorbed by the TPSs which are located at the head and end of the section. However, there are many TPSs in the DC metro system, as shown in Fig.1. In the Chinese DC metro system, catenaries are connected through the positive bus of the rectifier. Hence, the power supply system of the Chinese DC metro system is a through-type structure. Meanwhile, the rails along the metro line are electrically connected by fishplates or welding. The stray current collection mats for collecting stray current are laid under the track bed. The connection terminals are soldered to the stray current collection mats and electrically connected by cable to the adjacent stray current collection mats. The tunnel reinforcement in the tunnel constructed by the cut-and-cover method is electrically connected by welding. It can be used as an auxiliary stray current collection mat. TPS of the Chinese metro system usually adopts the unearthed scheme. Stray current in the Stray current collection mat and tunnel reinforcement can pass through the soil back to the negative pole of the TPS. That is, the return current system of the Chinese DC metro system consisting of rails-stray current collection mat-tunnel reinforcement-ground is also a through-type structure.
Because the power supply system and the return current system are both the through-type structure, the train current can raise up from all TPSs along the metro line when one train operates in one section. Meanwhile, the traction return current can return to all TPSs. Therefore, the calculation of stray current based on a single power supply section (PPS) is not in line with the actual DC metro system.
Besides, the currents in the power supply system are injected into the rails by trains and flow in the return current system. At the same time, the currents in the return current system are absorbed by TPSs and regenerative braking trains along the metro line. Different power flow distribution in the power supply system can cause the trend of the current distribution in the return current system to change, and vice versa. The current distribution in the power supply system and the return current system can interact with each other. That is, there is a coupling relationship between the power supply system and the return current system. In conclusion, it is necessary to study the calculation method of stray current in multiple PPSs and consider the coupling relationship between the power supply system and the return current system.
The infinitesimal method is to calculate the current and potential by establishing differential equations based on the infinitesimal model. Usually, the position of trains and TPSs are used as boundary points in the infinitesimal model. Any two adjacent boundary points can be considered as a calculation section. The current and potential in each calculation section can be expressed as the functions of rail position based on Kirchhoff's law. Then the current and potential 96582 VOLUME 8, 2020 at any position can be obtained. When the train position changes, the new functions for current and potential can be obtained based on the new boundary conditions. Therefore, the infinitesimal method could be applied to the calculation of dynamic stray current in multiple PPSs with multiple trains running. Fortunately, parameters corresponding to different calculation sections can be set to different values. Furthermore, by adding the virtual boundary points (VBPs), current sources whose value are set to 0A in a section, so that the calculation section can be divided into any number of calculation sections. As a result, stray current distribution in the return current system with the nonuniform electrical parameters can be solved. By constructing the return current system model with the multi-layer soil model, the effect of soil resistivity variability in the vertical direction on current distribution can be analyzed by the infinitesimal method too. This is in line with the fact that the return current system has different electrical performance and the geological environment is variable along the metro line.

C. TECHNICAL CONTRIBUTION AND PAPER STRUCTURE
For the DC metro system with the through-type structure, the infinitesimal method based calculation of metro stray current in multiple PPSs is proposed in this paper. The calculation method is verified by comparison with the simulation results in CDEGS [10]. Moreover, stray current and rail potential are calculated by considering the return current system with different electrical performance along the metro line, and the train operation. In comparison to the existing methods, the main contributions of the method in this paper are the following: 1) The calculation method of stray current is proposed by considering the DC metro system is the through-type structure. Base on this method, the stray current distribution in multiple PPSs can be evaluated. The method is more suitable for analyzing stray current in the actual DC metro system.
2) The coupling relationship between the power supply system and the return current system is considered in the calculation model of stray current. Therefore, the stray current distribution under this coupling relationship can be calculated in this paper.
3) The stray current calculation based on the infinitesimal method can be applied to evaluate dynamic stray current caused by trains running in multiple PPSs. Moreover, the return current system with nonuniform electrical performance and variable geological environment could be considered based on the infinitesimal method.
This paper is organized as follows. In Section II and III, the detailed calculation process of the method in this paper is explained. In Section IV, to verify the method and analyze the effect of electrical performance of the return current system on the stray current distribution, stray current and rail potentials are calculated. In Section V, dynamic stray current and dynamic rail potentials are obtained by considering that train runs along the Shijiazhuang Metro Line 3 in China. Finally, the conclusions are presented in Section VI.

II. MODELING OF STRAY CURRENT IN MULTIPLE METRO POWER SUPPLY SECTIONS
Based on the characteristic of the DC traction power system as shown in Fig. 1, the calculation model for stray current has been established. The calculation model is divided into two parts: the traction power supply system sub-model with multiple TPSs, trains; the return current system sub-model with complex stray current return path.

A. TRACTION POWER SUPPLY SYSTEM SUB-MODEL
To calculate feeder currents, the TPS can be represented by voltage source V s , and resistance r [23], [24]. Considering that the stray current is influenced by the train current and the accuracy of the train model has only a slight influence on stray current, the train can be equivalent to current source C, while the train is generally equivalent to a current source model [25] or a power model [26]. The rail is equivalent to ladder network cells consisting of rail resistance R r and leakage conductance G to model the ground leakage current, each rail cell being represented by an equivalent π -circuit [27].
If there are n TPSs represented by equivalent voltage sources of TPSs V ss (s = 1, 2, 3, . . . , n) with equivalent resistances of TPSs r s and m trains represented by equivalent current sources of trains C t (t = 1, 2, 3, . . . , m) in the metro line. Due to the different corrosion degree of rail along the metro line, the electrical parameters of RSTG may be different in different regions of the metro line. The VBP is added to the model at each division point of the adjacent regions of these regions except the positions of TPS and trains, and the number of VBPs is p. The power supply system can be divided into (m + n + p − 1) sections. For each section c(c = 1, 2, 3, . . . , m + n + p − 1), considering the overhead line resistance R ci , rail resistance R ri , and equivalent leakage conductance of rail G i in section i, the traction power supply system sub-model can be obtained as shown in Fig. 2.

B. RETURN CURRENT SYSTEM SUB-MODEL
The return current system sub-model of RSTG can be divided into several regions according to the stray current return path conditions, and the electrical parameters of the RSTG in each region can be regarded as uniform distribution. The feeder current of each TPS and the train current of each train are regarded as equivalent current sources, I sw (w = 1, 2, 3, . . . , m + n). In addition, the virtual equivalent current VOLUME 8, 2020 sources, I sv (v = 1, 2, 3, . . . , p), at these VBPs can be set to 0A. According to the geographical relationship between I sw and I sv , I sq (q = 1, 2, 3, . . . , m + n + p) is obtained by combining I sw and I sv in the sequence from the starting point to the endpoint of the line. The return current system sub-model is established as shown in Fig. 3. In some cases, the uniform soil model may lead to an underestimate of the current in the buried structures such as stray current collection mat and tunnel reinforcement, compared to multi-layer soil model. To account for the impact caused by the multi-layer soil structure, the return current system sub-model with horizontal multiple layers of soil could be established by increasing the number of soil layers, as shown in Fig. 4. This is meaningful for evaluating the corrosion of buried structures. R r , R p , R j can be calculated using the cross-sectional area, length, and resistivity of each conductor's material. Since stray current can only propagate within a certain range of soil, the equivalent soil model needs to be established to represent the area where stray current exists, as shown in Fig. 5. R sm (m = 1, . . . , u) is the resistance of the soil layer m(m = 1, . . . , u) in the x-direction. u is the number of soil layers, when u = 1, it means that the soil model is a uniform structure. According to the law of resistance, R sm can be estimated by the formula (1). Since stray current is a direct-current phenomenon, it is reasonable that the soil within a certain range is equivalent to the lumped parameter.
The value of l could be set based on reasonable assumptions and could be modified by comparing the calculation results, which are obtained by the method in this paper and the simulation in software CDEGS [10]. dx represents the unit length in the x-direction. h m is the thickness of the soil layer m. ρ m (m = 1, . . . , u) represents the resistivity of the soil layer m. The value of ρ m and h m can be determined by the inversion of parameters, which are obtained by the Wenner four-probe test method [31]. The conductance G c and G e could be estimated by the moment method [32] or the software CDEGS [10]. G a , G d , G f is mainly determined by rail insulation performance and not the soil structure. The stray current collection mat and the tunnel reinforcement are covered with concrete, so G b is also hardly affected by the soil structure. G a , G b , G d and G f could refer to [29].
The influence of the soil model on the current distribution in the return current system is complicated. The current distribution depending on the comprehensive influence of the parameters l, ρ m , h m , u of the soil model on the return current system parameters. In some cases, the smaller the value of h 1 , l and the larger the value of ρ 1 , u, the more charge will be accumulated in the buried structure of the return current system. The more accurate the parameter evaluation, the more accurate the effect of analyzing the sensitivity of the soil model on the current distribution will be in the return current system.

A. SOLUTION OF TRACTION POWER SUPPLY SYSTEM SUB-MODEL
In traction power supply system model, for the section i, the conductance between rail and stray current collection mat G ai , the conductance between rail and tunnel reinforcements G di , and the conductance between rail and ground G fi are considered to calculate G i as follow: At a specific moment, each train current I T = [I T1 , I T2 , . . . , I Tm ] and each train location S T = [S T1 , S T2 , · · · , S Tm ] can be viewed as invariant. The Mesh current of the traction power supply system sub-model, I nτ (τ = 1, 2, 3, . . . , 2m + 2n + 2p − 2), is shown in Fig. 4. The feeder current of each TPS can be solved by the mesh current method as follows: 1) The initial voltage of each train, V Tm ], is assumed as the rated voltage V e .
2) According to S T , the equivalent resistances of all meshes except the train branches are calculated, such as R c2 , R r2 , and r 2 , in Fig. 6.
4) Each I nτ is calculated by the mesh current method. The currents of each overhead line, rail and the branch current of TPS are calculated through I nτ .
5) The V (k+1) T can be recalculated by the loop constructed based on train branches, rail branches, TPS branches, and overhead line branches. Then the convergence condition is judged.
If the convergence condition is not satisfied, return to step 3) to continue the iteration. Otherwise, if the convergence condition is satisfied, the feeder currents of TPS and the train currents (i.e. I sw (w = 1, 2, 3, . . . , m + n)) can be obtained by the mesh currents.

B. SOLUTION OF RETURN CURRENT SYSTEM SUB-MODEL
The current without backtracking to the TPS through the rail is called stray current [1]. Therefore, stray current I z  is equivalent to the sum of I p , I j , and I d . The positions of the rail at I sq (q = 1, 2, 3, . . . , m + n + p) injected into are regarded as the dividing points Q q , and each section between adjacent points is regarded as a calculation section. As shown in Fig. 7, the return current system sub-model has (m + n + p − 1) calculation sections built as an infinitesimal model. When the soil structure is assumed to be uniform, R d represents the resistance of the pseudo uniform soil model R s . When the soil structure is assumed to be horizontally layered, R sm of each layer is connected in parallel with each other. R d is the value of R s1 . . . , R su in parallel. For each calculation section, if the resistance and conductance of RSTG are uniformly distributed, the infinitesimal unit of each calculation section can be represented and shown in Fig. 8. Then the following equations can be obtained by Kirchhoff's law.
Equation (6) is a first-order linear homogeneous differential equation with constant coefficients. The formula of the potential differences (PD) including U 1 , U 2 , and U 3 and of the VOLUME 8, 2020 current intensity (CI) including I r , I p , I j , and I d in RSTG at point x in each calculation section can be written in the form of (7). where Based on the solution of the first-order linear homogeneous differential equations with constant coefficients, the general solution of (7) can be written as (9).
The current direction is defined as positive from calculation section q to calculation section q + 1. According to the resistivity and conductance parameters in each calculation section, the relationships of PD and CI in return current system at point x can be obtained as the formula (8). Therefore, for the (m+n+p−1) calculation sections, there are 7(m+n+ p−1) pending coefficients K k (k = 1, 2, . . . , 7(m+n+p−1)), which can be determined by 7(m + n + p − 1) boundary conditions as following: 1) For the boundary points of nonterminal dividing points, Q q (q = 2, . . . , m + n + p − 1), the electrical quantities including I p , I j , I d and each PD calculated in the two adjacent calculation sections of boundary point Q q should be equal, and the current burst of I r at Q q equals to I sq respectively. Therefore, 7 boundary conditions can be determined at each of these boundary points, and thus 7 (m + n + p − 2) boundary conditions can be determined. The boundary conditions determined at each boundary point are: where U 1,q (Q q ), U 2,q (Q q ), and U 3,q (Q q ) are the PD at Q q in calculation section q, respectively. Similarly, I r,q (Q q ), I p,q (Q q ), I j,q (Q q ), and I d,q (Q q ) are the CI at Q q in calculation section q, respectively.
2) For the boundary point of the first TPS position and the last TPS position, CI in the RSTG at Q q (q = 1, m + n + p) all return to the TPS, and the current in the external ground remains conserved. Therefore, 7 boundary conditions can be determined: Based on the above analysis, the 7(m + n + p − 1) equations of PD and CI can be obtained, and the pending coefficients can be solved by the 7(m + n + p − 1) boundary conditions. Then, the PD and CI of the RSTG can be calculated. Prediction of the efficiency of stray current collection mats before construction of the metro system is essential to prevent stray current. Based on return current system sub-model, the efficiency η of the stray current collection mats could be assessed by [9].
where I T P is the accumulated stray current by the stray current collection mats, and I T z is the total stray current.

IV. CALCULATION AND VERIFICATION
Take the section from Xisanzhuang Station to Shijiazhuang Station of Shijiazhuang Metro Line 3 as an example. Between two adjacent TPS is called a PPS, so there are four PPSs called PPS #1, PPS #2, PPS #3, PPS #4, respectively, as shown in Fig. 9. In addition, assuming the soil is uniformly distributed. Both rails are used as the traction current return rails and the stray current direction is defined as positive from Xisanzhuang Station to Shijiazhuang Station. The rated voltage of each TPS is set to 1500V, based on the calculation of the model parameter settings [24], [28], [29]. It is assumed that r is 0.02240 and the electrical parameters of RSTG in PPS #2 are different from that of the other PPSs, as shown in Table 1.

A. METHOD VERIFICATION
Take 4 trains running on the line as an example, and the trains' operation mode is shown in Table 2.  As shown in Fig. 10, the distribution of rail potential and stray current along the metro line are calculated by the above method and verified by CDEGS simulation. CDEGS modeling techniques can be found in [5]. Meanwhile, the distributions of rail potential and stray current in a single PPS are obtained by selecting the section from Xisanzhuang station to Berlin Zhuang station only considering the operation of train 1.
According to Fig. 10, The rail potential along the line shows inflection points at the positions where train 1, train 2, and train 4 are located. Since the amplitude of the traction current of train 1 is the largest, the amplitude of the rail potential here is the largest. At the same time, due to the return current system is the through-type structure, the maximum amplitude of the stray current appears between train1 with the maximum positive traction current and train 2 with the maximum negative traction current. In addition, the efficiency of the stray current collection mats is calculated, and the value is 35.52%. TPSs adopt the unearthed scheme, currents in the stray current collection mats are only returned to the negative pole of TPSs through the ground. Therefore, the efficiency of the stray current collection mats is low.
More importantly, it is clear that the rail potential and stray current obtained by the infinitesimal based calculation method are similar to the result by the CDEGS simulation. Actually, the single PPS method neglects the influence of multiple TPSs, multiple trains and underground continuous stray current return paths of other PPSs on rail potential and stray current along the metro line.
Part of traction current flows through the electrically connected RSTG and returns to other PPSs outside the PPS where the train is located. In other words, the trains operating in different PPSs have an interactive effect on rail potential and stray current. That is accounted for in the calculation method for multiple power sections. Therefore, the rail potential and the stray current obtained by the method in multiple PPSs are obviously different from those obtained by the method in a single PPS. The result indicates that it is necessary to calculate the stray current based on the method in multiple PPSs.

B. INFLUENCE OF RETURN CURRENT SYSTEM ELECTRICAL PERFORMANCE ON STRAY CURRENT
To study the influence of the electrical performance of the rail on rail potential and stray current, the trains' operation mode is set as shown in Table 2. Using the parameters listed in Table 3 and Table 4, respectively, the rail potential and stray current along the line are calculated and shown in Fig. 11.  In Fig. 11(a), according to Ohm's law, as the rail resistance Rr increases, the amplitude of the rail potential increases because the traction current is constant. At the same time, the stray current along the metro line increases due to the amplitude of rail potential increases as shown in Fig. 11(b). VOLUME 8, 2020 The stray current is much smaller than the rail current, so the rail potential is mainly affected by the rail current and the rail resistance. In Fig. 11(c), when the rail leakage conductance changes, the rail potential does not change significantly. Based on the rail potential is hardly changed, the stray current increases significantly due to the reduced insulation between the rail to the ground. The influence of the nonuniform electrical performance of the return current system on the stray current is studied, based on the simulation parameter settings in Table 5. RP refers to the resistance parameters of the metro stray current calculation model, which includes R c , R r , R p , R j , R d . GP refers to the conductivity parameters of the metro stray current calculation model, which includes G a , G b , G c , G d , G e , G f . RP and GP are set to equal to the parameters in PPS #1 shown in Table 1. Stray current distribution and rail potential distribution in the return current system with nonuniform parameters and uniform parameters are shown in Fig. 12.
C1# indicates that the return current system parameters are uniform. C2# and C3# indicate that the return current system parameters are nonuniform. In C1#, the rail potential at the positions where trains are located is larger. At the same time, the amplitude of stray current in the middle of the metro line is relatively larger than that in the head and end of the line. In C2#, the rail potential in PPS #2 is smaller compared to C1# which the return current system with uniform parameters. The rail potential difference mainly depends on the rail resistance and rail current. The change of the GP has a small effect on the rail current. Therefore, the rail potential difference has hardly changed as shown in C2# of Fig. 12(a). In other words, the trend of the rail potential is hardly changed along the metro line. Based on the change of the amplitude of the rail potential, the stray current distribution has changed accordingly in C2#. Compared with C1#, the maximum amplitude of the stray current in C2# is larger. Since the rail potential in PPS #1 of C2# is larger, the stray current in PPS #1 is larger in C2#. Similarly, the stray current in PPS #3 of C2# and PPS #4 of C2# is smaller. In C3#, the rail potential is increased and the stray current is reduced in PPS #4 compared to C2#. The stray current in other PPSs of C3# is also reduced because the return current system is the through-type structure. Actually, metro systems often exhibit different electrical performance along the metro line. Nonuniform metro system parameters will affect the rail potential distribution and the stray current distribution along the metro line and is complicated.

V. DYNAMIC STRAY CURRENT IN MULTIPLE METRO POWER SUPPLY SECTIONS
Setting the train dwell time at 45s to obtain the typical curve of train current and train position of a single train running from Xisanzhuang station to Shijiazhuang station by train performance simulation [30], as shown in Fig. 13.
When the exit time interval is set to 300s during the steady operation period, the time interval between adjacent trains leaving Berlin Zhuang Station is 300s, the train diagram is shown in Fig. 14. In the calculation process, the calculation time window is set from 900s to 1200s. Based on this, the spatiotemporal distributions of rail potential and stray current along the section from Xisanzhuang Station to Shijiazhuang Station are solved and shown in Fig. 15. According to Fig. 15, the maximum positive rail potential is 60.9329 V, which occurs in 82 s at 1.063 km. That is, when the train is running in the traction state to 1.063km,  the maximum positive rail potential appears. At this moment, the maximum positive stray current appears, and the value is 10.37A at 3.048km. However, stray current at the position of maximum positive rail potential is 4.3561A. The maximum negative rail potential is -52.0572 V, which occurs in 205 s at 2.2 km. That is, the maximum negative rail potential appears when the train is running in the braking state to 2.2 km. Similarly, the maximum negative stray current appears at this time, which is -12.3596A at 4.78km. Stray current at the position of maximum negative rail potential is −7.1596A.
When trains are in the traction state or braking state, the rail potential is larger in the area where trains are located and the amplitude of the stray current along the line is at a high level. Moreover, the amplitude of stray current in the middle of the metro line is relatively larger than that in the head and end of the line. Since most stray current at both the head and end region of the line returns to the traction power supply system. By calculating the spatiotemporal distribution of rail potential and stray current along the line, it is valuable for the protection against the rail overvoltage problem and for the mitigation of stray current electrolytic corrosion.

VI. CONCLUSION
In this study, the calculation model of stray current in multiple PPSs is established. Based on the ''RSTG'' stray current return path, rail potential and stray current along the metro line are calculated using the infinitesimal based calculation method. The main conclusions are as follows: 1) The accuracy of the infinitesimal based calculation method is verified by comparison with the simulation results in CDEGS. The simulation result indicates that it is necessary to calculate and evaluate stray current based on the model which includes multiple PPSs.
2) The infinitesimal based calculation method can be used to solve the metro stray current with a variable number of PPSs and different return current system electrical performance in different regions. And it provides ideas for the calculation of stray current under horizontal multiple layers of soil structure.
3) The infinitesimal based calculation method has great adaptability to assess the degree of stray current during metro operation. The method can be used for the prediction and assessment of the stray current distribution along the metro line. Therefore, it is valuable for preventing overvoltage on rails and reducing stray current electrolytic corrosion. Her research interest includes analysis and mitigation of metro stray current. VOLUME 8, 2020