電力高壓設備高性能有限元電磁仿真方法綜述與展望

鄒軍，楊帥，程啟問

（清華大學電機系， 北京 100084）

0 引言

構建以新能源為主體的新型電力系統是實現中國“30·60”“雙碳”戰略目標的基礎支撐。新型電力系統包含交流線路、直流線路、高壓電力設備和大量新能源接入的電力電子裝置。由于新型電力設備和電力電子裝置的接入，電力系統中將存在大量時變的和非線性的“元件”，電力系統運行特性變得更為復雜［1-3］。為模擬電力系統的特性，亟須逐步從器件級、變換器級、系統級開展不同尺度建模和仿真，這些仿真的基礎之一是對高壓電氣設備和裝置的電磁行為特性進行高精度建模。

有限元方法是一種通用的偏微分方程求解技術，其對于任意幾何形狀和材料都能夠穩定獲得數值解，為實際工程實踐提供定量化的指導［4］。科學和技術領域中對電磁建模的需求不斷增長，催生了以有限元方法為基礎的商業軟件。商用電磁仿真軟件具有友好的人機接口，可與計算機輔助設計（computer aided design，CAD和計算機輔助工程computer aided engineering，CAE）的軟件結合，在實踐中獲得了巨大的成功。在電磁仿真領域，借助各種變換、因子分解和對角化技術，求解數十億自由度的稠密矩陣已經成為可能。新型寬帶求解器、等效方法和區域分解法等計算技術可以對多尺度問題進行全波分析，而不必對實際問題做過多的簡化。有人認為數值計算機技術已經“完結”，無需研究新方法和技術，所有的電磁問題可以采用商業軟件和集成更高算力的計算機/工作站最終得以解決。但是，由于新應用場景不斷出現，例如納米電磁學、手征材料和強非線性多物理耦合等場合，商業軟件功能局限性日益突出，對于計算規模在數千萬單元、瞬態、非線性和多物理場耦合問題，仿真能力和手段還非常有限。對于所仿真的系統往往通過增加硬件開銷的方式采用蠻力進行計算，計算效率無法保證。因此，對有限元方法和加速技術的研究重新被學術和工業界人士所重視。

本文聚焦于有限元方法，從高性能有限元技術的發展路徑出發，討論有限元技術發展的方向和當前的瓶頸問題。同時，結合電力系統電磁仿真的特點，簡要討論了計算電磁學的關鍵技術問題和可能的發展路線圖。

1 工程電磁場有限元方法回顧與挑戰

在20世紀初，不同研究領域的科學家都曾獨立提出了對微分方程的變分形式進行離散的求解思路［5］，其中德裔美國數學家Courant的工作被認為是有限元方法的起源：Courant在分析軸的扭轉問題時將軸的橫截面離散成若干三角形單元，并用每個節點上未知量的線性插值表示截面上的應力分布。在20世紀50年代初，波音公司的工程師獨立發展出了基于能量極小原理的“直接剛性方法”，并廣泛應用于機翼結構的建模仿真［6］。20世紀年代初，為了突出在有限單元上進行積分運算的特性，Clough正式提出了“有限元方法”這一名詞。作為Rayleigh-Ritz方法的一種形式，有限元分析開始被學術界廣泛接受，不同領域的學者嘗試在本領域中使用有限元方法，應用數學家也開始著手建立有限元方法的數學理論。1962年至1972年被稱為有限元方法發展的黃金十年［7］：有限元的求解對象從橢圓方程拓展到雙曲方程和拋物方程，應用場景從結構力學拓展到熱傳導和流體分析；經典的有限元技術諸如分片測試、協調單元和非協調單元、等參變換在這一時期被提出；針對有限元方法的誤差分析和收斂性分析理論被建立；在這一時期也涌現了大量經典有限元計算程序，如NASA發展的結構有限元分析軟件NASTRAN、John Swanon發展的通用有限元分析軟件ANSYS等。為了解決傳統變分理論的有限元方法在求解時出現的數值振蕩問題，Zienkiewicz等人提出了基于迎風格式的Petrov Galerkin有限元方法［8］；在此基礎上針對暫態問題發展出了與波傳播特征相結合的Taylor Galerkin有限元方法［9］；1990年前后，Cockburn和舒其望等人提出了具有通量守恒特性的離散Galerkin有限元方法［10］。

