13.174. CONTA174 - 3-D 8-Node Surface-to-Surface Contact

Matrix or VectorGeometryShape Functions Integration Points
Stiffness and Stress Stiffness MatricesQuadIf KEYOPT(4) = 0 (has midside nodes) Equation 11–84, Equation 11–85, and Equation 11–86 2 x 2
TriangleIf KEYOPT(4) = 0 (has midside nodes) Equation 11–114 3

The following CONTA174 topics are available:

13.174.1. Introduction

CONTA174 is an 8-node element that is intended for general rigid-flexible and flexible-flexible contact analysis. In a general contact analysis, the area of contact between two (or more) bodies is generally not known in advance. CONTA174 is applicable to 3-D geometries. It may be applied for contact between solid bodies or shells.

13.174.2. Contact Kinematics

Contact Pair

In studying the contact between two bodies, the surface of one body is conventionally taken as a contact surface and the surface of the other body as a target surface. For rigid-flexible contact, the contact surface is associated with the deformable body; and the target surface must be the rigid surface. For flexible-flexible contact, both contact and target surfaces are associated with deformable bodies. The contact and target surfaces constitute a "Contact Pair".

The CONTA174 contact element is associated with the 3-D target segment elements (TARGE170) using a shared real constant set number. This element is located on the surface of 3-D solid, shell elements (called underlying element). It has the same geometric characteristics as the underlying elements. The contact surface can be either side or both sides of the shell or beam elements.

Location of Contact Detection

Figure 13.26:  Contact Detection Point Location at Gauss Point

Contact Detection Point Location at Gauss Point

CONTA174 is a surface-to-surface contact element. The contact detection points are the integration points and are located either at nodal points or Gauss points. The contact element is constrained against penetration into the target surface at its integration points. However, the target surface can, in principle, penetrate through into the contact surface. See Figure 13.26: Contact Detection Point Location at Gauss Point. CONTA174 uses Gauss integration points by default (Cescotto and Charlier([214]), Cescotto and Zhu([215])), which generally provides more accurate results than when using the nodes themselves as the integration points. A disadvantage of using nodal contact points is that, for a uniform pressure, the kinematically equivalent forces at the nodes are unrepresentative and indicate release at corners.

The surface projection method (based on Puso and Laursen ([383] and [384])) is available for the 3-D elements, CONTA173 and CONTA174. This option enforces contact constraints on an overlapping region of contact and target surfaces (see Figure 13.27: Surface-Projection-Based Contact) rather than on individual contact nodes or Gauss points. This method passes patch tests and generally offers smoother stress distributions (for both contact elements and underlying elements) than the other contact detection options. The contact forces vary smoothly during large sliding, and the forces do not jump when contact nodes slide off edges of target surfaces. The surface projection method also satisfies moment equilibrium for frictional or bonded contact. Other contact detection options may not satisfy moment equilibrium exactly when an offset exists along the contact normal direction under certain circumstances (for example, when using a penalty based formulation, when accounting for shell thickness effects, or when using a user-defined contact offset).

Figure 13.27:  Surface-Projection-Based Contact

Surface-Projection-Based Contact

Penetration Distance

The penetration distance is measured along the normal direction of the contact surface (at integration points) to the target surface (Cescotto and Charlier([215])). See Figure 13.28: Penetration Distance. It is uniquely defined even when the geometry of the target surface is not smooth. Such discontinuities may be due to physical corners on the target surface, or may be introduced by a numerical discretization process (e.g. finite elements). Based on the present way of calculating penetration distance, there is no restriction on the shape of the rigid target surface. Smoothing is not always necessary for the concave corner. For the convex corner, it is still recommended to smooth out the region of abrupt curvature changes (see Figure 13.29: Smoothing Convex Corner).

Figure 13.28:  Penetration Distance

Penetration Distance

Figure 13.29:  Smoothing Convex Corner

Smoothing Convex Corner

For the surface-projection-based contact method, the contact detection remains at contact nodal points and the contact normal is perpendicular to the contact surface. However, the contact penetration distance is computed over the overlapping region (see Figure 13.27: Surface-Projection-Based Contact) in an average sense.

Pinball Algorithm

