Tangent contact stiffness

Figure 132: Visualization of the tangential differential displacements

To find the tangent contact stiffness matrix, please look at Figure 132, part a). At the beginning of a concrete time increment, characterized by time $ t_n$, the slave node at position $ \boldsymbol{p}_n$ corresponds to the projection vector $ \boldsymbol{q}_n$ on the master side. At the end of the time increment, characterized by time $ t_{n+1}$ both have moved to positions $ \boldsymbol{p}_{n+1}$ and $ \boldsymbol{q}_{n+1}$, respectively. The differential displacement between slave and master surface changed during this increment by the vector $ \boldsymbol{s}$ satisfying:

$\displaystyle \boldsymbol{s}=(\boldsymbol{p}_{n+1} - \boldsymbol{q}_{n+1}) -(\boldsymbol{p}_n - \boldsymbol{q}_n).$ (228)

Here, $ \boldsymbol{q_{n+1}}$ satisfies

$\displaystyle \boldsymbol{q}_{n+1}= \sum_j \varphi_j(\xi^m _n, \eta^m_n) \boldsymbol{q_j}_{n+1} .$ (229)

Since (the dependency of $ \varphi_j$ on $ \xi$ and $ \eta$ is dropped to simplify the notation)

$\displaystyle \boldsymbol{p}_{n}= \boldsymbol{X} + \boldsymbol{u}_{n},$ (230)

$\displaystyle \boldsymbol{p}_{n+1}= \boldsymbol{X} + \boldsymbol{u}_{n+1},$ (231)

$\displaystyle \boldsymbol{q}_{n}= \sum_j \varphi_j \left[ \boldsymbol{X_j} + (\boldsymbol{u_j})_{n} \right ],$ (232)

$\displaystyle \boldsymbol{q}_{n+1}= \sum_j \varphi_j \left[ \boldsymbol{X_j} + (\boldsymbol{u_j})_{n+1} \right ],$ (233)

this also amounts to

$\displaystyle \boldsymbol{s}= \boldsymbol{u}_{n+1} - \sum_j \varphi_j (\boldsym...
...- \left [ \boldsymbol{u}_{n} - \sum_j \varphi_j (\boldsymbol{u_j})_{n} \right].$ (234)

Notice that the local coordinates take the values of time $ t_n$ (the superscript m denotes iteration m within the increment). The differential tangential displacement now amounts to:

$\displaystyle \boldsymbol{t} \equiv \boldsymbol{t_{n+1}}=\boldsymbol{t_{n}} + \boldsymbol{s}- s \boldsymbol{n},$ (235)


$\displaystyle s=\boldsymbol{s} \cdot \boldsymbol{n}.$ (236)

Derivation w.r.t. $ \boldsymbol{u_i}$ satisfies (straightforward differentiation):

$\displaystyle \frac{\partial \boldsymbol{t} }{\partial \boldsymbol{u_i} } = \fr...
...ldsymbol{u_i} } - s \frac{\partial \boldsymbol{n} }{\partial \boldsymbol{u_i} }$ (237)

$\displaystyle \frac{\partial s}{\partial \boldsymbol{u_i} } = \frac{\boldsymbol...
...bol{m}\Vert } \cdot \frac{\partial \boldsymbol{s} }{\partial \boldsymbol{u_i} }$ (238)


$\displaystyle \frac{\partial \boldsymbol{n} }{\partial \boldsymbol{u_i} } = \fr...
...bol{m } \otimes \frac{\partial \Vert \boldsymbol{m} \Vert }{\boldsymbol{u_i} }.$ (239)

The derivative of $ \boldsymbol{s}$ w.r.t. $ \boldsymbol{u_i}$ can be obtained from the derivative of $ \boldsymbol{r}$ w.r.t. $ \boldsymbol{u_i}$ by keeping $ \xi$ and $ \eta$ fixed (notice that the derivative is taken at $ t_{n+1}$, consequently, all derivatives of values at time $ t_n$ disappear):