隨著計算電磁學研究對象越來越復雜，模型越來越精細，有限元離散后生成的代數方程組規模愈來愈大，這對現有有限元方法提出了新的挑戰。以高電壓工程中的流注發展問題為例，受限于CFL條件，迭代的時間步僅為數皮秒，實際仿真中需要進行上千步迭代，而每一個時間步需要求解一個涉及上億自由度的問題［11］；在仿真計算高壓絕緣子表面電位時，絕緣子的幾何尺度遠小于桿塔的幾何尺度，導致整個計算區域離散后產生數百萬個節點［12］；大地電阻率反演問題中，采用精細的像素級剖分后將產生數百萬待反演參數，在反演過程中需要使用有限元方法反復求解正問題［13］。由于特高壓工程的推進，大量超常規的高壓設備也需要重新進行設計：大容量電力電子器件中電-熱-力三種物理場相互耦合，共同決定了開關性能；針對高壓套管中電場優化需要解決分層介質導致的多重尺度問題；在電力變壓器中各類軟磁材料引入的強非線性對數值仿真的穩定性和精度提出挑戰。總結來看，電力系統高壓設備電磁仿真中存在大量三維、時域以及非線性問題亟待解決，為此，中國學者提出了“計算高壓學”的概念，而有限元是其核心支撐技術之一。在新型電力系統中，有大量電力電子設備的應用，考慮到整個系統的電磁特性，計算空間尺度可能會跨越3個數量級以上，目前還不能完全采用場的方法仿真。

2 高性能計算有限元方法的應用進展

2.1 高性能計算硬件發展

高性能計算，即使用盡可能多的硬件資源完成特定問題的計算。將一個串行計算程序轉化成并行程序并移植到高性能計算機上運行，涉及兩類操作：分解和映射［14］：將計算任務分解成一組按照任務依賴圖并行執行的任務，然后映射到各個計算單元上執行。為了在最短運行時間內執行完相關任務，需要根據實際計算單元特點將任務分解成特定粒度的子任務，同時最小化單元間通信開銷，減少不同計算單元之間的交互。由此可見，高性能計算技術的應用依托于計算平臺本身。

在處理器發展早期，計算性能的提升主要靠提升主頻和處理器內部的隱式并行［15］：提升主頻，即縮短處理器的時鐘周期，從而提升單位時間內浮點數運算能力；處理器內部隱式并行則是通過流水線技術結合超長指令實現指令級并行，提高硬件利用率，進而提升計算性能。由于單個處理器計算性能低，內存小，為了獲得更高算力，需要將大量處理器通過高速通信網絡連接，構成大規模并行處理機（mas parallel processor， MPP）［16］。在這一階段，主要是通過區域分解法將大規模問題分解成若干小規模的子問題并分配給各個結點進行計算，再通過結點間的通信傳遞不同子問題間的相互影響，最終實現對原問題的求解。區域分解法本質上是通過犧牲計算時間，使用有限計算資源求解大規模計算問題，當問題規模較小時，使用區域分解法將帶來近十倍的耗時增長［17］。

