并行計算合成絕緣子串電壓分布及金具表面電場強度

厲天威 阮江軍 杜志葉 黃道春

（1. 廣西電力工業勘察設計研究院 南寧 530023 2. 武漢大學電氣工程學院 武漢 430072）

1 引言

合成絕緣子相對于玻璃和瓷絕緣子具有重量輕、強度高、耐污性能好、便于運輸和安裝、運行維護方便等優點，在國內外得到了廣泛的應用[1]。但在運行過程中有出現芯棒脆斷的情況，分析表明主要是高壓端部電場強度較高，電暈嚴重，產生電腐蝕并加劇酸腐作用的結果[2]，所以必須降低高壓端的電壓分布。隨著電網超高壓的發展，在線路上發生污閃事故的絕緣子幾乎都損壞了均壓環、端部金具及金具附近的傘裙、護套。這些被損壞的絕緣子有很多是由于絕緣子端部區域的場強過高引起的，給電網的安全運行帶來極大的危害。計算合成絕緣子串的電位分布以及金具表面的電場強度并以此來采取防護措施勢在必行。文獻[1]和文獻[3]用三維有限元法計算高壓輸電線路的絕緣子串的電位分布，由于單機計算速度和內存大小的限制，網格剖分數量有限，因此計算結果會存在一定的誤差。本文通過并行計算進行分析及對比，來得到復合絕緣子串電位分布和金具表電場強度，以確定絕緣子串屏蔽環的最優位置并降低絕緣子串端部電場強度。

有限元并行計算已經成為解決大規模數值計算問題的最重要手段之一[4-6]。有限元并行最直接和最初的實現是利用數值并行求解器對傳統有限元最耗時的線性方程組求解過程進行并行化求解以加快整個分析過程[7-11]，即線性方程組的求解是并行過程，而單元分析和線性方程組的形成是串行過程。為進一步提高并行效率，本文提出采用區域分解法[12-15]“分而治之”的并行策略在分布式存儲的并行環境下實現有限元并行計算，即將有限元求解區域劃分成若干大小近似相等的分區，然后將這些分區信息映射到各個處理器上并行地進行計算。每個處理器負責該分區的全部計算，通過消息傳遞方式與其他處理器交換分區邊界數據信息。該并行策略在單元級上實施有限元并行求解，避免形成總體剛度矩陣，其優點在于計算高度并行化和數據通信量少。

本文建立了500kV高壓輸電線路合成絕緣子串的三維有限元模型，運用基于非重疊型區域分解的并行計算技術對絕緣子串的電位分布以及金具表面的電場強度進行計算來優化屏蔽環位置，確定并降低絕緣子端部電場強度。

2 非重疊型區域分解法

根據分區是否重疊，區域分解法分為重疊型區域分解法和非重疊型區域分解法。兩種分區方法如圖1所示。

圖1 重疊型和非重疊型區域分解Fig.1 The overlapping and non-overlapping domain decomposition

以Schwarz交替法[12]為基礎的重疊型區域分解算法的基本思想是把求解區域Ω 分解為兩個或多個較小的相互重疊的分區Ω1，Ω2，…，Ωn，然后在這兩個或多個分區之間交替求解，最終得到整個求解區域的解。Schwarz交替法的邊值問題為

其中，aij為對稱正定矩陣，a0≥0。把求解區域Ω 分解為 M個分區，當Ωk與Ωm相鄰時，Ωk∩ Ωm≠ ? 。選擇初始近似 u0∈(Ω )，n=0，則并行計算各分區上邊值問題為

再令n=n+1繼續迭代，直至收斂。

基于Schwarz交替法的重疊型區域分解法，各分區與其他分區有重疊，重疊區域中的未知量與多個分區相關聯，因而計算量相對較大。其計算效率取決于重疊區域的大小，當分區個數增多時迭代收斂速度減慢，特別對于三維問題，增加重疊區域嚴重制約了計算速度的提高。

