4.8.2.15. TRACK_SOIL_FORCE
This subroutine generates a user - defined soil force.
Language type |
Subroutine |
FORTRAN |
track_soil_force (TIME, INFO, UPAR, NPAR, DIRD, DIRV, DISP, LGORI, NGPOS, NGORI, NGVEL, LENGTH, WIDTH, IFLAG, RESULT) |
C/C++ |
track_soil_force (double time, int info[], double upar[], int npar, double dird[], double dirv[], double disp[], double lgori[], double ngpos[], double ngori[], double ngvel[], double length, double width, int iflag, double result[3]) |
Variable Name |
Size |
Description |
time |
double |
Current simulation time of RecurDyn/Solver. |
info |
int[3] |
ID of body & node. |
upar |
double[30] |
Parameters defined by the user. The maximum size of array is 30. |
npar |
int |
Number of user parameters. |
dird |
double[3] |
Direction vector of a ground patch. |
dirv |
double[3] |
Returned force vector and six-dimensional array. |
lgori |
double[9] |
Global orientation of a track link. |
ngpos |
dou ble[3] |
Global position of a contact node. |
ngvel |
double[3] |
Global velocity of a contact node |
length |
double |
Length of meshed rectangle at a track link. |
width |
double |
Width of meshed rectangle at a track link. |
iflag |
int |
When RecurDyn/Solver makes its initial call to this routine, the flag is true. Otherwise, the flag is false. |
result |
double[3] |
Returned force vector. |
Example
Figure 4.118 Track vehicle model using track soil user subroutine
Note
If the user wants to run this model using a user subroutine, the user can refer it in the directory (<install dir>\Help\usub\**).
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
#include "stdafx.h"
#include "DllFunc.h"
#include "math.h"
TrackSoilForce_API void __cdecl track_soil_force
(double time, int info[], double upar[], int npar, double dird[], double dirv[],
double disp[], double lgori[], double ngpos[], double ngori[], double ngvel[],
double length, double width, int iflag, double result[3])
{
using namespace rd_syscall;
// Parameter Information
// time : Simulation time of RD/Solver. (Input)
// info : ID of body & node & end_flag. (Input, 1 : TRUE)
// upar : parameters defined by user. (Input)
// npar : Number of user parameters. (Input)
// dird : Direction vector of ground patch. (Input)
// dirv : Velocity of a contact node relative to direction vector. (Input)
// disp : Sinkage & shear displacements. (Input)
// lgori : Global orientation of a track link. (Input)
// ngpos : Global position of a contact node. (Input)
// ngori : Gloabal orientation of a contact node. (Input)
// ngvel : Gloabal velocity of a contact node. (Input)
// length : Length of meshed rectangle at a track link. (Input)
// width : Width of meshed rectangle at a track link. (Input)
// iflag : When RD/Solver initializes arrays, the flag is true. (Input, -1)
// result : Returned track soil force vector. (Output, Size : 3)
// User Statement
int body_id, node_id, end_flag, i;
double 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;
//--- ASSIGN INFORMATION
body_id = info[0];
node_id = info[1];
end_flag = info[2];
//--- UNIT VECTOR OF DIRECTION
for(i=0; i<3; i++){
u_x[i] = dird[i];
u_y[i] = dird[3+i];
u_z[i] = dird[6+i];
}
//--- VEOCITY OF A CONTACT NODE RELATIVE TO DIRECTION VECTOR
vel_x = dirv[0];
vel_y = dirv[1];
vel_z = dirv[2];
//-------------------------------------------------------
//--- 1. PRESSURE - SINKAGE RELATIONSHIP
//-------------------------------------------------------
//--- SINKAGE & SHEAR DISPLACEMENTS
z = disp[0];
delta_sx = disp[1];
delta_sz = disp[2];
z_max = disp[3];
//--- DEFINE SOIL PARAMETERS
k_c = 0.00047613;
k_phi = 0.00076603;
n = 1.1;
cohesion = 0.00104;
phi = 28.0*acos(-1.0)/180.0;
m_k = 25.0;
//--- CALCULATE AREA RATIO RELATIVE TO A NORMAL VECTOR OF TRIANGULAR PATHCH
ratio_n = fabs( ngori[3]*u_y[0] + ngori[4]*u_y[1]
+ ngori[5]*u_y[2] );
area = length*width*ratio_n;
//--- FOR REPETITIVE LOAD
ratio_zmax = 0.01;
//--- CHECK IF SINKAGE IS IN THE LOADING OR UNLOADING/RELOADING CONDITION
//--- PLASTIC EFFECT
if( (z>=0.0) && (z>=z_max) ) {
pre = ( k_c/width + k_phi ) *pow(z,n);
}
//--- ELASTIC EFFECT
else if ( (z>=0.0) && (z<z_max) && (z>=(z_max*(1.0-ratio_zmax))) ) {
pre_max = ( k_c/width + k_phi ) * pow(z_max,n);
ku = pre_max/(z_max*ratio_zmax);
pre = pre_max - ku*(z_max-z);
}
else {
pre = 0.0;
}
force = pre * area;
//--- GET GLOBAL FORCE
for(i=0; i<3; i++){
result[i] = force * u_y[i];
}
//-------------------------------------------------------
//--- 2. SHEAR STRESS - SHEAR DISPLACEMENT RELATIONSHIP
//-------------------------------------------------------
if( vel_x>=0.0 ) {
sgx = -1.0;
}
else {
sgx = 1.0;
}
if( vel_z>=0.0 ) {
sgz = -1.0;
}
else {
sgz = 1.0;
}
//----/ LONGITUDINAL FORCE /
pre_sx = ( cohesion + pre*tan(phi) )
*( 1.0 - exp(-delta_sx/m_k) );
area = length*width*ratio_n;
force_x = pre_sx * area;
for(i=0; i<3; i++){
result[i] = result[i] + fabs(force_x)*sgx*u_x[i];
}
//----/ LATERAL FORCE /
pre_sz = ( cohesion + pre*tan(phi) )
*( 1.0 - exp(-delta_sz/m_k) );
area = length*width*ratio_n;
force_z = pre_sz * area;
for(i=0; i<3; i++){
result[i] = result[i] + fabs(force_z)*sgz*u_z[i];
}
}