故障切除時(shí)間與三相快速重合閘對(duì)電力系統(tǒng)穩(wěn)定性的影響仿真教材

1、影響電力系統(tǒng)暫態(tài)穩(wěn)定性的因素i研究對(duì)象ii系統(tǒng)模型研究影響系統(tǒng)暫態(tài)穩(wěn)定性的因素，首先要確定一個(gè)系統(tǒng)模型，本例選取參數(shù)可靠的 美國(guó)西部電網(wǎng)等值模型SWCC-9,該系統(tǒng)為三機(jī)九節(jié)點(diǎn)系統(tǒng)，如圖1.1所示100MW18kV 18/2300.0085 + j0.072B/2 = j0.07459 j0.05853230/13.813.8kV35MVar:0.0119 + j0.1008B/2 = j0.10458230kVlol o J + C2D o230kVnuj + 930 o6125MW90MW50MVar0080 nur 一一 2fe30MVar230kVOTD o16.5kV1圖1.1 WS

2、CC-9系統(tǒng)模型1.2系統(tǒng)參數(shù)1.2.1節(jié)點(diǎn)參數(shù)表1.1節(jié)點(diǎn)已知參數(shù)節(jié)點(diǎn)類型電壓幅值電壓角度發(fā)電機(jī)有功發(fā)電機(jī)無(wú)功負(fù)荷有功負(fù)荷無(wú)功1V 91.04000.71640.2705002PV1.0251.63000.0665003PV1.0250.8500-0.1086004PQ00005PQ001.25000.50006PQ000.90000.30007PQ00008PQ001.00000.35009PQ0000上表中發(fā)電機(jī)有功、無(wú)功出力和負(fù)荷的有功無(wú)功功率均為以 100MVA為基準(zhǔn)時(shí)的標(biāo)幺值。1.2.2支路參數(shù)表1.2支路參數(shù)首節(jié)點(diǎn)末節(jié)點(diǎn)電阻電抗電納一半450.01000.08500.088046

3、0.01700.09200.0790570.03200.16100.1530690.03900.17000.1790780.00850.07200.0745890.01190.10080.1045140.00000.05760.0000270.00000.06250.0000390.00000.05860.0000上表中所有的參數(shù)均為標(biāo)幺值，對(duì)于變壓器支路。最后三行表示三臺(tái)變壓器參數(shù), 已經(jīng)計(jì)算出變壓器的等效電抗并直接在表格中給出。1.2.3發(fā)電機(jī)參數(shù)對(duì)于發(fā)電機(jī)，采用二階經(jīng)典模型，并對(duì)系統(tǒng)作如下假設(shè):(1) 輸入的機(jī)械功率保持恒定；(2) 忽略阻尼效應(yīng)；(3) 負(fù)荷采用恒阻抗模型。表1.3發(fā)電

4、機(jī)參數(shù)發(fā)電機(jī)XdXH10.14600.060823.6420.89580.11986.4031.31250.18133.01以上阻抗參數(shù)均以標(biāo)幺值表示，額定轉(zhuǎn)速下存儲(chǔ)的能量（H）也轉(zhuǎn)化為以100MVA為基準(zhǔn)的標(biāo)幺值。2待研究的影響因素及其仿真流程設(shè)計(jì)2.1故障切除時(shí)間的影響由等面積法則可在理論上分析得出如下結(jié)論：系統(tǒng)發(fā)生短路故障后，故障切除的越快，越有利于提高系統(tǒng)的暫態(tài)穩(wěn)定性，且存在臨界切除時(shí)間。假設(shè)系統(tǒng)0時(shí)刻以前處于穩(wěn)定運(yùn)行狀態(tài)，0時(shí)刻在線路5-7上靠近7的母線出口處 發(fā)生三相短路，故障切除時(shí)間為tc，仿真三臺(tái)發(fā)電機(jī)的轉(zhuǎn)速與功角變化曲線，通過(guò)觀察 曲線判斷系統(tǒng)是否失去穩(wěn)定。不斷改變故障時(shí)間，

