Massless Node-to-Face Contact

The massless node-to-face contact formulation has been introduced in CalculiX in order to model perfectly hard contact within explicit dynamic calculations. The implementation closely follows [69].

For the sake of simplicity, frictionless contact is taken as example. Let $g$ denote the gap, $\lambda$ is the contact force.

Figure 135: Normal hard contact

For a zero contact force the gap can be any positive real number including zero ( $\mathbb{R}^+_0$), for a strictly positive force the gap is zero and a negative contact force is not plausible (i.e. no adhesion is assumed). These dependencies can be represented in the form of $g(\lambda)$ as in Figure 135. For dynamic calculations with friction the gap velocity $\gamma=\dot{g}$ is used in the contact laws. In principle, the gap velocity can be positive or negative. Indeed, if an object is approaching an obstacle the gap decreases and the velocity is negative. After collision the gap is increasing and the velocity is positive. However, at collision (i.e. $g=0$) $\gamma \ge 0$ is satisfied (the gap remains zero for static equilibrium or is increasing after a collision) and now Figure 135 also applies for $\gamma(\lambda)$. In order for the velocity to be strictly positive, the force has to be zero.

Now, a couple of mathematical concepts are introduced. For further details the reader is referred to [52] and [102].

• A function is set-valued if it can have more than one value for a given argument. For instance, the gap velocity function satisfies

  
  

  $$\lambda = 0 \Rightarrow 0 \le \gamma \le \infty$$ (322)

  $$\lambda > 0 \Rightarrow \gamma = 0$$ (323)

  
  

  The admissible set on which $\gamma$ is defined is $\mathbb{R}^+_0$.

• A set $\mathcal{C}$ is convex if any linear combination of elements $\boldsymbol{x_1}$ and $\boldsymbol{x_2}$ of the form $(1-\alpha) \boldsymbol{x_1}+ \alpha \boldsymbol{x_2}, 0 \le \alpha \le 1$ also belongs to $\mathcal{C}$. For instance, the interior of a circle in two-dimensional space is convex. This also applies to the interior supplemented by the boundary of the circle. The interior of a ring is not convex. $\mathbb{R}^+_0$ is convex.

