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 3D geometries. It may be applied for contact between solid bodies or shells.

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 3D target segment elements (TARGE170) using a shared real constant set number. This element is located on the surface of 3D 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.24: 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.24: Contact Detection Point Location at Gauss Point. CONTA174 uses Gauss integration points by default (Cescotto and Charlier([215]), Cescotto and Zhu([216])), 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 standard surface projection method (based on Puso and Laursen ([379] and [380])) is available for the 3D element, CONTA174. This option enforces contact constraints on an overlapping region of contact and target surfaces (see Figure 13.25: 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.25: Surface-Projection-Based Contact

The dual shape function projection method (based on Popp et al. [443]) is an extension of the standard surface projection method. It uses a dual shape function for the discretization of contact traction instead of using a displacement shape function. The dual shape function is derived based on the biorthogonality relationship from the standard displacement shape function. Therefore, this method reduces dependent terms from the adjacent nodes for contact constraints built on each contact node. The resulting gap function and slip increments are closer to the local geometric gap/slip at the contact nodes.

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([216])). See Figure 13.26: 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.27: Smoothing Convex Corner).

Figure 13.26: Penetration Distance

Figure 13.27: 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.25: 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:

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.

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–222)

(13–223)

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.28: 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 2D and 3D contact elements (CONTA172, CONTA174, CONTA175, CONTA177, and CONTA178).

Orthotropic Friction

The orthotropic friction model is based on two coefficients of friction, μ₁ and μ₂, to model different stick-slip behavior in different directions. Use orthotropic friction model in 3D contact only (CONTA174, CONTA175, 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–224)

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([358])):

(13–225)

(13–226)

(13–227)

(13–228)

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–225.

The coefficient of friction (μ₁ and μ₂) 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([317])):

(13–229)

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([136]).

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

(13–230)

where:

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

(13–231)

(13–232)

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

(13–233)

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 Programmer's Reference for a detailed description on writing a USERFRIC subroutine. You can use it with any 2D or 3D contact element (CONTA172, CONTA174, CONTA175, 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–234)

where:

= equivalent slip predictor over the current substep (Laursen and Simo[218])

= 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([218]) 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–235)

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[394]):

(13–236)

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[394]). To study the effects of squeal damping on the complex eigenvalue extraction analysis, the above equation is rewritten:

(13–237)

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)

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

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–238)

where:

The contact pressure is:

(13–239)

where:

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

(13–240)

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–241)

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–242)

As a special case, when contact friction is used in a steady-state rolling analysis, the frictional stress calculation is changed to a viscous formulation. The frictional stress for isotropic friction given in Equation 13–240 is modified as:

(13–243)

where is the slip velocity in direction i. All other symbols have the same meaning as in Equation 13–240.

Similarly, the orthotropic friction calculation given in Equation 13–241 is modified as:

(13–244)

where is the slip velocity in principle direction i.

Augmented Lagrangian Method

The augmented Lagrangian method is an iterative series of penalty updates to find the Lagrange multipliers (that is, 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–245)

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 (for example, 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 (for example, PCG) can not be used.

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

(13–246)

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–240).

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.

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–247)

(13–248)

(13–249)

where:

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–250)

(13–251)

where:

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

(13–252)

where:

Ac = area domain of the contact surface

= time increment

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 ([373]), and Armero and Petőcz ([374])).

Impact Constraints and Contact Forces

In Mechanical APDL, 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–253)

where:

K_n = contact normal stiffness

u_n = contact gap size

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

P^n-1 = normal contact pressure at the end of previous substep

and:

(13–254)

where:

K_s = 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–245 and Equation 13–246 but with additional variables as shown above in Equation 13–253 and Equation 13–254.

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 Mechanical APDL 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 Mechanical APDL, 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.

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 CONTA172, CONTA174, CONTA175, and CONTA177.

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

(13–255)

(13–256)

and

(13–257)

where:

The debonding parameter is defined as:

(13–258)

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–259)

and

(13–260)

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–261)

(13–262)

and

(13–263)

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–264)

and

(13–265)

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–266)

with

(13–267)

(13–268)

(13–269)

(13–270)

where:

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

The debonding modes are based on input data:

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:

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.

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 [415]). Mechanical APDL 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:

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

Archard Wear Model

The Archard wear model (see Archard [416]) 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–271)

where:

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–272)

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.

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([280])) are:

Thermal Contact Conduction

(13–273)

where:

Heat Convection

(13–274)

where:

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

Heat Radiation

(13–275)

where:

Heat Generation Due to Frictional Sliding

(13–276)

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

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

(13–277)

where:

Ac = area domain of the contact surface

= time increment

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–278)

where:

Electrostatic (KEYOPT(1) = 5 or 6)

(13–279)

where:

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

(13–280)

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–281)

where:

The magnetic contact permeance coefficient is defined as:

(13–282)

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

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–283)

(13–284)

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–285)

(13–286)

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–287)

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–288)

(13–289)

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.

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

Diffusion Flow for closed contact and near-field contact :

(13–290)

(13–291)

where:

Diffusion Flow for far-field contact :

(13–292)

where: