方 錦,閻秀恪,鐘立國,任自艷,張殿海
(沈陽工業大學電氣工程學院,遼寧 沈陽 110870)
在使用傳統有限元法對大型電氣設備進行電磁場數值分析時,為了提高計算精度,需要對結構模型進行較為精細的剖分,剖分產生的大量網格使得計算規模龐大,需要大量的計算資源和計算時間,有時甚至無法求解。為了提升有限元計算效率,學者們在網格剖分、方程組求解、并行計算等方面做了很多研究。
在有限元計算中非線性問題的求解會導致一些耗時問題。迭代法是非線性問題求解常用的方法,如直接迭代法、增量法、牛頓-拉夫遜方法等,其中牛頓-拉夫遜算法常常被使用,并不斷得到改進。牛頓-拉夫遜方法在求解時需要通過修改全局Jacobi矩陣進行迭代更新。隨著計算規模增大,全局Jacobi矩陣也變得十分巨大,修改全局Jacobi矩陣消耗大量的時間。
最初,傳輸線方法(transmission line model,TLM)被應用于研究波的傳播理論,波的傳播過程中入射和反射過程類似于非線性迭代過程,因此也被研究用于求解非線性問題。TLM法中電路網絡節點電壓方程中的導納矩陣與FEM的系數矩陣具有相同的性質,根據對偶原理建立與有限元網絡相對應的等效電路網絡,因此用于求解電路網絡的TLM方法,也可以用于求解有限元網絡。TLM方法與有限元法相結合從而形成傳輸線有限元法(transmission line model-finite element method,TLM-FEM)。由于傳輸線的入射電壓與反射電壓之間存在時間延遲,混合電路網絡可以解耦為線性部分和非線性部分。應用TLM方法的解耦特性,可將有限元網絡轉化為1個線性系統和若干個獨立的非線性部分。
TLM-FEM起初應用于求解非線性的靜磁場問題[1-8]。目前已經應用在電氣設備電磁場和溫度場的分析中[9-15],文獻[16]在使用傳輸線有限元等效模型求解軸對稱磁場問題時,在入射階段采用并行求解器加快單步迭代計算速度。文獻[17]結合傳輸線方法提出了一種黑盒電路模型提高了復雜問題解決能力。目前對于TLM-FEM的研究在入射階段的線性方程組求解均使用LU分解法,未能充分利用全局矩陣對稱稀疏的特點加速分解過程。同時反射階段非線性負載兩端電壓的計算一般使用牛頓迭代法,迭代過程也消耗了很多時間。因此在入射階段線性方程組的高效求解與反射階段計算方法的研究可以提升整體求解效率。
本文推導了TLM-FEM的數學模型,將分段線性化方法引入到求解反射電壓的松弛方法中,用于更新反射階段的參數,以提高求解的效率和精度。根據單元系數矩陣計算及全局矩陣LDLT分解計算的特點,在GPU平臺上實現了2個過程的并行計算。數學模型和算法均由C++語言自編程實現。將上述算法和程序應用于單相變壓器的磁場計算,并與ANSYS仿真結果進行了比較和分析。
傳統傳輸線迭代方法包含入射階段和反射階段。在含有非線性元件的電路中,無損傳輸線將線性系統和非線性元件連接在一起,根據傳輸線理論由于傳輸線特性導納與終端導納間存在不匹配現象,導致電信號在傳輸線上的傳播會在兩端不斷反射。只有當傳輸線兩端電壓相等時,電路系統達到穩定狀態能求出負載導納兩端的電壓值,從而實現電路問題的求解。圖1為非線性電路中第e個非線性負載的TLM方法的等效電路。

圖1 非線性負載的TLM方法

TLM法中電路網絡節點電壓方程的導納矩陣與有限元方法的系數矩陣具有相同性質,以三角形單元為例,根據對偶原理,有限元單元系數矩陣對偶于電路單元導納矩陣,節點磁位對偶于節點電壓,右端激勵項對偶于電路的電流源,因此可以得到如圖2所示的傳輸線有限元的等效模型。

圖2 傳輸線有限元等效模型的建立過程

傳輸線有限元等效過程主要分為3個步驟:一是先將在有限元三角形單元的邊緣分別添加負載導納,得到等效電路拓撲網絡;二是在單元邊緣與負載導納之間添加1組傳輸線;三是將傳輸線有限元等效電路解耦得到入射和反射等效電路,最終建立傳輸線有限元等效模型。
當等效電路中有激勵產生時,首先電壓信號從各個單元的非線性負載入射進線性系統。此時傳輸線迭代方法的入射階段可以被描述為
YV=I+Ii
(1)
式中:Y為入射階段的全局系數矩陣;I為電流源激勵向量;Ii為入射階段等效電流源向量。
當傳輸線一端的負載導納與傳輸線特性導納之間存在不匹配現象,電壓信號將會從線性系統反射回各個單元的負載導納。此時TLM方法的反射階段在各個負載導納上獨立進行,以端口ij的反射階段的計算為例。
(2)

(3)
(4)
負載導納與對應特性導納間的不匹配,可以通過式(3)、式(4)修正。
后續迭代重復上述過程,只有當傳輸線兩端電壓差達到收斂誤差則傳輸線迭代過程停止,同時輸出式(1)中求解出的每個節點電位。

圖3 結合優化松弛方法的TLM迭代收斂過程