$\displaystyle \frac{\partial \boldsymbol{s} }{\partial \boldsymbol{u_i} } = \delta _{ip} \boldsymbol{I} - \varphi_i (1 - \delta _{ip}) \boldsymbol{I}.$ (240)

Physically, the tangential contact equations are as follows (written at the location of slave node p):

First, a difference form of the additive decomposition of the differential tangential displacement is derived. Starting from

$\displaystyle \boldsymbol{t}= \boldsymbol{t^e} + \boldsymbol{t^p},$ (245)

one obtains after taking the time derivative:

$\displaystyle \boldsymbol{\dot{t}}= \boldsymbol{\dot{t}^e} + \boldsymbol{\dot{t}^p}.$ (246)

Substituting the slip evolution equation leads to:

$\displaystyle \dot{\gamma} \frac{\boldsymbol{F_T} }{\Vert \boldsymbol{F_T} \Vert }=\boldsymbol{\dot{t}} - \boldsymbol{\dot{t}^e},$ (247)

and after multiplying with $ K_t$:

$\displaystyle K_t \dot{\gamma} \frac{\boldsymbol{F_T} }{\Vert \boldsymbol{F_T} \Vert }=K_t \boldsymbol{\dot{t}} - K_t \boldsymbol{\dot{t}^e}$ (248)

Writing this equation at $ t_{n+1}$, using finite differences (backwards Euler), one gets:

$\displaystyle K_t \Delta \gamma_{n+1} \frac{\boldsymbol{F_T}_{n+1} }{\Vert\bold...
...lta \boldsymbol{t}_{n+1} - K_t \boldsymbol{t^e}_{n+1} + K_t \boldsymbol{t^e}_n,$ (249)

where $ \Delta \gamma_{n+1} \equiv \dot{\gamma} _{n+1} \Delta t$ and $ \Delta \boldsymbol{t}_{n+1} \equiv \boldsymbol{t}_{n+1} -
\boldsymbol{t}_n$. The parameter $ K_t$ is assumed to be independent of time.

Now, the radial return algorithm will be described to solve the governing equations. Assume that the solution at time $ t_n$ is known, i.e. $ \boldsymbol{t^e}_n$ and $ \boldsymbol{t^p}_n$ are known. Using the stick law the tangential forc $ \boldsymbol{F_T}_n$ can be calculated. Now we would like to know these variables at time $ t_{n+1}$, given the total differential tangential displacement $ \boldsymbol{t_{n+1}}$. At first we construct a trial tangential force $ \boldsymbol{F_T}_{n+1}^{trial}$ which is the force which arises at time $ t_{n+1}$ assuming that no slip takes place from $ t_n$ till $ t_{n+1}$. This assumption is equivalent to $ \boldsymbol{t^p}_{n+1}=\boldsymbol{t^p}_n$. Therefore, the trial tangential force satisfies (cf. the stick law):

$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} = K_t(\boldsymbol{t}_{n+1} - \boldsymbol{t^p}_n).$ (250)

Now, this can also be written as:

$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} = K_t(\boldsymbol{t}_{n+1} - \boldsymbol{t}_n + \boldsymbol{t}_n - \boldsymbol{t^p}_n).$ (251)


$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} = K_t \Delta \boldsymbol{t}_{n+1} + K_t \boldsymbol{t^e}_n.$ (252)

Using Equation (249) this is equivalent to:

$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} = K_t \Delta \gamma_{n+1} \frac{\b...
...l{F_T}_{n+1} }{\Vert\boldsymbol{F_T}_{n+1} \Vert} + K_t \boldsymbol{t^e}_{n+1},$ (253)


$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} = (K_t \Delta \gamma_{n+1} \frac{1 }{\Vert\boldsymbol{F_T}_{n+1} \Vert} + 1) \boldsymbol{F_T}_{n+1} .$ (254)

From the last equation one obtains

$\displaystyle \boldsymbol{F_T}_{n+1}^{trial} \parallel \boldsymbol{F_T}_{n+1}$ (255)

and, since the terms in brackets in Equation (254) are both positive:

$\displaystyle \Vert \boldsymbol{F_T}_{n+1}^{trial}\Vert = K_t \Delta \gamma_{n+1} + \Vert \boldsymbol{F_T}_{n+1} \Vert.$ (256)