隨著芯片制造工藝發展和指令集的更新，單核處理器性能逐漸逼近理論上限，為進一步提升處理器性能，多核（multi-core）構架應運而生：處理器由多個獨立的計算核心構成，各個計算核心通過共享高速緩存與多線程技術相結合，實現更高計算性能。多核構架受限于高速緩存一致性，隨著核數增加，緩存信息的維護成本越來越高，難度越來越大。為進一步提升處理器計算性能，在多核構架基礎上發展出了眾核（many-core）構架，即放棄高速緩存一致性，通過堆疊更多計算核心來提升算力。典型的眾核處理器有通用圖形處理器（GPU）、志強協處理器（Xeon Phi）等。基于不同設計理念的構架具有完全不同的特性，相應的處理器適用于處理不同的問題［18］：1）多核處理器：適合順序執行復雜指令，串行計算性能高；由于線程切換開銷大，并行計算性能較弱。2）眾核處理器：適合執行單指令多數據（single instruction multi-data，SIMD）或單指令多線程（single instruction multi-thread，SIMT）操作，并行計算性能高；由于每個核上寄存器數量較少，不適合執行復雜指令。

隨著眾核處理器的出現，有限元方法能夠進行更細粒度的分解，實現單元級并行，即單個核心（或單個線程）處理一個單元。隨著問題規模增大，GPU的顯存成為限制求解規模的主要瓶頸［19］：以2020年最新發布的A100顯卡為例，其最大顯存僅為40 GB。為解決這一問題，研究人員提出了多GPU聯合仿真的解決方法。如圖1所示，將整體剛性矩陣分塊，每塊GPU負責一塊子矩陣的生成和裝配；求解矩陣方程時，各個GPU并行完成矩陣-向量乘法，將計算結果傳遞給CPU，再由CPU進行組裝并分發給各個GPU。多GPU聯合仿真時，由于引入了CPU進行收據收集和分發，需要考慮由此產生的CPU-GPU通信開銷：隨著并行GPU的數量增加，從CPU向GPU分發數據時會造成數據總線阻塞，相應地通訊開銷會成倍增長。如表1—2所示，當GPU數量從2塊增長為4塊時，從CPU向GPU分發數據的時間從4 ms增長為8 ms［20］。

表1 塊GPU聯合仿真時各類操作的耗時 (GPU數量為2)Tab.1 Time-consuming of various operations in block GPU co-simulation (the number of GPU is 2)

圖1 多GPU聯合仿真模式Fig.1 Multi-GPU co-simulation mode

目前主流服務器、工作站往往同時提供多核處理器和眾核處理器，而在上述提到的研究中，計算任務完全由GPU執行，CPU只起到居中協調的作用。為了進一步提升算力，需要讓閑置的CPU也參與計算。為實現CPU和GPU聯合的異構計算，需要設計合理的任務分配機制。研究人員提出了兩類任務分配策略：靜態分配——在任務開始前將任務分配給各個計算單元；動態分配——在任務執行過程中將任務分配給各個計算單元。

表2 塊GPU聯合仿真時各類操作耗時 (GPU數量為4)Tab.2 Time-consuming of various operations in block GPU co-simulation (the number of GPU is 4)

如圖2所示，由于CPU向GPU的通信開銷以及GPU的準備時間，在問題規模較小時使用CPU執行的耗時要小于GPU執行的耗時。基于此，文獻［21］制定了靜態的任務分配策略：當問題規模小于設定閾值時，只使用CPU；隨著問題規模增大，將大部分任務轉移到GPU上執行，余下少量任務留在CPU上執行。在特定問題規模之下，這種靜態分配策略比僅使用GPU的策略快1.5倍。

圖2 不同問題規模下CPU和GPU計算時間對比Fig.2 Comparison of CPU and GPU computing time under different problem scales

文獻［22］在利用時域有限元仿真時采用了動態任務分配策略：根據上一個時間步的各個計算單元的計算耗時，調整下一個時間步各個計算單元的任務量，從而實現動態負載均衡。最終計算結果如表3所示，可以看到在動態任務分配機制下，CPU+GPU的異構計算性能優于僅使用CPU或GPU的計算模式。

表3 不同硬件平臺上，GPU-CPU聯合仿真加速比Table 3 GPU-CPU co-simulation speedup ratio on different hardware platforms

