Tangent contact stiffness

Figure: Visualization of the tangential differential displacements
\begin{figure}\epsfig{file=tangentcontact.eps,width=12cm}\end{figure}

To find the tangent contact stiffness matrix, please look at Figure [*], 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:

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

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

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

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

\boldsymbol{p}_{n}= \boldsymbol{X} + \boldsymbol{u}_{n}, (212)

\boldsymbol{p}_{n+1}= \boldsymbol{X} + \boldsymbol{u}_{n+1}, (213)

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

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

this also amounts to

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

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:

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

where

s=\boldsymbol{s} \cdot \boldsymbol{n}. (218)

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

\frac{\partial \boldsymbol{t} }{\partial \boldsymbol{u_i} } = \frac{\partial \boldsymbol{s} }{\partial \boldsymbol{u_i} } - \boldsymbol{n} \otimes\frac{\partial s}{\partial \boldsymbol{u_i} } - s \frac{\partial \boldsymbol{n}}{\partial \boldsymbol{u_i} } (219)

\frac{\partial s}{\partial \boldsymbol{u_i} } = \frac{\boldsymbol{s} }{\| \boldsymbol{m}\|} \cdot \frac{\partial \boldsymbol{m} }{\partial\boldsymbol{u_i} } - \frac{s}{\| \boldsymbol{m}\| } \frac{\partial \|\boldsymbol{m} \| }{\boldsymbol{u_i} } + \frac{\boldsymbol{m} }{\|\boldsymbol{m}\| } \cdot \frac{\partial \boldsymbol{s} }{\partial\boldsymbol{u_i} } (220)

and

\frac{\partial \boldsymbol{n} }{\partial \boldsymbol{u_i} } = \frac{1}{\| \boldsymbol{m}\| } \frac{\partial \boldsymbol{m} }{\partial \boldsymbol{u_i}}- \frac{1}{\| \boldsymbol{m}\|^2 } \boldsymbol{m } \otimes \frac{\partial \| \boldsymbol{m} \| }{\boldsymbol{u_i} }. (221)

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):

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

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

\boldsymbol{t}= \boldsymbol{t^e} + \boldsymbol{t^p}, (227)

one obtains after taking the time derivative:

\boldsymbol{\dot{t}}= \boldsymbol{\dot{t}^e} + \boldsymbol{\dot{t}^p}. (228)

Substituting the slip evolution equation leads to:

\dot{\gamma} \frac{\boldsymbol{F_T} }{\|\boldsymbol{F_T} \| }=\boldsymbol{\dot{t}} - \boldsymbol{\dot{t}^e}, (229)

and after multiplying with K_t:

K_t \dot{\gamma} \frac{\boldsymbol{F_T} }{\|\boldsymbol{F_T} \| }=K_t \boldsymbol{\dot{t}} - K_t \boldsymbol{\dot{t}^e} (230)

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

K_t \Delta \gamma_{n+1} \frac{\boldsymbol{F_T}_{n+1}}{\|\boldsymbol{F_T}_{n+1} \|} = K_t \Delta \boldsymbol{t}_{n+1} - K_t\boldsymbol{t^e}_{n+1} + K_t \boldsymbol{t^e}_n, (231)

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):

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

Now, this can also be written as:

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

or

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

Using Equation ([*]) this is equivalent to:

\boldsymbol{F_T}_{n+1}^{trial} = K_t \Delta \gamma_{n+1} \frac{\boldsymbol{F_T}_{n+1}}{\|\boldsymbol{F_T}_{n+1} \|} + K_t\boldsymbol{t^e}_{n+1}, (235)

or

\boldsymbol{F_T}_{n+1}^{trial} = (K_t \Delta \gamma_{n+1} \frac{1}{\|\boldsymbol{F_T}_{n+1} \|} + 1) \boldsymbol{F_T}_{n+1} . (236)

From the last equation one obtains

