MATLAB試驗(yàn)電力系統(tǒng)暫態(tài)穩(wěn)定分析報(bào)告

1、標(biāo)準(zhǔn)文案實(shí)驗(yàn)三電力系統(tǒng)暫態(tài)穩(wěn)定分析電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算實(shí)際上就是求解發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程的初值問題，從而得出8 -t和co-t的關(guān)系曲線。每臺(tái)發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程是兩個(gè)一階非線性的常微分方程。因此，首先介紹常微分方程的初值問題的數(shù)值解法。-、常微分方程的初值問題（一）問題及求解公式的構(gòu)造方法我們討論形如式（3-1 ）的一階微分方程的初值問題寸（x）=f（x, y）, axb金（3-1 ）y（xo） =y（）設(shè)初值問題（3-1）的解為y（x）,為了求其數(shù)值解而采取離散化方法，在求解區(qū)間 a,b上取一組節(jié)點(diǎn)a =&：:：xi ：- ： xi :二 xi 1 ： : xn =b稱hi =xrxi （

2、i =0,1，n1 ）為步長。在等步長的情況下，步長為b _ ah 二n用yi表示在節(jié)點(diǎn)xi處解的準(zhǔn)確值y（x）的近似值。設(shè)法構(gòu)造序列 ）所滿足的一個(gè)方程（稱為差分方程）yi+ =y +h 陰（xi, yi ,h）（3-2 ）作為求解公式，這是一個(gè)遞推公式，從（x, y）出發(fā)，采用步進(jìn)方式，自左相右逐步算出y（x）在所有節(jié)點(diǎn)X上的近似值yi （i=1,2，n）。在公式（3-2 ）中，為求yi+只用到前面一步的值 yi ,這種方法稱為單步法。 在公式（3-2 ）中的yi +由yi明顯表示出，稱為顯式公式。而形如（3-3）y+=yi+h 甲3，yi, yi+, h）（3-3）的公式稱為隱式公式，因

3、為其右端中中還包括y書。如果由公式求時(shí)，不止用到前一個(gè)節(jié)點(diǎn)的值，則稱為多步法。(3-4)由式（3-1 ）可得dy = f (x, y)dx兩邊在x , xf上積分，得大全標(biāo)準(zhǔn)文案xi . 1-y（Xi+） =y（x。f （x, y（x） dx （3-5）Xi由此可以看出，如果想構(gòu)造求解公式，就要對右端的積分項(xiàng)作某種數(shù)值處理。這種求解公式的構(gòu)造方法叫做數(shù)值積分法。（二）一般的初值問題的解法.歐拉法和改進(jìn)歐拉法對于初值問題（3-1）,采用數(shù)值積分法，從而得到（3-5）。對于（3-5）右端的積分用矩形公式（取左端點(diǎn)），則得到Xi 1f （x, y（x）dx : h f （Xi, y（xj Xi進(jìn)而得

