2.5. Subroutines for Customizing Contact Interfacial Behavior

The subroutines enable you to:

2.5.1. Subroutine usercnprop (Programming Your Own Contact Properties)

This subroutine applies to the CONTA17x contact elements.

*deck,usercnprop                      USERDISTRIB
      subroutine usercnprop (ndim,coor,nkeyopt,keyopt,nrl,rlconst,
     x  nintIn,intIn,nrealIn,realIn,kupdhis,localr,nuval,nintp,usvr,
     x  ncomp,stress,strain0,strain,kstat,mu,kcnprop,cnprop,keyerr)
c
c *** primary function:   Allow users to define their own contact properties
c                         in real constant table
c                         This logic is accessed with real constant defined
c                         by table name: %_CNPROP% 
c                         (e.g. rmod,cid,kcnprop,%_CNPROP%)
c
c *** Notice - This file contains ANSYS Confidential information ***
c
c
c         Copyright ANSYS.  All Rights Reserved.
c         *** ansys, inc.
c
c  input arguments:
c     variable (type,sze,intent)    description
c
c     elem     (int,sc,in)         - element number
c     intpt    (int,sc,in)         - element integration point number
c     ndim     (int,sc,in)         - number of dimensions of the problem
c                                    = 2 2D
c                                    = 3 3D
c     nintp    (int,sc,in)         - the total number of integration points of
c                                    an element to be used with this routine
c     nuval    (int)               - number of additional state variables per
c                                    integration point
c         note:  nuval x nintp = nstv(on nsvr command); cannot exceed 840!
c
c     intIn    (int,ar(*),in)      - integer variables passed in
c     intIn    (int,ar(*),in)      - integer variables passed in
c                                    intIn(1) = element number
c                                    intIn(2) = element integration point number
c                                    intIn(3) = material reference number
c                                    intIn(4) = element type ID number (absolute value)
c                                               > 0 for CONTA171-CONTA177
c                                               < 0 for CONTA178
c                                    intIn(5) = real constant ID number
c                                    intIn(6) = associated contact nodal number
c                                    intIn(7) = contact indicator
c                                               0: intersection is found
c                                               otherwise: no intersection
c                                    intIn(8) = target element number
c                                    intIn(9) = flag for forcing sliding 
c                                               frictional case
c                                               0 : not forcing
c                                               1 : forcing (Slip direction is 
c                                                   defined through CMROT command)
c                                    intIn(10) = 1 first pass through 
c                                                  (1st iteration)
c                                                  (useful for initializing state 
c                                                  variables to a non-zero value)
c                                              = 2 first pass through key of
c                                                  a restart
c                                              = 3 first pass through key of
c                                                  a rezoning
c                                    intIn(11) = current load step number
c                                    intIn(12) = current substep number
c                                    intIn(13) = current equilibrium iteration 
c                                                number
c                                    intIn(14) = flag for using unsymmetric 
c                                                matrices (nropt,unsym)
c                                                0 : symmetric
c                                                1 : unsymmetric
c                                    intIn(15) = Linear perturbation flag
c                                                0 : a general load step
c                                                1 : a linear perturbation step
c                                    intIn(16) = key to indicate output pass
c                                                0 : not a output pass
c                                                1 : output pass
c                                    intIn(17) = key to indicate if history-
c                                                dependent variables 
c                                                (user defined) need to be 
c                                                updated after the substep has 
c                                                converged
c                                                1 : update (converged)
c                                                0 : do not update (not converged)
c                                    intIn(18) = key to indicate transient effects
c                                                1 : transient is active
c                                                0 : transient is not active
c                                    intIn(19) = large deformation key [nlgeom cmd]
c                                                1 : on
c                                                0 : off
c                                    intIn(20) = analysis type (derived from 
c                                                antype cmd)
c                                                0 : a static analysis
c                                                1 : a buckling analysis
c                                                2 : a modal analysis
c                                                3 : a harmonic analysis
c                                                4 : a transient analysis
c                                                7 : a substructure analysis
c                                                8 : a spectrum analysis
c                                    intIn(21) = key for displacement & force 
c                                                convergence
c                                                1 : converged
c                                                0 : not converged
c     realIn   (dp,ar(*),in)       - real variables passed in
c                                    realIn(1) = contact element length
c                                    realIn(2) = contact element depth
c                                    realIn(3) = area associated with the contact 
c                                                detection point
c                                    realIn(4) = pinball radius
c                                    realIn(5) = un-scaled normal penalty stiffness
c                                    realIn(6) = time (or frequency for a harmonic 
c                                                analysis) at the beginning of this 
c                                                load step
c                                    realIn(7) = time (or frequency for a harmonic 
c                                                analysis) at the end of this load step 
c                                    realIn(8) = current time value (or frequency value
c                                                for a harmonic analysis)
c                                    realIn(9) = time increment (or frequency increment 
c                                                for a harmonic analysis) over this
c                                                substep
c                                    realIn(10)= temperature offset from absolute
c                                                 zero
c                                    realIn(11)= geometric penetration/gap 
c                                                (current substep) 
c                                                > 0 : gap
c                                                < 0 : penetration
c                                    realIn(12)= time increment scaling factor to
c                                                be used for structural transient 
c                                                dynamics
c     nkeyopt  (int,sc,in)         - number of key options
c     keyopt   (int,ar(nkeyopt),in)- array containing key options
c                                    keyopt(1) : Select degree of freedom
c                                    keyopt(2) : Contact algorithm
c                                    ... so on (see ANSYS documentation)
c     nrl      (int,sc,in)         - number of real constants
c     rlconst  (dp,ar(nrl),in)     - array containing real constants 
c                                    Elements CONTA171 to CONTA177
c                                    rlconst(1) : R1
c                                    rlconst(2) : R2
c                                    rlconst(3) : FKN
c                                    rlconst(4) : FTOLN
c                                    ... so on (see ANSYS documentation)
c                                    Element CONTA178
c                                    rlconst(1) : FKN
c                                    rlconst(2) : GAP
c                                    ... so on (see ANSYS documentation)
c     
c     kcnprop  (int,sc,in)         - the position of constant in the real set
c                                    
c                                    (see ANSYS contact element manual)
c     ncomp    (int,sc,in)         - number of stress/force component
c                                    = 9 for CONTA171-CONTA177
c                                    = 7 for CONTA178

