張 萍 劉 寧 聶鑫鵬 吉增強

(1.省部共建電工裝備可靠性與智能化國家重點實驗室(河北工業(yè)大學) 天津 300130；2.河北工業(yè)大學人工智能與數據科學學院 天津 300130；3.中電建(?？?新能源有限公司 ?？?570100；4.河北燃氣有限公司 石家莊 050000)

1 引言

單/多導體傳輸線的數學模型—傳輸線方程，作為一階雙曲型偏微分方程組，并沒有固定的求解方式，利用現有的數學手段并不能得到其精確的解析解[1-2]，但是可以利用數值計算方式，得到其較為準確的數值解。

用于求解傳輸線方程的數值方法分為兩類：頻域法和時域法。對于頻域分析方法，就是利用拉普拉斯變換將傳輸線方程轉化為純代數方程，在頻域內求得電壓、電流的頻率響應，再通過拉氏反變換得到其時域解[3-4]。

對于大規(guī)模集成電路電壓、電流的計算，頻域分析法存在計算量較大的問題。因此一般考慮在時域內求解傳輸線方程。比較常用的時域分析方法有時域有限差分法、Lax 差分法、Upwind 差分法和Lax-Wendroff 差分法等。上述方法均需要利用一定的離散方法將傳輸線方程轉換成一階擬線性偏微分方程組，然后進行計算，其推導過程較為復雜，并且在計算過程中涉及到大量狀態(tài)矩陣求逆運算，數值計算精度以及計算效率較低，并且對離散的時間步長和空間步長有一定的穩(wěn)定性要求[5-7]。

為解決傳統(tǒng)數值運算中出現的推導過程復雜、大量的狀態(tài)矩陣求逆運算以及計算時間、空間步長影響數值穩(wěn)定性等問題，本文提出利用求解一階線性常微分方程組的數值方法對單/多傳輸線的數學模型-傳輸線方程進行數值分析[8]。求解過程主要有兩個步驟：一是將傳輸線方程利用緊致有限差分法進行空間四階離散，將得到的狀態(tài)方程組利用矩陣轉換得到一階線性常微分方程組[9]；二是利用精細積分法與微分求積法(Differential quadrature method,DQM)相結合的高精度單步積分法對一階線性常微分方程組進行求解[10-11]。通過理論分析可以得出，本文算法有固定的計算格式，且沒有涉及到狀態(tài)矩陣的求逆運算，其穩(wěn)定性和計算時間步長無關等優(yōu)勢，可提高計算精度以及計算效率，最后通過仿真實例驗證本文算法的優(yōu)越性。

2 高精度單步積分法

針對一階常系數微分方程的初值問題[12]

式中，t、tf分別為時間變量和截止時間；f(t,X(t))為時間變量t對狀態(tài)變量X的函數；X0為狀態(tài)變量X在初始時間t0的函數值。

經過一系列的數值變換可將式(1)變成常系數微分方程形式

式中，H為n階定常數矩陣；r(X,t)為關于X、t的非線性函數。

在一個時間積分區(qū)間[tk,tk+1]，運用矩陣指數函數和卷積運算可得到式(3)精確的時域動態(tài)響應[13]。

式中，τ為時間積分區(qū)間[tk,tk+1]內任意時刻；Xk為時間tk時刻函數值，為給定值。Xk+1為tk+1時刻函數值，為待求值。

將式(3)等號右邊第一項用精細積分法計算[10]，可得Xk+1計算公式為

其中，T=exp[(tk+1－tk)H]=exp(ΔtH)可用精細積分法求得，等式右邊第一項計算精度可以達到計算機精度，故數值誤差主要來源于等式右邊的第二項積分，此積分又被稱為Duhamel 積分。本文采用時域微分求積法對Duhamel 積分項進行近似求解。

2.1 Duhamel 積分項計算格式

對于式(1)中一階常微分方程的初值問題，s階時域微分求積法可以表示為[11]

式中，Δt是積分時間步長；cj是網格點；mij和nj是與網格點相關的積分系數[14]。

故式(4)中Duhamel 積分項的計算公式可寫為

式中，i=1,2,…,s；t i=t k+i×Δt/s。

對于式(6)中的X?k+i/s，采用4 階Runge-Kutta進行預估，此時s=4，故i=1,2,3,4。X?k+i/s計算公式為[15]

式中，

根據指數矩陣乘法可得以下關系式

再次利用精細單步求積法計算指數矩陣T3=exp(HΔt/4)，然后利用式(9)的方法求得T1、T2。當求出Duhamel 積分項的近似值后，將其代入式(4)得到Xk+1。

尤為注意的是，上述算法不涉及到矩陣的求逆，不會因多次求逆而導致數值誤差，并且對于線性常微分方程組，式(7)的預估過程可以省略，因此在很大程度上提高了計算效率。

2.2 算法穩(wěn)定性分析

根據傅里葉方法可得迭代公式的穩(wěn)定性條件，故本文方法穩(wěn)定性條件為[16-17]

從式(10)可知只要Re( Δt?λi) ≤ 0即能滿足算法的穩(wěn)定性要求。其中λi為矩陣H的特征值；Δt為計算時間步長，是正實數。所以本文算法穩(wěn)定只需要Re(λi) ≤ 0，即時間步長Δt的取值并不影響算法的穩(wěn)定性。

3 電報方程空間離散格式

描述傳輸線的電報方程時域形式為[18]

式中，L0、G0、C0、R0分別為變壓器繞組的單位長度電感、電導、電容、電阻參數；t為時間變量；x為空間變量；u(x,t)、i(x,t)為傳輸線上關于空間、時間的電壓和電流列矢量。

為了獲得線路沿線的電壓和電流變量，將整條傳輸線均勻分成N段。假設傳輸線長度為L，則每段長度Δx=L/N，設un為x=n?Δx處的電壓值，n∈處的電流值，n∈[1,N]，如圖1 所示。

圖1 傳輸線分段模型

本文采用緊致有限差分法(Compact finite difference, CFD)對電報方程進行離散，基于CFD 的空間四階差分公式為

式中，Y(x)為u(x)或i(x)，k1= 1/24，k2= 1 -2k1。

電報方程的空間四階離散格式為

對傳輸線的首尾端采用空間二階插值公式

傳輸線首尾的離散格式為

式中，0i、ni分別為首端和尾端電流。

根據式(13)、(15)、(16)可得2N+1 個狀態(tài)方程，故可得其矩陣形式的線性常微分方程組

式中，H1、H2、s(t) 參數分別為

矩陣H1、H2中的參數分別為

將式(17)轉化為以下標準形式

式中，r(t) 為離散后的源項；，H為常數矩陣，H=-H1-1，r(t) =H1-1s(t)。

以上推導沒有結合實際線路的邊界條件，下文將給出具體線路的邊界條件。

4 仿真算例

4.1 單導體傳輸線模型

如圖2 所示，單相高壓長輸電線路系統(tǒng)在t=0.0 s 時斷路器接通，實際輸電線路空載，在空載線路首尾兩端分別連接電阻Rs和RL，電阻值為10 ??？蛰d輸電線路長度為l，假設輸電線路為均勻線路，其參數R0、L0、C0分別為輸電線路單位長度的電阻、電感和電容值。為了計算空載輸電線路的末端電壓uN，采用傳輸線方程對輸電線路系統(tǒng)進行建模。

圖2 單相高壓長輸電線路系統(tǒng)示意圖

當進行合閘操作時，線路系統(tǒng)兩端滿足以下條件

將整個線路平均分成30 段，即N= 30。系統(tǒng)輸入正弦電壓信號US，其初始相位ω0= π/2。

采用本文方法對上面空載線路的末端電壓uN進行計算，計算時間步長為Δt，并在不同的時間步長內與時域有限差分法(Finite difference time domain，FDTD)結果做對比，具體仿真結果如圖3～6 所示。

圖3 空載線路末端電壓uN 仿真結果(Δt=1 μs)

圖4 空載線路末端電壓uN 仿真結果(Δt=40 μs)

圖5 本文方法在不同時間步長下空載線路末端電壓uN 仿真結果

圖6 兩種方法在不同時間步長下空載線路末端電壓uN仿真結果(本文方法Δt=40 μs，FDTD Δt=1 μs)

由圖3～6 仿真結果可得如下結論：① 本文方法所得結果與FDTD 所得結果基本一致，可驗證本文方法的有效性；② 無論使用本文方法還是FDTD計算線路的末端電壓，都會產生一定的數值振蕩。不同的是，本文所用方法產生數值振蕩較小，計算精度較高；③ 本文方法使用于FDTD 的40 倍仿真時間步長時，本文所用方法的仿真結果基本沒有變化，但是采用FDTD 得出的計算結果，其振蕩會加劇，驗證了本文算法的穩(wěn)定性并沒有受到仿真步長的影響，以及可利用大步長進行仿真計算，提高計算效率。

4.2 多導體傳輸線模型

為擴展本文算法適用性，將本文算法用于實際大型變壓器繞組暫態(tài)過電壓的計算。

特快速暫態(tài)過電壓(Very fast transient overvoltage，VFTO)進入大型電力變壓器會對其繞組產生絕緣損壞，因此計算在VFTO 下大型電力變壓器繞組的匝間過電壓能有效預防繞組的絕緣損壞。

在不考慮變壓器組件頻變效應情況下，將變壓器繞組等效成一根根首尾相連的傳輸線，如圖7 所示，故可用描述多導體傳輸線的電報方程進行建模研究[19-20]。

圖7 變壓器繞組的多導體傳輸線模型

圖7 中，uS和uR表示每匝繞組的輸入端電壓與輸出端電壓；iS和iR表示每匝繞組的輸入端電流與輸出端電流。

大型電力變壓器繞組的多導體傳輸線模型對應的電報方程中R0、L0、C0、G0分別是單位長度的電阻、電感、電容、電導矩陣。變壓器模型及其基本參數如表1 所示。

本文所用大型電力變壓器模型是由保定天威集團提供的江西豐城SFP-800000/500 自耦變壓器，雙線圈高低高結構，自鐵心向外繞組排列：內高壓(Hi)、低壓和外高壓(Ho)；單相接線方式為YNd11 三相牽引；其鐵心形式為單相三柱式結構，如圖8 所示。其中A到X1 為外高壓繞組，采用內屏蔽連續(xù)式，上下對稱，具體線段數及段內匝數如表1 所示。A1到X為內高壓繞組，分為4 餅，每餅匝。

變壓器繞組首端接入50 ? 電阻，即RS=50 ?，尾端接地?？傻美@組的邊界條件為

式中，e0為沖擊電壓源。

沖擊電壓源e0為上升時間4 ns，峰值639 kV的沖擊電壓源，如圖9 所示。利用本文方法與FDTD計算在同一變壓器模型下A 段8 匝繞組的末端過電壓，如圖10～12 所示。

圖9 沖擊電壓源

圖10 第1 匝末端過電壓

圖11 第2 匝末端過電壓

圖12 前8 匝末端過電壓

由圖10～12 仿真結果可得如下結論。

(1) 本文方法仿真波形與FDTD 仿真波形基本一致，可驗證本文算法在大型電力變壓器繞組電磁暫態(tài)計算模型的實用性。然而由于計算過程中出現大量的矩陣求逆運算，影響了數值計算精度，FDTD仿真結果局部出現了較為輕微的數值振蕩。

(2) 從大型變壓器繞組自身分析：① 前幾匝繞組的末端過電壓的最大值出現在第一匝末端，為沖擊電壓源幅值的47%；② 隨著匝數的增加，繞組末端過電壓幅值在不斷減小。并且在每一匝仿真時間內，電壓值隨著時間的增加而逐漸衰減。

5 結論

本文提出結合精細積分法與時域微分求積法對傳輸線方程進行求解。通過對算法的推導以及仿真實例驗證，可得出以下結論。

(1) 由算法穩(wěn)定性條件可知，本文方法穩(wěn)定性條件與時間步長無關。FDTD 的計算時間步長則要滿足穩(wěn)定條件：Δt≤Δx/c(c=3×108m/s)，因此本文算法可利用大步長進行仿真計算，提高計算效率。

(2) 本文方法避免了因狀態(tài)矩陣求逆而帶來的數值誤差，抑制了數值振蕩，提高了計算精度。

(3) 首次將本文方法用于計算VFTO 下大型電力變壓器繞組的過電壓，由仿真結果可知，變壓器第一匝繞組中VFTO 的幅值最大，并隨著匝數的增加，VFTO 的幅值呈遞減狀態(tài)，在實際工程中，應注重變壓器繞組首端的絕緣保護。