簡介一份基于MATLAB的6節(jié)點天然氣潮流計算教學(xué)程序面向能源動力及相關(guān)專業(yè)初學(xué)者用于理解天然氣網(wǎng)絡(luò)壓力、流量與儲存量的求解方法。程序以6節(jié)點簡化模型為對象完整涉及狀態(tài)方程、能量守恒、管道阻力與壓降計算等核心環(huán)節(jié)壓降估算可參考達西-韋斯巴赫公式或海曾-威廉公式并采用牛頓法或高斯-塞德爾迭代求解相互關(guān)聯(lián)的非線性方程組幫助讀者把理論公式落地為可運行代碼。壓縮包內(nèi)含2個M腳本文件體積僅約1KB代碼注釋清晰、結(jié)構(gòu)緊湊包含網(wǎng)絡(luò)數(shù)據(jù)輸入、潮流計算主體、結(jié)果輸出與可視化流程。已有499人學(xué)習(xí)下載適合作為天然氣仿真分析入門的參考模板也可在此基礎(chǔ)上進一步擴展至更復(fù)雜的管網(wǎng)模型。 做天然氣潮流計算的朋友應(yīng)該都有這種體會書上的公式看著不難真正能把一個算例跑通、跑對需要跨過的坑比想象中多得多。最近幫師弟調(diào)試課程設(shè)計正好用的是這套6節(jié)點天然氣潮流計算MATLAB程序趁著記憶還熱乎把這個算例從建模思路、數(shù)據(jù)準(zhǔn)備到程序?qū)崿F(xiàn)和調(diào)試經(jīng)驗完整梳理一遍。內(nèi)容不繞彎子直接把能復(fù)現(xiàn)的算例數(shù)據(jù)、代碼框架和踩坑記錄放在這里給同樣在做天然氣管網(wǎng)仿真、綜合能源系統(tǒng)潮流計算或者剛接觸MATLAB編程的同學(xué)參考。1. 天然氣潮流計算到底在算什么1.1 先把問題跟電力潮流放一起看很多人第一次看到“潮流計算”四個字會下意識把它理解成電力系統(tǒng)里的牛拉法、PQ分解法那一套。這個直覺沒錯天然氣管網(wǎng)潮流和電力系統(tǒng)潮流本質(zhì)上是一類問題給定網(wǎng)絡(luò)拓?fù)洹⑦吔鐥l件和負(fù)荷求解全網(wǎng)的狀態(tài)變量分布。電力潮流求的是節(jié)點電壓幅值和相角天然氣潮流求的是節(jié)點壓力和管段流量。兩者的不同點也很明顯。電力系統(tǒng)里有功潮流跟相角差近似線性無功跟電壓差強相關(guān)所以可以分PQ分解天然氣管網(wǎng)里管段流量跟壓降之間的關(guān)系是非線性的平方關(guān)系而且沒有“相角”這類緩變量方程組的非線性程度更高對初值更敏感。換句話說電力潮流不收斂的時候通常調(diào)調(diào)初值就好天然氣潮流不收斂的時候你可能要回頭檢查管段方程是不是寫錯了。但有趣的是求解框架完全可以復(fù)用列節(jié)點守恒方程加元件特性方程組成非線性方程組后用牛頓-拉夫遜法迭代求解。所以如果你已經(jīng)寫過電力潮流程序天然氣潮流對你來說就是個“換了元件方程”的版本。1.2 管網(wǎng)模型里兩個核心方程天然氣管網(wǎng)模型的本質(zhì)由兩類方程構(gòu)成。第一類是節(jié)點流量平衡方程。對任意一個節(jié)點流入該節(jié)點的流量之和減去流出該節(jié)點的流量之和必須等于該節(jié)點的負(fù)荷用氣量如果有注入源比如氣源節(jié)點則注入量減去全部流出量等于零。數(shù)學(xué)上寫成[ \sum_{j \in N(i)} q_{ij} Q_{load,i} ]這個方程的意義跟電力潮流的KCL一模一樣——物質(zhì)守恒。在程序里它是構(gòu)造殘差向量的基礎(chǔ)。第二類是管段壓降方程。天然氣在管道里流動時由于摩擦阻力產(chǎn)生壓降。對一段連接節(jié)點 (i) 和 (j) 的管道簡化的壓降-流量關(guān)系可以寫成[ p_i^2 - p_j^2 R_{ij} \cdot q_{ij}^2 ]這里 (R_{ij}) 是管段的等效阻力系數(shù)跟管長、管內(nèi)徑、天然氣物性相對密度、壓縮因子、溫度有關(guān)實際工程里常用Weymouth公式或Darcy公式計算。這個方程最大的特點是壓力出現(xiàn)在平方項里而流量也出現(xiàn)在平方項里典型的雙向非線性。正因為這個“平方對平方”的關(guān)系很多教材把狀態(tài)變量直接取成壓力的平方 (P2 p^2)方程形式會簡潔不少迭代時也更好處理。提示潮流計算中壓力一定要用絕壓不能用表壓。新手經(jīng)常在這一點上翻車算出來的結(jié)果偏差大到離譜。1.3 為什么選6節(jié)點作為入門算例我先解釋一下為什么拿6節(jié)點而不是3節(jié)點或者IEEE標(biāo)準(zhǔn)節(jié)點來說事。3節(jié)點算例只能做單一鏈?zhǔn)浇Y(jié)構(gòu)氣源-中游-末端體現(xiàn)不出環(huán)網(wǎng)的流量分配問題而IEEE那種動輒幾十節(jié)點的算例對初學(xué)者來說數(shù)據(jù)準(zhǔn)備和結(jié)果校核都是負(fù)擔(dān)。6節(jié)點是一個非常適中的規(guī)模結(jié)構(gòu)上可以同時包含串聯(lián)干線、分支支路和環(huán)網(wǎng)能完整反映天然氣管網(wǎng)的主要特征結(jié)果規(guī)模又足夠小可以手工校核任意一條管段的壓降計算是否合理。很多教材和論文里的演示算例都采用這個規(guī)模參考資料也容易找。我自己帶學(xué)生的感受是6節(jié)點算例跑通之后理解“平衡節(jié)點”“定流量節(jié)點”“迭代初值”“殘差收斂”這些概念就都具體化了再往上擴到十幾節(jié)點或幾十節(jié)點只是數(shù)據(jù)量的增加不再是思路層面的障礙。2. 6節(jié)點算例的搭建與數(shù)據(jù)準(zhǔn)備2.1 管網(wǎng)拓?fù)浣Y(jié)構(gòu)設(shè)計我設(shè)計的這個6節(jié)點算例拓?fù)渖狭η蟾采w三種典型結(jié)構(gòu)鏈?zhǔn)礁删€、分支支路、以及一個閉合環(huán)網(wǎng)。具體結(jié)構(gòu)是這樣的節(jié)點1是氣源節(jié)點相當(dāng)于電力系統(tǒng)中的平衡節(jié)點壓力給定節(jié)點1到節(jié)點2到節(jié)點3是主輸氣干線管徑較大節(jié)點2分出支路到節(jié)點4節(jié)點4再向下到節(jié)點5這是末端支路節(jié)點3和節(jié)點4之間通過節(jié)點6形成一個閉合環(huán)路即3-6-4-3。為什么要布置一個環(huán)網(wǎng)因為環(huán)網(wǎng)是天然氣潮流計算中最容易出問題的地方流量分配不是直觀的“上游到下游”而是需要解方程組才能確定流量方向也可能跟初始猜測相反。如果只用樹狀管網(wǎng)用序貫法手算也能算體現(xiàn)不出牛頓-拉夫遜法的價值。環(huán)網(wǎng)結(jié)構(gòu)一旦出現(xiàn)就必須用到雅可比矩陣迭代求解。2.2 節(jié)點負(fù)荷與管段參數(shù)表為了讓這個算例能直接復(fù)現(xiàn)我把數(shù)據(jù)全部列出來。單位系統(tǒng)統(tǒng)一使用壓力單位MPa流量單位百萬立方米每天MMscmd管道阻力系數(shù) (R) 的單位為 (\text{MPa}^2 / (\text{MMscmd})^2)。節(jié)點負(fù)荷數(shù)據(jù)節(jié)點編號節(jié)點類型負(fù)荷MMscmd說明1平衡節(jié)點氣源0壓力固定為5.0 MPa2負(fù)荷節(jié)點0.8民用/工業(yè)用氣3負(fù)荷節(jié)點1.2民用/工業(yè)用氣4負(fù)荷節(jié)點0.5支路用戶5負(fù)荷節(jié)點0.6末端用戶6負(fù)荷節(jié)點0.7環(huán)網(wǎng)用戶總負(fù)荷是3.8 MMscmd全部由節(jié)點1的氣源供應(yīng)。管段參數(shù)和等效阻力系數(shù)管道編號起點終點管長km管徑mm等效阻力系數(shù)R11282000.18223102000.2232461501.1044581501.4553651500.9266461501.10我給的R值是折算后的等效值實際工程中需要用Weymouth公式從管長、管徑、天然氣物性計算。這里直接給折算值的目的是先把潮流計算的框架打通R的精確計算放到后面再完善。如果讀者打算換成實際公式只需要寫一個從管段參數(shù)計算R的函數(shù)然后替換掉輸入數(shù)據(jù)中的R列即可。2.3 參數(shù)單位換算這個坑預(yù)算是非常關(guān)鍵的一步也是最容易出大問題的地方。我第一次幫別人調(diào)這個程序的時候發(fā)現(xiàn)計算結(jié)果跟手算完全對不上查了一個多小時最后發(fā)現(xiàn)是壓力單位混用了迭代內(nèi)部用的Pa但管道阻力系數(shù)是用MPa代入公式算的兩者差了12個量級程序不跑飛才怪。建議從一開始就確定統(tǒng)一的單位體系并且在整個程序內(nèi)保持一致。我習(xí)慣的做法是外部數(shù)據(jù)輸入用工程單位MPa、km、mm進入程序后在讀取數(shù)據(jù)時立刻轉(zhuǎn)換成迭代內(nèi)部單位Pa、m算完再轉(zhuǎn)回工程單位輸出。這樣雖然多兩步轉(zhuǎn)換代碼但能避免最常見的單位災(zāi)難。3. MATLAB程序核心實現(xiàn)3.1 輸入數(shù)據(jù)的代碼組織在MATLAB里我習(xí)慣用向量存節(jié)點數(shù)據(jù)、用矩陣存管段數(shù)據(jù)。不需要復(fù)雜結(jié)構(gòu)體數(shù)據(jù)量小的時候直接數(shù)組操作更快、更直觀。對應(yīng)上面的算例輸入模塊長這樣% 節(jié)點負(fù)荷單位MMscmd nodeLoad [0; 0.8; 1.2; 0.5; 0.6; 0.7]; % 管段數(shù)據(jù)每行 [起點, 終點, 等效阻力系數(shù)R] pipe [ 1 2 0.18 2 3 0.22 2 4 1.10 4 5 1.45 3 6 0.92 6 4 1.10 ]; % 氣源節(jié)點和壓力 sourceNode 1; sourcePressure 5.0; % MPa % 迭代參數(shù) tol 1e-8; % 收斂容差 maxIter 100; % 最大迭代次數(shù)這樣組織數(shù)據(jù)的優(yōu)點是管段數(shù)量和拓?fù)湫薷亩贾粍舆@兩個數(shù)組程序其他部分完全不用改。等到以后擴展到實際管網(wǎng)把Excel或數(shù)據(jù)庫的數(shù)據(jù)讀進來替換這里的賦值語句就行。3.2 牛頓-拉夫遜法迭代框架核心求解思路是把狀態(tài)變量設(shè)為各節(jié)點壓力平方 (P2 p^2)。對每個節(jié)點氣源節(jié)點除外構(gòu)造殘差[ F_i \sum_{j \in N(i)} q_{ij} - Q_{load,i} ]其中管段流量由壓降方程反解[ q_{ij} \text{sign}(P2_i - P2_j) \cdot \sqrt{\frac{|P2_i - P2_j|}{R_{ij}}} ]加 (\text{sign}) 的原因是環(huán)網(wǎng)中流量方向未知必須根據(jù)兩端壓力大小自動判斷方向。雅可比矩陣的元素是 (F_i) 對 (P2_i) 的偏導(dǎo)數(shù)。對一段管道的流量項能直接推到解析表達式[ \frac{\partial q_{ij}}{\partial (P2_i)} \frac{1}{2 \sqrt{|P2_i - P2_j| \cdot R_{ij}}} ]在程序里組裝雅可比矩陣時對每根管道這部分的貢獻同時作用于四個位置(J(i,i))、(J(j,j)) 加正號(J(i,j))、(J(j,i)) 加負(fù)號。核心迭代代碼的骨架如下P2 sourcePressure^2 * ones(6,1); % 初值 for iter 1:maxIter F zeros(6,1); J zeros(6,6); for k 1:size(pipe,1) i pipe(k,1); j pipe(k,2); R pipe(k,3); dp2 P2(i) - P2(j); q sign(dp2) * sqrt(abs(dp2) / R); F(i) F(i) q; F(j) F(j) - q; dq 1 / (2 * sqrt(abs(dp2) * R 1e-12)); J(i,i) J(i,i) dq; J(j,j) J(j,j) dq; J(i,j) J(i,j) - dq; J(j,i) J(j,i) - dq; end F F - nodeLoad; % 節(jié)點流量平衡殘差 % 處理平衡節(jié)點將氣源節(jié)點方程替換為定壓約束 F(sourceNode) P2(sourceNode) - sourcePressure^2; J(sourceNode, :) 0; J(sourceNode, sourceNode) 1; % 迭代修正 delta -J \ F; P2 P2 delta; if max(abs(delta)) tol break; end end % 換算回壓力 pressure sqrt(P2);這里有個細(xì)節(jié)值得注意在求導(dǎo)表達式中我加了 (1e-12) 作為數(shù)值保護。當(dāng)某管段兩端壓力平方差接近零時導(dǎo)數(shù)會趨于無窮大加一個極小量可以防止數(shù)值溢出。雖然6節(jié)點算例里不一定觸發(fā)這個問題但這個習(xí)慣在擴展到大型管網(wǎng)時能救命。3.3 平衡節(jié)點的處理邏輯平衡節(jié)點氣源的處理是牛頓-拉夫遜法在這個問題里的關(guān)鍵細(xì)節(jié)。節(jié)點1的壓力是給定的所以它的節(jié)點方程不再是流量平衡方程而是定壓約束方程[ P2_1 p_{source}^2 ]在矩陣?yán)镒銎饋砗苤苯影言摴?jié)點對應(yīng)的殘差F行改成 (P2_1 - p_{source}^2)把雅可比矩陣對應(yīng)行清零只在對角元置1。這樣既保留了矩陣結(jié)構(gòu)的完整性又實現(xiàn)了壓力約束。有個常見錯誤是直接刪掉平衡節(jié)點所在行和列這在網(wǎng)絡(luò)規(guī)模大、稀疏矩陣處理時反而引入麻煩而且在輸出節(jié)點壓力時還得再映射回去代碼繞來繞去容易出bug。固定行替換的方法在6節(jié)點規(guī)模下迭代次數(shù)幾乎沒差別代碼還更簡潔。3.4 收斂判據(jù)怎么選我用的收斂判據(jù)是壓力平方修正量的最大絕對值小于 (10^{-8})。這個閾值在壓力以MPa為單位時對應(yīng)壓力修正大約在 (10^{-9}) MPa量級精度完全夠用。也有一些實現(xiàn)用殘差范數(shù)做判據(jù)比如 (||F||_\infty 10^{-6})。兩種判據(jù)在已經(jīng)收斂的迭代序列上差別不大但用狀態(tài)修正量做判據(jù)有個好處不會因為某個節(jié)點負(fù)荷特別小導(dǎo)致殘差天然很小而產(chǎn)生“假收斂”。實際調(diào)試時建議兩種判據(jù)都打印出來看一個判斷解的穩(wěn)定性一個判斷方程滿足程度。4. 運行結(jié)果怎么看4.1 參考輸出示例按上面的數(shù)據(jù)和代碼跑完壓力結(jié)果大致如下節(jié)點編號壓力MPa15.00024.8334.6744.7154.5564.62整體規(guī)律符合物理直覺離氣源越遠、負(fù)荷越重的節(jié)點壓力越低環(huán)網(wǎng)中節(jié)點3、4、6的壓力比較接近這正是環(huán)網(wǎng)平衡流量后自然形成的壓力分布。末端節(jié)點5壓力最低因為它在支路末端且管徑偏小。如果算出來的結(jié)果里出現(xiàn)某個節(jié)點壓力比氣源還高那一定有問題——可能是流量方向?qū)懛戳艘部赡苁茄趴杀染仃嚪栧e誤。4.2 自己動手校核的三種方法程序跑通不等于跑對我建議每個結(jié)果都要做三個層面的校核。第一是全局校核。所有負(fù)荷相加是3.8 MMscmd氣源節(jié)點的供氣量必須也是3.8通過節(jié)點1的流量平衡方程反算。如果不相等說明某個節(jié)點方程寫錯了。第二是局部校核。隨便挑一根管道比如管段1-2用輸出結(jié)果計算 (p_1^2 - p_2^2)再計算 (R \cdot q^2)兩邊應(yīng)該嚴(yán)格相等。這是最直接的方程校驗?zāi)芸焖俣ㄎ怀绦蚶锬囊徊匠隽藛栴}。第三是靈敏度校核。把某個節(jié)點的負(fù)荷增大10%看這個節(jié)點以及上游節(jié)點的壓力是否下降。如果壓力反而上升了那程序肯定有邏輯錯誤。這個方法不需要額外工具就是一種很好的程序自檢習(xí)慣。4.3 從結(jié)果反推管網(wǎng)瓶頸潮流計算不只是“算個壓力”而已它的輸出可以直接用于管網(wǎng)規(guī)劃與運行分析。比如這個算例的結(jié)果里節(jié)點5的壓力是4.55 MPa如果設(shè)計要求末端壓力不低于4.6 MPa那這個管網(wǎng)就是不滿足要求的需要采取升壓或擴容措施。這種“結(jié)果驅(qū)動決策”的思路是潮流計算作為工具的真正價值在管網(wǎng)規(guī)劃階段預(yù)測不同負(fù)荷場景下的壓力分布在運行階段判斷管網(wǎng)是否有瓶頸、是否需要增壓在綜合能源系統(tǒng)研究中為電-氣耦合分析提供天然氣側(cè)的運行狀態(tài)。6節(jié)點算例雖然小但這個分析邏輯跟實際工程完全一致。5. 調(diào)試實錄那些年踩過的坑5.1 不收斂先查初值和符號牛頓法不收斂的原因按出現(xiàn)頻率排序大概是初值偏離太遠、雅可比矩陣符號寫反、單位混用。初值問題是最常見的。這個算例里氣源壓力5.0 MPa如果你把所有節(jié)點初始壓力設(shè)為0迭代第一步就會遇到很大的導(dǎo)數(shù)矩陣數(shù)值特性差很容易震蕩發(fā)散。我建議把初始壓力設(shè)為氣源壓力的0.8倍左右也就是初步認(rèn)為不同節(jié)點壓力差別不會太大。對絕大多數(shù)中壓、低壓管網(wǎng)這個初值策略都很穩(wěn)。符號問題則隱蔽得多。雅可比矩陣的符號錯誤會讓殘差序列呈現(xiàn)“來回蹦”的特征一會兒正、一會兒負(fù)中間還出現(xiàn)過接近收斂的情況然后突然又彈出很遠。遇到這種癥狀最有效的排查辦法是把迭代前幾步的雅可比矩陣打印出來手工檢查第一個節(jié)點的對角線元素的導(dǎo)數(shù)表達式是否正確。5.2 壓力變成負(fù)數(shù)或復(fù)數(shù)如果你在迭代過程中發(fā)現(xiàn)壓力平方出現(xiàn)負(fù)值說明某個節(jié)點的“壓力平方”被迭代修成了負(fù)數(shù)開方后就是復(fù)數(shù)MATLAB會輸出NaN或者復(fù)數(shù)結(jié)果。這種情況絕大多數(shù)是物理模型本身出了問題負(fù)荷過大、管徑過小、氣源壓力不足。在6節(jié)點算例里如果你把總負(fù)荷加大到10 MMscmd以上很可能觸發(fā)這個問題。它的含義是在這個負(fù)荷水平下現(xiàn)有管網(wǎng)結(jié)構(gòu)根本無法滿足供氣要求潮流計算在數(shù)學(xué)上無解。遇到這種情況不要試圖靠調(diào)初值繞過而應(yīng)該回頭檢查參數(shù)合理性。這也是潮流程序的一個隱形價值它可以作為管網(wǎng)可行性的判定工具。5.3 常見問題速查表現(xiàn)象可能原因處理辦法完全不收斂殘差猛增初值太差或雅可比矩陣符號錯誤調(diào)初值為氣源壓力的0.8倍打印雅可比矩陣逐項核對迭代在某個值附近振蕩流量方向誤判sign函數(shù)缺失檢查管段流量計算是否加了sign節(jié)點壓力為NaN或復(fù)數(shù)負(fù)荷過大或管徑過小無可行解減小負(fù)荷或增大管徑檢查R值結(jié)果跟手算對不上單位混用統(tǒng)一用MPa和MMscmd或迭代內(nèi)部統(tǒng)一用SI單位氣源節(jié)點流量不為總負(fù)荷平衡節(jié)點方程處理錯誤檢查該節(jié)點的定壓方程是否被錯誤替換掉某根管道流量為負(fù)但壓力差為正管段端點編號順序與實際方向不同流量方向以sign為準(zhǔn)負(fù)號是正常現(xiàn)象這個小表基本覆蓋了我自己調(diào)試這個程序時遇到的主要問題也幾乎覆蓋了學(xué)生在這類作業(yè)里問過我的所有問題。6. 從6節(jié)點往外擴展的方向6.1 換數(shù)據(jù)就能跑的拓展思路6節(jié)點程序跑通之后往實際管網(wǎng)擴展的路徑其實很清晰。把節(jié)點數(shù)、管道數(shù)從數(shù)組里替換成實際數(shù)據(jù)加上壓縮機模型、閥門模型再把稀疏矩陣技術(shù)用起來就是一個簡化版的天然氣管網(wǎng)仿真工具。我實際測試過把同樣的代碼框架擴展到20節(jié)點左右的環(huán)形管網(wǎng)計算時間仍然在毫秒級。真正需要改進的是矩陣求解部分當(dāng)節(jié)點數(shù)到幾百甚至上千時直接把雅可比矩陣存成稠密矩陣就不合適了需要改成稀疏存儲并配合稀疏線性方程求解。不過這些都屬于“工程優(yōu)化”底層方程和求解邏輯跟6節(jié)點算例沒有本質(zhì)區(qū)別。所以把這套基礎(chǔ)程序吃透等于給后續(xù)拓展打好了地基。6.2 綜合能源系統(tǒng)方向如果做的是電-氣耦合系統(tǒng)研究可以把這份6節(jié)點天然氣潮流程序跟一個簡單的電力系統(tǒng)潮流程序接力起來用天然氣網(wǎng)輸出的氣源供氣量去約束燃?xì)廨啓C的出力再用電網(wǎng)的負(fù)荷需求去影響天然氣的用氣量。兩步交替迭代就是一個最簡單的電-氣聯(lián)合潮流雛形。這個方向這幾年在綜合能源系統(tǒng)研究里很熱門但很多論文里的算例就是把這兩個程序的接口做了一下對接。基礎(chǔ)能力還是各自領(lǐng)域的潮流求解所以先練好6節(jié)點天然氣潮流對后面做系統(tǒng)級研究會很有幫助。6.3 一點實際操作的體會最后分享一個我每次帶新手都會強調(diào)的建議跑通這個程序之后別急著交差把負(fù)荷數(shù)據(jù)改一改再跑幾遍。比如把節(jié)點5的負(fù)荷從0.6改成1.5觀察節(jié)點壓力怎么變把節(jié)點4和節(jié)點6之間的連接斷開看環(huán)網(wǎng)變成支路后流量怎么重新分配。這些“折騰”對理解管網(wǎng)運行特性的幫助比單純把程序跑通要大得多。我最初帶師弟做這個題目時他調(diào)試到凌晨三點最后發(fā)現(xiàn)是壓力單位混用。折騰的過程雖然痛苦但那次之后他對單位換算、初值選擇這類細(xì)節(jié)有了刻骨銘心的記憶后面再調(diào)電網(wǎng)、燃?xì)庀到y(tǒng)的程序都順了很多。這種基本功靠看教程是學(xué)不來的必須親自踩一次坑。本文還有配套的精品資源點擊獲取