分類號： 學校代碼：10079 密級：華北電力大學碩 士 學 位 論 文題 目：陡波前過電壓下頻變變壓器繞組過電壓時域算法研究英文題目：ime-domain Calculation of ransformer indingsOvervoltage under VFTO Considering Frequency-dependentParameters研究生姓名：馬宜軍 專業：電工理論與新技術研究方向：電網絡理論及其在電力系統中的應用導師姓名：梁貴書、董華英 職稱：教授2006年12月22日聲 明本人鄭重聲明：此處所提交的碩士學位論文《陡波前過電壓下頻變變壓器繞組過電壓時域算法研究，是本人在華北電力大學攻讀碩士學位期間，在導師指導下進行的研究工作和取得的研究成果據本人所知除了文中特別加以標注和致謝之處外論文中不包含其他人已經發表或撰寫過的研究成果也不包含為獲得華北電力大學或其他教育機構的學位或證書而使用過的材料與我一同工作的同志對本研究所做的任何貢獻均已在論文中作了明確的說明并表示了謝意。學位論文作者簽名： 日 期： 本人完全了解華北電力大學有關保留使用學位論文的規定即①學校有權保管、并向有關部門送交學位論文的原件與復印件②學校可以采用影印縮印或其它復制手段復制并保存學位論文③學校可允許學位論文被查閱或借閱④學校可以學術交流為目的,復制贈送和交換學位論文⑤同意學校可以用不同方式在不同媒體上發表傳播學位論文的全部或部分內容。(涉密的學位論文在解密后遵守此規定)作者簽名： 導師簽名 日 期： 日 期： 華北電力大學碩士學位論文目錄華北電力大學碩士學位論文目錄目 錄中文摘要英文摘要第一章緒論11.1陡波前過電壓（VFTO）11.1.1VFTO的產生11.1.2VFTO的波形特點與產生機理11.1.3VFTO對變壓器絕緣的影響21.2本課題國內外研究概況31.2.1變壓器模型的建立31.2.2多導體傳輸線方程的數值解法31.3本文完成的主要工作5第二章預備知識62.1多導體傳輸線的基本知識62.2變壓器分布參數模型的建立62.3變壓器線圈MTL模型基本參數的計算8第三章VFTO下變壓器線圈的空間離散模型103.1緊湊有限差分法103.2變壓器多導體傳輸線模型的空間離散123.3矢量匹配和遞歸卷積133.4稀疏矩陣技術16第四章頻變變壓器的線圈時域仿真184.1直接求解微分方程仿真電位184.1.1精細積分法184.1.2仿真與實驗驗證194.2微分方程離散為代數方程仿真電位254.2.1反向差分法254.2.2BDF應用于變壓器線圈284.2.3仿真驗證與實例計算294.3變壓器線圈餅間并聯MOV限制VFTO幅值36第五章基于krylov子空間的時域降階算法40iiii5.1基于krylov子空間的代數方程的降階405.1.1子空間簡述405.1.2Krylov子空間415.1.3Arnoldi算法425.1.4Lanczos算法445.2仿真驗證46第六章結論與展望47參考文獻48致謝52攻讀學位期間發表的學術論文和參加科研情況53華北電力大學碩士學位論文華北電力大學碩士學位論文PAGEPAGE1第一章緒論1.1陡前過電壓（VFO）1.1.1VFO氣體絕緣變電站(GasInsulatedSubstation簡稱GIS),又稱為全封閉組合電器[1]，GIS是將高壓電器放于接地的金屬殼內，以高壓氣體為主要絕緣的電站。所用的絕緣氣體主要是SF6。GIS具有結構緊湊、體積小、占地少、運行可靠、維護工作量少、對環境污染小等優點，為城市和人口稠密地區以及大型水電工程建設節省用地創造了條件，所以在電力系統中得到越來越廣泛的應用。在GIS變電站中由于斷路器隔離開關以及接地開關操作或帶電線路對地閃絡，甚至雷電波入侵，都可能會在GIS內部產生一個上升速度極快(幾到幾十ns)的電壓陡波這個電壓陡波沿著GIS管道傳播遇到波阻抗發生改變就會發生反射和折射所有在GIS中產生的多次反射和折射的各行波分量疊加在一起就會形成波頭很陡、頻率高達幾到幾十MHz甚至上百MHz數量級的陡波前過電壓，又稱特快速暫態過電壓(eryFastransientOvervoltage，簡稱VFTO)[2]。陡波前過電壓主要包括內部陡波前過電壓、瞬態外殼電壓和外部陡波前過電壓三種在GIS內高壓導體或管道和外殼之間產生的VFFO稱內部陡波過電壓出現在GIS外部及GIS以外設備上或設備內的陡波前過電壓稱之為外部陡波前過電壓；在GIS外殼和地之間產生的陡波前過電壓稱為瞬態外殼電壓[3]。1.1.2VFO的波形特點與產生機理VFTO波形包括多個頻率段，波形比較復雜，通常由四個分量組成[4]：(1)階躍電壓；(2)由于GIS內母線管道(即電暈屏蔽彎管等)波阻抗的多次微弱變化形成的極高頻范圍(最高達100MHz)；(3)由于GIS母線管道和電纜末端或架空線終端處波阻抗的顯著變化而引起的反射形成的高頻范圍(最高達30MHz)；(4)由于外部設備的大電容(如電容式互感器或輸電線載波系統的耦合電容)引起的諧振產生的低頻范圍(0.1～5MHz)。因此，內部陡波前過電壓的波形取決于GIS的內部結構和外部配置。此外，由于陡波前過電壓的行波特性，其波形隨位置不同可能有很大的變化(在某些情況下，1米的距離就會造成顯著的變化[4])通常情況下并不是因為VFTO幅值特別高而是其高頻振蕩的電壓波對變壓器的絕緣和系統的正常運行造成危害[5]。1.1.3VFO對變壓器絕緣的影響當GIS進行斷路器隔離開關等開關操作時產生的陡波前過電壓極有可能引起GIS內部設備及與GIS相連的電力系統高壓設備的絕緣事故[6]。決定該類絕緣事故的主要因素有[7,8]：(1)開關在操作期間的動態絕緣特性；(2)高壓設備的絕緣介質在VFTO作用下的絕緣耐受能力。通過對VFTO的產生機理、波形及傳播途徑的分析可知，VFTO的特性與高壓設備耐壓試驗所采用的雷電沖擊波和操作沖擊波的特性有很大不同。盡管制造商已掌握了對高壓設備(如變壓器)在標準雷電沖擊波和操作波作用下的各種絕緣數據，但是，由于VFTO的高陡度及其高頻振蕩特性，它會怎樣影響絕緣強度，絕緣的擊穿強度是多少，絕緣裕度應當怎樣選取等等問題，制造者和設計者目前還不十分清楚[9]。在各種電力設備中，變壓器受VFTO影響非常大，因為VFTO電壓波頻率高，當幅值不是特別高時，甚至不足以使變壓器入口避雷器動作的情況下，就可能侵入變壓器造成匝間絕緣的破壞因此不論是與GIS直接連接的還是非直接連接的變壓器，其內部的絕緣都會受到很大的威脅。但是當變壓器不是直接和開關相連接的時候，VFTO會有一定削弱。由陡波前過電壓引起的絕緣擊穿事故已在許多國家電網中出現，如美國電力系統(ACP)、加拿大及中國廣東核電站均發生過由于此類因素引起的超高壓變壓器絕緣擊穿事故[3]對于VFTO造成的變壓器內部的絕緣擊穿，其原因就目前的研究而言，可以歸納為以下兩個方面[10,1]：(1)在開關觸頭擊穿瞬間產生的VFTO到達變壓器時，在變壓器端部加上了一個陡波頭波，對直接連接的變壓器，其上升時間可能只有數十ns，遠遠低于雷電沖擊波試驗時的波頭上升時間(約為1.2μs)，該陡波的波頭在變壓器繞組上造成極不均勻的匝間電壓分布，大部分電壓降落在靠近入波端的一小部分線圈或導體上，電位梯度極大，危害較大。對非直接連接的變壓器，因為經過了其它的設備，陡波頭趨于平緩，其作用與雷電沖擊波相近。(2)VFTO中含有的振蕩諧波的頻率與變壓器中的若干固有振蕩頻率匹配從而引起諧振產生幅值很高的高頻諧振過電壓導致繞組與鐵芯以及匝間的絕緣破壞。由于變壓器的絕緣設計只考慮了雷電沖擊波以及操作波作用下引起的沖擊情況，故匝間絕緣在高頻振蕩情況下顯得很脆弱。1.2本題國內外研究概況1.2.1變壓器模型的建立變壓器繞組建模的方法有場的方法[12]、場-路結合的方法[13]和路的方法[14-17]三大類。但由于變壓器線圈的纏繞結構復雜，采用前兩種方法計算量龐大。將“場”問題簡化為“路”來求解，可有效地減少計算量。實驗和仿真結果的吻合表明，從工程應用的角度來說這種簡化是合理、可以接受的。因此，實際中主要采用路的方法。采用路的方法建立的模型又分為以下幾類：(1)低頻和中頻下繞組模型低頻和中頻下一般采用集中參數模型。文獻[3]基于變壓器單餅或雙餅為單元建立了等值集中電路模型，利用該模型仿真雷電波侵入變壓器的暫態過程的時候可以滿足需要。如果用此模型仿真變壓器VFTO作用下暫態電壓分布存在以下問題，首先仿真精度無法滿足要求；其次只能計算線餅首末端的過電壓值，對于線匝之間的電壓分布，不能通過該電路模型計算獲得[8]。此外，還有很多集中參數的變壓器模型，例如：EMTP下的BCTRAN和TRELEG這兩個標準的變壓器模型[14]。(2)高頻下繞組模型頻率高的時候，線匝長度相對電磁波的波長已經不能忽略，因此需要采用分布參數理論建立變壓器繞組的高頻模型[14-21]建立變壓器繞組的電路高頻模型主要有兩種方法：基于多端口網絡理論建立變壓器線圈的集中參數電路模型[14-16]；通過分割變壓器線圈為若干單元來建立詳細的內部模型[17-21]。文獻[17]采用集中參數元件構造模型在最容易損壞的前幾匝線圈,每匝線圈為一單元構造模型,其余線圈以兩餅為一單元但是文中實驗驗證的輸入波采用波前時間為μs級的雷電波,而VFTO下波前時間達到ns級；文獻[18]采用分布參數與集中參數相結合的混合模型，模型中集中電路部分的元件參數不易確定；文獻[19]對殼式變壓器線圈采用單導體傳輸線和多導體傳輸線相結合的模型，該模型僅考慮了餅內互感，未考慮餅間互感；實踐表明，該模型應用于其它類型變壓器線圈時，仿真計算結果與實際測量存在較大誤差[17]文獻[20-21]將變壓器的每一匝線圈看成一條傳輸線建立了繞組的全多導體傳輸線模型可以詳細反映變壓器線圈的電磁過程。1.2.2多導體傳輸線方程的數值解法由于VFTO頻率很高變壓器線圈參數呈現明顯的頻變效應[19]因此一般采用頻域分析方法對多導體傳輸線模型進行分析[19-21]即借助快速傅立葉變換首先求其頻域響應，然后轉換到時域。這種頻域分析方法對于大型變壓器計算量龐大，且難以處理非線性問題。為了解決這一問題，本文采用時域方法分析頻變線圈。其基本思想是直接在時域用數值方法將傳輸線的偏微分方程通過空間變量的離散，轉化為常微分方程求解。求解多導體傳輸線方程的時域數值分析方法主要有以下幾種：(1)時域有限差分法(FDTD)[22]，采用具有二階精度的中心差分近似替代偏微分方程中的空間和時間導數，即dpl)=dl

