APDL梁单元转子建模及临界转速计算命令流

*STEP,1,1e-3,1,10
*SET,ANGULAR_VELOCITY,1
*SET,DENSITY,7850
*SET,MODULUS_OF_ELASTICITY,210e9
*SET,POISSON_RATIO,0.3
*SET,LENGTH,0.5
*SET,WIDTH,0.1
*SET,THICKNESS,0.01
*SET,NUMBER_OF_ELEMENTS,20
*SET,RADIUS,0.2
*SET,NUMBER_OF_NODES,NUMBER_OF_ELEMENTS+1
*SET,NUMBER_OF_DOFS,2*NUMBER_OF_NODES
*SET,YOUNGS_MODULUS,MODULUS_OF_ELASTICITY/(1-POISSON_RATIO^2)
*SET,INERTIA,WIDTH*THICKNESS^3/12
*SET,MASS,DENSITY*LENGTH*WIDTH*THICKNESS
*SET,MASS_MATRIX,MASS/6*[2,1;1,2]
*SET,STIFFNESS_MATRIX,YOUNGS_MODULUS*INERTIA/LENGTH*[1,-1;-1,1]
*SET,ROTATION_MATRIX,[cos(ANGULAR_VELOCITY*TIME),-sin(ANGULAR_VELOCITY*TIME);sin(ANGULAR_VELOCITY*TIME),cos(ANGULAR_VELOCITY*TIME)]
*SET,INITIAL_DISPLACEMENTS,[0;0;0;0]
*SET,INITIAL_VELOCITIES,[0;0;0;0]
*SET,TIME_STEP,1e-5
*DIM,COORDINATES,NUMBER_OF_NODES,3
*DIM,DOF_LIST,NUMBER_OF_DOFS
*DIM,NODE_LIST,NUMBER_OF_NODES
*VGEN,NODE_LIST,1,NUMBER_OF_NODES
*DO,I,1,NUMBER_OF_NODES
COORDINATES(I,1)=RADIUS*cos((I-1)/(NUMBER_OF_NODES-1)*pi)
COORDINATES(I,2)=RADIUS*sin((I-1)/(NUMBER_OF_NODES-1)*pi)
*ENDDO
*DO,I,1,NUMBER_OF_NODES
DOF_LIST(2*I-1)=2*I-1
DOF_LIST(2*I)=2*I
*ENDDO
*K,M,DOF_LIST,NUMBER_OF_DOFS,NUMBER_OF_NODES,COORDINATES,NUMBER_OF_ELEMENTS,1,STIFFNESS_MATRIX,MASS_MATRIX
*SET,INITIAL_DISPLACEMENTS,[0;0;0;0;0;0;0;0]
*SET,INITIAL_VELOCITIES,[0;0;0;0;0;0;0;0]
*SET,TIME,0
*DO WHILE (TIME<1)
*SOLVE,TRAN
*SET,ROTATION_MATRIX,[cos(ANGULAR_VELOCITY*(TIME+TIME_STEP)), -sin(ANGULAR_VELOCITY*(TIME+TIME_STEP)); sin(ANGULAR_VELOCITY*(TIME+TIME_STEP)), cos(ANGULAR_VELOCITY*(TIME+TIME_STEP))]
*SET,INITIAL_DISPLACEMENTS,ROTATION_MATRIX*INITIAL_DISPLACEMENTS
*SET,INITIAL_VELOCITIES,ROTATION_MATRIX*INITIAL_VELOCITIES+[0;0;0;0;0;0;ANGULAR_VELOCITY;0]
*SET,TIME,TIME+TIME_STEP
*ENDDO
*SET,EIGENVALUE,0
*SET,EIGENVECTOR,0
*SET,NUMBER_OF_MODES,2
*SET,EIGENVALUE_MATRIX,NUMBER_OF_DOFS,NUMBER_OF_DOFS
*SET,EIGENVECTOR_MATRIX,NUMBER_OF_DOFS,NUMBER_OF_DOFS
*SET,COMPLEX_MATRIX,NUMBER_OF_DOFS,2*NUMBER_OF_DOFS
*SOLVE,EIGEN
*DO,I,1,NUMBER_OF_MODES
EIGENVALUE(I)=EIGENVALUE_MATRIX(2*I-1,2*I-1)
EIGENVECTOR(:,I)=EIGENVECTOR_MATRIX(:,2*I-1)
*ENDDO
*SET,UNBALANCE_VECTOR,[0;0;0;0;0;0;1;0]
*SET,CORRECTION_FACTOR,0.001
*SET,ITERATION_NUMBER,0
*SET,TOLERANCE,1e-6
*SET,MAXIMUM_ITERATION_NUMBER,100
*SET,ANGULAR_VELOCITY_INCREMENT,0.1
*SET,CRITICAL_ANGULAR_VELOCITY,0
*DO WHILE (CRITICAL_ANGULAR_VELOCITY==0 AND ITERATION_NUMBER<MAXIMUM_ITERATION_NUMBER)
*SET,ITERATION_NUMBER,ITERATION_NUMBER+1
*SET,ANGULAR_VELOCITY,ANGULAR_VELOCITY+ANGULAR_VELOCITY_INCREMENT
*SET,ROTATION_MATRIX,[cos(ANGULAR_VELOCITY*TIME),-sin(ANGULAR_VELOCITY*TIME);sin(ANGULAR_VELOCITY*TIME),cos(ANGULAR_VELOCITY*TIME)]
*SET,INITIAL_DISPLACEMENTS,ROTATION_MATRIX*INITIAL_DISPLACEMENTS
*SET,INITIAL_VELOCITIES,ROTATION_MATRIX*INITIAL_VELOCITIES+[0;0;0;0;0;0;ANGULAR_VELOCITY;0]
*SET,UNBALANCE_FORCE,UNBALANCE_VECTOR'*ROTATION_MATRIX*EIGENVECTOR*EIGENVECTOR'*ROTATION_MATRIX'*UNBALANCE_VECTOR
*SET,CRITICAL_ANGULAR_VELOCITY,ANGULAR_VELOCITY*(1+CORRECTION_FACTOR*UNBALANCE_FORCE/EIGENVALUE(1))
*SET,ANGULAR_VELOCITY_INCREMENT,0.01*ANGULAR_VELOCITY_INCREMENT*abs((CRITICAL_ANGULAR_VELOCITY-ANGULAR_VELOCITY)/ANGULAR_VELOCITY)
*ENDDO
*PRINT,CRITICAL_ANGULAR_VELOCITY
*FINISH