`***contact`#

Description#

The ***contact section defines zones of contact in static and dynamic mechanical problems. Typically, applying contact is one of the more troublesome aspects of nonlinear finite element analysis, because of increased discontinuities and sometimes unavoidable over- and underconstraints. This is above the obvious increase in modeling and input file complexity. After Z-set version 8.2, major updates to the contact system have been made which hopefully increase the robustness and quality of contact mechanics. The algorithms of Z-set are particularly good at avoiding global convergence issues with chatter, and in performing well when there are large differences in the relative stiffnesses of the contacting parts.

Contact is enforced with a group of contact zones. Zones are modeled by a node set acting as the impactor (slave) surface which may come in contact with a target (master) surface area. Each zone may have a different contact model (law) with or without frictional behaviors, and can include a variety of control parameters. The contact condition is defined in terms of a contact stress, but can alternatively be defined directly in terms of nodal reaction forces. ¹ ¹Note that for second order interpolation functions the magnitude and even direction of the nodal reactions is not intuitive. The fact is that treating nodal reactions directly assumes a discrete force acting over a zero area, which of course results in a singular stress in a continuum. Thus for these higher interpolations the only real choice is to deal with stresses distributed over an area.

To define contact zones, one must identify target/master surface areas, and corresponding impact/slave sets of nodes. Under normal use, i.e. stress-based contact, a node set (nset) is required for the impactor side, and a boundary set (liset or faset) is required for both target and impactor. ² ²This is in contrast to older versions where a bset was required for the target only. A facility is available to automatically generate the impactor bset from its nset, but in general it is preferable to set this up ahead of time. One can use the ***mesher **bset command (page ) with the *use_nset sub-command. The boundary sets must be set up so that the normals point outward (see page ), which can be verified graphically in Zmaster using the Bset normals option of the Mesh Bset... command. An example view of a properly set up two-dimensional example is given in the above image, where the target nset is not visible because it coincides with its corresponding liset/faset (as it should).

Note that if contact is occurring on only one node, or if there is a zero radius corner in contact, the contact stress will still be distributed over the adjacent boundaries. If the purpose of this modeling is to impose a boundary condition, or if the frictional behavior is not important, it is possible to use the **force_basis option, or to define the impactor boundary set with a dot set. In the latter case, each node has unit area and does not have any additional connectivity.

Solution method#

Contact introduces an additional system of constraints on the displacement field. These are implicitly defined in terms of the contact stresses. The contact strategy in Z-set consists in solving a sub-problem for the contact stresses and relative displacements between surfaces, from which the additional global constraints are computed (a process very similar to a sub-structuring problem). The nonlinear contact solution is thus done by running “sub-space iterations” depending on the chosen contact law, and a flexibility matrix found from the global stiffness. The flexibility matrix is a (possibly) full matrix describing the deformation of each node on a contacting surface due to a unit force on each other node on the contact surface. Note that these terms will only be zero if the nodes are on separate bodies. The contact-modified solution intervenes in the global iterative procedure in several steps:

An estimate of the global field is made without enforcing the contact in the current iteration, by applying the classic Newton method using the last converged residual. This residual consists of both the internal reactions and the added contact reactions from the last iteration.
The flexibility matrix is computed from the global stiffness matrix by using repeated back solves of unit loads over all contact nodes. Because back solves are being made, it is very advantageous to have all the contact nodes numbered at the end of the global unknowns. However, this will in general increase the total bandwidth of the problem. Very often the computation of the flexibility matrix will be the most important part of the solution time of a contact problem. Currently the sparse matrix solvers include a compromise renumbering scheme for best performance.
The next step is then proceeding with the “inner” iterations by using this flexibility matrix to solve for the contact stresses and displacements. There are both direct solvers (involving the inverse of the flexibility matrix) and block diagonal iterative solvers in case that direct solution becomes too costly. The solver can be chosen with the **solve_method option. By default, the currently implemented contact laws use a direct solver. Running in verbose mode (either using the -v command line switch, the **verbose output option, or by using the **verbose_contact contact option) will show a great deal of additional output on the progress of the local contact solution. This option should be enabled if there are any difficulties with the contact solution.
The last stage updates the set of global residual equations to include the contact reactions, which will in turn be used for the next dof prediction in the global Newton iterations.

