|
SUBROUTINE TRACK_SOIL_FORCE & (TIME,INFO,UPAR,NPAR,DIRD,DIRV,DISP, & LGORI,NGPOS,NGORI,NGVEL,LENGTH,WIDTH,IFLAG, & RESULT) C---- TO EXPORT * SUBROUTINE !DEC$ ATTRIBUTES DLLEXPORT,C::TRACK_SOIL_FORCE
c---- INCLUDE SYSTEM CALL INCLUDE 'SYSCAL.F'
c---- DEFINE VARIABLES c Parameter Information c TIME : Simulation time of RD/Solver. (Input) c INFO : ID of body & node & end_flag. (Input, 1 : TRUE) c UPAR : parameters defined by user. (Input) c NPAR : Number of user parameters. (Input) c DIRD : Direction vector of ground patch. (Input) c DIRV : Velocity of a contact node relative to direction vector. (Input) c DISP : Sinkage & shear displacements. (Input) c LGORI : Global orientation of a track link. (Input) c NGPOS : Global position of a contact node. (Input) c NGORI : Global orientation of a contact node. (Input) c NGVEL : Global velocity of a contact node. (Input) c LENGTH : Length of meshed rectangle at a track link. (Input) c WIDTH : Width of meshed rectangle at a track link. (Input) c IFLAG : When RD/Solver initializes arrays, the flag is true. (Input, -1) c RESULT : Returned track soil force vector. (Output, Size : 3)
DOUBLE PRECISION TIME, UPAR(*) INTEGER INFO(3), NPAR DOUBLE PRECISION DIRD(9), DIRV(9), DISP(4) DOUBLE PRECISION LGORI(9), NGPOS(3), NGORI(9), NGVEL(3) DOUBLE PRECISION LENGTH, WIDTH LOGICAL IFLAG DOUBLE PRECISION RESULT[REFERENCE](3)
c---- USER STATEMENT INTEGER BODY_ID, NODE_ID, END_FLAG DOUBLE PRECISION K_C, K_PHI, N, COHESION, PHI, M_K, AREA, PRE_MAX, & Z, Z_MAX, PRE, FORCE, SGX, SGZ, PRE_SX, PRE_SZ, & DELTA_SX, DELTA_SZ, FORCE_X, FORCE_Z, & U_X(3), U_Y(3), U_Z(3), VEL_X, VEL_Y, VEL_Z, & KU, RATIO_ZMAX, RATIO_N
c---- ASSIGN INFORMATION BODY_ID = INFO(1) NODE_ID = INFO(2) END_FLAG = INFO(3)
c---- UNIT VECTOR OF DIRECTION DO I = 1, 3 U_X(I) = DIRD(I) U_Y(I) = DIRD(3+I) U_Z(I) = DIRD(6+I) ENDDO
c---- VEOCITY OF A CONTACT NODE RELATIVE TO DIRECTION VECTOR VEL_X = DIRV(1) VEL_Y = DIRV(2) VEL_Z = DIRV(3)
C-------------------------------------------------------- C---- 1. PRESSURE - SINKAGE RELATIONSHIP C-------------------------------------------------------- C---- SINKAGE & SHEAR DISPLACEMENTS Z = DISP(1) DELTA_SX = DISP(2) DELTA_SZ = DISP(3) Z_MAX = DISP(4)
c---- DEFINE SOIL PARAMETERS K_C = 0.00047613D0 K_PHI = 0.00076603D0 N = 1.1D0 COHESION = 0.00104d0 PHI = 28.0D0*DACOS(-1.0D0)/180.0D0 M_K = 25.0D0
C---- CALCULATE AREA RATIO RELATIVE TO A NORMAL VECTOR OF TRIANGULAR PATHCH RATIO_N = DABS( NGORI(4)*U_Y(1) + NGORI(5)*U_Y(2) & + NGORI(6)*U_Y(3) ) AREA = LENGTH*WIDTH*RATIO_N
C---- FOR REPETITIVE LOAD RATIO_ZMAX = 0.01D0
C---- CHECK IF SINKAGE IS IN THE LOADING OR UNLOADING/RELOADING CONDITION c---- PLASTIC EFFECT IF( Z.GE.0.0D0 .AND. Z.GE.Z_MAX ) THEN PRE = ( K_C/WIDTH + K_PHI ) *Z**N
c---- ELASTIC EFFECT ELSE IF( Z.GE.0.0D0 .AND. Z.LT.Z_MAX & .AND. Z.GE.Z_MAX*(1.0D0-RATIO_ZMAX) ) THEN PRE_MAX = ( K_C/WIDTH + K_PHI ) * Z_MAX**N KU = PRE_MAX/(Z_MAX*RATIO_ZMAX) PRE = PRE_MAX - KU*(Z_MAX-Z) ELSE PRE = 0.0D0 ENDIF
FORCE = PRE * AREA
c---- GET GLOBAL FORCE DO I = 1, 3 RESULT(I) = FORCE * U_Y(I) ENDDO
C-------------------------------------------------------- C---- 2. SHEAR STRESS - SHEAR DISPLACEMENT RELATIONSHIP C-------------------------------------------------------- SGX = - DSIGN(1.0D0,VEL_X) SGZ = - DSIGN(1.0D0,VEL_Z)
c-----/ LONGITUDINAL FORCE / PRE_SX = ( COHESION + PRE*DTAN(PHI) ) & *( 1.0D0 - DEXP(-DELTA_SX/M_K) ) AREA = LENGTH*WIDTH*RATIO_N FORCE_X = PRE_SX * AREA DO I = 1, 3 RESULT(I) = RESULT(I) + DABS(FORCE_X)*SGX*U_X(I) ENDDO
c-----/ LATERAL FORCE / PRE_SZ = ( COHESION + PRE*DTAN(PHI) ) & *( 1.0D0 - DEXP(-DELTA_SZ/M_K) ) AREA = LENGTH*WIDTH*RATIO_N FORCE_Z = PRE_SZ * AREA DO I = 1, 3 RESULT(I) = RESULT(I) + DABS(FORCE_Z)*SGZ*U_Z(I) ENDDO
RETURN END |