擾動法在風網解算中的應用研究*

趙自豪,李鵬慧,呂海建

(內蒙古科技大學 礦業與煤炭學院，內蒙古 包頭 014010)

0 引言

礦井分風系統解算的本質是解非線性方程組，解算方法可以分為回路分析法和節點分析法2類。回路分析法利用節點風量連續方程和回路風壓連續方程構建方程組，在處理進風井和出風井的開口節點時，需要引入虛擬回路的概念[1-2]，該方法多用牛頓迭代法進行求解，其難點在于獨立回路矩陣的求取、解算結果的初值依賴和如何提高解算效率上；節點分析法利用節點風量連續方程和節點壓降方程構建方程組，在處理進風井和出風井的開口節點時，往往引入固定等級節點概念[3-4]，該方法多用線性逼近法進行解算，其難點在于如何處理計算結果的擺動收斂問題。

回路分析法按照解算方程組的方法不同，分為牛頓法、斯考特-恒斯雷法、連續法(同倫法)等。牛頓法是目前最常用的風網解算方法，該方法雖然對初值有依賴，但理論上嚴格收斂，且有唯一解，其缺點是每次迭代都需要重新計算雅克比矩陣的逆矩陣，計算量偏大。文獻[5-8]在對雅克比矩陣的對稱性進行分析的基礎上，提出利用LDLT分解法求解雅克比矩陣的逆矩陣，并用計算機并行運算的方式進行實驗驗證，結果表明利用LDLT分解法相較于傳統解法可節約50%的工作量。斯考特-恒斯雷法是基于在牛頓法不嚴密假設基礎上的1種風網快速解法，該方法較牛頓法計算量小，但容易出現計算結果不收斂的情況。文獻[9-11]基于傳統的最小生成樹原理對斯考特-恒斯雷法進行改進，通過改進的BFS生成樹創建法和雙通路快速回路圈劃法，有效避免斯考特-恒斯雷法的初值依賴問題，加大收斂速度，避免計算結果發散的情況發生。擾動法是在工程和科學計算中常用的1種解析近似方法，該方法基于漸近級數理論，對擾動級數在工程允許的范圍內進行截斷，從而獲得1個近似的表達式，該表達式的解根據截斷項的多少，在一定精度內逼近原方程的解。利用該方法，無需執行循環運算，無需設定初值，計算時間可預估，在工程允許的誤差范圍內，簡單有效。Basha等[12]首先在供水網絡解算中引入擾動分析法，通過對含多級加壓水泵和出水點的復雜水網進行模擬計算，表明在工程允許的精度范圍內，該方法相對于傳統方法快捷高效；Bush[13]從數學的角度論證擾動法在解算多元二次方程組時的精度問題和收斂問題，從理論上證明該方法解算流體網絡的可行性；王樹剛等[14]將擾動法的原理應用于流體網絡的靈敏度分析問題，提出流體網絡中靈敏分支的特征值計算公式和靈敏分支的識別方法，并進行理論驗證。

本文在深入研究礦井風網解算傳統方法的優缺點及發展趨勢的基礎上，將擾動法引入風網解算領域；并針對風網解算的特點，演化出擾動法的矩陣表示法，針對3種常見的風機工況點，提出擾動法的擾動方程系數矩陣的構成原理和MATLAB擾動解算法的算法步驟；就擾動法中出現的數學公式定義域問題給出解決辦法；通過2個具體實例，驗證該方法的可靠性和高效性。

1 通風基本關系式的擾動法表示

擾動法是1種基于級數的近似算法，本文用到的2個級數展開如式(1)～(2)所示：

(1)

(2)

式中：z為擾動法展開式里的1個未知數。

1.1 通風阻力定律的擾動法表示

在不考慮方向性的情況下，通風阻力為正，則反映風量與通風阻力之間關系的通風阻力定律可改寫為式(3)[15]：

Q=αhx

(3)

式中：Q為風量，m3/s；h為通風阻力，Pa；α為風量和通風阻力的關系系數；x為常數，取0.5。

α可以用風阻表示，如式(4)所示：

(4)

式中：R為風阻，Pa。

對式(3)利用δ展開，可以獲得1個關于變量δ的多項式擾動級數。將h表示為式(5)：

h=h0+h1δ+h2δ2+h3δ3+h4δ4+O(δ5),δ=x-1

(5)

式中：δ為用擾動法展開式的變量；h0～h4均為展開式各項的系數；O(δ5)為展開式的余項。

則得式(6)：

hx=hhδ=heδlnh

(6)

將式(5)代入lnh，并取前5項，得式(7)：

(7)

利用式(1)～(2)和式(5)，則式(6)最終展開為δ擾動級數，得式(8)：

(8)

則通風阻力定律可表示為式(9)：

(9)

1.2 風機特性曲線的擾動展開表達式

對于任意采用機械通風的礦井，至少存在1個風機。風機的風量和風壓之間的關系可表示為式(10)：

(10)

式中：hp為風機的風壓,Pa；Qp為風機的風量,m3/s；p,q,r為各項系數；η為風機效率。

1)若p>0

