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

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

1 引言

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

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

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

2 非重疊型區(qū)域分解法

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

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

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

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

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

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

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

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

3 并行計算實現(xiàn)

3.1 網(wǎng)格劃分

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

圖2 非結構化網(wǎng)格分區(qū)示意圖Fig.2 The sketch map of unstructured mesh partition

3.2 剛度矩陣的形成

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

圖3 鄰接圖Fig.3 The adjacency graph

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

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

3.3 數(shù)據(jù)通信

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

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

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

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

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

4 算例及分析

4.1 控制方程

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

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

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

其復數(shù)形式為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 六個分區(qū)網(wǎng)格圖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機構建的并行計算機系統(tǒng)，基于非重疊型區(qū)域分解法，采用消息傳遞的并行編程對等模式，對三維有限元計算的數(shù)學模型開發(fā)了SPMD的并行程序，實現(xiàn)了對500kV高壓輸電線路的絕緣子串電位分布和金具表面電場強度的并行計算。計算結果表明并行計算系統(tǒng)相對于單機計算可以加快計算速度，獲得良好的加速比。以合成絕緣子串電位分布的結果為依據(jù)來優(yōu)化屏蔽環(huán)位置，通過計算復合絕緣子串金具表面電場強度，以確定并降低絕緣子端部電場強度。

[1]司馬文霞, 楊慶, 孫才新, 等. 基于有限元和神經(jīng)網(wǎng)絡方法對超高壓合成絕緣子均壓環(huán)結構優(yōu)化的研究[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]. 電網(wǎng)技術, 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]沈鼎申, 張孝軍, 萬啟發(fā), 等. 750kV線路絕緣子串電壓分布的有限元計算[J]. 電網(wǎng)技術, 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]謝德馨, 唐任遠, 王爾智, 等. 計算電磁學的現(xiàn)狀與發(fā)展趨勢——第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 集群的大規(guī)模電力系統(tǒng)方程并行 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]吉興全, 王成山. 電力系統(tǒng)并行計算方法比較研究[J]. 電網(wǎng)技術, 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]呂濤, 石濟民, 林振寶, 等. 區(qū)域分解算法—偏微分方程數(shù)值解新技術[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.