特高壓直流輸電線路離子流場快速穩定計算方法

修連成，杜志葉，岳國華，何靖萱，蔡泓威，易凡

(1.武漢大學電氣與自動化學院，武漢430072；2.空軍預警學院雷達士官學校，武漢430345)

0 引言

近年來，隨著我國電能需求的不斷增加，特高壓直流輸電技術在電力系統中得到廣泛運用[1-3]。當特高壓直流輸電線路的導線和金具表面電場強度超過起暈電場強度時，其表面會發生電暈[4-7]，導線和金具周圍產生了正、負電荷。因為直流輸電線路導線極性不變，電荷受到電場力的作用向空間中擴散，與標稱電場作用形成了合成電場，提高了地面電場強度[8-10]。由于輸電電壓等級的不斷提高，線路附近的電磁環境問題更加嚴重。因此，特高壓直流輸電線路下方的離子流場研究具有非常重要的意義。

解析求解是一種可以精確計算離子流場的方法，其原理是根據場域及其內部物理量分布對稱特性，從而獲得與邊界條件對應的唯一解[11]。但是，該方法只能求解相對簡單的結構，如同軸圓柱模型、平行板模型等。為了解決復雜模型離子流場求解問題，文獻[12-13]基于Deutsch假設提出了一種離子流場計算方法，實現了復雜模型的離子流場快速求解。然而，該假設認為空間電荷只會影響電場數值，而不會改變電場方向，這與實際情況有所區別。而且在含有風速影響時，采用Deutsch假設的離子流場求解精度較低。文獻[14]提出一種利用Kaptzov假設求解離子流場的有限元方法。針對不同結構的直流輸電線路問題，文獻[15-18]利用基于Kaptzov假設的有限元法對其離子流場進行了求解和分析。文獻[19-20]通過有限元法計算了單極和雙極直流輸電線路的離子流場。有限元法求解過程主要分為兩個部分：首先，通過假設的初始電荷密度求解泊松方程，實現獲取空間電場信息。隨后，利用空間電場和電流連續性方程計算空間電荷密度。將上述兩部分反復迭代即可獲得地面合成電場和離子流場分布。但是該方法在求解復雜模型時經常會發生迭代過程不收斂的情況。

為了快速準確地求解特高壓直流線路離子流場，文獻[21]提出了一種上流有限元法來求解直流線路的離子流場。上流有限元法的原理是搜尋和確定上流有限單元，從而實現由電荷源點向外擴散求解的目標。基于上流有限元法和Kaptzov假設對輸電線路的合成電場和離子流密度進行了求解[22-23]。文獻[24]基于上流有限元法提出了一種減少剖分網格數量的方法。該方法將二維上流有限元法求解的離子流場當作三維求解時輸電導線的求解邊界。但是上述文獻對初始導體表面電荷敏感度較高，迭代過程容易無法收斂。文獻[25]提出了一種偏差估算法來計算初始導線表面電荷密度。但是該方法的迭代步數較多，迭代計算時間較長。

本文通過基于牛頓-拉夫遜的上流有限元法求解特高壓直流線路下方的地面合成電場與離子流場。該方法不僅無需調整和優化算法參數，還提高了收斂能力。該算法應用簡單，魯棒性強，可求解復雜情況下的地面合成電場和離子流密度。

1 離子流場模型

特高壓直流輸電線路離子流場模型中的正、負電荷主要受到了風力和電場力的影響。因此，離子流場的控制方程可以表示為：

(1)

j+=ρ+(k+E+W)

(2)

j-=ρ-(-k-E+W)

(3)

(4)

(5)

(6)

式中：φ為電位，V；ρ+、ρ-分別為正、負電荷密度，C/m3；ε0為空氣介電常數；j+、j-分別為正、負離子流密度，A/m2；k+、k-分別為正、負離子遷移率，m2/Vs；E為電場強度，V/m；W為風速，m/s；Rn為復合系數。

將式(1)—(6)代入式(4)和式(5)，可得：

