薄殼單元法在接觸物體間電磁力計算中的應用

劉慧娟， 傅為農

(1.北京交通大學電氣工程學院，北京 100044;2.香港理工大學電氣工程學院，香港特別行政區)

薄殼單元法在接觸物體間電磁力計算中的應用

劉慧娟1， 傅為農2

(1.北京交通大學電氣工程學院，北京 100044;2.香港理工大學電氣工程學院，香港特別行政區)

針對在利用有限元法計算接觸物體間電磁力時在物體接觸表面沒有積分單元可用的問題，首先，基于薄殼單元的思想，提出了一種在體積分和面積分之間確定附加薄殼單元的方法和一種在3D有限元法中利用附加薄殼單元精確計算接觸物體間電磁力的方法，然后，基于虛功原理，推導了計算薄殼單元節點電磁力的T－Ω公式和雅可比矩陣表達式，并給出了計算物體接觸表面分布電磁力的算法及其實現方法，最后將所提出的方法用于實例計算，并將3D有限元法計算結果與解析計算結果進行比較，證實所提出方法的有效性和正確性。

薄殼單元法;3D有限元法;接觸面電磁力;虛功;T－Ω公式;雅可比矩陣;電磁設備;解析計算

0 引言

在電磁設備的設計與性能分析中，常采用有限元法(finite element method，FEM)，為了能更精確反映電磁設備內容的電磁過程，通常需要采用3D FEM進行瞬態電磁場或渦流場的計算［1］。3D有限元公式有兩類:一類是采用矢量磁位(magnetic vector potential)A和標量電位(electric scalar potential)V的A－V公式，另一類是采用矢量電位T和標量磁位Ω的T－Ω公式，在低頻電磁設備的設計與分析中，通常采用 T － Ω 公式［2－3］。

在低頻電磁設備的有限元計算中，通常需要進行電磁力的計算［4－6］，而利用FEM計算某物體所受電磁力時，是在FEM的后處理中，通過沿該物體表面的一層單元進行體積分而計算得到的［7－9］。但是，在許多電磁設備中，常存在物體間的相互接觸，此時，若需要利用FEM進行接觸物體間電磁力的計算，就會出現在物體接觸表面沒有積分單元可用的問題，因此，利用通常的FEM后處理方法就不能計算接觸物體間的電磁力。為此，必須在兩物體的接觸面構造出新的積分單元，才能進行電磁力計算。文獻［10］首次提出在接觸物體的表面添加附加薄殼單元來計算電磁力的思想，卻并沒有給出具體3D FEM的實現算法。

本文基于薄殼單元的思想，提出了一種在體積分和面積分之間確定附加薄殼單元的方法，然后，基于虛功原理，推導了利用局部磁能偏導數來計算3D有限元電磁力的T－Ω公式和雅可比矩陣表達式，并給出了計算物體接觸表面分布電磁力的算法及其實現方法，該方法可用于計算接觸表面局部和全局的電磁力，特別適合于永磁體(permanent magnet，PM)材料與其它物體接觸的表面。最后將該方法用于計算兩個PM物體接觸面的電磁力值，并將FEM計算結果與文獻［11］的解析計算結果進行了比較，證實了本文方法的有效性和正確性。

1 節點電磁力計算公式

在每個單元中，利用與位移s相關的磁共能的導數來計算各節點電磁力的公式為

其中:B是磁密，H是磁場強度，Hc是永磁體(PM)的矯頑力，Re是自然坐標中的單元面積;R^e是參考坐標中的單元面積，J表示坐標變換之間的雅可比矩陣。

其中Hp是場源，為固體導體或多匝線圈中的已知全部電流或電流密度。

可將磁場效應Heff表示為

其中Nm是具有10個節點的二階節點單元的形狀函數。于是Ω為

用一種新穎的四面體單元來計算磁標量Ω，即

在每個邊界單元上的磁場強度H為

其中τm是具有6條邊的邊界單元的形狀函數。

將式(2)、式(3)代入式(1)，有

全局坐標與局部坐標之間的關系為

其中ζm是坐標的形狀函數。雅可比矩陣為

其中|J|=6V，V 是單元的體積，xij=xi－ xj，yij=yiyj，zij=zi－ zj，且 dxdydz=|J(x，y，z)|dζ0dζ1dζ2。假設節點i沿方向的位移為

2 接觸物體間薄殼單元的處理

