本發(fā)明涉及一種數(shù)值仿真方法,具體涉及一種應(yīng)用于開關(guān)電路電磁暫態(tài)分析的數(shù)值仿真方法。
背景技術(shù):
電磁暫態(tài)分析是電力系統(tǒng)仿真的重要組成部分,是了解電力系統(tǒng)暫態(tài)復(fù)雜行為的必要工具。當(dāng)前常用的電磁暫態(tài)仿真軟件大多采用節(jié)點(diǎn)電壓法作為底層算法。節(jié)點(diǎn)電壓法是先采用數(shù)值積分方法對(duì)系統(tǒng)中動(dòng)態(tài)原件的特性方程差分化,得到等效的計(jì)算電導(dǎo)與歷史項(xiàng)電流源并聯(lián)形式的諾頓等效電路,進(jìn)而獲得整個(gè)電路系統(tǒng)的節(jié)點(diǎn)電導(dǎo)矩陣,求解得到系統(tǒng)中各節(jié)點(diǎn)電壓的瞬時(shí)值。使用節(jié)點(diǎn)分析法分析開關(guān)電路時(shí),存在的主要問題是在開關(guān)動(dòng)作時(shí)刻出現(xiàn)強(qiáng)烈的數(shù)值振蕩。雖然各個(gè)節(jié)點(diǎn)分析法軟件均采用一定的數(shù)值技術(shù)來消除振蕩,而這些消除數(shù)值振蕩的技術(shù)往往是一階的線性插值技術(shù),會(huì)帶來數(shù)值精度的下降。因此仍然要求采用非常小的步長來保證足夠的計(jì)算精度。
技術(shù)實(shí)現(xiàn)要素:
本發(fā)明的目的是提供一種應(yīng)用于開關(guān)電路電磁暫態(tài)分析的數(shù)值仿真方法,此方法在開關(guān)動(dòng)作時(shí)刻不會(huì)產(chǎn)生數(shù)值振蕩,具有數(shù)值穩(wěn)定性好,數(shù)值精度高(二階精度)的特點(diǎn)。
本發(fā)明的目的是采用下述技術(shù)方案實(shí)現(xiàn)的:
本發(fā)明提供一種應(yīng)用于開關(guān)電路電磁暫態(tài)分析的數(shù)值仿真方法,其改進(jìn)之處在于,所述方法包括下述步驟:
(1)列寫電路的系統(tǒng)狀態(tài)方程,并通過變量代換,將其變換一階形式的狀態(tài)方程;
(2)在開關(guān)動(dòng)作時(shí)刻,對(duì)電路的系統(tǒng)狀態(tài)方程進(jìn)行修改;
(3)利用直接積分法對(duì)一階狀態(tài)方程進(jìn)行迭代求解;
(4)在開關(guān)動(dòng)作時(shí)刻對(duì)電路節(jié)點(diǎn)電位進(jìn)行重新初始化操作;
(5)在開關(guān)狀態(tài)切換時(shí)刻使用極小步長后向歐拉運(yùn)算消除數(shù)值振蕩;
(6)開關(guān)動(dòng)作時(shí)刻的等價(jià)性證明。
進(jìn)一步地,所述步驟(1)中,電路的系統(tǒng)狀態(tài)方程表達(dá)式如下:
其中:
通過變量代換,將系統(tǒng)狀態(tài)方程①變換為一階狀態(tài)方程,如下:
其中:
k1、k2均表示系數(shù)矩陣;x為狀態(tài)變量;e是單位對(duì)角矩陣。
進(jìn)一步地,所述步驟(2)中,若連接在兩個(gè)節(jié)點(diǎn)之間的理想開關(guān)在開關(guān)動(dòng)作時(shí)刻發(fā)生動(dòng)作,則在開關(guān)動(dòng)作時(shí)刻td修改式①中的電阻系數(shù)矩陣kr;使用式④對(duì)電路的系統(tǒng)狀態(tài)方程進(jìn)行修改,式④表達(dá)式如下:
其中,turnon表示理想開關(guān)在td時(shí)刻開通;turnoff表示理想開關(guān)在td時(shí)刻斷開;除kr(i,i)=rs-1,kr(i,j)=-rs-1,kr(j,i)=-rs-1和kr(j,j)=rs-1以外,kr中的其他元素均為零,0v理想電壓源內(nèi)阻rs趨于零;kr,+是開關(guān)動(dòng)作前的電阻系數(shù)矩陣,kr,-是開關(guān)動(dòng)作后的電阻系數(shù)矩陣;kr是一個(gè)n階矩陣,kr(i,i)表示矩陣中第i行,第i列的元素;kr(i,j)表示矩陣中第i行,第j列的元素;kr(j,i)表示矩陣中第j行,第i列的元素;kr(j,j)表示矩陣中第j行,第j列的元素;
進(jìn)一步地,所述步驟(3)中,對(duì)于標(biāo)準(zhǔn)形式的一階動(dòng)態(tài)方程②,用直接積分法進(jìn)行迭代求解,直接積分法的迭代格式如下式⑤所示:
其中:δt是第n個(gè)時(shí)刻與n+1時(shí)刻之間的時(shí)間步長;β是一個(gè)可供選擇的系數(shù),當(dāng)β取值為0時(shí),稱為前差法,穩(wěn)定性為條件穩(wěn)定,數(shù)值精度為一階精度;當(dāng)β取值為0.5時(shí),稱為crank-nicolson法,穩(wěn)定性為無條件穩(wěn)定,數(shù)值精度為二階精度;當(dāng)β取值為1時(shí),稱為后差法,穩(wěn)定性為無條件穩(wěn)定,數(shù)值精度為一階精度;xn代表了第n步,即t=nδt時(shí)刻的系統(tǒng)狀態(tài)變量;xn+1代表了第n+1步,即t=(n+1)δt時(shí)刻的系統(tǒng)狀態(tài)變量;rn和rn+1分別代表了第n步,即t=nδt時(shí)刻和第n+1步,即t=(n+1)δt時(shí)刻的系統(tǒng)輸入向量。
進(jìn)一步地,所述步驟(4)包括:
在開關(guān)動(dòng)作時(shí)刻td,將電容等效為電壓源,將電感等效為電流源,根據(jù)td-時(shí)刻的電容電壓和電感電流計(jì)算td+時(shí)刻的電位分布;td-表示td時(shí)刻前的某一個(gè)時(shí)刻,且td-td-趨于零;td+表示td時(shí)刻后的某一個(gè)時(shí)刻,且td+-td趨于零;
若將電路中的電容、電感動(dòng)態(tài)元件均做等效處理后,電路中僅存在電阻支路和電流源支路,獲得靜態(tài)電路方程,求解式⑥中的
其中:
k=kr,d++krc,其中
對(duì)于某一個(gè)連接在i號(hào)節(jié)點(diǎn)與j號(hào)節(jié)點(diǎn)之間、電容為c的支路來講,系數(shù)矩陣krc中包含其等效成理想電壓源后,關(guān)于其內(nèi)阻rs趨于零的支路信息;矩陣krc中除krc(i,i)=rs-1,krc(i,j)=-rs-1,krc(j,i)=-rs-1,krc(j,j)=rs-1以外,其他元素均為零;krc是一個(gè)n階矩陣,krc(i,i)表示矩陣krc中第i行,第i列的元素;krc(i,j)表示矩陣krc中第i行,第j列的元素;krc(j,i)表示矩陣krc中第j行,第i列的元素;krc(j,j)表示矩陣krc中第j行,第j列的元素;
向量i中除包含原電路中所有電流源以及理想電壓源的諾頓等效電流源信息is,d+以外,還包含電感等效電流源的信息il,d-以及電容等效理想電壓源的信息isc,d-,如式⑧所示:
i=is,d++il,d-+isc,d-⑧
其中:
將式⑦和⑧代入⑥,并整理成式⑨的格式:
其中:
進(jìn)一步地,所述步驟(5)包括:
在開關(guān)動(dòng)作時(shí)刻td使用后向歐拉運(yùn)算:
如果在開關(guān)動(dòng)作時(shí)刻td,以一個(gè)極小的時(shí)間步長δt0對(duì)系統(tǒng)狀態(tài)方程①施加后向歐拉運(yùn)算,其中δt0=td+-td-,使用td-時(shí)刻的狀態(tài)量xd-,求得td+時(shí)刻的狀態(tài)量xd+,計(jì)算表達(dá)式如下:
其中,
考慮到δt0趨于0,式
其中:ψd+和ψd-分別表示開關(guān)動(dòng)作后各個(gè)節(jié)點(diǎn)的電壓時(shí)間積分量和開關(guān)動(dòng)作前各個(gè)節(jié)點(diǎn)的電壓時(shí)間積分量;
將式
如果式
進(jìn)一步地,所述步驟(6)包括下述步驟:
<1>電感相關(guān)項(xiàng)的等價(jià)性證明:
對(duì)于電感支路來講,存在式
或者寫成矩陣格式,如式
繼續(xù)擴(kuò)展為下述
由于每一個(gè)電感支路都有式
<2>電容相關(guān)項(xiàng)的等價(jià)性證明:
考慮電容支路的諾頓等效電流源,式⑨中的電容相關(guān)項(xiàng)展開為式
將式
將式
由于
開關(guān)動(dòng)作時(shí)刻使用極小步長后向歐拉運(yùn)算,等價(jià)于根據(jù)電路中的電容電壓和電感電流對(duì)節(jié)點(diǎn)電位重初始化的操作。
與最接近的現(xiàn)有技術(shù)相比,本發(fā)明提供的技術(shù)方案具有的優(yōu)異效果是:
由于電感元件上的電流方向始終為正,當(dāng)電流源的輸出電流過零時(shí),電感電流的變化率將發(fā)生由負(fù)到正的一次跳變,從而導(dǎo)致電感兩端的電壓
而本發(fā)明所述的狀態(tài)方程法,由于在開關(guān)動(dòng)作時(shí)刻使用了極小步長的后差運(yùn)算來對(duì)節(jié)點(diǎn)電位進(jìn)行重初始化操作,這樣在線性區(qū)段再使用梯形法時(shí),總體上仍然可以在保持二階數(shù)值精度。因此,本發(fā)明所述方法與基于梯形積分法建立方程的節(jié)點(diǎn)分析法相比,具有數(shù)值精度高的優(yōu)點(diǎn)。
附圖說明
圖1是本發(fā)明提供的開關(guān)動(dòng)作時(shí)刻電容的等效電路圖;
圖2是本發(fā)明提供的開關(guān)動(dòng)作時(shí)刻電感的等效電路圖;
圖3是本發(fā)明提供的電源頻率50hz的橋式電流源整流電路圖;
圖4是本發(fā)明提供的pscad軟件的仿真結(jié)果與本發(fā)明所述方法仿真結(jié)果的對(duì)比示意圖。
具體實(shí)施方式
下面結(jié)合附圖對(duì)本發(fā)明的具體實(shí)施方式作進(jìn)一步的詳細(xì)說明。
以下描述和附圖充分地示出本發(fā)明的具體實(shí)施方案,以使本領(lǐng)域的技術(shù)人員能夠?qū)嵺`它們。其他實(shí)施方案可以包括結(jié)構(gòu)的、邏輯的、電氣的、過程的以及其他的改變。實(shí)施例僅代表可能的變化。除非明確要求,否則單獨(dú)的組件和功能是可選的,并且操作的順序可以變化。一些實(shí)施方案的部分和特征可以被包括在或替換其他實(shí)施方案的部分和特征。本發(fā)明的實(shí)施方案的范圍包括權(quán)利要求書的整個(gè)范圍,以及權(quán)利要求書的所有可獲得的等同物。在本文中,本發(fā)明的這些實(shí)施方案可以被單獨(dú)地或總地用術(shù)語“發(fā)明”來表示,這僅僅是為了方便,并且如果事實(shí)上公開了超過一個(gè)的發(fā)明,不是要自動(dòng)地限制該應(yīng)用的范圍為任何單個(gè)發(fā)明或發(fā)明構(gòu)思。
本發(fā)明涉及一種應(yīng)用于開關(guān)電路電磁暫態(tài)分析的數(shù)值仿真方法,包括下述步驟:
(1)列寫電路的系統(tǒng)狀態(tài)方程,并通過變量代換,將其變換一階形式的狀態(tài)方程;
一般線性電路(由電阻、電容、電感以及電源等基本元件組成的,具有n+1個(gè)節(jié)點(diǎn)的電路,其中0號(hào)節(jié)點(diǎn)定義為參考地電位節(jié)點(diǎn))的系統(tǒng)狀態(tài)方程表達(dá)式如下:
其中:
通過變量代換,將系統(tǒng)狀態(tài)方程①變換為一階狀態(tài)方程,如下:
其中:
k1、k2均表示系數(shù)矩陣;x為狀態(tài)變量;e是單位對(duì)角矩陣。
(2)在開關(guān)動(dòng)作時(shí)刻,對(duì)電路的系統(tǒng)狀態(tài)方程進(jìn)行修改;
當(dāng)電路網(wǎng)絡(luò)中的理想開關(guān)元件發(fā)生動(dòng)作時(shí)(導(dǎo)通或關(guān)斷),在開關(guān)動(dòng)作時(shí)刻快速修改狀態(tài)方程系數(shù)矩陣的方法:
若連接在兩個(gè)節(jié)點(diǎn)之間的理想開關(guān)在開關(guān)動(dòng)作時(shí)刻發(fā)生動(dòng)作,則在開關(guān)動(dòng)作時(shí)刻t0修改式①中的電阻系數(shù)矩陣kr;使用式④對(duì)電路的系統(tǒng)狀態(tài)方程進(jìn)行修改,式④表達(dá)式如下:
其中,turnon表示理想開關(guān)在td時(shí)刻開通;turnoff表示理想開關(guān)在td時(shí)刻斷開;除kr(i,i)=rs-1,kr(i,j)=-rs-1,kr(j,i)=-rs-1和kr(j,j)=rs-1以外,kr中的其他元素均為零,0v理想電壓源內(nèi)阻rs趨于零;kr,+是開關(guān)動(dòng)作前的電阻系數(shù)矩陣,kr,-是開關(guān)動(dòng)作后的電阻系數(shù)矩陣;kr是一個(gè)n階矩陣,kr(i,i)表示矩陣中第i行,第i列的元素;kr(i,j)表示矩陣中第i行,第j列的元素;kr(j,i)表示矩陣中第j行,第i列的元素;kr(j,j)表示矩陣中第j行,第j列的元素;
(3)利用直接積分法對(duì)一階狀態(tài)方程進(jìn)行迭代求解;
對(duì)于標(biāo)準(zhǔn)形式的一階動(dòng)態(tài)方程②,用直接積分法進(jìn)行迭代求解,直接積分法的迭代格式如下式⑤所示:
其中:δt是第n個(gè)時(shí)刻與n+1時(shí)刻之間的時(shí)間步長;xn代表了第n步,即t=nδt時(shí)刻的系統(tǒng)狀態(tài)變量;xn+1代表了第n+1步,即t=(n+1)δt時(shí)刻的系統(tǒng)狀態(tài)變量;rn和rn+1分別代表了第n步,即t=nδt時(shí)刻和第n+1步,即t=(n+1)δt時(shí)刻的系統(tǒng)輸入向量;β是一個(gè)可供選擇的系數(shù),當(dāng)β取不同的值時(shí),相關(guān)的名稱如下表1所述:
表1β的選擇與數(shù)值穩(wěn)定性
其中,梯形法具有二階的數(shù)值精度,同時(shí)是無條件穩(wěn)定的,所以在線性問題分析中被廣泛采用。但在開關(guān)電路分析中,梯形法在開關(guān)動(dòng)作時(shí)刻會(huì)產(chǎn)生強(qiáng)烈的數(shù)值振蕩。
(4)在開關(guān)動(dòng)作時(shí)刻對(duì)電路節(jié)點(diǎn)電位進(jìn)行重新初始化操作;
當(dāng)開關(guān)動(dòng)作引起電路拓?fù)湓趖d時(shí)刻發(fā)生改變時(shí),其中的節(jié)點(diǎn)電位有可能在td時(shí)刻前后發(fā)生跳變(如圖3中,在二極管開通或關(guān)斷時(shí)刻,3號(hào)節(jié)點(diǎn)的電位會(huì)發(fā)生跳變),使用梯形法計(jì)算時(shí),會(huì)引起相關(guān)節(jié)點(diǎn)電位的數(shù)值振蕩。這就需要在開關(guān)時(shí)刻對(duì)電路各節(jié)點(diǎn)電位進(jìn)行重新初始化操作。
根據(jù)電路理論,在開關(guān)動(dòng)作前后,電容支路有保持兩端電位差不變的特性,電感支路有保持支路電流不變的特性。因此可以在td時(shí)刻,將電容等效為電壓源(如圖1),將電感等效為電流源(如圖2),根據(jù)td-時(shí)刻的電容電壓和電感電流計(jì)算td+時(shí)刻的電位分布。在開關(guān)動(dòng)作時(shí)刻td,將電容等效為電壓源,將電感等效為電流源,根據(jù)td-時(shí)刻的電位分布;td-表示td時(shí)刻前的某一個(gè)時(shí)刻,且td-td-趨于零;td+表示td時(shí)刻后的某一個(gè)時(shí)刻,且td+-td趨于零;
若將電路中的電容、電感等動(dòng)態(tài)元件都做了等效處理后,電路中僅存在電阻支路和電流源支路,獲得靜態(tài)電路方程。求解式⑥中的
其中:
(5)在開關(guān)狀態(tài)切換時(shí)刻使用極小步長后向歐拉運(yùn)算消除數(shù)值振蕩;
系數(shù)矩陣k中含電阻支路、電壓源內(nèi)阻支路的信息kr,d+和電容等效電壓源的內(nèi)阻信息krc;
對(duì)于某一個(gè)連接在i號(hào)節(jié)點(diǎn)與j號(hào)節(jié)點(diǎn)之間、電容為c的支路來講,系數(shù)矩陣krc中包含其等效成理想電壓源后,關(guān)于其內(nèi)阻rs趨于零的支路信息;矩陣krc中除krc(i,i)=rs-1,krc(i,j)=-rs-1,krc(j,i)=-rs-1,krc(j,j)=rs-1以外,其他元素均為零;krc是一個(gè)n階矩陣,krc(i,i)表示矩陣krc中第i行,第i列的元素;krc(i,j)表示矩陣krc中第i行,第j列的元素;krc(j,i)表示矩陣krc中第j行,第i列的元素;krc(j,j)表示矩陣krc中第j行,第j列的元素;
向量i中除包含原電路中所有電流源以及理想電壓源的諾頓等效電流源信息is,d+以外,還包含電感等效電流源的信息il,d-以及電容等效理想電壓源的信息isc,d-,如式⑧所示:
i=is,d++il,d-+isc,d-⑧
其中:
將式⑦和⑧代入⑥,并整理成式⑨的格式:
其中:
在開關(guān)動(dòng)作時(shí)刻td使用后向歐拉運(yùn)算:
如果在開關(guān)動(dòng)作時(shí)刻td,以一個(gè)極小的時(shí)間步長δt0對(duì)系統(tǒng)狀態(tài)方程①施加后向歐拉運(yùn)算,其中δt0=td+-td-,使用td-時(shí)刻的狀態(tài)量xd-,求得td+時(shí)刻的狀態(tài)量xd+,計(jì)算表達(dá)式如下:
其中,
將式⑩進(jìn)一步展開為式
考慮到δt0趨于0,式
其中:ψd+和ψd+分別表示開關(guān)動(dòng)作后各個(gè)節(jié)點(diǎn)的電壓時(shí)間積分量和開關(guān)動(dòng)作前各個(gè)節(jié)點(diǎn)的電壓時(shí)間積分量;
將式
如果式
(6)開關(guān)動(dòng)作時(shí)刻的等價(jià)性證明:
<1>電感相關(guān)項(xiàng)的等價(jià)性證明:
對(duì)于電感支路來講,存在式
或者寫成矩陣格式,如式
繼續(xù)擴(kuò)展為下述
由于每一個(gè)電感支路都有式
<2>電容相關(guān)項(xiàng)的等價(jià)性證明:
考慮電容支路的諾頓等效電流源,式⑨中的電容相關(guān)項(xiàng)展開為式
將式
將式
由于
這就證明了在開關(guān)動(dòng)作時(shí)刻:使用極小步長后差法,等價(jià)于根據(jù)電路中的電容電壓和電感電流對(duì)節(jié)點(diǎn)電位重初始化的操作。
實(shí)施例
圖3示出了一個(gè)橋式電流源整流電路,表2是其電路拓?fù)浯鎯?chǔ)在計(jì)算機(jī)中的形式。使用專利方法進(jìn)行仿真分析,基本步長設(shè)定為2ms。
表2算例的支路信息
首先對(duì)表3中的1~5號(hào)支路執(zhí)行電路網(wǎng)絡(luò)狀態(tài)方程的列寫方法的代碼,得到系數(shù)矩陣如式
從而得到電路的系統(tǒng)方程
其中:
定義零時(shí)刻的狀態(tài)量為x0=(0,0,0,0,900,0,0,0,0,0)t。
5.1前半個(gè)周期
在0+時(shí)刻,由于理想二極管的單向?qū)ㄌ匦裕琩2和d3將開通,同時(shí)d1和d4保持關(guān)斷。參照式所述的規(guī)則,將初始狀態(tài)的電阻系數(shù)矩陣kr修改為kr1,如式
其中:
5.1.1單步后差計(jì)算
由于d2和d3在零時(shí)刻開通,為了避免出現(xiàn)數(shù)值振蕩,這是需要施加一次后向差分法,來獲取0+時(shí)刻的狀態(tài)量。這個(gè)后向差分運(yùn)算的時(shí)間步長取δt0=10-7s,并將δt0和β=1代入式,得到t=10-7s時(shí)刻的系統(tǒng)狀態(tài)量
其中
5.1.2梯形迭代計(jì)算
隨后調(diào)整時(shí)間步長為δt=0.002-10-7s,根據(jù)t=10-7s時(shí)刻的系統(tǒng)狀態(tài)量
5.2后半個(gè)周期
在t=10ms時(shí)刻,d2和d3關(guān)斷,d1和d4導(dǎo)通,參照式所述的規(guī)則,將前半個(gè)周期的系數(shù)矩陣k21修改為后半個(gè)周期的k22,如式
其中
與前半個(gè)周期一樣,在電路拓?fù)淝袚Q的t=10ms時(shí)刻,取δt0=10-7s,得到t=0.01+10-7s時(shí)刻的系統(tǒng)狀態(tài)量
結(jié)果分析:
仿真得到的電感電壓如表3和圖4所示。同時(shí)給出了pscad軟件的仿真結(jié)果進(jìn)行對(duì)比,其中的精確解是使用pscad軟件在1μs步長下獲得的。
表3電感電壓
以上實(shí)施例僅用以說明本發(fā)明的技術(shù)方案而非對(duì)其限制,盡管參照上述實(shí)施例對(duì)本發(fā)明進(jìn)行了詳細(xì)的說明,所屬領(lǐng)域的普通技術(shù)人員依然可以對(duì)本發(fā)明的具體實(shí)施方式進(jìn)行修改或者等同替換,這些未脫離本發(fā)明精神和范圍的任何修改或者等同替換,均在申請(qǐng)待批的本發(fā)明的權(quán)利要求保護(hù)范圍之內(nèi)。