c     stress   (dp,ar(ncomp),in)   - stress components at the beginning of 
c                                    the current iteaation/substep.
c                                    stress(1) = frictional stress in direction 1
c                                    stress(2) = frictional stress in direction 2 
c                                                (3D only)
c                                    stress(3) = contact normal pressure
c                                    > 0 : compression
c                                    < 0 : tension
c                                    the above contact traction must be defined in 
c                                    a local coordinate system (see localr)
c                                    stress(4) = heat flux (per area) 
c                                                flowing into contact
c                                    stress(5) = heat flux (per area) 
c                                                flowing into target
c                                    < 0 heat flowing into a surface
c                                    > 0 heat flowing out of a surface
c                                    stress(6) = electrical current density 
c                                                (or pore fluid flux density) 
c                                                (per area) flowing into contact
c                                    stress(7) = electrical current density 
c                                                (or pore fluid flux density) 
c                                                (per area) flowing into target
c                                    > 0 current flowing out of a surface
c                                    < 0 current flowing into a surface
c                                    stress(8) = diffusion flux density 
c                                                (per area) flowing into contact
c                                    stress(9) = diffusion flux density 
c                                                (per area) flowing into target
c                                    > 0 flux flowing out of a surface
c                                    < 0 flux flowing into a surface
c     strain0  (dp,ar(ncomp),in) - strain components in the end of the previous 
c                                    substep
c                                    (see strain for each component definition)
c     strain  (dp,ar(ncomp),in)  - current strain components
c                                    strain(1) = slip increment in direction 1
c                                    strain(2) = slip increment in direction 2
c                                               (3D only)
c                                    strain(3) = contact normal gap/penetration
c                                    < 0 : gap
c                                    > 0 : penetration
c                                    strain(4) = temperature at the contact point
c                                                (from TEMP DOF or temperature load)
c                                    strain(5) = temperature at the target point
c                                                (only from TEMP DOF)
c                                    strain(6) = voltage (or pore pressure) 
c                                                at the contact point
c                                    strain(7) = voltage (or pore pressure) 
c                                                at the target point
c                                    strain(8) = concentrationat the contact point
c                                    strain(9) = concentrationat the target point
c     kstat    (int,sc,in)         - contact status at the end of the previous 
c                                    substep
c                                    3 : stick
c                                    2 : sliding
c                                    1 : open contact (near)
c                                    0 : open contact (far)
c     mu       (dp,sc,in)          - The frictional coef at the end of previous 
c                                    substep
c     coor     (dp,ar(6),in)       - Coordinates of the contact detection point
c                                    coor(1) current x
c                                    coor(2) current y
c                                    coor(3) current z
c                                    coor(4) initial x
c                                    coor(5) initial y
c                                    coor(6) initial z
c     localr   (dp,ar(3,3),in)     - the direction cosines of the local surface 
c                                    coordinate system at contact detection
c                                    localr(1,1),localr(1,2),localr(1,3) in slip 
c                                                                        direction 1
c                                    localr(2,1),localr(2,2),localr(2,3) in slip 
c                                                                        direction 2
c                                    localr(3,1),localr(3,2),localr(3,3) in normal 
c                                                                        direction
c                                    
c     usvr     (dp,ar(nuval,nintp),inout)- additional state variables from
c                                    previous equilibrium iteration (saved
c                                    if the nsvr command is used)
c     kupdhis  (int,sc,in)         - key to indicate if history-dependent
c                                    variables (user defined) need to be 
c                                    updated after the substep has converged
c                                    1 : update (converged)
c                                    0 : do not update (not converged)
c
c  output arguments:
c     variable (type,sze,intent)    description

c     cnprop   (dp,ar(5),out)      - user defined real constant value and 
c                                    derivatives w.r.t. kcnprop position
c                                    cnprop(1) = user defined real constant value
c                                    (e.g. kcnprop = 3 for normal contact 
c                                      stiffness FKN.
c                                      positive as scaling factor; 
c                                      negative value as the absolute value)
c                                    cnprop(2) = derivative of the real constant 
c                                                w.r.t. geometric penetration/gap
c                                    cnprop(3) = derivative of the real constant 
c                                                w.r.t. contact normal pressure
c                                    cnprop(4) = derivative of the real constant 
c                                                w.r.t. temperature at contact
c                                    cnprop(5) = derivative of the real constant 
c                                                w.r.t. temperature at target
c     usvr     (dp,ar(nuval,nintp),inout)- updated additional state variables
c                                    They are passed in as the values at the 
c                                    beginning of this substep.
c                                    They are updated to be the values at the 
c                                    end of this substep
c                                    Use NSVR command to size usvr array and
c                                    set nuval to same value as number of
c                                    variables on NSVR commands
c                                    Use userou.F to save these values
c                                    on NMISC record for output purposes.
c                                    The number of user defined output items on 
c                                    NMISC should be equal or less than NSTV
c                                    on nsvr command). It cannot exceed 120.
c
c      keyerr  (int,sc,inout)      - key to indicate if there is any element 
c                                    formulation error.
c                                    The error could be caused by too
c                                    large incremental step, illegal model.
c                                    = 0 no error (present value before calling)
c                                    = 1 some error happens. ANSYS will
c                                    decide to stop the analysis or cutback
c                                    the substep (bi-section) based on other
c                                    user input and information at higher
c                                    level.
c

