一種電磁場(chǎng)-傳輸線組合的時(shí)域有限差分方法

胡 榕，崔 翔，陳維江，張衛(wèi)東

(1.華北電力大學(xué)新能源電力系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京102206;2.國(guó)家電網(wǎng)公司，北京100031)

0 引言

電力系統(tǒng)中的氣體絕緣開(kāi)關(guān)設(shè)備(Gas Insulated Switchgear，GIS)在隔離開(kāi)關(guān)操作過(guò)程中產(chǎn)生的快速瞬態(tài)過(guò)電壓(Very Fast Transient Overvoltage，VFTO)會(huì)影響GIS 及相連設(shè)備的安全運(yùn)行［1-2］。在設(shè)計(jì)建造GIS 設(shè)備之初，需要通過(guò)建模仿真計(jì)算獲得GIS 設(shè)備的VFTO 水平及其分布規(guī)律，保障設(shè)備的安全運(yùn)行。

通常可以采用等效電路模型［3-5］或三維電磁場(chǎng)模型［6-7］對(duì)GIS 設(shè)備進(jìn)行建模和時(shí)域仿真計(jì)算。

等效電路模型因模型簡(jiǎn)單、參數(shù)容易提取等特點(diǎn)被廣泛使用，特別是國(guó)際大電網(wǎng)組織(Conseil International des Grands Réseaux électriques，CIGRE)工作組WG33/13.09 推薦了GIS 設(shè)備各元件的等效電路模型后，眾多學(xué)者在此基礎(chǔ)上開(kāi)展了大量的建模仿真研究工作［3-5］。在這些研究中，GIS 設(shè)備的一些簡(jiǎn)單傳輸線結(jié)構(gòu)如GIS 管母線、分支管母線等元件，若采用等效的傳輸線電路模型則可以獲得較高的計(jì)算精度。但對(duì)于GIS 設(shè)備的一些復(fù)雜結(jié)構(gòu)如隔離開(kāi)關(guān)、套管、端部等元件，受等效電路模型及其參數(shù)的限制，等效電路模型無(wú)法反映高頻瞬態(tài)電磁波在其內(nèi)部的傳播過(guò)程，使得等效電路模型適用于低頻情況，與實(shí)測(cè)的VFTO 波形相比，采用等效電路模型進(jìn)行仿真計(jì)算得到的波形只能夠大致反映VFTO的波形走向，在波形細(xì)節(jié)上不能很好反映實(shí)際的VFTO 特性。

采用三維電磁場(chǎng)模型計(jì)算GIS 設(shè)備中的VFTO，能夠很好地模擬GIS 設(shè)備中復(fù)雜結(jié)構(gòu)及其對(duì)VFTO 的影響［6-7］，且具有較高的計(jì)算精度，但存在模型復(fù)雜、剖分困難、消耗大量計(jì)算機(jī)內(nèi)存及計(jì)算時(shí)間，特別是對(duì)于多分支、大尺度的GIS 設(shè)備的建模，普通的計(jì)算機(jī)是難以勝任的，常需要高性能的服務(wù)器進(jìn)行模型建立和計(jì)算。

縱觀GIS 設(shè)備，具有傳輸線結(jié)構(gòu)的元件占設(shè)備的絕大部分，隔離開(kāi)關(guān)、套管和端部等具有三維非傳輸線結(jié)構(gòu)的元件僅占設(shè)備的少數(shù)部分。對(duì)于具有傳輸線結(jié)構(gòu)的元件，采用傳輸線等效電路模型具有較高的計(jì)算精度，沒(méi)有必要采用三維電磁場(chǎng)模型;對(duì)于三維非傳輸線結(jié)構(gòu)的元件，其內(nèi)部波過(guò)程復(fù)雜，傳輸線等效電路模型會(huì)帶來(lái)較大的計(jì)算誤差，只有采用三維電磁場(chǎng)模型，才能獲得具有較高計(jì)算精度的仿真結(jié)果。均衡計(jì)算精度和計(jì)算代價(jià)，并充分發(fā)揮等效電路模型和三維電磁場(chǎng)模型各自的優(yōu)點(diǎn)，本文針對(duì)GIS 設(shè)備的VFTO 的建模仿真需求，提出了一種場(chǎng)-路組合的建模方法，即將三維電磁場(chǎng)時(shí)域有限差分方法與一維傳輸線時(shí)域有限差分方法相組合的建模方法(finite difference time domain combined with 3D and 1D，3D-1D FDTD)。

