基于等效磁荷法用蒙特卡洛法计算永磁轴承磁力

2013-07-21张钢张坚张海龙孟庆涛樊曼

轴承 2013年10期

张钢，张坚，张海龙，孟庆涛，樊曼

(上海大学机电工程与自动化学院，上海 200072)

永磁悬浮轴承(以下简称永磁轴承)是利用永磁体在气隙中产生的磁场实现转轴稳定悬浮于空间的一种新型高技术轴承。由于其具有寿命长、无磨损、无需润滑等显著优点，因此得到了广泛关注和研究[1-7]。承载能力和刚度计算是永磁轴承研究与设计过程中的重要环节。而永磁轴承的磁力解析计算，目前常用的方法是利用等效磁荷法或分子电流法建立磁力模型，通过求解一个复杂的四重积分得到[4-7]。由于四重积分求解比较复杂，有学者直接用Matlab中的符号积分函数求解，不仅计算周期长，而且往往得不到正确的数值解；也有学者用C++进行编程求解，但是编程过程较为复杂[7]。

基于以上原因，提出应用蒙特卡洛法对基于等效磁荷法建立的永磁轴承磁力模型的四重积分方程组进行求解；并与ANSYS分析得到的磁力值进行比对。

1 永磁轴承磁力模型

永磁轴承的基本结构如图1所示，其中A为动磁环，B为静磁环，c0为两磁环间的初始气隙。当与转子固定的动磁环受到z轴正方向外力Fz作用时，将在沿z轴正方向上产生一轴向位移z(0

图1 永磁轴承结构图

图2 动磁环A受力示意图

依据等效磁荷理论[1]，可建立面1和面3之间的磁力F13,z的解析式

dr1dr3dαdβ/{[(c0-z+2L)2+(r3cosβ-r1cosα)2+

(r3sinβ-r1sinα)2]3/2}；

(1)

同理可建立面2和面4、面1和面4以及面2和面3的磁力F24,z，F14,z，F23,z的解析式

dαdβ/{[(c0-z)2+ (r4cosβ-r2cosα)2+

(r4sinβ-r2sinα)2]3/2}，

(2)

dαdβ/{[(c0-z+L)2+ (r4cosβ-r1cosα)2+(r4sinβ-r1sinα)2]3/2}，

(3)

dr2dr3dαdβ/{[(c0-z+L)2+(r3cosβ-r2cosα)2+ (r3sinβ-r2sinα)2]3/2}，

(4)

式中：Br为磁环剩磁(假设两磁环剩磁相同)；μ0为真空磁导率；(r1,α)，(r2,α)，(r3,β),(r4,β)分别为面1，2，3，4上任一点的极坐标。

规定磁力正方向为坐标轴正方向，则轴向外载荷Fz为

Fz=F13,z+F24,z-F14,z-F23,z。

(5)

永磁轴承的刚度是指动磁环在受到外扰力时，抵抗偏离平衡位置的能力。因此可以定义永磁轴承轴向刚度Kz为

(6)

2 蒙特卡洛法对磁力的求解

蒙特卡洛法求解多重积分的思想为任意一个积分都可看作是某个随机变量的数学期望[8]。因此，在利用蒙特卡洛法计算多重积分时，采用随机变量的算术平均值作为其近似值。

2.1 基本蒙特卡洛积分法

设D为n维空间Rn的一个区域，f(x)∈D⊂Rn→R，区域D上的n重积分用下式表示

(7)

I可以被认为等于区域D的测度乘以函数f(p)的期望。计算时首先要选择概率密度函数g(p1,p2,…，pn)，并且在区域D内满足g(p1,p2,…，pn)>0。

则(7)式可改写为

dp1dp2…dpn。

(8)

无论区域D的形状如何，根据概率密度函数g(p1,p2,…,pn)进行抽样时(抽样次数为N)，均可以给出积分的估计为

(9)

标准差为

(10)

(11)

积分误差为

(12)

式中：Xα为与置信度1-α有关的参数，当α已知，可以查到相应的Xα；σ为标准差。

基本的蒙特卡洛法就是找一个包含区域D的超立方体(测度已知，为Mc)，在D内随机生成N(N一般足够大)个均匀分布的点，统计落入区域D的点，假设有m个，则区域D的测度为

(13)

函数f的期望为

(14)

从而有

(15)

将(15)式应用到四重积分，可得

(16)

式中：bj为四重积分的4个积分上限；aj为对应的4个积分下限；j=1,2,3,4。

2.2 磁力的求解流程

如图3所示，按照图示流程即可将(16)式用程序实现，求解磁力非常方便。

图3 蒙特卡洛法计算磁力流程图

采用Matlab进行编程，生成超立方体内的随机数可以调用函数unifrnd，调用格式为[9]：unifrnd(积分下限，积分上限，1，N)。

值得注意的是，蒙特卡洛法是一种近似方法，是概率统计的方法，定义不同的N值，积分值不会相同。下面将通过定义不同的N值和多次求解磁力来探究N值对结果精度以及可重复性的影响。

3 永磁轴承磁力计算分析

对照图1的轴承结构，文中计算所采用的尺寸为R1=10 mm，R2=20 mm，R3=10 mm，R4=20 mm，L=15 mm，初始轴向气隙c0=6 mm。定义永磁环性能参数：Br=1.231 T，Hc=917 530 A/m，则磁环的相对磁导率μr为

(17)

将上述参数代入到等效磁荷法公式中，编制程序进行求解。

3.1 样本数量N对磁力值的影响

其他参数均不变，改变N值得出的永磁轴承磁力曲线如图4、图5所示。

图4 改变N值时磁力曲线变化示意图

图5 不同N值时z=4 mm处的磁力曲线(重复运行程序6次)

从图4可以看出，随着N值的减小，磁力曲线波动越来越大，曲线越来越不平滑，计算数值准确性随之大大降低；并且积分区域越大，蒙特卡洛法的计算精度越低，因此，当磁环面积较大时，应该提高抽样数量N。

从图5可以看出，随着N值的减小，数据的可重复性差；蒙特卡洛积分误差与抽样数的平方根成反比，在满足工程应用(误差在10%以内)的前提下，N应取一合适值(106)。

3.2 蒙特卡洛法与数值解法对比

研究表明，在计算永磁轴承的承载能力时，有限元数值解是较为准确的[2-3，10]。用ANSYS求解时，若只考虑轴向气隙的变化，采用轴对称二维模型，单元类型选择plane53，磁环网格尺寸为0.1 mm，空气网格尺寸为2 mm，给最外层空气边界施加通量平行条件(即Az=0)，将动磁环定义为一个组件，求解即可得到磁力线分布和轴向力的值。

图6和图7分别为蒙特卡洛法(N=1 000 000)与有限元数值解法的轴向力和轴向刚度计算结果对比。