(7)

式中V+(-)為正、負離子的遷移速度。

特高壓直流輸電線路下方的離子流場求解屬于強非線性問題，很難實現直接準確求解，因此需要簡化離子流場模型[26]。同時，假設條件應對地面離子流場計算結果影響較小，即

1)忽略特高壓直流輸電導線的電暈層厚度；

2)電暈已到達穩態，不考慮暫態過程；

3)導線表面場強滿足Kaptzov假設；

4)正、負離子的遷移率、復合系數和風速與電場無關，即設置為恒定常數。

2 基于牛頓-拉夫遜的離子流場求解算法

求解電流連續性方程是離子流場計算的重要環節。因此，本文通過上流有限元法求解電流連續性方程。上流有限元法的基本原理是通過三角單元的兩個已知電荷密度節點去計算另一未知電荷密度節點[5,7]。空間正電荷密度可描述為：

(8)

式中：下標i,j,m分別代表不同節點；bijm和cijm均為三角形單元的節點坐標函數；SΔ為三角形單元面積。

由式(7)—(8)可得特高壓直流輸電線路的正電荷密度。正、負電荷密度求解思路一致。因此，節點i處的正、負電荷密度為：

(9)

式中：

(10)

(11)

(12)

式(9)的求解結果實現了從源點電荷發出并擴散到整個場域，與實際情況相符，保證了電流連續方程的收斂性。

在迭代求解過程中，導線表面電場場強變化較大時會導致更新的導體表面電荷密度發生大幅度振蕩，進而使迭代過程不收斂。同時，導線表面電荷密度的初值設置對計算泊松方程和電流連續性方程時的收斂性影響較大。當電荷密度的初值選擇不恰當時，迭代過程容易發生振蕩不收斂的情況。為此，本文提出了一種基于牛頓-拉夫遜的迭代收斂算法，具體證明如下。

假設迭代電荷密度方程屬于單變量非線性方程，即可得：

f(ρi)=0

(13)

將式(13)改寫成式(14)。

f[ρi(k-1)+Δρi(k-1)]=0

(14)

式中：ρi(k-1)為節點i在第k-1次迭代時的電荷密度；ρi(k-1)和Δρi(k-1)分別為迭代電荷密度方程的近似解和真解誤差。

根據泰勒公式原理，式(14)可改寫為：

f[ρi(k-1)+Δρi(k-1)]=

f[ρi(k-1)]+f′[ρi(k-1)]Δρi(k-1)+

(15)

迭代電荷密度方程的真解誤差Δρi(k-1)相對較小，因此可以忽略其二階以上的高次項，即：

f[ρi(k-1)+Δρi(k-1)]=

f[ρi(k-1)]+f′[ρi(k-1)]Δρi(k-1)

(16)

根據式(14)和式(16)可知，真解誤差Δρi(k-1)的修正解可表示為：

ρi(k)=ρi(k-1)+Δρi(k-1)

=ρi(k-1)-f(k-1)/f′(k-1)

(17)

根據式(17)可知，改進算法依據迭代過程的前2步導線表面電荷密度來求解新的電荷密度。該算法在導線表面場強與起暈場強相差較遠(較近)時，電荷密度更新較快(較慢)，保證了算法的快速穩定收斂。因此，本算法對導體表面電荷密度的初始值不敏感。具體控制過程如下。

ρi(k)=ρi(k-1)+Δρi(k-1)=

(18)

式中：Emaxk-1、Emaxk-2分別為第k-1、k-2次迭代時導線表面的最大電場強度；E0為導線表面起暈場強。

式(18)解決了導體表面電荷密度快速收斂到真值附近問題。但是當計算空間電荷密度時，前后兩次迭代的空間電荷密度大幅度變化，可能會導致計算不收斂。因此，本文引入了迭代變化系數到空間電荷密度更新中，確保迭代穩定收斂。空間電荷密度可以表示為：

ρis(k)=(1-v)ρis(k-1)+ρisc(k)v