5、求出恰好能使系統(tǒng)維持穩(wěn)定的故障切 除時(shí)間，該時(shí)間就是臨界切除時(shí)間。2.2三相快速重合閘的影響由等面積法則分析可知，在系統(tǒng)發(fā)生短路故障后，在減速過(guò)程中自動(dòng)重合閘動(dòng)作將 會(huì)增大加速面積，可能使原來(lái)無(wú)法維持暫態(tài)穩(wěn)定的系統(tǒng)保持穩(wěn)定。但是重合閘要求的時(shí) 間比較苛刻，如果重合閘在系統(tǒng)已經(jīng)失去穩(wěn)定才動(dòng)作則對(duì)系統(tǒng)的暫態(tài)穩(wěn)定起不到作用， 因此要求時(shí)間通常很短。這時(shí)如果采用單相重合閘由于有潛供電流的影響，單相重合閘 的動(dòng)作時(shí)間不能太快，因此本例只考慮三相重合閘。在2.1的基礎(chǔ)上，tar時(shí)刻自動(dòng)重合閘動(dòng)作，調(diào)整參數(shù)，做對(duì)比仿真實(shí)驗(yàn)，一次不投入重合閘，另一次投入重合閘，通過(guò)判 斷系統(tǒng)是否失去穩(wěn)定來(lái)得出重合閘對(duì)系統(tǒng)暫態(tài)

6、穩(wěn)定性是否有影響。然后，固定故障切除 時(shí)間tc為某一個(gè)定值，改變tar的值并觀察系統(tǒng)是否失去穩(wěn)定來(lái)確定自動(dòng)重合閘的最長(zhǎng)動(dòng) 作延時(shí)。2.3仿真流程設(shè)計(jì)首先要求出系統(tǒng)的穩(wěn)態(tài)運(yùn)行參數(shù)，即系統(tǒng)潮流分布，這些參數(shù)作為暫態(tài)過(guò)程的初始值。潮流計(jì)算采用Newton-Raphson迭代法，求出各個(gè)母線的電壓和角度。2.3.1發(fā)電機(jī)初態(tài)在動(dòng)態(tài)分析時(shí)，首先應(yīng)將發(fā)電機(jī)和負(fù)荷用相應(yīng)的模型等效。本文將發(fā)電機(jī)等效為二階經(jīng)典模型，將負(fù)荷等效為恒阻抗負(fù)荷。本系統(tǒng)中具有三臺(tái)發(fā)電機(jī)，因此會(huì)引入三個(gè)內(nèi)節(jié)點(diǎn)。因此節(jié)點(diǎn)導(dǎo)納矩陣丫將增廣到12階。可表示為:丫12 12Y 3 3Y 93Y39Y9 9在 Y33 中:在 Y39 和 Y93

7、中:在 Y99 中:Yii 二丄,I =1,2,3JXdiIIYi 丄,I =1,2,3JXdi1Yii =Yi,i =1,2,3Jxdi丫 | =丫|論怕皿8(2.1)(2.2)(2.3)(2.4)將網(wǎng)絡(luò)等效在發(fā)電機(jī)內(nèi)節(jié)點(diǎn)構(gòu)成的網(wǎng)絡(luò)中，可以得到一個(gè)3階的降階節(jié)點(diǎn)導(dǎo)納矩陣, 該矩陣可由以下關(guān)系解出:對(duì)于故障前和故障后的降階節(jié)點(diǎn)導(dǎo)納矩陣，利用式 （2.5）可以計(jì)算出對(duì)應(yīng)的降階節(jié)點(diǎn) 導(dǎo)納矩陣；對(duì)于故障中的節(jié)點(diǎn)導(dǎo)納矩陣，在 12階的增廣矩陣中去掉故障母線所在的那 一行和那一列，利用下式計(jì)算降階矩陣。丫 =丫3 3 - 丫 38 丫代 丫8 3(2.6)發(fā)電機(jī)初態(tài)中包含初始的電壓幅值和功角，根據(jù)圖2.

8、1，可以求解發(fā)電機(jī)1、2、3的初始狀態(tài)。LEi R + jQ * ，二Vi j( 11)Xdi,i =1,2,3V(2.7)LEi 二 、i(2.8)2.3.2列寫(xiě)發(fā)電機(jī)動(dòng)態(tài)方程發(fā)電機(jī)采用經(jīng)典模型，其動(dòng)態(tài)方程為:d idtdtPm-巳r2H,21,2,3(2.9)二臺(tái)發(fā)電機(jī)共有六個(gè)狀態(tài)量，定義:得到:(2.10)2.3.3求解發(fā)電機(jī)狀態(tài)利用龍格一一庫(kù)塔數(shù)值積分方法，對(duì)發(fā)電機(jī)狀態(tài)方程積分，求解發(fā)電機(jī)的狀態(tài)。認(rèn)為在故障過(guò)程中，發(fā)電機(jī)內(nèi)電壓的幅值是不變的，只有功角改變；認(rèn)為發(fā)電機(jī)的機(jī)械功率是不變的，其值等于故障前的發(fā)電機(jī)電磁功率3仿真過(guò)程與結(jié)果3.1潮流計(jì)算結(jié)果表3.1潮流計(jì)算結(jié)果節(jié)點(diǎn)類型電壓幅值電

