基于階躍級數逼近法的城市軌道交通鋼軌短路電流趨膚效應研究

喻 樂 和敬涵 王小君 薄志謙

（1.北京交通大學電氣工程學院 北京 100044 2.ALSTOM電網公司 英國斯塔福德 ST17 4LX）

1 引言

城市軌道交通牽引網大量采用第三軌（接觸軌道）供電方式。第三軌容易因絕緣老化、受損發生對走行軌或接地短路故障，嚴重威脅著軌道交通安全運行，造成不同程度的停運事故。

目前，國內城市軌道交通直流饋線保護設備主保護原理為電流上升率與電流增量保護，該方法能區分大部分故障情況，但隨著城市軌道交通的快速發展，該原理對于運行區段內多機車用電與高阻短路等情況存在難以區分的問題。因此，大量研究針對該問題展開：文獻[1-5]利用小波等方法以區分短路電流趨膚效應特征與機車正常用電電流時間常數特征；文獻[6-8]采用穩態分析方法對城市軌道交通供電系統饋線短路現象進行了分析。

當鋼軌發生短路故障時，暫態初始時刻鋼軌中電流會產生趨膚效應，使得暫態初始時刻電阻大、內電感小，但隨著時間的推移電阻逐漸減小而內電感逐漸增大，由趨膚效應引起暫態過程時間常數變化規律與MTRC地鐵短路試驗結果基本一致[9]。對鋼軌短路趨膚效應的準確計算是研究針對城市軌道交通供電系統的配套直流饋線保護、監測及故障定位等技術開發的基礎，且對直流饋線保護定值整定有著重要作用。

目前鋼軌短路電流趨膚效應的仿真計算所采用的方法主要有以下三類：①時域計算方法：由H. S.Carslaw與J. C. Jaeger[10]為解決導體熱傳導問題最先提出，并由Edward J. Tuohy等人[11]引入趨膚效應計算。文獻[12]將該算法應用于鋼軌趨膚效應計算，取得了良好的效果；②頻域計算方法：由R. J. Hill等人提出[13,14]，該方法將貝塞爾方程轉換至頻域，將電路整體進行頻域計算，最后將結果采用傅里葉反變換法轉至時域，獲得其時域解。Lawrence J.Giacoletto[15]采用Laplace變換方式進行頻域與時域的轉換，其思路與R. J. Hill較為類似，對于交流暫態計算較為實用；③導體細分法：由P.Silvester[16-18]提出，該方法的基本思路就是將整個導體劃分為有限個子部分，通過分析細分導體間的電磁關系，建立電感、電阻計算矩陣，獲得導體整體阻抗表達式從而進行計算，該方法計算結果準確，在短路計算中應用較多。

本文采用階躍級數逼近的方法對鋼軌短路趨膚效應建立了數學模型，并將該方法采用 S-Function建立鋼軌模型與牽引變電站模型在Matlab/Simulink環境下聯合仿真，提高了仿真計算速度及結果的準確性。

2 鋼軌趨膚效應計算

2.1 Matlab/Simulink仿真模型

本文基于 Matlab/Simulink建立城市軌道交通牽引供電系統模型，包含直流牽引變電站（10kV交流三相電源、二十四脈波牽引整流機組）模塊及第三軌及走行軌模塊[19]，如圖1所示。由于牽引變電站可由Matlab/Simulink自帶模塊組成，對此不再贅述。本文采用S-Function對鋼軌進行建模，將鋼軌等效為可控電流源，通過對每個計算步長獲得鋼軌電壓并計算其電流返回模型，以實現對鋼軌趨膚效應特性仿真計算。

圖1 牽引供電系統Matlab/Simulink模型圖Fig.1 Matlab/Simulink model of DC railway power supply system

由于鋼軌具有較大的不規則橫截面，其截面示意圖如圖2所示。趨膚效應使得鋼軌電阻和內電感的計算變得非常復雜。本文采用等效截面積法將鋼軌等效為等截面積的圓柱形導體，截面半徑為