2.5.2. Subroutine userfric (Writing Your Own Friction Laws)

This subroutine applies to the CONTA17x contact elements.

*deck,userfric                      USERDISTRIB
      subroutine userfric (elem,mat,intpt,nkeyopt,keyopt,nrl,rlconst,
     x ncomp,npropu,uprop,kfirst,kfsteq,kn,kt,elen,kstat,timval,
     x timinc,tcont,ttarg,toffst,dslip,slip,pres,tau,dt,usvr,
     x fdiss,elener,kupdhis,mu,dtdp,coor)
c
c *** primary function:   Allow users to write their own friction laws.
c                         This logic is accessed with tb,fric with tbopt=user.
c                         The below demonstration logic is the same as using
c                         tb,fric for isotropic Coulomb friction.
c                         Other friction laws may require more general 
c                         definition of friction forces.
c *** secondary function: demonstrate the use of user-written friction laws
c                         in this routine:
c                         a. update history variables
c                         b. compute consistent tangent matrix
c
c *** Notice - This file contains ANSYS Confidential information ***
c
c
c         Copyright ANSYS.  All Rights Reserved.
c         *** ansys, inc.
c
c  input arguments:
c     variable (type,sze,intent)    description
c
c     elem     (int,sc,in)         - element number (label)
c     mat      (int,sc,in)         - material reference number
c     intpt    (int,sc,in)         - element integration point number
c     nkeyopt  (int,sc,in)         - number of key options
c     keyopt   (int,ar(nkeyopt),in)- array containing key options
c                                    keyopt(1) : Select degree of freedom
c                                    keyopt(2) : Contact algorithm
c                                    ... so on (see ANSYS documentation)
c     nrl      (int,sc,in)         - number of real constants
c     rlconst  (dp,ar(nrl),in)     - array containing real constants
c                                    Elements CONTA171 to CONTA177
c                                    rlconst(1) : R1
c                                    rlconst(2) : R2
c                                    rlconst(3) : FKN
c                                    rlconst(4) : FTOLN
c                                    ... so on (see ANSYS documentation)
c                                    Element CONTA178
c                                    rlconst(1) : FKN
c                                    rlconst(2) : GAP
c                                    ... so on (see ANSYS documentation)
c     ncomp    (int,sc,in)         - no. of friction stress components (1 or 2)
c     npropu   (int,sc,in)         - no. of user-defined friction properties
c     uprop    (dp,ar(npropu),in)  - user-defined material properties
c     kfirst   (int,sc,in)         - 1 if first time through, 0 otherwise
c                                    (useful for initializing state variables
c                                    to a non-zero value)
c     kfsteq   (int,sc,in)         - 1 if first equilibrium iteration of a
c                                    substep, 0 otherwise
c     kn       (dp,sc,in)          - normal penalty stiffness
c     kt       (dp,sc,in)          - tangential penalty stiffness
c                                    (an initial guess is provided but
c                                     the user must pick a suitable
c                                     value that allows minimal tangential 
c                                     slip during sticking without 
c                                     adversely affecting the convergence;
c                                     a possible choice could be kt=mu*kn).
c                                     For Lagrange multiplier method (keyopt(2)=4
c                                     use a small kt (several orders of magnitude 
c                                     smaller than mu*pres).
c     elen     (dp,sc,in)          - length of contact element
c     kstat    (int,sc,inout)      - contact status
c                                    3 : stick
c                                    2 : sliding
c                                    1 : open contact (near)
c                                    0 : open contact (far)  
c     timval   (dp,sc,in)          - current time value
c     timinc   (dp,sc,in)          - time increment over this substep
c     tcont    (dp,sc,in)          - contact surface temperature
c                                    (from temperature DOF or temperature load)
c     ttarg    (dp,sc,in)          - target surface temperature 
c                                    (only from temperature DOF)
c     toffst   (dp,sc,in)          - temperature offset from absolute zero
c     dslip    (dp,ar(ncomp),in)   - slip increment (current substep)
c     slip     (dp,ar(ncomp),inout)- accumulated slip (previous substep)
c     pres     (dp,sc,in)          - normal pressure/force (current substep)
c                                    > 0 : compression
c                                    < 0 : tension 
c     tau      (dp,ar(ncomp),inout)- frictional stress (previous substep)
c                                    Lagrange multiplier contribution is added 
c                                    if keyopt(2)=4  
c     usvr     (dp,ar(nuval,nintp),inout)- additional state variables from
c                                    previous equilibrium iteration (saved
c                                    if the nsvr command is used)
c     kupdhis  (int,sc,in)         - key to indicate if history-dependent
c                                    variables (user defined) need to be 
c                                    updated after the substep has converged
c                                    1 : update (converged)
c                                    0 : do not update (not converged)
c     coor     (dp,ar(6),in)       - Coordinates of the contact detection point
c                                    coor(1) current x
c                                    coor(2) current y
c                                    coor(3) current z
c                                    coor(4) initial x
c                                    coor(5) initial y
c                                    coor(6) initial z
c
c  output arguments:
c     variable (type,sze,intent)    description
c
c     kstat    (int,sc,inout)      - updated contact status
c     mu       (dp,sc  inout)      - updated friction coefficient
c     slip     (dp,ar(ncomp),inout)- updated accumulated slip
c     tau      (dp,ar(ncomp),inout)- updated frictional stress
c     dt       (dp,ar(5,5),out)    - material tangent modulus
c                                    rows and columns of dt matrix are
c                                    associated to:
c                                    row 1 : frictional stress in direction 1 
c                                    row 2 : frictional stress in direction 2
c                                    row 3 : normal pressure 
c                                    row 4 : blank
c                                    row 5 : blank
c                                    col 1 : sliding in direction 1 
c                                    col 2 : sliding in direction 2
c                                    col 3 : normal gap 
c                                    col 4 : temperature at contact 
c                                    col 5 : temperature at targte
c                                    relevant components to be filled in are:
c                                    dt(1,1): d(tau1)/d(slip1)
c                                    dt(1,2): d(tau1)/d(slip2)
c                                    dt(1,3): d(tau1)/d(normal gap)
c                                    dt(2,1): d(tau2)/d(slip1)
c                                    dt(2,2): d(tau2)/d(slip2)
c                                    dt(2,3): d(tau2)/d(normal gap)
c                                    dt(3,3): d(pres)/d(normal gap)
c                                    dt(3,3) set to kn internally
c                                    dt(1,4) : d(tau1)/d(tcont)
c                                    dt(1,5) : d(tau1)/d(ttarg)
c                                    dt(2,4) : d(tau2)/d(tcont)
c                                    dt(2,5) : d(tau2)/d(ttarg)
c     dtdp     (dp,ar(ncomp),out)  - partial derivative of the frictional 
c                                    stress in direction 1/2 w.r.t. normal
c                                    pressure used in Lagrange multiplier 
c                                    method (keyopt(2)=3,4). 
c     usvr     (dp,ar(nuval,nintp),inout)- updated additional state variables
c                                    For example, mu value and absolute
c                                    accumulated slip could be output as follows:
c                                    usvr(1,intpt) : mu
c                                    usvr(2,intpt) : abs. acc. slip in dir1
c                                    usvr(3,intpt) : abs. acc. slip in dir2
c                                    Use NSVR command to size usvr array and
c                                    set nuval to same value as number of
c                                    variables on NSVR commands
c                                    Use userou.F to save these values
c                                    on NMISC record for output purposes.
c                                    The number of user defined output items on 
c                                    NMISC should be equal or less than NSTV
c                                    on nsvr command). It cannot exceed 120.
c
c     fdiss    (dp,sc,out)         - incremental frictional dissipation
c                                    per unit area
c     elener   (dp,sc,out)         - incremental elastic stored energy 
c                                    per unit area
c
c  fortran parameters (to be defined by the user):
c     variable (type)              description
c     nuval    (int)               - number of additional state variables per
c                                    integration point
c     nintp    (int)               - maximum number of integration points of
c                                    an element to be used with this routine
c                                    (14 is the maximum)
c         note:  nuval x nintp = nstv(on nsvr command); cannot exceed 840!
c
c  internal variables:
c     variable (type,sze)      description
c     dtfac    (dp,sc)         - temporary variable
c     taulim   (dp,sc)         - limit frictional stress
c     taueq    (dp,sc)         - equivalent frictional stress
c     dir1     (dp,sc)         - slip increment direction 1
c     dir2     (dp,sc)         - slip increment direction 2
c     dslipeq  (dp,sc)         - equivalent slip increment
c     oldt1    (dp,sc)         - frictional stress 1 from prev substep
c     oldt2    (dp,sc)         - frictional stress 2 from prev substep
c     err      (dp,ar(2))      - data array for diagnostic message
c

2.5.3. Subroutine userinter (Writing Your Own Contact Interactions)

This subroutine applies to the CONTA17x contact elements.

*deck,userinter                      USERDISTRIB
      subroutine userinter (ndim,coor,nkeyopt,keyopt,nrl,rlconst,
     x  npropu,uprop,nintIn,intIn,nrealIn,realIn,kupdhis,localr,
     x  nuval,nintp,usvr,ncomp,stress,strain0,strain,
     x  kstat,mu,dt,dtdp,kdamp,damp,fdiss,elener,keyerr,keycnv)
c
c *** primary function:   Allow users to write their own interaction behavior.
c                         This logic is accessed with tb,inter with tbopt=user.
c *** secondary function: demonstrate the use of user-written interface laws
c                         in this routine:
c                         a. update history variables
c                         b. compute consistent tangent matrix
c
c *** Notice - This file contains ANSYS Confidential information ***
c
c
c         *** Copyright ANSYS.  All Rights Reserved.
c         *** ansys, inc.
c
c  input arguments:
c     variable (type,sze,intent)    description
c
c     elem     (int,sc,in)         - element number
c     intpt    (int,sc,in)         - element integration point number
c     ndim     (int,sc,in)         - number of dimensions of the problem
c                                    = 2 2D
c                                    = 3 3D
c     nintp    (int,sc,in)         - the total number of integration points of
c                                    an element to be used with this routine
c     nuval    (int)               - number of additional state variables per
c                                    integration point
c         note:  nuval x nintp = nstv(on nsvr command); cannot exceed 840!
c
c     intIn    (int,ar(*),in)      - integer variables passed in
c                                    intIn(1) = element number
c                                    intIn(2) = element integration point number
c                                    intIn(3) = material reference number
c                                    intIn(4) = element type ID number (absolute value)
c                                               > 0 for CONTA171-CONTA177
c                                               < 0 for CONTA178
c                                    intIn(5) = real constant ID number
c                                    intIn(6) = associated contact nodal number
c                                    intIn(7) = contact indicator
c                                               0: intersection is found
c                                               otherwise: no intersection
c                                    intIn(8) = target element number
c                                    intIn(9) = flag for forcing sliding 
c                                               frictional case
c                                               0: - not forcing
c                                               1: - forcing (Slip direction is 
c                                                    defined through CMROT command)
c                                    intIn(10) = 1 first pass through 
c                                                  (1st iteration)
c                                                  (useful for initializing state 
c                                                   variables to a non-zero value)
c                                              = 2 first pass through key of 
c                                                  a restart
c                                              = 3 first pass through key of 
c                                                  a rezoning
c                                    intIn(11) = current load step number
c                                    intIn(12) = current substep number
c                                    intIn(13) = current equilibrium iteration number
c                                    intIn(14) = flag for using unsymmetric matrices 
c                                                (nropt,unsym)
c                                                0: - symmetric
c                                                1: - unsymmetric
c                                    intIn(15) = Linear perturbation flag
c                                                0: - a general load step
c                                                1: - a linear perturbation step
c                                    intIn(16) = key to indicate output pass
c                                                0: not a output pass
c                                                1: output pass
c                                    intIn(17) = key for displacement & force 
c                                                convergence
c                                                1: converged
c                                                0: not converged
c                                    intIn(18) = key to indicate transient effects 
c                                                1 : transient is active
c                                                0 : transient is not active
c                                    intIn(19) = large deformation key [nlgeom cmd]
c                                                1 : on
c                                                0 : off
c                                    intIn(20) = analysis type (derived from antype)
c                                                0 : a static analysis
c                                                1 : a buckling analysis
c                                                2 : a modal analysis
c                                                3 : a harmonic analysis
c                                                4 : a transient analysis
c                                                7 : a substructure analysis
c                                                8 : a spectrum analysis
c     realIn   (dp,ar(*),in)       - real variables passed in
c                                    realIn(1) = contact element length
c                                    realIn(2) = contact element depth
c                                    realIn(3) = area associated with the contact 
c                                                detection point
c                                    realIn(4) = pinball radius
c                                    realIn(5) = unscaled normal penalty stiffness
c                                    realIn(6) = time (or frequency for a harmonic 
c                                                analysis) at the beginning of this 
c                                                load step
c                                    realIn(7) = time (or frequency for a harmonic 
c                                                analysis) at the end of this load step 
c                                    realIn(8) = current time value (or frequency value
c                                                for a harmonic analysis)
c                                    realIn(9) = time increment (or frequency increment 
c                                                for a harmonic analysis) over this
c                                                substep
c                                    realIn(10) = temperature offset from absolute 
c                                                 zero
c                                    realIn(11) = geometric penetration/gap 
c                                                 (current substep) 
c                                                 > 0 : gap
c                                                 < 0 : penetration
c                                    realIn(12) = time increment scaling factor to
c                                                 be used for structural transient 
c                                                 dynamics
c                                    realIn(13) = convection coefficient (SFE command)
c                                    realIn(14) = bulk temp (SFE command)
c     nkeyopt  (int,sc,in)         - number of key options
c     keyopt   (int,ar(nkeyopt),in)- array containing key options
c                                    keyopt(1) : Select degree of freedom
c                                    keyopt(2) : Contact algorithm
c                                    ... so on (see ANSYS documentation)
c     nrl      (int,sc,in)         - number of real constants
c     rlconst  (dp,ar(nrl),in)     - array containing real constants 
c                                    Elements CONTA171 to CONTA177
c                                    rlconst(1) : R1
c                                    rlconst(2) : R2
c                                    rlconst(3) : FKN
c                                    rlconst(4) : FTOLN
c                                    ... so on (see ANSYS documentation)
c                                    Element CONTA178
c                                    rlconst(1) : FKN
c                                    rlconst(2) : GAP
c                                    ... so on (see ANSYS documentation)
c     ncomp    (int,sc,in)         - number of stress/force component
c                                    = 9 for CONTA171-CONTA177
c                                    = 7 for CONTA178
c     stress   (dp,ar(ncomp),inout)- stress components (current substep)
c                                    It is passed in as the stress at the beginning 
c                                    of the current substep. It is updated to be 
c                                    the stress at the end of this current substep
c                                    stress(1) = frictional stress in direction 1
c                                    stress(2) = frictional stress in direction 2 
c                                                (3D only)
c                                    stress(3) = contact normal pressure
c                                    > 0 : compression
c                                    < 0 : tension
c                                    the above contact traction must be defined in
c                                    a local coordinate system (see localr)
c                                    Lagrange multiplier contribution is added 
c                                    if keyopt(2)=3,4  
c                                    stress(4) = heat flux (per area) 
c                                                flowing into contact
c                                    stress(5) = heat flux (per area) 
c                                                flowing into target
c                                    < 0 heat flowing into a surface
c                                    > 0 heat flowing out of a surface
c                                    stress(6) = electrical current density 
c                                                (or pore fluid flux density) 
c                                                (per area) flowing into contact
c                                    stress(7) = electrical current density 
c                                                (or pore fluid flux density) 
c                                                (per area) flowing into target
c                                    > 0 current flowing out of a surface
c                                    < 0 current flowing into a surface
c                                    stress(8) = diffusion flux density 
c                                                (per area) flowing into contact
c                                    stress(9) = diffusion flux density 
c                                                (per area) flowing into target
c                                    > 0 flux flowing out of a surface
c                                    < 0 flux flowing into a surface
c     strain0  (dp,ar(ncomp),in)   - strain components in the end of the previous 
c                                    substep
c                                    (see strain for each component definition)
c     strain   (dp,ar(ncomp),in)   - current strain components
c                                    strain(1) = slip increment in direction 1
c                                    strain(2) = slip increment in direction 2 
c                                                (3D only)
c                                    strain(3) = contact normal gap/penetration
c                                    < 0 : gap
c                                    > 0 : penetration
c                                    strain(4) = temperature at the contact point
c                                               (from TEMP DOF or temperature load)
c                                    strain(5) = temperature at the target point
c                                               (only from TEMP DOF)
c                                    strain(6) = voltage (or pore pressure) 
c                                                at the contact point
c                                    strain(7) = voltage (or pore pressure) 
c                                                at the target point
c                                    strain(8) = concentrationat the contact point
c                                    strain(9) = concentrationat the target point
c                                    
c     kstat    (int,sc,inout)      - contact status (current substep)
c                                    It is passed in as the status at the 
c                                    beginning of the current substep. 
c                                    It is updated to be the status at the 
c                                    end of the current substep
c                                    3 : stick
c                                    2 : sliding
c                                    1 : open contact (near)
c                                    0 : open contact (far)  
c     coor     (dp,ar(6),in)       - Coordinates of the contact detection point
c                                    coor(1) current x
c                                    coor(2) current y
c                                    coor(3) current z
c                                    coor(4) initial x
c                                    coor(5) initial y
c                                    coor(6) initial z
c     localr   (dp,ar(3,3),in)     - the direction cosines of the local surface 
c                                    coordinate system at contact detection
c                                    localr(1,1),localr(1,2),localr(1,3) in 
c                                                          slip direction 1
c                                    localr(2,1),localr(2,2),localr(2,3) in
c                                                          slip direction 2
c                                    localr(3,1),localr(3,2),localr(3,3) in 
c                                                           normal direction
c                                    
c     npropu   (int,sc,in)         - number of user-defined interaction properties
c     uprop    (dp,ar(npropu),in)  - user-defined material properties
c
c     usvr     (dp,ar(nuval,nintp),inout)- additional state variables from
c                                     previous equilibrium iteration (saved
c                                     if the nsvr command is used)
c     kupdhis  (int,sc,in)         - key to indicate if history-dependent
c                                    variables (user defined) need to be 
c                                    updated after the substep has converged
c                                    1 : update (converged)
c                                    0 : do not update (not converged)
c
c  output arguments:
c     variable (type,sze,intent)   description
c
c     kstat    (int,sc,inout)      - updated contact status
c     stress   (dp,ar(ncomp),inout)- updated stress components
c                                    
c     dt   (dp,ar(ncomp,ncomp),out)- interface stiffness matrix:
c                                    dt(i,j) defines the partial derivative of 
c                                    the ith stress component at the current 
c                                    substep w.r.t. the jth component of the 
c                                    relative strain increment array.
c                                    If symmetric solver option used, ANSYS will 
c                                    symmetrize the matrix bu averaging the 
c                                    off-diagonal terms.
c                                    rows and columns of dt matrix are
c                                    associated to:
c                                    row 1 : frictional stress in direction 1 
c                                    row 2 : frictional stress in direction 2
c                                    row 3 : normal pressure 
c                                       > 0 : compression
c                                       < 0 : tension
c                                    row 4 : heat flux out the contact surface
c                                       < 0 heat flowing into contact
c                                       > 0 heat flowing out of target
c                                    row 5 : heat flux out the target surface
c                                       < 0 heat flowing into target
c                                       > 0 heat flowing out of target
c                                    row 6 : electrical current density
c                                            (or pore prssure) 
c                                            flowing out the contact surface
c                                       > 0 current flowing out of contact
c                                       < 0 current flowing into contact
c                                    row 7 : electrical current density
c                                            (or pore prssure) 
c                                       > 0 current flowing out of target
c                                       < 0 current flowing into target
c                                    row 8 : diffusion flux density
c                                            flowing out the contact surface
c                                       > 0 flux flowing out of contact
c                                       < 0 flux flowing into contact
c                                    row 9 : diffusion flux density
c                                       > 0 flux flowing out of target
c                                       < 0 flux flowing into target
c                                    col 1 : sliding in direction 1 
c                                    col 2 : sliding in direction 2
c                                    col 3 : normal gap 
c                                       < 0 : gap
c                                       > 0 : penetration
c                                    col 4 : temperature at the contact surface
c                                    col 5 : temperature at the target surface
c                                    col 6 : voltage at the contact surface
c                                    col 7 : voltage at the target surface
c                                    col 8 : concentration at the contact surface
c                                    col 9 : concentration at the target surface
c                                    relevant components to be filled in are:
c                                    dt(1,1): d(tau1)/d(slip1)
c                                    dt(1,2): d(tau1)/d(slip2)
c                                    dt(1,3): d(tau1)/d(normal gap)
c                                    dt(1,4): d(tau1)/d(tempC)
c                                    dt(1,5): d(tau1)/d(tempT)
c                                    dt(1,6): d(tau1)/d(voltC)
c                                    dt(1,7): d(tau1)/d(voltT)
c                                    dt(1,8): d(tau1)/d(concC)
c                                    dt(1,9): d(tau1)/d(concT)
c                                    dt(2,1): d(tau2)/d(slip1)
c                                    dt(2,2): d(tau2)/d(slip2)
c                                    dt(2,3): d(tau2)/d(normal gap)
c                                    ...
c                                    dt(3,1): d(pres)/d(slip 1)
c                                    dt(3,2): d(pres)/d(slip 2)
c                                    dt(3,3): d(pres)/d(normal gap)
c                                    ...
c                                    dt(4,1): d(fluxC)/d(slip 1)
c                                    dt(4,2): d(fluxC)/d(slip 2)
c                                    dt(4,3): d(fluxC)/d(normal gap)
c                                    dt(4,4): d(fluxC)/d(tempC)
c                                    dt(4,5): d(fluxC)/d(tempT)
c                                    dt(4,6): d(fluxC)/d(voltC)
c                                    dt(4,7): d(fluxC)/d(voltT)
c                                    dt(4,8): d(fluxC)/d(concC)
c                                    dt(4,9): d(fluxC)/d(concT)
c                                    ...
c                                    dt(5,4): d(fluxT)/d(tempC)
c                                    dt(5,5): d(fluxT)/d(tempT)
c                                    dt(5,6): d(fluxT)/d(voltC)
c                                    dt(5,7): d(fluxT)/d(voltT)
c                                    dt(5,8): d(fluxT)/d(concC)
c                                    dt(5,9): d(fluxT)/d(concT)
c                                    ...
c                                    dt(6,4): d(eleC)/d(tempC)
c                                    dt(6,5): d(eleC)/d(tempT)
c                                    dt(6,6): d(eleC)/d(voltC)
c                                    dt(6,7): d(eleC)/d(voltT)
c                                    dt(6,8): d(eleC)/d(concC)
c                                    dt(6,9): d(eleC)/d(concT)
c                                    ...
c                                    dt(7,4): d(eleT)/d(tempC)
c                                    dt(7,5): d(eleT)/d(tempT)
c                                    dt(7,6): d(eleT)/d(voltC)
c                                    dt(7,7): d(eleT)/d(voltT)
c                                    dt(7,8): d(eleT)/d(concC)
c                                    dt(7,9): d(eleT)/d(concT)
c                                    ...
c                                    dt(8,4): d(diffC)/d(tempC)
c                                    dt(8,5): d(diffC)/d(tempT)
c                                    dt(8,6): d(diffC)/d(voltC)
c                                    dt(8,7): d(diffC)/d(voltT)
c                                    dt(8,8): d(diffC)/d(concC)
c                                    dt(8,9): d(diffC)/d(concT)
c                                    ...
c                                    dt(9,4): d(diffT)/d(tempC)
c                                    dt(9,5): d(diffT)/d(tempT)
c                                    dt(9,6): d(diffT)/d(voltC)
c                                    dt(9,7): d(diffT)/d(voltT)
c                                    dt(9,8): d(diffT)/d(concC)
c                                    dt(9,9): d(diffT)/d(concT)
c     dtdp     (dp,ar(ncomp),out)  - partial derivative of the frictional stress
c                                    in direction 1,2 w.r.t. normal pressure 
c                                    used in Lagrange multiplier method 
c                                    (keyopt(2)=3,4). 
c     damp     (dp,ar(3,3),out)    - interface damping matrix (structure only)
c                                    it can be used only in Linear perturbation
c                                    modal analysis or transient analysis or
c                                    harmonic analysis in frequence domain.
c                                    damp(i,j) defines the partial derivative of 
c                                    the ith stress component at the current 
c                                    substep w.r.t. the jth component of the 
c                                    strain increment rate array.
c                                    rows and columns of dt matrix are
c                                    associated to:
c                                    row 1 : frictional stress in direction 1 
c                                    row 2 : frictional stress in direction 2
c                                    row 3 : normal pressure 
c                                    col 1 : sliding rate in direction 1 
c                                    col 2 : sliding rate in direction 2
c                                    col 3 : normal gap rate
c     kdamp    (int,sr,out)        - damping matrix index
c                                    0 : no damping matrix 
c                                    1 : taking damping matrix into account
c     usvr     (dp,ar(nuval,nintp),inout)- updated additional state variables
c                                    For example, mu value and absolute/relative
c                                    accumulated slip could be output as follows:
c                                    usvr(1,intpt) : mu
c                                    usvr(2,intpt) : abs. acc. slip in dir1
c                                    usvr(3,intpt) : abs. acc. slip in dir2
c                                    usvr(4,intpt) : acc. slip in dir1
c                                    usvr(5,intpt) : acc. slip in dir2
c                                    They are passed in as the values at the 
c                                    beginning of this substep. They are updated
c                                    to be the values at the end of this substep.
c                                    Use NSVR command to size usvr array and
c                                    set nuval to same value as number of
c                                    variables on NSVR commands
c                                    Use userou.F to save these values
c                                    on NMISC record for output purposes.
c                                    The number of user defined output items on 
c                                    NMISC should be equal or less than NSTV
c                                    on nsvr command). It cannot exceed 120.
c
c     mu       (dp,sc,inout)       - The current frictional coefficient
c     fdiss    (dp,sc,out)         - incremental frictional dissipation
c                                    per unit area
c     elener   (dp,sc,inout)       - Total elastic stored energy 
c                                    per unit area. 
c                                    Previous converged value is passed in and
c                                    current total should be the output 
c                                    
c     keyerr (int,sc,out)          - key to indicate if there is any element 
c                                    formulation error, like 
c                                    contact status changes abruptly,
c                                    too much penetration.
c                                    The error could be caused by too
c                                    large incremental step, illegal model.
c                                    = 0 no error (preset value before calling)
c                                    = 1 some error happens. ANSYS will
c                                    decide to stop the analysis or cutback
c                                    the substep (bi-section) based on other
c                                    user input and information at higher
c                                    level.
c     keycnv (int,sc,inout)        - key to flag if this element satisfies
c                                    the user defined element convergence
c                                    criterion. 
c                                    = 1, yes, the criterion is satisfied
c                                      or don't have any criterion at all
c                                      it is preset value before calling
c                                    = 0, no, the element doesn't satisfy
c                                      element convergence criterion. If
c                                      this is the case, the iteration will
c                                      not converge even when both force
c                                      and displacement converge 
c
c  internal variables:
c     variable (type,sze)      description
c

