邊光滑有限元-邊界元耦合法計算二維瞬態渦流場

王洋洋，蔣興良

(輸配電裝備及系統安全與新技術國家重點實驗室(重慶大學)，重慶 400044)

在瞬態渦流場研究中，常用的數值計算方法有積分法[1-2]、有限元法[3]、邊界元法以及有限元邊界元耦合法[4-7]. 雖然傳統有限元法(FEM)[8]作為最成功的數值方法之一，已被廣泛用于解決工程和科學的各種實際問題，但是因為其有固有的缺點，如剖分網格密度以及計算階數的復雜程度都對有限元計算精度影響非常大. 單元數量的增加，及其插值函數精度的提高，才能得到更精確的數值解. 特別是在某些復雜的設計優化實例中，有限元法需要重復多次進行網格剖分以及重要部分網格加密操作，這不僅增加了設計的難度，而且還降低了設計計算的效率. 因此國內外學者一直致力于研究如何克服傳統有限元法的缺點.

近年來，在努力開發先進數值方法時，文獻[9]將無網格中的應變光滑技術與傳統有限元法相結合，引入梯度光滑的操作概念，提出一系列廣義梯度光滑有限元法. 其包括用梯度光滑技術分別對構造的節點光滑域、邊光滑域以及面光滑域進行光滑操作而形成的點光滑有限元法(NS-FEM)[10]、邊光滑有限元法(ES-FEM)[11]以及面光滑有限元法(FS-FEM)[12]. 3種方法各有優缺點. 相對于傳統有限元法，NS-FEM在精度和效率方面有所提高. 但是NS-FEM在空間上是穩定的在時間上可能是不穩定的，并且由于“過軟”特性而不能直接應用于解決動態問題[13]. 然而，ES-FEM可以克服時間不穩定性[14]，由于本文研究的是時變渦流場，因此克服時間不穩定性對于時變渦流場分析尤其重要. 此外，對于固體力學問題，具有線性積分的ES-FEM比線性有限元法更加準確和有效[15]. 至此，ES-FEM近年來得到了廣泛的應用. 對于FS-FEM更適用于三維問題.

本文提出一種用ES-FEM與邊界元耦合法計算二維瞬態渦流場. 該耦合法同時具有ES-FEM計算精度高以及邊界元法占用計算機內存少等優點. 首先將問題域離散化為一組三角形單元，然后進一步構造與三角形單元邊相關的積分域，采用伽遼金加權余量法，推導出基于邊光滑有限元-邊界元耦合算法的離散公式，在導體區域內采用ES-FEM，非導體區域內采用邊界元法. 研究導體區域的磁場分布，并與測量結果和有限元-邊界元耦合法計算結果進行對比. 計算結果表明，在相同網格密度條件下，基于ES-FEM邊界元耦合法計算精度更高，該方法具有更好的應用前景.

1 電磁場的分析模型

電磁分析的問題實際上是在給定邊界條件下求解一組麥克斯韋方程組的問題. 如圖1所示，為開域渦流場求解區域示意圖. 整個電磁場求解域分為渦流區V1和非渦流區V2，Γ1是V1和V2的內部邊界，非渦流區V2包括無源電流的非渦流區和有源電流的非渦流區.V1區域內采用ES-FEM，V2區域內采用BEM, 在內部邊界上滿足的邊界條件為

(1)

(2)

ν

(3)

圖1 開域渦流場求解示意圖

2 邊光滑有限元-邊界元耦合法介紹

2.1 邊光滑有限元法

邊光滑有限元法是在求解有限元的基礎上對其求解域進一步光滑，形成多個光滑子單元，然后在新形成的光滑子單元內引入梯度光滑操作，梯度光滑操作是指在新形成的光滑域上對有限元求解的系數矩陣進行光滑平均運算，其數學表達式為

其中Ax為光滑域的面積. 圖2為邊光滑有限元法的計算流程圖.

圖2 計算流程圖

對于二維求解域，先用劃分軟件Hypermesh對其進行三角形劃分，其中生成N個節點和p條邊. 傳統有限元的離散矩陣[16]為

(4)

(5)

式中AL為光滑域的面積，在二維求解域單元內，AL的表達式可根據參考文獻[9]計算得到.

圖3 基于邊的光滑域

對式(7)進行格林公式換算，那么對光滑域內位函數的面積分可以轉化為光滑域邊界上的線積分，即

(6)

將式(6)代入式(5)中，可得

(7)

其中：

(8)

(9)

對光滑域ΩL采用高斯積分，式(12)的光滑域邊界積分形式為

(10)

式中：NA為電磁場求解域邊界L的個數，NB為高斯點的個數，δq為高斯點積分的權函數.

由式(10)可知，dij是由光滑域各個邊中點處的形函數求得，而各個中點處的形函數可通過場節點插值得到，如圖4所示，點1、2、3、4分別為場節點，A、B分別為中心節點，g1、g2、g3、g4分別為高斯積分點[17].

基于邊L的光滑有限元系統的局部光滑域的系數矩陣可以寫成：