pl+Δ/2)?pl?Δ/2)+ο(Δ2)Δ其中，Δ是步長間隔并且要構造差分網格用以實現中心差分近似這是FDTD的關鍵步驟。該方法簡單直接，但是為了得到理想的精度，差分點數非常密集，降低了計算效率。(2)特征法[24]，將具有時間t和空間坐標x的偏微分方程轉換到兩組特征線族上以常微分方程的形式求解，具有較高的計算效率。該方法最基本特點是將電壓電流分解為入射波和反射波，且各具有波速Φ,這對于單根傳輸線的方程容易做到。對于單根線，傳輸線越簡單，該方法優勢越明顯，計算均勻無耗線時，其計算效率遠遠高于其他方法。當傳輸線較為復雜時，例如對于單根有耗線或非均勻線，相對優勢就減少了。當傳輸線參數高度頻變時，該方法就難以處理了。在實際的系統中，大多數都是多導體耦合線，因為特征法本質是將電壓和電流分解為入射波和反射波，因此需要對多導體耦合線解耦，對無耗耦合線容易做到。對于有耗多導體耦合線就比較復雜，電報方程右側有兩項，兩個分布參數矩陣不相同，方程組去耦時要求對兩個分布參數矩陣同時對角化，在一般情況下比較困難，只有某些特定的多導體耦合系統能做到，這就在很大程度上限制了特征法的應用范圍。(3)微分求積法(DQ)[23]，主要思想是將某點對座標的微分算子以該座標全部定義域中一系列離散點的函數值加權逼近，將微分方程化為常微分方程或代數方程求解。其主要環節在于利用一組試驗函數確定微分求積逼近式中的線性加權全系數。該方法得到的方程組所包含的微分方程數可明顯少于前述由差分而轉化成的微分方程組的方程數，但是DQ方法中試驗函數的選擇至關重要，試驗函數有很多種類，沒有任何試驗函數適用所有仿真對象，當所選試驗函數不適合仿真對象時，就會嚴重影響精度。(4)緊湊有限差分法[2],主要思想是采用四階精度的插值公式αf(x)

+αf(x)

+αf(x)

=fi12?fi121 2 1xi1

xi

xi1 x對時域電報方程空間離散化為微分方程式中f代表電壓或者電流變量α1,α2為根據采用的插值公式決定的常數。該方法離散點僅是傳統離散方法的1/3，極大地提高了計算效率。1.3本完成的主要工作本文主要研究VFTO下變壓器繞組的傳輸線高頻電路模型的時域快速求解方法，確定變壓器在VFTO作用下繞組內部的電位分布。這對于改進變壓器的絕緣設計從而提高電力系統的穩定運行將起到重要的作用。本文的具體工作如下：(1)在時域計算頻變參數是難點，本文采用矢量匹配結合遞歸卷積在處理時域頻變問題，在保證高的計算精度情況下，避免了常規卷積積分的巨大計算量。使用緊湊有限差分法在空間離散電報方程使用精細積分法求解微分方程得到了VFTO下頻變變壓器繞組的時域響應。(2)基于(1)中的理論，采用了反向差分法求解微分－代數方程組，仿真了內屏蔽式變壓器線圈VFTO下的響應，極大地減少了計算量。利用時域算法可以處理非線性的優點，仿真了MOV做變壓器的內保護時線圈的電位分布,仿真結果表明，線圈諧振電壓幅值得到了有效的抑制。(3)研究了Krylov子空間降階方法，以降低求解大規模代數方程的計算量，并且減少內存的使用仿真結果表明在保證精度的基礎上方程的階數大幅降低。第二章預備知識2.1多體傳輸線的基本知識當信號頻率提高時，導體的長度相對電磁波的波長不能忽略，就需要應用傳輸線理論[25]研究問題傳輸線是約束電磁波沿著規定方向傳輸能量和信息的系統傳輸線的幾何長度l與電磁波的工作波長λ之比值lλ稱為傳輸線的電長度，通常把lλ>.05的傳輸線稱為長線。這時傳輸線導體上存在的損耗電阻，兩導體間介質損耗產生的電導、傳輸線的自感以及兩導線間的互電容，這些量沿線分布，因而這些量稱為分布參數。在傳輸線方程的建立以及分析中，最基本的假設是所分析的傳輸線為TEM結構,即電場和磁場的方向為橫向。傳輸線的數學表達形式為電報方程，其時域形式為??u=?i?L?i??x ?t?? ??i=?u?Cu??x ?t其中，R,G,L,C分別為傳輸線單位長度的電阻電導電感電容如果初始電壓、電流值為零，頻域電報方程可寫成下列形式?d?dx?