The contact algorithm outlined above has many advantages, at the cost of computing the flexibility matrix. Firstly the method isolates the severe nonlinearities involved in solving the contact state from the global material and geometry nonlinearities. In particular, the effect of “chattering” is almost totally isolated to the local contact iterations, where it is very much more efficiently treated. In addition, the local system defining the contact is nonsymmetric in the case of sliding friction, but because these terms are solved only in the local iterations, the global matrix can remain symmetric even under high friction (as long as other equations have led to a symmetric tangent before the contact has been considered).

Contact behavior#

The contact behavior is defined in terms of the quantities schematically represented in the following figure:

As shown, there is a normal mode of contact (to be enforced to null within a given limit), and a shear mode with a displacement vector defined in local coordinate axes \(\left( \xi_1,\xi_2 \right)\) of the face element. The following variable definitions are used in the descriptions of the contact laws:

\(\sigma\) the normal contact stress (scalar), defined with compressive values being positive.
\(\vect \tau\) the tangent stress vector.
\(u_n\) the normal displacement (scalar) between potentially contacting surfaces. Values greater than zero have a gap between the surfaces and values less than zero indicate a penetration.
\(d\vect u_t\) the incremental tangent displacement (vector), measuring the projected slip on the element face during the increment.
\(k_n, k_t\) intrinsic contact stiffnesses in the normal and tangent directions (units of stress/length), as estimated by the code from the global matrix.
\(\mu\) (not drawn, dimensionless scalar) the friction coefficient as a fraction of the normal stress.

With these variable definitions we describe the following contact models:

normal contact with the following auxiliary variables \(\hat{\sigma} = \sigma - k_n u_n\) and \(\hat{\vect \tau} = \vect \tau - k_t d \vect u_t\), the normal contact law is expressed as

(55)#\[\begin{split}\begin{aligned} &\hbox{if~~} \hat{\sigma}>0 \\ &\hskip.5cm \hbox{then} \\ &\hskip1cm k_n u_n = 0 \\ &\hskip1cm \hbox{if~~}\vert\hat{\vect \tau}\vert < \mu \sigma \\ &\hskip1.5cm \hbox{then~~} k_t d\vect u_t = 0 \\ &\hskip1.5cm \hbox{else~~} \vect \tau=\mu\sigma\vect n_t~~~\hbox{with}~~~ \vect n_t=\hat{\vect \tau}/\vert\hat{\vect \tau}\vert \\ &\hskip.5cm \hbox{else~~} \sigma=0, \vect \tau=0 \end{aligned}\end{split}\]

The normal and tangential stiffnesses \(k_n>0\) and \(k_t>0\) are normalizing parameters which are automatically set by the program.

penalty contact is identical to normal contact, except that it allows for the contact constraints to be violated using the penalty parameters in the normal and tangential directions \(\gamma_n\) and \(\gamma_t\). /The larger the penalty, the more the constraint is satisfied./ Numerically, these penalty parameters tend to make the local contact matrix diagonally dominant and can condition the matrix under situations where other constraints conflict with the contact constraints. When the penalty parameters approach infinity, the penalty contact model behaves exactly as the standard normal contact. In addition, in problems with dynamic impact, it is useful to include a term \(\beta\) which effectively damps any oscillations introduced by severe impact, while introducing only a small penalty under less impulsive conditions. Defining the auxiliary forces \(\hat{\sigma} = \sigma - \left[\gamma_n^{-1}\sigma + k_n(u_n + \beta\gamma_n^{-1}\,du_n)\right]\) and \(\hat{\vect \tau} = \vect \tau - \left[\gamma_t^{-1}\,d\vect \tau + k_t\,d\vect u_t\right]\), the penalty contact model is expressed as

(56)#\[\begin{split}\begin{aligned} &\hbox{if~~} \hat{\sigma}>0 \\ &\hskip.5cm \hbox{then} \\ &\hskip1cm \gamma_n^{-1}\sigma + k_n(u_n + \beta\gamma_n^{-1}\,du_n) = 0 \\ &\hskip1cm \hbox{if~~}\vert\hat{\vect \tau}\vert < \mu\sigma \\ &\hskip1.5cm \hbox{then~~} \gamma_t^{-1}\,d\vect \tau + k_t\,d\vect u_t = 0 \\ &\hskip1.5cm \hbox{else~~} \vect \tau=\mu\sigma\vect n_t~~~ \hbox{with}~~~\vect n_t=\hat{\vect \tau}/\vert\hat{\vect \tau}\vert \\ &\hskip.5cm \hbox{else~~} \sigma=0, \vect \tau=0 \\ \end{aligned}\end{split}\]

is almost identical to penalty contact, with the exception that the friction stress is limited by a user-defined maximum shear strength \(\tau_{max}\). Thus, the normal contact model is a subset of the penalty contact model, and the penalty contact model is a subset of the isotropic Coulomb contact model. Defining the auxiliary forces \(\hat{\sigma} = \sigma - \left[\gamma_n^{-1}\sigma + k_n(u_n + \beta\gamma_n^{-1}\,du_n)\right]\) and \(\hat{\vect \tau} = \vect \tau - \left[\gamma_t^{-1}\,d\vect \tau + k_t\,d\vect u_t\right]\), the isotropic Coulomb contact model is expressed as:

(57)#\[\begin{split}\begin{aligned} &\hbox{if~~} \hat{\sigma}>0 \\ &\hskip.5cm \mbox{then} \\ &\hskip1cm \gamma_n^{-1}\sigma + k_n(u_n + \beta\gamma_n^{-1}\,du_n) = 0 \\ &\hskip1cm \hbox{if~~} \vert\hat{\vect \tau}\vert < \hbox{min}(\mu\sigma, \tau_{max}) \\ &\hskip1.5cm \hbox{then~~} \gamma_t^{-1}\,d\vect \tau + k_t\,d\vect u_t = 0 \\ &\hskip1.5cm \hbox{else~~} \vect \tau= \hbox{min}\left(\mu\sigma,\tau_{max}\right)\vect n_t~~~ \hbox{with}~~~\vect n_t=\hat{\vect \tau}/\vert\hat{\vect \tau}\vert \\ &\hskip.5cm \hbox{else~~} \sigma=0, \vect \tau=0 \\ \end{aligned}\end{split}\]

is a generalization of the isotropic coulomb model.

Syntax#

Contact takes a number of main control commands, and sub-blocks defining the different zones and behavior models. One can alternately define the contact model to apply to all zones, or zone-by-zone.

The **soft_param and **penalty_param commands are no longer supported in versions 8.2 and newer. This is now specified using the **zone *behavior_coef command, as described below for each specific contact model. A warning will be issued, but these commands will be ignored.

***contact [ contact model ]: \(~\,\) **zone [ contact model ] \(~\,~\,\) zone subcommands [ **conv precision max_iter [ force_precision ] ] [ **stable eps [ precondition``|``damp ] ] [ **init_d_stress [ sequence ] ] [ **force_basis ] [ **solve_method direct``|``iterative ] [ **limit_activity t1 t2 ] [ **spy_node node [ node ] … ] [ **verbose_contact [ filename ] ]

**zone

defines a contact zone, as well as the law defining its contact model. Current (as of version 8.2) contact models are normal, penalty, coulomb and ortho_coulomb. The keyword soft is deprecated and automatically converted to penalty instead. More than one zone may be specified. If multiple zones exist which all use the same contact model, then the contact model may be specified immediately after the ***contact command. Available sub-options are described on page .

**conv