本文首先介紹3D-1D FDTD 方法中一維時(shí)域有限差分方法(1D FDTD)和三維時(shí)域有限差分方法(3D FDTD)的計(jì)算格式。然后重點(diǎn)給出1D FDTD 和3D FDTD 在連接截面上的耦合邊界條件以及迭代計(jì)算格式，給出3D-1D FDTD 方法時(shí)域計(jì)算的穩(wěn)定條件。最后，以一段GIS 設(shè)備長(zhǎng)直管母線為例，分別采用本文提出的3D-1D FDTD 計(jì)算方法和純一維傳輸線等效電路模型計(jì)算方法計(jì)算長(zhǎng)直管母線不同位置的瞬態(tài)電壓和瞬態(tài)電流，對(duì)比該兩種計(jì)算方法的計(jì)算結(jié)果，驗(yàn)證3D-1D FDTD 方法的正確性和有效性。

1 時(shí)域有限差分方法

本文提出的3D-1D FDTD 建模方法，是建立在1D FDTD 方法和3D FDTD 方法基礎(chǔ)上的。為后續(xù)討論方便，下面，簡(jiǎn)要介紹這兩種方法，并給出相應(yīng)的計(jì)算公式。

1.1 1D FDTD 方法

1D FDTD 方法是將傳輸線方程進(jìn)行時(shí)域有限差分形式離散。假設(shè)入射電磁波在傳輸線上產(chǎn)生的分布電壓源和分布電流源分別為UF和IF，則頻域的傳輸線方程為［8-10］:

式中:Z 為單位長(zhǎng)度的串聯(lián)阻抗矩陣;Y 為單位長(zhǎng)度的并聯(lián)導(dǎo)納矩陣。若Z 和Y 中包含頻變參數(shù)，則在復(fù)頻域中采用矢量匹配法將矩陣Z 和Y 按極點(diǎn)進(jìn)行部分分式展開(kāi)［8-9］，得到

式中:L 和C 分別為電感和電容矩陣，該矩陣為常數(shù)矩陣;bk和qk為留數(shù)矩陣;ak和pk為極點(diǎn);Ns為展開(kāi)的項(xiàng)數(shù)。將式(2)代入式(1)中，初始條件為0 時(shí)，進(jìn)行單邊拉普拉斯變換，得到

對(duì)式(3)中的卷積項(xiàng)采用一階中心差分形式進(jìn)行離散展開(kāi)，式(3)中第一個(gè)等式令l = (m +1/2)Δl，t = (n +1)Δt;第二個(gè)等式令l = mΔl，t =(n +1/2)Δt;Δl、Δt 分別為空間步長(zhǎng)和時(shí)間步長(zhǎng);m、n 分別為空間第m 個(gè)節(jié)點(diǎn)和時(shí)間第n 個(gè)節(jié)點(diǎn);若I(l，(n +1/2)Δt)記為In+1/2l，則式(3)中的卷積項(xiàng)整理得到

將式(4)代入式(3)，應(yīng)用一階中心差分公式對(duì)式(3)進(jìn)行離散展開(kāi)，整理獲得