• The indicator function $\psi$ for a set $\mathcal{C}$ is defined by:

  $$\psi_{\mathcal{C}}(\boldsymbol{x}) = \left \{ \begin{array}{r l} ... ...\infty & \text{if} \quad \boldsymbol{x} \not\in \mathcal{C} \end{array} \right.$$ (324)

  
  

  So the indicator function is zero for all elements included in the set and infinity else. The indicator function for $\mathbb{R}^+_0$ is shown in Figure 136 (bold line)

  Figure 136: Indicator function for $\mathbb{R}^+_0$

  

• A function $\mathcal{C} \rightarrow \mathbb{R}$ is convex when $\{(\boldsymbol{x},y) \vert y \ge f(\boldsymbol{x}); \boldsymbol{x} \in \mathcal{C}\}$ is a convex set. For instance, for a function on a set in the two-dimensional domain $\mathbb{R}^2$ this means that the three-dimensional volume above the function values has to be a convex set. Since the indicator function $\psi_\mathcal{C}$ is zero everywhere in the set $\mathcal{C}$ it is clear that the indicator function of a convex set is convex.

• Now, the concept of subdifferential can be explained, which is a generalization of the differential for functions with kinks. The subdifferential $\partial f(\boldsymbol{x})$ of a convex function $f(\boldsymbol{x})$ is defined by:

  $$\partial f(\boldsymbol{x}) = \{ \boldsymbol{y} \vert f(\boldsymbo... ...(\boldsymbol{x^*}-\boldsymbol{x}); \forall \boldsymbol{x^*} \in \mathcal{C} \},$$ (325)

  
  

  Figure 137: Subdifferential at several positions

  

  for $f(\boldsymbol{x}) < +\infty$, where “$\cdot$” is the inner product. This means that all function values have to exceed a “tangent” straight line with the subdifferential as slope. As shown in Figure 137 the subdifferential at point $b$, where the function is continuous differentiable, coincides with the derivative. In point $a$, where the function is continuous but not differentiable, the subdifferential consists of all tangent lines in between the left and right derivative at that point. Thus, the subdifferential in $a$ is multivalued and a set-valued function. The same applies to the origin in Figure 136 (dashed lines). Indeed, by comparing Figures 135 and 136 one observes that the subdifferential of the indicator function of $\mathbb{R}^+_0$ coincides with $-\gamma$:

  $$\partial \psi _{\mathbb{R}^+_0} = -\gamma.$$ (326)

  
  

  From analysis it is well known that for a continuous differentiable function a minimum requires that the derivative is zero. From Figure 137 it is obvious that for a minimum of a $C_0$-continuous function

  $$0 \in \partial f(\boldsymbol{x})$$ (327)

  
  

  has to be satisfied. For the indicator function $\psi$ the definition of subdifferential reduces to:

  $$\partial \psi_{\mathcal{C}} (\boldsymbol{x}) = \{ \boldsymbol{y} ... ...(\boldsymbol{x^*}-\boldsymbol{x}); \forall \boldsymbol{x^*} \in \mathcal{C} \},$$ (328)

  
  

  since the indicator function is zero in $\mathcal{C}$.

• Now, the normal cone of $\mathcal{C}$ in $\boldsymbol{x}$ is defined by:

  $$\mathcal{N}_{\mathcal{C}} (\boldsymbol{x}) = \{ \boldsymbol{y} \v... ...0; \forall \boldsymbol{x^*} \in \mathcal{C}; \boldsymbol{x} \in \mathcal{C} \},$$ (329)

  
  

  i.e. it is the set of all vectors which make an angle $\ge 90 ^{\circ}$ with all vectors connecting $\boldsymbol{x} \in \mathcal{C}$ with any other point $\boldsymbol{x^*} \in \mathcal{C}$.

  Figure 138: Normal cone at several positions

  

  Looking at Figure 138 the normal cone in $a$ is $\{\boldsymbol{0}\}$ , in $b$ it is (the vectors on) a straight line locally orthogonal to $\mathcal{C}$ and in $c$ it is (the vectors within) a cone bordered by the dashed lines. By comparing Equations (328) and (329) it is clear that:

  $$\partial \psi_{\mathcal{C}} (\boldsymbol{x}) = \mathcal{N}_{\mathcal{C}} (\boldsymbol{x}) .$$ (330)

  
  

  This establishes a relationship between the subdifferential and the normal cone concept.

• Next, the concept of proximal point is introduced. The proximal point of $\boldsymbol{z}$ w.r.t. $\mathcal{C}$ is defined by:

  $$_{\mathcal{C}} (\boldsymbol{z}) = _{\{\boldsymbol{x^*} \in \; \mathcal{C} \} } \Vert \boldsymbol{z} - \boldsymbol{x^*} \Vert,$$ (331)

  
  

  i.e. it is the point within $\mathcal{C}$ which is closest to $\boldsymbol{z}$. For instance, looking at Figure 138 the proximal point of $\boldsymbol{b^*}$ is $\boldsymbol{b}$ and the proximal point of $\boldsymbol{c^*}$ is $\boldsymbol{c}$. Now, the concept of proximal point can be linked to the normal cone and subdifferential definitions. Indeed:

  
  

  $$\boldsymbol{x} = _{\mathcal{C}} (\boldsymbol{z})$$ (332)

  $$\Leftrightarrow \boldsymbol{x} = _{\{\boldsymbol{x^*} \in \; \mathcal{C}\} } \Vert \boldsymbol{x^*} - \boldsymbol{z} \Vert$$ (333)

  $$\Leftrightarrow \boldsymbol{x} = _{\{\boldsymbol{x^*} \in \; \mathcal{C} \} } \left [ \frac{1}{2} \Vert \boldsymbol{x^*} - \boldsymbol{z} \Vert ^2 \right ]$$ (334)

  $$\Leftrightarrow \boldsymbol{x} = \left [ \frac{1}{2} \Vert \boldsymbol{x^*} - \boldsymbol{z} \Vert ^2 + \psi_{\mathcal{C}} (\boldsymbol{x^*}) \right ]$$ (335)

  
  

  Notice that by adding the indicator function the constraint $\{\boldsymbol{x^*} \in \mathcal{C}\}$ was removed and a convex constrained minimization problem is turned into a convex unconstrained minimization problem. A minimum is obtained if zero belongs to the subdifferential, i.e.:

  
  

  $$\boldsymbol{0} \in \partial \left [ \frac{1}{2} \Vert \boldsymbol{x} - \boldsymbol{z} \Vert ^2 + \psi_{\mathcal{C}} (\boldsymbol{x}) \right ]$$ (336)

  $$\Leftrightarrow \boldsymbol{0} \in \boldsymbol{x}- \boldsymbol{z} + \partial \psi_{\mathcal{C}} (\boldsymbol{x})$$ (337)

  $$\Leftrightarrow \boldsymbol{z} - \boldsymbol{x} \in \partial \psi_{\mathcal{C}} (\boldsymbol{x})$$ (338)

  $$\Leftrightarrow \boldsymbol{z} - \boldsymbol{x} \in \mathcal{N}_{\mathcal{C}} (\boldsymbol{x})$$ (339)

  
  

  This means that one can write:

  $$\boldsymbol{\xi } \in \mathcal{N}_{\mathcal{C}} (\boldsymbol{x}) \Leftrightarrow \boldsymbol{x} = _{\mathcal{C}} (\boldsymbol{\xi }+ \boldsymbol{x} ).$$ (340)

  
  

  This finishes the mathematical excursion.

Applying this to the relationship between the gap velocity $\gamma$ and the normal contact force $\lambda$ one can write:

$$-\gamma \in \partial \psi _ {\mathbb{R}^+_0} (\lambda) \Leftright... ...w - \gamma \in \mathcal{N}_{\mathbb{R}^+_0} (\lambda) \Leftrightarrow \lambda = _{\mathbb{R}^+_0} (-\gamma + \lambda),$$ (341)

which turns a set relationship into a nonlinear equation, which is easier to solve. Notice that solving for a proximal point usually boils down to a projection on the admissable set.

If friction is present the feasible domain consists of the positive real numbers (including zero) for the normal contact, and a disk with a radius equal to the friction coefficient times the normal force:

$$-\gamma_n \in \mathcal{N}_{\mathbb{R}^+_0} (\lambda_n)$$ (342)

$$-\boldsymbol{\gamma_t} \in \mathcal{N}_{S_p} (\boldsymbol{\lambda_t })$$ (343)

$$S_p = \{ \boldsymbol{\lambda_t} \; \vert \; \Vert \boldsymbol{\lambda_t} \Vert \le \mu \lambda_n\}$$ (344)

Notice that the feasible domain is split into a normal and a tangential domain. Therefore, also the projection is split: in normal direction the projection is on $\mathbb{R}^+_0$, in tangential direction on $S_p$.

Figure 139: Contact area

In the following a node-to-face contact definition is assumed (cf. Section 6.7.6) . For the gap definition one can start from MPC's connecting a slave node with the opposite master face, e.g. for node $a$ in Figure 139:

$$\boldsymbol{u_a}= 0.5 \boldsymbol{u_b} + 0.5 \boldsymbol{u_c},$$ (345)

where $\boldsymbol{u}$ represents the displacement field. In the more general case the coefficients have to be determined from the shape functions of the master face. The gap vector can now be written as (assuming the gap is at the beginning of the calculation zero):

$$\boldsymbol{g} = \boldsymbol{u_a}- 0.5 \boldsymbol{u_b} - 0.5 \boldsymbol{u_c}.$$ (346)

This vector can be transformed in a $\{ n, \boldsymbol{t} \}$ coordinate system leading to a normal component and two tangential components for the gap. This can be done for all slave nodes leading to the relationship:

$$\{g\} = [W_b]^T \{U_b\} + \{g_0(t)\}.$$ (347)

$\{g\}$ contains the gap in local coordinates. Its length is 3 times the number of slave nodes in the model (for three-dimensional applications). $\{U_b\}$ is a vector containing the displacements of all contact boundary nodes, i.e. slave nodes and master nodes, in the global orthogonal coordinate system. $\{ g_0(t) \}$ represents the initial gap. Assuming small sliding and small deformations, the matrix $[W_b]$ can be assumed constant. Taking the derivative w.r.t. time of the above equation yields:

$$\{\gamma \} = [W_b]^T \{V_b\} + \{ \dot{g}_0(t) \},$$ (348)

where $\{V_b\}$ denotes the velocity of all contact boundary nodes. Now, the equations of motion are written thereby neglecting the mass and damping terms for the contact boundary nodes (the index i denotes the internal nodes):

$$\begin{bmatrix}0 & 0 \\ 0 & M_{ii} \end{bmatrix} \begin{Bmatrix}\... ...\end{Bmatrix} = \begin{Bmatrix}{F}_b(t) +W_b \lambda \\ {F}_i(t) \end{Bmatrix},$$ (349)

or

$$[K_{bb}] \{U_b\} + [K_{bi}] \{U_i\} = \{F_b(t)\}+[W_b] \{ \lambda \}$$ (350)

and

$$[M_{ii}] \{\dot{V}_i\} + [D_{ii}] \{V_i\}+ [K_{ii}] \{U_i\} + [K_{ib}] \{U_b\} = \{F_i(t) \}.$$ (351)

Here, $\{\lambda \}$ represent the contact forces in the slave nodes in a local coordinate system. The size of this vector is three times the number of slave nodes. $[W_b] \{ \lambda \}$ are the contact forces in all nodes (slave and master) in global coordinates. The motivation for neglecting the inertia terms at the contact boundary comes from the fact that these forces were observed to lead to substantial scatter in the solution in the contact area.

From Equation (350) one obtains for the displacements at the boundary nodes for time increment $k$:

$$\{U_b\}^k = [K_{bb}]^{-1} \left ( \{F_b(t^k)\}-[K_{bi}] \{U_i\}^k+[W_b] \{ \lambda \}^k \right )$$ (352)

Now, a Verlet scheme in its leapfrog form [69] is proposed based on the trapezoidal rule and a shifted grid. This means that displacements are evaluated at times $t^k$ and $t^{k+1}$, while velocities are evaluated at times $t^{k-\frac{1}{2}}$ and $t^{k+\frac{1}{2}}$. Evaluating Equation (348) at time step $k$ and using $\boldsymbol{v} ^{k-\frac{1}{2}}$ one obtains:

$$\{\gamma \}^k = [W_b]^T \frac{(\{U_b\}^k-\{U_b\}^{k-1})}{\Delta t} + \{ \dot{g}_0(t^k) \}.$$ (353)

and after substituting $\{U_b\}^k$:

$$\{\gamma \}^k = \frac{1}{\Delta t}[W_b]^T [K_{bb}]^{-1}[W_b] \{ \lambda \}^k$$

$$+ \frac{1}{\Delta t}[W_b]^T \left ( [K_{bb}]^{-1} ( \{F_b(t^k)\}-[K_{bi}] \{U_i\}^k) -\{U_b\}^{k-1} \right )$$

$$+ \{ \dot{g}_0(t^k) \},$$ (354)

or, introducing the symmetric matrix $[G]$ and vector $\{c\}$:

$$\{\gamma \}^k= [G] \{ \lambda \}^k+\{c\}.$$ (355)

So the contact constraint at time $t^k$ (also called an inclusion problem):

$$-\boldsymbol{\gamma }^k \in \mathcal{N}_{\mathcal{C}} (\boldsymbol{\lambda }^k)$$ (356)

now amounts to:

$$- ([G] \{ \lambda \}^k+\{c\}) \in \mathcal{N}_{\mathcal{C}} (\{ \lambda \}^k),$$ (357)

or

$$\{ \lambda \}^k = _{\mathcal{C}} \left [ \{ \lambda \}^k - ([G] \{ \lambda \}^k + \{c\} ) \right ].$$ (358)

The latter equation can be solved in an iterative way [102]:

$$\{ \lambda \}^k_{n+1} = _{\mathcal{C}} \left [ \{ \lambda \}^k_n - [r]([G] \{ \lambda \}^k_n + \{c\} ) \right ],$$ (359)

where $[r]$ is a diagonal matrix with relaxation factors and $n$ is the iteration index. The relaxation factors are taken to be

$$r_{ii} = \frac{ \omega}{ G_{ii}} \quad G_{ii} > \sum_{j,j\ne i} \vert G_{ij} \vert$$ (360)

$$r_{ii} = \frac{\omega}{ \sum_{j,j\ne i} \vert G_{ij} \vert } \quad G_{ii} \le \sum_{j,j\ne i} \vert G_{ij} \vert ,$$ (361)

where $0 \le \omega \le 2$ is a relaxation parameter.

Summarizing, the contact displacements in time step $t^k$ are solved using the following steps (knowing $\{F_b(t^k)\}$ , $\dot{g}_0(t^k)$, $\{U_i\}^k$ and $\{U_b\}^{k-1}$):

• determine $[G]$ and $\{c\}$ using Equations (354) and (355)
• solve for $\{ \lambda \}^k$ using Equation (359) with $\{ \lambda \}^{k-1}$ as initial guess.
• solve for $\{U_b\}^k$ using Equation (352).

The solution of the contact problem, however, is restricted to the active contact degrees of freedom. Indeed, only for the nodes in contact the gap velocity is positive and Figure 135 applies. The procedure to determine these active contact degrees of freedom is called an active set strategy. Two cases are considered. For the initially open case (i.e. no contact at $t^{k-1}$) the gap is calculated at $t^k$ by substituting Equation (352) for $\{ \lambda \}^k=0$ into Equation (347). The active degrees of freedom are those for which the gap is nonpositive. For the preloaded case (i.e. the contact forces at $t^{k-1}$ are positive) the preload is calculated from Equation (352) assuming sticking contact, i.e. $\{U_b\}^k=\{U_b\}^{k-1}$. Then, the active degrees of freedom are those, for which this sticking contact is changing, i.e. either the normal preload is negative or the tangential preload exceeds the sticking range and slip occurs. For details the reader is referred to [69].

Once $\{U_b\}^k$ is calculated $\{V_i\}^{k+\frac{1}{2}}$ can be calculated by substituting $\{\dot{V}_i\}^k=(\{V_i\}^{k+\frac{1}{2}}-\{V_i\}^{k-\frac{1}{2}} )/\Delta t$ and $\{V_i\}^k=(\{V_i\}^{k+\frac{1}{2}}+\{V_i\}^{k-\frac{1}{2}} )/2$ in Equation (351) expressed at time $t^k$ leading to:

$$\left( \frac{[M_{ii}]}{\Delta t} + \frac{[D_{ii}]}{2} \right) \{V_i\}^{k+\frac{1}{2}} = \{F_i(t^k) \}-[K_{ii}] \{U_i\}^k-[K_{ib}]\{U_b\}^k$$

$$+ \left( \frac{[M_{ii}]}{\Delta t} - \frac{[D_{ii}]}{2} \right)\{V_i\}^{k-\frac{1}{2}}.$$ (362)

The solution of this set of equations requires a linear equation solver. The mass matrix does not change during the calculation. If the damping matrix does not change either, the factorization step in the linear equation solver can be done just once at the start of the calculation. This drastically reduces the computation time. Knowing $\{V_i\}^{k+\frac{1}{2}}$ the value of the displacements in the internal nodes can be obtained from:

$$\{U_i\}^{k+1} = \{U_i\}^k + \{V_i\}^{k+\frac{1}{2}} \Delta t.$$ (363)

Consequently, the overall algorithm can be summarized as follows, knowing $\{F_b(t^k)\}$ , $\dot{g}_0(t^k)$, $\{U_i\}^k$, $\{U_b\}^{k-1}$, $\{ \lambda \}^{k-1}$ and $\{V_i\}^{k-\frac{1}{2}}$:

• Determine the active set
  • If the active set is empty $\{\lambda\}^k=\{0\}$ for an initially open contact, and $\{\lambda\}^k=\{\lambda ^{\text{pre}}\}$ for a preloaded contact.
  • If the active set is not empty, determine $\{ \lambda \}^k$ from Equation (357) with $\{ \lambda \}^{k-1}$ as starting value.
• Determine $\{U_b\}^k$ from Equation (352).
• Determine $\{V_i\}^{k+\frac{1}{2}}$ from Equation (362).
• Determine $\{U_i\}^{k+1}$ from Equation (363).

In the first increment (k=1) $\{U_i\}^1$, $\{U_b\}^{0}$ and $\{V_i\}^{\frac{1}{2}}$ have to be known.

Right now, the following limitations apply to the massless contact method:

• If NLGEOM is active on the *STEP card selective mass scaling is not reliable and should not be applied, i.e. the user should not specify a minimum increment size underneath the *DYNAMIC card.
• If NLGEOM is active on the *STEP card large rotations will lead to wrong results.