單個計算機的算力有限，為了進一步提升計算能力，基于多機互聯的分布式計算應運而生：大量由通信網絡連接的計算機構成一張網絡，通過任務分配調用網絡中閑置的計算資源，實現巨大算力的使用，這便是網格計算；而由特殊的高速通信網絡和同構計算機互聯形成的則是計算集群如常見的超級計算機。隨著虛擬化技術的發展，網格計算進一步發展為如今的云計算，“計算”作為一種資源以越來越便捷的方式提供給客戶。

2.2 有限元方法加速技術

對于非線性問題和時域問題，通過線性化、時域離散后都可化簡為若干穩態問題的求解。使用有限元方法求解穩態問題可分為3步：單元矩陣生成、整體矩陣裝配以及線性方程組求解。在實際計算中這3步都可能成為瓶頸。如圖3所示，對于非線性色散介質中電場仿真，隨著計算規模增長，單元矩陣分析和裝配耗時占總計算時長的99%［23］；而在電磁輻射場仿真中，表4結果指出線性方程組的求解耗時在總計算耗時中占主導地位［20］。下面，將分別介紹單元矩陣生成、裝配和線性方程組求解這3步的加速技術。

圖3 不同計算規模下整體剛性矩陣生成耗時在總計算時間中占比Fig.3 Proportion of the overall rigid matrix generation time in different calculation scales in the total calculation time

表4 矩陣裝配和求解環節的計算時間對比Tab.4 Comparison of calculation time between matrix assembly and solution

2.2.1 單元矩陣生成

單元分析時僅用到各個單元自身信息，可以進行逐單元并行：對于低階單元，往往使用單個線程（或單個核心）來生成對應的單元剛性矩陣［24］；對于高階單元，單元剛性矩陣規模較大，需要分配更多計算資源來生成一個單元矩陣［25］。由于不同單元在分析過程中都需要讀入相關單元的節點信息，大量的內存讀寫操作使得這類傳統的并行方法效率提升有限。另一種加速思路是將單元分析過程進行張量化改寫。

2.2.2 整體矩陣裝配

獲得各個單元的局部剛性矩陣后，需要裝配形成整體剛性矩陣。由于一個節點的所有鄰接單元矩陣對這一結點均有貢獻，在裝配時的歸并操作將導致內存的沖突讀寫。為解決這一問題，一種常見的解決思路是染色法：對計算區域內單元進行染色，使得相鄰單元不同色。每次將具有相同顏色單元的局部剛性矩陣裝配至整體剛性矩陣中，從而避免了競爭沖突；另一種解決方法是將矩陣裝配環節和后續的線性方程組求解環節相結合，如逐單元（element by element）有限元［26］、節點分解（nodal domain decomposition）有限元［27］等，這些免裝配方法實際上都可視為多波前方法的變種［28］。

2.2.3 稀疏線性方程組求解

通過單元矩陣生成和裝配后，原偏微分方程被離散成大型稀疏線性方程組。如何求解稀疏線性代數方程組一直是高性能計算領域的核心議題。針對大型稀疏線性方程組的解法可以分為兩類。

第一類是基于高斯消去法的直接求解器。直接計算方法可分為填充元優化、符號分解、數值分解、消元回代4步。其中填充元優化是對矩陣進行重排序，在保證數值穩定性的同時，使得后續矩陣分解時產生盡可能少的填充元，從而提高計算效率。常見的方法如最小度算法、重復剖分算法以及多重剖分算法。而符號分解則是現代快速直接解法特有的技術，通過使用消去樹這一數據結構，確定分解因子矩陣的非零元位置，指導分解路徑，如基于一般消去樹的多波前方法和基于超節點消去樹的多元最小度方法等。直接解法的優點是計算穩定可靠，因此也成為大量通用有限元計算軟件中的默認求解器；其缺點是計算過程中產生的分解因子、波前矩陣等中間變量規模不定，隨著問題規模增大，超出內存容量后超出部分只能存儲在硬盤中，硬盤的I/O操作成為計算瓶頸制約了直接解法的性能。