(5)
(6)

(7)

(8)
(9)
綜上所述,在進行TLM-FEM迭代時,反射階段的求解可以通過直接計算式(9)得到,計算結果用于更新式(1)中等效電流源。結合優化松弛方法的TLM-FEM求解流程的偽代碼如圖4所示。

圖4 結合優化松弛方法的TLM-FEM程序

在近年來,GPU由于其出色的并行計算能力被經常應用于大規模的數值計算,GPU目前采用單指令多線程執行模型,可以將計算任務分配在多個線程上同時處理,非常適合執行邏輯簡單的大規模密集型計算任務。在TLM-FEM計算中,各個單元系數矩陣的計算是獨立進行的,計算過程適合采用GPU單指令多線程的計算模式。以三角形單元為例,單元系數矩陣元素包括每個傳輸線導納和線性部分導納,單元系數矩陣表示為
(10)
式中:b、c的值與三角形節點坐標i、j、m有關;Δ為三角形單元的面積;線性單元的磁阻率ve是固定值,非線性單元會提前輸入1個預設值。
由式(10)可知,每個單元的計算主要是利用有限元剖分信息執行乘法及加法操作,計算需要的內存少、邏輯相關弱,可設置1個或多個線程對每個單元執行計算操作。在CUDA平臺上利用GPU上進行單元系數矩陣的并行計算,其過程分為3個步驟:
a. cudaMalloc申請計算所需要的內存空間,并由cudaMemcpy將參與計算的數據從主機端復制至設備端;
b. 根據剖分網格單元數量分配線程塊的個數以及線程數量,調用內核函數在GPU上計算;
c. cudaMemcpy將計算結果從設備端復制回主機端,進行傳輸線有限元后續計算。
在TLM迭代中,全局矩陣Y被組裝以后,在后續迭代中始終保持不變。因此,使用直接法求解式(1)更加有效。通常,LU分解方法被廣泛應用于求解這一類線性方程。Y只在初次迭代中被分解,其分解結果可以在后續迭代中被重復使用。
由于全局矩陣Y具有對稱性,因此分解產生的矩陣U能進一步被分解為對角矩陣D和上三角矩陣LT。
Y=LU=LDLT
(11)
式(11)的矩陣形式為
(12)
在LDLT分解中,由于全局系數矩陣也具有稀疏性,因此計算范圍可以由所有元素縮小到非零元素,這樣可以減小計算規模和節省內存。
以矩陣Y第k列分解為例,全局矩陣的LDLT分解過程如圖5所示。

圖5 全局矩陣的分解過程
圖5中,將該列中所有的非零元素提取出來。提取出來的元素組成的向量A表示為
A=[a0a1…am]T
(13)
在Y中被提取出的非零元素對應位置信息由一個新的數組記錄。
row=[r0r1…rm]
(14)
每次提取均在同一列進行,因此只需要記錄行號,A中各元素的行號記錄在row中。
其主要分為2個過程:
a. 使用向量標量除法分解當前列中所有非零元素。
ai=ai/a0
(15)
b. 更新Y中那些對第k列有依賴的元素。
(16)
A與AT的乘積所指示的Y中相同位置的元素會被更新為
yrirj=yrirj-aiaja0(i,j (17) 由上述過程可知,全局矩陣中每1列元素的分解操作和更新操作是2個并行度較高的流程,適合在GPUs上并行執行,2個流程間采用串行方式。單元分析及全局矩陣LDLT分解的并行過程見圖6。 圖6 全局矩陣的LDLT并行分解過程 為了驗證優化后的TLM-FEM在電磁場計算中的求解精度和效率,本文將算法和程序應用于1臺型號為DSP-241000 kVA/500 kV的單相變壓器磁場計算,并將計算結果與采用商用軟件ANSYS計算得到的結果進行對比分析。變壓器的結構模型和網格剖分見圖7。變壓器一次側施加空載電流,二次側開路,計算得到的磁場分布見圖8。 圖7 變壓器幾何模型和網格剖分 圖8 變壓器磁通密度分布 在求解域內隨機選擇6個位置,提取ANSYS和本文算法的計算結果如表1所示。 表1 ANSYS與本文算法結果對比 如表1所示,本文算法誤差是在可以接受范圍,位置1處的誤差大于5%,分析原因是由于該處的值本身較小,其余各點的值均在3%以內,證明本文算法正確。 為考察本文開發的并行計算對計算時間的影響,對結合LU分解的TLM(TLM-LU)方法和本文算法在不同單元、節點數目下進行求解,計算時間見表2。本文算法程序在帶有NVIDIA 940MX GPU和Inter i7-6700HQ CPU的計算機上執行。 表2 TLM-LU和本文算法計算時間的對比 由表2可知,本文算法相較于TLM-LU方法計算時間減少,隨著網格數量增多、計算規模增大,加速效果更為明顯。 本文提出優化松弛法用于計算TLM-FEM反射階段的參數,以提高求解的效率和精度,同時將單元系數矩陣計算及全局矩陣LDLT分解計算在GPU平臺上并行實現。針對于單相變壓器電磁場問題的求解,該方法可以達到與商業軟件一致的求解精度。本文算法能通過縮短全局矩陣分解時間和反射過程計算時間從而提升整體求解效率,模型數量規模越大,加速效果越明顯。
3 算法驗證與應用




4 結語