MATLAB实验 电力系统暂态稳定分析

实验三电力系统暂态稳定分析电力系统暂态稳定计算实际求解发电机转子运动方程初始值问题，得到-t与-t的关系曲线. 各发电机转子运动方程为两个一阶非线性常微分方程。 因此，首先介绍常微分方程初值问题的数值解法。一、常微分方程的初值问题(1)求解问题和计算公式的结构方法研究了式(3-1)一阶微分方程的初值问题(3-1)将初始值问题(3-1)的解作为为了求出其数值解而离散化的方法，在求解区间中取一系列的节点()称为步骤。 对于等步骤，步骤为使用表示节点上解的正确值的近似值。构筑顺序满足的方程式(称为差分方程式)(3-2)作为求解式，这是递归式，将各节点的近似值()从左相向右相进行阶段性计算。在式(3-2)中，为了只求出前一步骤的值，将该方法称为单步法。 式(3-2)中的数字明显表示，称为显示式。 形状如(3-3)(3-3)的公式称为隐式公式。 因为右边还包括在内。通过公式求出时，不仅是前一个节点的值，被称为多阶段法。可以由式(3-1)得到=(3-4)两侧用，积分，得到(3-5)由此可知，构筑求解式需要对右端的积分项进行某种数值处理。 求解该公式的结构方法称为数值积分法。(2)一般初始值问题的解法1 .欧拉法和改良欧拉法对于初始值问题(3-1)，通过采用数值积分法，得到(3-5)。 (3-5)用矩形式(取左端点)求出右端积分时进而，得到(3-1)求解式(3-2)(=0，1，2，n-1) (3-6)这个公式叫做Euler (Euler )格式。如果将式(3-5)右端的积分设为梯形初始值问题(3-1)的梯形求解式如式(3-7)那样得到(=0，1，2，n-1) (3-7)式(3-7)是隐式式. 以Euler形式求出一个初步的近似值，记为预报值后，将式(3-7)的右端的值替换为预报值进行计算，将其称为修正值，这样制作出的预报-修正方法可以称为改良Euler形式(3-8)2 .龙格-库塔方法单步法中使用最广泛的是龙吸虫法(runge kutta )法、简称rk法。 接着，直接给出4次伦格-库塔法的计算式(3-9)(3-9)这也称为标准(古典)龙格-库塔法。研究例3-1以下微分方程的初值问题解答：这是一个特殊的微分方程，其解的解析表达式应用龙格-库塔法，设=0.25，根据公式(3-9)编制程序，从零向左向右在所有节点上依次计算近似值。 计算结果如表3-1所示。 计算结果表明，四次龙格-库塔法精度较高。表3-12.00.399956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-8实际上，MATLAB给常微分方程式提供了很好的解题指令，使得解常微分方程式变得容易，可以用图表表示问题和解答。 因此，我们不是根据式(3-9)写复杂的程序，而是应用MATLAB提供的常微分方程式的解题器来解决问题。 用MATLAB编制的问题解决程序如下所示。首先创建描述常微分方程的ODE文件，文件名为22222222222222222222函数dy=myfun (x，y )dy=zeros (1，1 )dy=1/(1 x2)-2*y2；用解题器命令编制解y的程序。clearx0=0；for i=1:4xm=2*i；y0=0；x，y=ode45(myfun，x0 xm，y0 )format long (格式长)y (长度(y ) )结束plot(x，y，- )执行上述程序以获得若干点的函数值，同时也可以获得函数y的曲线。从图3-1运算结果绘制y曲线二、简单电力系统的暂态稳定性(1)物理过程分析某个单纯的电力系统，如图3-2(a )所示，正常运转时发电机通过变压器和双工电路向无限大系统供电。 若将发电机用电位设为其等效电位，则电位与无限大系统间电抗(3-10 )这时发电机产生的电磁力(3-11 )突然，在一次输电线的始端发生非对称短路时，如图3-2(b )所示。 故障期间发电机电位与无限大系统间的相关电抗(3-12 )故障时发电机输出的电磁功率(3-13 )发生短路故障后，线路继电保护装置如图3-2(c )所示，快速切断故障线路两端的断路器。 此时发电机电位与无限大系统的关联电抗(3-14 )发电机输出功率(3-15 )图3-2简单电力系统及其等效电路(a )正常运转方式及其等效电路(b )故障状况及其等效电路(c )故障切除后及其等效电路如果发电机在正常时向无限大系统输送有效电力，则原动机输出的机械电力相等。 假定调速器在故障后几秒内的作用，即机械功率始终保持不变。 因此，该简单的电力系统的正常运转、故障期间及故障切除后的电力特性曲线如图3-3所示。图3-3简单系统的正常运转、故障期间及故障切除后的电力特性曲线在上述简单的电力系统中，可以根据等面积法则求出极限切除角。 但是，在实际工作中，有必要知道多长时间切除故障线路，也就是说，有必要知道与极限切除角相对应的极限切除时间。 解决这个问题需要求解发电机转子运动方程。(2)求解发电机转子运动方程求解发电机转子的运动方程，得到-t和-t的关系曲线。 其中t曲线一般称为摇摆曲线。 上述简单电力系统中故障期间的转子运动方程式(3-16 )式中，的功率角为转子角速度(单位为弧度)，标度值为转子的同步角速度，即，=314.16，单位为弧度/秒的发电机的惯性时间常数为秒，其单位分别为机械、电磁力和标准值。这是两个一阶非线性常微分方程，其起始条件是已知的=0；=1.0；=故障切除后，系统参数发生变化，发电机功率特性发生变化，因此必须求出另一个微分方程(3-17 )式中变量的含义与前述相同，其中为标记或值。 这个方程的开始条件是=这里，是规定的切除时间，因为与时刻对应的和，所以能够根据故障期间的-t与-t的关系曲线求出(和均不发生突变)。 通常，在计算类型在发生故障几秒后的过程中，在不总是超过180且振幅值较小时，系统是过渡稳定的。当发电机与无限大系统之间失去振动和同步时，感应电流流过发电机的转子电路，特别是阻尼绕组，形成阻尼转矩(也称为非同步转矩)。 作为微小振动，衰减力可以表现如下=(3-18 )式中，称为衰减力系数转子角速度的偏移量、标尺值转子角速度、标尺值。 阻尼系数除了与发电机参数有关外，还与原功角、振荡频率有关。 通常是正数。 当原始功角较小或定子电路中存在串联电容器，定子电路的总电阻大于总电抗时，d可能为负数。 考虑到阻尼力的影响，故障后的转子运动方程(3-19 )在电力系统的暂态稳定计算中，使用数值计算法求出故障期间的曲线后，从曲线中找出与极限切除角对应的极限切除时间，仅求出微分方程式(3-16 )即可，另一个问题是故障切除时间已知，需要求出摇摆曲线来判断系统的稳定性，微分方程式(3-16 )和考虑到衰减转矩的影响，此时需要分别阶段性地求解微分方程式(3-16 )和(3-19 )。三、例题例3-2某个简单的电力系统如图3-4所示，设定基准值=220MVA、=209KV。 已经在图中示出了转换参数，其中，一个线路的电抗=0.486，=8.18秒。 假设在电力线的某一次始端发生了两相接地故障。 假设=常数。 (1)计算保持过渡稳定所要求的极限切除角。 (2)计算极限切除时间，建立0.15秒切除故障时的t曲线。图3-4某简单电力系统的配线图解：计算系统正常运行的方法，确定和。可从3-3(a )的顺序网络得到，在这种情况下，系统的总电抗为=0.295 0.138 0.243 0.122=0.798发电机的过渡电位=1.41=34.53(2)故障后的功率特性另外，3-3(b )反相、零相网络能够得到故障点的反相、零相等电抗=0.222=0.123正序网络故障点的附加电抗是因此，故障时等效电路如图3-3(c )所示因此，故障发电机的最大功率如下：(3)故障切除后的电力特性故障切除后的等效电路如图3-3(d )所示这种情况下最大功率图3-5例题7-12的等效电路(a )正常动作等效电路(b )反相和零相的等效电路(c )故障时等效电路(d )故障切除后的等效电路(4)计算极限切除角=0.458(5)找到极限切除时间根据(3-16 )，首先计算初始值，y (1)=-y (2)=。 制作记述故障期间中的转子运动方程式的ODE文件，文件名为myequ。函数dy=my equ (t，y )dy=zeros (2，1 )f=50； w1=2*pi*f；dy(1)=(y(2)-1)*w1；dy (2)=(1/8.18 ) * (1.0-0.504 * sin (y (1) ) )；用解题器命令编制解y的程序。cleart0=0； tm=0.25；d0=(34.53/180)*pi； w0=1；T，Y=ode45(myequ，t0 tm，d0 w0 )。plot (t ) (y (:1 )/pi ) * 180，-，0.194，62.76，* .等)文本(0.194，60，delta_cmax=62.76circ，FontSize，10 )text (0. 194，56，t_cmax=0.194s，FontSize，10 )图3-6例题7-12的-t曲线图3-6表示短路发生后0秒到0.25秒之间的-t计算曲线，从最大切除角()发现极限切除时间为0.194秒。 由图3-6可知，故障切除时间超过0.194秒时，发电机的工作角度越来越大，最终失去过渡性稳定。 到临界切除时间切除故障，计算分析发电机摆动曲线的情况如下。(6)不考虑阻尼转矩影响，故障切除时间为0.15秒时，通过计算求出-t曲线首先，创建描述故障中转子运动方程式的ODE文件，文件名为“myfun01”。函数dy=myfun 01 (t，y )f=50； w1=2*pi*f；TJ=8.18； Pt=1.0； P2m=0.504；dy=zeros (2，1 )dy(1)=(y(2)-1)*w1；dy(2)=(1/TJ)*(Pt-P2m*sin(y(1) ) )；创建描述故障切除后转子运动方程的ODE文件。 文件名为“myfun02”。函数dy=myfun 02 (t，y )f=50； w1=2*pi*f；TJ=8.18； Pt=1.0； P3m=1.35；dy=zeros (2，1 )dy(1)=(y(2)-1)*w1；dy(2)=(1/TJ)*(Pt-P3m*sin(y(1) ) )；用解题器命令创建解y的小程序。cleart0=0； tc=0.15； tm=2.0；d0=(34.53/180)*pi； w0=1.0；T1，Y1=ode45(myfun01，t0 tc，d0 w0 )DC=y1(长度(y1)、1 )WC=y1 (长度(y1 )、2 )T2，Y2=ode45(myfun02，tc tm，dc wc )plot(T1 )、(y1 (:1 )/pi ) * 180、-、T2、(y2 (:1 )/pi ) * 180、-、tc、(dc/pi)*180、*text (0. 28，50it t _ c =0. 15 s，FontSize，8 )文本(0. 28，43it delta _ c =51.71CIRC，FontSize，8 )xlabel(itt )标签(it delta )计算结果表明，功角沿故障切除后的功角特性曲线按等面积规律等幅振动。 实际上，由于衰减扭矩的影响，振幅逐渐衰减，功角最终以=47.8运行。 因此，发电机能暂时保持稳定状态。图3-7不考虑衰减转矩的影响，在0.15秒内切除了故障时的发电机曲线(7)考虑到阻尼转矩的影响，故障切除时间为0.15秒时，通过计算求出-t曲线描述故障中转子运动方程式的ODE文件与(6)相同，文件名也为“myfun01”。重写记述了故障切除后的转子运动方程式的