當兩個物體相互直接接觸，電磁力計算就需要沿接觸表面進行面積分。為此，在接觸面添加如圖1(a)和圖1(b)所示的三棱柱薄殼單元，且薄殼單元的高δ接近于0，薄殼單元的材料設為線性材料，則其磁導率等于真空磁導率，即μ=μ0。

磁場強度可表示為

于是薄殼單元中的磁共能為

圖1 三棱柱薄殼單元Fig.1 A triangular prism shell element

據文獻［10］，將式(19)對位移s求偏導數，可得到單個薄殼單元的節點力為

式(20)表明:通過保持邊界上Hi的恒定，即電流恒定，可實現對電磁力的計算。

在式(21)中，只有雅克比矩陣與位移相關，參考單元的函數和面積均與位移無關，因此有

從自然坐標系(x，y，z)到參考坐標系(ξ，η，ζ)變換的雅可比矩陣為

式(24)和式(25)中五面體單元的J和|J|可根據式(23)的基本雅克比矩陣推導得到。

3 算法的具體實現

在算法的實現過程中，需要在垂直方向設定兩個標志位nf和nt，參考圖2，將其定義為:nf=1表明該節點在需要計算電磁力的物體上，nf=0表明該節點不在需要計算電磁力的物體上。nt=－1表明該節點在需要計算電磁力物體的外表面，nt=1表明該節點在接觸面上，nt=0表明該節點不在需要計算電磁力的物體上。

對每個單元定義

當ef＞0，且至少有一個節點的nt=－1時，則需要對該單元進行體積分。當et=3時，則需要對該單元進行面積分。

圖2 各節點的標志位nt值示例Fig.2 An example of the flag nt

若將體積分改為面積分，就需要附加一個面積分單元。為簡單起見，利用圖2所示的3D示例來進行說明，并將薄殼單元按比例放大為如圖3所示。當薄殼單元的高度δ不等于0時，則需要進行電磁力計算的物體將被一層體單元包圍，薄殼單元3也包含在其中。于是，當薄殼單元的高度δ接近于0時，在體單元1和薄殼單元2之間就需要有另一個附加的薄殼單元3(如圖3所示)。

圖3 附加表面單元的示例Fig.3 An example of an additional surface element

在3D問題中，如圖4所示，當四面體在體積分和面積分之間的分界面上(圖4中每個四面體頂點上標注的數字為nt標志位值)時，存在3種情況，且只有情況1和情況2需要面積分。根據4個節點的標志位nt的數值，將很容易確定需要在哪個面上進行面積分。

圖4 四面體位于體積分和面積分之間的3種情況Fig.4 Three cases when the tetrahedron is in the interface between volume integration and surface integration

4 應用實例

為驗證本文方法的可行性，將該方法用于實例計算，為進一步驗證該方法的正確性，選用文獻［11］中的兩個相互接觸的矩形永磁體的實例模型，并將本文的3D FEM計算結果與該文獻中的解析計算結果進行比較。計算實例如圖5所示，兩永磁體沿z軸放置，具有相同的x和y坐標，其尺寸均為40 mm ×20 mm ×10 mm，且 Br=1.1 T，μ =μ0，當兩永磁體在z軸方向的距離為0時，表明兩個物體相互接觸。

表1為兩永磁體之間距離變化時，利用3D FEM計算的電磁力和文獻［11］中解析計算的電磁力的數值，其中FEM的剖分單元數為83 666。從表1可見，有限元數值計算結果與解析計算結果有很好的一致性，特別是當兩永磁體相互接觸時，有限元計算結果與解析計算結果的差為0.041N，誤差為0.022%，表明利用本文的方法能正確計算兩永磁體接觸面之間的電磁力。

圖5 兩永磁體模型Fig.5 The model of two PM objects

表2為兩物體相互接觸時，采用不同的有限元剖分方法，計算的全部電磁力和附加薄殼單元電磁力的值。可見，當剖分網格較大時，附加薄殼單元的引入非常重要，并且，隨著有限元剖分單元數量的增加，附加薄殼單元的電磁力逐漸減小，這是因為剖分單元數增多時，每個薄殼單元的尺寸都逐漸變小，因而對積分的貢獻也就逐漸變小。

表1 3D FEM與解析計算結果的比較Table 1 Comparison of forces between two methods

表2 不同剖分情況全部電磁力和附加薄殼單元電磁力值Table 2 The total force and the force of additional shell elements vs number of FEM elements

5 結論

