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):
\boldsymbol{t}= \boldsymbol{t^e} + \boldsymbol{t^p}. | (223) |
\boldsymbol{F_T} = K_t \boldsymbol{t^e}. | (224) |
\| \boldsymbol{F_T} \| \le \mu \| \boldsymbol{F_N} \| | (225) |
\boldsymbol{\boldsymbol{\dot{t}^p} }= \dot{\gamma} \frac{\boldsymbol{F_T} }{\|\boldsymbol{F_T} \| } | (226) |
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:
The only equation which is left to be satisfied is the Coulomb slip limit. Two possibilities arise:
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
No extra slip occurred from t_n to t_{n+1}.
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.