defines the convergence parameters. The first (real) parameter precision defines a non-dimensional tolerance defining the reduction of the contact residual relative to the residual at the beginning of the inner contact iterations (i.e. the global problem residual without the additional contact influences of this global iteration included). A moderate relative value is therefore quite acceptable, even for problems needing a high absolute final residual.

The default value of 1.e-5 is quite small, so in most cases this command will be useful. A recommended value of 1.e-2 indicates that the contact convergence is two orders of magnitude smaller than the current global convergence state. Unless it is believed that the influence of contact is very very significantly altering the global behavior, this value should be ok.

The second (integer) parameter is the number of iterations before contact divergence will be signaled (default 200). The last optional (real) parameter force_precision indicates the desired magnitude of the contact residual relative to the magnitude of contact stress (default 1.e-5). Only one of the two convergence criteria needs to be satisfied in order to achieve convergence.

Unlike previous versions, the contact solution algorithm in 8.2 always assigns the initial guess of the contact stresses based on the contact stress determined from the previous global iteration. As the global iteration converges, the initial absolute error in the local contact iteration will also decrease. It is recommended that force_precision be set only as tight as the expected convergence rate of the global iteration. When the global iterations converge, the absolute error in the contact law will be satisfied to the same degree as the global iterations.

**stable

takes one real value eps (default 1.e-2) and the choice between precondition (default) or damping. If the static model contains multiple bodies and at least one of these bodies has a rigid body mode (if the contact constraints are not considered), then this command is absolutely required. For dynamic problems it is not needed. The command instructs Z-set to add a term on the diagonal of the stiffness matrix of the contact nodes, in effect adding a grounded incremental spring to each node involved in contact. The value of eps then denotes the relative magnitude of the spring constant.

The specifier precondition instructs Z-set to add only the incremental spring to the stiffness matrix. It does not add the force of the spring to the external forces. This option stabilizes the contact solution algorithm, but in no way affects the converged solution.

With the specifier damping, not only the incremental spring is added to the stiffness matrix, but also the force of the incremental spring is added to the external forces. This actually affects the converged solution. This specifier can be conveniently used to avoid quasi-static singularity problems when a free body is not initially in contact in the initial loading stages. A small value relative to the problem stiffness is recommended, as is a convergence study to ensure that this numerical technique does not significantly alter the structural behavior.

**init_d_stress

causes the algorithm to make an initial guess for the stress increments based on the previous solution, thereby possibly accelerating convergence. This has advantages especially for monotonic loading paths. For non-monotonic loading paths, the option sequence disables this guess for the first increment of each sequence, so that the contact algorithm is not confused too much by sudden changes in loading direction. This command is somewhat similar in spirit to the ***resolution **init_d_dof command (see page ). Without this command, the initial guess for the stress increment is set to zero, that is, the initial guess of the contact stress at the end of the current increment is assigned a value equal to the contact stress at the end of the previous increment.

**force_basis

sets the algorithm based on forces instead of stresses. Synonyms: **lumped and **forced_based.

**solve_method

specifies the type of solver used for the contact solution. The option direct (default value for all currently implemented contact laws) uses the recommended direct solver. A second option iterative uses less memory, but requires more iterations to converge.

**limit_activity

only activates the contact conditions between times \(t_1\) and \(t_2\) (taken as two real parameters).

**spy_node

causes Z-set to print out more detailed information of the contact stresses at the specified list of nodes. These nodes must lie in the node sets given under **zone impactor.

**verbose_contact

gives additional information about the current state of all contact nodes involved. This information is written to an output file filename, or in absence of a user-supplied filename to a file called problem.contact_info.

Advice when problems occur#

The following hints may help out with common issues of setting up contact problems.