式中，A為鋼軌截面積；R為鋼軌等效半徑。

圖2 鋼軌等效圓柱體截面示意圖Fig.2 Equivalent circular cylinder conductor of rail

分析和計算任何電磁場問題的基本出發點都是麥克斯韋的電磁場方程組。將麥克斯韋電磁場基本方程轉換為相量形式

式中，μ 為磁導率；σ 為電導率；J 為導體電流密度。

將式（2）在如圖3所示單位圓柱體模型下展開可得

圖3 單位圓柱形導體模型Fig.3 Model of unit circular cylinder conductor

式（3）為圓柱型導體在時變電磁場中考慮趨膚效應時的求解方程，是典型的貝塞爾方程。對于此類方程的求解過程需要對實際應用中的邊界條件展開討論并采用貝塞爾函數求解獲得。

2.2 導體表面疊加階躍電壓信號的趨膚效應計算

考慮以下情況：階躍電壓信號疊加在圓柱形導體表面：

（1）在 t=0時刻：根據階躍信號定義可知：dU/dt→∞，即信號角頻率趨于無窮大。根據趨膚效應定義

式中，Js為導體表面電流密度；ρ為導體電阻率；ω為電流角頻率，ω=2πf；d為電流與導體表面間的距離。

此時導體除表面電流密度以外，內部電流密度為零。因此，此時導體表面電流密度可采用Js=σEs=σU；計算獲得由于導體內部電流趨于無窮小，因此導體電流為無窮小。在實際計算中，可以近似認為It=0=0。

在 t→∞時刻：導體內部任意位置電流密度等于Js。因此，此時電流大小為：It=∞=πR2Js。

以上兩種情況可以視為電流計算的邊界條件。

計算導體表面疊加單位階躍電場強度E時導體電流密度J(r，t)關于時間t與半徑r的解為

式中，J0為零階第一類貝塞爾函數；J1為一階第一類貝塞爾函數；αn為零階第一類貝塞爾函數的根。

圓導體內電流可由式（5）積分獲得

將式（5）代入式（6）可得圓柱形導體中的時域電流解為

在t =0與t =∞間的時間電流大小可以通過暫態內電阻與內電感計算獲得。

一種暫態內電阻定義方法為通過功率損失定義[11]

另一種適用于時間較短時的電阻定義為

導體內電感可以采用內部磁鏈計算

式中，ψint為導體內部磁鏈；Ik為半徑為r的圓形積分路徑內的電流；IR為圓導體內電流。

將式（7）~式（10）代入式（11）計算獲得導體電壓，并令 Js=Uσ =1·σ，即穩態時電壓為 1V時所產生導體表面電流密度；為使計算結果準確，本文在0～0.000 1s內，采用Rint2作為內電感計算，部分結果見表1；在0.000 1～1s內，采用Rint1作為內電感計算，部分結果見表2。

表1 時間較短時電壓計算結果Tab.1 Voltage calculation results with short time

表2 時間較長時電壓計算結果Tab.2 Voltage calculation results with long time

由表1與表2可知：計算結果除在初始時刻與1V有一定誤差，其他時刻基本符合。該誤差可能由內電感定義造成。并由式（11）可知：導體疊加電壓信號幅值與電流響應幅值為正比關系。因此，導體疊加階躍電壓信號所產生的電流響應可以采用式（7）進行計算。

2.3 階躍電壓級數分析方法

根據上節分析，可知對于鋼軌疊加單位階躍電壓信號所產生電流響應可用式（7）表示。因此，對電壓信號的電流響應，本文采用信號的階躍電壓級數逼近方法，其步驟如下：

（1）將對鋼軌施加的電壓信號采用階躍函數的級數逼近方法進行表示。設在離散采樣的形式下，電壓信號 f(t)可以用下式離散表示：f[n]=f(nT)，T為采樣時間間隔。

將電壓離散信號 f[n]用單位階躍信號 U(t)級數進行逼近

