VFRIC

User subroutine to define frictional behavior for contact surfaces.

User subroutine VFRIC:

can be used to define the frictional behavior between contact pair surfaces;
can be used when the classical Coulomb friction model is too restrictive and a more complex definition of shear transmission between contacting surfaces is required;
must provide the entire definition of shear interaction between the contacting surfaces;
can use and update solution-dependent state variables;
cannot be used in conjunction with softened tangential surface behavior; and
cannot be used with the general contact algorithm.

This page discusses:

Terminology
Local Coordinate System
Frictional Forces
User Subroutine Interface
Variables to Be Defined
Variables That Can Be Updated
Variables Passed in for Information

Terminology

The use of user subroutine VFRIC requires familiarity with the following terminology.

Surface Node Numbers

The “surface node number” refers to the position of a particular node in the list of nodes on the surface. For example, there are nSecNod nodes on the secondary surface. Number $n, n = 1, 2...$ nSecNod, is the surface node number of the nth node in this list; jSecUid $(n)$ is the user-defined global number of this node. An Abaqus/Explicit model can be defined in terms of an assembly of part instances (see Assembly Definition). In such models a node number in jSecUid is an internally generated node number. If the original node number and part instance name are required, call the utility routine VGETPARTINFO (see Obtaining Part Information).

Contact Points

The nodes on the secondary surface that are in contact in the current time increment are defined as “contact points.” The number of contact points is passed into this subroutine as nContact. The array jConSecid(nContact) gives the surface node numbers for the contact points.

Local Coordinate System

A local coordinate system is defined for each contact point to facilitate specification of frictional forces and incremental slips. The local 1-direction for both two-dimensional and three-dimensional contact is tangential to the main surface, and it is defined by $t_{1} = d s / | d s |$ , where $d s$ is the incremental slip vector. The incremental slip vector used to define $t_{1}$ corresponds to the incremental slip in the current time increment for penalty contact and the predicted incremental slip for kinematic contact. The main surface normal direction, $n$ , is the local 2-direction for two-dimensional contact and the local 3-direction for three-dimensional contact. The local 2-direction for three-dimensional contact is given by $t_{2} = n \times t_{1}$ , which is also tangent to the main surface. The vectors are shown in Figure 1 and Figure 2. The direction cosines for $t_{1}$ and $n$ with respect to the global coordinate system are available in dirCosT1 and dirCosN, respectively. In the case of zero incremental slip ( $| d s | = 0$ ) we choose an arbitrary direction for $t_{1}$ that is orthogonal to the normal direction, $n$ .

Frictional Forces

You specify the frictional force, fTangential, at each contact point in local coordinates in this subroutine. The array fTangential is dimensioned such that only the tangential components can be specified. Any components of the frictional force that are not specified will remain equal to zero. For three-dimensional contact with isotropic friction, only the first component of the frictional force need be specified since the second component should be zero. A “stick force” at each contact point is provided in the array fStickForce to assist you in setting the appropriate frictional force values. The stick force is the force required to prevent additional “plastic” slipping. The stick force at each contact point is provided as a scalar value as it would act in the direction opposite to $t_{1}$ . The stick force is computed prior to calling user subroutine VFRIC by either the kinematic or the penalty contact algorithm. See Contact Constraint Enforcement Methods in Abaqus/Explicit for descriptions of the kinematic and penalty contact algorithms and the user interface for choosing between them. The first component of the frictional force should be in the range between zero and minus the stick force value. Typically, the stick force will be positive and the first component of the applied frictional force will be negative, opposing the incremental slip. Penalty contact includes an elastic slip regime due to finite penalty stiffness, so occasionally, during recovery of elastic slip, the stick force will be negative, indicating that it is appropriate for the first component of the frictional force to be positive (i.e., acting in the same direction as the incremental slip). A noisy or unstable solution is likely to result if the first component of fTangential is set outside of the range between zero and negative the value of the stick force.