(11)

對上述局部系數矩陣進行整體組裝，得到域內整體系數矩陣為

(12)

圖4 ES-FEM形函數分布

2.2 邊界元法

對于非渦流區應用格林公式，可得對應邊界積分方程[5]為

(13)

采用伽遼金余量法對式(13)離散，得

(14)

2.3 基于邊光滑有限元邊界元耦合算法

將式(4)和式(14)聯立，可得

(15)

對式(15)中的導數項，采用向后差分，即

(16)

(17)

(18)

3 計算算例與分析

3.1 脈沖線圈-鋁板表面磁場分布

脈沖線圈-鋁板的結構參數：矩形截面為0.6 mm×4.8 mm的銅絲纏繞成內徑、外徑分別為3.2、25.5 mm,匝數為30的圓柱形線圈[19]. 脈沖線圈流過的瞬時電流波形，如圖5所示. 在線圈正上方放置0.8 mm厚的鋁板，電導率δ=3.48×107s/m. 線圈與鋁板的間隙為2 mm. 采用本文所提的基于邊有限元邊界元耦合法求得鋁板表面的磁感應強度，其結果分別與測量值以及有限元-邊界元法的計算結果進行對比，如圖6～8所示.

圖5 脈沖電流隨時間變化

圖6 磁矢量位A線的分布云圖

圖7 切向磁感應強度隨時間的變化(r=10 mm)

由圖6可知，離線圈越近，磁矢量位A值越大，與實際情況一致，說明本文的計算方法合理. 為了驗證本文所提方法的計算精度，設相對誤差Bi為

(19)

式中：B0為測量值，B1是計算值，i=x、y分別代表磁場沿切向、垂直方向的磁場分量.

圖8 垂直磁感應強度隨時間的變化(r=7.6 mm)

由圖7、8可知，基于邊光滑有限元邊界元耦合法與測量結果一致，因此可用本文所提的方法計算渦流場. 在網格單元劃分密度相同的條件下，運行兩種方法效率對比見表1，雖然本文所提方法的運行時間稍長于有限元邊界元耦合法，但其計算精度遠遠高于有限元邊界元耦合法，另外，對于目前計算機的計算能力而言，本文通過改變模型的大小，對比兩種方法的運行時間. 由表2可知，兩種方法運行時間的差別可以忽略. 由此可知，在計算渦流場時基于邊光滑有限元邊界元法能大大提高計算精度.

表1 效率對比

表2 運行時間對比

綜上所述，在進行渦流場分析時，應優先考慮基于邊光滑有限元邊界元耦合法.

3.2 鋁板表面磁場分布的影響因素

3.2.1 電流密度變化的影響

脈沖電流-鋁板模型參數如上所述，在進行電流密度對其磁場影響分析時，保持模型的其他參數不變，改變流過脈沖線圈的電流密度Js從10×106A/m2到180×106A/m2變化，間隔為30×106A/m2. 則上述電流密度變化的情況下鋁板表面磁場分量分布如圖9所示.

圖9 磁場強度隨渦流密度的變化

如圖9所示，磁場強度隨渦流密度的增大而增大，這是因為渦流密度Js與系數矩陣F成正相關，即當Js增大時，F增大. 由式(20)可知,隨著F增大，所求磁矢量位A增大，因此磁場強度增大.

3.2.2 脈沖線圈-鋁板間隙變化的影響

脈沖電流-鋁板模型參數如上所述，在進行脈沖線圈-鋁板間隙對其磁場影響分析時，保持模型的其他參數不變，改變間隙h從2 mm到8 mm之間變化，間隔為1 mm. 則上述間隙變化的情況下鋁板表面磁場分量分布如圖10所示.

圖10 磁場強度隨間隙的變化

由圖10可知，隨著脈沖線圈與鋁板之間間隙h的增大，磁場強度逐漸減小. 這說明如果為了增大鋁板中的渦流作用，即產生大的脈沖力，那么脈沖線圈需盡可能近的接近鋁板.

4 結 論

1)提出了一種計算渦流場分布的基于邊光滑有限元邊界元耦合法，該方法即具有ES-FEM計算精度高和邊界元法占用內存少的優點，并詳細給出了該方法的推導過程. 同時，利用本文所提的基于邊光滑有限元邊界元耦合法計算分析脈沖線圈上方鋁板表面的磁場分布，與測量結果對比表明，該方法是正確合理的.

2)在單元網格密度劃分相同的條件下，將本文得到的計算結果與傳統有限元邊界元耦合法計算得到的結果對比，雖然本文所提方法的運行時間稍長于有限元邊界元耦合法，但其計算精度遠遠高于有限元邊界元耦合法，另外，對于目前計算機的計算能力而言，兩種方法運行時間的差別可以忽略.

3)在計算模型的其他參數保持不變，只改變渦流密度的情況下，磁場強度隨渦流密度的增大而增大且曲線增大率幾乎保持不變，同理，磁感應強度隨著脈沖線圈與鋁板之間間隙h的增大而減小，且其曲線降低率也幾乎保持不變.