式中，當n=0時，f[(n-1)T]=0。

由f[n]的級數表達式可知n-1，f[n]可表示為一系列階躍電壓分量的疊加形式，而單位階躍電壓U(nT)的系數為前后兩次采樣電壓值的差，即{f(nT)-f[(n-1)T]}。

（2）對鋼軌施加單位階躍電壓信號U(t)時所產生的電流響應為 I(t)，即式（7），所以單位階躍電壓U(nT)所對應的電流響應為I(t-nT)。因此，階躍電壓分量{f(nT)-f((n-1)T)}U(nT)所對應的電流響應可表示為{f(nT)-f((n-1)T)}I(t-nΤ)。

（3）計算電壓信號用階躍級數表示的逼近式中的每個階躍電壓分量所對應的電流響應，并將這些電流響應進行時域內的組合，即電流響應的時域解：

該方法具體計算過程如圖4所示。

圖4 階躍電壓級數分解Fig.4 Analysis of step series of voltage method

2.4 頻域電流計算驗證

本文采用穩態頻域計算方法對該算法進行驗證。驗證方法為假設圓導體上疊加高頻電壓，計算其頻域阻抗，從而獲得其電流響應。圓柱形導體頻域阻抗[1,9,20,21]為

其中

因此導體電流可以通過下式計算獲得

采用不同頻率單位交流電壓對階躍級數算法進行計算。測試中，采樣率為 10kHz，計算結果如圖5～圖8所示，與頻域電流計算結果對比見表3，可以發現：在高頻情況下，該算法與頻域計算方法基本符合，這證明了該算法的有效性。

圖5 100Hz電流響應Fig.5 Current response of 100Hz excitement

圖6 300Hz電流響應Fig.6 Current response of 300Hz excitement

圖7 700Hz電流響應Fig.7 Current response of 700Hz excitement

圖8 1kHz電流響應Fig.8 Current response of 1kHz excitement

表3 電壓階躍級數計算與頻域阻抗計算結果比較Tab.3 Comparison between calculation results of step series voltage method and frequency domain impedance method

3 地鐵供電系統仿真及分析

本文采用文獻[9]所提供實驗參數對 Matlab/Simulink模型進行測試。該實驗對倫敦地鐵進行了短路測試，短路點距離牽引變電站 2.882km，采樣率為10kHz。計算結果與測試結果對比見表4。

表4 電壓階躍級數計算與實測結果比較Tab.4 Comparison between calculation results of step series voltage method and test results

從以上結果可知，該模型很好地反應了鋼軌在短路時趨膚效應的暫態特性，計算初始時刻存在一定誤差，為仿真計算與實際實驗中開關暫態特性有所不同導致。

4 結論

本文首先對鋼軌等效圓柱形導體疊加階躍電壓信號這一情況展開分析，采用信號的階躍級數分析方法建立鋼軌趨膚效應模型，并用頻域計算對其進行了驗證。最后，將該方法采用S-Function模塊在Matlab/Simulink仿真實現，并與城市軌道交通短路試驗實測數據進行對比，取得了良好的效果。

[1]衛志農, 鄭玉平, 何樺. 基于小波變換的地鐵遠方短路電流分析[J]. 電力系統及其自動化學報. 2003,15(6): 5-8.Wei Zhinong, Zheng Yuping, He Hua. The analyzing of distant short-circuit fault current in rapid transit system s based on wavelet transform[J]. Proceedings of the CSU-EPSA, 2003, 15(6): 5-8.

[2]Chang C S, Kumar S, Liu B, Khambadkone A,Real-time detection using wavelet transform and neural network of short-circuit faults within a train in DC transit systems[J]. IEE Procedings-Elects. Power Applications, 2001, 148(3): 251-256.

[3]Yu L, He J H, Yip H T, et al. A new DC feeder protection based on wavelet transform[C]. Powertech,2009 IEEE Bucharest, 2009.