對于橢圓形偏微分方程的邊值問題，非重疊型區域分解算法的邊值問題為

如何在減少重疊區域的同時不降低迭代收斂速度成為本文所關注的問題。對于非重疊型區域分解，只有區域之間的邊界上的變量與多個分區相關聯，迭代計算速度更快。因此，本文采用非重疊型區域分解法來進行并行求解。

3 并行計算實現

3.1 網格劃分

對于有限元網格如圖2a所示，利用METIS（分區軟件）對該網格劃分了3個分區。分區結果如圖2b所示。其中各分區的單元的尖點（部分節點）都在與之相鄰的分區上，在實行有限元并行求解時需要利用消息傳遞的方式來更新這些節點的值。在圖2c中，處理器Proc 0在進行計算時，除了包括該分區里的所有節點和單元外，還必須包括與這些節點相關聯的外邊界單元的信息，如虛線所示的單元。因此，在圖2b的單元分區結果的基礎上，進行的分區后處理的結果如圖2d所示。對于圖2d中的三個分區，各處理器需要計算的節點在圖中的分界線以內，各分區相關聯的相鄰分區中的節點在分界線以外并靠近分界線，各分區的單元為標有同樣標記的節點構成的所有單元之和，如圖打點標記。

圖2 非結構化網格分區示意圖Fig.2 The sketch map of unstructured mesh partition

3.2 剛度矩陣的形成

有限元分區的鄰接圖如圖3所示，整個求解區域被劃分成3個分區，分別分配給3個處理器Proc 0、Proc 1和Proc 2進行并行計算。

圖3 鄰接圖Fig.3 The adjacency graph

若用單處理器對整個求解區域進行求解，需要形成總體剛度矩陣。而在各分區內單獨形成剛度矩陣時，各分區剛度矩陣分別對應于總體剛度矩陣的多行，如圖4所示。如Proc 0、Proc 1和Proc 2分別對應總體剛度矩陣的{0,1,3}行、{4}行和{2,5}行。

圖4 分布式系數矩陣的形成Fig.4 The formation of distributed coefficient matrices

3.3 數據通信

MPI[16]（消息傳遞函數庫）提供了與C、Fortran和C++語言的綁定，適用于分布內存的并行計算機系統。MPI程序采用的是SPMD（單程序多數據流）執行模式，即一個程序同時形成多個獨立的進程，各進程由獨立操作系統調度，享有獨立的CPU和內存資源。進程間通過消息傳遞互相交換信息，消息傳遞快可以減少各處理器間的等待時間，從而提高并行度和并行計算效率。

為了充分利用主節點的高計算性能，并且考慮到利用相同的數值計算方法在單機上計算結果有助于對比分析并行計算的效率，本文采用對等模式來替代主從模式，對等模式基本框架如圖5所示。對等模式即程序中各進程（包括控制部分）地位相同，只是處理的數據不同，且程序可以利用單個進程運行。由圖可知，進程0除了發送和接受數據外還承擔著計算的任務。

圖5 對等模式并行程序基本框架Fig.5 The basic frame of equal mode parallel program

各處理器間的數據交換主要表現為：在前處理階段，主處理器把各分區的單元節點信息發送到各從處理器，數據量較大，只需要傳遞一次；在迭代開始后，各主、從處理器把新的邊界節點值傳遞到相鄰處理器來更新邊界節點值，數據量相對較小。這些處理器之間的數據交換通過調用MPI庫函數來實現，采用非阻塞通信的方式來減少數據交換時間，以提高數據傳輸效率。網格分區后的各個分區的單元和節點數據均存儲在各處理器所在存儲器中，實現存儲局部化，這樣每個分區可以很方便存取需要的數據。各分區的剛度矩陣、荷載向量、分區的邊界節點值和中間計算結果都存儲在處理器所在存儲器中，這樣大大降低各處理器之間通信。

