Contact formulations

Contact formulations are used in Solfec-1.0 to facilitate solution of the constraints equation

\[\mathbf{C}\left(\mathbf{U},\mathbf{R}\right)=\mathbf{0}\]

once contact points are detected between bodies and included into the constraints. Depending on a solver type a contact formulation can be either used locally (on a single contact point level) or globally (on the level of a all contact points simultaneously). This however does not affect the mathematics of a contact formulation itself. Sections below describe the formulations employed by Solfec-1.0 and point to references, source code sections, and descriptions of solvers that utilize them.

The frictional contact law

The frictional contact law in Solfec-1.0 employs the so called Signorini–Coulomb conditions. The velocity based Signorini non–penetration condition reads

(74)\[\bar{U}_{N}\ge0\,\,\,R_{N}\ge0\,\,\,\bar{U}_{N}R_{N}=0\]

where \(\bar{U}_{N}=U_{N}^{t+h}+e\min\left(0,U_{N}^{t}\right)\), \(e\) is the velocity restitution coefficient, \(U_{N}\) is the the normal relative velocity, and \(R_{N}\) is the normal reaction. The normal direction is consistent with the positive gap velocity so that (74) states, that any violation of the non–penetration results in a reaction force or velocity driving at the penetration–free configuration. While using \(\bar{U}_{N}\) allows to account for the Newton impact law, for models where multiple impacts occur during one time step, using \(e>0\) cannot be justified from a theoretical standpoint. The Coulomb’s friction law reads