=?(R+jωL)=?Z?d=?(G+jωC)=?dx其中Z=R+ωL為單位長度的串聯阻抗，Y=G+ωC為單位長度的并聯導納。2.2變器分布參數模型的建立為了建立變壓器線圈的多導體傳輸線(MTL,MulticondutorransmissionLines)模型，對變壓器線圈進行如下理想化假設[26]：(1)認為線圈的平均直徑遠大于其幅向尺寸(繞組的徑向寬度)因而可以忽略線匝彎曲的影響，且所有線匝的長度近似相等(取平均值，即平均匝長)。(2)認為線餅間連線及電壓源引線充分短，因而可以忽略它們對電磁場分布的影響。(3)認為線圈的平均匝長大于所分析線圈的截面尺寸，因而認為電磁波沿線匝傳播過程在同一子午面上是瞬時建立起的電磁場分布，即忽略電磁波沿軸向、徑向的延時效應。根據上面的假設，可將線圈在線端處沿子午面剖開，將線匝展成直線，每一線匝成一“傳輸線對于連續式線圈模型這些傳輸線按線圈繞制關系首尾相連(如圖2-1所示)。為了方便，以下仍將“傳輸線”稱為線匝。規定線匝按電氣聯接順序編號。當線匝連接順序改變時（線圈類型改變，線匝編號隨之改變。這樣，線匝的邊界條件為第i根線末端的電壓電流分別等于第i＋1根線首端的電壓電（i＝1，2，…，N－1；第一根線首端接電壓源，第N根線末端或接地、或接負載阻抗、或懸空，邊界條件總數為2N個。US

IS

IR)

URUS()

IS(2)

IR(2)

UR(2)US(N?)

IS(N?)

IR(N?)

UR(N?)US(N)

IS(N)

IR(N)

UR(N)圖2-1線圈首位相連示意圖沿用傳統的建立MTL方程的方法，可得?U(x,t)=?I(x,t)?LI(x,t)? x t?

(2-1)?I(x,t)=GU(x,t)?CU(x,t)x t其中U(x,t)和I(x,t)是沿線分布的N×1的電壓、電流列向量．R，L，C，G分別是單位長度上的電阻、電感、電容和電導矩陣。將上式變換到頻域可得頻域MTL方程為?dU(x,s)=?(R+sL)I(x,s)=?Z(s)I(x,s)? dx?

(2-2)?dI(x,s)=?G+s)U(x,s)=?Y(s)U(x,s)? dxZ和Y分別是N×N串聯阻抗矩陣和并聯導納矩陣。對于圖2.1所示的連續式線圈的MTL模型，其邊界條件可以描述為iri)=isi+)uri)=usi+)

i=,,",N?1us)=ust)（變壓器入口VFTO波形）上式中下標r代表線匝末端，s代表線匝首端。最后加上線圈最后一匝末端的邊界條件，根據線圈末端邊界條件的不同(可能是直接接地，懸空，經過阻抗接地等)，添加不同的邊界條件方程(可能是電壓為零電流為零或電壓和電流滿足一定的約束條件)。求解MTL上的電壓電流分布情況，通常可分為下列三個步驟[27]：(1)從靜電場的角度出發，計算傳輸線分布參數模型的單位長度參數。(2)從電路角度建立MTL方程并添加邊界條件；(3)求解第二步中在邊界條件約束下的MTL方程。2.3變器線圈MTL模型基本參數的計算(1)電容參數的確定文獻[27]中使用有限元軟件FEMM計算靜電場中變壓器線匝間儲存的能量計算時候激勵導體加1V電壓，其余導體電勢為零，可以由如下公式計算電容1N12∑ijΔij2

=Wji12其中Wj為第j個激勵導體產生的能量，i為與激勵導體之間有電容的導體編號。因為Δij值為1所以當有足夠數目的類似上述方程聯立就可以計算出匝間餅間電容，從而形成線圈的電容矩陣C。(2)電阻參數的確定計算電阻參數，必須考慮高頻下的集膚效應。高頻時單位長度電阻為[25]R=0+S其中，0表示靜態情況下的電阻，S是由于集膚效應導致的頻變電阻。fs=式中

δ(1+d2)

(2-3)δ=1

1d22(1+d2)1，d2分別為矩形導體橫截面的長和寬，σ為導體的電導率，f為對應的頻率。為了把頻域的線圈模型轉化到時域，需要對頻變電阻做如下處理，把f代入s表達式(2-3)中

=ωπs sR= f = 1s s

ωR′ωσδ(1+d2) σδ(1+d2) π其中把s=jω代入上式

′=

1 1δ(1+d2) πR=R′

ω=R′ jω=

2R′

jω=

2R′

s(1?j)s s s j

1+js 2 s(3)電感參數的確定同電阻參數類似，高頻時電感參數也是由兩部分構成的[25]。L=e+i其中，e是低頻情況下的電感，i是在高頻時由于集膚效應造成的電感。e的計算公式為L=εrC1e c2i其中，C為(1)中所計算的電容矩陣，εr為絕緣材料的相對介電常數，c為光在真空中的速度，取c=3×8/s。L其表達式如下[28]：iL=si ω其中，s為電阻的頻變部分，ω為相應的頻率。線圈的串聯阻抗為Z=R+jωL=0+S+jω(e+i)=0+jωe+S+jωi=0+jωe+S+jS把S的表達式(2-3)代入上式，可以得到Z=0+se+

2′s上式右側前兩項可以使用反拉氏變換直接化到時域，第三項通過矢量匹配和遞歸卷積處理，具體方法將在以后章節詳細敘述。(4)電導參數的確定電導G可以由下式決定[28]G=σCε其中，σ是導體電導率；ε是周圍介質介電常數，C是電容。根據文獻[25]，當頻率低于吉赫茲的時候，電導受頻變因素的影響可以忽略，VFTO信號最高在百z，所以本文不考慮電導的頻變。第三章VFTO下變壓器線圈的空間離散模型本章以連續式線圈為例，詳細敘述了VFTO下頻變線圈的空間離散。本章首先使用緊湊有限差分法(CompactFiniteDifference,CFD)[29-30]對變壓器繞組的多導體傳輸線模型進行空間離散運用矢量匹配法[31]處理頻變參數并用遞歸卷積[32-34]處理由此引起的卷積項，從而把帶有卷積積分的偏微分方程化為狀態方程；由于矢量匹配可以較為準確地逼近頻變參數緊湊有限差分法只需取傳統分段方法1/3的段數，空間離散點少；遞歸卷積只需考慮前一步的值，避免了普通卷積積分大量的計算，因此，基于本章方法的時域仿真算法不僅可以方便地處理非線性，而且可有效地減少計算量。3.1緊有限差分法變壓器多導體模型在時域的數學形式是偏微分方程組CFD就是用來空間離散這個偏微分方程組，把模型數學形式化作微分方程組。為了敘述的簡便，以單導體情況為例進行方法的說明。設單根傳輸線的長度為l，如圖3-1所示。沿空間將其等分為M段，每段長度為Δx=l