Use the **verbose_contact contact option to examine the progress details.
Check the target (master) surface normals to make sure that they point outward from the target body. This can be easily checked in Zmaster’s Mesh mode, under the Bsets.. popup to observe bset normals. In the case that some normals are inverted, one can use the inverse_liset (page ) or faset_align mesher commands (page ).
Check to see that the target is the stiffer and/or less refined surface.
Check the contact behavior. All contact models from versions before 8.3 are deprecated. The preferred “standard” model is coulomb which includes penalty parameters for “soft” behavior in the case of overconstraint, can optionally capped friction behavior.
Check if any of the bodies are unconstrained in any of the degrees of freedom if not for contact. In this case, the global matrix will be singular and will no doubt cause problems during the first prediction of new displacement fields. In this case use the **stable precondition command with a small value (e.g. 1.e-4). This command adds a small value to the global matrix pivots for the contact nodes, stabilizing the global matrix conditioning, but not actually altering the residual computation (i.e. not changing the physics at all). Of course a large stabilization value makes the global matrix less correct, and eventually convergence will be inhibited by this parameter. This method is only needed for static problems.
Check to see if a body is underconstrained in some dofs, and may not always be in contact. If this is the case one can use the **stable damp command which will include incrementally defined springs on the contact nodes with a given stiffness. Please be certain (by doing an appropriate convergence study) that the value of this parameter does not actually alter the results.
Check to see if there are impactor nodes near the point of expected contact. If the slave surface is too coarse or the warning distance is too large, the contact may not be properly detected. In particular, if the impactor is too coarse there may be cases where the impactor does not penetrate the target, but the target does penetrate the impactor. The definition of the master/slave relationship of contact zones allows the master nodes to penetrate the slave.
Check to see if there are sharp convex points in the contact zone. This can cause problems because the surface normals vary tremendously on the master surface, or by redistributing the contact stress erroneously around corners of slave surfaces.
Check the integration order of the elements. Generally underintegrated elements will produce less oscillations (spatial and temporal) in the solution. This is especially the case under highly hydrostatic constraint or in locations where there is a rapid change in the contact condition (e.g. Hertzian contact problems). This will be especially important for isochoric plasticity and viscoplasticity problems, or otherwise close to incompressible elements.
Check for odd convergence and ‘chattering’’ when the contact residual rocks back and forth between two residual values. Sometimes there is a problem with a node alternating between sticking and sliding contact states. In such an instance increase the slipping penalty parameter.

`**zone`#

Description#

Contact zones define different contact pairs (impactor/target), zone by zone control parameters, and contact material behaviors.

Syntax#

The zone syntax is summarized below:

**zone [ contact model ] \(~\,\) *impactor impactor_zone \(~\,\) *target target_zone \(~\,\) *behavior_coef \(~\,~\,\) key-value-pairs [ *dont_check_bcs ] [ *variable_friction file [ position ] ] [ *warning_distance warning_distance ]

**zone: specifies a contact zone with its specific contact model. Details of each model can be found in the following pages. Currently (as of Z-set 8.2), valid choices are normal, penalty, coulomb and ortho_coulomb. Multiple **zone commands may exist, each with its own contact model. If all contact zones have the same contact model, this model may be specified immediately after the ***contact command, instead of here. There is no default behavior. Parameters of the contact model are given with the *behavior_coef command.
*impactor (or alternatively *slave): needs the nset impactor_zone. The associated bset (liset or faset) has to exist as well. The bset has to be oriented with normals pointing away from the interior of the associated body.
*target (or alternatively *master): takes the bset target_zone. Again, this bset must be oriented properly.
*behavior_coef: sets the parameters of the specific contact model. Details can be found in the following pages.
*dont_check_bcs: specifies that the code will not check whether contact conditions are compatible with other boundary conditions that may exist at the same node, as it does by default. With the default behavior, the code will try to fix any incompatibilities, or if that is not possible, the contact condition will be ignored. A warning message will be given if this is the case. This option disables the verification for the current zone.
*variable_friction (or *file): specifies a variable friction coefficient \(\mu\), depending on the current and/or cumulated value of tangential slip. This variable friction cannot (yet) depend on any external parameter, such as temperature. A precise description of this behavior can be found in the material manual under ***behavior variable_friction. Not implemented for the ortho_coulomb contact model.
*warning_distance: specifies the maximum distance warning_distance for which the algorithm will try to detect if contact occurs between target and impactor of the current zone. In other words, if the distance between target and impactor is larger than this warning distance, the algorithm will not try to detect contact. This distance should be chosen larger than the decrease in distance between target and impactor during the increment of interest, otherwise there is a risk of penetration without the contact algorithm coming into action.

