How can I compute harmonics for my magnetic solution?


The magnetic field values can be stored in an APDL array using a
path operation and then fourier coefficients can be computed. The sample
input and macro below demonstrate the procedure on a dummy 30 degree model.

! Create a dummy 30 degree emag model.

/prep7
csys,1
k,1,1
k,2,1,10
k,3,1,20
k,4,1,30

k,5,1.5
k,6,1.5,10
k,7,1.5,20
k,8,1.5,30

k,9,2
k,10,2,10
k,11,2,20
k,12,2,30

k,13,2.2
k,14,2.2,10
k,15,2.2,20
k,16,2.2,30

k,17,3
k,18,3,10
k,19,3,20
k,20,3,30

a,1,5,6,2
*repe,3,1,1,1,1

a,5,9,10,6
*repe,3,1,1,1,1

a,9,13,14,10
*repe,3,1,1,1,1

a,13,17,18,14
*repe,3,1,1,1,1

et,1,53

mp,murx,1,1

mp,murx,2,1 ! magnet
mp,mgxx,2,1/12.57e-7

mp,murx,3,1000 ! iron

esize,,10

mat,1
amesh,7,9

mat,2
amesh,5

mat,3
amesh,all

nsel,s,loc,x,1
nsel,a,loc,x,3
d,all,az,0
nsel,all

fini

/solu
solve
fini

/post1
set,last

lpath,node(2.1,0,0),node(2.1,5,0),node(2.1,10,0),node(2.1,15,0),node(2.1,20,0),node(2.1,25,0),node(2.1,30,0)

pdef,B_Sum,B,Sum

plpath,b_sum
paget,b_vals,table

parsav,all ! save array to file.parm

PARRES,NEW,file,parm,
*dim,btab,table,289,1
*SET,pi,acos(-1)
*vfact,180/pi ! convert path x variable to degrees from 0 to 360
*vfun,btab(1,0),copy,b_vals(1,4)
*vfun,btab(1,1),copy,b_vals(1,5)

fourier,'btab','fout',20 !

*stat,fout ! fourier coefficients

###############
FOURIER.MAC
!
! arg1 is the table array name containing the force/angle data
! angle as the index column and the force in column 1
! arg2 is array type array to be created and to contain the
! mode, isym, and coefficients of the evaluaion points
! in columns 1, 2, and 3 respectively.
! arg3 is the quantity of coefficients to consider - minimum of 2
! and the input may be increased by 1 in the macro if applicable
! arg4 is the number of curve data points to coefficients ratio
! (minimum is 2)
! arg5 is a switch touse: 0 for sine and cosine terms, -1 for sine
! only, and1 for cosine only.
! arg6 is the initial mode to be considered (default is 0)
! arg7 is the evaluation increment (default 1 degree)
! arg8 is the quantity of coefficients to evaluate (defaults to
! all calcualted coefficients - arg3)
! arg9 is the parameter file save/delete key - 0 delete (default)
! any other value will be the file extension on file 'parmset'
! ar18, ar19, ar23, ar24, ar25, ar26, ar27, ar28 are internally used
!
ar23=nint(2>arg3)
ar24=(2>nint(arg4))
*if,arg5,eq,0,then
ar25=0
*else
ar25=(arg5/abs(arg5))
*endif
ar26=(0>nint(arg6))
*IF,ar25,EQ,0,THEN ! sine and cosine terms
*if,ar26,eq,0,then
*IF,NINT(ar23/2-0.1),EQ,(ar23/2),THEN
ar23=(ar23+1)
*ENDIF
*elseif,nint(ar23/2-.1),ne,(ar23/2)
ar23=ar23+1
*endif
*elseif,ar25,eq,-1 ! sine terms only
*if,ar26,eq,0,then
ar26=1
*endif
*ENDIF
*if,arg7,eq,0,then
ar27=1
*elseif,arg7,gt,30
ar27=30
*else
ar27=arg7
*endif
ar28=nint(abs(arg8))
*if,ar28,eq,0,then
ar28=ar23
*elseif,ar28,gt,ar23
ar28=ar23
*endif

