4.8.2.2. CONTACT_FORCE
This subroutine generates a user - defined contact force.
Language type |
Subroutine |
FORTRAN |
contact force (time,upar,npar,pen,rvel,jflag,iflag,result) |
C/C++ |
contact force (double time, double upar[], int npar, double pen, double rvel[], int jflag, int iflag, double* result) |
Variable Name |
Size |
Description |
time |
double |
Current simulation time of RecurDyn/Solver. |
upar |
double[30] |
Parameters defined by the user. The maximum size of array is 30. |
npar |
int |
Number of user parameters. |
pen |
double |
Penetration at the contact point. |
rvel |
dou ble[3] |
Relative velocity between contact pair w.r.t. contact reference frame. |
jflag |
int |
When RecurDyn/Solver evaluates Jacobian, the flag is true. Otherwise, the flag is false. |
iflag |
int |
When RecurDyn/Solver makes its initial call to this routine, the flag is true. Otherwise, the flag is false |
result |
double[6] |
Returned force vector and six dimensional arrays. Though the returned force has six arrays, only three terms of the translational force are available. The reference frame of reference is contact point reference frame. |
Example
Figure 4.106 Contact force model using Contact force 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\**).
C---- SUB. CONTACT_FORCE
SUBROUTINE CONTACT_FORCE
& (TIME,UPAR,NPAR,PEN,RVEL,JFLAG,IFLAG,RESULT)
C---- TO EXPORT * SUBROUTINE
!DEC$ ATTRIBUTES DLLEXPORT,C::CONTACT_FORCE
C---- INCLUDE SYSTEM CALL
INCLUDE 'SYSCAL.F'
C---- DEFINE VARIABLES
C Parameter Information
C TIME : Simulation time of RD/Solver. (Input)
C UPAR : Parameters defined by user. (Input)
C NPAR : Number of user parameters. (Input)
C PEN : Penetration at the contact point. (Input)
C RVEL : Relative velocity between contact pair w.r.t. contact reference frame. (Input, Size : 3)
C JFLAG : When RD/Solver evaluates a Jacobian, the flag is true. (Input)
C IFLAG : When RD/Solver initializes arrays, the flag is true. (Input)
C RESULT : Returned contact force vector. (Output, Size : 6)
DOUBLE PRECISION TIME, UPAR(*), PEN, RVEL(*)
INTEGER NPAR
LOGICAL JFLAG, IFLAG
DOUBLE PRECISION RESULT[REFERENCE](6)
C---- USER STATEMENT
DOUBLE PRECISION DPEN, K, C, M1, M2, M3, STV, DTV, SFC, DFC
DOUBLE PRECISION SFO, DFO, TVEL, FC(3), ZERO, MU
DOUBLE PRECISION FRICTIONFORCE
INTEGER I
LOGICAL ERRFLG
DPEN = RVEL(3)
K = UPAR(1)
C = UPAR(2)
M1 = UPAR(3)
M2 = UPAR(4)
M3 = UPAR(5)
STV = UPAR(6)
DTV = UPAR(7)
SFC = UPAR(8)
DFC = UPAR(9)
ZERO = 0.0D0
C---- COMPUTE CONTACT NORMAL FORCE
IF ( DABS(PEN) .LT. 1.0D-12 ) THEN
SFO = 0.0D0
ELSE
SFO = K * (DABS(PEN)**M1)
ENDIF
IF ( DABS(DPEN) .LT. 1.0D-12 ) THEN
DFO = 0.0D0
ELSE
DFO = C * (DABS(PEN)**M3) * (DABS(DPEN)**M2)
ENDIF
IF ( DPEN .GT. ZERO ) DFO = - DFO
FC(3) = SFO + DFO
C---- COMPUTE FRICTION FORCE
TVEL = SQRT(RVEL(1)*RVEL(1) + RVEL(2)*RVEL(2))
FC(1) = 0.0D0
FC(2) = 0.0D0
IF ( TVEL .GE. 1.0D-12 ) THEN
IF ( TVEL .LT. STV ) THEN
CALL RD_HAVSIN(DABS(TVEL),-STV,-SFC,STV,SFC,0,MU,ERRFLG)
ELSE
CALL RD_HAVSIN(DABS(TVEL),STV,SFC,DTV,DFC,0,MU,ERRFLG)
ENDIF
FRICTIONFORCE = MU * DABS(FC(3))
IF ( TVEL == 0.0d0 ) THEN
FC(1) = 0.0D0
FC(2) = 0.0D0
ELSE
FC(1) = -FRICTIONFORCE*RVEL(1)/TVEL
FC(2) = -FRICTIONFORCE*RVEL(2)/TVEL
ENDIF
ENDIF
DO I = 1, 3
RESULT(I) = FC(I)
ENDDO
RETURN
END
#include "stdafx.h"
#include "DllFunc.h"
#include "math.h"
ContactForce_API void __cdecl contact_force(double time, double upar[], int npar, double pen, double rvel[], int jflag, int iflag, double* result)
{
using namespace rd_syscall;
// Parameter Information
// time : Simulation time of RD/Solver. (Input)
// upar : Parameters defined by user. (Input)
// npar : Number of user parameters. (Input)
// pen : Penetration at the contact point. (Input)
// rvel : Relative velocity between contact pair w.r.t. contact reference frame. (Input, Size : 3)
// jflag : When RD/Solver evaluates a Jacobian, the flag is true. (Input)
// iflag : When RD/Solver initializes arrays, the flag is true. (Input)
// result : Returned contact force vector. (Output, Size : 6)
// User Statement
double dpen, K, C, M1, M2, M3, STV, DTV, SFC, DFC;
double SFO, DFO, TVEL, FC[3], ZERO, MU;
int i,ERRFLG;
dpen = rvel[2];
K = upar[0];
C = upar[1];
M1 = upar[2];
M2 = upar[3];
M3 = upar[4];
STV = upar[5];
DTV = upar[6];
SFC = upar[7];
DFC = upar[8];
ZERO = 0.0;
//---- COMPUTE CONTACT NORMAL FORCE
if ( fabs(pen) < 1.0e-12 ) {
SFO = 0.0;
}
else {
SFO = K * pow(fabs(pen),M1);
}
if ( fabs(dpen) < 1.0e-12 ) {
DFO = 0.0;
}
else {
DFO = C * pow(fabs(pen),M3) * pow(fabs(dpen),M2);
}
if ( dpen > ZERO ) DFO = -DFO;
FC[2] = SFO + DFO;
//---- COMPUTE FRICTION FORCE
TVEL = sqrt(rvel[0]*rvel[0] + rvel[1]*rvel[1]);
FC[0] = 0.0;
FC[1] = 0.0;
if ( TVEL >= 1.0e-12 ) {
if ( TVEL < STV ) {
rd_havsin(TVEL,-STV,-SFC,STV,SFC,0,&MU,&ERRFLG);
}
else {
rd_havsin(TVEL,STV,SFC,DTV,DFC,0,&MU,&ERRFLG);
}
double frictionForce = MU * fabs(FC[2]);
if (TVEL == 0.0) {
FC[0] = 0.0;
FC[1] = 0.0;
}
else {
FC[0] = -frictionForce*rvel[0]/TVEL;
FC[1] = -frictionForce*rvel[1]/TVEL;
}
}
for(i=0;i<3;i++) result[i] = FC[i];
}