MATLAB實驗電力系統(tǒng)暫態(tài)穩(wěn)定分析Word版

整理為word格式整理為word格式整理為word格式實驗三電力系統(tǒng)暫態(tài)穩(wěn)定分析電力系統(tǒng)暫態(tài)穩(wěn)定計算實際上就是求解發(fā)電機轉(zhuǎn)子運動方程的初值問題，從而得出δ-t和ω-t的關系曲線。每臺發(fā)電機的轉(zhuǎn)子運動方程是兩個一階非線性的常微分方程。因此，首先介紹常微分方程的初值問題的數(shù)值解法。常微分方程的初值問題（一）問題及求解公式的構(gòu)造方法我們討論形如式（3-1）的一階微分方程的初值問題（3-1）設初值問題（3-1）的解為，為了求其數(shù)值解而采取離散化方法，在求解區(qū)間[]上取一組節(jié)點稱（）為步長。在等步長的情況下，步長為用表示在節(jié)點處解的準確值的近似值。設法構(gòu)造序列所滿足的一個方程（稱為差分方程）（3-2）作為求解公式，這是一個遞推公式，從（，）出發(fā)，采用步進方式，自左相右逐步算出在所有節(jié)點上的近似值（）。在公式（3-2）中，為求只用到前面一步的值，這種方法稱為單步法。在公式（3-2）中的由明顯表示出，稱為顯式公式。而形如（3-3）（3-3）的公式稱為隱式公式，因為其右端中還包括。如果由公式求時，不止用到前一個節(jié)點的值，則稱為多步法。整理為word格式整理為word格式整理為word格式由式（3-1）可得 =（3-4）兩邊在[，]上積分，得（3-5）由此可以看出，如果想構(gòu)造求解公式，就要對右端的積分項作某種數(shù)值處理。這種求解公式的構(gòu)造方法叫做數(shù)值積分法。（二）一般的初值問題的解法歐拉法和改進歐拉法對于初值問題（3-1），采用數(shù)值積分法，從而得到（3-5）。對于（3-5）右端的積分用矩形公式(取左端點)，則得到進而得到（3-1）的求解公式（3-2）（=0，1，2，n-1）（3-6）此公式稱為歐拉（Euler）格式。如果對式（3-5）右端的積分用梯形公式則可以得到初值問題（3-1）的梯形求解公式如式（3-7）（=0，1，2，n-1）（3-7）式（3-7）是個隱式公式?？梢圆扇∠扔脷W拉格式求一個的初步近似值，記作，稱之為預報值，然后用預報值替代式（3-7）右端的，再計算得到，稱之為校正值，這樣建立起來的預報－校正方法稱為改進歐拉格式（3-8）龍格—庫塔方法整理為word格式整理為word格式整理為word格式在單步法中，應用最廣泛的是龍格－庫塔（Runge-kutta）法，簡稱R－K法。下面直接給出一種四階的龍格－庫塔法的計算公式（3-9）（3-9）它也稱為標準（古典）龍格－庫塔法。例3-1研究下列微分方程的初值問題解：這是一個特殊的微分方程，其解的解析式可以給出，為應用龍格－庫塔法，取=0.25，根據(jù)式（3-9）編寫一段程序，由零開始自左相右逐步算出在所有節(jié)點上的近似值。計算結(jié)果見表3-1。計算結(jié)果表明，四階龍格－庫塔方法的精度是較高的。表3-12.00.399956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-8實際上，MATLAB為常微分方程提供了很好的解題指令，使得求解常微分方程變得很容易，并且能將問題及解答表現(xiàn)在圖形上。因此，我們可以不用根據(jù)式（3-9）編寫較復雜的程序，而只需應用MATLAB提供的常微分方程解題器來解決問題。下面給出用MATLAB編寫的解題程序。首先編寫描述常微分方程的ODE文件，文件名為ˊmyfunˊ，便于解題器調(diào)用它。functiondy=myfun(x,y)dy=zeros(1,1);dy=1/(1+x^2)-2*y^2;再編寫利用解題器指令求解y的程序。整理為word格式整理為word格式整理為word格式clearx0=0;fori=1:4xm=2*i;y0=0;[x,y]=ode45('myfun',[x0xm],[y0]);formatlongy(length(y))endplot(x,y,'-')運行上述程序，在得到幾個點的函數(shù)值的同時，也得到函數(shù)y的曲線，如圖3-1所示。圖3-1根據(jù)運算結(jié)果畫出y的曲線簡單電力系統(tǒng)的暫態(tài)穩(wěn)定性（一）物理過程分析某簡單電力系統(tǒng)如圖3-2(a)所示，正常運行時發(fā)電機經(jīng)過變壓器和雙回線路向無限大系統(tǒng)供電。發(fā)電機用電勢作為其等值電勢，則電勢與無限大系統(tǒng)間的電抗為（3-10）這時發(fā)電機發(fā)出的電磁功率可表示為（3-11）如果突然在一回輸電線路始端發(fā)生不對稱短路，如圖3-2(b)所示。故障期間發(fā)電機電勢與無限大系統(tǒng)之間的聯(lián)系電抗為（3-12）在故障情況下發(fā)電機輸出的電磁功率為整理為word格式整理為word格式整理為word格式（3-13）在短路故障發(fā)生之后，線路繼電保護裝置將迅速斷開故障線路兩端的斷路器，如圖3-2(c)所示。此時發(fā)電機電勢與無限大系統(tǒng)間的聯(lián)系電抗為（3-14）發(fā)電機輸出的功率為（3-15）圖3-2簡單電力系統(tǒng)及其等值電路（a）正常運行方式及其等值電路；（b）故障情況及其等值電路；（c）故障切除后及其等值電路如果正常時發(fā)電機向無限大系統(tǒng)輸送的有功功率為，則原動機輸出的機械功率等于。假定不計故障后幾秒種之內(nèi)調(diào)速器的作用，即認為機械功率始終保持。因此，可以得到此簡單電力系統(tǒng)正常運行、故障期間及故障切除后的功率特性曲線如圖3-3所示。圖3-3簡單系統(tǒng)正常運行、故障期間及故障切除后的功率特性曲線整理為word格式整理為word格式整理為word格式對于上述簡單電力系統(tǒng)，我們可以根據(jù)等面積定則求得極限切除角。但是，實際工作需要知道在多少時間之內(nèi)切除故障線路，也就是要知道與極限切除角對應的極限切除時間。要解決這個問題，必須求解發(fā)電機的轉(zhuǎn)子運動方程。（二）求解發(fā)電機的轉(zhuǎn)子運動方程求解發(fā)電機轉(zhuǎn)子運動方程可以得出δ-t和ω-t的關系曲線。其中δ-t曲線一般稱為搖擺曲線。在上述簡單電力系統(tǒng)中故障期間的轉(zhuǎn)子運動方程為（3-16）式中，——功率角，其單位為弧度；——轉(zhuǎn)子角速度，標幺值；——轉(zhuǎn)子的同步角速度，即==314.16，其單位為弧度/秒；——發(fā)電機的慣性時間常數(shù)，其單位為秒；、――分別為機械和電磁功率，標幺值。這是兩個一階的非線性常微分方程，它的起始條件是已知的，即==0；==1.0；==故障切除后，由于系統(tǒng)參數(shù)改變，以致發(fā)電機功率特性發(fā)生變化，必須開始求解另一組微分方程：（3-17）式中變量含義同前述，其中也為標幺值。這組方程的起始條件為=；=；=其中為給定的切除時間；、為與時刻對應的和，它們可由故障期間的δ-t和ω-t的關系曲線求得（和都是不突變的）。一般來說，在計算故障發(fā)生后幾秒種的過程中，如果δ始終不超過180o，而且振蕩幅值越來越小，則系統(tǒng)是暫態(tài)穩(wěn)定的。當發(fā)電機與無限大系統(tǒng)之間發(fā)生振蕩或失去同步時，在發(fā)電機的轉(zhuǎn)子回路中，特別是阻尼繞組中將有感應電流而形成阻尼轉(zhuǎn)矩（也稱為異步轉(zhuǎn)矩）。當作微小振蕩時，阻尼功率可表達為：整理為word格式整理為word格式整理為word格式==（3-18）式中，稱為阻尼功率系數(shù)；為轉(zhuǎn)子角速度的偏移量，標幺值；為轉(zhuǎn)子角速度，標幺值。阻尼功率系數(shù)除了與發(fā)電機的參數(shù)有關外，還和原始功角、的振蕩頻率有關。在一般情況下它是正數(shù)。在原始功角較小，或者定子回路中有串聯(lián)電容使定子回路總電阻相對于總電抗較大時，D可能為負數(shù)。如果考慮阻尼功率的影響，則故障后的轉(zhuǎn)子運動方程又可表達為（3-19）電力系統(tǒng)暫態(tài)穩(wěn)定計算包括兩類問題，一類是應用數(shù)值計算法得出故障期間的曲線后，根據(jù)曲線找到與極限切除角對應的極限切除時間，此時只需要求解微分方程（3-16）；另一類是已知故障切除時間，需要求出搖擺曲線來判斷系統(tǒng)的穩(wěn)定性，此時需要分段分別求解微分方程（3-16）和（3-17）。如果考慮阻尼轉(zhuǎn)矩的影響，則此時需要分段分別求解微分方程（3-16）和（3-19）。例題例3-2某簡單電力系統(tǒng)如圖3-4所示，取基準值=220MVA，=209KV。換算后的參數(shù)已經(jīng)標在圖中，其中一回線的電抗=0.486，=8.18秒。設電力線路某一回的始端發(fā)生兩相接地短路。假定=常數(shù)。（1）計算保持暫態(tài)穩(wěn)定而要求的極限切除角。（2）計算極限切除時間，并且作出在0.15秒切除故障時的δ-t曲線。圖3-4某簡單電力系統(tǒng)的接線圖解：計算系統(tǒng)正常運行方式，決定和。由3-3(a)的正序網(wǎng)絡可得，此時系統(tǒng)的總電抗為=0.295+0.138+0.243+0.122=0.798整理為word格式整理為word格式整理為word格式發(fā)電機的暫態(tài)電勢為：==1.41==34.53o（2）故障后的功率特性又由3-3(b)的負序、零序網(wǎng)絡可得故障點的負序、零序等值電抗為==0.222==0.123所以在正序網(wǎng)絡故障點上的附加電抗為：于是故障時等值電路如圖3-3(c)所示，則因此，故障期間發(fā)電機的最大功率為：（3）故障切除后的功率特性故障切除后的等值電路如圖3-3(d)所示此時最大功率為整理為word格式整理為word格式整理為word格式圖3-5例題7-12的等值電路（a）正常運行等值電路；（b）負序和零序等值電路；（c）故障時等值電路；（d）故障切除后等值電路（4）計算極限切除角==0.458（5）找出極限切除時間根據(jù)（3-16），首先計算初值，令y(1)=，y(2)=。編寫描述故障期間轉(zhuǎn)子運動方程的ODE文件，文件名為ˊmyequˊ。functiondy=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)));整理為word格式整理為word格式整理為word格式再編寫利用解題器指令求解y的程序。cleart0=0;tm=0.25;d0=(34.53/180)*pi;w0=1;[T,Y]=ode45('myequ',[t0tm],[d0w0]);plot(T,(Y(:,1)/pi)*180,'-',0.194,62.76,'*')text(0.194,60,'\delta_{cmax}=62.76\circ','FontSize',10)text(0.194,56,'t_{cmax}=0.194s','FontSize',10)圖3-6例題7-12的δ-t曲線圖3-6給出短路發(fā)生后0秒到0.25秒期間的δ-t計算曲線，根據(jù)最大切除角（）找到極限切除時間為0.194秒。由圖3-6可見，如果故障切除時間大于0.194秒，則發(fā)電機的功角將不斷地增大，最終失去暫態(tài)穩(wěn)定。在極限切除時間之前切除故障，發(fā)電機的搖擺曲線的狀況將在下面作計算、分析。（6）不考慮阻尼轉(zhuǎn)矩的影響，當故障切除時間為0.15秒時通過計算得出δ-t曲線首先編寫描述故障期間轉(zhuǎn)子運動方程的ODE文件，文件名為”myfun01”functiondy=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)));再編寫描述故障切除后轉(zhuǎn)子運動方程的ODE文件，文件名為”myfun02functiondy=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;整理為word格式整理為word格式整理為word格式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',[t0tc],[d0w0]);dc=Y1(length(Y1),1);wc=Y1(length(Y1),2);[T2,Y2]=ode45('myfun02',[tctm],[dcwc]);plot(T1,(Y1(:,1)/pi)*180,'-',T2,(Y2(:,1)/pi)*180,'-',tc,(dc/pi)*180,'*')text(0.28,50,'\it{t}_{c}=0.15s','FontSize',8)text(0.28,43,'\it{\delta}_{c}=51.71\circ','FontSize',8)xlabel('\it{t}')ylabel('\it{\delta}')計算結(jié)果表明，功角沿著故障切除后的功角特性曲線根據(jù)等面積定則作等幅振蕩，如圖3-7所示。實際上，由于阻尼轉(zhuǎn)矩的影響，振蕩的幅度是逐漸衰減的，功角最終運行在=47.8o。因此，發(fā)電機能夠保持暫態(tài)穩(wěn)定。圖3-7不考慮阻尼轉(zhuǎn)矩影響，當0.15秒切除故障時發(fā)電機的曲線（7）考慮阻尼轉(zhuǎn)矩的影響，當故障切除時間為0.15秒時通過計算得出δ-t曲線描述故障期間轉(zhuǎn)子運動方程的ODE文件與（6）相同，文件名也為”myfun01重新編寫描述故障切除后轉(zhuǎn)子運動方程的ODE文件，文件名為”myfun03functiondy=myfun03(t,y)f=50;w1=2*pi*f;TJ=8.18;Pt=1.0;P3m=1.35;D=15;dy=zeros(2,1);d