本文采用并行預處理共軛梯度法等對各分區線性方程組進行迭代求解，每一步迭代都要與相鄰的分區交換更新數據，調用 MPI函數庫中提供的MPI_Send，MPI_Recv等函數可以方便地實現數據交換。

4 算例及分析

4.1 控制方程

輸電線路通電情況下其周圍工頻電場變化頻率很低，屬準靜態場范疇。當電場的變化頻率很低時，庫侖電場遠大于感應電場，即0Bt??≈，于是有

若只考慮電場，則用電流連續性方程代入式（5），且不考慮式（7），則準靜態電場的控制方程可表示為

而J=σE，D= ε E，且有E =-?φ ，故

其復數形式為

假設電介質都為均勻介質，于是本文計算靜態電磁場的泊松方程可表示為

4.2 并行計算合成絕緣子串電位分布

為了給合成絕緣子的設計運行提供參考，本文運用有限元并行計算方法來計算 500kV高壓輸電線路合成絕緣子串上的電壓分布。絕緣子型號為FXBW2—500/160，大傘43片，小傘42片，結構高度為 4450mm，絕緣子介電常數為 3.5，兩端金具長度分別為150mm，屏蔽環外徑280mm；采用4分裂導線LGJ—400/50，分裂間距為450mm，導線長 20m；橫擔寬、高各 0.8m，長 16m。考慮到實際線路結構的對稱性和計算機的計算能力有限，以導線側芯棒端部中心為坐標原點，只對絕緣子串、金具、導線和橫擔的1/4部分建立三維模型進行分析，最后用空氣包圍起來，模擬無窮大區域。圖6a為合成絕緣子串的三維模型，圖6b為該模型外包空氣層后的1/4模型。為驗證本文算法正確性，對圖6b的三維模型采用4節點的四面體單元分別剖分為27 868和115 291個節點，基于區域分解法的思想，對該單元網格進行分區，來進行并行求解。

圖6 合成絕緣子串三維分析模型Fig.6 The 3D model of composite insulator string under the transmission line and the rung

在本研究室的機群環境下（6節點，Fedora core 2.0操作系統），并行計算程序采用MPI消息傳遞的并行編程的對等模式，數值求解方法用C++語言編寫實現，本文在芯棒中間位置均勻選取40個有代表性的點，精確計算得到各點的值作為參考。圖7為各絕緣子的耐受電壓百分比，橫坐標1～40是絕緣子串上依次從高壓端到低壓端上點的編號，可以看出并行計算和單機計算的結果能夠很好地吻合，從而驗證了本文運用以上并行求解算法計算的正確性。

圖7 并行計算與單機計算結果比較Fig.7 Comparison between sequential and parallel computation

運用各種線性方程組并行求解器PCG（預處理共軛梯度法）和 GMRES（廣義最小余量法）等分別對各分區內形成的線性方程組進行并行求解，求解器的計算時間與處理器數目之間的對比如圖8所示。可以看出，隨著處理器數目增加，各求解器計算時間在不斷減少。在求解精度相同的情況下，BICGSTAB（穩定雙共軛梯度法）法計算時間最少，CGS（復共軛梯度平方法）法和 GMRES法次之，CG（共軛梯度法）法計算時間最長。

圖8 各求解器計算時間與處理器數目的對比Fig.8 Time comparison for the number of processors

并行計算加速比和效率是衡量并行計算機系統和并行算法優劣的兩個基本指標。圖9a為利用GMRES計算得到的并行計算加速比和效率。可以看出，當計算量較小時，六處理器加速比為2.2。因為數據量較少，數據交換時間在整個計算過程中占的比重較大，導致加速比不高。計算量增大時，隨著處理器數目的增多，計算時間不斷地減少，出現了超線性加速情況，即 s(n)＞n（s(n)為加速比，n為處理器數）。超線性加速是由于在機群中每個處理器可以訪問的內存容量與單處理器擁有的內存容量相同，因此，在任何時刻總能保存更多的數據，減少了對工作速度相對較慢的硬盤的訪問。從圖9b可以看出，并行計算效率隨著處理器的增加有減小的趨勢。根據以上分析，適當增大處理器數目可以減少計算時間，并保證較高的求解效率。因算法特點不同，各求解器迭代次數也不同，迭代次數越多，處理器間數據交換越頻繁，對于27 868自由度規模的問題，各求解器迭代次數見下表。

