! FEA computer laboratory NL6: POST BUCKLING OF A PLATE WITH CENTRAL HOLE
!! Initial setup
! White canvas
/RGB,INDEX,100,100,100,0
/RGB,INDEX,0,0,0,15
/ESHAPE,1.0 
! [N], [Nm], [N/m^2=Pa], [N/m]
!
! Define constant parameters
A=0.25
R=0.04
t1=0.25E-3
t=3E-3
!
F=180E3
!
Exx=127.6E9
Eyy=11.03E9
Ezz=11.03E9
nuxy=0.27
nuyz=0.03
nuxz=0.03
Gxy=5.74E9
Gyz=5.74E9
Gxz=5.74E9
!
PI=4*ATAN(1)
!! SETUP: element types, materials, layup
/PREP7 
!
ET,1,SHELL181   
KEYOPT,1,8,1
!
ESYS,0
!  
MPTEMP,1,0  
MPDATA,EX,1,,Exx
MPDATA,EY,1,,Eyy
MPDATA,EZ,1,,Ezz
MPDATA,NUXY,1,,nuxy 
MPDATA,NUYZ,1,,nuyz 
MPDATA,NUXZ,1,,nuxz 
MPDATA,GXY,1,,Gxy   
MPDATA,GYZ,1,,Gyz   
MPDATA,GXZ,1,,Gxz
!
! Define shell layup
SECTYPE,1,SHELL
SECDATA, t1,1,30,3  
SECDATA, t1,1,-30,3 
SECDATA, t1,1,30,3  
SECDATA, t1,1,-30,3 
SECDATA, t1,1,30,3  
SECDATA, t1,1,-30,3
!
SECDATA, t1,1,-30,3 
SECDATA, t1,1,30,3  
SECDATA, t1,1,-30,3 
SECDATA, t1,1,30,3  
SECDATA, t1,1,-30,3 
SECDATA, t1,1,30,3  
SECOFFSET,MID   
!*
SECPLOT, 1 !<a href="images/layup.png" target="_blank">Click here to view image</a>
!! Modeling
K,1,0,0
K,2,A,0
K,3,A,A
K,4,0,A
!
A,1,2,3,4
CYL4,A/2,A/2,R  
ASBA,1,2
!
NUMCMP,ALL
!! Meshing
*DO,III,1,4
LESIZE,III,,,20
*ENDDO
!
*DO,III,5,8
LESIZE,III,,,10
*ENDDO
!
AMESH,1 !<a href="images/mesh.png" target="_blank">Click here to view image</a>
!
! detach mesh
MODMSH,DETACH
! Move nodes to 'help' the formation of buckling 
*GET,NNUMBER,NODE,0,NUM,MAX
!
*DO,I,1,NNUMBER
TMPZ = SIN(PI*NX(I)/A)*SIN(PI*NY(I)/A)*(3E-5)
N,I,NX(I),NY(I),TMPZ
*ENDDO
FINISH
!! Solution
/SOLU   
! Setting up the solution  
ANTYPE,0
NLGEOM,1
NSUBST,300,1000,100
OUTRES,ALL,ALL
! Setting the kinematic constraints on lines
NSEL,S,LOC,X,0
D,ALL,ALL,0
!
NSEL,S,LOC,Y,0
NSEL,A,LOC,Y,A
D,ALL,UZ,0
D,ALL,ROTY,0
!
NSEL,S,LOC,X,A
D,ALL,UZ,0
D,ALL,ROTX,0
D,ALL,ROTY,0
!
ALLSEL,ALL
D,ALL,ROTZ,0
!
! Setting the loading on 
F,NODE(A,A/2,0),FX,-F
FINISH
!
! Coupling the x displacements
/PREP7
NSEL,S,LOC,X,A
CP,1,UX,ALL
ALLSEL,ALL !<a href="images/constraints.png" target="_blank">Click here to view image</a>
FINISH
/SOL
SOLVE
FINISH
! Time history postprocessing
/POST26 
NSOL,2,NODE(A/2,A/2+R,0),U,Z, UZ_A
NSOL,3,NODE(A/2+R,A/2,0),U,Z, UZ_B
NSOL,4,NODE(A,A/2,0),U,X, UX_C
/AXLAB,X,'Time [s]'
/AXLAB,Y,'Displacement [m]'
XVAR,1
PLVAR,2,3,4 !<a href="images/displacement.png" target="_blank">Click here to view image</a>
FINISH
! Display the displaced structure at the last timestep
/POST1
SET,LAST
/VIEW,1,1,1,1
/DSCALE,ALL,4
PLNSOL, U,Z, 0,1.0 !<a href="images/disp_struct_uz.png" target="_blank">Click here to view image</a>
!<a href="images/ud.mp4" target="_blank">Click here to view animation</a>
! Load the macro for saving the displacement values at a given load magnitude
/INPUT,LOAD_DISP,txt !<a href="LOAD_DISP.txt" target="_blank">LOAD DISP MACRO</a>
!<a href="AZD.OUT" target="_blank">Z displacement of point A</a>
!<a href="BZD.OUT" target="_blank">Z displacement of point B</a>
!<a href="CXD.OUT" target="_blank">X displacement of point C</a>
!<a href="images/LOAD_DISP.png" target="_blank">Load value as a function of displacement</a>
!<a href="LOAD_DISP_PLOT.m" target="_blank">MATLAB file for plotting the data in files</a>