3D永磁体对另一物体的磁力计算

zhangxw023

目前有个任务需要计算一永磁体对另一物体的吸引力，但是不知道计算结果对不对，所以建了个模型进行分析，发现结果有问题，但是不知道问题出在哪，请各位大侠指教，谢谢！
模型说明：两圆柱体，其中一个为永磁体，另一个定义为任何物体，可修改相对磁导率，两个圆柱体由稍大一点的圆柱形空气包1包裹，然后再由更大球形空气包2包裹，空气包2外表元设为无限远边界，设空气包1的目的是细化网格，无特殊意义。我在计算时设定试验物体的相对磁导率为1，也就是空气，理论上永磁体对其磁力为0，但结果计算出来力却很大，麻烦各位帮忙看一下


命令流如下，网格划分较细，各位可以酌情修改：
/prep7
emunit,mks
et,1,solid96         !磁体单元
et,2,infin111         !无限边界单元
!空气磁导率
mp,murx,1,1
!b-h曲线 定义永磁体
HC=895000
TB,BH,2,,26
TBPT,, 130.000000 , 0.100000000
TBPT,, 170.000000 , 0.200000000
TBPT,, 197.000000 , 0.300000000
TBPT,, 218.000000 , 0.400000000
TBPT,, 250.000000 , 0.500000000
TBPT,, 290.000000 , 0.600000000
TBPT,, 338.000000 , 0.700000000
TBPT,, 400.000000 , 0.800000000
TBPT,, 472.000000 , 0.900000000
TBPT,, 570.000000 , 1.00000000
TBPT,, 682.000000 , 1.10000000
TBPT,, 810.000000 , 1.20000000
TBPT,, 975.000000 , 1.30000000
TBPT,, 1600.00000 , 1.40000000
TBPT,, 2520.00000 , 1.50000000
TBPT,, 3520.00000 , 1.60000000
TBPT,, 4760.00000 , 1.70000000
TBPT,, 8300.00000 , 1.80000000
TBPT,, 12000.0000 , 1.90000000
TBPT,, 17000.0000 , 2.00000000
TBPT,, 23000.0000 , 2.10000000
TBPT,, 32000.0000 , 2.20000000
TBPT,, 42500.0000 , 2.30000000
TBPT,, 44500.0000 , 2.32000000
TBPT,, 48200.0000 , 2.35000000
TBPT,, 52200.0000 , 2.37000000
TBPLOT,BH,2,,,
mp,mgxx,2,0
mp,mgyy,2,0
mp,mgzz,2,hc   !极化方向定义为Z方向
mp,murx,3,1    ！定义试验材料3为空气，铁磁材料等，更改相对磁导率即可
!建模
radus=0.4
height=0.4
CYLIND,radus,0,0.1*height,height,0,360,
CYLIND,radus,0,-0.1*height,-1*height,0,360,
CYLIND,radus+0.05,0,-1.1*height,1.1*height,0,360, ！建立空气包1将永磁体和试验件包裹
SPHERE,0,0.8,0,360,                               ！建立空气包2
vovlap,all
numcmp,volu
vplot
!空气划分网格
！定义材料属性
VSEL, , , ,       1
VATT,       2, ,   1,       0
VSEL, , , ,       2
VATT,       3, ,   1,       0
VSEL, , , ,       3,4
VATT,       1, ,   1,       0
mshk,0               !设置为自由网格划分
msha,1,3d            !划分网格采用四面体单元
esize,0.04
!永磁体磁体、试验件、空气包1划分网格
ALLSEL,ALL
vmesh,1
vmesh,2
vmesh,4
esize,0.1
vmesh,3              !空气包2划分网格
！定义远场单元
ASEL,S,,,13,14 !选择空气包2外表面
CM,A_AIR,AREA ! 定义组件 A_AIR
CMSEL,S,A_AIR !选择组件 A_AIR
EXTOPT,ESIZE,1 !生成单元层
TYPE,2
MAT,1
LOCAL,12,2
VEXT,ALL,,,0.001
NUMCMP,ALL
CSYS,0
ASLV
CMSEL,U,A_AIR
NSLA,S,1
SF,ALL,INF !Define the infinite node
ESEL,ALL
NSEL,ALL
FINISH
！对需要进行力计算的部件施加边界条件
ALLSEL,ALL
ESEL,S,MAT,,2                ！定义永磁体为部件aaa
ALLSEL,BELOW,ELEM   
EPLOT   
CM,aaa,ELEM
ESEL,S,MAT,,3                ！定义试验体为部件bbb
ALLSEL,BELOW,ELEM   
EPLOT   
CM,bbb,ELEM  
/PREP7!
FMAGBC,'AAA'
FMAGBC,'BBB'
！求解
ALLSEL,ALL
MAGSOLV,2, , ,0.001,25,0

zhangxw023

分析结果： 圆柱形物体在磁场中受力，理论应该为0
              SUMMARY OF FORCES BY VIRTUAL WORK
Load Step Number:       1.
Substep Number:         1.
Time:             0.1000E+01
Units of Force:        ( N )   
Component        Force-X       Force-Y        Force-Z
    BBB       -0.56049E+03   0.13539E+04    0.11444E+06 ___________________________________________________                     
           SUMMARY OF FORCES BY MAXWELL STRESS TENSOR
Units of Force:  ( N )   
Component        Force-X       Force-Y        Force-Z
    BBB        0.21498E+03   0.87481E+02    0.16922E+06
_____________________________________________________                  
Note: Maxwell forces are in the Global Cartesian coordinate system.     
Virtual work forces are in the element ESYS coordinate system.         
The following element table items are available for printing and        
plotting of the forces obtained by Virtual Work and the Maxwell Stress  
Tensor methods.                                                         
     Element Item Name      Method      Direction
          FVW_X          Virtual Work       X
          FMX_X          Maxwell Stress     X
          FVW_Y          Virtual Work       Y
          FMX_Y          Maxwell Stress     Y
          FVW_Z          Virtual Work       Z
          FMX_Z          Maxwell Stress     Z
___________________________________________________________________