The only equation which is left to be satisfied is the Coulomb slip limit. Two possibilities arise:

  1. $ \Vert \boldsymbol{F_T}_{n+1}^{trial}\Vert \le \mu \Vert
\boldsymbol{F_N}_{n+1} \Vert$.

    In that case the Coulomb slip limit is satisfied and we have found the solution:

    $\displaystyle \boldsymbol{F_T}_{n+1} = \boldsymbol{F_T}_{n+1}^{trial} = K_t( \boldsymbol{t_{n+1}} - \boldsymbol{t^p}_n)$ (257)


    $\displaystyle \frac{\partial \boldsymbol{F_T}_{n+1}}{\partial \boldsymbol{t}_{n+1} } = K_t \boldsymbol{I}.$ (258)

    No extra slip occurred from $ t_n$ to $ t_{n+1}$.

  2. $ \Vert \boldsymbol{F_T}_{n+1}^{trial}\Vert > \mu \Vert
\boldsymbol{F_N}_{n+1} \Vert$.

    In that case we project the solution back onto the slip surface and require $ \Vert \boldsymbol{F_T}_{n+1}\Vert = \mu \Vert
\boldsymbol{F_N}_{n+1} \Vert$. Using Equation (256) this leads to the following expression for the increase of the consistency parameter $ \gamma$:

    $\displaystyle \Delta \gamma _{n+1}= \frac{\Vert\boldsymbol{F_T}_{n+1}^{trial}\Vert- \mu \Vert \boldsymbol{F_N}_{n+1} \Vert }{K_t},$ (259)

    which can be used to update $ \boldsymbol{t^p}$ (by using the slip evolution equation):

    $\displaystyle \Delta \boldsymbol{t^p}= \Delta \gamma_{n+1} \frac{\boldsymbol{F_...
...rac{\boldsymbol{F_T}_{n+1}^{trial} }{\Vert\boldsymbol{F_T}_{n+1}^{trial} \Vert}$ (260)

    The tangential force can be written as:

    $\displaystyle \boldsymbol{F_T}_{n+1}= \Vert \boldsymbol{F_T}_{n+1} \Vert \frac{...
...c{\boldsymbol{F_T}_{n+1}^{trial}}{\Vert \boldsymbol{F_T}_{n+1}^{trial} \Vert }.$ (261)

    Now since

    $\displaystyle \frac{\partial \Vert \boldsymbol{a} \Vert}{\partial \boldsymbol{b...
...ymbol{a}\Vert } \cdot \frac{\partial \boldsymbol{a} }{\partial \boldsymbol{b} }$ (262)


    $\displaystyle \frac{\partial }{\partial \boldsymbol{b} } \left (\frac{\boldsymb...
...ght ) \right ] \cdot \frac{\partial \boldsymbol{a} }{\partial \boldsymbol{b} },$ (263)

    where $ \boldsymbol{a}$ and $ \boldsymbol{b}$ are vectors, one obtains for the derivative of the tangential force:

    $\displaystyle \frac{\partial \boldsymbol{F_T}_{n+1} }{\partial \boldsymbol{t}_{n+1} }$ $\displaystyle = \mu \boldsymbol{\xi }_{n+1} \otimes \left [ \frac{\boldsymbol{F...
...t \frac{\partial \boldsymbol{F_N}_{n+1}}{\partial \boldsymbol{t}_{n+1}} \right]$    
      $\displaystyle + \mu \frac{\Vert\boldsymbol{F_N}_{n+1} \Vert}{\Vert\boldsymbol{F...
...frac{\partial \boldsymbol{F_T}_{n+1}^{trial} }{\partial \boldsymbol{t}_{n+1}, }$ (264)


    $\displaystyle \boldsymbol{\xi }_{n+1} \equiv \frac{\boldsymbol{F_T}_{n+1}^{trial}}{\Vert \boldsymbol{F_T}_{n+1}^{trial} \Vert }.$ (265)

    One finally arrives at (using Equation (258)

    $\displaystyle \frac{\partial \boldsymbol{F_T}_{n+1} }{\partial \boldsymbol{u_i}_{n+1} }$ $\displaystyle = \mu \boldsymbol{\xi }_{n+1} \otimes \left [ -\boldsymbol{n} \cdot \frac{\partial \boldsymbol{F_N}_{n+1}}{\partial \boldsymbol{u_i}_{n+1}} \right]$    
      $\displaystyle + \mu \frac{\Vert\boldsymbol{F_N}_{n+1} \Vert}{\Vert\boldsymbol{F...
...ot K_t \frac{\partial \boldsymbol{t}_{n+1} }{\partial \boldsymbol{u_i}_{n+1}. }$ (266)

    All quantities on the right hand side are known now (cf. Equation (213) and Equation (237)).

    In CalculiX, for node-to-face contact, Equation (228) is reformulated and simplified. It is reformulated in the sense that $ \boldsymbol{q}_{n+1}$ is assumed to be the projection of $ \boldsymbol{p}_{n+1}$ and $ \boldsymbol{q}_n$ is written as (cf. Figure 132, part b))

    $\displaystyle \boldsymbol{q}_{n}= \sum_j \varphi_j(\xi^m _{n+1}, \eta^m_{n+1}) \boldsymbol{q_j}_{n} .$ (267)

    Part a) and part b) of the figure are really equivalent, they just represent the same facts from a different point of view. In part a) the projection on the master surface is performed at time $ t_n$, and the differential displacement is calculated at time $ t_{n+1}$, in part b) the projection is done at time $ t_{n+1}$ and the differential displacement is calculated at time $ t_n$. Now, the actual position can be written as the sum of the undeformed position and the deformation, i.e. $ \boldsymbol{p}=(\boldsymbol{X}+\boldsymbol{v})^s$ and $ \boldsymbol{q}=(\boldsymbol{X}+\boldsymbol{v})^m$ leading to:

    $\displaystyle \boldsymbol{s}=(\boldsymbol{X}+\boldsymbol{v})^s_{n+1} - (\boldsy...
...mbol{v})^s_n + (\boldsymbol{X}+\boldsymbol{v})^m_n(\xi^m _{n+1}, \eta^m_{n+1}).$ (268)

    Since the undeformed position is no function of time it drops out:

    $\displaystyle \boldsymbol{s}=\boldsymbol{v}^s_{n+1} - \boldsymbol{v}^m_{n+1}(\x...
...a^m_{n+1}) -\boldsymbol{v}^s_n + \boldsymbol{v}^m_n(\xi^m _{n+1}, \eta^m_{n+1})$ (269)


    $\displaystyle \boldsymbol{s}$ $\displaystyle =\boldsymbol{v}^s_{n+1} - \boldsymbol{v}^m_{n+1}(\xi^m _{n+1}, \e...
...n+1}) -\boldsymbol{v}^s_n + \boldsymbol{v}^m_n(\xi^m _{n}, \eta^m_{n})\noindent$ (270)
      $\displaystyle + \boldsymbol{v}^m_n(\xi^m _{n+1}, \eta^m_{n+1})- \boldsymbol{v}^m_n(\xi^m _{n}, \eta^m_{n})$ (271)

    Now, the last two terms are dropped, i.e. it is assumed that the differential deformation at time $ t_n$ between positions $ (\xi^m _{n},
\eta^m_{n})$ and $ (\xi^m _{n+1}, \eta^m_{n+1})$ is neglegible compared to the differential motion from $ t_n$ to $ t_{n+1}$. Then the expression for $ \boldsymbol{s}$ simplifies to:

    $\displaystyle \boldsymbol{s}=\boldsymbol{v}^s_{n+1} - \boldsymbol{v}^m_{n+1}(\x...
...\eta^m_{n+1}) -\boldsymbol{v}^s_n + \boldsymbol{v}^m_n(\xi^m _{n}, \eta^m_{n}),$ (272)

    and the only quantity to be stored is the difference in deformation between $ \boldsymbol{p} $ and $ \boldsymbol{q}$ at the actual time and at the time of the beginning of the increment.