ar18=ar23*ar24
ar19=nint(360/ar27)
ar19=(ar19>(ar28+1))

*DIM,MODE,,ar23
*DIM,ISYM,,ar23
*DIM,COEFF,,ar23
*DIM,THETA,,ar18
*DIM,FORCE_IN,,ar18
!
! Fill in mode/isym vector values
!
*if,ar25,eq,0,then ! sine and cosine terms
*if,ar26,eq,0,then
*vlen,((ar23-1)/2)
*vfill,mode(1),ramp,ar26,1
*vfill,mode((ar23-1)/2+1),ramp,(1>ar26),1
*vlen,((ar23-1)/2)
*vfill,isym(1),ramp,1
*vfill,isym((ar23-1)/2+1),ramp,-1
*else
*vlen,(ar23/2)
*vfill,mode(1),ramp,ar26,1
*vfill,mode(ar23/2+1),ramp,(1>ar26),1
*vlen,((ar23)/2)
*vfill,isym(1),ramp,1
*vfill,isym(ar23/2+1),ramp,-1
*endif
*else ! sine only or cosine only terms
*vfill,mode(1),ramp,ar26,1
*vfill,isym(1),ramp,ar25
*endif
!
*VFILL,THETA(1),RAMP,0,(360/ar18)
*vitrp,force_in(1),arg1(1),theta(1)
!
! Calculate coefficients
!
*MFOURI,FIT,COEFF(1),MODE(1),ISYM(1),THETA(1),FORCE_IN(1)
!
*dim,%arg2%,,ar23,3
*VFUN,arg2(1,1),COPY,MODE(1)
*VFUN,arg2(1,2),COPY,ISYM(1)
*VFUN,arg2(1,3),COPY,COEFF(1)
*DIM,OLDORDER,,ar23
*VABS,,1 ! USE ABSOLUTE VALUE
*MOPER,OLDORDER(1),arg2(1,1),SORT,arg2(1,3)
*DIM,DESCEND,,ar23
*VFILL,DESCEND(1),RAMP,ar23,-1
*MOPER,OLDORDER(1),arg2(1,1),SORT,DESCEND(1)
! print out mode/isym/coefficients
/nopr
*vwrite
('mode isym coefficient')
*vlen,ar28
*vwrite,arg2(1,1),arg2(1,2),arg2(1,3)
(f5.0,3x,f5.0,3x,g16.9)
/gopr
!
! Evaluate coefficients
*dim,EV_MODE,,ar28
*dim,EV_ISYM,,ar28
*dim,EV_COEFF,,ar28
*VFUN,EV_MODE(1),COPY,arg2(1,1)
*VFUN,EV_ISYM(1),COPY,arg2(1,2)
*VFUN,EV_COEFF(1),COPY,arg2(1,3)

!
*DIM,EVAL_TH,,ar19
*DIM,EVAL_EV,,ar19
*DIM,EVAL,TABLE,ar19,2
!
*VFILL,EVAL_TH(1),RAMP,0,(360/ar19)
*vitrp,EVAL(1,1),arg1(1),EVAL_TH(1)
!
*MFOURI,EVAL,EV_COEFF(1),EV_MODE(1),EV_ISYM(1),EVAL_TH(1),EVAL_EV(1)
!
! Graphically compare results to input
!
EVAL(0,0)=1
EVAL(0,1)=1
EVAL(0,2)=2
*MFUN,EVAL(1,0),COPY,EVAL_TH(1)
*MFUN,EVAL(1,2),COPY,EVAL_EV(1)
!
*VPLOT,EVAL(1,0),EVAL(1,1),2
!
*IF,arg9,ne,0,THEN
! save results
PARSAV,ALL,parmset,%arg9%
*endif
*set,mode(1)
*set,isym(1)
*set,theta(1)
*set,force_in(1)
*set,coeff(1)
*set,EVAL(1,1)
*SET,EVAL_TH(1)
*SET,EVAL_EV(1)
*SET,EV_MODE(1)
*SET,EV_ISYM(1)
*SET,EV_COEFF(1)
*SET,OLDORDER(1)
*SET,DESCEND(1)





Show Form
No comments yet. Be the first to add a comment!