9、壓角度發(fā)電機(jī)有功發(fā)電機(jī)無(wú)功負(fù)荷有功負(fù)荷無(wú)功1V 91.04000.00000.71640.2705002PV1.02509.28001.63000.0665003PV1.02504.66480.8500-0.1086004PQ1.0258-2.216800005PQ0.9956-3.9888001.25000.50006PQ1.0127-3.6874000.90000.30007PQ1.02583.719700008PQ1.01590.7275001.00000.35009PQ1.03241.966700003.2發(fā)電機(jī)初態(tài)表3.2發(fā)電機(jī)初態(tài)發(fā)電機(jī)內(nèi)電壓功角11.05662.271621.05

10、0219.731631.017013.16643.3故障切除時(shí)間影響暫態(tài)穩(wěn)定性仿真假設(shè)0時(shí)刻發(fā)生故障，經(jīng)過(guò)5工頻周期后切除故障，即tc = 0.083 s,重合閘不投入得到三臺(tái)發(fā)電機(jī)的轉(zhuǎn)速和功角變化曲線。3863843823803783760圖3.1發(fā)電機(jī)轉(zhuǎn)速變化曲線140012001000800600400200圖3.2發(fā)電機(jī)功角變化曲線圖3.1和圖3.2分別為該過(guò)程中三臺(tái)發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真總時(shí)間為4 s。從曲線中可以看出系統(tǒng)沒(méi)有失去暫態(tài)穩(wěn)定?，F(xiàn)在改變故障切除時(shí)間tc為12個(gè)工頻周期，即0.2 s,三臺(tái)發(fā)電機(jī)的動(dòng)態(tài)曲線分別如下。圖3.3發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.4發(fā)電機(jī)功

11、角變化曲線圖3.3和圖3.4分別為該過(guò)程中三臺(tái)發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真 總時(shí)間為4 s。從曲線中可以明顯看出系統(tǒng)已經(jīng)失去暫態(tài)穩(wěn)定。因此可以得出故障切除時(shí)間會(huì)影響系統(tǒng)的暫態(tài)穩(wěn)定性，并且時(shí)間越短越有利于提高系統(tǒng)的暫態(tài)穩(wěn)定性。通過(guò)改變故障切除時(shí)間，可以求出臨界故障切除時(shí)間tcc = 0.16 s，這時(shí)的發(fā)電機(jī)暫態(tài)曲線如下：4000.511.522.533.54t (s)5933958338537圖3.5臨界切除時(shí)間下的發(fā)電機(jī)轉(zhuǎn)速變化曲線O502O002O50000500000.511.52t (s)2.533.543.4重合閘影響暫態(tài)穩(wěn)定性仿真0050.511.522.533.54t

12、 (s)4204800044604043803600圖3.7 tc = 0.17 s未投入重合閘時(shí)發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.8 tc = 0.17 s未投入重合閘時(shí)發(fā)電機(jī)功角變化曲線取故障切除時(shí)間tc = 0.17 s,由3.3節(jié)可知，該時(shí)間大于臨界故障切除時(shí)間，因此暫 態(tài)穩(wěn)定不能維持，發(fā)電機(jī)暫態(tài)曲線如圖3.7和圖3.8所示。圖3.9 tc = 0.17 s投入重合閘時(shí)發(fā)電機(jī)轉(zhuǎn)速變化曲線4O1L圖3.10 tc = 0.17 s投入重合閘時(shí)發(fā)電機(jī)功角變化曲線22現(xiàn)將自動(dòng)重合閘投入，重合閘動(dòng)作延時(shí)為0.25 S,即動(dòng)作時(shí)刻為tar = 0.42 S,這時(shí)發(fā) 電機(jī)的轉(zhuǎn)速和功角變化曲線分別如圖 3.9