2.5.4. Subroutine userwear (Writing Your Own Wear Law)

This subroutine applies to CONTA17x contact elements.

*deck,userwear                      USERDISTRIB
      subroutine userwear(WearInc,WearDir,TotWearOld,strain,
     x     stress,temperature,dtime,YieldStress,nTbprop,Tbprop,
     x     coor,kstat,elem,intpt,ndim,localr,intIn,realIn,usvr,
     x     keyopt,rlconst)
c *** Primary Function:   Calculates the Wear Increment and (optionally) wear direction
c *** Notice - This file contains ANSYS Confidential information ***
c
c
c         *** Copyright ANSYS.  All Rights Reserved.
c         *** ansys, inc.
c input arguments:
c     variable (type,sze,intent)    description
c
c     ndim     (int,sc,in)         - number of dimensions of the problem
c                                    = 2 2D
c                                    = 3 3D
c     TotWearOld(dp,ar(ndim),in)   - Total Wear at the contact point at the previous substep
c     strain   (dp,ar(3),in)   - current strain components in contact surface coordinate system
c                                    strain(1) = slip increment in direction 1
c                                    strain(2) = slip increment in direction 2 
c                                                (3D only)
c                                    strain(3) = contact normal gap/penetration
c                                    < 0 : gap
c                                    > 0 : penetration
c     stress   (dp,ar(3),in)   - stress components  in contact surface coordinate system
c                                    stress(1) = frictional stress in direction 1
c                                    stress(2) = frictional stress in direction 2 
c                                                (3D only)
c                                    stress(3) = contact normal pressure
c                                    > 0 : compression
c                                    < 0 : tension
c     temperature  (dp,sc,in)      - temperature
c     dtime (dp,sc,in)             - time increment
c     YieldStress (dp,sc,in)       - Yield stress of underlying element (defined only for Plastic material-see doc for restrictions)
c     nTbprop (int,sc,in)          - Number of TBdata for Tb,Wear per field
c     Tbprop (dp,ar(nTbprop,in)    - TB data for the the Tb,Wear option at the given temperature       
c     coor     (dp,ar(6),in)       - Coordinates of the contact detection point
c                                     coor(1) current x
c                                     coor(2) current y
c                                     coor(3) current z
c                                     coor(4) initial x
c                                     coor(5) initial y
c                                     coor(6) initial z
c     kstat    (int,sc,in)          - contact status (current substep)
c                                     3 : stick
c                                     2 : sliding
c                                     1 : open contact (near)
c                                     0 : open contact (far) 
c     elem     (int,sc,in)         - element number
c     intpt    (int,sc,in)         - element integration point number
c     localr   (dp,ar(3,3),in)     - the direction cosines of the local surface 
c                                    coordinate system at contact detection
c                                    localr(1,1),localr(1,2),localr(1,3) in 
c                                                          slip direction 1
c                                    localr(2,1),localr(2,2),localr(2,3) in
c                                                          slip direction 2
c                                    localr(3,1),localr(3,2),localr(3,3) in 
c                                                           normal direction
c     intIn    (int,ar(*),in)      - integer variables passed in
c                                     intIn(1) = target element number the contact element
c                                                is in contat with (or is closest)
c                                     intIn(2) = Attached contact element number (if any)
c                                                to the target element in passed in intIn(1)
c                                     intIn(3) = number of additional state variables per
c                                                integration point (nuval)
c                                     intIn(4) = the total number of integration points of
c                                                an element to be used with this routine (nintp)
c                                     intIn(5) = key to indicate if history-dependent
c                                                variables (user defined) need to be 
c                                                updated after the substep has converged
c                                                1 : update (converged)
c                                                0 : do not update (not converged)
c     realIn   (dp,ar(*),in)       - real variables passed in
c                                    realIn(1) = contact element length
c                                    realIn(2) = contact element depth
c                                    realIn(3) = area associated with the contact 
c                                                detection point
c     
c     usvr     (dp,ar(nuval*nintp),inout)- additional state variables from
c                                     previous equilibrium iteration (saved
c                                     if the nsvr command is used  
c     keyopt   (int,ar(*),in)     - array containing key options for the element
c     rlconst  (dp,ar(*),in)      - array containing real constants for the element

c Output Arguments:
c     variable (type,sze,intent)   description
c
c     WearInc  (dp,sc)             - Increment in the Wear (magnitude)- User must define
c     WearDir  (dp,ar(ndim),inout) - Direction cosines in global coordinate system
c                                    in which wear increment will be applied- Optional  
c                                    default coming in -Contact normal direction