The position and the motion of a contact element relative to its associated target surface determine the contact element status. The program monitors each contact element and assigns a status:

STAT = 0 Open far-field contact
STAT = 1 Open near-field contact
STAT = 2 Sliding contact
STAT = 3 Sticking contact

A contact element is considered to be in near-field contact when the element enters a pinball region, which is centered on the integration point of the contact element. The computational cost of searching for contact depends on the size of the pinball region. Far-field contact element calculations are simple and add few computational demands. The near-field calculations (for contact elements that are nearly or actually in contact) are slower and more complex. The most complex calculations occur the elements are in actual contact.

Setting a proper pinball region is useful to overcome spurious contact definitions if the target surface has several convex regions. The current default setting should be appropriate for most contact problems.

13.174.3. Frictional Model

Coulomb's Law

In the basic Coulomb friction model, two contacting surfaces can carry shear stresses. When the equivalent shear stress is less than a limit frictional stress (τlim), no motion occurs between the two surfaces. This state is known as sticking. The Coulomb friction model is defined as:

(13–259)

(13–260)

where:

τlim = limit frictional stress
μ = coefficient of friction for isotropic friction (input as MU using either TB command with Lab = FRIC or MP command; orthotropic friction defined below
P = contact normal pressure
b = contact cohesion (input as COHE on R command)

Once the equivalent frictional stress exceeds τlim, the contact and target surfaces will slide relative to each other. This state is known as sliding. The sticking/sliding calculations determine when a point transitions from sticking to sliding or vice versa. The contact cohesion provides sliding resistance even with zero normal pressure.

CONTA174 provides an option for defining a maximum equivalent frictional stress τmax (input as TAUMAX on RMORE command) so that, regardless of the magnitude of the contact pressure, sliding will occur if the magnitude of the equivalent frictional stress reaches this value.

Figure 13.30:  Friction Model

Friction Model

Contact elements offer two models for Coulomb friction: isotropic friction and orthotropic friction.

Isotropic Friction

The isotropic friction model uses a single coefficient of friction μiso based on the assumption of uniform stick-slip behavior in all directions. It is available with all 2-D and 3-D contact elements (CONTA171, CONTA172, CONTA173, CONTA174, CONTA175, CONTA176, CONTA177, and CONTA178).

Orthotropic Friction

The orthotropic friction model is based on two coefficients of friction, μ1 and μ2, to model different stick-slip behavior in different directions. Use orthotropic friction model in 3-D contact only (CONTA173, CONTA174, CONTA175, CONTA176, and CONTA177). The two coefficients are defined in two orthogonal sliding directions called the principal directions (see Element Reference for more details). The frictional stress in principal direction i=1,2 is given by:

(13–261)

By appropriately scaling the frictional stresses in principal directions the expressions for scaled limit frictional stress ( ) and scaled equivalent frictional stress ( ) for orthotropic friction can be written in a form similar to isotropic friction (Michalowski and Mroz([362])):

(13–262)

(13–263)

(13–264)

(13–265)

where:

= scaled frictional stress in direction i = 1,2
= scaled equivalent frictional stress
= scaled limit frictional stress
μeq = equivalent coefficient of friction for orthotropic friction
μ1, μ2 = coefficients of friction in first and second principal directions (input as MU1 and MU2 using TB command with Lab = FRIC)

While scaled frictional stresses are used for friction computations, actual frictional stresses are output after applying the inverse scaling in Equation 13–262.

The coefficient of friction (μ1 and μ2) can have dependence on time, temperature, normal pressure, sliding distance, or sliding relative velocity (defined as fields on the TBFIELD command). Suitable combinations of up to two fields can be used to define dependency, for example, temperature or temperature and sliding distance; see Contact Friction in the Material Reference for details.

Static and Dynamic Friction

CONTA174 provides the exponential friction model, which is used to smooth the transition between the static coefficient of friction and the dynamic coefficient of friction according to the formula (Benson and Hallquist([318])):

(13–266)

where:

v = slip rate
μs= Rf μd = static coefficient of friction
Rf = ratio of static and dynamic friction (input as FACT on RMORE command)
c = decay coefficient (input as DC on RMORE command)

Integration of Frictional Law

The integration of the frictional mode is similar to that of nonassociated theory of plasticity (see Rate-Independent Plasticity). In each substep that sliding friction occurs, an elastic predictor is computed in contact traction space. The predictor is modified with a radial return mapping function, providing both a small elastic deformation along with a sliding response as developed by Giannakopoulos([135]).

The flow rule giving the slip increment for orthotropic friction can be written as:

(13–267)

where:

dui = slip increment in principal direction i = 1, 2
λ = Lagrange multiplier for friction

By appropriately scaling the slip increment, it can be shown that the Lagrange multiplier is equal to the scaled equivalent slip increment:

(13–268)

(13–269)

and the direction of scaled slip increment is same as that of scaled frictional stress.

(13–270)

Thus, computations for orthotropic friction use the same framework as isotropic friction except for scaled slip increments and scaled frictional stresses which are converted to actual values for output.

User-defined Friction

For friction models that do not follow Coulomb’s law, you can write a USERFRIC subroutine. Refer to the Guide to User-Programmable Features in the Mechanical APDL Programmer's Reference for a detailed description on writing a USERFRIC subroutine. You can use it with any 2-D or 3-D contact element (CONTA171, CONTA172, CONTA173, CONTA174, CONTA175, CONTA176, CONTA177, and CONTA178) with penalty method for tangential contact (select KEYOPT(2) = 0, 1, or 3). Use the TB,FRIC command with TBOPT = USER to choose the user define friction option, and specify the friction properties on the TBDATA command.

Friction models involve nonlinear material behavior, so only experienced users who have a good understanding of the theory and finite element programming should attempt to write a USERFRIC subroutine.

Constitutive Relationship for Frictional Sliding

Regarding frictional sliding, the linearized constitutive relationship between tangential stresses and sliding increments can be expressed:

(13–271)

where:

= equivalent slip predictor over the current substep (Laursen and Simo[217])
= current equivalent slip rate
ti = the unit vector in the direction of slip (i = 1, 2)

can be obtained from the tabular input (define the friction coefficient to be a function of pressure via the TB,FRIC command).

can be obtained from the tabular input (define friction coefficient to be a function of sliding velocity via TB,FRIC command) or from the definition of Static and Dynamic Friction via real constants FACT and DC.

Algorithmic Symmetrization

Contact problems involving friction produce non-symmetric stiffness. Using an unsymmetric solver (NROPT,UNSYM) is more computationally expensive than a symmetric solver for each iteration. For this reason, a symmetrization algorithm developed by Laursen and Simo([217]) is used by which most frictional contact problems can be solved using solvers for symmetric systems. If frictional stresses have a substantial influence on the overall displacement field and the magnitude of the frictional stresses is highly solution dependent, any symmetric approximation to the stiffness matrix may provide a low rate of convergence. In such cases, the use of an unsymmetric stiffness matrix is more computationally efficient.

Forced Frictional Sliding Using Velocity Input

In a static analysis, you can model steady-state frictional sliding between two flexible bodies or between a flexible and a rigid body with different velocities. In this case the sliding velocities no longer follow the nodal displacements, and they are predefined through the CMROTATE command. This command sets the velocities on the nodes of the element component as an initial condition at the start of a load step. The program determines the sliding direction based on the direction of the sliding velocities. The frictional stresses then can be expressed by:

(13–272)

where is obtained from the CMROTATE command.

The linearized frictional stress with respect to the sliding velocity can be expressed in the following form (Bajer, Belsky, and Kung[398]):

(13–273)

In the complex eigenvalue extraction analysis using the QRDAMP or DAMP method, the above terms will be included in the damping matrix. The first term is, in general, negative in the case of a negative friction-sliding velocity gradient, and it can destabilize the system. Therefore, it is referred to as destabilizing squeal damping. The destabilizing squeal damping is proportional to the contact pressure, but it follows the sliding direction. The second term is proportional to the equivalent friction stress but orthogonal to the sliding direction. It may be needed to avoid false unstable modes and is, therefore, referred to as stabilizing squeal damping (Bajer, Belsky, and Kung[398]). To study the effects of squeal damping on the complex eigenvalue extraction analysis, the above equation is rewritten:

(13–274)

where:

cd = a scaling factor for destabilizing squeal damping (real constant FDMD when KEYOPT(16) = 0)
cs = a scaling factor for stabilizing squeal damping (real constant FDMS when KEYOPT(16) = 0)
= friction-sliding velocity gradient (real constant FDMD when KEYOPT(16) = 1)
= destabilizing damping coefficient (real constant FDMD when KEYOPT(16) = 2)
= stabilizing damping coefficient (real constant FDMS when KEYOPT(16) = 1 or 2)

13.174.4. Contact Algorithm

Four different contact algorithms are implemented in this element (selected by KEYOPT(2)).

  • Pure penalty method

  • Augmented Lagrangian method (Simo and Laursens([216]))

  • Pure Lagrange multiplier method (Bathe([2]))

  • Lagrange multiplier on contact normal and penalty on frictional direction

Pure Penalty Method

This method requires both contact normal and tangential stiffness. The main drawback is that the amount penetration between the two surfaces depends on this stiffness. Higher stiffness values decrease the amount of penetration but can lead to ill-conditioning of the global stiffness matrix and to convergence difficulties. Ideally, you want a high enough stiffness that contact penetration is acceptably small, but a low enough stiffness that the problem will be well-behaved in terms of convergence or matrix ill-conditioning.

The contact traction vector is:

(13–275)

where:

P = normal contact pressure
τ1 = frictional stress in direction 1
τ2 = frictional stress in direction 2

The contact pressure is:

(13–276)

where:

Kn = contact normal stiffness
un = contact gap size

The frictional stress for isotropic friction is obtained by Coulomb's law:

(13–277)

where:

Ks = tangential contact stiffness (input as FKT on R command)
Δui = slip increment in direction i over the current substep; In the case of sticking, Δui is the elastic slip increment, which represents the reversible tangential motion within the current substep.
= equivalent slip increment over the current substep
μiso = coefficient of friction
= frictional stress in direction i = 1,2 at the end of previous substep

For orthotropic friction, slip increment and frictional stress are scaled so that

(13–278)

where:

= scaled tangential contact stiffness in principal direction i = 1, 2
= slip increment in principal direction i = 1, 2 over the current substep
= scaled equivalent slip increment over the current substep
μeq = equivalent coefficient of friction
= scaled frictional stress in principal direction i = 1, 2 at the end of previous substep

For consistency between scaled friction stress and scaled slip increment, the scaled tangential contact stiffness in principal direction i = 1, 2 must be defined as:

(13–279)

Augmented Lagrangian Method

The augmented Lagrangian method is an iterative series of penalty updates to find the Lagrange multipliers (i.e., contact tractions). Compared to the penalty method, the augmented Lagrangian method usually leads to better conditioning and is less sensitive to the magnitude of the contact stiffness coefficient. However, in some analyses, the augmented Lagrangian method may require additional iterations, especially if the deformed mesh becomes excessively distorted.

The contact pressure is defined by:

(13–280)

where:

ε = compatibility tolerance (input as FTOLN on R command)
λi = Lagrange multiplier component at iteration i

The Lagrange multiplier component λi is computed locally (for each element) and iteratively.

Pure Lagrange Multiplier Method

The pure Lagrange multiplier method does not require contact stiffness. Instead it requires chattering control parameters. Theoretically, the pure Lagrange multiplier method enforces zero penetration when contact is closed and "zero slip" when sticking contact occurs. However the pure Lagrange multiplier method adds additional degrees of freedom to the model and requires additional iterations to stabilize contact conditions. This will increase the computational cost. This algorithm has chattering problems due to contact status changes between open and closed or between sliding and sticking. The other main drawback of the Lagrange multiplier method is overconstraint in the model. The model is overconstrained when a contact constraint condition at a node conflicts with a prescribed boundary condition on that degree of freedom (e.g., D command) at the same node. Overconstraints can lead to convergence difficulties and/or inaccurate results. The Lagrange multiplier method also introduces zero diagonal terms in the stiffness matrix, so that iterative solvers (e.g., PCG) can not be used.

The contact traction components (i.e., Lagrange multiplier parameters) become unknown DOFs for each element. The associated Newton-Raphson load vector is:

(13–281)

Lagrange Multiplier on Contact Normal and Penalty on Frictional Direction

In this method only the contact normal pressure is treated as a Lagrange multiplier. The tangential contact stresses are calculated based on the penalty method (see Equation 13–277).

This method allows only a very small amount of slip for a sticking contact condition. It overcomes chattering problems due to contact status change between sliding and sticking which often occurs in the pure Lagrange Multiplier method. Therefore this algorithm treats frictional sliding contact problems much better than the pure Lagrange method.

13.174.5. Viscous Damping

Viscous damping is primarily used to damp relative motions between the contact and target surfaces for open contact. It provides a certain resistance to reduce the risk of rigid body motion. The viscous damping traction is proportional but opposite to the relative velocities along contact normal or tangential directions:

(13–282)

(13–283)

(13–284)

where:

dn = damping coefficient in normal direction
dt = damping coefficient in tangential direction

The total damping force is calculated by the integration of damping traction over the contact elements.

The damping coefficient is specified as a function of the opening gap:

(13–285)

(13–286)

where:

upinb = pinball radius (real constant PINB)
dn0 = normal contact damping coefficient (real constant FDMN)
dt0 = tangential contact damping coefficient (real constant FDMT)

The contact damping energy can be obtained by integrating damping traction over the contacting area within the time increment:

(13–287)

where:

Ac = area domain of the contact surface
= time increment

13.174.6. Energy and Momentum Conserving Contact

To correctly model the physical interaction between contact and target surfaces in a transient dynamic analysis, the contact forces must maintain force and energy balance, and ensure proper transfer of linear momentum. This requires imposing additional constraints on relative velocities between contact and target surfaces (see Laursen and Chawla ([377]), and Armero and Petőcz ([378])).

Impact Constraints and Contact Forces

In ANSYS the penetration constraints and the relative velocity constraints between contact and target surfaces are collectively referred to as impact constraints. These constraints can be selected by setting KEYOPT(7) = 4 for any of the 2D or 3D contact elements and are valid for all types of contact interactions (flexible-to-flexible, flexible-to-rigid, and rigid-to-rigid) with and without friction.

An automatic time stepping scheme is used to predict the time of impact and adjust the size of the time increment to minimize penetration. The resulting penetration may not be within the penetration tolerance specified using real constant FTOLN. Rather, the allowable penetration is decided such that energy and momentum are conserved. When contact is detected, the relative velocity constraints are imposed using one of the four contact algorithms: pure penalty method, augmented Lagrangian method, pure Lagrange multiplier method, or Lagrange multiplier in contact normal and penalty in frictional direction method. In the case of rough contact (KEYOPT(12) = 1) the relative velocity constraint is imposed in the tangential direction also to prevent slip. In the case of standard contact (KEYOPT(12) = 0) with friction, the slip increment and frictional stress are computed by taking the relative velocity constraint into consideration.

For the pure penalty method, contact pressure P and friction stresses τi for isotropic friction are defined as:

(13–288)

where:

Kn = contact normal stiffness

un = contact gap size

= algorithmic contact gap size (based on the relative velocity constraint)

Pn-1 = normal contact pressure at the end of previous substep

and:

(13–289)

where:

Ks = tangential contact stiffness (input as FKT on R command)

= slip increment in direction i over the current substep

= equivalent slip increment over the current substep

= algorithmic slip increment in direction i over the current substep

= algorithmic equivalent slip increment over the current substep

μiso = coefficient of friction

= frictional stress in direction i = 1,2 at the end of previous substep

For other contact algorithms, the expressions for contact pressure and frictional stresses are defined in a similar manner as shown in Equation 13–280 and Equation 13–281 but with additional variables as shown above in Equation 13–288 and Equation 13–289.

Energy and Momentum Balance

Imposition of the impact constraints at Gauss points of contact elements ensures satisfaction of momentum and energy balance in a finite element sense. Since the impact constraints act only on the contact/target interface, energy balance is not enforced for the underlying finite elements used to model the interior of the contact and target bodies. Total energy at the contact/target interface is conserved for frictionless or rough contact when relative velocity constraints are satisfied exactly. If the relative velocity constraints are not satisfied to a tight tolerance there may be some loss of kinetic energy.

When friction is specified for contact elements, energy is conserved when the contact and target surfaces are not slipping (STICK) with respect to each other, and energy equal to the work done by frictional forces is dissipated when the contact and target surfaces are slipping (SLIP) with respect to each other.

Energy is also lost when numerical damping is used for the time integration scheme.

As per the classical theory of impact, exact conservation of energy during impact between rigid bodies is identified with elastic impact. It corresponds to a coefficient of restitution (COR) of 1. The impact constraints in ANSYS for impact between rigid bodies satisfy the conditions of elastic impact when the constraints are satisfied exactly and no numerical damping or friction is specified. Energy loss during impact between rigid bodies can be modeled by specifying a coefficient of restitution value (input as real constant COR) smaller than 1.

Time Integration Scheme

The impact constraints are formulated such that they can be used with both methods available for implicit transient dynamic analysis in ANSYS, the Newmark method and the HHT method. An important reason for using the impact constraints is that they make the time integration scheme numerically more stable without using large numerical damping. A small amount of numerical damping may still be needed to suppress high frequency noise.

13.174.7. Debonding

Debonding refers to separation of bonded contact (KEYOPT(12) = 2, 3, 4, 5 or 6). It is activated by associating a cohesive zone material model (input with TB,CZM) with contact elements. Debonding is available only for pure penalty method and augmented Lagrangian method (KEYOPT(2) = 0,1) with contact elements CONTA171, CONTA172, CONTA173, CONTA174, CONTA175, CONTA176, and CONTA177.

A cohesive zone material model is provided with bilinear behavior (Alfano and Crisfield([366])) for debonding. The model defines contact stresses as:

(13–290)

(13–291)

and

(13–292)

where:

P = normal contact stress (tension)
τ1 = tangential contact stress in direction 1
τ2 = tangential contact stress in direction 2
Kn = normal contact stiffness
Kt = tangential contact stiffness
un = contact gap
u1 = contact slip distance in direction 1
u2 = contact slip distance in direction 2
d = debonding parameter
direction 1 and direction 2 = principal directions in tangent plane

The debonding parameter is defined as:

(13–293)

with d = 0 for Δ 1 and 0 < d 1 for Δ > 1, and Δ and χ are defined below.

Debonding allows three modes of separation: mode I, mode II and mixed mode.

Mode I debonding is defined by setting

(13–294)

and

(13–295)

where:

= contact gap at the maximum normal contact traction (tension)
= contact gap at the completion of debonding (input on TBDATA command as C2 using TB,CZM)

Mode II debonding is defined by setting

(13–296)

(13–297)

and

(13–298)

where:

= equivalent tangential slip distance at the maximum equivalent tangential stress,
= equivalent tangential slip distance at the completion of debonding (input on TBDATA command as C4 using TB,CZM)

Mixed mode debonding is defined by setting

(13–299)

and

(13–300)

The constraint on χ that the ratio of the contact gap distances be same as the ratio of tangential slip distances is enforced automatically by appropriately scaling the contact stiffness values.

For mixed mode, debonding is complete when the energy criterion is satisfied:

(13–301)

with

(13–302)

(13–303)

(13–304)

(13–305)

where:

σmax = maximum normal contact stress (input on TBDATA command as C1 using TB,CZM)
τmax = maximum equivalent tangential contact stress (input on TBDATA command as C3 using TB,CZM)

Verification of satisfaction of energy criterion can be done during postprocessing of results.

The debonding modes are based on input data:

  1. Mode I for normal data (input on TBDATA command as C1, C2, and C5).

  2. Mode II for tangential data (input on TBDATA command as C3, C4, and C5).

  3. Mixed mode for normal and tangential data (input on TBDATA command as C1, C2, C3, C4, C5 and C6).

Artificial damping can be used to overcome convergence difficulties associated with debonding. It is activated by specifying the damping coefficient η (input on TBDATA command as C5 using TB,CZM).

Tangential slip under compressive normal contact stress for mixed mode debonding is controlled by appropriately setting the flag β (input on TBDATA command as C6 using TB,CZM). Settings on β are:

β = 0 (default) indicates no tangential slip
β = 1 indicates tangential slip is allowed

After debonding is completed the surface interaction is governed by standard contact constraints for normal and tangential directions. Frictional contact is used if friction is specified for contact elements.

13.174.8. Contact Surface Wear

Progressive loss of material from the surface of a solid body when in contact with another body (in any phase) is called wear. At a micro-scale, the process of wear is quite complex with various mechanical and chemical processes resulting in material failure and loss. At a continuum scale, wear is approximated by phenomenological models that relate various quantities at the contact surface with material loss (see Meng and Ludema [419]). ANSYS models wear using such a phenomenological approach. Material loss due to wear is approximated by repositioning the contact nodes at the contact surface. The new coordinates of the nodes are determined by a wear model. Since the contact nodes are moved to new positions, the contact variables (for example, contact pressure) change. The underlying continuum elements also experience a loss in material (and volume), thus simulating wear. After the nodes are moved to a new position, equilibrium may be lost and the program continues to iterate until equilibrium is achieved.

The following options are available for defining a wear model:

  • Archard wear model (see Archard [420])

  • User-defined wear model

Wear is activated using the TB command with Lab = WEAR for the material assigned to contact elements. The following contact elements support modeling wear: CONTA171, CONTA172, CONTA173, CONTA174, and CONTA175.

Archard Wear Model

The Archard wear model (see Archard [420]) is a popular sliding wear model with fairly good results for simulating wear. The original model proposed by Archard assumed that the rate of volume loss due to wear is linearly proportional to the contact pressure and sliding velocity at the contact surface. A generalized version of this model allowing proper law dependence on contact pressure and velocity is implemented in the program. The Archard model is defined using the TB command with TBOPT = ARCD and Lab = WEAR.

Wear is assumed to occur in the inward normal direction of the surface, taken as the direction opposite the contact normal direction. The rate of wear at a contact node, , is given by:

(13–306)

where:

K = wear coefficient
H = material hardness
P = contact pressure
vrel = the relative sliding velocity
m = pressure exponent
n = velocity exponent

In the above equation, K, H, m, and n are user-specified values from the TBDATA command. If a value of 1 is input for the fifth constant (C5) on the TBDATA command, the calculation is based on the nodal stresses of the underlying solid element instead of the contact pressure. The stress evaluated at the gauss point of the solid element is first copied to its nearest node. The contact element overlaid on the solid element to which the gauss point belongs uses this stress, σ, and the contact normal, n, to calculate the component of the traction vector in the contact normal direction: . This normal component of the traction vector is used instead of contact pressure in the wear rate equation.

If C5 on the TBDATA command is set to 10 or 11, the wear increment is averaged over the contact area of the contact pair such that the total volume lost due to wear is the same as the total volume lost due to wear when each node wears with a different amount (C5 = 0 or 1). The average wear increment is calculated as:

(13–307)

where A represents the contact area at the point and the summation is over all points in contact in a contact pair.

User-Defined Wear Model (USERWEAR)

If the Archard wear model does not fit your needs, you can define your own wear model by programming the subroutine USERWEAR. Most of the relevant contact results and properties are passed in the USERWEAR subroutine. You can define the wear increment and direction of the wear increment. The user-defined wear model is defined using the TB command with TBOPT = USER and Lab = WEAR.

13.174.9. Thermal/Structural Contact

Combined structural and thermal contact is specified if KEYOPT(1) = 1, which indicates that structural and thermal DOFs are active. Pure thermal contact is specified if KEYOPT(1) = 2. The thermal contact features (Zhu and Cescotto([281])) are:

Thermal Contact Conduction

(13–308)

where:

q = heat flux (heat flow rate per area)
Kc = thermal contact conductance coefficient (input as TCC on R command)
TT = temperature on target surface
TC = temperature on contact surface

Heat Convection

(13–309)

where:

hf = convection coefficient (input on SFE command with Lab = CONV and KVAL = 1)

Heat Radiation

(13–310)

where:

σ = Stefan-Boltzmann constant (input as SBCT on R command)
ε = emissivity (input using EMIS on MP command)
F = radiation view factor (input as RDVF on R command)
To = temperature offset (input as VALUE on TOFFST command)

Heat Generation Due to Frictional Sliding

(13–311)

where:

qc = amount of frictional dissipation on contact surface
qT = amount of frictional dissipation on target surface
Fw = weight factor for the distribution of heat between two contact and target surfaces (input as FWGT on R command)
Ff = fractional dissipated energy converted into heat (input on FHTG on R command)
= equivalent frictional stress
v = sliding rate


Note:  When KEYOPT(1) = 2, heat generation due to friction is ignored.


The frictional dissipation energy can be obtained by integrating dissipation density over the contacting area within the time increment

(13–312)

where:

Ac = area domain of the contact surface
= time increment

13.174.10. Electric Contact

Combined structural, thermal, and electric contact is specified if KEYOPT(1) = 3. Combined thermal and electric contact is specified if KEYOPT(1) = 4. Combined structural and electric contact is specified if KEYOPT(1) = 5. Pure electric contact is specified if KEYOPT(1) = 6. The electric contact features are:

Electric Current Conduction (KEYOPT(1) = 3 or 4)

(13–313)

where:

J = current density
σ/L = electric conductivity per unit length (input as ECC on R command)
VT = voltage on target surface
VC = voltage on contact surface

Electrostatic (KEYOPT(1) = 5 or 6)

(13–314)

where:

13.174.11. Magnetic Contact

The magnetic contact is specified if KEYOPT(1) = 7. Using the magnetic scalar potential approach, the 3-D magnetic flux across the contacting interface is defined by:

(13–315)

where:

ψn = magnetic flux
ϕt = magnetic potential at target surface (MAG degree of freedom)
ϕc = magnetic potential at contact surface (MAG degree of freedom)
CM = magnetic contact permeance coefficient
μo = free space permeability
A = contact area
= normal component of the "guess" magnetic field (See Equation 5–16)

The gap permeance is defined as the ratio of the magnetic flux in the gap to the total magnetic potential difference across the gap. The equation for gap permeance is:

(13–316)

where:

t = gap thickness

The magnetic contact permeance coefficient is defined as:

(13–317)

The above equations are only valid for 3-D analysis using the Magnetic Scalar Potential approach.

13.174.12. Pore Fluid Contact

Combined structural and pore fluid contact is specified if KEYOPT(1) = 8. Combined structural, pore fluid, and thermal contact is specified if KEYOPT(1) = 9.

Pore Fluid Flow Permeability (for closed contact and near-field contact):

(13–318)

(13–319)

where:

= pore fluid flux density out of contact surface
= pore fluid flux density out of target surface
= pore fluid contact permeability coefficient (input as PCC on R command)
= pore pressure on target surface
= pore pressure on contact surface

Pore Fluid Flow Due to Gap Rate Change (for near-field contact):

(13–320)

(13–321)

where:

= gap distance rate between the contact and target
= weight factor for the distribution of gap pore fluid flow between two contact and target surfaces (input as FPWT on R command)
= gap pore fluid flow participation factor (input as FPFT on R command)

This pore fluid flow is associated with removing or adding fluid for near field contact. For steady-state analysis, the changing rate for gap distance is zero so that there is no contribution from this gap-rate-dependent fluid flow.

Pore Fluid Flow (for far-field contact):

(13–322)

where:

= the pore pressure seepage coefficient (input as PSEE on R command)
= the ambient pore pressure (input as ABPP on R command)

Drainage-Only Flow (for near-field contact):

(13–323)

(13–324)

Drainage-only flow indicates that normal pore fluid flow occurs only from the contact surface to the environment when the contact surface pore pressures are positive. No fluid flow can enter the contact surface when the contact surface pore pressures are negative.

13.174.13. Diffusive Contact

Diffusive contact combined with other multiphysics fields is specified if KEYOPT(1) > 10.

Diffusion Flow for closed contact and near-field contact :

(13–325)

(13–326)

where:

qc = diffusion flux density out of contact surface
qt = diffusion flux density out of target surface
Kd = contact diffusivity coefficient (input as DCC on R command)
Cc = concentration on contact surface
Ct = concentration on target surface

Diffusion Flow for far-field contact :

(13–327)

where:

Kc = diffusive convection coefficient (input as DCON on R command)
C = ambient concentration (input as ABDC on R command)


Release 18.2 - © ANSYS, Inc. All rights reserved.