马宜军-VFTO下频变变压器绕组过电压时域算法研究

分类号： 学校代码：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)