圖9 并行計算加速比和求解效率Fig.9 The efficiency of parallel computation

表 各求解器迭代次數比較Tab. The iterative numbers of the solvers

4.3 并行計算金具表面電場強度

基于4.2節并行算法正確性的驗證，現以500kV高壓輸電線路復合絕緣子為例，考慮桿塔、導線和屏蔽環等因數的影響，對金具表面的電場強度進行并行計算。根據對稱性建立有桿塔的四分之一模型，采用4分裂導線，如圖10a所示，網格剖分如圖10b～圖10c所示，高壓端電位分布如圖10d所示。本文首先對整體網格模型進行分區，分區數目為2～6個，然后分別利用 2～6個處理器來進行計算。其中，6個分區的分區結果如圖11所示，分區過程綜合考慮了負載平衡等因素，各分區計算量大致相等。

首先對絕緣子周圍電場分布進行計算，研究對象的取值點如圖12右，在每片絕緣子的大小兩個傘裙之間、小傘的下方處取點，共42個。由圖12左可知單機和并行求解的結果能夠很好地吻合，在上、下屏蔽環附近，電場強度有急劇增大的趨勢，可知絕緣子周圍電場強度的最大值在下屏蔽環附近，大約為920kV/m，未達到均勻電場中擊穿場強的峰值(3000kV/m)。隨后，本文依次對絕緣子上、下鐵帽和導線周圍線夾表面電場強度分別進行了計算，計算結果和取值點分別如圖13所示，為進一步分析電暈特性等提供了依據。在上鐵帽和下鐵帽周圍0.09m處取21個點。在離線夾所在平面0.02m處分裂導線周圍取21個點，這些點均勻分布在以分裂導線中心為圓心，半徑為1.2R的半圓上，R為四分裂導線分裂半徑。

圖10 500kV輸電線路模型及網格圖Fig.10 The model of 500kV transmission lines and the mesh

圖11 六個分區網格圖Fig.11 The six partitions mesh graph

圖12 絕緣子串周圍電場分布Fig.12 The electric field around the insulator

圖13 金具表面電場分布Fig.13 The electric field intensity distribution on the surface of hardware fittings

5 結論

本文利用多臺 PC機構建的并行計算機系統，基于非重疊型區域分解法，采用消息傳遞的并行編程對等模式，對三維有限元計算的數學模型開發了SPMD的并行程序，實現了對500kV高壓輸電線路的絕緣子串電位分布和金具表面電場強度的并行計算。計算結果表明并行計算系統相對于單機計算可以加快計算速度，獲得良好的加速比。以合成絕緣子串電位分布的結果為依據來優化屏蔽環位置，通過計算復合絕緣子串金具表面電場強度，以確定并降低絕緣子端部電場強度。

[1]司馬文霞, 楊慶, 孫才新, 等. 基于有限元和神經網絡方法對超高壓合成絕緣子均壓環結構優化的研究[J]. 中國電機工程學報, 2005, 25(17): 115-120.Sima Wenxia, Yang Qing, Sun Caixin, et al.Optimization of corona ring design for EHV composite insulator using finite element and neural network method[J]. Proceedings of the CSEE, 2005,25(17): 115-120.

[2]張鳴, 陳勉. 500kV羅北甲線合成絕緣子芯棒脆斷原因分析[J]. 電網技術, 2003, 27(12): 51-53.Zhang Ming, Chen Mian. Anlysis on core brittle fracture of composite insulators in a certain 500kV transmission line[J]. Power System Technology, 2003,27(12): 51-53.