4、到（3-1 ）的求解公式（3-2 ） TOC o 1-5 h z yi + =yi +h f（Xi,yi）（ i =0, 1, 2, n-1 ）（3-6）此公式稱為歐拉（Euler ）格式。如果對式（3-5）右端的積分用梯形公式X X1hf (x, y(x)dx : - (f(Xi, y(Xi) f(Xi 1, y(Xi 1)則可以得到初值問題（3-1 ）的梯形求解公式如式（3-7 ）yi+=yi +2 區(qū)）+f （為十，yq （ i =0, 1, 2, n-1 ）（3-7）式（3-7）是個(gè)隱式公式?？梢圆扇∠扔脷W拉格式求一個(gè)y（x+）的初步近似值，記作 與書，稱之為預(yù)報(bào)值，然后用預(yù)報(bào)值yi+

5、替代式（3-7）右端的yi+,再計(jì)算得到y(tǒng) ,稱之為校正值，這樣建立起來的預(yù)報(bào)-校正方法稱為改進(jìn)歐拉格式(3-8)y 1 =yi - h f (Xi, yi)yi 1 =yi h f (Xi, yi) - f (x 1 ,yi 1) I22.龍格一庫塔方法在單步法中，應(yīng)用最廣泛的是龍格庫塔（ Runge-kutta ）法，簡稱R- K法。下面直接給出一種四階的龍格一庫塔法的計(jì)算公式（3-9 ）大全標(biāo)準(zhǔn)文案1Va =yi +(Ki +2K2 +2K3 +K4) 一 6(3-9) TOC o 1-5 h z Ki =h f (Xi ,yi) h1水2 =h f(Xi +-, yi +- Ki) 2

6、2. h 1K3 =h f (x +-,yi +-K2) 22K4 =h f(Xi - h,yi - K3)它也稱為標(biāo)準(zhǔn)(古典)龍格-庫塔法。例3-1研究下列微分方程的初值問題 HYPERLINK l bookmark36 o Current Document *1C 2jy =2y1 +x y(0) =o解：這是一個(gè)特殊的微分方程，其解的解析式可以給出，為應(yīng)用龍格庫塔法，取 h =0.25 ,根據(jù)式(3-9)編寫一段程序，由零開始自左相右逐步算出y(x)在所有節(jié)點(diǎn)X上的近似值y。計(jì)算結(jié)果見表3-1。計(jì)算結(jié)果表明，四階龍格庫塔方法的精度是較高的。表3-1Xnyny(Xn ) -yn2.00.3

7、99956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-83-9)編寫較復(fù)雜的實(shí)際上，MATLAB為常微分方程提供了很好的解題指令，使得求解常微分方程變得很容易，并且能將問題及解答表現(xiàn)在圖形上。因此，我們可以不用根據(jù)式(程序，而只需應(yīng)用MATLABt供的常微分方程解題器來解決問題。下面給出用MATLA面寫的解題程序。首先編寫描述常微分方程的OD或件，文彳名為 myfun，便于解題器調(diào)用它。function dy = myfun(x,y) dy = zeros(1,1);dy=1/(1+xA2)-2*yA2;再編寫

8、利用解題器指令求解y的程序。clearx0=0;for i=1:4xm=2*i;大全標(biāo)準(zhǔn)文案y0=0;x,y = ode45( myfun ,x0 xm,y0);format longy(length(y) endplot(x,y,-)運(yùn)行上述程序，在得到幾個(gè)點(diǎn)的函數(shù)值的同時(shí)，也得到函數(shù)y的曲線，如圖3-1所示。圖3-1 根據(jù)運(yùn)算結(jié)果畫出y的曲線簡單電力系統(tǒng)的暫態(tài)穩(wěn)定性(一)物理過程分析某簡單電力系統(tǒng)如圖3-2(a)所示，正常運(yùn)行時(shí)發(fā)電機(jī)經(jīng)過變壓器和雙回線路向無限大系統(tǒng)供電。發(fā)電機(jī)用電勢E作為其等值電勢，則電勢E 與無限大系統(tǒng)間的電抗為x t=x； +xT1+xT2(3-10 )2這時(shí)發(fā)電機(jī)發(fā)

9、出的電磁功率可表示為P=EUsin6 =Pp sin6(3-11)x如果突然在一回輸電線路始端發(fā)生不對稱短路，如圖3-2(b)所示。故障期間發(fā)電機(jī)電勢E ,與無限大系統(tǒng)之間的聯(lián)系電抗為Xlx(XdXti)( Xt2)xq=(Xd +x)十(二+Xt2)+2 (3-12)2x在故障情況下發(fā)電機(jī)輸出的電磁功率為P=EUsin6 =Pp sin 6(3-13)x在短路故障發(fā)生之后，線路繼電保護(hù)裝置將迅速斷開故障線路兩端的斷路器，如圖3-2(c)所示。此時(shí)發(fā)電機(jī)電勢E與無限大系統(tǒng)間的聯(lián)系電抗為大全標(biāo)準(zhǔn)文案x = Xd XT1 Xl XT2(3-14 )發(fā)電機(jī)輸出的功率為sin :. pPmM sin、

10、.(3-15)Plll EU X IE jXd jXT1j jXLjXT2 UjXLjXT2jXLjX:(b)(c)圖3-2簡單電力系統(tǒng)及其等值電路（a）正常運(yùn)行方式及其等值電路；（b）故障情況及其等值電路;(c)故障切除后及其等值電路如果正常時(shí)發(fā)電機(jī)向無限大系統(tǒng)輸送的有功功率為P0,則原動(dòng)機(jī)輸出的機(jī)械功率Pt等P0。因此，可3-3所示。于R。假定不計(jì)故障后幾秒種之內(nèi)調(diào)速器的作用，即認(rèn)為機(jī)械功率始終保持 以得到此簡單電力系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線如圖PrnpP.Pt =P0、0Qk:;c0m圖3-3 簡單系統(tǒng)正常運(yùn)行、故障期間及故障切除后的功率特性曲線但是，實(shí)際工作需對于

11、上述簡單電力系統(tǒng)， 我們可以根據(jù)等面積定則求得極限切除角。要知道在多少時(shí)間之內(nèi)切除故障線路，也就是要知道與極限切除角對應(yīng)的極限切除時(shí)間。解決這個(gè)問題，必須求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程。（二）求解發(fā)電機(jī)的轉(zhuǎn)子運(yùn)動(dòng)方程大全標(biāo)準(zhǔn)文案求解發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)方程可以得出8-t和co-t的關(guān)系曲線。其中8 -t曲線一般稱為搖擺曲線。在上述簡單電力系統(tǒng)中故障期間的轉(zhuǎn)子運(yùn)動(dòng)方程為5：,、(1-1);，1.(3-16)d ,1 ,=(Pt -P M sin、.)dt Tj式中，二一一功率角，其單位為弧度；.一一轉(zhuǎn)子角速度，標(biāo)幺值；i 轉(zhuǎn)子的同步角速度，即缶1= 2日=314.16 ,其單位為弧度/秒；Tj 發(fā)電機(jī)的慣性

12、時(shí)間常數(shù)，其單位為秒;R、p 分別為機(jī)械和電磁功率，標(biāo)幺值。這是兩個(gè)一階的非線性常微分方程，它的起始條件是已知的，即1 Ptt =to =0,= = coo =1.0 . = = 60 = sin Pm故障切除后，由于系統(tǒng)參數(shù)改變，以致發(fā)電機(jī)功率特性發(fā)生變化，必須開始求解另一組 微分方程：d、.( -1),1.(3-17)d 1愉二T；(Pt -p m sin、)式中變量含義同前述，其中 P型也為標(biāo)幺值。這組方程的起始條件為t=tc； 0 = 0c； CO = Oc其中tc為給定的切除時(shí)間；&、a為與tc時(shí)刻對應(yīng)的6和6 ,它們可由故障期間的8 -t和CO -t的關(guān)系曲線求得(6和8都是不突變

13、的)。一般來說，在計(jì)算故障發(fā)生后幾秒種的過程 中，如果8始終不超過 180o ,而且振蕩幅值越來越小，則系統(tǒng)是暫態(tài)穩(wěn)定的。當(dāng)發(fā)電機(jī)與無限大系統(tǒng)之間發(fā)生振蕩或失去同步時(shí)，在發(fā)電機(jī)的轉(zhuǎn)子回路中， 特別是阻尼繞組中將有感應(yīng)電流而形成阻尼轉(zhuǎn)矩(也稱為異步轉(zhuǎn)矩)。當(dāng)作微小振蕩時(shí)，阻尼功率可表達(dá)為：Pd = D/ = DQ1)(3-18)式中，d稱為阻尼功率系數(shù)；4為轉(zhuǎn)子角速度的偏移量，標(biāo)幺值；CO為轉(zhuǎn)子角速度，標(biāo)幺值。阻尼功率系數(shù) d除了與發(fā)電機(jī)的參數(shù)有關(guān)外，還和原始功角、尷的振蕩頻率有關(guān)。在一般情況下它是正數(shù)。 在原始功角較小，或者定子回路中有串聯(lián)電容使定子回路總電阻相對 于總電抗較大時(shí)，D可能為負(fù)數(shù)

14、。如果考慮阻尼功率的影響，則故障后的轉(zhuǎn)子運(yùn)動(dòng)方程又可 表達(dá)為大全標(biāo)準(zhǔn)文案5：二(.：. 一1)，1(3-19 )dtPt DQ，_1)_P m sin c.dt Tj3-16 );另電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算包括兩類問題，一類是應(yīng)用數(shù)值計(jì)算法得出故障期間的曲線后, 根據(jù)曲線找到與極限切除角對應(yīng)的極限切除時(shí)間，此時(shí)只需要求解微分方程(類是已知故障切除時(shí)間，需要求出搖擺曲線來判斷系統(tǒng)的穩(wěn)定性，此時(shí)需要分段分別求解微 分方程(3-16)和(3-17)。如果考慮阻尼轉(zhuǎn)矩的影響，則此時(shí)需要分段分別求解微分方程(3-16 )和(3-19 )。三、例題例3-2 某簡單電力系統(tǒng)如圖3-4所示，取基準(zhǔn)值Sb=220M

15、VA Ub=209KV。換算后的參數(shù) 已經(jīng)標(biāo)在圖中，其中一回線的電抗Xl =0.486, Tj =8.18秒。設(shè)電力線路某一回的始端發(fā)生兩相接地短路。假定 E=常數(shù)。(1)計(jì)算保持暫態(tài)穩(wěn)定而要求的極限切除角。(2)計(jì)算極限切除時(shí)間，并且作出在 0.15秒切除故障時(shí)的8 -t曲線。GT1Lx2 =0.432x =0.486XT1 =0.138LXd =0.295Xlo =4XL圖3-4某簡單電力系統(tǒng)的接線圖解：計(jì)算系統(tǒng)正常運(yùn)行方式，決定 E 和瓦。由3-3(a)的正序網(wǎng)絡(luò)可得，此時(shí)系統(tǒng)的總電抗為x =0.295+0.138+0.243+0.122=0.798發(fā)電機(jī)的暫態(tài)電勢為：0.2 0.798

16、 21.0 0.798 2 E = (1.0 .)() =1.41,1.01.0、.=tg一口798=34.53 o1.0 0.2 0.798(2)故障后的功率特性又由3-3(b)的負(fù)序、零序網(wǎng)絡(luò)可得故障點(diǎn)的負(fù)序、零序等值電抗為(0.432 0.138) (0.243 0.122)X2 丁= =0.222一 (0.432 0.138) (0.243 0.122)大全標(biāo)準(zhǔn)文案_ 0.138(0.972 +0.122) =0 1230.138 (0.972 0.122)所以在正序網(wǎng)絡(luò)故障點(diǎn)上的附加電抗為:=0.0790.222x0.123必一0.222 +0.123于是故障時(shí)等值電路如圖3-3(c

17、)所示，則x =0.433 0.3650.433 0.365_0.079= 2.80因此，故障期間發(fā)電機(jī)的最大功率為:EUx .1.41 1.02.8= 0.504(3)故障切除后的功率特性故障切除后的等值電路如圖3-3(d)所示x -0.295 0.138 0.486 0.122=1.041此時(shí)最大功率為P MEU1.41 1.0 =1.35XI1.0410. j, 1.0，、，h -180 -sin1.35=132.20E=1.41 j 0.295 j0.138 j 0.243 j 0.122U =1.0、0 =34.53f(1)(a)j 0.432 j 0.138 j 0.243F0 =

18、1.0Q0 =0.2 j 0.122f(2)j 0.138 j0.243 4 j 0.122f(0)U =1.0 j 0.122I(b)EV.41 j 0.295 j 0.138 j 0.243,j 0.079U = 1.0(c)El.41j 0.295 j 0.138 j 0.486 j 0.122q(d)圖3-5 例題7-12的等值電路(a)正常運(yùn)行等值電路；(b)負(fù)序和零序等值電路；(c)故障時(shí)等值電路；(d)故障切除后等值電路(4)計(jì)算極限切除角大全標(biāo)準(zhǔn)文案COS CcmPt ( ：0 -、)P，m COS.：,h -Pm COS%P |MP M1.0 (132.2 -34.53) 1

19、.35cos132.2 _0.504cos34.531801.35 -0.504=0.458、.cm =62.740(5)找出極限切除時(shí)間tcm根據(jù)(3-16),首先計(jì)算初值34.536。=* n=0.6027,切。=1.0180令y(1)= S, y(2)= 0。編寫描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的OD政件，文彳名為 myequ，。function dy = myequ(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的程序。cleart

20、0=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,*)text(0.194,60, delta_cmax=62.76circ , FontSize ,10)text(0.194,56,t_cmax=0.194s, FontSize ,10)圖3-6 例題7-12的S-t曲線圖3-6給出短路發(fā)生后0秒到0.25秒期間的8 -t計(jì)算曲線，根據(jù)最大切除角6cm大全標(biāo)準(zhǔn)文案(=62.740)找到極限切除時(shí)間tcm為0.194秒。由圖3-6可見

21、，如果故障切除時(shí)間大于0.194秒，則發(fā)電機(jī)的功角將不斷地增大，最終失去暫態(tài)穩(wěn)定。在極限切除時(shí)間之前切除故障，發(fā)電機(jī)的搖擺曲線的狀況將在下面作計(jì)算、分析。(6)不考慮阻尼轉(zhuǎn)矩的影響，當(dāng)故障切除時(shí)間為 0.15秒時(shí)通過計(jì)算得出8 -t曲線首先編寫描述故障期間轉(zhuǎn)子運(yùn)動(dòng)方程的OD或件，文彳名為“ myfun01”。function dy = myfun01(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);再編寫描述故障切

22、除后轉(zhuǎn)子運(yùn)動(dòng)方程的OD或件，文彳名為“ myfun02”。function dy = myfun02(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(length(Y1),1);wc=Y1(le

23、ngth(Y1),2);T2,Y2 = ode45( myfun02 ,tc tm,dc wc); TOC o 1-5 h z plot(T1,(Y1(:,1)/pi)*180,- ,T2,(Y2(:,1)/pi)*180,- ,tc,(dc/pi)*180,*)text(0.28,50,itt_c=0.15s, FontSize ,8)text(0.28,43,itdelta_c=51.71circ, FontSize ,8)xlabel( itt)ylabel( itdelta )計(jì)算結(jié)果表明，功角6沿著故障切除后的功角特性曲線根據(jù)等面積定則作等幅振蕩，如圖3-7所示。實(shí)際上，由于阻尼轉(zhuǎn)矩的影響，振蕩的幅度是逐漸衰減的，功角 6最終運(yùn)行在6k =47.8 o 。因此，發(fā)電機(jī)能夠保持暫態(tài)穩(wěn)定。大全標(biāo)準(zhǔn)文案圖3-7不考慮阻尼轉(zhuǎn)矩