Massless contact

Massless contact is coded uniquely for explicit dynamics. The underlying contact definition uses a node-to-surface approach and fields such as islavnode and imastnode are common with the penalty contact formulation.

Within a new step some of the fields for massless contact have to be set up only once in the step. Therefore, they are calculated at the beginning of nonlingeo.c in the section in which the initial acceleration for the $\alpha$-method is calculated (just after setting up the matrices in mafillsmmain.c).

At first the degrees of freedom needed to set up the $[W_b]$ matrix (cf. Section 6.7.9) are collected in create_contactdofs.f. These are:

• the slave node degrees of freedom. Each of them is characterized by a node number “node” and a global direction “j”. The value of nactdof(j,node) is stored in field kslav, the value 10*node+j in lslav. Subsequently kslav is sorted in ascending order taking lslav along. The size of these fields is 3*nslavs, since no SPC's or MPC's are allowed in the contact area for massless contact.
• the slave and master nodes degrees of freedom. They are stored in a similar way in fields ktot and ltot. The size of these fields is neqtot:=3*nslavs+3*nmasts.

Furthermore, in the same routine the friction coefficient in each slave node is stored in field fric(1...nslavs).

Then, the matrices in between the big brackets in Equation (362) are calculated. They consist of a linear combination of the mass and damping matrix, or, since Rayleigh damping is assumed, a linear combination of the mass and stiffness matrix. Subsequently, these matrices are reduced to the internal degrees of freedom only (recall that $[M_{ii}]$ and $[D_{ii}]$ is used in Equation (362)) in subroutine reducematrix.f. Finally, the matrix corresponding to the left hand side is factorized. It will be repeatedly used for solving the system in each increment.

Now, the increment loop can start. The following actions are performed:

• At the beginning of each increment the $[W_b]$ matrix is formed. This is performed in subroutine gencontelem_n2f.f. Each set of three columns in the matrix corresponds to the dependence of the degrees of freedom of a slave node in a local coordinate system (splitting the normal direction from the tangential directions) on the degrees of freedom of the nodes of the opposite master face in global coordinates. The matrix is stored using fields auw, jqw, iroww and nzsw. Vector auw contains the nonzero values column by column, the corresponding row number is stored in iroww, in total there are nzsw entries. The first entry in each column is stored in jqw(i), i=1,..3*nslavs+1. [for linear calculations (i.e. without geometric and material nonlinearities; characterized by masslesslinear $>$ 0) $[W_b]$ is calculated only once at the start of the step.]
• contrary to penalty dynamic calculations, no extrapolation from the previous result is done, i.e. prediction.c is skipped.
• contrary to penalty dynamic calculations, the stiffness matrix is assembled in mafillsmmain.c. [For linear calculations the stiffness matrix is calculated only once at the start of the step.]
• contrary to penalty dynamic calculations, no residual is calculated in calcresidual.c. Instead, routine massless.c is called with the following actions:
  • The rows in $[W_b]$ are rearranged; in gencontelem_n2f.f the order is such that first the slave nodes are treated, followed by the master nodes. In routine expand_auw the new row order is the one dictated by field nactdof. [For linear calculations this is done only once at the start of the step.]
  • In routine extract_matrices $[K_{bb}]$, $[K_{bi}]$,$[K_{ib}]$ and $[K_{ii}]$ are extracted from the global stiffness matrix. Matrix $[K_{bb}]$ is stored using fields jqbb, aubb, adbb, irowbb and icolbb (the latter contains the number of nonzero entries in each column), because of symmetry only half the matrix is stored. The subdiagional terms of matrix $[K_{bi}]$ are stored using fields jqbi, aubi, irowbi and nzsbi. The subdiagonal terms of matrix $[K_{ib}]$ are stored using fields jqib, auib, irowib and nzsib. To obtain matrix $[K_{bi}]$ it has to be complemented by the transpose of the $[K_{ib}]$ matrix. Notice that matrix $[K_{ib}]$ is assumed to have as many rows as $[K]$ and three times nslavs columns (any entries not belonging to the contact nodes are set to zero; this does not blow up the storage, since only non-zeros are stored). Mutatis mutandis this also applies to $[K_{bi}]$. Matrix $[K_{ii}]$ is actually not extracted, the entries of matrices $[K_{bi}]$, $[K_{ib}]$ and $[K_{bb}]$ are set to zero (off-diagonal) or one (diagonal). [For linear calculations this is done only once at the start of the step.]
  • Calculation of $\{F_b(t^k)\} - [K_{bi}]\{U_i\}^k$ in resforccont.f. This is the residual force due to external forces and forces from the internal nodes in the absence of contact forces $\{ \lambda \}^k$. For linear calculations the external force due to distributed loads (e.g. centrifugal forces) is calculated only once at the start of the step (only if the distributed load does not change in the step). Concentrated loads are added on an increment basis.
  • factorizing $[K_{bb}]$ (for linear calculations this is done only once at the start of the step) and premultiplying the residual force with $[K_{bb}]^{-1}$ by solving a linear equation system. The result is stored in field gapdisp and is the displacement field at the boundary nodes in global coordinates assuming zero contact force.
  • calculating the gap by premultiplying gapdisp by $[W_b]^T$ and adding $\{g_0\}$ in detectactivecont.f; determining the active contact degrees of freedom by identifying negative gap values. The number of active degrees of freedom is nacti, the corresponding entries in iacti(j),j=1...3*nlavs are numbered in ascending order, the nonactive degrees in iacti are set to zero.
  • calculate $[G]$ and $\{c\}$ (using gapdisp) and storing them in gmatrix(1..nacti,1..nacti) and cvec(1..nacti). For linear calculations a maximal $[G]$ matrix (stored in fullgmatrix) is calculated once at the start of the step assuming all contact degrees of freedom do be active. it is appropriately condensed to the right size in each increment based on the truly active degrees of freedom. In the same way a maximal relaxation vector fullr is calculated.
  • solving the inclusion problem in auglag_inclusion.f. Field alglob(1..neqtot) contains $[W_b] \{ \lambda \}$, i.e. the contact forces in slave and master nodes in global coordinates, field al(1..3*nslavs) contains $\{\lambda \}$, i.e. the contact forces in the slave nodes in a local coordinate system.
  • calculate $\{U_b\}^k$ using Equation (352).
  • calculate the right hand side of Equation (362).

• Solve the equation system corresponding to Equation (362). Since the left hand side is stored in fields adb and aub this is covered by the regular penalty dynamic calculations and does not require special treatment.
• Calculate a pseudo velocity in the contact boundary nodes by calculating $(\{U_b\}^k-\{U_b\}^{k-1})/\Delta t$; this is needed for the next step.
• Use routines resultsini.c and iniparll.c to calculate $\{U_i\}^{k+1}$ and $\{U_b\}^{k}$.