(19)

式中：ρis(k)、ρis(k-1)分別為上一步迭代過程中的空間電荷密度；ρisc(k)為當前迭代步中根據上一步有限元法計算得到的空間電荷密度；ν為迭代變化系數。

綜述所述，在迭代過程前期，輸電導線表面的電荷密度距離實際電荷密度相差較遠。因此，本算法加大了調整電荷量，確保迭代過程快速進行，提高迭代過程計算效率，使得導線表面電場強度可以快速接近起暈場強。在迭代過程后期，本算法減少了更新的電荷量，使得本算法可以穩定收斂。離子流場的收斂條件為：

(20)

(21)

式中：δE、δρ分別為導線表面的電場強度和電荷密度相對誤差；Emax為導線表面的最大電場強度。

當導線表面未發生電暈時，通過求解拉普拉斯方程獲得空間電場分布。而求解離子流場主要分為四個步驟：通過空間電荷密度求解泊松方程得到空間電場分布信息；利用空間電場信息和上流有限元法求解電流連續性方程；判斷是否滿足Kaptzov假設和電荷密度結果穩定的收斂條件；如果不滿足穩態收斂條件，則采用基于牛頓-拉夫遜的迭代收斂算法來修改電荷密度。將上述過程重復迭代到達收斂條件。離子流場求解的流程圖如圖1所示。

圖1 離子流場計算流程圖

3 算法驗證

3.1 試驗模型驗證

通過單極直流試驗線路模型來驗證本文所提出的新型迭代收斂算法的正確性和有效性。該模型的直流導線施加電壓為30 kV，導線半徑為0.8 mm，直流導線高度為0.394 m，試驗平臺如圖2所示。在導體模型下方放置場磨探頭和威爾遜板，用于測量地面電場強度和離子流密度。圖中左側的高壓發生器控制臺用于控制導線電壓等級。送風機用于控制試驗風速。

圖2 試驗平臺

參考國家標準GB/T 37543—2019中的評價依據對試驗數據進行了統計。圖3給出了不同風速下的本文算法與試驗測量數據的對比結果。圖3(a)和3(b)中的紫(深)色點和黃(淺)色點分別表示在風速為0 m/s和3 m/s時的測量結果。

圖3 計算和測量結果對比

根據圖3(a)，本算法求解的地面合成電場結果與試驗線路的測量結果相差較小。特別是在風速值為0 m/s時，計算值與測量值基本一致。本算法的離子流密度計算值與實際測量值對比如圖3(b)所示。實際測量峰值下降速度比本文方法的計算值下降速度稍快。但計算值與實際測量值基本一致，證明了本算法的正確性。需要說明的是，風速對離子流密度和合成電場影響明顯。電荷由逆風側被吹至順風側，合成電場最大值位置也向順風側發生偏移。

3.2 算法收斂性與穩定性分析

本文利用云南—廣東± 800 kV特高壓直流輸電線路模型來驗證本算法的收斂性和穩定性。該線路導線數據詳見參考文獻[27]。當初始導線表面電荷不同時，本算法求解離子流場的迭代收斂結果如圖4所示。圖中導線表面初始電荷密度分別為1×10-8、5×10-8、1×10-7、1.5×10-7和2×10-7C/m3。顯然，本算法在不同初始導線表面電荷密情況下都能夠快速穩定收斂，且迭代步數都少于20步。

圖4 收斂效果驗證

當初始導線表面電荷為2×10-7C/m3時，將本文算法、收斂控制因子法以及傳統方法進行了對比，如圖5所示。傳統方法對電荷敏感度較高，容易出現迭代振蕩，從而導致計算無法收斂。迭代收斂控制因子法和本文算法都能夠實現穩定的迭代收斂。

圖5 不同方法收斂速度對比