式(5)為式(1)傳輸線方程的時(shí)域有限差分迭代計(jì)算格式，考慮了阻抗和導(dǎo)納矩陣中參數(shù)的頻變，通過(guò)起始條件及入射的分布電壓源與電流源激勵(lì)，經(jīng)過(guò)時(shí)間和空間上的迭代，可依次計(jì)算出傳輸線沿線各個(gè)時(shí)刻的電壓和電流。

1.2 3D FDTD 方法

3D FDTD 方法是將Maxwell 方程組中的旋度方程進(jìn)行時(shí)域有限差分形式離散。若考慮有損耗介質(zhì)，且介質(zhì)為各向同性線性介質(zhì)，則Maxwell 旋度方程和本構(gòu)關(guān)系為

式中:Je、Jm分別是電流密度和磁流密度;ε 為介質(zhì)的介電常數(shù)，(F/m);μ 為磁導(dǎo)率，(H/m);δ 為電導(dǎo)率，(S/m);δm是磁電阻率，(Ω/m)，對(duì)應(yīng)磁場(chǎng)的損耗;對(duì)于無(wú)損耗的介質(zhì):δ = 0，δm= 0。

設(shè)計(jì)算的區(qū)域是無(wú)源空間，對(duì)計(jì)算區(qū)域以Yee晶胞作為剖分單元進(jìn)行剖分，將式(6)進(jìn)行時(shí)間和空間上的一階中心差分形式離散，整理得到的6個(gè)分量的時(shí)域有限差分格式為［11-12］

式中:Δx、Δy 和Δz 分別為x、y、z 方向上的空間步長(zhǎng)。C、Cˊ、D 和Dˊ的表達(dá)式如式(8)，式中η 取x、y、z，Δt 為時(shí)間步長(zhǎng)。

對(duì)于計(jì)算區(qū)域的吸收邊界采用單軸各向異性介質(zhì)完全匹配層(uniaxial perfectly matched layer，UPML)。UPML 是Sacks (1995年)和Gedney(1996年)提出的，3D FDTD 的UPML(含棱邊和角頂區(qū))時(shí)域有限差分格式的計(jì)算公式見(jiàn)文獻(xiàn)［11-12］，本文的UPML 的最外側(cè)截?cái)噙吔邕x擇為理想導(dǎo)體邊界。結(jié)合計(jì)算區(qū)域的初始條件和UPML 吸收邊界條件，通過(guò)式(7)和式(8)在時(shí)間和空間上依次迭代，可以計(jì)算出區(qū)域中任何一點(diǎn)的電場(chǎng)和磁場(chǎng)。

2 3D-1D FDTD 方法

2.1 耦合邊界

1D FDTD 方法和3D FDTD 方法在連接截面上的耦合邊界需要滿足的條件為:邊界處必須滿足橫向電磁波(transverse electromagnetic wave，TEM)或準(zhǔn)TEM 波傳播的條件即電磁波傳播的方向上沒(méi)有電場(chǎng)和磁場(chǎng)分量或者電磁波傳播的方向上有電場(chǎng)和磁場(chǎng)的分量，但它們必須遠(yuǎn)遠(yuǎn)小于橫向分量。以圓導(dǎo)體結(jié)構(gòu)的傳輸線為例說(shuō)明耦合邊界處的處理情況，傳輸線填充介質(zhì)中的剖分單元為Yee 晶胞。

圖1(a)中展示了耦合邊界處填充介質(zhì)中不同時(shí)刻的電流(I)、電壓(U)、電場(chǎng)強(qiáng)度(E)、磁場(chǎng)強(qiáng)度(H)之間的處理關(guān)系。可知，在ti+1/2時(shí)刻，依據(jù)式(7)，邊界處的H 由前一時(shí)刻的值和相鄰的電場(chǎng)強(qiáng)度計(jì)算得到。而I 則由H 的環(huán)線積分計(jì)算得到，計(jì)算公式為