第二類是迭代解法。預計于Gauss消去法的直接解法不同，迭代解法是在特定的解空間中不斷尋找更接近真解的元素，其中最常見的是基于krylov子空間的迭代方法：如先驗共軛梯度迭代算法（preconditioned conjugate gradient，PCG）方法、雙共軛梯度（biconjugate gradient，BiCG）方法等。在迭代求解過程中需要不斷進行稀疏矩陣向量乘法操作（sparse matrix-vector multiplication，SpMV）來獲得當前解對應的殘差，并基于此進行修正。如何高效地執行SpMV操作成為該類解法的研究重點。研究人員通過調整稀疏矩陣的存儲格式來實現內存的合并讀寫，優化SpMV操作的性能。常見的存儲格式可分為兩類：行主元的COO格式、CSR格式、CTO格式以及列主元的ELL格式，Sliced ELL格式等［29］。針對硬件配置選擇特定的存儲格式，在存儲空間最小化同時通過緩存的合并讀寫實現時空效率的最大化。

迭代求解器中另一個關鍵技術是預條件處理。隨著問題規模增大，剛性矩陣的條件數越來越大，通過預條件處理能夠顯著提升迭代算法的收斂性能。在并行求解器中，剛性矩陣往往是分塊存儲，如何根據局部矩陣信息生成預條件矩陣是并行求解器研究的熱點，目前發展的并行預條件矩陣有：不完全LU分解預條件子，分塊Jacobi預條件子，最小二乘多項式預條件子以及多重網格預條件子等。相較于直接求解器，迭代求解器計算過程中對存儲的需求固定，借助區域分解方法容易實現稀疏矩陣-向量乘法操作的并行化，適合求解較大規模問題；而迭代解法的缺點也非常明顯：無法保證收斂、生成性態優良的預條件矩陣較為耗時。

縱觀有限元發展歷史，求解器經歷了迭代方法（Gauss-Seidel方法）、直接方法（變帶寬方法）、PCG方法再到多波前等直接解法的螺旋式發展［30］。如今直接解法能夠有效求解百萬自由度的有限元方程，而迭代解法多用于更大規模問題［31］。無論是直接解法還是迭代解法，性能的提升都是通過發掘稀疏矩陣自身結構并與計算機硬件相結合來實現，隨著計算機硬件發展，如何利用多核和眾核處理器實現高效求解成為新的問題。

3 高性能有限元技術發展現狀與展望

3.1 從工業軟件角度重新認識電磁仿真軟件

隨著物聯網（internet of things，IoT）技術的成熟，實時獲取設備運行過程中的各種參數成為可能。為了有效利用這些參數對設備全生命周期進行預測和控制，數字孿生（digital twin）這一概念被提出。由于實際設備的物理模型過于復雜，傳統數值仿真往往采用簡化模型：如使用集總參數模型、將非線性時變的材料特性線性化等。仿真模型和實際模型間存在較大的差異，極大地限制了仿真結果的準確性，因此傳統數值仿真僅作為優化設計的輔助工具。大量設備內部參數的出現使得數值仿真模型能夠深入設備內部，盡可能真實地還原實際運行工況，在此基礎上得到的仿真計算結果準確度大大提高，基于此能夠進行對設備運行過程中出現的故障進行檢測、定位以及預測。

數字孿生概念的出現極大地拓寬了電磁場數值仿真的應用場景，同時也提出新的挑戰：精細的模型往往兼具時域、三維、非線性等特征，某些復雜的問題還涉及熱場、流場、應力場等多物理場耦合。如何利用邊緣計算、云計算等新興技術實現穩定高效地求解是一個亟待解決的問題。

3.2 借力國際開放資源