本文提出了一種在3D FEM中利用附加薄殼單元計算接觸物體間電磁力的方法。通過對應用實例的計算表明:

1)在3D FEM計算的后處理中，在體積分和面積分之間添加薄殼單元積分，可以精確計算接觸物體間附加薄殼單元的節點力，其計算結果的精度與解析計算結果精度一致。

2)在利用3D FEM計算接觸物體間電磁力時，通過附加薄殼單元，既可實現對接觸物體全部電磁力的精確計算，也可實現對其局部電磁力的精確計算。并且，FEM網格剖分的數量和質量直接影響計算結果的精確度。

3)所提出的利用薄殼單元法計算接觸物體間電磁力的T－Ω公式適用于其它材料類型的3D接觸物體模型。當然，也可推導相應的3D A－V公式應用于FEM的后處理，進行接觸物體間電磁力的計算。

［1］ STANCHEVA R D，IATCHEVA I I.3-D electromagnetic force distribution in the end region of turbogenerator［J］.IEEE Transactions on Magnetics，2009，45(3):1000 －1003.

［2］ ZHOU P，FU W N，LIN D，et al.Numerical modeling of magnetic devices［J］.IEEE Transactions on Magnetics，2004，40(4):1803－1809.

［3］ ZHOU P，LIN D，FU W N，et al.A general cosimulation approach for coupled field-circuit problems［J］.IEEE Transactions on Magnetics，2006，42(4):1051 －1054.

［4］ ZHU Z Q，WU L J，XIA Z P.An accurate subdomain model for magnetic field computation in slotted surface-mounted permanentmagnet machines［J］.IEEE Transactions on Magnetics，2010，46(4):1100－1115.

［5］ JIANG Hao，HUANG Xueliang，ZHOU Gan，et al.Analytical force calculations for high-precision planar actuator with Halbach magnet array［J］.IEEE Transactions on Magnetics，2009，45(10):4543－4546.

［6］ RAVAUD R，LEMARQUAND G，LEMARQUAND V.Force and stiffness of passive magnetic bearings using permanent magnets.part 1:axial magnetization［J］.IEEE Transactions on Magnetics，2009，45(7):2996 －3002.

［7］ COULOMB J.A methodology for the determination of global electromechanical quantities from a finite element analysis and its application to the evaluation of magnetic forces，torques and stiffness［J］.IEEE Transactions on Magnetics，1983，19(6):2514－2519.

［8］ REN Z，RAZEK A.Local force computation in deformable bodies using edge elements［J］，IEEE Transactions on Magnetics，1992，28(2):1212－1215.

［9］ FU W N，HO S L.Error estimation for the computation of force using the virtual work method on finite element models［J］.IEEE Transactions on Magnetics，2009，45(3):1388 －1391.

［10］ REN Z，CENDES Z.Shell elements for the computation of magnetic forces［J］.IEEE Transactions on Magnetics，2001，37(5):3171－3174.

［11］ AKOUN G，YONNET J-P.3D analytical calculation of the forces exerted between two cuboidal magnets［J］.IEEE Transactions on Magnetics，1984，20(5):1962 －1964.

(編輯:張詩閣)

Computation of electromagnetic force on the interface using shell element method

LIU Hui-juan1，FU Wei-nong2
(1.School of Electrical Engineering，Beijing Jiaotong University，Beijing 100044，China;
2.Department of Electrical Engineering，The Hong Kong Polytechnic University，Hong Kong，China)

In order to compute the electromagnetic force between two objects touching each other，an additional shell element is proposed and a simple method to determine the additional shell elements on the interface between volume integration and surface integration is also presented.A general formulation of nodal virtual work electromagnetic force computation between two objects touching each other for three-dimensional(3D)finite element analysis was presented.The T-Ω formulation and the Jacobian matrix were deduced and the implementation methodology was presented.An example with analytical solutions is used to verify the computed results.

shell element method;three-dimensional finite element method;interface electromagnetic force;virtual work;T-Ω formulation;Jacobian matrix;electromagnetic devices;analytical computation

O 242.21;TM 351

A

1007－449X(2012)08－0101－06

2012－04－01

香港理工大學科研項目(GU489)

劉慧娟(1967—)，女，博士，副教授，研究方向為電機及其控制、電機電磁場計算等;

傅為農(1962—)，男，博士，副教授，研究方向為電磁場數值計算的理論與應用、電機及特種電機的優化設計。

劉慧娟