圖1(b)中展示出環(huán)線積分的路徑c，要求積分路徑c 為最靠近中心金屬導(dǎo)體的環(huán)線。依照式(5)，在ti+1時(shí)刻，電壓U 是由前一時(shí)刻的值及相鄰的電流計(jì)算得出。電場(chǎng)強(qiáng)度E 則由U 給出，在耦合邊界處忽略傳播方向上的電場(chǎng)強(qiáng)度分量;假設(shè)中心金屬導(dǎo)體和金屬外殼之間的電壓為標(biāo)幺值1，由靜態(tài)電磁場(chǎng)理論計(jì)算獲得填充介質(zhì)中各個(gè)剖分單元的節(jié)點(diǎn)電位，電位分布如圖1(c)所示，剖分單元上的電場(chǎng)強(qiáng)度標(biāo)幺值等于相鄰的節(jié)點(diǎn)電位除以節(jié)點(diǎn)之間的距離，電場(chǎng)強(qiáng)度的方向?yàn)閅ee 單元規(guī)定的方向，因而邊界的電場(chǎng)強(qiáng)度E 等于電壓U 乘以電場(chǎng)強(qiáng)度標(biāo)幺值。采用上述給出的耦合邊界條件和處理關(guān)系，建立邊界處I、U、E、H 之間的計(jì)算關(guān)系，實(shí)現(xiàn)1D FDTD 方法與3D FDTD 方法之間的銜接。

圖1 3D-1D FDTD 的耦合邊界Fig.1 Coupled boundary of 3D-1D FDTD

2.2 穩(wěn)定條件

3D-1D FDTD 方法的穩(wěn)定運(yùn)行取決于時(shí)域有限差分格式的穩(wěn)定運(yùn)行。穩(wěn)定運(yùn)行主要集中表現(xiàn)在方法的時(shí)間步長(zhǎng)、信號(hào)的最小波長(zhǎng)與空間步長(zhǎng)之間的關(guān)系上。式(5)中1D FDTD 的空間步長(zhǎng)Δl，式(7)中3D FDTD 的空間步長(zhǎng)Δx、Δy、Δz 滿足Courant 穩(wěn)定條件時(shí)［11-12］:

式中:υ 為信號(hào)的最大相速度;Δt 為計(jì)算時(shí)間步長(zhǎng)。若取Δx = Δy = Δz = Δl = Δs，則3D-1D FDTD 方法的穩(wěn)定運(yùn)行條件為

式中:λmin為信號(hào)的最小波長(zhǎng);離散系數(shù)β 一般選取20～100。依據(jù)式(11)和式(12)合理選取α 和β 使得3D-1D FDTD 方法在合理的數(shù)值色散范圍內(nèi)穩(wěn)定運(yùn)行。

2.3 計(jì)算流程

圖2 展示了3D-1D FDTD 方法的計(jì)算流程。首先建立計(jì)算區(qū)域的3D-1D FDTD 計(jì)算模型，給定激勵(lì)源和計(jì)算時(shí)間，依據(jù)式(11)和式(12)設(shè)置空間步長(zhǎng)和時(shí)間步長(zhǎng);然后依據(jù)式(5)、式(7)及耦合邊界上的處理關(guān)系在空間和時(shí)間上進(jìn)行迭代計(jì)算區(qū)域內(nèi)的電壓、電流、電場(chǎng)強(qiáng)度和磁場(chǎng)強(qiáng)度;最后輸出計(jì)算結(jié)果，完成區(qū)域內(nèi)計(jì)算。

圖2 3D-1D FDTD 方法的計(jì)算流程圖Fig.2 Calculation diagram of 3D-1D FDTD method

3 計(jì)算示例