After user subroutine VFRIC is called, frictional forces that oppose the forces specified at the contact points are distributed to the main nodes. For balanced main-secondary contact we then compute weighted averages of the frictional forces for both main-secondary orientations. These forces are directly applied if the penalty contact algorithm is being used. If the kinematic contact algorithm is being used, the frictional forces are converted to acceleration corrections by dividing by the nodal masses.

User Subroutine Interface

      subroutine vfric(
C Write only - 
     1 fTangential, 
C Read/Write - 
     2 statev,
C Read  only - 
     3 kStep, kInc, nContact, nFacNod, nSecNod, nMainNod,
     4 nFricDir, nDir, nStateVar, nProps, nTemp, nPred, numDefTfv,
     5 jSecUid,jMainUid, jConSecid, jConMainid, timStep, timGlb,
     6 dTimCur, surfInt, surfSec, surfMain, lContType,
     7 dSlipFric, fStickForce, fTangPrev, fNormal, frictionWork,
     8 shape, coordSec, coordMain, dirCosSl, dircosN, props,
     9 areaSec, tempSec, preDefSec, tempMain, preDefMain) 
C
      include `vaba_param.inc'
C
      character*80 surfInt, surfSec, surfMain
C
      dimension props(nProps), statev(nStateVar,nSecNod), 
     1 dSlipFric(nDir,nContact),
     2 fTangential(nFricDir,nContact),
     3 fTangPrev(nDir,nContact),
     4 fStickForce(nContact), areaSec(nSecNod),
     5 fNormal(nContact), shape(nFacNod,nContact),
     6 coordSec(nDir,nSecNod), coordMain(nDir,nMainNod), 
     7 dirCosSl(nDir,nContact), dircosN(nDir,nContact),
     8 jSecUid(nSecNod), jMainUid(nMainNod),
     9 jConSecid(nContact), jConMainid(nFacNod,nContact)
     1 tempSec(nContact), preDefSec(nContact,nPred),
     2 tempMain(numDefTfv), preDefMain(numDefTfv,nPred)

      user coding to define fTangential
      and, optionally, statev

      return
      end

Variables to Be Defined

fTangential(nFricDir, nContact): This array must be updated to the current values of the frictional force components for all contact points in the local tangent directions. See Figure 1 and Figure 2 for definition of the local coordinate system. This array will be zero (no friction force) until you reset it.

Variables That Can Be Updated

statev(nstateVar, nSecNod): This array contains the user-defined solution-dependent state variables for all the nodes on the secondary surface. You define the size of this array (see Frictional Behavior for more information). This array will be passed in containing the values of these variables prior to the call to user subroutine VFRIC. If any of the solution-dependent state variables is being used in conjunction with the friction behavior, it must be updated in this subroutine. The state variables are available even for secondary nodes that are not in contact. This may be useful when, for example, the state variables need to be reset for secondary nodes that are not in contact.

Variables Passed in for Information

kStep: Step number.
kInc: Increment number.
nContact: Number of contacting secondary nodes.
nFacNod: Number of nodes on each main surface facet (nFacNod is 2 for two-dimensional surfaces, nFacNod is 4 for three-dimensional surfaces). If the main surface is an analytical rigid surface, this variable is passed in as 0.
nSecNod: Number of secondary nodes.
nMainNod: Number of main surface nodes, if the main surface is made up of facets. If the main surface is an analytical rigid surface, this variable is passed in as 0.
nFricDir: Number of tangent directions at the contact points (nFricDir = nDir - 1).
nDir: Number of coordinate directions at the contact points. (In a three-dimensional model nDir will be two if the surfaces in the contact pair are two-dimensional analytical rigid surfaces or are formed by two-dimensional elements.)
nStateVar: Number of user-defined state variables.
nProps: User-specified number of property values associated with this friction model.
nTemp: 1 if the temperature is defined and 0 if the temperature is not defined.
nPred: Number of predefined field variables.
numDefTfv: Equal to nContact if the main surface is made up of facets. If the main surface is an analytical rigid surface, this variable is passed in as 1.
jSecUid(nSecNod): This array lists the user-defined global node numbers (or internal node numbers for models defined in terms of an assembly of part instances) of the nodes on the secondary surface.
jMainUid(nMainNod): This array lists the user-defined global node numbers (or internal node numbers for models defined in terms of an assembly of part instances) of the nodes on the main surface. If the main surface is an analytical rigid surface, this array is passed in as a dummy array.
jConSecid(nContact): This array lists the surface node numbers of the secondary surface nodes that are in contact.
jConMainid(nFacNod, nContact): This array lists the surface node numbers of the main surface nodes that make up the facet with which each contact point is in contact. If the main surface is an analytical rigid surface, this array is passed in as a dummy array.
timStep: Value of step time.
timGlb: Value of total time.
dtimCur: Current increment in time from $t = t_{c u r r} - Δ t$ to $t = t_{c u r r}$ .
surfInt: User-specified surface interaction name, left justified.
surfSec: Secondary surface name.
surfMain: Main surface name.
lContType: Contact type flag. This flag is set based on the type of constraint enforcement method (see Contact Constraint Enforcement Methods in Abaqus/Explicit) being used: 1 for kinematic contact and 2 for penalty contact. Stick conditions are satisfied exactly with the kinematic contact algorithm; they are satisfied only approximately (subject to an automatically chosen penalty stiffness value) with the penalty contact algorithm.
dSlipFric(nDir, nContact): This array contains the incremental frictional slip during the current time increment for each contact point in the current local coordinate system. These incremental slips correspond to tangential motion in the time increment from $t = t_{c u r r} - Δ t$ to $t = t_{c u r r}$ . For penalty contact this incremental slip is used to define the local coordinate system at each contact point (see Figure 1 and Figure 2) so that only the first component of dSlipFric can be nonzero in the local system. The contact points for kinematic contact are determined based on penetrations detected in the predicted configuration (at $t = t_{c u r r} + Δ t$ ), and the predicted incremental slip direction is used to define the local coordinate system at each contact point. If the slip direction changes between increments, dSlipFric may have a nonzero component in the local 2-direction and, if the surface is faceted and the contact point moves from one facet to another, in the local 3-direction.
fStickForce(nContact): This array contains the magnitude of frictional force required to enforce stick conditions at each contact point. For kinematic contact this force corresponds to no slip; for penalty contact this force depends on the previous frictional force, the value of the penalty stiffness, and the previous incremental slip. The penalty stiffness is assigned automatically. Occasionally, during recovery of elastic slip associated with the penalty method, the stick force will be assigned a negative value.
fTangPrev(nDir, nContact): This array contains the values of the frictional force components calculated in the previous increment but provided in the current local coordinate system (zero for nodes that were not in contact).
fNormal(nContact): This array contains the magnitude of the normal force for the contact points applied at the end of current time increment; i.e., at time $t = t_{c u r r}$ .
frictionWork: This variable contains the value of the total frictional dissipation in the entire model from the beginning of the analysis. The units are energy per unit area.
shape(nFacNod, nContact): For each contact point this array contains the shape functions of the nodes of its main surface facet, evaluated at the location of the contact point. If the main surface is an analytical rigid surface, this array is passed in as a dummy array.
coordSec(nDir, nSecNod): Array containing the nDir components of the current coordinates of the secondary nodes.
coordMain(nDir, nMainNod): Array containing the nDir components of the current coordinates of the main nodes. If the main surface is an analytical rigid surface, this array is passed in as a dummy array.
dirCosSl(nDir, nContact): Direction cosines of the incremental slip at the contact points.
dircosN(nDir, nContact): Direction cosines of the normals to the main surface at the contact points.
props(nProps): User-specified vector of property values to define the frictional behavior between the contacting surfaces.
areaSec(nSecNod): Area associated with the secondary nodes (equal to 1 for node-based surface nodes).
tempSec(nContact): Current temperature at the secondary nodes.
preDefSec(nContact,nPred): Current user-specified predefined field variables at the secondary nodes (initial values at the beginning of the analysis and current values during the analysis).
tempMain(numDefTfv): Current temperature at the nearest points on the main surface.
preDefMain(numDefTfv,nPred): Current user-specified predefined field variables at the nearest points on the main surface (initial values at the beginning of the analysis and current values during the analysis).