```prolog parameter( iTimStep = 1, * iTimGlb = 2, * iDTimCur = 3, * iTrackThic = 4, * nData = 4 ) C dimension jFlags(nFlags), rData(nData) C character*80 surfInt, surfSlv, surfMst C user coding to define stress, and, optionally, fluxSlv, fluxMst, state, sed, sfd, scd, spd, and svd C return end ``` # Variable to be defined stress(nDir, nBlock) On entry this array contains the stress defined in the local system at the proximity points during the previous time increment. It must be updated to the stress at the interface in the current time increment. # Variables that can be updated fluxSlv(nBlock) On entry this array contains the flux entering the slave surface during the previous time increment. It must be updated to the flux entering the slave surface during the current increment. fluxMst(nBlock) On entry this array contains the flux entering the master surface during the previous time increment. It must be updated to the flux entering the master surface during the current time increment. state(nStates,nNodState,nBlock) This array contains the user-defined solution-dependent state variables for the proximity points. The use of the state variables is applicable only for node-to-face contact. See “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide, for more information on the size of this array. This array will be passed in containing the values of these variables prior to the call to user subroutine VUINTERACTION. If any of the solution-dependent state variables is being used in conjunction with the interaction, it must be updated in this subroutine. These state variables need to be updated with care: outside the user subroutine these state variables are single-valued per slave node, but multiple proximity points may refer to the same slave node. Each proximity point may be passed into the user subroutine independently in a given increment, possibly on separate calls to the user subroutine; therefore, you may end up advancing the state variables for the associated node multiple times for a single time increment. To keep track of whether or not a node state is advanced, you may want to use one of the state variables exclusively for this purpose. You could set that selected state variable to the current increment number and update the state only if it is not already set to the current increment number. # sed(nBlock) On entry this array contains the elastic energy density at the proximity points at the beginning of the increment. It can be updated to contain the elastic energy density at the end of the current time increment. These values contribute to the output variable ALLSE and have no effect on other solution variables. The use of this variable is applicable only for node-to-face contact. # sfd(nBlock) This array can be updated to contain the increment in frictional dissipation at each proximity point (units of energy per unit area). These values contribute to the output variables SFDR and ALLFD and have no effect on other solution variables. The use of this variable is applicable only for node-to-face contact. # scd(nBlock) This array can be updated to contain the increment in creep dissipation at each proximity point (units of energy per unit area). These values contribute to the output variables SFDR and ALLCD and have no effect on other solution variables. The use of this variable is applicable only for node-to-face contact. # spd(nBlock) This array can be updated to contain the increment in plastic dissipation at each proximity point (units of energy per unit area). These values contribute to the output variables SFDR and ALLPD and have no effect on other solution variables. The use of this variable is applicable only for node-to-face contact. # svd(nBlock) This array can be updated to contain the increment in viscous dissipation at each proximity point (units of energy per unit area). These values contribute to the output variables SFDR and ALLVD and have no effect on other solution variables. The use of this variable is applicable only for node-to-face contact. # Variables passed in for information # nBlock Number of proximity points to be processed in this call to VUINTERACTION. # nBlockAnal 1 for analytical rigid master surface; nBlock otherwise. # nBlockEdge nBlock for edge type slave surface; 1 otherwise. # nNodState 1 for node-to-face and node-to-analytical rigid surface contact; not applicable for edge-to-edge contact. nNodSlv 1 for node-to-face and node-to-analytical rigid surface contact; 2 for edge-to-edge contact. nNodMst 1 for analytical rigid master surface; 2 for edge-type master surface; 4 for facet-type master surface. nDir Number of coordinate directions at the proximity points (equal to 3). nStates Number of user-defined state variables. nProps User-specified number of property values associated with this interaction model. nTemp 1 if the temperature is defined and 0 if the temperature is not defined. nFields Number of predefined field variables. jFlag(1) Step number. jFlag(2) Increment number. jFlag(3) 1 for node-to-face contact, 2 for edge-to-edge contact, and 3 for node-to-analytical rigid surface contact. rData(1) Value of step time. rData(2) Value of total time. rData(3) Current increment in time from $t = t_{curr} - \Delta t$ to $t = t_{curr}$ . rData(4) This variable contains the value of the tracking thickness specified for the surface interaction. surfInt User-specified surface interaction name, left justified. surfSlv Slave surface name, currently set to a blank. # surfMst Master surface name, currently set to a blank. # jSlvUid(nNodSlv,nBlock) This array lists the surface node numbers of the slave surface nodes associated with each proximity point. # jMstUid(nNodMst,nBlockAnal) This array lists the surface node numbers of the master surface nodes that make up the facet, edge, or analytical rigid surface associated with each proximity point. # props(nProps) User-specified vector of property values to define the interaction between the tracking surfaces. # penetration(nBlock) The relative position of the proximity points, with respect to the master surface, in the normal direction (a positive value indicates a penetration, and a negative value indicates a gap) during the current time increment. # drDisp(nDir,nBlock) An array containing the increments in relative positions of the proximity points with respect to the associated master surfaces during the current time increment. # dRot(2,2,nBlock) This argument is currently undefined. # stiffDef(nBlock) Values of the default penalty stiffnesses (stress per unit penetration, units of $\mathrm { F L } ^ { - 3 } )$ . # conductDef(nBlock) Values of the default penalty conductances (units of $\boldsymbol { J } \boldsymbol { \theta } ^ { - 1 } \boldsymbol { \mathrm { T } } ^ { - 1 } )$ . # dircos(nDir,nDir,nBlock) Direction cosines of the local surface coordinate system. The first array index corresponds to the components of the local directions, and the second array index corresponds to the local direction number. The first direction (dircos(1..nDir,1,...)) is the normal to the surface. The second direction (dircos(1..nDir,2,...)) is the first surface tangent. For a three-dimensional surface, the third direction (dircos(1..3,3,...)) is the second surface tangent. If the master surface is an analytical rigid surface, the numbers in dircos are valid only if the corresponding parts in penetration are valid (i.e., not equal to r\_MaxVal). # coordSlv(nDir,nNodSlv,nBlock) Array containing the nDir components of the current coordinates of the proximity points. # coordMst(nDir,nNodMst,nBlockAnal) Array containing the nDir components of the current coordinates of the nodes on the master surface. If the master surface is an analytical rigid surface, this array is passed in as a dummy array. # areaProx(nBlock) Contact area associated with a proximity point. The sum of the contact areas among all proximity points associated with a single slave node equals the surface area associated with that slave node (equal to 1 for node-based surface nodes). Therefore, the contact area at a proximity point depends on the number of other proximity points currently associated with the same slave node. A proximity point contributes a contact normal force to the associated slave node that is equal to stress(1,k) multiplied by areaProx(k). # shapeSlv(nNodSlv,nBlockEdge) For edge-to-edge contact this array contains the shape functions of the nodes of its slave edge, evaluated at the location of the contact point. If the contact is not edge-to-edge, this array is passed in as a dummy array. # shapeMst(nNodMst,nBlockAnal) For node-to-face and edge-to-edge contact this array contains the shape functions of the nodes of its master surface, evaluated at the location of the contact point. If the master surface is an analytical rigid surface, this array is passed in as a dummy array. # tempSlv(nBlock) Current temperature at the proximity points on the slave surface. # tempMst(nBlockAnal) Current temperature at the points on the master surface closest to the proximity points. # dTempSlv(nBlock) Increment in the temperature during the previous time increment at the proximity points on the slave surface. # dTempMst(nBlockAnal) Increment in the temperature during the previous time increment at the points on the master surface closest to the proximity points. # fieldSlv(nFields,nBlock) Current user-specified predefined field variables at the proximity points on the slave surface (initial values at the beginning of the analysis and current values during the analysis). # fieldMst(nFields,nBlockAnal) Current user-specified predefined field variables at the points on the master surface closest to the proximity points (initial values at the beginning of the analysis and current values during the analysis). dFieldSlv(nFields,nBlock) Increment in the user-specified predefined field variables during the previous time increment at the proximity points on the slave surface. dFieldMst(nFields,nBlockAnal) Increment in the user-specified predefined field variables during the previous time increment at the points on the master surface closest to the proximity points. # 1.2.22 VUMAT: User subroutine to define material behavior. # Product: Abaqus/Explicit WARNING: The use of this user subroutine generally requires considerable expertise. You are cautioned that the implementation of any realistic constitutive model requires extensive development and testing. Initial testing on a single-element model with prescribed traction loading is strongly recommended. The component ordering of the symmetric and nonsymmetric tensors for the three-dimensional case using C3D8R elements is different from the ordering specified in “Three-dimensional solid element library,” Section 28.1.4 of the Abaqus Analysis User’s Guide, and the ordering used in Abaqus/Standard. # References • “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide • \*USER MATERIAL # Overview # User subroutine VUMAT: • is used to define the mechanical constitutive behavior of a material; • will be called for blocks of material calculation points for which the material is defined in a user subroutine (“Material data definition,” Section 21.1.2 of the Abaqus Analysis User’s Guide); • can use and update solution-dependent state variables; • can use any field variables that are passed in; and • can be used in an adiabatic analysis, provided you define both the inelastic heat fraction and the specific heat for the appropriate material definitions and you store the temperatures and integrate them as user-defined state variables. # Component ordering in tensors The component ordering depends upon whether the tensor is symmetric or nonsymmetric. # Symmetric tensors For symmetric tensors such as the stress and strain tensors, there are ndir+nshr components, and the component order is given as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, a stress tensor contains ndir direct stress components and nshr shear stress components, which are passed in as

