!
=======================================
! USER INPUT
! =======================================
! - define a Named Selection called CM_Section_Forces scoped to the body of
interest
! - define a Local Coordinate System at the starting point, with it's local Z axis
in the direction
! that you want to create the bending moment diagram. Moment are calculated about
this origin, which moves along the body with each increment.
! ARG1 - coordinate system (number) at the starting point, with it's local Z axis
in the direction that you want to create the bending moment diagram
! ARG2 - distance, measured from the local coordinate system (in ARG1) over which
you want to calculate the bending moment
! ARG3 - number of increments to split the distance into to calculate results
! Note: a text file of results will automatically open when the calculation is
complete, but you'll need to close this to complete the update
! =======================================
! AUTHOR AND DATE
! =======================================
! Nick Stevens - Tor Engineering Ltd
! September 2016
! http://tor-eng.com
! Thanks to SimuTech Group, in particular this article:
http://www.simutechgroup.com/FEA/fea-tips-tricks-stress-linearization.html
! Stress resultants article: https://en.wikipedia.org/wiki/Stress_resultants
! =======================================
! INITIALISE MACRO
! =======================================
resume
!########################################################################
! If you've already solved the model, then any subsequently created coordinate
! systems or named selections won't be available, so instead you can create
! them here as a fudge to avoid resolving the model
!########################################################################
! create coordinate system - the number (12 in this case) is populated in ARG1
LOCAL, 12, 0, 675, 30, 80, , 90
! create named selection, which needs to be called CM_Section_Forces
esel,s,mat,,11,13
cm,CM_Section_Forces,elem
allsel
!########################################################################
/show,png
/GFILE,500 !image size
/CMAP,_TEMPCMAP_,CMP,,SAVE
/RGB,INDEX,100,100,100,0
/RGB,INDEX,0,0,0,15
increment=ARG2/(ARG3-1)
! dimension results arrays
*dim,asurf_Mx,table,ARG3
*dim,asurf_My,table,ARG3
*dim,asurf_Fz,table,ARG3
*dim,asurf_Vx,table,ARG3
*dim,asurf_Vy,table,ARG3
*dim,asurf_Tz,table,ARG3
*dim,asurf_Dist,table,ARG3
/post1
set,last
Local_CS=ARG1 !coordinate system
cmsel,s,CM_Section_Forces ! select elements of body/bodies of interest
nsle,s,active ! select nodes
WPCSYS,,Local_CS ! place working plane at origin of coord sys ARG1
/CPLANE,1 ! cutting plane is working plane
WPSTYL,,,,,,0,2,1 ! cartesian WP with triad shown
! =======================================
! START LOOP
! =======================================
*do,i,1,ARG3
! get origin and orientation of "Local_CS" coordinate system
*GET,x_orig,CDSY,Local_CS,LOC,X ! get the origin
*GET,y_orig,CDSY,Local_CS,LOC,Y
*GET,z_orig,CDSY,Local_CS,LOC,Z
*GET,THXY,CDSY,Local_CS,ANG,XY ! get the orientation
*GET,THYZ,CDSY,Local_CS,ANG,YZ
*GET,THZX,CDSY,Local_CS,ANG,ZX
! create surface
SUCR,mysurf,CPLANE,3
SUSEL,S,mysurf
! get info needed to rotate results in to local coordinate system
*get,thxy,CDSY,4,ANG,XY ! First rotation, about the Global Z axis, that put
working plane where it is
*get,thyz,CDSY,4,ANG,YZ ! Second rotation, about the new X axis, that put
working plane where it is
*get,thzx,CDSY,4,ANG,ZX ! Third rotation, about the newest Y axis, that put
working plane where it is
*get,origx,CDSY,4,LOC,X ! location of the origin of the working plane, in
Global Cartesian
*get,origy,CDSY,4,LOC,Y ! location of the origin of the working plane, in
Global Cartesian
*get,origz,CDSY,4,LOC,Z ! location of the origin of the working plane, in
Global Cartesian
! Form cosine and sin terms as used in a rotation matrix
*afun,deg
c1=cos(thxy) ! Cosine of THXY, first rotation about Global Z
s1=sin(thxy)
c2=cos(thyz) ! Cosine of THYZ, second rotation about the new X
s2=sin(thyz)
c3=cos(thzx) ! Cosine of THZX, third rotation about the newest Y
s3=sin(thzx)
*afun,rad ! return to default radians
! Form the 3x3 rotation matrix
R11=c3*c1-s3*s2*s1
R12=c3*s1+s3*s2*c1
R13=-s3*c2
R21=-c2*s1
R22=c2*c1
R23=s2
R31=s3*c1+c3*s2*s1
R32=s3*s1-c3*s2*c1
R33=c3*c2
SUCALC,MyZero,DA,ZERO,,,,0.0 ! a dummy surface result (of zeroes) that is used in
calculations later
SUCALC,GCXX,GCX,SUB,MyZero, , ,-origx
SUCALC,GCYY,GCY,SUB,MyZero, , ,-origy
SUCALC,GCZZ,GCZ,SUB,MyZero, , ,-origz
SUCALC,sr11,DA,ZERO,,,,r11 ! rotation matrix terms
SUCALC,sr12,DA,ZERO,,,,r12
SUCALC,sr13,DA,ZERO,,,,r13
SUCALC,sr21,DA,ZERO,,,,r21
SUCALC,sr22,DA,ZERO,,,,r22
SUCALC,sr23,DA,ZERO,,,,r23
SUCALC,sr31,DA,ZERO,,,,r31
SUCALC,sr32,DA,ZERO,,,,r32
SUCALC,sr33,DA,ZERO,,,,r33
! Form working plane coordinates
! First do X
SUCALC,WPX1,sr11,MULT,GCXX ! R11*X
SUCALC,WPX2,sr12,MULT,GCYY ! R12*Y
SUCALC,WPX3,sr13,MULT,GCZZ ! R13*Z
SUCALC,WPX,WPX1,ADD,WPX2
SUCALC,WPX,WPX,ADD,WPX3 ! Form X ON WORKING PLANE = R11*X+R12*Y+R13*Z
! Then do Y
SUCALC,WPY1,sr21,MULT,GCXX ! R21*X
SUCALC,WPY2,sr22,MULT,GCYY ! R22*Y
SUCALC,WPY3,sr23,MULT,GCZZ ! R23*Z
SUCALC,WPY,WPY1,ADD,WPY2
SUCALC,WPY,WPY,ADD,WPY3 ! Form Y ON WORKING PLANE = R21*X+R22*Y+R23*Z
! Then do Z (for checking - all points should be on the same plane)
SUCALC,WPZ1,sr31,MULT,GCXX ! R31*X
SUCALC,WPZ2,sr32,MULT,GCYY ! R32*Y
SUCALC,WPZ3,sr33,MULT,GCZZ ! R33*Z
SUCALC,WPZ,WPZ1,ADD,WPZ2
SUCALC,WPZ,WPZ,ADD,WPZ3 ! Form Z ON WORKING PLANE = R31*X+R32*Y+R33*Z
RSYS,ARG1 !results coordinate system using local CS
! Mapping onto the surface
SUMAP,szz,S,Z ! map SZ (axial) onto the surface
SUMAP,sxz,S,XZ ! map shear stress SXY
SUMAP,syz,S,YZ ! map shear stress SYZ
! Form products in anticipation of calculating stress resultants
sucalc,szz_y,szz,MULT,WPY
sucalc,szz_x,szz,MULT,WPX
sucalc,s_tors1,syz,MULT,WPX
sucalc,s_tors2,sxz,MULT,WPY
SUCALC,s_tors,s_tors1,SUB,s_tors2
! Form integrals. This produces stress resultants about the origin of the local
coordinate system
sueval,Mx,szz_y,INTG
sueval,My,szz_x,INTG
sueval,Fz,szz,INTG
sueval,Vx,sxz,INTG
sueval,Vy,syz,INTG
sueval,Tz,s_tors,INTG
! Finally map results to output arrays
asurf_Mx(i)=Mx
asurf_My(i)=-My !negative sign required for My - see Wikipedia on Stress
Resultants
asurf_Fz(i)=Fz
asurf_Vx(i)=Vx
asurf_Vy(i)=Vy
asurf_Tz(i)=Tz
asurf_Dist(i)=(i-1)*increment
! Offset the working plane and local coordinate system ready to go round the loop
again
WPOFFS,,,increment
CSWPLAN,ARG1,0
*enddo
! =======================================
! END LOOP
! =======================================
! =======================================
! WRITE AND OPEN RESULTS
! =======================================
*cfopen,Section_Results,txt ! creates file called
Section_Results.txt, which will be in Solver Files Directory
*vwrite,'Dist','Vx','Vy','Fz','Mx','My','Tz'
(7(2X,A16))
*vwrite,asurf_Dist(1),asurf_Vx(1),asurf_Vy(1),asurf_Fz(1),asurf_Mx(1),asurf_My(1),a
surf_Tz(1)
(7(2X,F16.2))
*cfclos
/title,Vx against Distance
/GCOLUMN, 1, Vx
/AXLAB, X, Distance
/AXLAB, Y, Shear Force Vx
*VPLOT, asurf_Dist, asurf_Vx
/title,Vy against Distance
/GCOLUMN, 1, Vy
/AXLAB, X, Distance
/AXLAB, Y, Shear Force Vy
*VPLOT, asurf_Dist, asurf_Vy
/title,Fz against Distance
/GCOLUMN, 1, Fz
/AXLAB, X, Distance
/AXLAB, Y, Axial Force Fz
*VPLOT, asurf_Dist, asurf_Fz
/title,Mx against Distance
/GCOLUMN, 1, Mx
/AXLAB, X, Distance
/AXLAB, Y, Moment Mx
*VPLOT, asurf_Dist, asurf_Mx
/title,My against Distance
/GCOLUMN, 1, My
/AXLAB, X, Distance
/AXLAB, Y, Moment My
*VPLOT, asurf_Dist, asurf_My
/title,Tz against Distance
/GCOLUMN, 1, Tz
/AXLAB, X, Distance
/AXLAB, Y, Torsion Moment Tz
*VPLOT, asurf_Dist, asurf_Tz
/VIEW,1,1,1,1
/TRIA,LBOT
/PSYMB,CS,1
/title,Orientation System and Surface
WPSTYL,,,,,,0,0,1 ! cartesian WP with triad shown
nplot
!supl,mysurf,WPZ,1
! opens the created results file in Notepad:
/sys,notepad Section_Results.txt
/CMAP,_TEMPCMAP_,CMP
/show,close
rsys,0
csys,0
allsel
