title, Forces between two permanent magnets
/nopr
!npts=16
!*dim,dist,table,npts
!*dim,F,table,npts,2
d=0
d_air=0
!*do,i,1,npts
!parsave,all
!/clear,nostart             
!parres
! Problem parameters
Pi=acos(-1)
mu0=4*Pi*1e-7                     ! permeability of free space, H/m
mu0r=1                            ! relative permeability of free space
Br=0.38                           ! residual induction, tesla
a=10 $ b=6 $ c=3                  ! lower magnet dimensions
aa=6 $ bb=10 $ cc=3               ! upper magnet dimensions
xo=-4 $ yo=-4 $ zo=8              ! upper magnet center offset
/com, *** Upper Magnet Displacement d = %d% (mm)
/PREP7
! Solid model
block,-a,a,-b,b,-c,c              ! lower magnet volume
wpoff,xo,yo,zo
block,-aa+d,aa+d,-bb,bb,-cc,cc    ! upper (moving) magnet volume
vplot
wpoff,-xo,-yo,-zo
block,-30,30+d_air,-34,26,-23,31  ! air box
vovlap,all  
! Finite element model
esize,1
et,1,SOLID236                     ! magnetic solid
keyop,1,7,1                       ! condense forces to the corner nodes
vmesh,1,2  
esize,4
numcmp,volu
msha,1,3D 
vmesh,3 
! Material properties
mp,mgzz,2,Br/mu0
mp,murx,2,mu0r
mp,murx,1,mu0r
vsel,s,,,1,2,,1                 
emod,all,mat,2 
alls
nsel,s,ext             
d,all,az,0                        ! flux parallel boundary conditions
alls
vlscale,all,,,1e-3,1e-3,1e-3,,,1  ! scale the volumes to mm
WPSTYLE,,,,,,,,0
/ANG,1
/VIEW,1,1,1,1
/VUP,1,Z
esel,s,mat,,2
eplot
alls
fini
/SOLU
solve  
fini
/POST1
vsel,s,,,2,,,1                    ! select lower magnet along with the
                                  ! associated elements and nodes
esln
EMFT                              ! sum up magnetic forces
alls
esel,s,mat,,2
/VSCALE,,0.3
plvect,b,,,,vect,elem,on,on       ! plot magnetic field
plvect,fmag,,,,vect,node,on,on    ! plot magnetic forces
alls
fini
!/com, i = %i%
!dist(i)=d
!F(i,1)=_fxsum                    ! FX sum calculated by EMFT
!F(i,2)=_fysum                    ! FY sum calculated by EMFT
!d=d+2                            ! upper magnet displacement update
!*if,d,gt,8,then
!d_air=d-8
!*endif
!*enddo
!/axlab,x,Distance d (mm)
!/axlab,y,Forces acting the magnet (N)
!/gcol,1,Fx
!/gcol,2,Fy
!*vplot,dist(1),F(1,1),2