傳輸線方程高精度直接積分的數(shù)值求解方法*

張 萍 劉 寧 聶鑫鵬 吉增強(qiáng)

(1.省部共建電工裝備可靠性與智能化國(guó)家重點(diǎn)實(shí)驗(yàn)室(河北工業(yè)大學(xué)) 天津 300130；2.河北工業(yè)大學(xué)人工智能與數(shù)據(jù)科學(xué)學(xué)院 天津 300130；3.中電建(?？?新能源有限公司 ?？?570100；4.河北燃?xì)庥邢薰?石家莊 050000)

1 引言

單/多導(dǎo)體傳輸線的數(shù)學(xué)模型—傳輸線方程，作為一階雙曲型偏微分方程組，并沒有固定的求解方式，利用現(xiàn)有的數(shù)學(xué)手段并不能得到其精確的解析解[1-2]，但是可以利用數(shù)值計(jì)算方式，得到其較為準(zhǔn)確的數(shù)值解。

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

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

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

2 高精度單步積分法

針對(duì)一階常系數(shù)微分方程的初值問題[12]

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

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

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

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

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

將式(3)等號(hào)右邊第一項(xiàng)用精細(xì)積分法計(jì)算[10]，可得Xk+1計(jì)算公式為

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

2.1 Duhamel 積分項(xiàng)計(jì)算格式

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

式中，Δt是積分時(shí)間步長(zhǎng)；cj是網(wǎng)格點(diǎn)；mij和nj是與網(wǎng)格點(diǎn)相關(guān)的積分系數(shù)[14]。

故式(4)中Duhamel 積分項(xiàng)的計(jì)算公式可寫為

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

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

式中，

根據(jù)指數(shù)矩陣乘法可得以下關(guān)系式

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

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

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

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

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

3 電報(bào)方程空間離散格式

描述傳輸線的電報(bào)方程時(shí)域形式為[18]

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

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

圖1 傳輸線分段模型

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

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

電報(bào)方程的空間四階離散格式為

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

傳輸線首尾的離散格式為

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

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

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

矩陣H1、H2中的參數(shù)分別為

將式(17)轉(zhuǎn)化為以下標(biāo)準(zhǔn)形式

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

以上推導(dǎo)沒有結(jié)合實(shí)際線路的邊界條件，下文將給出具體線路的邊界條件。

4 仿真算例

4.1 單導(dǎo)體傳輸線模型

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

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

當(dāng)進(jìn)行合閘操作時(shí)，線路系統(tǒng)兩端滿足以下條件

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

采用本文方法對(duì)上面空載線路的末端電壓uN進(jìn)行計(jì)算，計(jì)算時(shí)間步長(zhǎng)為Δt，并在不同的時(shí)間步長(zhǎng)內(nèi)與時(shí)域有限差分法(Finite difference time domain，F(xiàn)DTD)結(jié)果做對(duì)比，具體仿真結(jié)果如圖3～6 所示。

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

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

圖5 本文方法在不同時(shí)間步長(zhǎng)下空載線路末端電壓uN 仿真結(jié)果

圖6 兩種方法在不同時(shí)間步長(zhǎng)下空載線路末端電壓uN仿真結(jié)果(本文方法Δt=40 μs，F(xiàn)DTD Δt=1 μs)

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

4.2 多導(dǎo)體傳輸線模型

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

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

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

圖7 變壓器繞組的多導(dǎo)體傳輸線模型

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

大型電力變壓器繞組的多導(dǎo)體傳輸線模型對(duì)應(yīng)的電報(bào)方程中R0、L0、C0、G0分別是單位長(zhǎng)度的電阻、電感、電容、電導(dǎo)矩陣。變壓器模型及其基本參數(shù)如表1 所示。

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

變壓器繞組首端接入50 ? 電阻，即RS=50 ?，尾端接地。可得繞組的邊界條件為

式中，e0為沖擊電壓源。

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

圖9 沖擊電壓源

圖10 第1 匝末端過電壓

圖11 第2 匝末端過電壓

圖12 前8 匝末端過電壓

由圖10～12 仿真結(jié)果可得如下結(jié)論。

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

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

5 結(jié)論

本文提出結(jié)合精細(xì)積分法與時(shí)域微分求積法對(duì)傳輸線方程進(jìn)行求解。通過對(duì)算法的推導(dǎo)以及仿真實(shí)例驗(yàn)證，可得出以下結(jié)論。

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

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

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