M。u(x)表示在x=(n?12(Δx)(n=,,??,M)處的電壓i(x)表示在x=n(Δx)(n=,,??,M)處的電流。i0u0

xu12

i1u212 i2

iM?2uM?32iM1uM12iMuM圖3-1傳輸線的分段傳輸線的時域電報方程為?u(x,t)=?Li(x,t)?i(x,t)? x t?

(3-1)?i(x,t)=Cu(x,t)?u(x,t)x t式中，R,G,L和C分別為傳輸線單位長度的電阻、電導、電感和電容,u和i分別為線上的電壓和電流,x為長度上空間坐標，t為時間變量。使用下式對電報方程進行空間離散[29]αf(x)

+αf(x)

+αf(x)

=fn12?

fn121 2 1

(3-2)x n1

x n

x n1 Δx[30]式中，f代表u(x)或者i(x)；α1,α2為根據采用的插值公式決定的常數

。當采用四階精度的插值公式時，α1=1/24，α2=1/12。使用式(3-2)離散式(3-1)中第一個方程，令f(x)=u(x,t)，式(3-2)可化為αu(x,t)

+αu(x,t)

+αu(x,t)

u ?u=n12 n121 2 1x n1

x n

x n1 x把式(3-1)中第一個方程等號右側代入上式等號左側，經過化簡α

+Ldn1＋α

i+Ldn＋α

+Ldn1

u ?u=－n12 n12

(3-3a)1( n1

dt) 2( n

dt) 1(

n1 dt) x同理，令f(x)=i(x,t)，可以利用式(3-2)離散方程式(3-1)中第二個方程α

du+C n+32＋α

u

du+C n12＋α

du+C n12)1 n+32 dt

2 n12 dt

1 n12 dt

(3-3b)=－n1?n(n=,??,M?2)x由于式(3-3)中，每點的值和其前后點的值相關，但是首端沒有前一點，末端沒有后一點因此空間離散方程(3-3)對傳輸線兩端不適用本文對兩端采用二階精度的離散公式[30]。對于首端有?α(Ri

+Ld0)?α(RI

+Ld1)=12?u0? 4 0 dt?

1 1 dt Δx

(3-4)?α

Gu

du+C 12)?αGu

du i?i+C 32)=1 0對于末端有

3 12 dt

1 32 dt Δx?α

(Ri

+LdM)?α(Ri

+LdM1)=uM?uM1/2? 4 M dt

1 M1 dt Δx? (3-5)??α

Gu

+C M12)?αu

u i+C M?32)=M

M13 M12 dt

1 M?32 dt x其中，α3=1?α1,α4=1/2?α1。經過上述空間離散得到2M+1個方程，其矩陣形式為Pdx+Qx+f=0dt

(3-6)其中 x=?12 u32

"uM1/2 0 1

TT"M ;Tf=[0 " 0

u0

" uM];?3 1 ?? ??1

a2 % ?? %%% 0 ?? ?? % a2 1 ?? 1 3 ?P=? ?? 4? 0 1

1 ?2 % ?? ?? % % % ?? ?? % 2? 1

1?4?3 1

?1 1 ?? ??1

2 % %% 0 ?? %%% %% ?? ?? %2 1

%% ?? 1 3

?1 1?? ?Q=?1??1% 0

4 1 0 ?1 2 % ?? ?? %% %%% ?? %% %%% ?? ?? 0 % 1

% 2

1?? ?? ?1

1 4?式中，i=C?Δx?αi，i=L?Δx?αi,i=G?Δx?αi,i=R?x?αi由于矩陣P逆矩陣存在[35]，則由方程(3-6)得

(i=1,2,3,4)。dx=x+f (3-7)P其中A=?P?Q，fP

dt P=?Pf。CFD和傳統差分方法相比，采用的點數僅為傳統方法的1/3；并且可以保證系統降階時的無源性。3.2變器多導體傳輸線模型的空間離散連續式變壓器線圈模型是由首尾相連的N根多導體傳輸線構成[27]如圖2-1所示。多導體傳輸線模型的方程為?u(x,t)=?Li(x,t)?Ri(x,t)? x t?

(3-8)?i(x,t)=?Cu(x,t)?u(x,t)x t上式只是把式(3-1)中的