Component	2D Case	3D Case
1	$\sigma_{11}$	$\sigma_{11}$
2	$\sigma_{22}$	$\sigma_{22}$
3	$\sigma_{33}$	$\sigma_{33}$
4	$\sigma_{12}$	$\sigma_{12}$
5		$\sigma_{23}$
6		$\sigma_{31}$

The shear strain components in user subroutine VUMAT are stored as tensor components and not as engineering components; this is different from user subroutine UMAT in Abaqus/Standard, which uses engineering components. # Nonsymmetric tensors For nonsymmetric tensors there are ndir+2\*nshr components, and the component order is given as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, the deformation gradient is passed as

Component	2D Case	3D Case
1	$F_{11}$	$F_{11}$
2	$F_{22}$	$F_{22}$
3	$F_{33}$	$F_{33}$
4	$F_{12}$	$F_{12}$
5	$F_{21}$	$F_{23}$
6		$F_{31}$
7		$F_{21}$
8		$F_{32}$
9		$F_{13}$

# Initial calculations and checks In the data check phase of the analysis Abaqus/Explicit calls user subroutine VUMAT with a set of fictitious strains and a totalTime and stepTime both equal to 0.0. This is done as a check on your constitutive relation and to calculate the equivalent initial material properties, based upon which the initial elastic wave speeds are computed. # Defining local orientations All stresses, strains, stretches, and state variables are in the orientation of the local material axes. These local material axes form a basis system in which stress and strain components are stored. This represents a corotational coordinate system in which the basis system rotates with the material. If a user-specified coordinate system (“Orientations,” Section 2.2.5 of the Abaqus Analysis User’s Guide) is used, it defines the local material axes in the undeformed configuration. # Special considerations for various element types The use of user subroutine VUMAT requires special consideration for various element types. # Shell and plane stress elements You must define the stresses and internal state variables. In the case of shell or plane stress elements, NDIR=3 and NSHR=1; you must define strainInc(\*,3), the thickness strain increment. The internal energies can be defined if desired. If they are not defined, the energy balance provided by Abaqus/Explicit will not be meaningful. # Shell elements When VUMAT is used to define the material response of shell elements, Abaqus/Explicit cannot calculate a default value for the transverse shear stiffness of the element. Hence, you must define the element’s transverse shear stiffness. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User’s Guide, for guidelines on choosing this stiffness. # Beam elements For beam elements the stretch tensor and the deformation gradient tensor are not available. For beams in space you must define the thickness strains, strainInc(\*,2) and strainInc(\*,3). strainInc(\*,4) is the shear strain associated with twist. Thickness stresses, stressNew(\*,2) and stressNew(\*,3), are assumed to be zero, and any values you assign are ignored. # Pipe elements For pipe elements the stretch tensor and the deformation gradient tensor are not available. The axial strain, strainInc(\*,1), and the shear strain, strainInc(\*,4), associated with twist are provided along with the hoop stress, stressNew(\*,2). The hoop stress is predefined based on your pipe internal and external pressure load definitions (PE, PI, HPE, HPI, PENU, and PINU), and it should not be modified here. The thickness stress, stressNew(\*,3), is assumed to be zero and any value you assign is ignored. You must define the axial stress, stressNew(\*,1), and the shear stress, stressNew(\*,4). You must also define hoop strain, strainInc(\*,2), and the pipe thickness strain, strainInc(\*,3). # Deformation gradient The polar decomposition of the deformation gradient is written as , where and are the right and left symmetric stretch tensors, respectively. The constitutive model is defined in a corotational coordinate system in which the basis system rotates with the material. All stress and strain tensor quantities are defined with respect to the corotational basis system. The right stretch tensor, , is used. The relative spin tensor represents the spin (the antisymmetric part of the velocity gradient) defined with respect to the corotational basis system. # Special considerations for hyperelasticity Hyperelastic constitutive models in VUMAT should be defined in a corotational coordinate system in which the basis system rotates with the material. This is most effectively accomplished by formulating the hyperelastic constitutive model in terms of the stretch tensor, , instead of in terms of the deformation gradient, . Using the deformation gradient can present some difficulties because the deformation gradient includes the rotation tensor and the resulting stresses would need to be rotated back to the corotational basis. # Objective stress rates The Green-Naghdi stress rate is used when the mechanical behavior of the material is defined using user subroutine VUMAT. The stress rate obtained with user subroutine VUMAT may differ from that obtained with a built-in Abaqus material model. For example, most material models used with solid (continuum) elements in Abaqus/Explicit employ the Jaumann stress rate. This difference in the formulation will cause significant differences in the results only if finite rotation of a material point is accompanied by finite shear. For a discussion of the objective stress rates used in Abaqus, see “Stress rates,” Section 1.5.3 of the Abaqus Theory Guide. # Material point deletion Material points that satisfy a user-defined failure criterion can be deleted from the model (see “Userdefined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide). You must specify the state variable number controlling the element deletion flag when you allocate space for the solution-dependent state variables, as explained in “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide. The deletion state variable should be set to a value of one or zero in VUMAT. A value of one indicates that the material point is active, while a value of zero indicates that Abaqus/Explicit should delete the material point from the model by setting the stresses to zero. The structure of the block of material points passed to user subroutine VUMAT remains unchanged during the analysis; deleted material points are not removed from the block. Abaqus/Explicit will pass zero stresses and strain increments for all deleted material points. Once a material point has been flagged as deleted, it cannot be reactivated.