通過數學變換，式(10)可以表示成式(11)：

Qp=ap(hp+cp)x-bp

(11)

2)若p<0

同理，得式(12)：

Qp=ap(-hp+cp)x-bp

(12)

3)若p=0

此時風機特性曲線為1條直線。

類似式(8)，風機性能函數方程可展開為式(13)：

(13)

2 風網中風流流動的基本定律

在1個分風網絡中存在如式(14)所示的關系：

NP=NJ+NL+NF-1

(14)

式中：NP為分支數；NJ為連接節點數；NL為獨立回路數；NF為固定等級節點數。

此處，固定等級節點是在水力網絡解算時的常用說法，指的是固定水頭壓力或固定流量的節點。在風網解算中，也就是進風井或出風井的開口節點。為表達簡潔僅討論NF=2的情形。

對于連接節點，滿足風量連續方程，如式(15)所示：

BAHx=0

(15)

式中：A=(α)NP×NP為對角線矩陣，對角線上按順序存儲每個分支的α值；B=(b)NJ×NP為連接節點的基本關聯矩陣；H=(h)NP×1為分支風壓向量[14-15]。

對于獨立回路，滿足回路風壓連續方程，如式(16)所示：

CH=0

(16)

式中：C=(c)NL×NP為僅與連接節點相聯分支的獨立回路矩陣。

在風網解算中，可以用虛擬回路的方式將固定等級節點轉變為連接節點。本文不假設虛擬回路，則有NF-1個通路壓降方程或固定等級節點風量連續方程。

常見的3種風機工況下的擾動方程矩陣表達式如下。

1)若風機處于固定風壓工況點Hp

根據固定等級節點的分布，選定NF-1個連接2個固定等級節點的獨立風路組合，得出通路壓降方程，如式(17)所示：

FH=Hp

(17)

式中：F=(f)(NF-1)×NP為通路分支關聯方程系數矩陣；Hp=(hp)(NF-1)×1為每個通路對應的總壓降。

2)若風機處于固定風量工況點Qp

選定NF-1個固定等級節點，寫出針對固定等級節點的風量連續方程，如式(18)所示：

BfAHx=Qp

(18)

式中：Bf=(bf)(NF-1)×NP為固定等級節點的基本關聯矩陣；Qp=(qp)(NF-1)×1為對應固定等級節點的出入風量。

3)若給定風機的性能函數

此時需要同時確定風機的風量Qp和風壓Hp，相當于增加NF-1個未知量。

①通路壓降方程如式(19)所示：

FH-hp=0

(19)

②固定等級節點風量連續方程如式(20)所示：

BfAHx-ap(hp+cp)x=bp

(20)

3 分風網絡擾動法的矩陣表示法及求解公式

3.1 固定風壓工況分風網絡求解

對式(15)～(17)聯立方程，得式(21)：

(21)

令M,N,P如式(22)所示：

(22)

得式(23)：

Mh+N(hx-h)=P

(23)

將式(5)，式(8)代入式(23)，可得固定風壓工況點時的通風網絡方程擾動展開式。

根據擾動理論，首先忽略擾動展開式的擾動項，得關于h0的線性方程組，如式(24)所示：

(24)

然后，將h0代入通風方程擾動展開式，且只保留δ擾動項，同時令P=0，得出關于h1的線性方程組，如式(25)所示：

(25)

類似的，可求得h2，h3，h4的解。h2，h3的解如式(26)～(27)所示：

(26)

(27)

式中：h0，h1，h2，h3均為NP×1的列向量，如h0=(hoi)NP×1代表分支1～NP的通風阻力關于h的主擾動項系數。h1～h3分別為1～3次擾動項系數。

從以上各擾動系數的表達式形式來看，除擾動主項外，其他擾動系數括號外的部分均相同，其結果可重復利用。只有括號內的算式部分有差異。但該部分執行的是向量元素之間的代數運算，相比矩陣運算，其復雜度非常小。因此，隨著要求計算精度增加，其增加的運算量不明顯。

3.2 固定風量工況分風網絡求解

聯立式(15)～(16)，式(18)，得式(28)：

(28)

記M，N，P如式(29)所示：

(29)

可寫出同式(23)形式完全相同的方程式，其求解過程與式(24)～(27)完全相同。

3.3 給定風機性能函數的分風網絡求解

假定礦山風機為抽出式通風，對于連接風機的固定等級節點，其風機分支的基本關聯系數為1。在風機性能函數p<0時，聯立式(12)，式(15)～(16)，式(19)～(20)，得式(30)：

(30)

(31)

仿照式(23)，可寫出給定風機性能函數下的通風網絡方程組，如式(32)所示：

(32)

將式(5)、式(8)、式(13)代入式(32)，獲得通風網絡方程擾動展開式。

(33)

(34)

4 擾動法解算通風網絡的流程

式(3)成立的基本條件是h>0，但在擾動法解算風網時需要計算lnh0，hx./h0，所以在計算過程中，當h0≤0時，會出現計算錯誤，從而導致解算失敗。下面分別討論如何處理此種情形。