R,L,G,C變為N×N矩陣，i(x)，u(x)變為N維列向量。T1 2 N1 N其中 u(x)=[1(x,u2(x)???uN1(x,uNT1 2 N1 N

，i(x)=[i(x,i(x)???i (x,i(x)T。首尾相連的邊界條件為ijM

=(j)0,

ujM

=(j)0

(j=,2???N?);運用緊湊有限差分法,由式(3-8)可得與式(3-6)相似的下列方程?d?+??+?=0 (3-9)dt上式中各量的意義與式(3-6)中相同只是，?和?均為(2M+)N維方陣，?和?均為(2M+)N維列向量。將電流邊界條件代入式(3-9)，合并同類項，消去未知數1M,2M,??,(N)M；將電壓邊界條件代入式(3-9)，把第2MN+1～(2M+)N?1個方程對應加到第MN+2~(M+)N個方程，然后劃去第2MN+1～(2M+)N?1個方程。這樣，?向量中未知數1M,u2M???，(N)M被消去。至此，式(3-9)化為2MN+1維的方程組。?中10,uNM為已知的邊界條件。上述對方程的變換，相當于對矩陣?和?把第2MN+1～(2M+)N?1列（行）對應加到第MN+2~(M+)N（行然后劃去第2MN+1～(2M+)N?1（行。列向量?和x的處理方式與?和?中對行的處理類似。由于?可逆[35]，故方程(3-9)可變為下列標準狀態方程由于?可逆，故方程(3-9)可變為下列標準狀態方程d?=??+?

(3-10)dt P式中各項的意義與方程(3-7)相同。上述處理方法同樣適用于其它類型的變壓器線圈。3.3矢匹配和遞歸卷積VFTO作用下變壓器繞組傳輸線模型的分布參數呈現明顯的頻變特性。引入矢量匹配[1](ectorFitting)和遞歸卷積[2,33](ecursiveConvolution)相結合的方法處理頻變參數下面以集膚效應導致的頻變效應[25]為例進行說明該方法對其他因素引起的參數隨頻率變化同樣適用。對頻變參數傳輸線，頻域中的電報方程為?dU(x,s)=?Z(s)I(x,s)? dx?

(3-1)?dI(x,s)=Y(s)U(x,s)? dx其中Z(s)=R(s)+sL(s)，Y(s)=G(s)+sC(s)。考慮集膚效應時，單位長度電阻可表示為式(3-1)中第一式可寫為

R(s)=0+1sdU(x,s)=?R

+R s+sL

+sL(s))I(x,s)

(3-12)(dx 0 1 0 1(使用矢量匹配，上式中 s，1(s)化為如下形式N kf(s)≈∑ i i1s?i式中ki,i分別為第i個極點和相應的留數。將上式代入式(3-12)得dU(x,s) ?R R

kR L s

kj

?I(x,s)n m∑ ∑=?? +

i + + ?

(3-13)dx ?

0 1 s?p 0 s?p?? i1 Ri

j1

j?式中n,m分別代表

s,L(s)匹配的階數，k,p分別表示相應的極點和留數。對式(3-13)取拉氏反變換u(x,t)=?Ri(x,t)?L

i(x,t)?λ?λ

(3-14)x其中1,2均為卷積項

0 0 t 1 2?λ=λ(x,t)=

n tp(tτ)k e ix, d?1 1??

∑i1

Ri0

(3-15)?λ=λ

(x,t)=

m t p(tτ)Ri ( τ)τj ( Ri ( τ)τj ( )2 2 ∑j1

j0j對式(3-1)中第二個方程可類似處理。將具有N匝的變壓器線圈使用N根耦合頻變的多導體傳輸線系統建模把每根導體分為M段，得到形式如式(3-9)，傳輸線上分布電壓和電流空間離散方程?d?+??+?+?+?=0

(3-16)dt 1 212其中，?,?是由形如(3-14)的附加項（卷積項）形成的2MN+1維列向量12意義和(3-9)中相同。對任意空間點x，式(3-15)可以做如下變化：n m1(x,t)=∑i(x,t),2=∑2j(x,t)i1若時間步長為τ，在時刻n有如下遞推關系

j1?i(x,n+)=e?

pRτ

i(x,n)+di(x,n)+di2i(x,n+)

(3-17)2j(x,n+)=e

pLτ

2j(x,n)+dji(x,n)+dj2i(x,n+)k ? epRτ?1?

k ?epRτ?1 ?其中 di1=

Ri

?eRτ?

?，di2=

Ri?

?1?pi?

pτ ?

pi?

pτ ?d =k

?pτ?eLj?

epLτ?1??，d

?epτ?1=k?

??1?,j1 Lj?

pLτ ?

j2 Lj? pτ ?? ? ? Lj ?遞歸卷積關鍵是式(3-17)，對(3-17)中一式進行證明，二式類似可證明。由式(3-15)n1 Ri( )n1 Ri( )λ(x,t

)=k

e i(x,t)dti n1 Ri0(=kRi(

t∫np(tτ)( )∫np(tτ)( )0

t∫n1p(tτ∫n1p(tτ)( ) )tn

(3-18)n1 Rin1 Ri( )

ptτ=eiλ

(x,t

)+k

e i(x,t)dtni nn令

Rit?(x,t

)=tn1eRi(tτ)i(x,t)dt

(3-19)ni n1 tn利用線性插值逼近傳輸線上分布電流關于時間t的函數，即假定電流在各時間點之間是線性的，則有如下結論i(x,t)=i(x,n)+i(x,n+)?i(x,n)(t?t)τ n

(3-20)把式(3-20)代入式(3-19)t+

ptτ?(x,t )=

e i(x,τ)dti n1

t?t+ pt?tn1n1

tn+

?tdt??ixt

+?epτ

n1t+ ptn1tet

t?tndt??ixtn1 Ri( )nn1 Ri( )n

Ri(n1) 1

? (,n) ?

Ri(n

) ? (

,n1)∫Ri?n τ ? ? n∫Ri

τ ?(3-21)使用分步積分法，式(3-21)可化為? 1?

epRτ?1?

1?epRτ?1 ?i(x,tn1)=

epRτ

? ?i(x,tn)+ ?1

?i(x,tn1)

(3-22)+ ? ? ? ?pi?

pτ ?

pi?

pτ ?將式(3-22)代入式(3-18)，即可得到式(3-17)中第一式。式(3-17)等號右邊的第三項可與方程中的同類未知量項合并，其余兩項可看作傳輸線上隨時間、空間變化的已知附加分布電源項，可與方程中的已知項合并。這樣，方程(3-16)可轉化為標準的狀態方程。3.4稀矩陣技術稀疏矩陣是指含有大量零元素，非零元素在全部矩陣元素中占比例很小的矩陣。可以考慮利用稀疏矩陣內含有大量零元素的特點，采取措施只存儲非零元素和只對非零元素進行運算減少存儲量和運算量這一類方法稱之為稀疏矩陣技術[36]。電路方程中，除了個別的，大多數是稀疏矩陣。而且隨著電路規模的增加，矩陣的階數增高稀疏程度也越大例如10個節點的電路系數矩陣的非零元素可能占50%；而100個節點的典型電路中，非零元素僅占5%左右。據統計，對多數實際電路，其系數矩陣中的非零元素數目在4n~6n之間(n為矩陣的階數)。稀疏矩陣技術應用主要優點：(1)減少存儲量只存儲非零元素節省存儲空間如果矩陣是對稱的還可以進一步節省存儲空間。(2)提高計算速度只進行非零元素的運算避免無效的零運算采用合理的數據存儲技術，使運算過程中對數據信息的分類、檢索、插入和刪除等處理盡可能方便，以提高數據檢索效率，減少運算時間和存儲量。以上兩個問題的解決是相互依存的又可能是相互矛盾的往往需要綜合權衡，在存儲量、計算時間之間采取適當的折衷方案。為了感性認識稀疏矩陣技術，在這里簡單介紹一種比較簡單的稀疏矩陣存儲方法：線性表存儲方法，對于稀疏矩陣T，可以使用3個向量α,β,γ分別存儲T中非零元素的橫坐標，縱坐標，數值。這樣就節省了大量的存儲空間。選擇存儲方案不僅考慮節省存儲空間，還要考慮存取時間，檢索、查找方便等因素，如果運算過程中，稀疏矩陣是變化的，就要考慮刪除、添加、改變非零元素的效率，即存取效率問題。本文算法中，式(3-6),(3-16)中矩陣P、Q，基本的電阻矩陣R,電容矩陣C，電導矩陣G，均為稀疏矩陣。使用稀疏矩陣存儲這些矩陣，并且涉及這些空間離散以后稀疏矩陣非零元素比例為15%~20%采用稀疏矩陣以后稀疏矩陣存儲空間節省了70%。稀疏矩陣是一個計算方法，在本文以后章節，涉及大型稀疏矩陣，均要采用稀疏矩陣技術，用來節省存儲空間和計算量。應用稀疏矩陣應注意以下問題：(1)對一個稀疏矩陣進行某些運算，例如求逆、平方等。結果很可能不再是稀疏矩陣，應立即轉為普通矩陣存儲格式。否則占用存儲空間不但不會減少，還會增加很多，并且由于存取數據的原因，計算速度也會大幅下降。(2)稀疏矩陣只對非零元素進行計算，但是某些計算，零元素也會產生數值，例如：對稀疏矩陣T,求eT,其中零元素會產生數值，應避免錯誤。第四章頻變變壓器的線圈時域仿真基于第三章的理論，本章分別使用了精細積分法和反向差分法，在時域仿真了VFTO作用下變壓器線圈的電位分布仿真與實驗結果接近證明了理論的正確性。然后對拉西瓦750kV晉東南1000kV變壓器進行了VFTO作用下的電位分布計算，對其絕緣設計進行指導。文中還計算了不同上升沿時間最大匝間電壓的變化規律，以找出VFTO對絕緣影響的最主要因素本章最后仿真了用氧化鋅非線性電阻做變壓器內保護時變壓器的電壓分布，發揮了時域方法處理非線性問題的優點，取得了良好的效果。4.1直求解微分方程仿真電位4.1.1精細積分法CFD把高頻變壓器模型數學形式化作微分方程組，而精細積分法[37,39](PreciseIntegrationMethod,PIM)則是用來求解這個微分方程組，由此得出變壓器線圈在VFTO下過電壓分布。式(3-7)的解析解為[39]xt)=et

x(0)+

ftf0e

A(t?ξ)P

ξ)dξ以τ為時間步長，則t時刻的解為τ

τ ξx(n+)=e

x(n)+0

e fPt?ξ)ξ

(4-1)在時間步長內，設非齊次項fp是線性的或可以線性化的，即可以表示為fPt)=0+1t?tn)

(4-2)其中0=

fPtn,1=

fP;tn將式(4-2)代入式(4-1)可得迭代公式x(n+)=eτ?x(n)+τ?f

(n)?+τ?f

(n+)

(4-3)2 p 2 p該算法的關鍵在于計算eτ。令T=eτ，m=2N，N=2,Δt=τ/m，則eA?Δt≈1+A?Δt+1[A?Δt2+1[A?Δt3+1[A?Δt42 3! 4!=1+a因此2T=e2

Aτ

=?e

A?t

N2N=1+a2N

2N1[ ]=1+[ ]

2N1[ ]×1+[ ]這種分解共進行N次。而

1+a]×1+a]=1+2a+a×a精細計算可用下列語句實現

for（j＝0；j<N，j++）a=2a+a×aendT=1+a這里最后一步再加單位矩陣l，能有效防止嚴重喪失有效數字。PIM的優點為(1)Δt的值極小，eA?Δt近似到了10-30絕大多數計算機已經作為舍入誤差，保證了計算的高精度；(2)eτ只需要計算一次，可以一直使用；時間步長τ可以取較大的數值，提高了計算效率。4.1.2仿真與實驗驗證為了驗證所提方法的可行性和準確性，首先對連續式變壓器線圈模型進行了計算和測量。線圈模型的基本參數見下表[27]表4-1變壓器線圈基本參數參數量參數值線餅數18每餅匝數10導體寬度(mm)6.95導體高度(mm)1.2平均匝長(m)1.4828絕緣紙厚度(mm)3.00絕緣相對介電常數2.2使用本文敘述方法，VFTO下線圈數學模型在頻域為如下電報方程[27]?dU(x,s) ?R=? +RR

s+sL

+sL

1?I(x,s)? dx

?0 1 0 1 s?? ? ???dI(x,s)=?(sC? dx

0+0

)U(x,s)其中，0,1,0,1,C0,0為由線圈模型單位長度的電氣參數常數矩陣，β為常量。合并同類項，上式可化為?dU(x,s)=?R

+sL

+?

s)I(x,s)? dx?

0 0 1(1 1 1?dI(x,s)=?(s?+G(1 1 1

)U(x,s)dx 0式中 ?＝R＋L，本文在H～0MHz頻率范圍內采用對數平均分布的104個頻點使用10階實數極點和留數進行矢量匹配，匹配前后均方根誤差為5×3。 s匹配后的幅頻特性曲線如圖4-1所示。圖4-1 s匹配前、后的幅頻特性曲線在仿真計算及試驗中，線圈首端直接輸入電壓，末端直接接地。輸入電壓波形的上升沿時間為20ns，如圖4-2。計算中每根傳輸線分為5段，時間步長為5ns。圖4-3分別給出了第2、4、6、8、10、1214、16餅末端電壓的仿真計算波形和測量波形；表4-2為測量與仿真輸出電壓峰值的對比；表4-3為測量與仿真輸出電壓的諧振頻率及對應的幅值由圖4-3表4-2和表4-3可見本文仿真的結果與實際測量結果基本吻合。誤差主要是由于變壓器線圈模型電氣參數的誤差所致。在采用完全傳輸線模型和相同計算（CPU為PⅢ533的條件下本文方法的仿真時間為1807s，比文獻[20]中頻域方法的2592s減少了30%的計算量。圖4-2仿真及實驗中輸入電壓波形(a)第2餅輸出波形 (b)第4餅輸出波形(c)第6餅輸出波形 (d)第8餅輸出波形(e)第10餅輸出波形 (f)第12餅輸出波形(g)第14餅輸出波形 (h)第16餅輸出波形圖4-3計算和測量時域結果的比較表4-2計算和測量結果的峰值比較餅號測量(電壓/V)仿真(電壓/V)2100.0199.40468.3165.84664.6260.34856.8155.491048.4451.171246.3150.951442.5641.251640.6341.05表4-3測量與仿真結果頻譜比較餅號測量仿真諧振頻率(MHz)頻譜×10-5V/Hz諧振頻率(MHz)頻譜×10-5V/Hz20.503.280.493.640.903.550.894.411.303.521.313.241.702.981.752.152.101.812.211.3810.00.8110.040.6610.81.051.10.8512.070.5812.050.5740.605.10.585.280.974.560.915.251.312.881.313.471.701.351.750.9810.000.5510.040.491.900.4912.050.3813.200.3413.120.3360.586.280.566.170.953.140.923.291.600.711.590.782.040.732.010.769.400.639.310.7110.810.5210.900.581.200.331.230.2980.608.310.587.430.971.320.911.701.302.151.312.761.701.131.771.1410.100.5410.080.4710.800.5310.950.561.900.4712.050.47100.607.150.606.940.971.810.931.981.291.301.271.201．711.261.761.1710.000.5410.130.481.890.451.940.46120.606.150.606.430.972.880.933.101.301.431.311.411.711.291.761.2410.000.5510.090.481.900.461.940.48140.605.10.585.290.974.560.915.301.312.871.323.451.701.351.770.9810.00.5510.040.491.900.4912.070.38160.605.510.585.590.974.620.914.881.312.771.322.921.701.191.741.0510.00.5310.090.491.900.4812.030.41本文還仿真了不考慮頻變時的結果，在圖4-4中，虛線是不考慮頻變時的各餅輸出電壓，實線是考慮頻變時的各餅輸出電壓。從圖中可以看到，未考慮頻變時，電壓峰值明顯比考慮頻變時高。主要原因是，考慮頻變時，由于集膚效應，電流集中在導線表層，阻抗要大的多，所以電壓波衰減比較快。(a)第2餅輸出波形 (b)第4餅輸出波形(c)第8餅輸出波形 (d)第14餅輸出波形圖4-4考慮頻變前、后計算結果的比較4.2微方程離散為代數方程仿真電位由于內屏蔽式線圈和連續式線圈邊界條件的不同，導致PIM的數值精度變差，因此引入了基于反向差分代替微分代數方程仿真內屏蔽式變壓器線圈的方法。由于內屏蔽式變壓器線圈中屏蔽線末端都是懸空的，所以屏蔽線末端邊界條件是Ib=0(pb代表屏蔽線末端的編號)。從第三章方程(3-6)可以看到，對于懸空屏蔽線附加向量f中的電壓邊界條件是未知的；如果在方程(3-6)中添加懸空線電流邊界條件，P矩陣不可逆就無法使用PIM求解狀態方程為了解決上述問題本節把方程(3-6)首先改寫為另一種形式然后使用反向差分法(ackwardDierentiationForulas,BDF)進行仿真。由于本方法不再對邊界條件限制，故可以應用到任何形式的線圈。4.2.1反向差分法(1)CFD方程的改寫方程組(3-3)、(3-4)、(3-5)可以改寫做如下形式Pdx+Qx+f=0dt方程(4-1)中各矩陣和向量的意義和(3-6)中已經不同。其中

(4-4)?0 3 1 0 ?? ?? 1

a2 % ?? %%% ?? ?? % a2 1 0 ?? 1 3 ?P=? ??# # 4 1 ?? 0 1

2 % ?? ?? % % % ?? ???0 0

% 21

1?4? 3 1

1 1 ?? ?? 1

2 % %% 0 ?? %%% 0

%% ?? ?? %2 1

%% ?? 1 3

1 1?? ?Q=??1 1 0

4 1 ?? %%

1 2 % ?? ?? %% %%% ?? %% %%% ?? ?? 0 %% %2

1?? ?? ?11

1 4?x=?u0f=[0

12"

320],

" uM1/2 uM

I0 1

T" IM ,矩陣PQ的維數是(2M+)N×(2M+3)N,也就是方程數目為(2M+)N未知數數目為(2M+3)N。注意方程(4-4)中并未添加邊界條件，因此附加向量f為零向量。對于首尾相連的線匝，有如下邊界條件ijM

=(j+)0，ujM

=(j+)0(j=,2???N?)對于末端懸空的線匝，邊界條件為添加邊界條件有兩種方法：

ipb=0。(a)使用類似第三章中處理連續線圈邊界條件的方法進行矩陣變換添加首尾相連的線匝的邊界條件；對末端懸空的線匝，寫做新的方程，列入方程組。(b)把首尾相連的線和末端懸空線均寫做新的方程，列入方程組。無論采用那種方法，都要把變壓器線圈首、末端邊界條件列方程寫入方程組。第一種方法可以節省存儲空間和計算量。本文采用第一種方法進行仿真。添加邊界條件以后方程和未知數個數均為(2M+)N+2(Npb+)其中Npb為屏蔽線的懸空端數目。(2)化微分方程為代數方程[39]對于方程(4-4)，添加邊界條件以后，微分項的系數矩陣是不可逆的，因此，需要化作代數方程求解，以避免求逆問題。對于形式如下的方程稱為微分代數方程f(x,,t)=0； (4-5)其中未知數個數和方程個數相等。相對傳統所述顯式狀態方程=h(x,t)

(4-6)避免了把方程首先轉化為(4-6)式的形式，這樣就可以避免直接求逆。假設方程(4-5)在t=tn，t=tn1，…，t=tn?k時的解x(t)均已經求得，其中步長i=ti1?ti不一定一致。使用i表示x(i)，則方程(4-5)在t=tn1時的解n1必須滿足f(n1,(n+),tn1)=0

(4-7)使用反向差分公式,在任何給定精度范圍在t=tn1時，用n1和k個過去值n,n1,",n?k1表示時間導數的現在值(tn1)如果設n1表示(tn1)的近似值則k階BDF表示為1k1n1=? ∑αixn?ig(n1)hi=0

(4-8)其中α0,α1,",αk為常數h為時間步長式(4-8)中只有n1是未知的因此表示為其函數g(n1)。將(4-8)代入式(4-7)f(n1,g(n1),tn1)=0

(4-9)上式成為以未知變量n1表示的一個代數方程組，可以是線性或者非線性的，本文中高頻變壓器線圈模型的方程是線性的，但是本文方法也可以計算非線性模型。BDF的階數k取值范圍為3-6k取6的時候精度已經非常高了足以滿足絕大多數問題的需要。本文仿真過程階數可以設定，仿真結果證明，采用4階精度時已經可以滿足精度要求當變步長時BDF系數是與步長有關的一個方程組本文中步長是固定的，式(4-8)中系數α可由如下公式求解?1 1 1 1 "

1??α0?

?0?? ?? ? ???01 2 3 "

k??α1?

?1?2?01 222

2 "

k2??α?

?0?? ??

?=??

(4-10)?# # # # " #??#??# # # # " #??#?

?#??#?? ?? ? ??01 2k

k " kk α 0? ??k? ??上式左邊系數矩陣是范德蒙德(andermonde)矩陣定步長情況下α0,α1,",αk的值與離散時間和步長無關。(3)反向差分表示反向差分公式使用以上算法，執行(4-8)式時，必須存儲前k步的值n,n1,",n?k1，文獻[40]證明存儲反向差分有利于減少舍入誤差，令則式(4-8)化為

n?in?i?n?i,i=,,",k1k1

(4-1)n1=?

∑?iΔxn?i

(4-12)hi=0式中式(4-8)化為

i?i∑αij=0kpn1=n+∑iΔxn?ip

(4-13)i1式中ii1+∑γjj1文獻[39]中證明了反向差分表示反向差分公式是等價的。PAGEPAGE284.2.2BDF應用于變壓器線圈本章4.1中得到的方程組，為了敘述簡便沒有添加頻變因素，現在使用矢量匹和配遞歸卷積，可以得到考慮頻變時變壓器線圈的數學模型。頻變的變壓器線圈可表示為??+??+?

+?=0

(4-14)t 00其中?，?為相關系數矩陣，?為未知電壓電流列向量，?為邊界條件，?為頻變0產生的附加項，式中各項和第三章中相比，意義相同，但是內容已經不同。應用改進后的BDF使用式(4-12)逼近式(4-14)中的導數項化簡后在時刻n+1,可以表示為?(

??a???

n1+1)

+??

n1+f0

+?=0

(4-15)0 1上式中，a=α/τ，?為前k部值的加權項，其余各相意義和式(4-14)中相同。0 1對式(4-15)進一步化簡? ?(a?+?)?n1=?(??1+f0+?)? ?

(4-16)顯然，方程(4-16)的解為

? 1

?? ? ?n1=?(aP+?)

?(P?1+f0+λ)

(4-17)對于方程(4-16)，如果階數比較低，可以直接求逆。如果階數比較高，本文使用了LU分解法[36]來求解。方程(4-16)可以簡寫為代數方程形式Ax=bL、U分解后有

A=LU其中，L為下三角矩陣，U為上三角矩陣。則原方程可以化為令則原方程化為

LUx=bx=y

(4-18)Ly=b (4-19)使用追趕三角法求解方程(4-19)得到y值，再使用一次追趕三角法求解方程(4-18)得到x值，即解得了原方程的解。本文使用LU分解法主要原因在于，分解后可以保證矩陣L、U的稀疏性，使用稀疏矩陣技術，節省計算量和存儲空間。4.2.3仿真驗證與實例計算(1)仿真驗證由于本節的算法同時可以求解本章第一節試驗中的連續式線圈，所以首先用上一節的實驗檢驗本節算法的正確性。仿真與測量結果比較見圖4-5。由圖可以看到仿真與測量結果在幅值與波形都很接近，造成誤差的主要原因是變壓器線圈的復雜性和建立模型做了一些簡化。相對本章第一節的方法，本節方法使矩陣維數增大，計算機資源占有增多。圖4-5仿真與測量結果比較(2)實例驗證天威集團使用變壓器線圈沖擊電壓計算軟件仿真雷電波下線圈過電壓分布經過與實測比較，該軟件結果與實際接近，精度足夠指導實際生產。該軟件對每餅線圈使用集中參數模型進行建模無法仿真VFTO高頻信號作用下的線圈中過電壓分布應用本文算法計算拉西瓦750kV主變雷電波作用下線圈過電壓分布和已有軟件結果比較，以證明本文算法的正確性。輸入波為雙指數函數疊加成的標準實驗雷電波，波前時間為1.2μs，波尾時間為50μs，峰值為100kV。即??t

?t?()ut=A?eτ1?eτ2()? ?? ?

()其中，τ1為波尾時間，τ2為波前時間，A=104kV.本文算法和軟件得到的一些餅的末端對地電壓如表4-4所示表4-4軟件和本文方法各餅峰值比較餅號軟件電壓峰值(電壓/V)本文方法電壓峰值(電壓/V)兩者誤差百分比293.5292.101.5184％487.6086.121.6895％1071.0782.2115.6747％1871.6073.192.2207％2469.8361.921.3275％3458.3655.435.0206％4433.6635.164.4563％5414.7114.014.7587％表4-4給出了前幾餅和中間兩者誤差比較大的一些餅的末端電位從表中數據可以看出在VFO最容易造成損害的前幾餅兩者非常接近中間的一些餅誤差比較大，但是O一般不會對這些餅造成損害。誤差主要原因是本文算法未考慮副邊的影響。在O作用時間比較短副邊影響小雷電波作用時間比較長副邊影響就會造成比較大的誤差。(3)上升沿時間對匝間電壓的影響由文獻[5]VFTO頻率對變壓器的匝間絕緣影響相當大本文在不改變輸入波幅值的情況下，僅改變電壓波的頻率。當頻率提高時候，上升沿時間相應變短。以拉西瓦750kV主變為例子作用在變壓器入口處的VFTO電壓保持幅值不變波前時間分別為20ns,60ns,100ns.圖4-6中給出了不同波前時間情況下最大匝間電壓的波形最大匝間電壓峰值的絕對值分別為19.2kV18.7kV16.1kV由圖可以看到，保持幅值不變上升沿時間相應變短，匝間電壓絕對值越大，對匝間絕緣影響越大。因此頻率越高，對匝間絕緣的影響也越大。圖4-6波前時間為20ns、60ns、100ns時最大匝間電壓的波形(4)實例計算本文對設計中的晉東南1000kV主變高壓線圈VFTO下電壓分布進行了計算。該變壓器高壓線圈屏蔽段示意圖如圖4-7所示。總匝數為842。按照設計要求，GIS開關投切時可能產生波頭很陡的暫態過電壓(VFTO)(暫按2.5倍設備最高運行電壓考慮)，因此施加在變壓器高壓繞組首端的暫態過電壓最大幅值按2.5*100/3＝1588kV考慮本文輸入波的波前時間為20ns是情況比較嚴重的時候，如圖4-8(a)所示。圖4-8(b)-4-8(m)給出了變壓器線圈匝間電壓，可以看出最高匝間電壓出現在第1餅1號線與屏蔽線之間峰值為384.8kV從圖中還可以看出電壓差比較大的匝主要出現在線圈前幾匝,也就是電壓降主要在前幾匝上是絕緣保護的重點而以后的線匝壓降就比較低了。圖4-8(n)-4-8(r)給出了變壓器前幾匝線圈對地電壓從圖中可以看出電壓變化比較大的匝主要出現在線圈前幾匝，這也是前幾匝線圈匝間電壓比較大的原因。圖4-8(s)-4-8(v)給出了變壓器線圈前幾餅的餅間電壓從圖中可以看出電壓差比較大的餅主要出現在線圈前幾餅，峰值為445.3kV。圖4-7 高壓線圈屏蔽段結構示意圖(a)線圈入口處VFO電壓波形波前時間為2ns,峰值為5kV

(b)最大匝間電壓峰值為348k,出現在第1線圈與屏蔽線之間(c)第1餅2號線與其前一匝屏蔽線的電壓差，峰值絕對值為224kV

(d)第1餅2號線與其后面屏蔽線的電壓差，峰值絕對值為2.5kV(e)第1餅3號線與其前面屏蔽線的電壓差，峰值絕對值為2.3kV

(f)第1餅3號線與其后面屏蔽線的電壓差，峰值絕對值為1.3kV(g)第1餅4號線與其前面屏蔽線的電壓差，峰值絕對值為1.1kV

(h)第1餅6號工作線與其前面屏蔽線的電壓差，峰值絕對值為103kV(i)第1餅6號工作線與7號工作線的電壓差，峰值絕對值為83.1kV

(j)第2餅8號線與其前面屏蔽線的電壓差，峰值絕對值為78kV(k)第3餅6號工作線與7號工作線的電壓差，峰值絕對值為88.4kV

(l)第5餅7號工作線與8號工作線的電壓差，峰值絕對值為26.1kV()第1餅3號工作與4號工線的電壓差，峰值絕對值為1.kV

(n)第1第1匝屏線的電壓電壓峰值為186kV(o)第1餅2號工作線的電壓電壓峰值為18.V

(p)第1第2匝屏線的電壓電壓峰值為148kV(q)第1餅6號工作線的電壓電壓峰值為16.V

(r)第2餅2號工作線電壓電壓峰值為10.V(s)第1餅第2餅最電壓值，位置在第1匝作線，峰值為453kV

(t)第2餅第3餅最電壓值，位置在3餅8工作線，峰值為350kV(u)第34餅間最大電壓值位置在4餅134工作線，峰值為319kV

(v)第45餅間最大電壓值位置在5餅5號工線，峰值為787kV圖4-8變壓器入口電壓、匝間電壓、對地電壓、餅間電壓4.3變器線圈餅間并聯MOV限制VFO幅值文獻[41][42]分別通過實驗和計算研究了在變壓器每餅上跨接金屬氧化物可變電阻(MentalOxidearisto,MOV)來限制變壓器內部的電壓諧振，取得了令人滿意的效果。但是文中把每餅線圈用集中參數元件建模，這種模型只適用于雷電波侵入的情況，無法滿足研究VFTO作用下的電壓分布的需要。本文使用多導體傳輸線對變壓器線圈建模，在拉西瓦水電站750kV主變高壓線圈餅間跨接MOV。仿真VFTO作用下的線圈電位分布因為MOV為非線性元件頻域方法無法處理而本文使用的時域方法可以很容易處理含非線性元件的問題。本文中的MOV是氧化鋅非線性電阻，MOV的U-I特性[43]可以用下式表示U=Iα式中α為非線性系數，取值為0.01-0.05。其中k=CH/Sα

(4-21)H為厚度，S為截面積，C為材料常數。也就是k與幾何尺寸，材料常數，非線性系數有關。MOV做變壓器線圈內保護原理可以用如圖4-9的U-I特性圖說明當作用在MOV兩端的電壓差小于閥值時，MOV工作在小電流區，如圖4-9所示的1區，流過MOV的電流小于1mA,這時MOV呈現很大的電阻，幾百兆歐級，連接MOV的支路近似于開路。在過電壓作用下，當作用在MOV兩端的電壓差超過閥值時，MOV的工作點要發生躍變，進入大電流區．如圖4-9中2區所示．在大電流區，MOV呈現很小的電阻，流過電流可達百安或千安級，并伴隨吸收過電壓能量，實際上其電壓變化不大這樣對VFTO過電壓來說。MOV可以將過電壓限制在一定的范圍內，對共振過電壓．MOV能起阻尼作用，抑制局部共振的幅值。當電流持續增大，MOV過流能力達到飽和時，電流值隨電壓變化比2區顯著變小，如圖4-9所示3區，進入3區臨界點電流大約是10kA左右。U1 2 3I圖4-9 MOV的U-I特性曲線圖4-10給出了拉西瓦水電站主變高壓線圈前幾餅線匝分布簡圖以及MOV的連接方式.該變壓器線圈詳細結構類似圖4-7中所示，也是內屏蔽式，為了清楚明了，圖4-10中略去了屏蔽線文獻[44]使用EMTP仿真了操作過電壓,在變壓器入口處的波形如圖4-1所示圖4-10變壓器線圈中MOV連接簡圖圖4-1變壓器入口處操作過電壓本文使用公式(4-21)構造了一個壓控非線性電阻(MOV)，MOV閥值電壓均為1.5kV,非線性系數取0.03。在狀態方程(4-14)中添加MOV的VCR方程，代入接MOV的線匝端部的KCL方程這樣就得到了接MOV時的狀態方程.新狀態方程未知變量中增加了三個MOV的電流為未知量。計算過程中，在未知量里把MOV電流寫在一起，以便使用矩陣分塊計算，加快計算速度。圖4-12給出了變壓器第10匝入口處對地電壓,該處是由諧振引起的對地電壓最高的線匝,由圖可見,MOV抑制諧振電壓效果明顯.由圖4-13可以看到餅間電壓差比較大主要在前四餅之間,接MOV前后相比餅間電壓有了明顯的降低電壓峰值限制在了MOV閥值電壓內.圖4-14所示為各MOV中的電流曲線,電流峰值為2.1kA，出現在編號為1的非線性電阻處.圖4-12加MOV前后第10匝入口端電壓圖4-13加避雷器前后1-2、2-3、3-4餅間電壓圖4-14避雷器上的電流第五章基于krylov子空間的時域降階算法求解大規模矩陣問題包括的線性方程組求解是科學工程計算中的重大課題。本文涉及了高頻下的變壓器線圈模型，按照第四章的方法，多導體耦合線模型最后形成一個高階的代數方程組最高達到七千多階計算量非常龐大幾乎已經到了PC機的處理極限。因此引入了階數縮減技術(OrderReduction)[45-48]，來節省存儲空間和降低計算量。最近幾年,數學領域關于階數縮減方面的研究工作取得了許多重大進展，許多成果已經應用到了工程方面，例如微電子學科。在過去一段時間里，電子和電力學科研究中也使用了很多階數縮減的方法。第一類是漸進波形估計法(WE)又稱Pade逼近[45]來實現對大型線性集中網絡的階數縮減原理是尋求一種形式簡單而特性又和原網絡函數逼近的函數。這種方法適用范圍廣，但是，它不能保證階數縮減模型的無源性第二類是基于Krylov子空間的階數縮減方法[46-50]這類方法是使用低階矩陣的特征值來等效高階矩陣的特征值可以保證系統的無源性例如Arnoldi,Lanczos算法。本文使用Arnoldi,Lanczos算法來處理本文的高階模型。求解線性方程組問題的Krylov子空間方法可以追溯到上世紀50年代初，此后Lanczos,Arnoldi,Paige,Saad[51-52]等人做了許多的工作。人們做了大量的理論分析和數值試驗,充分認識到Krylov子空間方法是求解大型線性方程組和大型矩陣特征值問題的一類最有效的方法。5.1基于krylov子空間的代數方程的降階5.1.1空間簡述為便于理解，首先簡述子空間映射在解大型線性方程中的應用。設一個線性方程組以下面的矩陣方程表示式中A∈Rn×n

為系數矩陣、x∈Rn

Ax=b (5-1)為待求變量列向量、b∈Rn為常系數列向量。令m是n維空間中m個線性無關的n維向量m=[1,v2,"vm]×mm∈RKm=span(1,v2,"vm)則稱Km是以上述m個n維向量作為基底在n維空間中展成的m維子空間

[24]利用向子空間投影的概念可求上述方程(5-1)的近似解。設x(m)表示利用子空間投影的近似解，此時可得到如下的等式m其中X(m)∈K，因此可以得到m

x(m)=b (5-2)令 (m) (m)

(Ax

(m)

?b)

⊥vj

j=,,",m

(5-3)x =mym其中V為上述子空間的一個基，式中y(m)為一個m維向量。m則可將方程(5-3)改寫為mVT(m

my

(m)

?b)