迭代收斂控制因子法對初始導線表面電荷密度不敏感，但該方法需要設置控制因子，且控制因子值的大小與收斂性和計算速度直接相關。當控制因子較小時，迭代過程中的電荷密度變化不大，進而實現穩定收斂。但是，其迭代速度較慢，計算時間變長。當控制因子較大時，迭代收斂速度變快，然而導體表面電荷變化較大，可能導致無法收斂。由圖5可知，迭代收斂控制因子法雖然能夠較快速到達電荷真值附近，但是仍需要多次迭代才能達到收斂條件。而本文所提出的算法迭代速度較快，且無需復雜的算法參數調整與優化。特高壓直流輸電線路的導線表面電荷和表面電場強度變化情況如圖6所示。

圖6 導線表面電荷和表面電場強度變化曲線

本文算法在迭代過程初期的導線表面電荷密度變化程度較大，而迭代后期的電荷變化較小，使得本文算法可以快速穩定迭代收斂。迭代初期(3～10步)，導線附近合成電場強度變化劇烈，使得導體表面電荷變化率較大，從而實現快速地迭代到起暈場強附近。迭代后期，導線表面電荷變化率較小，確保離子流場迭代穩定收斂，防止在起暈場強附近發生振蕩。因此，本文算法實現了快速穩定地計算特高壓直流線路下方離子流場。

4 風速對±800 kV特高壓直流輸電線路離子流場影響分析

本文以± 800 kV特高壓直流輸電線路為例進行了地面離子流場求解，得到了地面合成場強和地面離子流密度在不同風速下的分布規律。具體線路參數和測量數據見參考文獻[27]。圖7和圖8分別表示了不同風速時地面合成電場強度和離子流密度。

圖7 地面合成場強求解結果

圖8 地面離子流密度求解結果

圖7—8中風向從左側向右側為正，風速從-3 m/s增加到3 m/s。當風速從-3 m/s增加時，合成場強和離子流密度的最大值位置逐漸向右側移動，且其最大值先減小后增加。本文計算的合成電場和離子流密度與測量結果基本一致，驗證了本算法的有效性和正確性。

隨著風速向左增加，正極性導線下方地面合成場強和離子流密度逐漸增大，而負極性導線下方的離子流場有所減小。當風速值從0 m/s增加時，負極性下方的合成電場和離子流密度逐漸減小。這是由于負極性導線下方的空間電荷都被吹到了正極性導線下方，從而導致負極性導線下方的合成電場和離子流密度減小。當風速值為3 m/s時，最大合成場強為36.58 kV/m，且最大離子流密度為83.6 nA/m2。同時，對比于風速為0 m/s時的最大合成電場，其位置發生了約5 m的偏移。綜上所述，隨著風速向右(向左)增大，合成場強和離子流密度的分布曲線和最大值都逐漸向右(向左)移動，而且順風側導線的合成場強和離子流密度數值均會顯著變大。

本文算法在不同風速情況下求解合成電場和離子流場時，迭代過程一般少于20步，相比于傳統方法縮短了計算時間，且滿足工程精度要求。同時，本文所提出算法不但對導體表面電荷密度初值選擇不敏感，而且無須調整和優化算法參數，具有較好的魯棒性。

5 結語

為了解決特高壓輸電線路離子流場計算的快速穩定收斂問題，本文基于牛頓-拉夫遜法提出了一種快速穩定求解離子流場算法。該算法在導線表面場強距離起暈電場強度相差較大(較小)時，電荷密度更新較多(較少)，從而實現了快速穩定收斂。通過對實驗室模型線路的地面合成電場和離子流場計算，驗證了本算法的有效性和正確性。同時，本文對±800 kV 特高壓輸電線路下方離子流場進行了求解，證實了該算法可以提高求解速度并保證收斂的穩定性。隨著向右風速的增加，正極性導線下方的地面合成場強和離子流密度逐漸增大，而負極性導線下方的離子流場數值逐漸減小。該算法不但對導體表面電荷密度初值選擇不敏感，而且無需調整算法參數，具有較好的魯棒性。研究結論有助于快速準確計算特高壓直流輸電線路下方離子流場。