13、和圖3.10所示。從圖中課明顯看出此時(shí)系統(tǒng)沒(méi) 有失去暫態(tài)穩(wěn)定性，說(shuō)明三相快速重合閘能夠提高系統(tǒng)的暫態(tài)穩(wěn)定性。MATLAB仿真程序代碼:%生成節(jié)點(diǎn)導(dǎo)納矩陣Ypofl = zeros ( 9);for n =1:9Ypofl(BRANCH (n, 1), BRANCH (n ,1) = Ypofl ( BRANCH( n,1), BRANCH (n , 1) +1j * BRANCH(n,5) + 1/( BRANCH(n ,3) + 1 j * BRANCH(n ,4);Ypofl( BRANCH (n ,2), BRANCH (n ,2) = Ypofl ( BRANCH( n,2), BRA

14、NCH (n ,2) +1j * BRANCH(n,5) + 1/( BRANCH(n ,3) + 1 j * BRANCH(n ,4);Ypofl(BRANCH( n, 1), BRANCH( n, 2)=Ypofl(BRANCH( n, 1),BRANCH(n ,2)-1/( BRANCH(n,3) + 1 j * BRANCH(n ,4);Ypofl(BRANCH( n, 2), BRANCH( n, 1)=Ypofl(BRANCH( n, 2),BRANCH( n , 1)-1/( BRANCH(n,3) + 1 j * BRANCH(n ,4);endG = real(Ypofl )

15、;B = imag(Ypofl );tolera nee=1e-10 ;tol= 1;LoopCou nt=0;delta_y=zeros ( 14 , 1);H = zeros(8);N = zeros(8, 6);M = zeros(6, 8);L = zeros(6);%計(jì)算潮流while tol tolera nee% Delta_yfor mi=2:9sum=0;for n =1:9sum=sum + BUS (m3) * BUS (n ,3) *(G( m, n) * eos(BUS( m 4) - BUS ( n ,4)endfor m4: 9sumfor0;= 1:9sum-B

16、 ( m, n) *endcos=sum + BUS( m, 3) * BUS(BUS( m 4) - BUS ( n ,4);n, 3) *(G(m, n) * sin ( BUS( m 4) - BUS(n ,4)delta_yend(m+5) = BUS ( m 6) - BUSsum ;% Hfor m = 2: 9for n =2:9if m = nsum= 0;for k = 1:9sum= sum + BUS( m 3) * BUS( k, 3) * ( G( m, k) * sin ( BUS( m, 4) - BUS( k, 4)-B ( m k) * cos ( BUS(

17、m 4) - BUS ( k, 4);endH( m- 1 , m- 1)= BUS ( m, 3)A 2 * B ( m, m) + sum ;BUS( n, 4)elseH(m- 1 , n- 1)=-BUS ( m, 3) * BUS ( n,3) * ( G( m n) * sin ( BUS( m, 4)BUS( n, 4)-B(m, n) * cos(BUS( m, 4) - BUS ( n, 4);endendend% Nfor m =2:9forn= 4:9ifm = n%行列不等時(shí)N(m- 1 , n- 3)=-BUS ( m, 3) * BUS ( n,3) * ( G(

18、m n) * cos ( BUS( m, 4)BUS( n, 4)+ B(m, n) * sin(BUS( m, 4) - BUS ( n, 4);endendendfor m =2:7sum=0;forn= 1:9sum=sum+ BUS (m3) * BUS (n ,3) * ( G( m, n)cos ( BUS( m, 4)endfor m4: 9for2: 9ifm = nsumfor k0;= 1:9sumBUS( k, 4)(m, k) *endsin=sum + BUS ( m, 3) * BUS ( k, 3)(BUS( m, 4) - BUS (k, 4);* (G(m,

19、k)* cos ( BUS( m 4)(m- 3, m- 3)=BUS ( m 3)A 2 * G ( m m)sum ;BUS( n, 4)endelse(m- 3, n- 1) + B ( m, n) *endsin=BUS (m3) * BUS ( n ,3) *(BUS( m, 4) - BUS ( n, 4);(G( m, n) *cos ( BUS( m, 4)endfor m4: 9for4: 9ifm = nsum0;sumfor kBUS( k, 4)(m, k) *end=sum + BUS ( m, 3) * BUS ( k, 3)(BUS( m, 4) - BUS (k

20、, 4);* (G(m, k)* sin ( BUS( m 4)elseBUS( n, 4)cos(m- 3, m- 3)=BUS ( m 3)a 2 * B ( m, m)=-BUS ( m, 3) * BUS ( n, 3)sum ;(m- 3, n- 3)-B ( m, n) * cos ( BUS( m, 4) - BUS ( n, 4);end* (G( m n)* sin ( BUS( m, 4)endend%雅可比矩陣JACOB = H N ; M L ;%修正量delta_x= - JACOB delta_y%修正電壓幅值和相位form = 4: 9endform = 2: 9

21、BUS(m 4) = BUS ( m, 4) + delta_x(m- 1);end%計(jì)算精度tol = abs ( max ( delta_y );LoopCou nt= LoopCou nt + 1;endBUS (:,4)= BUS (:,4)disp disp(POWERFLOW *180 / pi ;,date );(BUS);Y_pf33=zeros ( 3);Y_pf39=zeros ( 3, 9);Y_pf99=Ypofl ;Y_pf33(1,1)= 1/( 1j *GEN( 1, 3);Y_pf33(2,2)= 1/( 1j *GEN( 2, 3);Y_pf33(3, 3)=

22、 1/( 1j *GEN( 3, 3);Y_pf39(1,1)= - 1/( 1j * GEN( 1 , 3);Y_pf39(2,2)= - 1/( 1j *GEN(2, 3);Y_pf39(3, 3)= - 1/( 1j *GEN(3, 3);Y_pf93=transpose( Y_pf39 );Y_pf99(1,1)= Y_pf99 (1,1) + 1/(1j *GEN( 1, 3);Y_pf99(2,2)= Y_pf99 ( 2, 2) + 1/(1j *GEN(2, 3);Y_pf99(3,3)= Y_pf99 ( 3, 3) + 1/(1j *GEN(3, 3);Y_pf99(5,5

23、)= Y_pf99 ( 5, 5) + BUS(5, 7)/ BUS( 5,Y_pf99(6,6)= Y_pf99 ( 6, 6) + BUS(6, 7)/ BUS( 6,Y_pf99(8,8)= Y_pf99 ( 8, 8) + BUS(8, 7)/ BUS( 8,Yex =Y_pf33 Y_pf39; Y_pf93 Y_pf99;Yrpf=Y_pf33- Y_pf39/ Y_pf99* Y_pf93 ;%故障中Y_df33=Y_pf33 ;Y_df38=Y_pf39 ;Y_df38(:,7)=;Y_df83=transpose( Y_df38 );Y_df88=Y_pf99 ;%故障分析簡(jiǎn)

24、化的Y陣 故障前計(jì)算增廣丫陣3)a3)A3)A1j1j1j* BUS( 5, 8)/* BUS( 6, 8)/* BUS( 8, 8)/BUS( 5, 3)ABUS( 6, 3)aBUS( 8, 3)a2；2；2;Y_df88(7,:)=；Y_df88(:, 7)=;Yrdf = Y_df33 - Y_df38 / Y_df88 * Y_df83%故障后Y_af33=Y_pf33Y_af39=Y_pf39Y_af93=Y_pf93Y_af99=Y_pf99Y_af99(5,5) = Y_af99(5, 5)-1/(BRANCH( 3, 3)+ 1j * BRANCH( 3, 4)-1j * B

25、RANCH( 3 , 5);Y_af99(7,7) = Y_af99(7, 7)-1/(BRANCH( 3, 3)+ 1j * BRANCH( 3, 4)-1j * BRANCH( 3 , 5);Y_af99(5, 7)= 0；Y_af99(7, 5)= 0；Yraf=Y_af33- Y_af39/Y_af99* Y_af93 ;disp ( Before fault ,date);disp (Yrpf );disp (Yrdf );disp (Yraf );%故障前發(fā)電機(jī)狀態(tài)tempAng = BUS (1:3, 4)* pi / 180 ;tempV1 = BUS ( 1 : 3, 3).

26、* exp (1 j .* tempAng );tempV2 = tempV1 + conj ( BUS( 1:3,5)+ 1j .* BUS( 1 : 3, 6)./ tempV1 ).*1j .* GEN( 1 : 3, 3);E_pf = sqrt ( real (tempV2 ).A 2 + imag (tempV2 ).A 2);Delta_pf = atan (imag (tempV2 )./ real (tempV2 ).* 180. / pi ;GEN (:, 5) = real(tempV2 .* conj ( Yrpf * tempV2 );disp(Genlnit ,da

27、tedisp(曰);disp(E_pf );disp(delta);disp(Delta_pf );%故障中發(fā)電機(jī)狀態(tài)% 故障切除時(shí)間%仿真結(jié)束時(shí)間%重合閘時(shí)間t_cut= 0.17 ;t_end= 4;t_ar = 0.42 ;odetol optio ns1e-10 ;=odeset ( RelTol , 1e-10 );Delta_i nit=Delta_pf.* pi ./180 ;x0_df=ones (3,1).*2. *pi .*60 ;Delta_i nitTdf_out,Xdf_out =ode45 ( t,x) Gen_Fun (t,x,GEN,Yrdf,E_pfDelta_df=Xdf_out (:,4:6).*180./ pi ;Omega_