[4]張秀峰, 茍紹波. 利用小波變換檢測地鐵饋線故障[J]. 西南交通大學學報, 2002, 27(5): 548-552.Zhang Xiufeng, Gou Shaobo. Fault detection of subway feed line by means of wavelet transform[J].Journal of Southwest Jiaotong University, 2002, 27(5):548-552.

[5]孔瑋, 崔運海, 秦立軍,等. 直流牽引系統饋線微機保護裝置[J]. 電力系統自動化, 2005, 29(1): 92-96.Kong Wei, Cui Yunhai, Qin Lijun, et al. Research and development of microcomputer-based feeder protection for dc traction power supply system[J].Automation of Electric Power Systems, 2005, 29(1):92-96.

[6]張勛, 李夏青. 直流牽引供電系統的建模以及參數分析[J]. 北京石油化工學院學報. 2009, 19 (1): 5-9.Zhang Xun, Li Xiaqing. Modeling and parameter analysis of DC traction power supply system[J].Journal of Beijing Istitute of Petro-Chemical Technology, 2009,19 (1): 5-9.

[7]張海波, 張穎. 直流牽引供電系統穩態短路分析[J].城市軌道交通研究, 2010(3): 60-65, 68.Zhang Haibo, Zhang Ying. Analysis of steady-state short-circuit in DC traction power supply system[J].Urban Mass Transit, 2010(3): 60-65, 68.

[8]黃維軍, 周密, 王瀅. 直流牽引網穩態短路電流計算[J]. 電氣化鐵道, 2010(1):39-41.Huang Weijun, Zhou Mi, Wang Ying. DC traction net steady state short circuit current calculation[J].Electric Railway, 2010(1): 39-41.

[9]Brown J C, Allan J, Mellitt B. Calculation of remote short circuit fault currents for DC railways[J]. IEE Proceedings-B. I992, 139(4): 289-294.

[10]Carslaw and Jaeger. Conduction of heat in solids[M].1st ed. London: Oxford University Press, 1947.

[11]Edward J T, Thomas T H, Fullerton H P, Transient resistance of conductors[J]. IEEE Transactions on Power Apparatus and Systems, 1968, 87 (2): 455-462.

[12]孔瑋. 城市軌道交通直流牽引系統故障分析及若干問題的研究[D]. 北京: 華北電力大學, 2005.

[13]Hill R J, Carpenter D C. Modelling of nonlinear rail impedance in ac traction power systems[J]. IEEE Transactions on Power Delivery, 1991, 6(4): 1755-1761.

[14]Hill R J, Carpenter D C, Determination of rail internal impedance for electric railway traction system simulation[J]. IEE Proceedings-B 1991, 138 (6): 311-321.

[15]Giacoletto L J, Frequency-and-time-domain analysis of skin effects[J]. IEEE Transactions on Magnetics,1996, 32(1): 220-229.

[16]Silvester P, modal theory of skin effect in single and multiple turn coils[J]. IEEE Transactions on Power Apparatus and Systems, 1970, 89(1): 29-34.

[17]Silvester P. AC resistance and reactance of isolated rectangular conductors[J]. IEEE Transactions on Power Apparatus and Systems, 1967, 86 (6):770-774.

[18]Silvester P, The accurate calculations of skin effect in conductors of complicated shape[J]. IEEE Transactions on Power Apparatus and Systems, 1968, 87 (3):735-741.

[19]Yu L, He J H, Hu J, et al. Accurate track modeling for fault current on DC railways based on MATLAB/Simulink[C]. IEEE Power and Energy Society General Meeting, 2010.

[20]諸葛向彬. 高頻載流圓柱形導體的電阻計算[J].杭州大學學報, 1995(5): 74-77.Zhuge Xiangbing. The calculation of resistance of high frequency current cylinder conductor[J]. Journal of Hangzhou University, 1995(5): 74-77.

[21]崔運海. 直流牽引饋線微機保護裝置的研制[D]. 北京：華北電力大學, 2005.