在國際科學計算社區中不同研究團隊分別提出并發布了用于偏微分方程求解的高性能數值計算庫：Gmsh提供二維、三維結構的網格剖分；METIS則能夠對網格編號進行重排、分解；針對大型稀疏代數方程組，MUMPS提供了基于多波前方法的直接求解器，Hypre、AGMG等軟件則提供了基于代數多重網格方法的迭代求解器；可移植可拓展科學計算工具箱PETSc專門用于求解偏微分方程；Trillinos使用面向對象的軟件框架，通過調用各個組件實現大規模復雜多物理場求解；另外還存在大量基于上述計算組件二次開發的有限元計算軟件，如Deal.Ⅱ、FENICs、MFEM等。這些開源的高性能計算軟件是當前偏微分方程數值解研究的重要基礎，在計算流體力學、計算力學等領域都得到廣泛的應用，但在計算電磁學領域未見報道。

3.3 結合人工智能技術的數值計算方法

近年來，大數據和人工智能技術突飛猛進，深度神經網絡技術給傳統科學計算帶來新的詮釋。神經計算框架之下，科學計算問題可分為3類［32］：第1類問題是完全由數據驅動的問題，常見的如對風、光等自然資源的實時預測。由于實際模型未知，往往直接簡化為黑箱模型進行處理，采用神經網絡對數據進行擬合、挖掘；第2類問題則是基于部分物理機理和少量測量數據的不確定性模型，常見的如逆散射問題，即控制方程已知，根據測量數據反演出待求的模型參數。這一類問題往往可以寫為PDE約束下的優化問題，借助神經網絡這一優化求解器進行求解；第3類問題是基于已知模型的確定性問題，經典計算電磁學問題往往屬于這一類。有限元過程可以等價改寫為神經網絡的形式（FEA-net），即每個神經元的輸入為單元幾何信息，輸出為單元剛性矩陣［33］。應用神經網絡技術求解上述3類問題具有兩個顯著優勢：一是借助各類演化算法能夠處理復雜的非線性優化問題；二是神經網絡本身是對計算任務的細粒度分解，能夠實現高性能計算。將經典的計算電磁學問題與如今蓬勃發展的神經計算相結合，也是未來的發展方向之一。

3.4 發展異構計算和跨平臺的高性能計算方法

高性能計算技術的發展與高性能計算機的發展是密切相關的，基于高性能計算技術的有限元方法在計算電磁學中的應用經歷了從最開始只有單核CPU，發展到多核處理器，再到多核處理器和眾核處理器并存的異構計算環境。隨著硬件平臺在不斷更新換代，單位計算資源的成本越來越低：1997年，提供1TFlops計算能力的硬件需要花費4 600萬美元，而在2017年，使用相同的計算資源僅需要30美元［34］。相較于發展精巧的計算方法，高效利用“廉價”的硬件資源求解復雜的電磁學問題更為重要，也更加有效。

在目前高性能有限元計算方案中，還存在以下問題：其一，現有研究針對特定硬件的加速技術，這些技術往往受限于特定的問題或特定的計算平臺，尚缺乏一款具有可移植性、可拓展性以及較高硬件利用率的計算程序；其二，除多核與眾核處理器外，張量處理器 TPU［35］、神經計算單元 NPU［36］、量子處理器［37］等新興的計算單元也在不斷發展，有限元方法與這些硬件的結合也是一個亟待解決的問題。

4 結語

有限元方法作為一種通用的數值求解技術，以其強大的計算能力，在電磁仿真領域獲得巨大成功。由于新應用場景的出現，電力系統高壓設備電磁仿真對有限元方法提出了新的需求。從技術路徑發展上，高性能有限元技術是上述問題的一個可能的解決方案。本文介紹了高性能有限元方法的主要技術要點，討論其相關的發展途徑，希望對高性能有限元方法在電力系統高壓設備仿真的應用有所推動。