The subcommands *friction and *gap are now deprecated, as they are replaced with their equivalents in the *behavior_coef subcommand. Z-set will issue a warning, and in the case of *friction automatically convert to the new format. The *gap command will be entirely ignored.

`**zone normal`#

Description#

This is the normal contact model as described a few pages back. For the meaning of its coefficient the reader is referred to that description.

Syntax#

The normal behavior zones take all the standard **zone commands, with the following specific behavior coefficient:

**zone normal [ other zone subcommands ] \(~\,\) *behavior_coef \(~\,~\,\) friction mu

The default value is mu = 0.

`**zone penalty`#

Description#

This is the penalty contact model as described a few pages back. For the meaning of its coefficients the reader is referred to that description.

Syntax#

The penalty behavior zones take all the standard **zone commands, with the following specific behavior coefficients:

**zone penalty [ other zone subcommands ] \(~\,\) *behavior_coef \(~\,~\,\) friction mu \(~\,\) [ penalty_normal gamma_n ] \(~\,\) [ penalty_slip gamma_t ] \(~\,\) [ damp_normal beta ]

Default values are mu = 0., gamma_n = 1.e12, gamma_t = 1.e2 and beta = 1.

`**zone coulomb`#

Description#

This is the “standard” contact model for contact with Coulomb friction, as described a few pages back. For the meaning of its coefficients the reader is referred to that description. Note that this contact model effectively contains the normal and penalty contact models.

Syntax#

The Coulomb behavior zones take all the standard **zone commands, with the following specific behavior coefficients:

**zone coulomb [ other zone subcommands ] \(~\,\) *behavior_coef \(~\,~\,\) friction mu \(~\,\) [ max_shear tau_max ] \(~\,\) [ penalty_normal gamma_n ] \(~\,\) [ penalty_slip gamma_t ] \(~\,\) [ damp_normal beta ]

Default values are mu = 0., tau_max = 1.e12, gamma_n = 1.e12, gamma_t = 1.e2 and beta = 1.

`**zone ortho_coulomb`#

Description#

This is the orthotropic Coulomb contact model. It is similar to the isotropic Coulomb model, except that some parameters now have two be supplied with two values instead of one, and that the initial material orientation (vector) with respect to the master or slave surface now has to be specified.

Syntax#

The orthotropic Coulomb behavior zones take all the standard **zone commands, with the following behavior coefficient key-value pairs:

**zone ortho_coulomb [ other zone subcommands ] \(~\,\) *behavior_coef \(~\,~\,\) friction mu_1 mu_2 \(~\,\) [ max_shear taumax_1 taumax_2 ] \(~\,\) [ penalty_normal gamma_n ] \(~\,\) [ penalty_slip gamma_s1 gamma_s2 ] \(~\,\) [ damp_normal beta ] \(~\,\) [ direction nx ny nz ] \(~\,\) [ master | slave ]

Default values are mu_1 = mu_2 = 0., taumax_1 = taumax_2 = 1.e12, gamma_n = 1.e12, gamma_s1 = gamma_s2 = 1.e2, beta = 1., (nx ny nz) = (1. 0. 0.), and the last keyword has default value slave. Note that this behavior is currently not compatible with the *variable_friction zone subcommand.

***contact

On this page

***contact#

Description#

Solution method#

Contact behavior#

Syntax#

Advice when problems occur#

**zone#

Description#

Syntax#

**zone normal#

Description#

Syntax#

**zone penalty#

Description#

Syntax#

**zone coulomb#

Description#

Syntax#

**zone ortho_coulomb#

Description#

Syntax#

`***contact`#

`**zone`#

`**zone normal`#

`**zone penalty`#

`**zone coulomb`#

`**zone ortho_coulomb`#