(75)\[\begin{split}\left\{ \begin{array}{ll} \left\Vert \mathbf{R}_{T}\right\Vert \le\mu R_{N}\\ \left\Vert \mathbf{R}_{T}\right\Vert <\mu R_{N} & \Rightarrow\mathbf{U}_{T}=\mathbf{0}\\ \left\Vert \mathbf{R}_{T}\right\Vert =\mu R_{N} & \Rightarrow\exists_{\lambda\ge0}\mathbf{U}_{T}=-\lambda\mathbf{R}_{T} \end{array}\right.\end{split}\]

A friction force smaller than \(\mu R_{N}\) implies sticking, while sliding occurs with the force of value \(\mu R_{N}\) and direction opposite to the slip velocity.

Projected gradient formulation

This is an implicit formulation based on 1 and it is optionally used on an individual contact point level within the Gauss–Seidel solver. On the point level this formulation is implemented in dbs.c:35. The Signorini condition (74) is expressed as a projection

(76)\[R_{N}=\textrm{proj}_{R_{+}}\left(R_{N}-\rho\bar{U}_{N}\right)\]

where vector \(R_{N}-\rho\bar{U}_{N}\) is projected onto the set of positive real numbers \(R_{+}\). Similarly, the Coulomb law (75) is expressed as a projection as follows

(77)\[\mathbf{R}_{T}=\textrm{proj}_{D\left(\mu R_{N}\right)}\left(\mathbf{R}_{T}-\rho\mathbf{U}_{T}\right)\]

where \(D\left(F\right)\) is a two-dimensional \(\mathbf{0}\)–centred disc of radius \(\mu R_{N}\). In both cases above, \(\rho>0\). The name “projected gradient” refers to the above as resembling a gradient projection formula for an optimization problem, where \(\mathbf{U}\) expresses a derivative of an objective function with respect to \(\mathbf{R}\).

De Saxcé and Feng formulation

This is an implicit formulation based on 2 and it is optionally used on an individual contact point level within the Gauss–Seidel solver. On the point level this formulation is implemented in dbs.c:96. We express the Signorini–Coulomb law (74) and (75) as an inclusion. The friction cone \(K_{\alpha}\) is defined as

(78)\[K_{\alpha}=\left\{ \mathbf{R}_{\alpha}:\left\Vert \mathbf{R}_{\alpha T}\right\Vert \le\mu_{\alpha}R_{\alpha N},R_{\alpha N}\ge0\right\}\]

where \(\mu_{\alpha}\) is the coefficient of friction. It has been shown by De Saxcé and Feng 2, that the Signorini–Coulomb law can be expressed in a compact form

(79)\[\begin{split}-\left[\begin{array}{c} \mathbf{U}_{\alpha T}\\ \bar{U}_{\alpha N}+\mu_{\alpha}\left\Vert \mathbf{U}_{\alpha T}\right\Vert \end{array}\right]\in N_{K_{\alpha}}\left(\mathbf{R}_{\alpha}\right)\end{split}\]

where \(N_{K_{\alpha}}\) stands for the normal cone of the set \(K_{\alpha}\). For a convex set A the normal cone \(N_{A}\left(\mathbf{R}\right)\) at point \(\mathbf{R}\in A\) is defined as the set of all vectors \(\mathbf{V}\) such that \(\left\langle \mathbf{V},\mathbf{S}-\mathbf{R}\right\rangle \le0\) for all \(\mathbf{S}\in A\). Based on inclusion (79), the authors of 2 propose the following projection formula

(80)\[\begin{split}\mathbf{R}_{\alpha}=\mbox{proj}_{K_{\alpha}}\left(\mathbf{R}_{\alpha}-\rho\left[\begin{array}{c} \mathbf{U}_{\alpha T}\\ \bar{U}_{\alpha N}+\mu_{\alpha}\left\Vert \mathbf{U}_{\alpha T}\right\Vert \end{array}\right]\right)\end{split}\]

where \(\rho>0\). Formula (80) can be used instead of the projected gradient formulas (76) and (77). The appeal of this approach is in the separation of velocities on the left hand side of the inclusion (79) from forces on the right hand side, as well as in the constancy of the friction cone \(K_{\alpha}\), which together make this formulation seem even more like a statement of optimality for a constrained optimization problem. This may be helpful in formulating solution strategies based on already existing approaches.

Non–smooth force equation formulation

This is an implicit formulation based on 3 and it is used by default on an individual contact point level within the Gauss–Seidel solver. On the point level this formulation is implemented in dbs.c:142. The authors of 3 propose to express the Signorini and Coulomb conditions (74) and (75) as a non–smooth equation \(\mathbf{C}\left(\mathbf{U},\mathbf{R}\right)=\mathbf{0}\), where

(81)\[\begin{split}\mathbf{C}\left(\mathbf{U},\mathbf{R}\right)=\left[\begin{array}{c} \max\left(\mu d_{N},\left\Vert \mathbf{d}_{T}\right\Vert \right)\mathbf{R}_{T}-\mu\max\left(0,d_{N}\right)\mathbf{d}_{T}\\ R_{N}-\max\left(0,d_{N}\right) \end{array}\right]\end{split}\]

and

(82)\[d_{N}=R_{N}-\rho\bar{U}_{N}\]

(83)\[\mathbf{d}_{T}=\mathbf{R}_{T}-\rho\mathbf{U}_{T}\]

while \(\rho>0\). Equation (81) encapsulates the projection formulas (76) and (77) and it has been shown to work well as a basis for Newton type solution schemes in the finite–element context.

Non–smooth velocity equation formulation

This is an implicit formulation developed specifically for Solfec-1.0 based on formula (79) from 2. It is optionally used on an individual contact point level within the Gauss–Seidel solver. It is also the basis of contact linearization within the projected Newton solver. On the point level this formulation is implemented in scf.c. Using the local dynamics relationship

(84)\[\mathbf{U_{\alpha}}=\mathbf{B_{\alpha}}+\sum_{\beta}\mathbf{W}_{\alpha\beta}\mathbf{R}_{\beta}\]

let us define a function

(85)\[\begin{split}\mathbf{F}\left(\mathbf{R}\right)=\left[\begin{array}{c} ...\\ \mathbf{U}_{\alpha T}\left(\mathbf{R}\right)\\ \bar{U}_{\alpha N}\left(\mathbf{R}\right)+\mu_{\alpha}\left\Vert \mathbf{U}_{\alpha T}\left(\mathbf{R}\right)\right\Vert \\ ... \end{array}\right]\end{split}\]

and a total cone

\[K=\bigcup_{\alpha}K_{\alpha}\]

where \(\mu_{\alpha}\) is the coefficient of friction at a contact point \(\alpha\), \(K_{\alpha}\) is the corresponding friction cone (78), while the dependence \(\mathbf{U}_{\alpha}\left(\mathbf{R}\right)\) is given by (84). Formula (79) states, that the frictional contact constraints are satisfied if \(-\mathbf{F}\left(\mathbf{R}\right)\) belongs to the normal cone of the friction cone at \(\mathbf{R}\). Hence

\[-\mathbf{F}\left(\mathbf{R}\right)=\mathbf{R}-\mathbf{F}\left(\mathbf{R}\right)- \mbox{proj}_{K}\left(\mathbf{R}-\mathbf{F}\left(\mathbf{R}\right)\right)\]

which can be reduced to the usual projection formula \(\mathbf{R}=\mbox{proj}_{K}\left(\mathbf{R}-\mathbf{F}\left(\mathbf{R}\right)\right)\) or (80) with \(\rho=1\). Let us not do it though, but rather define a vector field

\[\mathbf{m}\left(\mathbf{S}\right)=\mathbf{S}-\mbox{proj}_{K}\left(\mathbf{S}\right)= \mathbf{n}\left(\mathbf{S}\right)\left\langle \mathbf{n}\left(\mathbf{S}\right),\mathbf{S}\right\rangle\]

where

(86)\[\begin{split}\mathbf{n}_{\alpha}\left(\mathbf{S}_{\alpha}\right)=\left\{ \begin{array}{lll} \mathbf{0} & \mbox{if} & \left\Vert \mathbf{S}_{\alpha T}\right\Vert -\mu_{\alpha}S_{\alpha N}\le0\\ \mathbf{S}_{\alpha}/\left\Vert \mathbf{S}_{\alpha}\right\Vert & \mbox{if} & \mu_{\alpha}\left\Vert \mathbf{S}_{\alpha T}\right\Vert +S_{\alpha N}<0\\ \frac{1}{\sqrt{1+\mu_{\alpha}^{2}}}\left[\begin{array}{c} \mathbf{S}_{\alpha T}/\left\Vert \mathbf{S}_{\alpha T}\right\Vert \\ -\mu_{\alpha} \end{array}\right] & \mbox{} & \mbox{otherwise} \end{array}\right.\end{split}\]

We can rewrite (79) as

(87)\[\mathbf{C}\left(\mathbf{R}\right)=\mathbf{F}\left(\mathbf{R}\right)+\mathbf{m}\left(\mathbf{R}-\mathbf{F}\left(\mathbf{R}\right)\right)=\mathbf{0}\mbox{ and }\mathbf{R}\in K\]

Note, that \(\mathbf{F}\left(\mathbf{R}\right)\) is expressed in terms of velocity, and so is \(\mathbf{C}\left(\mathbf{R}\right)\). Equation (87) expresses, in velocity form, the projection formula (80).

Semi–explicit penalty formulation

This is a simple penalty based formulation developed specifically for Solfec-1.0 and used within the penalty solver. On the point level this formulation is implemented in pes.c. Let

\[s=spring\mbox{ and }d=dashpot\mbox{ and }g=gap\mbox{ and }m=hpow\]

where \(hpow\) stands for the “Hertz power”. The normal reaction is computed as follows

(88)\[R_{N}=-s\cdot\frac{g^{t+h}+g^{t}}{2}-d\cdot\frac{U_{N}^{t+h}+U_{N}^{t}}{2}\]

where \(U_{N}\) is the normal relative velocity. Recall, that the gap function is computed for the configuration \(\mathbf{q}^{t}+\frac{h}{2}\mathbf{u}^{t}\), so that the gap function value computed during geometrical contact detection reads

\[g=g^{t}+\frac{h}{2}U_{N}^{t}\]

We then have

\[g^{t+h}=g^{t}+\frac{h}{2}\left(U_{N}^{t+h}+U_{N}^{t}\right)=g+\frac{h}{2}U_{N}^{t+h}\]

and since \(g^{t}=g-\frac{h}{2}U_{N}^{t}\) we can estimate

(89)\[R_{N}=-s\cdot\left(g+\frac{h}{4}\left(U_{N}^{t+h}-U_{N}^{t}\right)\right)-\frac{d}{2}\cdot\left(U_{N}^{t+h}+U_{N}^{t}\right)\]

We then use the diagonal block of local dynamics

\[\mathbf{U}^{t+h}=\mathbf{B}+\mathbf{W}\mathbf{R}\]

in order to estimate \(U_{N}^{t+h}\) as follows

\[U_{N}^{t+h}=B_{N}+\mathbf{W}_{NT}\mathbf{R}_{T}+W_{NN}R_{N}\]

where a previous tangential reaction \(\mathbf{R}_{T}\) is employed. Inserting this it into (89) results in

\[\bar{B}_{N}=B_{N}+\mathbf{W}_{NT}\mathbf{R}_{T}\]

(90)\[R_{N}=\left[-s\cdot\left(g+\frac{h}{4}\left(\bar{B}_{N}-U_{N}^{t}\right)\right)-\frac{d}{2}\cdot\left( \bar{B}_{N}+U_{N}^{t}\right)\right]/\left[1+\left(s\cdot\frac{h}{4}+\frac{d}{2}\right)\cdot W_{NN}\right]\]

The reason for using the above, rather than the classical \(R_{N}=-s\cdot g-d\cdot U_{N}^{t}\) is in an increased stability of the this approach. Since we aim at simplicity and want to avoid any nonlinear solve only at this stage we include the “Hertz power”

\[g_{1}=\mbox{min}\left(g+\frac{h}{4}\left(\bar{B}_{N}-U_{N}^{t}\right),0\right)\]

\[s_{1}=sm\left(-g_{1}\right)^{m-1}\]

\[R_{N}=\left[s\cdot\left(-g_{1}\right)^{m}-\frac{d}{2}\cdot\left(\bar{B}_{N}+U_{N}^{t}\right)\right] /\left[1+\left(s_{1}\cdot\frac{h}{4}+\frac{d}{2}\right)\cdot W_{NN}\right]\]

Again aiming at maximum simplicity and assuming \(\mathbf{U}_{T}^{t+h}=0\) we then estimate the tangential stick reaction

\[\mathbf{R}_{T}=-\mathbf{W}_{TT}^{-1}\left(\mathbf{B}_{T}+\mathbf{W}_{TN}R_{N}\right)\]

The complete interface law is expressed the below algorithm (where \(h\) is the time step, \(g\) is the contact gap, \(s\) is the spring constant, \(d\) is the damper constant, \(\mu\) refers there to the coefficient of friction, and \(m\) is the “Hertz power”).

SPRING_DASHPOT_CONTACT \(\left(h,g,s,d,\mu,cohesion,cohesive\right)\)
1 \(\,\,\) \(\bar{B}_{N}=B_{N}+\mathbf{W}_{NT}\mathbf{R}_{T}\)
2 \(\,\,\) if semi–explicit then
3 \(\,\,\,\,\,\,\) \(g_{1}=\mbox{min}\left(g+\frac{h}{4}\left(\bar{B}_{N}-U_{N}^{t}\right),0\right)\)
4 \(\,\,\,\,\,\,\) \(s_{1}=sm\left(-g_{1}\right)^{m-1}\)
5 \(\,\,\,\,\,\,\) \(R_{N}=\left[s\cdot\left(-g_{1}\right)^{m}-\frac{d}{2}\cdot\left(\bar{B}_{N}+U_{N}^{t}\right)\right]/\left[1+\left(s_{1}\cdot\frac{h}{4}+\frac{d}{2}\right)\cdot W_{NN}\right]\)
6 \(\,\,\) else \(R_{N}=s\cdot\left(-\min\left(g,0\right)\right)^{m}-d\cdot U_{N}^{t}\)
7 \(\,\,\) if not \(cohesive\) and \(R_{N}<0\) then \(\mathbf{R}=0\) return
8 \(\,\,\) \(\mathbf{R}_{T}=-\mathbf{W}_{TT}^{-1}\left(\mathbf{B}_{T}+\mathbf{W}_{TN}R_{N}\right)\)
9 \(\,\,\) if \(cohesive\) and \(R_{N}<-cohesion\) then \(cohesive=false\) and \(R_{N}=-cohesion\)
10 \(\,\) if \(\left\Vert \mathbf{R}_{T}\right\Vert >\mu\left|R_{N}\right|\) then
11 \(\,\,\,\,\,\) \(\mathbf{R}_{T}=\mu R_{N}\mathbf{R}_{T}/\left\Vert \mathbf{R}_{T}\right\Vert\)
12 \(\,\,\,\,\,\) if \(cohesive\) then \(cohesive=false\)

1

P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Computer Methods in Applied Mechanics and Engineering, 92 (3), 353-375, 1991.

2(1,2,3,4)

G. De Saxcé and Z. Q. Feng, The bipotential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms, Mathematical and Computer Modelling, 28, 225-245, 1998.

3(1,2)

S. Hüeber, G. Stadler, and B. I. Wohlmuth, A primal–dual active set algorithm for three–dimensional contact problems with Coulomb friction, SIAM Journal on Scientific Computing, 30 (2), 572-596, 2007.