(Difference between revisions)
 Revision as of 14:10, 2 May 2008 (view source) (→Initialisation of Contact Partners)← Older edit Latest revision as of 19:24, 1 November 2008 (view source)JMora (Talk | contribs) (50 intermediate revisions by one user not shown) Line 2: Line 2: This conditions represents the main part of the algorithm to solve contact problems. The algorithm is based on a surface to surface augmented lagrange/updated penalty formulation using a master and a slave contact surface. It is formulated to solve frictional contact problems in structural mechanics, the algorithm is formulated in a finite deformation context. This conditions represents the main part of the algorithm to solve contact problems. The algorithm is based on a surface to surface augmented lagrange/updated penalty formulation using a master and a slave contact surface. It is formulated to solve frictional contact problems in structural mechanics, the algorithm is formulated in a finite deformation context. {| border="1" cellpadding="5" cellspacing="0" align="center" {| border="1" cellpadding="5" cellspacing="0" align="center" − ! style="background:#ffdead;" | This condition and the underlying algorithm for contact problems has not been tested by means of a broader range of benchmark examples so far. It seem tio work quite fine for a number of different examples, but I do not give a guarantee that it is free of errors in its formulation or its implementation. + ! style="background:#ffdead;" | This condition and the underlying algorithm for contact problems has not yet been tested by means of a broader range of benchmark examples. It seems to work quite fine for a number of different examples, but I do not give any guarantee that it is free of errors in its formulation or its implementation. |} |} Line 23: Line 23: The constrained that defines the contact problem is now stated as: The constrained that defines the contact problem is now stated as: *$g(x)\leq 0$ *$g(x)\leq 0$ − Within the finite element formulation the problem is now solved by use of an updated penalty treatment of the contact constrained. To solve the problem the Energie of the system must be minimised by fulfilling the penetration constraint. After having introduced a contact stress + Within the finite element formulation the problem can be solved by use of an updated penalty treatment of the contact constrained. To solve the problem the energy of the system must be minimised by fulfilling the penetration constraint. After having introduced a contact stress *$t_N =<\lambda_N+\epsilon_Ng>$ *$t_N =<\lambda_N+\epsilon_Ng>$ − that acts on the surfaces if they are in contact the constraint can be rewritten in form of Kuhn-Tucker constraints: + that acts on the surfaces if they are in contact, the constraint can be rewritten in form of Kuhn-Tucker constraints: *$g\leq 0\;\;\;\;\;\;t_N\geq 0\;\;\;\;\;\;g t_N=0\;\;\;\forall\;\mathbf{x}\epsilon\Gamma^c$ *$g\leq 0\;\;\;\;\;\;t_N\geq 0\;\;\;\;\;\;g t_N=0\;\;\;\forall\;\mathbf{x}\epsilon\Gamma^c$ Whereas the contact stress $t_N$ consists of an (outer) Lagrange value $\lambda_N$, which is constant during the Newton iteration, and an inner penaltization of the gap $\epsilon_Ng$. Whereas the contact stress $t_N$ consists of an (outer) Lagrange value $\lambda_N$, which is constant during the Newton iteration, and an inner penaltization of the gap $\epsilon_Ng$. The updated penalty method now penalizes violation of the constraint, leading to the following potential: The updated penalty method now penalizes violation of the constraint, leading to the following potential: *$\Pi^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\frac{1}{2\epsilon_n}<\lambda_N+\epsilon_Ng>^2-\frac{1}{2\epsilon_n}\lambda_N^2d\Gamma^i$ *$\Pi^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\frac{1}{2\epsilon_n}<\lambda_N+\epsilon_Ng>^2-\frac{1}{2\epsilon_n}\lambda_N^2d\Gamma^i$ − The minimization of the potential is now been done by its variation: + which is minimized over by its variation: *$\delta W^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\delta g<\lambda_N+\epsilon_Ng>d\Gamma^i$ *$\delta W^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\delta g<\lambda_N+\epsilon_Ng>d\Gamma^i$ − Together with the minimisation of the system energy + Together with the minimization of the internal and external energies *$\delta W^{int, ext}=\sum_{i=1}^2\delta W_i^{int, ext}$ *$\delta W^{int, ext}=\sum_{i=1}^2\delta W_i^{int, ext}$ − the minimisation of the system energie with fulfillment of the penetration constraint can be written as + the minimization of the system energie with fulfillment of the penetration constraint can be written as *$\delta W=\delta W^{int, ext}+\delta W^c=0$ *$\delta W=\delta W^{int, ext}+\delta W^c=0$ This equation states a initial boundary value problem and can be solved by normal discretization methods in time and space. This equation states a initial boundary value problem and can be solved by normal discretization methods in time and space. Line 44: Line 44: ** next Uzawa iteration $uz\rightarrow uz+1$ ** next Uzawa iteration $uz\rightarrow uz+1$ The use of an updated penalty formulation allows to solve the problem exact even by using a (small) penalty value. It therefore combines the strength of the penalty method regarding simplicity with the accuracy of the Lagrange method. Nevertheless to increase the accuracy the penalized problem has to be solved a couple of times. The use of an updated penalty formulation allows to solve the problem exact even by using a (small) penalty value. It therefore combines the strength of the penalty method regarding simplicity with the accuracy of the Lagrange method. Nevertheless to increase the accuracy the penalized problem has to be solved a couple of times. + === Formulation with friction === === Formulation with friction === − If now friction does occur between the contacting surfaces a constrained to for the used friction law has to be considered. In case of Coulomb friction this friction law is stated as: + If friction occurs between the contacting surfaces a constrained for the stick and the slip of the surfaces has to be considered according to the used friction law. In case of Coulomb friction this friction law is stated as: − $|\boldsymbol{v}_T|=\left\{\begin{array}{c}=0\;\; \mbox{if}\;\;t_T\leq\mu t_N\\\geq 0 \;\; \mbox{if}\;\; t_T>\mu t_N\end{array}\right.$. This constrained is now treated in a similar way like it is described above. + $|\boldsymbol{v}_T|=\left\{\begin{array}{c}=0\;\; \mbox{if}\;\;t_T\leq\mu t_N\\\geq 0 \;\; \mbox{if}\;\; t_T>\mu t_N\end{array}\right.$. This constrained is treated within the present condition in a similar way like it is described above for contact constrained in normal direction. + == Algorithmic Formulation == == Algorithmic Formulation == === Contact Search === === Contact Search === − Before the start of the contact simulation a set of possible contact partners is initialised. Therefore for every quadrature point on the slave surface a partner on the master surface is search, e.g. the point on the master surface which is the closest to the slave quadrature point: + Before the start of the contact simulation a set of possible contact partners is initialised. Therefore for every quadrature point on the slave surface a partner on the master surface is searched, e.g. the point on the master surface which is the closest to the slave quadrature point: − *$\mbox{min}|\nu*\left(\boldsymbol{X}^1-\boldsymbol{X}^2\right)|$ where $\nu$ is the normal on the master surface + *$\mbox{min}|\boldsymbol{n}^2*\left(\boldsymbol{X}^1-\boldsymbol{X}^2\right)|$ where $\boldsymbol{n}^2$ is the normal on the master surface this search is performed by an outer and an inner search this search is performed by an outer and an inner search *outer search: search the node that is closest to the slave quadrature point *outer search: search the node that is closest to the slave quadrature point *inner search: for every master facet belonging to this node now a closest point projection is done *inner search: for every master facet belonging to this node now a closest point projection is done + === Initialisation of Contact Partners === === Initialisation of Contact Partners === − If a closest point has been found the slave surface and the master surface are stored as a ContactLink3D, e.g. one of the contact pairs where the constraints has to be fulfilled. + If a closest point has been found the slave surface and the master surface are stored as a ContactLink3D, e.g. one of the contact pairs where the constraints have to be fulfilled. Whereas now all the computations to compute the RHS and LHS are done in the ContactLink3D. The Constructor of the ContactLink3D gets the following informations to set up the pair of contact points: + * NewId: The Id of the contact link, given by the ContactUtility + * pGeometry: The contact link does not have a special geometry + * pProperties: The model part properties + * Master: A pointer on the master facet + * Slave:  A pointer on the slave facet + * MasterContactLocalPoint: The local coordinates of the point on the master surface that is a closest point projection + * SlaveContactLocalPoint: The local coordinates of the quadrature point on the slave surface + * SlaveIntegrationPointIndex: The integration point index of the quadrature point on the slave surface === Uzawa loop === === Uzawa loop === + * First step: After the contact links have been initialized the Lagrangian multipliers $[\lambda_N\;\;\lambda_T^1\;\;\lambda_T^2]$ are set to their initial values. + * Next steps: After each Newton loop the Multipliers are updated: + **$\lambda_n^{uz+1}=<\lambda_n^{uz}+\epsilon_ng>$ + * Is converged: The contact links are deleted, the Lagrange multipliers defined on the quadrature points of the slave surfaces are stored to be the initial values for the multipliers in the next time step. + === Inner loop === === Inner loop === + The fully discretized weak form including the penaltization of the constraint and the Lagrange Multipliers is solved iteratively by use of Newton's method. + == Element Formulation == == Element Formulation == + === Normal contact force === + $t_N= <\lambda_N+\epsilon_Ng>$ with $g=-\boldsymbol{n}^2\cdot\left(\boldsymbol{x}_1-\boldsymbol{x}_2\right)$ + + === Frictional contact force === + $t_T=\left\{\begin{array}{cc}\mu t_N\frac{\boldsymbol{t}_T^{trial}}{||\boldsymbol{t}_T^{trial}||}&if\;\;||\boldsymbol{t}_T^{trial}||\geq \mu t_N\\\boldsymbol{t}_T^{trial}&if\;\;||\boldsymbol{t}_T^{trial}||< \mu t_N\end{array}\right.$ + + whereas $t_T^{trial}=\boldsymbol{\lambda}_T+\epsilon_T\boldsymbol{v}^{\overline{12}}_{rel}$, having the relative tangential velocity $\boldsymbol{v}^{\overline{12}}_{rel}=\left(\dot{\boldsymbol{u}}^1-\dot{\boldsymbol{u}}^2\right)\left[\boldsymbol{\tau}^2_1\;\;\;\boldsymbol{\tau}^2_2\right]$ with $\boldsymbol{\tau}^2_i$ as the tangential vectors on the master surface + + $||\boldsymbol{t}_T^{trial}||=\sqrt{t^{trial}_{T_i}m_{ij}t^{trial}_{T_j}}$ with $m_{ij}=\frac{\boldsymbol{\tau}^2_i\cdot\boldsymbol{\tau}^2_j}{|\boldsymbol{\tau}^2_i||\boldsymbol{\tau}^2_j|}$ + + === Inner load vector === + $\boldsymbol{r}^{int}=\underbrace{\int_{\Gamma^1}\delta g t_N\boldsymbol{n}^2dA}_{\boldsymbol{r}^{int}_{normal}}+\underbrace{\int_{\Gamma}\delta\boldsymbol{\zeta}\boldsymbol{t}_tdA}_{\boldsymbol{r}^{int}_{friction}}$ + ::$=\sum_q^{NQ^1}\left\{\left[\boldsymbol{N}^1\;\;\left(-\boldsymbol{N}^2\right)\right]^Tt_N\left(\xi^q_1\right)\gamma^q\boldsymbol{n}^2|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|\right\} + +\sum_q^{NQ^1}\left\{\left[\boldsymbol{N}^1\;\;\left(-\boldsymbol{N}^2\right)\right]^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right]\boldsymbol{t}_T\left(\xi_1^q\right)\gamma^q\boldsymbol{n}^2|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|\right\}$ + + + with $\boldsymbol{N}^i=\left[\begin{array}{ccc|c|ccc}N_1^i&0&0&\cdots&N_n^i&0&0\\0&N_1^i&0&\cdots&0&N_n^i&0\\0&0&N_1^i&\cdots&0&0&N_n^i\end{array}\right]$ + + === LeftHandSide-Contribution === + $\boldsymbol{K}=\left[\begin{array}{cc}\boldsymbol{K}^1_1&\boldsymbol{K}^1_2\\\boldsymbol{K}^2_1&\boldsymbol{K}^2_2\end{array}\right]$ where the indices refer to the variation (top) and the linearization (bottom) regarding the degrees of freedom of the slave (1) and the master (2) facet. + + ====Contributions due to a linearization of the virtual work associated with the normal contact force==== + + $\boldsymbol{K}^{1}_{1}=\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^{1}_{2}=-\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|-\sum_q^{NQ}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\epsilon_N\left(\frac{\delta\boldsymbol{n}^2}{\delta\boldsymbol{u}^2}\left(\boldsymbol{x}^1-\boldsymbol{x}^2\right)\right)\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + ::$-\sum_q^{NQ^1}\left[\boldsymbol{N}^1\right]^T\frac{\delta\boldsymbol{n}^2}{\delta \boldsymbol{u}^2}t_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^{2}_{1}=-\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + ::$+\sum_q^{NQ^1}\left[\boldsymbol{N}^1\right]^T\frac{\delta\boldsymbol{n}^2}{\delta \boldsymbol{u}^2}t_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^{2}_{2}=\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|+\sum_q^{NQ}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\epsilon_N\left(\frac{\delta\boldsymbol{n}^2}{\delta\boldsymbol{u}^2}\left(\boldsymbol{x}^1-\boldsymbol{x}^2\right)\right)\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + ====Contributions due to a linearization of the virtual work associated with the frictional contact force==== + *$||\boldsymbol{t}_T^{trial}||< \mu t_N$: + + $\boldsymbol{K}_2^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \left[\frac{\partial\left(\boldsymbol{t}_T^{trial}\right)}{\partial \boldsymbol{u}^2}\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| + +\sum_q^{NQ}\left(\boldsymbol{N^1}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + $\boldsymbol{K}_2^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \left[\frac{\partial\left(\boldsymbol{t}_T^{trial}\right)}{\partial \boldsymbol{u}^2}\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| + -\sum_q^{NQ}\left(\boldsymbol{N^2}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + *$||\boldsymbol{t}_T^{trial}||\geq \mu t_N$: + $\boldsymbol{K}_1^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^1}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + $\boldsymbol{K}_2^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^2}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| + +\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu t_N\left[\frac{\partial\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}}{\partial \boldsymbol{u}^2}\right] + \gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + ::$+\sum_q^{NQ}\left(\boldsymbol{N^1}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + $\boldsymbol{K}_1^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^1}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + + $\boldsymbol{K}_2^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^2}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| + -\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] + \mu t_N\left[\frac{\partial\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}}{\partial \boldsymbol{u}^2}\right] + \gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + ::$-\sum_q^{NQ}\left(\boldsymbol{N^2}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$ + + === Damp-Matrix Contributions === + Only for frictional contact problems. + *$||\boldsymbol{t}_T^{trial}||< \mu t_N$: + + $\boldsymbol{K}^1_1=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^1_2=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^2_1=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^2_2=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + *$||\boldsymbol{t}_T^{trial}||\geq \mu t_N$: + + $\boldsymbol{K}^1_1=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_1}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^1_2=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_2}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^2_1=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_1}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + + $\boldsymbol{K}^2_2=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_2}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$ + ==References== ==References== * T.A. Laursen ''Computational Contact and Impact Mechanics'', Springer, 2002 * T.A. Laursen ''Computational Contact and Impact Mechanics'', Springer, 2002 Line 66: Line 173: *[[SlaveContactFace3D]] *[[SlaveContactFace3D]] *[[ContactUtility]] *[[ContactUtility]] + + + [[Category:Conditions]]

## General description of the condition

This conditions represents the main part of the algorithm to solve contact problems. The algorithm is based on a surface to surface augmented lagrange/updated penalty formulation using a master and a slave contact surface. It is formulated to solve frictional contact problems in structural mechanics, the algorithm is formulated in a finite deformation context.

This condition and the underlying algorithm for contact problems has not yet been tested by means of a broader range of benchmark examples. It seems to work quite fine for a number of different examples, but I do not give any guarantee that it is free of errors in its formulation or its implementation.

## Theoretical background of the used formulation

### Formulation without friction

Two bodies Ωi staying in contact over a contact surface $\Gamma^c_i$
Slave Γ1 and master Γ2 contact surface, and gap $g=\nu\cdot\left(\boldsymbol{y}_2-\boldsymbol{y}_1\right)$ between the surfaces
Three dimensional view on the contact surfaces and the gap between the two surfaces

The contact of two (or more) bodies in three dimension is a constrained continuum mechanical problem. The surface $\partial\Omega =\Gamma$ of the contacting continua Ω1 and Ω2 are constructed by the Neumann boundary Γσ, the Dirichlet boundary Γu and the contact boundary Γc. For the treated area Ω the following assumptions hold:

• $\Omega = \Omega_1 \cup \Omega_2$ and $\Omega_1 \cap \Omega_2=0$,
• having the boundaries $\partial\Omega=\Gamma^{\sigma}_i\cup\Gamma^{u}_i\cup\Gamma^{c}_i$ and $\Gamma^{\sigma}_i\cap \Gamma^{u}_i=0$, $\Gamma^{\sigma}_i\cap \Gamma^{c}_i=0$, $\Gamma^{u}_i\cap \Gamma^{c}_i=0$

The problem is constrained by the prescription of a penetration of the two bodies, e.g. of the penetration of the master surface Γ2 by the slave surface Γ1. This constrained is expressed by the gap function defined as the scalar projection of the distance vector between the point on the slave surface $\mathbf{x}\in \Gamma^1$ and its closest point projection onto the slave surface $\mathbf{y}(\mathbf{x})=\mbox{arg}\;\mbox{min}_{\mathbf{y}\epsilon\Gamma_2^c}\| \mathbf{x}-\mathbf{y}\left(\mathbf{x}\right)\|$ and the normal vector ν on the master surface:

• $g\left(\mathbf{x}\right)=-\mathbf{\nu}\left(\mathbf{x}-\mathbf{y}\right)$
$=-\mathbf\nu\left(\left[\mathbf{X}+\mathbf{u}_1\right]-\left[\mathbf{Y}-\mathbf{u}_2\right]\right)$
$=g_0-\mathbf\nu\left(\mathbf{u}_1-\mathbf{u}_2\right)$

The constrained that defines the contact problem is now stated as:

• $g(x)\leq 0$

Within the finite element formulation the problem can be solved by use of an updated penalty treatment of the contact constrained. To solve the problem the energy of the system must be minimised by fulfilling the penetration constraint. After having introduced a contact stress

• tN = < λN + εNg >

that acts on the surfaces if they are in contact, the constraint can be rewritten in form of Kuhn-Tucker constraints:

• $g\leq 0\;\;\;\;\;\;t_N\geq 0\;\;\;\;\;\;g t_N=0\;\;\;\forall\;\mathbf{x}\epsilon\Gamma^c$

Whereas the contact stress tN consists of an (outer) Lagrange value λN, which is constant during the Newton iteration, and an inner penaltization of the gap εNg. The updated penalty method now penalizes violation of the constraint, leading to the following potential:

• $\Pi^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\frac{1}{2\epsilon_n}<\lambda_N+\epsilon_Ng>^2-\frac{1}{2\epsilon_n}\lambda_N^2d\Gamma^i$

which is minimized over by its variation:

• $\delta W^c=\sum_{i=1}^2\int_{\Gamma_c^{i}}\delta g<\lambda_N+\epsilon_Ng>d\Gamma^i$

Together with the minimization of the internal and external energies

• $\delta W^{int, ext}=\sum_{i=1}^2\delta W_i^{int, ext}$

the minimization of the system energie with fulfillment of the penetration constraint can be written as

• δW = δWint,ext + δWc = 0

This equation states a initial boundary value problem and can be solved by normal discretization methods in time and space. The algorithm for the solution of the contact problem now consists of an inner iteration (Newton) and an outer iteration (Uzawa):

• start simulation with λ = 0
• Uzawa loop
• Solve $\delta W=\delta W^{int, ext}+\delta W^c\left(\lambda^{uz}_N=\mbox{const.}, \epsilon_N\right)=0$ by using Newtons loop
• update the Lagrange Multiplier $\lambda_N^{uz+1}=<\lambda_N^{uz}+\epsilon_Ng>$
• next Uzawa iteration $uz\rightarrow uz+1$

The use of an updated penalty formulation allows to solve the problem exact even by using a (small) penalty value. It therefore combines the strength of the penalty method regarding simplicity with the accuracy of the Lagrange method. Nevertheless to increase the accuracy the penalized problem has to be solved a couple of times.

### Formulation with friction

If friction occurs between the contacting surfaces a constrained for the stick and the slip of the surfaces has to be considered according to the used friction law. In case of Coulomb friction this friction law is stated as: $|\boldsymbol{v}_T|=\left\{\begin{array}{c}=0\;\; \mbox{if}\;\;t_T\leq\mu t_N\\\geq 0 \;\; \mbox{if}\;\; t_T>\mu t_N\end{array}\right.$. This constrained is treated within the present condition in a similar way like it is described above for contact constrained in normal direction.

## Algorithmic Formulation

### Contact Search

Before the start of the contact simulation a set of possible contact partners is initialised. Therefore for every quadrature point on the slave surface a partner on the master surface is searched, e.g. the point on the master surface which is the closest to the slave quadrature point:

• $\mbox{min}|\boldsymbol{n}^2*\left(\boldsymbol{X}^1-\boldsymbol{X}^2\right)|$ where $\boldsymbol{n}^2$ is the normal on the master surface

this search is performed by an outer and an inner search

• outer search: search the node that is closest to the slave quadrature point
• inner search: for every master facet belonging to this node now a closest point projection is done

### Initialisation of Contact Partners

If a closest point has been found the slave surface and the master surface are stored as a ContactLink3D, e.g. one of the contact pairs where the constraints have to be fulfilled. Whereas now all the computations to compute the RHS and LHS are done in the ContactLink3D. The Constructor of the ContactLink3D gets the following informations to set up the pair of contact points:

• NewId: The Id of the contact link, given by the ContactUtility
• pGeometry: The contact link does not have a special geometry
• pProperties: The model part properties
• Master: A pointer on the master facet
• Slave: A pointer on the slave facet
• MasterContactLocalPoint: The local coordinates of the point on the master surface that is a closest point projection
• SlaveContactLocalPoint: The local coordinates of the quadrature point on the slave surface
• SlaveIntegrationPointIndex: The integration point index of the quadrature point on the slave surface

### Uzawa loop

• First step: After the contact links have been initialized the Lagrangian multipliers $[\lambda_N\;\;\lambda_T^1\;\;\lambda_T^2]$ are set to their initial values.
• Next steps: After each Newton loop the Multipliers are updated:
• $\lambda_n^{uz+1}=<\lambda_n^{uz}+\epsilon_ng>$
• Is converged: The contact links are deleted, the Lagrange multipliers defined on the quadrature points of the slave surfaces are stored to be the initial values for the multipliers in the next time step.

### Inner loop

The fully discretized weak form including the penaltization of the constraint and the Lagrange Multipliers is solved iteratively by use of Newton's method.

## Element Formulation

### Normal contact force

tN = < λN + εNg > with $g=-\boldsymbol{n}^2\cdot\left(\boldsymbol{x}_1-\boldsymbol{x}_2\right)$

### Frictional contact force

$t_T=\left\{\begin{array}{cc}\mu t_N\frac{\boldsymbol{t}_T^{trial}}{||\boldsymbol{t}_T^{trial}||}&if\;\;||\boldsymbol{t}_T^{trial}||\geq \mu t_N\\\boldsymbol{t}_T^{trial}&if\;\;||\boldsymbol{t}_T^{trial}||< \mu t_N\end{array}\right.$

whereas $t_T^{trial}=\boldsymbol{\lambda}_T+\epsilon_T\boldsymbol{v}^{\overline{12}}_{rel}$, having the relative tangential velocity $\boldsymbol{v}^{\overline{12}}_{rel}=\left(\dot{\boldsymbol{u}}^1-\dot{\boldsymbol{u}}^2\right)\left[\boldsymbol{\tau}^2_1\;\;\;\boldsymbol{\tau}^2_2\right]$ with $\boldsymbol{\tau}^2_i$ as the tangential vectors on the master surface

$||\boldsymbol{t}_T^{trial}||=\sqrt{t^{trial}_{T_i}m_{ij}t^{trial}_{T_j}}$ with $m_{ij}=\frac{\boldsymbol{\tau}^2_i\cdot\boldsymbol{\tau}^2_j}{|\boldsymbol{\tau}^2_i||\boldsymbol{\tau}^2_j|}$

$\boldsymbol{r}^{int}=\underbrace{\int_{\Gamma^1}\delta g t_N\boldsymbol{n}^2dA}_{\boldsymbol{r}^{int}_{normal}}+\underbrace{\int_{\Gamma}\delta\boldsymbol{\zeta}\boldsymbol{t}_tdA}_{\boldsymbol{r}^{int}_{friction}}$

$=\sum_q^{NQ^1}\left\{\left[\boldsymbol{N}^1\;\;\left(-\boldsymbol{N}^2\right)\right]^Tt_N\left(\xi^q_1\right)\gamma^q\boldsymbol{n}^2|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|\right\} +\sum_q^{NQ^1}\left\{\left[\boldsymbol{N}^1\;\;\left(-\boldsymbol{N}^2\right)\right]^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right]\boldsymbol{t}_T\left(\xi_1^q\right)\gamma^q\boldsymbol{n}^2|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|\right\}$

with $\boldsymbol{N}^i=\left[\begin{array}{ccc|c|ccc}N_1^i&0&0&\cdots&N_n^i&0&0\\0&N_1^i&0&\cdots&0&N_n^i&0\\0&0&N_1^i&\cdots&0&0&N_n^i\end{array}\right]$

### LeftHandSide-Contribution

$\boldsymbol{K}=\left[\begin{array}{cc}\boldsymbol{K}^1_1&\boldsymbol{K}^1_2\\\boldsymbol{K}^2_1&\boldsymbol{K}^2_2\end{array}\right]$ where the indices refer to the variation (top) and the linearization (bottom) regarding the degrees of freedom of the slave (1) and the master (2) facet.

#### Contributions due to a linearization of the virtual work associated with the normal contact force

$\boldsymbol{K}^{1}_{1}=\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^{1}_{2}=-\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|-\sum_q^{NQ}\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)\epsilon_N\left(\frac{\delta\boldsymbol{n}^2}{\delta\boldsymbol{u}^2}\left(\boldsymbol{x}^1-\boldsymbol{x}^2\right)\right)\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$-\sum_q^{NQ^1}\left[\boldsymbol{N}^1\right]^T\frac{\delta\boldsymbol{n}^2}{\delta \boldsymbol{u}^2}t_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^{2}_{1}=-\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^1\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$+\sum_q^{NQ^1}\left[\boldsymbol{N}^1\right]^T\frac{\delta\boldsymbol{n}^2}{\delta \boldsymbol{u}^2}t_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^{2}_{2}=\sum_q^{NQ^1}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)^T\epsilon_N\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|+\sum_q^{NQ}\left(\left[\boldsymbol{N}^2\right]^T\boldsymbol{n}^2\right)\epsilon_N\left(\frac{\delta\boldsymbol{n}^2}{\delta\boldsymbol{u}^2}\left(\boldsymbol{x}^1-\boldsymbol{x}^2\right)\right)\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

#### Contributions due to a linearization of the virtual work associated with the frictional contact force

• $||\boldsymbol{t}_T^{trial}||< \mu t_N$:

$\boldsymbol{K}_2^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \left[\frac{\partial\left(\boldsymbol{t}_T^{trial}\right)}{\partial \boldsymbol{u}^2}\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| +\sum_q^{NQ}\left(\boldsymbol{N^1}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$\boldsymbol{K}_2^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \left[\frac{\partial\left(\boldsymbol{t}_T^{trial}\right)}{\partial \boldsymbol{u}^2}\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| -\sum_q^{NQ}\left(\boldsymbol{N^2}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

• $||\boldsymbol{t}_T^{trial}||\geq \mu t_N$:

$\boldsymbol{K}_1^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^1}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$\boldsymbol{K}_2^1=\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^2}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| +\sum_q^{NQ}\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu t_N\left[\frac{\partial\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}}{\partial \boldsymbol{u}^2}\right] \gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$+\sum_q^{NQ}\left(\boldsymbol{N^1}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$\boldsymbol{K}_1^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^1}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$\boldsymbol{K}_2^2=-\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu \left[\frac{\partial t_N}{\partial \boldsymbol{u}^2}\right]\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2| -\sum_q^{NQ}\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\;\;\;\;\frac{\boldsymbol{\tau}_2^2}{|\boldsymbol{\tau}_2^2|}\right] \mu t_N\left[\frac{\partial\frac{\boldsymbol{t}_T}{||\boldsymbol{t}_T||}}{\partial \boldsymbol{u}^2}\right] \gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

$-\sum_q^{NQ}\left(\boldsymbol{N^2}\right)^T\left[\frac{\partial\left(\frac{\boldsymbol{\tau}_1^2}{|\boldsymbol{\tau}_1^2|}\right)}{\partial\boldsymbol{u}^2}\right]\boldsymbol{t}_T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_1^2|$

### Damp-Matrix Contributions

Only for frictional contact problems.

• $||\boldsymbol{t}_T^{trial}||< \mu t_N$:

$\boldsymbol{K}^1_1=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^1_2=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^2_1=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^2_2=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\epsilon_T\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]^T\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

• $||\boldsymbol{t}_T^{trial}||\geq \mu t_N$:

$\boldsymbol{K}^1_1=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_1}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^1_2=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^1\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_2}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^2_1=-\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_1}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

$\boldsymbol{K}^2_2=\sum_{q}^{NQ^1}\left[\left(\boldsymbol{N}^2\right)^T\left[\frac{\boldsymbol{\tau}^2_1}{|\boldsymbol{\tau}^2_1|}\;\;\;\;\frac{\boldsymbol{\tau}^2_2}{|\boldsymbol{\tau}^2_2|}\right]\right]\mu t_N\frac{\partial\frac{\boldsymbol{t}_T}{|\boldsymbol{t}_T|}}{\partial\dot{\boldsymbol{u}}_2}\gamma^q|\boldsymbol{\tau}_1^1\times\boldsymbol{\tau}_2^1|$

## References

• T.A. Laursen Computational Contact and Impact Mechanics, Springer, 2002