以一段GIS 設(shè)備長(zhǎng)直管母線為例，分別采用純一維傳輸線等效電路模型和3D-1D FDTD 方法進(jìn)行建模計(jì)算。長(zhǎng)直管母線長(zhǎng)度為1.5 m，采用3D-1D FDTD 方法時(shí)，如圖3所示，將管母線分成3 個(gè)長(zhǎng)度為0.5 m 計(jì)算區(qū)域，分別為一維傳輸線等效電路模型段-三維電磁場(chǎng)模型段-一維傳輸線等效電路模型段，管母線兩端分別施加電壓源激勵(lì)和負(fù)載。計(jì)算區(qū)域的空間步長(zhǎng)為2.5 ×10-3m，時(shí)間步長(zhǎng)為4.8 ×10-12s，設(shè)置3 個(gè)觀測(cè)點(diǎn)，都位于計(jì)算區(qū)域的中心位置。另外采用純一維傳輸線等效電路模型建立的長(zhǎng)直管母線的仿真模型，時(shí)間步長(zhǎng)取值與3D-1D FDTD 方法一致。

圖3 GIS 管母線的3D-1D FDTD 計(jì)算模型Fig.3 3D-1D FDTD computational model of GIS bus

采用3D-1D FDTD 方法計(jì)算時(shí)，觀測(cè)點(diǎn)2 的輸出為電場(chǎng)強(qiáng)度E 和磁場(chǎng)強(qiáng)度H，為了比對(duì)方便，該點(diǎn)處電壓是由E 沿管母線徑向積分的平均值計(jì)算得到，而電流則由式(9)計(jì)算得到。下面，分別對(duì)方波、三角波和正弦三種不同形式的電壓源激勵(lì)產(chǎn)生的瞬態(tài)電壓和瞬態(tài)電流進(jìn)行時(shí)域計(jì)算。

當(dāng)電壓激勵(lì)源是上升時(shí)間為11 ns、幅值為1 V 的方波時(shí)，計(jì)算時(shí)間為20 000 個(gè)時(shí)間步長(zhǎng)。3D-1D FDTD 方法和純一維傳輸線等效電路模型的計(jì)算結(jié)果在3 個(gè)觀測(cè)點(diǎn)處的電壓、電流波形及頻譜圖如圖4所示。結(jié)果顯示，兩種方法計(jì)算得到的各觀測(cè)點(diǎn)處電壓、電流的波形走向一致。就同一觀測(cè)點(diǎn)而言，兩者的電壓波形吻合;采用純一維傳輸線等效電路模型計(jì)算得到的電流略高于3D-1D FDTD 方法計(jì)算的電流值。兩種方法計(jì)算得到的同一觀測(cè)點(diǎn)處電壓、電流波形的頻譜圖的頻率成分和幅值一致。

圖4 上升時(shí)間為11 ns 方波電壓源的計(jì)算結(jié)果Fig.4 Calculation results of square voltage source as 11 ns rise time

圖5 展示了電壓激勵(lì)源的上升時(shí)間為11 ns、下降時(shí)間為11 ns、幅值為1 V 的三角電壓波時(shí)，計(jì)算時(shí)間為20 000步，3 個(gè)觀測(cè)點(diǎn)處的電壓和電流的計(jì)算波形?？芍?，在同一觀測(cè)點(diǎn)處，兩種方法計(jì)算獲得的電壓波形吻合;電流波形初期，純一維傳輸線等效電路模型計(jì)算獲得的電流略大于3D-1D FDTD 方法計(jì)算的電流值;電流波形后期，兩種方法計(jì)算獲得的電流波形吻合。

電壓激勵(lì)源是頻率為100 MHz、幅值為1 V 的正弦電壓波時(shí)，計(jì)算時(shí)間為5 000步。3 個(gè)觀測(cè)點(diǎn)處計(jì)算得到的電壓和電流波形見(jiàn)圖6?？梢?jiàn)，在同一觀測(cè)點(diǎn)處，3D-1D FDTD 方法和純一維傳輸線等效電路模型分別計(jì)算得到的電壓波形和電流波形一致。

圖5 上升/下降時(shí)間為11 ns 三角波電壓源的計(jì)算結(jié)果Fig.5 Calculation results of triangular voltage source as 11 ns rise/fall time