[3]沈鼎申, 張孝軍, 萬啟發, 等. 750kV線路絕緣子串電壓分布的有限元計算[J]. 電網技術, 2003, 27(12):54-57.Shen Dingshen, Zhang Xiaojun, Wan Qifa, et al.Calculation of voltage distribution along insulator strings of 750kV by finite element method[J]. Power System Technology, 2003, 27(12): 54-57.

[4]謝德馨, 唐任遠, 王爾智, 等. 計算電磁學的現狀與發展趨勢——第14屆COMPUMAG會議綜述[J].電工技術學報, 2003, 18(5): 1-4.Xie Dexin, Tang Renyuan, Wang Erzhi, et al. Major achievements and future trends of computational electromagnetics—summary of the 14th conference of compumag[J]. Transactions of China Electrotechnical Society, 2003, 18(5): 1-4.

[5]鐘建英, 林莘, 何榮濤. 自能式 SF6斷路器熄弧性能的并行計算仿真[J]. 中國電機工程學報, 2006,26(20): 154-159.Zhong Jianying, Lin Xin, He Rongtao. The parallel computation simulation of the extinguishing arc performance in the self-extinguished SF6circuit breaker[J]. Proceedings of the CSEE, 2006, 26(20):154-159.

[6]Butrylo B, Musy F, Nicolas L, et al. A survey of parallel solvers for the finite element method in computational electromagnetics[J]. The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 2004, 23(2):531-546.

[7]劉洋, 周家啟, 謝開貴, 等. 基于 Beowulf 集群的大規模電力系統方程并行 PCG 求解[J]. 電工技術學報, 2006, 21(3): 105-111.Liu Yang, Zhou Jiaqi, Xie Kaigui, et al. Parallel PCG solution of large scale power system equations based on Beowulf cluster[J]. Transactions of China Electrotechnical Society, 2006, 21(3): 105-111.

[8]Watanabe K, Igarashi H. A parallel PCG method with overlapping domain decomposition[C]. 12th Biennial IEEE Conference on Electromagnetic Field Computation,2006: 49.

[9]Ito F, Amemiya N. Application of parallelized SOR method to electromagnetic field analysis of superconductors[J]. IEEE Transactions on Applied Superconductivity, 2004, 14(2): 1874-1877.

[10]Vollaire C, Nicolas L. Preconditioning techniques for the conjugate gradient solver on a parallel distributed memory computer[J]. IEEE Transactions on Magnetics, 1998, 34(5): 3347-3350.

[11]吉興全, 王成山. 電力系統并行計算方法比較研究[J]. 電網技術, 2003, 27(4):22-26.Ji Xingquan, Wang Chengshan. A comparative study on parallel processing applied in power system[J].Power System Technology, 2003, 27(4): 22-26.

[12]呂濤, 石濟民, 林振寶, 等. 區域分解算法—偏微分方程數值解新技術[M]. 北京: 科學出版社,1999.

[13]Imre Sebestyén. Electric-field calculation for HV insulators using domain-decomposition method[J].IEEE Transactions on Magnetics, 2002, 38(2): 1213-1216.

[14]Liu Peng, Jin Yaqiu. The finite-element method with domain decomposition for electromagnetic bistatic scattering from the comprehensive model of a ship on and a target above a large-scale rough sea surface[J].IEEE Transaction on Geoscience and Remote Sensing,2004, 42(5): 950-956.

[15]Liu Y Q, Yuan J S. A finite element domain decomposition combined with algebraic multigrid method for large-scale electromagnetic field computation[J]. IEEE Transactions on Magnetics, 2006,42(4): 655-658.

[16]都志輝. 高性能計算并行編程技術: MPI并行程序設計[M]. 北京:清華大學出版社, 2001.