1)h0<0

在通風網絡解算時，基本關聯矩陣、獨立回路矩陣等的選取依賴于事先假設的分支風流方向。如果事先假設的分支風流方向錯誤，會出現h<0的情況，進而h0<0。則式(23)～(25)中涉及到計算lnh0時會出現錯誤，導致計算失敗。出現這種情況時，修改原來假設的分支風流方向即可。反應在對應矩陣上，可以在式(22)解算結束后，根據h0的正負，利用代碼1更新矩陣M,N，h0。

>>Sig=diag((h0>=0)*2-1);

>>M=M*Sig;N=N*Sig;

>>h0=abs(h0);

代碼1

diag函數用來將向量展開為對角線矩陣。Sig為對角矩陣，當h0≥0時，對角線上相應位置為1，否則為-1。

2)h0=0

當風網的拓撲結構和分支參數處于某種特殊狀態(如對稱狀態)時，容易出現某些分支風量為0的情形。此時h=0，根據擾動理論，其所有擾動項系數均為0。在MATLAB中會導致無法計算lnh0，hx./h0，從而使解算失敗。這時，可以進行如代碼2所示處理，直接令擾動項系數為零，從而跳過對數和除法運算。

>>x0=h0,x0(～～x0)=1./x0(～～x0);

>>y=h0,y(～～y)=log(y(～～y));

代碼2

經過以上處理后，則hx./h0轉化為hx.×x0，需要計算lnh0的地方，替換為y。

綜上所述，利用擾動法解算通風網絡的流程如下：

1)分析通風網絡圖，找出連接節點、獨立回路和固定等級節點；

2)對網絡中所有的分支，根據R計算出α，寫出矩陣A。寫出對連接節點、固定等級節點的基本關聯矩陣B，Bf。寫出對獨立回路的獨立回路矩陣C和對固定等級節點的通路壓降矩陣F。寫出Hp，Qp等；

3)根據風機性能函數，計算參數ap，bp，cp；

4)將A，B，Bf，C，F，Hp，Qp，ap，bp，cp等構建矩陣M，N，P，T；

5)計算M的逆矩陣，并求出-M-1N；

8)根據代碼2處理后續的牽涉到計算lnh0，hx./h0的公式；

10)根據式(5)求解通風網絡各分支的通風阻力，根據式(3)求解各分支的風量。

5 計算實例

5.1 固定風壓工況點的計算實例

文獻[6]針對通風網絡1分風網絡利用擬牛頓法和直接解非線性方程法，在MATLAB中進行解算。此處風機工況點固定在500 Pa的風壓處。通風網絡1如圖1所示，其中，1～6表示支路，Ⅰ～Ⅴ表示節點。

圖1 通風網絡1

為驗證擾動法的有效性，根據事先擬定的分支風流流向，遵照風量連續方程和風壓連續方程，人為設定各分支風量Qr。指定部分分支風阻，推算出其他各分支風阻和所有分支風壓情況。各分支的風阻、風量和風壓見表1，單位取標準國際單位。此處的風量為對應風阻情況下的精確解。

表1 通風網絡1的相關參數

表2 擾動法和牛頓法解算效果比較

由表2可知，限定的精度下，擾動法的精度相當于牛頓法8次迭代時的精度，且用時遠小于牛頓法。隨著擾動項指數的增加，擾動法的計算誤差以0.5倍的速率遞減。但通過分析擾動法的計算原理可知，擾動項指數每增加1，需要計算的擾動項系數的個數增加1倍，理論上其計算時間也會增加。而牛頓法的每次迭代時間是相同的，并且迭代次數10次以內時誤差下降速率最大。從本文案例中可以看出，在分支誤差控制在1%左右時，擾動法用時遠小于牛頓法。

5.2 給定風機性能函數的求解實例

圖2 通風網絡2

根據分析，該風網有12個分支，6個連接節點，2個固定等級節點，5個獨立回路。故M，N均為13×13的矩陣，文獻[8]中計算出有2個分支風量為0，本文針對存在0風量分支的風網解算進行說明。

利用牛頓法進行驗算，證明用擾動法計算的精度同牛頓法計算20次時的精度相近，且所用時間遠小于牛頓法，具體驗算過程和產生的數據表格同5.1節類似。

6 結論

1)在擾動法系數求解的過程中，僅需要計算1次風網方程組系數矩陣M的逆，并且除擾動主項系數外，其他擾動項系數的-M-1N部分均相同，可反復利用，而括號內的部分是向量元素間的代數運算，故計算工作量遠小于其他傳統方法。

2)利用擾動法無需迭代運算，無需進行多次循環運算、判斷。并且其矩陣運算量相較于牛頓法的每次迭代都要小，在工程允許的精度范圍內，其耗時遠小于同精度下的牛頓法。

3)擾動法無需提供計算初值，克服了牛頓法的初值依賴性問題。

4)由于該方法中用到對數計算和除法運算，因此從數學角度講需要保證h≥0。本文通過增設MATLAB運算代碼來修改原來假設的分支風流方向，并對無風分支進行特殊處理，確保計算的成功。