\boldsymbol{F_T}_{n+1}^{trial} \parallel \boldsymbol{F_T}_{n+1} (237)

and, since the terms in brackets in Equation ([*]) are both positive:

\| \boldsymbol{F_T}_{n+1}^{trial}\| = K_t \Delta \gamma_{n+1} + \|\boldsymbol{F_T}_{n+1} \|. (238)

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:

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

    and

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

    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 ([*]) this leads to the following expression for the increase of the consistency parameter \gamma:

    \Delta \gamma _{n+1}= \frac{\|\boldsymbol{F_T}_{n+1}^{trial}\|- \mu \|\boldsymbol{F_N}_{n+1} \| }{K_t}, (241)

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

    \Delta \boldsymbol{t^p}= \Delta \gamma_{n+1} \frac{\boldsymbol{F_T}_{n+1}}{\|\boldsymbol{F_T}_{n+1} \|} = \Delta \gamma_{n+1} \frac{\boldsymbol{F_T}_{n+1}^{trial}}{\|\boldsymbol{F_T}_{n+1}^{trial} \|} (242)

    The tangential force can be written as:

    \boldsymbol{F_T}_{n+1}= \| \boldsymbol{F_T}_{n+1} \|\frac{\boldsymbol{F_T}_{n+1}}{\| \boldsymbol{F_T}_{n+1} \| } = \mu \|\boldsymbol{F_N}_{n+1} \| \frac{\boldsymbol{F_T}_{n+1}^{trial}}{\| \boldsymbol{F_T}_{n+1}^{trial} \| }. (243)

    Now since

    \frac{\partial \| \boldsymbol{a} \|}{\partial \boldsymbol{b} } =\frac{\boldsymbol{a} }{\| \boldsymbol{a}\| } \cdot \frac{\partial \boldsymbol{a} }{\partial \boldsymbol{b} } (244)

    and

    \frac{\partial }{\partial \boldsymbol{b} } \left (\frac{\boldsymbol{a} }{\| \boldsymbol{a} \|} \right ) = \frac{1}{\| \boldsymbol{a}\| } \left [\boldsymbol{I} - \left (\frac{\boldsymbol{a} }{\| \boldsymbol{a} \|} \right ) \otimes \left (\frac{\boldsymbol{a} }{\| \boldsymbol{a} \|} \right ) \right ] \cdot \frac{\partial \boldsymbol{a} }{\partial \boldsymbol{b} }, (245)

    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}, } (246)

    where

    \boldsymbol{\xi }_{n+1} \equiv \frac{\boldsymbol{F_T}_{n+1}^{trial}}{\| \boldsymbol{F_T}_{n+1}^{trial} \| }. (247)

    One finally arrives at (using Equation ([*])

    \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}. } (248)

    All quantities on the right hand side are known now (cf. Equation ([*]) and Equation ([*])).

    In CalculiX, for node-to-face contact, Equation ([*]) 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 [*], part b))

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

    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:

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

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

    \boldsymbol{s}=\boldsymbol{v}^s_{n+1} - \boldsymbol{v}^m_{n+1}(\xi^m _{n+1}, \eta^m_{n+1}) -\boldsymbol{v}^s_n+ \boldsymbol{v}^m_n(\xi^m _{n+1}, \eta^m_{n+1}) (251)

    or:

    \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 (252)
      \displaystyle + \boldsymbol{v}^m_n(\xi^m _{n+1}, \eta^m_{n+1})- \boldsymbol{v}^m_n(\xi^m _{n}, \eta^m_{n}) (253)

    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:

    \boldsymbol{s}=\boldsymbol{v}^s_{n+1} - \boldsymbol{v}^m_{n+1}(\xi^m _{n+1}, \eta^m_{n+1}) -\boldsymbol{v}^s_n+ \boldsymbol{v}^m_n(\xi^m _{n}, \eta^m_{n}), (254)

    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.