圖6 頻率為100 MHz 正弦波電壓源的計(jì)算結(jié)果Fig.6 Calculation results of sinusoidal voltage source as 100 MHz frequency

從圖4、5、6 中可知，波形無(wú)發(fā)散，在每個(gè)觀測(cè)點(diǎn)處，3D-1D FDTD 方法計(jì)算獲得的電壓和電流波形與純一維傳輸線等效電路模型仿真獲得的波形吻合。兩者波形的微小差異是由于3D-1D FDTD 方法是基于時(shí)域有限差分法，方法本身具有剖分誤差和數(shù)值色散，這些不可避免的誤差使得兩種方法計(jì)算獲得的波形有細(xì)微差異。與純一維傳輸線等效電路模型的計(jì)算結(jié)果比對(duì)后，說(shuō)明了3D-1D FDTD 方法是正確的、有效的。

4 結(jié)論

本文提出了一種場(chǎng)-路組合的建模方法，即將三維電磁場(chǎng)時(shí)域有限差分方法與一維傳輸線時(shí)域有限差分方法相組合的建模方法(3D-1D FDTD)。給出了3D-1D FDTD 方法在連接截面上的耦合邊界條件以及迭代計(jì)算格式，提出了時(shí)域計(jì)算的穩(wěn)定條件。以一段GIS 設(shè)備長(zhǎng)直管母線為例，計(jì)算了不同位置的瞬態(tài)電壓和瞬態(tài)電流，將計(jì)算結(jié)果與純一維傳輸線等效電路模型的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文方法的有效性。

［1］陳維江，李志兵，孫崗，等.特高壓氣體絕緣開(kāi)關(guān)設(shè)備中特快速瞬態(tài)過(guò)電壓特性的試驗(yàn)研究［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2011，31 (31):38-47.

［2］胡榕，崔翔，陳維江，等.特高壓氣體絕緣開(kāi)關(guān)設(shè)備瞬態(tài)外殼電壓特性的試驗(yàn)研究［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2014，34 (29):5244-5258.

［3］Grandl J，Eriksson A，Meppelink J，et al.Studies of very fast transient (VFT)in a 765 kV substation［R］.Paris:CIGRE，1988.

［4］Povh D，Schmitt H，Valcker O，et al.Modelling and analysis guidelines for very fast transients［J］.IEEE Transactions on Power Delivery，1996，11 (4):2028-2035.

［5］Okabe S，Takami J，Nojima K，et al.Circuit models of gas insulated switchgear elements for electromagnetic wave leakage analysis［J］.IEEE Transactions on Dielectrics and Electrical Insulation，2008，15 (5):1439-1448.

［6］Smajic J，Holaus W，Kostovic J，et al.3D Full-maxwell simulations of very fast transients in GIS［J］.IEEE Transactions on Magnetics，2011，47 (5):1514-1517.

［7］Hoshino T，Maruyama S，Sakakibara T.Simulation of propagating electromagnetic wave due to partial discharge in GIS Using FDTD［J］.IEEE Transactions on Power Delivery，2009，24 (1):153-159.

［8］Paul C R.Analysis of multiconductor transmission lines［M］.Hoboken，NJ:John Wiley，1994:252-358.

［9］齊磊.變電站瞬態(tài)電磁場(chǎng)對(duì)二次電纜的電磁耦合機(jī)理研究［D］.保定:華北電力大學(xué)，2006.

［10］高方平，姚纓英，季蘇蕾.用FDTD 法求解傳輸線方程［J］.華北電力大學(xué)學(xué)報(bào)，2012，39 (2):12-16.

［11］王長(zhǎng)清，祝西里.瞬變電磁場(chǎng)-理論和計(jì)算［M］.北京:北京大學(xué)出版社，2011.

［12］葛德彪，閆玉波.電磁波時(shí)域有限差分方法［M］.西安:西安電子科技大學(xué)出版社，2011.