5. Resolution strategy#

The resolution strategy is the same as that used by the continuous method in the classical finite element framework [bib 1]. The only difference is that with X- FEM, no pairing is necessary.

5.1. Resolution algorithm#

5.1.1. Contact-friction law#

With X- FEM, the points facing each other are known a priori and this intrinsic « pairing » does not change during the calculation (hypothesis of small displacements). Thus, it is not necessary to carry out a matching phase as in the classical framework. The geometric loop is also removed since there is no re-updating. For each time step, there are 3 interlocking loops left. The friction threshold loop solves the friction problem by looking for a fixed point on the Tresca friction thresholds. The contact status loop (similar to the active constraints method) makes it possible to determine the space of the effective contact points. At the deepest level, the Newton iteration loop makes it possible to solve the remaining nonlinearity, that due to the projection on the unit ball.

For a step of time:

Initialization of friction thresholds \({\lambda }_{s}\).

Loop over friction thresholdsInitialization of contact statuses \(\chi\)
- Loop over contact statuses
  Newton iterations Calculation of the tangent matrix and the second limb
  
  End of Newton iterations

Contact status update \(\chi\)

End of the active constraints loop

Update of friction thresholds \({\lambda }_{s}\)

End of the loop on friction thresholds

5.1.2. Cohesive law.#

In the implementation of cohesive laws, both mixed and regularized, we do not introduce status fields as we could do for contact: the different regimes are managed by the behavior routine itself, directly in Newton’s method. The only operation to be performed in addition to Newton’s iterations is therefore to update the internal variable.

For a step of time:
- Newton iterations Calculation of the tangent matrix and the second limb
- End of Newton iterations

Updating the internal variable \(\alpha\)

One could legitimately ask why the internal variable is not updated during Newton’s iterations. In fact, as it is a parameter measuring irreversibility, and determined by a maximum over time, it should only be updated at each converged time step. In fact, in the opposite case, if this parameter exceeds its equilibrium value during a Newton’s iteration, Newton’s algorithm will then be unable to reduce it to find the equilibrium value.

5.2. Criteria for stopping the loop on contact statuses#

For the supposed non-contacting points \((\chi =0)\), the condition of non-interpenetration \(({d}_{n}\le 0)\) is checked. The test is as follows: if there is interpenetration \(({d}_{n}>0)\) then the point is assumed to contact \((\chi =1)\) during the next iteration of active constraints. Numerically, the test is written \({d}_{n}>{\text{10}}^{-\text{16}}\).

For the points supposed to contact \((\chi =1)\), it is verified that the value of the contact reaction according to the normal is negative \(({\lambda }_{n}\le 0)\). The test is as follows: if there is detachment \(({\lambda }_{n}>0)\) then the point is assumed to be non-contacting \((\chi =0)\) during the next active stress iteration. Numerically, the test is written \({\lambda }_{n}>-{\text{10}}^{-3}\).

Note:

Contact status is defined independently for each integration point of each facet of contact, stop tests are carried out at each of these points.

Note:

The game at these integration points is calculated thanks to the interpolation of the displacement field on the parent element (3D) while the contact reaction is calculated thanks to the interpolation of the contact ddls on the contact facet (2D). We could also calculate the game by interpolating the displacement at the vertices of the contact facet (the latter being determined with the shape functions of the parent element) .

Note:

The numerical values of the 2 stop tests are difficult to determine. In some configurations, contact status oscillations occur and prevent the convergence of the algorithm. This phenomenon should be identified, and in the case where the values in the 2 cases (contacting and non-contacting) are similar, it could be considered that the convergence of the active constraints loop is reached.

5.3. Criteria for stopping the loop on friction thresholds#

The friction threshold loop is considered to have converged if the solution does not change too much from one iteration to the next. More precisely, let \(v\) be the mixed solution vector (displacement, contact multipliers, friction semi-multipliers), we converged at iteration \(i\) if:

\(\frac{\text{max}\mid {v}^{i}-{v}^{i-1}\mid }{\text{max}\mid {v}^{i-1}\mid }<{10}^{-3}\)

where \(\text{max}(u)\) means the max over all components of the \(u\) vector.

5.4. Writing the formulation during a Newton’s iteration#

Let’s rewrite the weak three-field formulation described in paragraph [§ 3] during a Newton iteration. It is necessary to take into account the loop on friction thresholds, the loop on active stresses. Thus, at this level, the friction thresholds become constants noted \({\lambda }_{s}\), the contact statuses \(\chi ({g}_{n})\) become constants \(\chi\). Moreover, the remaining problem being non-linear (because of the projection onto the unit ball), it is linearized by the Newton-Raphson method.

In the one-dimensional case, Newton’s method is an iterative process that makes it possible to approach the zeros of a continuous and differentiable function. We’re getting back to the \(F(x)=0\) resolution. We build a sequence of points \({x}^{k}\) by doing a Taylor expansion of \(F\) in the vicinity of \({x}^{k}\), which gives the first order:

\(F({x}^{k+1})\approx F({x}^{k})+{F}^{\text{'}}({x}^{k})({x}^{k+1}-{x}^{k})\)

By noting \({\mathrm{\delta x}}^{k}={x}^{k+1}-{x}^{k}\) the increment between two successive iterations, the equation linearized at iteration \(k+1\) is then as follows:

\({F}^{\text{'}}({x}^{k}){\mathrm{\delta x}}^{k}=-F({x}^{k})\)

In the case of the finite element method, \({F}^{\text{'}}({x}^{k})\) is similar to the tangent matrix, which can be calculated at each iteration if necessary, \({\mathrm{\delta x}}^{k}\) is the vector of the increments of the nodal unknowns, and \(F({x}^{k})\) is the second member. Note that \({F}^{\text{'}}({x}^{k})\) and \(F({x}^{k})\) only involve quantities from iteration \(k\), which are therefore known quantities.

The projection on the unit ball is written:

\({P}_{B(\mathrm{0,1})}(x)=\{\begin{array}{}x\text{si}x\in B(\mathrm{0,1})\\ \frac{x}{\parallel x\parallel }\text{sinon}\end{array}\)

The differential of this application, at any point not located on the edge of \(B(\mathrm{0,1})\), is defined by:

\({\partial }_{x}{P}_{B(\mathrm{0,1})}(x)\mathrm{\delta x}=K(x)\mathrm{\delta x}\)

with

\(K(x)=\{\begin{array}{ccc}{I}_{d}& \text{si}x\in B(\mathrm{0,1})& (\text{adhérence})\\ \frac{1}{\parallel x\parallel }({I}_{d}-\frac{x\cdot {x}^{T}}{{\parallel x\parallel }^{2}})& \text{sinon}& (\text{glissement})\end{array}\)

So the \({P}_{B(\mathrm{0,1})}({g}_{\tau })\) differential will be:

\(K({g}_{\tau }){\mathrm{\delta g}}_{\tau }=K(\Lambda +\frac{{\rho }_{\tau }}{\mathrm{\Delta t}}\Delta [[u]{]}_{\tau })(\mathrm{\delta \Lambda }+\frac{{\rho }_{\tau }}{\mathrm{\Delta t}}[[\mathrm{\delta u}]{]}_{\tau })\)

where \({g}_{\tau }\) is the friction semi-multiplier from the previous Newton iteration and \({\mathrm{\delta g}}_{\tau }\) is the increment of unknowns. Knowing the friction semi-multiplier resulting from the previous Newton iteration makes it easy to know whether the point is in the adherent or sliding state.

In the same way, in the case of the cohesive regularized law CZM_EXP_REG, for example, the differential of \({t}_{c}(\mathrm{〚}u\mathrm{〛})\) will be \(\frac{\mathrm{\partial }{t}_{c}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{\cdot }\mathrm{〚}\delta u\mathrm{〛}\) with:

Gold \(\frac{\mathrm{\partial }{t}_{c}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{=}H({\mathrm{〚}u\mathrm{〛}}_{\text{èq}}\mathrm{-}\alpha )\frac{\mathrm{\partial }{\sigma }_{\mathit{lin}}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}+(1\mathrm{-}H({\mathrm{〚}u\mathrm{〛}}_{\text{èq}}\mathrm{-}\alpha ))\frac{\mathrm{\partial }{\sigma }_{\mathit{dis}}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}+\frac{\mathrm{\partial }{\sigma }_{\mathit{pen}}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\)

We reuse the expressions of these three partial derivatives that are given in [R7.02.11], with \(\delta \mathrm{=}\mathrm{-}\mathrm{〚}u\mathrm{〛}\). In the expression for \({\sigma }_{\mathit{dis}}\), you must write \(\alpha \mathrm{=}{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\), which thus becomes a variable to be taken into account in the derivation. In practice, in the code, four cases are distinguished for clarity of reading:

- \({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\mathrm{\ge }\alpha \mathrm{et}\mathrm{〚}{u}_{n}\mathrm{〛}<0\) (non-contacting dissipative). So we have: \(\frac{\partial {t}_{c}}{\partial 〚u〛}=-{\sigma }_{c}\mathrm{exp}(-\frac{{\sigma }_{c}}{{G}_{c}}{〚u〛}_{\mathrm{èq}})(\frac{\text{Id}}{{〚u〛}_{\mathrm{èq}}}-\frac{〚u〛}{{〚u〛}_{\mathrm{èq}}}\otimes \frac{〚u〛}{{〚u〛}_{\mathrm{èq}}}(\frac{{\sigma }_{c}}{{G}_{c}}+\frac{1}{{〚u〛}_{\mathrm{èq}}}))\)
\({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}<\alpha \mathrm{et}\mathrm{〚}{u}_{n}\mathrm{〛}\mathrm{\ge }0\) (rubber band contacting). With \(({\tau }_{1},{\tau }_{2})\) a tangential plane base, we have: \(\frac{\mathrm{\partial }{t}_{c}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{=}\frac{\mathrm{\partial }{\sigma }_{\mathit{lin}}(\mathrm{〚}{u}_{\tau }\mathrm{〛})}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}+\frac{\mathrm{\partial }{\sigma }_{\mathit{pen}}(\mathrm{〚}{u}_{n}\mathrm{〛})}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{=}\mathrm{-}Cn\mathrm{\otimes }n\mathrm{-}\frac{{\sigma }_{c}}{\alpha }\mathrm{exp}(\frac{\mathrm{-}{\sigma }_{c}}{{G}_{c}}\alpha )({\tau }_{1}\mathrm{\otimes }{\tau }_{1}+{\tau }_{2}\mathrm{\otimes }{\tau }_{2})\)
\({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\mathrm{\ge }\alpha \mathrm{et}\mathrm{〚}{u}_{n}\mathrm{〛}\mathrm{\ge }0\) (dissipative contacting). By similar reasoning by replacing \({\sigma }_{\mathit{lin}}\) with \({\sigma }_{\mathit{dis}}\), we get: \(\frac{\partial {t}_{c}}{\partial 〚u〛}=-Cn\otimes n-\mathrm{exp}(\frac{-{\sigma }_{c}}{{G}_{c}}{〚u〛}_{\mathrm{èq}})\left[\frac{{\sigma }_{c}}{{〚u〛}_{\mathrm{èq}}}({\tau }_{1}\otimes {\tau }_{1}+{\tau }_{2}\otimes {\tau }_{2})-(\frac{{\sigma }_{c}^{2}}{{G}_{c}}+\frac{{\sigma }_{c}}{{〚u〛}_{\mathrm{èq}}})\frac{{〚u〛}_{\tau }\otimes {〚u〛}_{\tau }}{{{〚u〛}_{\mathrm{èq}}}^{2}}\right]\)
\({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}<\alpha \mathrm{et}\mathrm{〚}{u}_{n}\mathrm{〛}<0\) (non-contacting elastic). \(\frac{\mathrm{\partial }{t}_{c}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{=}\mathrm{-}\frac{{\sigma }_{c}}{\alpha }\mathrm{exp}(\frac{\mathrm{-}{\sigma }_{c}}{{G}_{c}}\alpha )\text{Id}\)

5.5. Linearization of the problem#

5.5.1. Integral writing with the Augmented Lagrangian method#

The linear system of the three equations at Newton’s iteration \(k+1\) is written as follows (to avoid making it heavier, the references to Newton’s iteration are omitted, because obviously, the unknowns are noted with a \(\delta\) in front, and the test fields are now noted with a star):

Find \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) such as:

\(\forall ({u}^{\text{*}},{\lambda }^{\text{*}},{\Lambda }^{\text{*}})\in {V}_{0}\times H\times H\)

Equation of balance	\(\begin{array}{c}{\mathrm{\int }}_{\Omega }\sigma (\mathit{\delta u})\mathrm{:}\epsilon ({u}^{\text{}})d\Omega \\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}\text{χδλ}\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{nd\Gamma }}_{c}+{\mathrm{\int }}_{{\Gamma }_{c}}{\text{χρ}}_{n}\mathrm{[}\mathrm{[}\mathit{\delta u}\mathrm{]}\mathrm{]}\mathrm{\cdot }n\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{nd\Gamma }}_{c}\\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}{\text{χμλ}}_{s}K({g}_{\tau }){\mathit{\delta g}}_{\tau }\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{\tau }{\mathit{d\Gamma }}_{c}\\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Omega }\sigma (u)\mathrm{:}\epsilon ({u}^{\text{}})d\Omega +{\mathrm{\int }}_{\Omega }f\mathrm{\cdot }{u}^{\text{}}d\Omega +{\mathrm{\int }}_{{\Gamma }_{t}}t\mathrm{\cdot }{u}^{\text{}}{\mathit{d\Gamma }}_{t}\\ +{\mathrm{\int }}_{{\Gamma }_{c}}\chi (\lambda \mathrm{-}{\rho }_{n}\mathrm{[}\mathrm{[}u\mathrm{]}\mathrm{]}\mathrm{\cdot }n)\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{nd\Gamma }}_{c}\\ +{\mathrm{\int }}_{{\Gamma }_{c}}{\text{χμλ}}_{s}{P}_{B(\mathrm{0,1})}({g}_{\tau })\mathrm{\cdot }\mathrm{[}\mathrm{[}{u}^{\text{*}}\mathrm{]}{\mathrm{]}}_{\tau }{\mathit{d\Gamma }}_{c}\end{array}\)
Law of contact	\(\begin{array}{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{(1\mathrm{-}\chi )}{{\rho }_{n}}\text{δλ}{\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\chi \mathrm{[}\mathrm{[}\mathit{\delta u}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{n\lambda }}^{\text{}}{\mathit{d\Gamma }}_{c}\\ \mathrm{=}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{(1\mathrm{-}\chi )}{{\rho }_{n}}\lambda {\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}+{\mathrm{\int }}_{\mathit{\Gamma c}}\chi \mathrm{[}\mathrm{[}u\mathrm{]}\mathrm{]}\mathrm{\cdot }n{\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\end{array}\)
Law of friction	\(\begin{array}{c}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}\mathit{\Delta t}}{{\rho }_{\tau }}\left[\mathit{\delta \Lambda }\mathrm{-}K({g}_{\tau }){\mathit{\delta g}}_{\tau }\right]{\Lambda }^{\text{}}{\mathit{d\Gamma }}_{c}+{\mathrm{\int }}_{\mathit{\Gamma c}}(1\mathrm{-}\chi ){\text{δΛΛ}}^{\text{}}{\mathit{d\Gamma }}_{c}\\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}\mathit{\Delta t}}{{\rho }_{\tau }}\left[\Lambda \mathrm{-}{P}_{B(\mathrm{0,1})}({g}_{\tau })\right]{\Lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}(1\mathrm{-}\chi )\text{Λ}{\text{Λ}}^{\text{}}{\mathit{d\Gamma }}_{c}\end{array}\)

5.5.2. Integral writing with the penalized method#

Find \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) such as:

\(\forall ({u}^{\text{*}},{\lambda }^{\text{*}},{\Lambda }^{\text{*}})\in {V}_{0}\times H\times H\)

Equation of balance	\(\begin{array}{c}{\mathrm{\int }}_{\Omega }\sigma (\mathit{\delta u})\mathrm{:}\epsilon ({u}^{\text{}})d\Omega \\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}\text{χ}\delta \lambda \mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{nd\Gamma }}_{c}\\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}{\text{χμλ}}_{s}{\kappa }_{\tau }K({\kappa }_{\tau }{\nu }_{\tau })\delta {\nu }_{\tau }\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{\tau }{\mathit{d\Gamma }}_{c}\\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Omega }\sigma (u)\mathrm{:}\epsilon ({u}^{\text{}})d\Omega +{\mathrm{\int }}_{\Omega }f\mathrm{\cdot }{u}^{\text{}}d\Omega +{\mathrm{\int }}_{{\Gamma }_{t}}t\mathrm{\cdot }{u}^{\text{}}{\mathit{d\Gamma }}_{t}\\ +{\mathrm{\int }}_{{\Gamma }_{c}}\chi \lambda \mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}\mathrm{]}\mathrm{\cdot }{\mathit{nd\Gamma }}_{c}\\ +{\mathrm{\int }}_{{\Gamma }_{c}}{\text{χμλ}}_{s}{P}_{B(\mathrm{0,1})}({\kappa }_{\tau }{\nu }_{\tau })\mathrm{\cdot }\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{\tau }{\mathit{d\Gamma }}_{c}\end{array}\)
Law of contact	\(\begin{array}{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{(1\mathrm{-}\chi )}{{\kappa }_{n}}\text{δλ}{\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\chi (\frac{\text{δλ}}{{\kappa }_{n}}+\mathrm{[}\mathrm{[}\mathit{\delta u}\mathrm{]}\mathrm{]}\mathrm{\cdot }n){\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\\ \mathrm{=}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{(1\mathrm{-}\chi )}{{\kappa }_{n}}\lambda {\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}+{\mathrm{\int }}_{\mathit{\Gamma c}}\chi (\frac{\text{λ}}{{\kappa }_{n}}+\mathrm{[}\mathrm{[}u\mathrm{]}\mathrm{]}\mathrm{\cdot }n){\lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\end{array}\)
Law of friction	\(\begin{array}{c}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}}{{\kappa }_{\tau }}\left[\mathit{\delta \Lambda }\mathrm{-}\textcolor[rgb]{1,0,0}{K({\kappa }_{\tau }{\nu }_{\tau }){\kappa }_{\tau }\delta {\nu }_{\tau }}\right]{\Lambda }^{\text{}}{\mathit{d\Gamma }}_{c}+{\mathrm{\int }}_{\mathit{\Gamma c}}(1\mathrm{-}\chi ){\text{δΛΛ}}^{\text{}}{\mathit{d\Gamma }}_{c}\\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}}{{\kappa }_{\tau }}\left[\Lambda \mathrm{-}{P}_{B(\mathrm{0,1})}({\kappa }_{\tau }{\nu }_{\tau })\right]{\Lambda }^{\text{}}{\mathit{d\Gamma }}_{c}\mathrm{-}{\mathrm{\int }}_{\mathit{\Gamma c}}(1\mathrm{-}\chi )\text{Λ}{\text{Λ}}^{\text{}}{\mathit{d\Gamma }}_{c}\end{array}\)

Here, we have chosen to solve the contact and friction problems implicitly (the contact and friction semi-multipliers are expressed as a function of the jump in displacement of the current Newton iteration). This choice makes the stiffness matrix non-symmetric, block \({B}_{r}\) (see [§ 5.6.1]) being non-zero because of the term in red when writing the law of friction while block \({{B}_{r}}^{T}\) is zero.

5.5.3. Integral writing for a formulation with regularized cohesive law#

Find \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) such as:

\(\forall ({u}^{\text{*}},{\lambda }^{\text{*}},{\Lambda }^{\text{*}})\in {V}_{0}\times H\times H\)

Equation of balance	\(\begin{array}{c}{\mathrm{\int }}_{\Omega }\sigma (\mathit{\delta u})\mathrm{:}\epsilon ({u}^{\text{}})d\Omega \\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}\left[\frac{\mathrm{\partial }{t}_{c,n}}{\mathrm{\partial }\mathrm{[}\mathrm{[}u\mathrm{]}{\mathrm{]}}_{n}}\mathrm{[}\mathrm{[}\delta u\mathrm{]}{\mathrm{]}}_{n}+\frac{\mathrm{\partial }{t}_{c,n}}{\mathrm{\partial }\mathrm{[}\mathrm{[}u\mathrm{]}{\mathrm{]}}_{\tau }}\mathrm{[}\mathrm{[}\delta u\mathrm{]}{\mathrm{]}}_{\tau }\right]\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{n}{\mathit{d\Gamma }}_{c}\\ \mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}\left[\frac{\mathrm{\partial }{t}_{c,\tau }}{\mathrm{\partial }\mathrm{[}\mathrm{[}u\mathrm{]}{\mathrm{]}}_{n}}\mathrm{[}\mathrm{[}\delta u\mathrm{]}{\mathrm{]}}_{n}+\frac{\mathrm{\partial }{t}_{c,\tau }}{\mathrm{\partial }\mathrm{[}\mathrm{[}u\mathrm{]}{\mathrm{]}}_{\tau }}\mathrm{[}\mathrm{[}\delta u\mathrm{]}{\mathrm{]}}_{\tau }\right]\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{\tau }{\mathit{d\Gamma }}_{c}\\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Omega }\sigma (u)\mathrm{:}\epsilon ({u}^{\text{}})d\Omega +{\mathrm{\int }}_{\Omega }f\mathrm{\cdot }{u}^{\text{}}d\Omega +{\mathrm{\int }}_{{\Gamma }_{t}}t\mathrm{\cdot }{u}^{\text{}}{\mathit{d\Gamma }}_{t}\\ +{\mathrm{\int }}_{{\Gamma }_{c}}({t}_{c,n}\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{n}+{t}_{c,\tau }\mathrm{[}\mathrm{[}{u}^{\text{}}\mathrm{]}{\mathrm{]}}_{\tau }){\mathit{d\Gamma }}_{c}\end{array}\)
Interface: normal part	\({\mathrm{\int }}_{\mathit{\Gamma c}}{\lambda }^{\text{*}}(\lambda +\text{δλ}\mathrm{-}{t}_{c,n}\mathrm{\cdot }n){\mathit{d\Gamma }}_{c}\mathrm{=}0\)
Interface: tangential part	\({\mathrm{\int }}_{\mathit{\Gamma c}}{\Lambda }^{\text{*}}(\Lambda +\delta \Lambda \mathrm{-}{t}_{c,\tau }){\mathit{d\Gamma }}_{c}\mathrm{=}0\)

5.5.4. Integral writing for a formulation with mixed cohesive law#

Find \((\mathit{\delta u},\mathit{\delta \lambda })\mathrm{\in }{V}_{0}\mathrm{\times }H\) such as:

\(\mathrm{\forall }({u}^{\text{*}},{\lambda }^{\text{*}})\mathrm{\in }{V}_{0}\mathrm{\times }H\)

Equation of balance

\(\begin{array}{c}{\mathrm{\int }}_{\Omega }\sigma (\delta u)\mathrm{:}\epsilon ({u}^{\text{*}})d\Omega \mathrm{-}{\mathrm{\int }}_{\Gamma }(\mathit{Id}\mathrm{-}r\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p})\mathrm{\cdot }\delta \lambda \mathrm{\cdot }\mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \\ +{\mathrm{\int }}_{\Gamma }r(\mathit{Id}\mathrm{-}r\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p})\mathrm{\cdot }\mathrm{〚}\delta u\mathrm{〛}\mathrm{\cdot }\mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \\ \mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Omega }\sigma (u)\mathrm{:}\epsilon ({u}^{\text{*}})d\Omega +{\mathrm{\int }}_{\Omega }f\mathrm{\cdot }{u}^{\text{*}}d\Omega +{\mathrm{\int }}_{{\Gamma }_{t}}t\mathrm{\cdot }{u}^{\text{*}}d\Gamma \\ +{\mathrm{\int }}_{{\Gamma }_{c}}\left[\lambda \mathrm{-}r(\mathrm{〚}u\mathrm{〛}+\delta (p))\right]\mathrm{\cdot }\mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \end{array}\)

Interface law

\(\begin{array}{c}\mathrm{-}{\mathrm{\int }}_{{\Gamma }_{c}}(1\mathrm{-}r\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p})\mathrm{\cdot }\mathrm{〚}u\mathrm{〛}\mathrm{\cdot }{\lambda }^{\text{*}}d\Gamma \mathrm{-}{\mathrm{\int }}_{\Gamma }\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p}\mathrm{\cdot }\delta \lambda \mathrm{\cdot }\lambda \text{*}d\Gamma \\ \mathrm{=}{\mathrm{\int }}_{\Gamma }(\mathrm{〚}u\mathrm{〛}+\delta )\mathrm{\cdot }{\lambda }^{\text{*}}d\Gamma \end{array}\)

5.6. Basic terms of rubbing contact#

5.6.1. Matrix writing of the linearized problem#

Using the notations from [bib 1], and considering the unified script adopted for the regularized friction and cohesive contact laws, the linearized system as solved in Newton’s \(k+1\) iteration can be put into matrix form:

Equation of balance	\(\begin{array}{}\left\{{u}^{\text{}}\right\}\left[{K}_{\text{méca}}\right](\mathrm{\delta u})\\ +\left\{{u}^{\text{}}\right\}{\left[A\right]}^{T}(\text{δλ})+\left\{{u}^{\text{}}\right\}\left[{A}_{u}\right](\mathrm{\delta u})\\ +\left\{{u}^{\text{}}\right\}{\left[{B}_{r}\right]}^{T}(\mathrm{\delta \Lambda })+\left\{{u}^{\text{}}\right\}\left[{B}_{u}\right](\mathrm{\delta u})\\ +\left\{{u}^{\text{}}\right\}\left[{D}_{u}\right](\mathrm{\delta u})\\ =\left\{{u}^{\text{}}\right\}({L}_{\text{méca}}^{1})\\ +\left\{{u}^{\text{}}\right\}({L}_{\text{cont}}^{1})\\ +\left\{{u}^{\text{}}\right\}({L}_{\text{frott}}^{1})\\ +\left\{{u}^{\text{}}\right\}({L}_{\text{coh}}^{1})\\ \end{array}\)
Law of contact	\(\begin{array}{}\left\{{\lambda }^{\text{}}\right\}\left[C\right](\text{δλ})+\left\{{\lambda }^{\text{}}\right\}\left[A\right](\mathrm{\delta u})\\ =\left\{{\lambda }^{\text{}}\right\}({L}_{\text{cont}}^{2})+\left\{{\lambda }^{\text{}}\right\}({L}_{\text{coh}}^{2})\\ \end{array}\)
Law of friction	\(\begin{array}{}\left\{{\Lambda }^{\text{}}\right\}\left[{F}_{r}\right](\mathrm{\delta \Lambda })+\left\{{\Lambda }^{\text{}}\right\}\left[{B}_{r}\right](\mathrm{\delta u})\\ =\left\{{\Lambda }^{\text{}}\right\}({L}_{\text{frott}}^{3})+\left\{{\Lambda }^{\text{}}\right\}({L}_{\text{coh}}^{3})\end{array}\)

where the column vectors are denoted \((x)\) and the row vectors \(\left\{x\right\}={(x)}^{T}\) .This system can be put in the following matrix form:

\(\left[\begin{array}{c}{{\rm K}}_{\text{méca}}+{{\rm A}}_{u}+{{\rm B}}_{u}+{D}_{u}\\ {\rm A}\\ {{\rm B}}_{r}\end{array}\begin{array}{c}{{\rm A}}^{T}\\ C\\ 0\end{array}\begin{array}{c}{{\rm B}}_{r}^{T}\\ 0\\ {F}_{r}\end{array}\right](\begin{array}{c}\mathrm{\delta u}\\ \text{δλ}\\ \mathrm{\delta \Lambda }\end{array})=(\begin{array}{c}{L}_{\text{méca}}^{1}+{L}_{\text{cont}}^{1}+{L}_{\text{frott}}^{1}+{L}_{\text{coh}}^{1}\\ {L}_{\text{cont}}^{2}+{L}_{\text{coh}}^{2}\\ {L}_{\text{frott}}^{3}+{L}_{\text{coh}}^{3}\end{array})\)

The unknown is the increment compared to the previous Newton iteration. The reference to Newton’s iteration number has been deliberately omitted.

Of course, the cohesion terms \({D}_{u}\), \({L}_{\text{coh}}^{1}\),, \({L}_{\text{coh}}^{2}\), and \({L}_{\text{coh}}^{3}\) only appear in the formulation for regularized cohesive law. If this is the case, all the other terms are null, except for the terms \({K}_{\text{méca}}\) and \({L}_{\text{méca}}^{1}\) of course, but also with the exception of \(C\), \({F}_{r}\), \({L}_{\text{cont}}^{2}\) and \({L}_{\text{frott}}^{3}\) which are used for post-processing, and which are deduced from the detailed expressions to be followed by considering \(\chi \mathrm{=}0\) and \({\rho }_{n}\mathrm{=}1\).

\({K}_{\text{méca}}\) is the mechanical stiffness matrix defined in paragraph [§3.2] of [R7.02.12].

\({A}_{u}\) is the matrix of increased stiffness due to contact.

\({B}_{u}\) is the stiffness matrix increased due to friction.

\({D}_{u}\) is the stiffness matrix due to cohesive forces.

\(A\) is the matrix linking the terms of displacement to those of contact (contact law matrix).

\({B}_{r}\) is the matrix linking the terms of displacement to those of friction (matrix of laws of friction). This matrix is noted \(B\) in bib 1, but in order not to confuse it with the matrix of derivatives of form functions, we’ll write it \({B}_{r}\).

\(C\) is the matrix for determining contact pressures in the non-contacting case.

\({F}_{r}\) is the matrix for determining the friction multipliers in the case of non-friction contact.

\({L}_{\text{méca}}^{1}\) is the second member representing internal forces and load increments.

\({L}_{\text{cont}}^{1}\) and \({L}_{\text{cont}}^{2}\) are the second members due to contact.

\({L}_{\text{frott}}^{1}\) and \({L}_{\text{frott}}^{3}\) are the second members due to friction.

\({L}_{\text{coh}}^{1}\), \({L}_{\text{coh}}^{2}\) and \({L}_{\text{coh}}^{3}\) are the second members due to cohesive forces.

Note:

We remind you that the system solved by Code_Aster is not of the type \(\mathrm{[}K\mathrm{]}\mathrm{[}U\mathrm{]}\mathrm{=}\mathrm{[}F\mathrm{]}\) but of the type \(\mathrm{[}K\mathrm{]}\mathrm{[}U\mathrm{]}+\mathrm{[}F\mathrm{]}\mathrm{=}0\) . There is therefore a minus sign between the second members given in this document and those coded in the Fortran files.

5.6.2. Expression of elementary contact matrices#

5.6.2.1. Augmented Lagrangian method#

Taking into account the discretizations of the fields mentioned in paragraphs [§3.2] of [R7.02.12] and [§ 4.1] of this document, the continuous matrix system above is replaced by a discrete system. More precisely, the matrix \(A\) has the following form:

\(\begin{array}{ccccccccccc}& & \phantom{\text{[}}& \delta {a}_{j}& \delta {b}_{j}& \delta {c}_{j}^{1}& \delta {c}_{j}^{2}& \delta {c}_{j}^{3}& \delta {c}_{j}^{4}& \phantom{\text{]}}& \\ {[A]}_{\mathrm{ij}}& \text{=}& [:ref:`& 0& x& x& 0& 0& 0& <& 0& x& x& 0& 0& 0& >\)] & {lambda} _ {i} ^ {text {*}}}end {array} `

Indeed, this is due to the fact that the contact terms cancel out for ddls whose associated form function is continuous. In fact,

\(\begin{array}{c}{\left\{{\lambda }^{\text{*}}\right\}}_{i}{\left[A\right]}_{\text{ij}}{\left(\mathit{\delta u}\right)}_{j}=-{\int }_{{\Gamma }^{1}}{\text{χψ}}_{i}{\lambda }_{i}^{\text{*}}{\mathrm{\varphi }}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}{\mathit{\delta b}}_{j}+{F}^{1}{\mathit{\delta c}}_{j}^{1}+{F}^{2}{\mathit{\delta c}}_{j}^{2}+{F}^{3}{\mathit{\delta c}}_{j}^{3}+{F}^{4}{\mathit{\delta c}}_{j}^{4}\right)\cdot n\mathit{d\Gamma }\\ +{\int }_{{\Gamma }^{2}}{\text{χψ}}_{i}{\lambda }_{i}^{\text{*}}{\mathrm{\varphi }}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}{\mathit{\delta b}}_{j}+{F}^{1}{\mathit{\delta c}}_{j}^{1}+{F}^{2}{\mathit{\delta c}}_{j}^{2}+{F}^{3}{\mathit{\delta c}}_{j}^{3}+{F}^{4}{\mathit{\delta c}}_{j}^{4}\right)\cdot n\mathit{d\Gamma }\\ =-{\int }_{{\Gamma }^{1}}{\text{χψ}}_{i}{\lambda }_{i}^{\text{*}}{\mathrm{\varphi }}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}({x}^{\text{-}}){\mathit{\delta b}}_{j}-\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right)\cdot n\mathit{d\Gamma }+{\int }_{{\Gamma }^{2}}{\text{χψ}}_{i}{\lambda }_{i}^{\text{*}}{\mathrm{\varphi }}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}({x}^{\text{+}}){\mathit{\delta b}}_{j}+\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right)\mathrm{.}n\mathit{d\Gamma }\\ ={\int }_{\Gamma }{\text{χψ}}_{i}{\lambda }_{i}^{\text{*}}{\mathrm{\varphi }}_{j}\left(2{\mathit{\delta b}}_{j}+2\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right)\cdot n\mathit{d\Gamma }\end{array}\)

Here, we can clearly see the product of the triangle’s \({\psi }_{i}\) shape functions with the triangle’s shape functions with the \({\phi }_{j}\) shape of the parent volume element.

Note:

Following this calculation, we note the expression for the jump in displacement as a function of the enriched degrees of freedom X- FEM: \(\begin{array}{c}{[[u]]}_{j}={\left({a}_{j}+{H}_{j}{b}_{j}+{F}^{1}{c}_{j}^{1}+{F}^{2}{c}_{j}^{2}+{F}^{3}{c}_{j}^{3}+{F}^{4}{c}_{j}^{4}\right)}_{1}\\ \phantom{\rule{14em}{0ex}}-{\left({a}_{j}+{H}_{j}{b}_{j}+{F}^{1}{c}_{j}^{1}+{F}^{2}{c}_{j}^{2}+{F}^{3}{c}_{j}^{3}+{F}^{4}{c}_{j}^{4}\right)}_{2}\\ \phantom{\rule{10em}{0ex}}=\left({a}_{j}+{H}_{j}({x}^{\text{-}}){b}_{j}-\sqrt{r}{c}_{j}^{1}\right)-\left({a}_{j}+{H}_{j}({x}^{\text{+}}){b}_{j}+\sqrt{r}{c}_{j}^{1}\right)\\ \phantom{\rule{10em}{0ex}}=-\left(2{b}_{j}+2\sqrt{r}{c}_{j}^{1}\right)\end{array}\)

Indeed, by construction (cf. [R7.02.12]), the jump coefficient for domain selection functions is equal to:

\({H}_{j}({x}^{\text{+}})-{H}_{j}({x}^{\text{-}})=2\)

Likewise, the matrix of increased stiffness due to contact has the following form:

\(\begin{array}{cccc}& & \begin{array}{cccccc}\phantom{\text{[}}\delta {a}_{j}& \delta {b}_{j}& \delta {c}_{j}^{1}& \delta {c}_{j}^{2}& \delta {c}_{j}^{3}& \delta {c}_{j}^{4}\phantom{\text{]}}\end{array}& \\ {[{A}_{u}]}_{\mathrm{ij}}& \text{=}& \left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& \times & \times & 0& 0& 0\\ 0& \times & \times & 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& \begin{array}{c}{{a}_{i}}^{\text{*}}\\ {{b}_{i}}^{\text{*}}\\ {{c}_{i}^{1}}^{\text{*}}\\ {{c}_{i}^{2}}^{\text{*}}\\ {{c}_{i}^{3}}^{\text{*}}\\ {{c}_{i}^{4}}^{\text{*}}\end{array}\end{array}\)

\({\left\{{u}^{\text{*}}\right\}}_{i}{\left[{A}_{u}\right]}_{\text{ij}}{(\text{δu})}_{j}={\int }_{\Gamma }{\text{χρ}}_{n}{\phi }_{i}\left\{({\mathrm{2b}}_{{i}^{\text{*}}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\cdot n\right\}{\phi }_{j}(2{\text{δb}}_{j}+2\sqrt{r}{\text{δc}}_{j}^{1})\cdot n\mathrm{d\Gamma }\)

The expression for the \(C\) matrix does not change compared to the classic case without X- FEM:

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{\left[C\right]}_{\text{ij}}{(\text{δλ})}_{j}=-{\int }_{\Gamma }\frac{1}{{\rho }_{n}}(1-\chi ){\psi }_{i}{\lambda }_{i}^{\text{*}}{\psi }_{j}{\text{δλ}}_{j}\mathrm{d\Gamma }\)

5.6.2.2. Penalized method#

The expression for matrix \(C\) is:

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{\left[C\right]}_{\text{ij}}{(\text{δλ})}_{j}=-{\int }_{\Gamma }\frac{1}{{\kappa }_{n}}{\psi }_{i}{\lambda }_{i}^{\text{*}}{\psi }_{j}{\text{δλ}}_{j}\mathrm{d\Gamma }\)

The matrix \({A}_{u}\) is zero.

The expression for matrix \(A\) is:

\(\begin{array}{}{\left\{{\lambda }^{\text{*}}\right\}}_{i}{\left[A\right]}_{\text{ij}}{(\mathrm{\delta u})}_{j}={\int }_{\Gamma }\chi {\text{ψ}}_{i}{\lambda }_{i}^{\text{*}}{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\cdot n\mathrm{d\Gamma }\end{array}\)

5.6.2.3. Formulation for a regularized cohesive law#

A \(C\) matrix necessary for post-processing is used, whose expression is that of augmented Lagrangian where \(\chi \mathrm{=}0\) and \({\rho }_{n}\mathrm{=}1\) have been written. So we get:

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{\left[C\right]}_{\text{ij}}{(\text{δλ})}_{j}\mathrm{=}{\mathrm{\int }}_{\Gamma }{\psi }_{i}{\lambda }_{i}^{\text{*}}{\psi }_{j}{\text{δλ}}_{j}\mathit{d\Gamma }\)

5.6.3. Expression of the second contact members#

5.6.3.1. Augmented Lagrangian method#

These expressions involve quantities from the previous Newton iteration (iteration \(k-1\)). That’s why the reference to the \(k-1\) index was explicitly included:

\({\left\{{u}^{\text{*}}\right\}}_{i}{({L}_{\text{cont}}^{1})}_{i}=-{\int }_{\Gamma }\chi {\phi }_{i}\left\{({\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\cdot n\right\}({\lambda }^{k-1}-{\rho }_{n}{d}_{n}^{k-1})\mathrm{d\Gamma }\)

The expression for vector \({L}_{\text{cont}}^{2}\) does not change compared to the classic case without X- FEM:

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{cont}}^{2})}_{i}={\int }_{\Gamma }{\psi }_{i}{\lambda }_{i}^{\text{*}}(\frac{1-\chi }{{\rho }_{n}}{\lambda }^{k-1}+{\mathrm{\chi d}}_{n}^{k-1})\mathrm{d\Gamma }\)

5.6.3.2. Penalized method#

These expressions involve quantities from the previous Newton iteration (iteration \(k-1\)). That’s why the reference to the \(k-1\) index was explicitly included:

\({\left\{{u}^{\text{*}}\right\}}_{i}{({L}_{\text{cont}}^{1})}_{i}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }\chi {\phi }_{i}\left\{({\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\mathrm{\cdot }n\right\}{\lambda }^{k\mathrm{-}1}\mathit{d\Gamma }\)

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{cont}}^{2})}_{i}={\int }_{\Gamma }\frac{1}{{\kappa }_{n}}{\psi }_{i}{\lambda }_{i}^{\text{*}}{\lambda }^{k-1}\mathrm{d\Gamma }+{\int }_{\Gamma }\chi {\psi }_{i}{\lambda }_{i}^{\text{*}}{d}_{n}^{k-1}\mathrm{d\Gamma }\)

5.6.3.3. Formulation for a regularized cohesive law#

A vector \({L}_{\text{cont}}^{2}\) necessary for post-processing is used, whose expression is that of augmented Lagrangian where \(\chi \mathrm{=}0\) and \({\rho }_{n}\mathrm{=}1\) have been written. So we get:

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{cont}}^{2})}_{i}\mathrm{=}{\mathrm{\int }}_{\Gamma }{\psi }_{i}{\lambda }_{i}^{\text{*}}{\lambda }^{k\mathrm{-}1}\mathit{d\Gamma }\)

5.6.4. Expression of friction matrices#

5.6.4.1. Augmented Lagrangian method#

In order to express quantities in the tangential plane, we use the expression in paragraph [§ 2.3], which is written in matrix form:

\({u}_{\tau }=(\text{Id}-n\otimes n)u=\left[P\right]u\)

In this expression, the matrix \(P\) designates the projection operator on the plane of normal \(n\). The expression for the matrix of this symmetric operator is:

\(\left[P\right]=\left[\begin{array}{ccc}1-{n}_{x}^{2}& -{n}_{x}{n}_{y}& -{n}_{x}{n}_{z}\\ -{n}_{x}{n}_{y}& 1-{n}_{y}^{2}& -{n}_{y}{n}_{z}\\ -{n}_{x}{n}_{z}& -{n}_{y}{n}_{z}& 1-{n}_{z}^{2}\end{array}\right]\)

where \({n}_{x}\), \({n}_{y}\), \({n}_{z}\) are the coordinates of normal \(n\) as defined in [§ 4.5]. With the choice of a constant normal per facet, this matrix, depending only on the normal, has the same value at each Gauss point and can be calculated only once for each facet.

The increased stiffness matrix due to friction is written as follows:

\({\left\{{u}^{\text{*}}\right\}}_{i}{\left[{B}_{u}\right]}_{\text{ij}}{(\mathrm{\delta u})}_{j}=-{\int }_{\Gamma }{\text{χμλ}}_{s}\frac{{\rho }_{\tau }}{\mathrm{\Delta t}}{\phi }_{i}\left\{{\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\}{\left[P\right]}^{T}\left[{K}_{n}\right]{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\left[P\right]\mathrm{d\Gamma }\)

where the matrix \({K}_{n}\) represents the tangent matrix of the projection onto the unit ball of the friction semi-multiplier increased at the previous Newton iteration: \({K}_{n}=K({g}_{\tau })\). It is a known matrix.

The expression for matrix \({B}_{r}\) is:

\({\left\{{\Lambda }^{\text{*}}\right\}}_{i}{\left[{B}_{r}\right]}_{\text{ij}}{(\mathrm{\delta u})}_{j}={\int }_{\Gamma }{\text{χμλ}}_{s}{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}\left[{K}_{n}\right]\left[P\right]{\phi }_{j}({\mathrm{2b}}_{j}+2\sqrt{r}{c}_{j}^{1})\mathrm{d\Gamma }\)

The expression for matrix \({F}_{r}\) is:

\(\begin{array}{}{\left\{{\Lambda }^{\text{*}}\right\}}_{i}{\left[{F}_{r}\right]}_{\text{ij}}{(\mathrm{\delta \Lambda })}_{j}={\int }_{\Gamma }\frac{{\text{χμλ}}_{s}\mathrm{\Delta t}}{{\rho }_{\tau }}{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}\left[{I}_{d}-{K}_{n}\right]{\psi }_{j}({\Lambda }_{j}^{1}{\tau }_{j}^{1}+{\Lambda }_{j}^{2}{\tau }_{j}^{2})\mathrm{d\Gamma }\\ +{\int }_{\Gamma }(1-\chi ){\psi }_{i}\left\{\begin{array}{cc}{\Lambda }_{i}^{1\text{*}}& {\Lambda }_{i}^{2\text{*}}\end{array}\right\}\cdot \left[\begin{array}{cc}{\tau }_{i}^{1}{\tau }_{j}^{1}& {\tau }_{i}^{1}{\tau }_{j}^{2}\\ {\tau }_{i}^{2}{\tau }_{j}^{1}& {\tau }_{i}^{2}{\tau }_{j}^{2}\end{array}\right]\cdot {\psi }_{j}(\begin{array}{c}{\Lambda }_{j}^{1}\\ {\Lambda }_{j}^{2}\end{array})\mathrm{d\Gamma }\end{array}\)

5.6.4.2. Penalized method#

The expression for matrix \({F}_{r}\) is:

\(\begin{array}{}{\left\{{\Lambda }^{\text{*}}\right\}}_{i}{\left[{F}_{r}\right]}_{\text{ij}}{(\mathrm{\delta \Lambda })}_{j}={\int }_{\Gamma }((1-\chi )+\chi \frac{\mu {\lambda }_{s}}{{\kappa }_{\tau }}){\psi }_{i}\left\{\begin{array}{cc}{\Lambda }_{i}^{1\text{*}}& {\Lambda }_{i}^{2\text{*}}\end{array}\right\}\cdot \left[\begin{array}{cc}{\tau }_{i}^{1}{\tau }_{j}^{1}& {\tau }_{i}^{1}{\tau }_{j}^{2}\\ {\tau }_{i}^{2}{\tau }_{j}^{1}& {\tau }_{i}^{2}{\tau }_{j}^{2}\end{array}\right]\cdot {\psi }_{j}(\begin{array}{c}{\Lambda }_{j}^{1}\\ {\Lambda }_{j}^{2}\end{array})\end{array}\mathrm{d\Gamma }\)

The expression for matrix \({B}_{u}\) is:

\({\left\{{u}^{\text{*}}\right\}}_{i}{\left[{B}_{u}\right]}_{\text{ij}}{(\mathrm{\delta u})}_{j}=-{\int }_{\Gamma }\chi {\text{μλ}}_{s}\frac{{\kappa }_{\tau }}{\Delta t}{\phi }_{i}\left\{{\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\}{\left[P\right]}^{T}\left[{K}_{n}\right]{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\left[P\right]\mathrm{d\Gamma }\)

The expression for matrix \({B}_{r}\) is:

In the penalized method, the stiffness matrix is not symmetric. We don’t have \({B}_{r}^{T}={B}_{r}\) but a zero matrix instead of \({B}_{r}^{T}\).

5.6.4.3. Formulation for a regularized cohesive law#

A \({F}_{r}\) matrix necessary for post-processing is used, whose expression is that of augmented Lagrangian where \(\chi \mathrm{=}0\) was written. So we get:

\({\left\{{\Lambda }^{\text{*}}\right\}}_{i}{\left[{F}_{r}\right]}_{\text{ij}}{(\mathit{\delta \Lambda })}_{j}\mathrm{=}{\mathrm{\int }}_{\Gamma }{\psi }_{i}\left\{\begin{array}{cc}{\Lambda }_{i}^{1\text{*}}& {\Lambda }_{i}^{2\text{*}}\end{array}\right\}\mathrm{\cdot }\left[\begin{array}{cc}{\tau }_{i}^{1}{\tau }_{j}^{1}& {\tau }_{i}^{1}{\tau }_{j}^{2}\\ {\tau }_{i}^{2}{\tau }_{j}^{1}& {\tau }_{i}^{2}{\tau }_{j}^{2}\end{array}\right]\mathrm{\cdot }{\psi }_{j}(\begin{array}{c}{\Lambda }_{j}^{1}\\ {\Lambda }_{j}^{2}\end{array})\mathit{d\Gamma }\)

5.6.5. Expression of the second members of friction#

5.6.5.1. Augmented Lagrangian method#

The second members of friction have the following expressions:

\(\begin{array}{}{\left\{{u}^{\text{*}}\right\}}_{i}{({L}_{\text{frott}}^{1})}_{i}=-{\int }_{\Gamma }{\text{χμλ}}_{s}{\phi }_{i}\left\{{\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{1\text{*}}\right\}{\left[P\right]}^{T}{P}_{B(\mathrm{0,1})}({g}_{\tau }^{k-1})\mathrm{d\Gamma }\\ {\left\{{\Lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{frott}}^{3})}_{i}=-{\int }_{\Gamma }\frac{{\text{χμλ}}_{s}\mathrm{\Delta t}}{{\rho }_{\tau }}{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}({\Lambda }^{k-1}-{P}_{B(\mathrm{0,1})}({g}_{\tau }^{k-1}))\mathrm{d\Gamma }\\ -{\int }_{\Gamma }(1-\chi ){\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}{\Lambda }^{k-1}\mathrm{d\Gamma }\end{array}\)

where \(k\mathrm{-}1\) represents the index of the previous Newton iteration.

5.6.5.2. Penalized method#

The second members of friction have the following expressions:

\(\begin{array}{}{\left\{{\Lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{frott}}^{3})}_{i}=-{\int }_{\Gamma }{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}((1-\chi ){\Lambda }^{k-1}+\chi \frac{\mu \mathrm{.}{\lambda }_{k}}{{\kappa }_{\tau }}({\Lambda }^{k-1}-{P}_{B(\mathrm{0,1})}({\kappa }_{\tau }{\nu }_{\tau }^{k-1})))\mathrm{d\Gamma }\end{array}\)

\({\left\{{u}^{\text{*}}\right\}}_{i}{({L}_{\text{frott}}^{1})}_{i}=-{\int }_{\Gamma }\chi {\text{μλ}}_{s}{\phi }_{i}\left\{{\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\}\left[{P}_{T}\right]{P}_{B(\mathrm{0,1})}({\kappa }_{\tau }{\nu }_{\tau }^{k-1})\mathrm{d\Gamma }\)

where \(k\mathrm{-}1\) represents the index of the previous Newton iteration.

5.6.5.3. Formulation for cohesive law#

A vector \({L}_{\text{frott}}^{3}\) necessary for post-processing is used, whose expression is that of the augmented Lagrangian where we wrote \(\chi \mathrm{=}0\). So we get:

\(\begin{array}{c}{\left\{{\Lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{frott}}^{3})}_{i}\mathrm{=}{\mathrm{\int }}_{\Gamma }{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}{\Lambda }^{k\mathrm{-}1}\mathit{d\Gamma }\end{array}\)

5.6.6. Expression of cohesion matrices#

Let two directions of the fixed base \(I\) and \(J\) have unit vectors \({e}_{I}\) and \({e}_{J}\). Let’s introduce the tangent matrix of the cohesive law into the fixed base \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\) of coefficients \({\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}}_{IJ}\mathrm{=}{e}_{I}\mathrm{\cdot }\frac{\mathrm{\partial }\delta }{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\mathrm{\cdot }{e}_{J}\). With the expression for \(\frac{\mathrm{\partial }{t}_{c}}{\mathrm{\partial }\mathrm{〚}u\mathrm{〛}}\) given in [§ 5.4], we have the tangent matrix of the cohesive law \(\mathrm{[}{K}_{\mathit{loc}}\mathrm{]}\) in the local base (see doc [R7.02.11]). We then get \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\) by \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\mathrm{=}{\mathrm{[}Q\mathrm{]}}^{T}\mathrm{[}{K}_{\mathit{loc}}\mathrm{]}\mathrm{[}Q\mathrm{]}\), where \(\mathrm{[}Q\mathrm{]}\) is an orthonormal transition matrix defined by:

\(\mathrm{[}Q\mathrm{]}\mathrm{=}\left[\begin{array}{ccc}{n}_{x}& {n}_{y}& {n}_{z}\\ {\tau }_{x}^{1}& {\tau }_{y}^{1}& {\tau }_{z}^{1}\\ {\tau }_{x}^{2}& {\tau }_{y}^{2}& {\tau }_{z}^{2}\end{array}\right]\)

Let \(i\) and \(j\) be two enriched nodes. Let \(\Gamma\) be the intersection of the supports of \(i\) and \(j\). The matrix \(\mathrm{[}{D}_{u}\mathrm{]}\) is then given by:

\({\mathrm{\{}{u}^{\text{*}}\mathrm{\}}}_{i}{\mathrm{[}{D}_{u}\mathrm{]}}_{\mathit{ij}}{\mathrm{\{}\delta u\mathrm{\}}}_{j}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }2{b}_{i}^{\text{*}}{\phi }_{i}\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}2{b}_{j}{\phi }_{j}d\Gamma\)

5.6.7. Expression of the second members of cohesion#

The second members of cohesion have the following expressions:

\({\mathrm{\{}{u}_{i}\mathrm{\}}}^{\text{*}}{({L}_{\mathit{coh}}^{1})}_{i}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }2{b}_{i}^{\text{*}}{\phi }_{i}{t}_{c}^{k\mathrm{-}1}d\Gamma\)

\({\left\{{\lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{coh}}^{2})}_{i}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }{\psi }_{i}{\lambda }_{i}^{\text{*}}({t}_{c,n}^{k\mathrm{-}1}\mathrm{\cdot }n)\mathit{d\Gamma }\)

\(\begin{array}{c}{\left\{{\Lambda }^{\text{*}}\right\}}_{i}{({L}_{\text{coh}}^{3})}_{i}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }{\psi }_{i}\left\{{\Lambda }_{i}^{1\text{*}}{\tau }_{i}^{1}+{\Lambda }_{i}^{2\text{*}}{\tau }_{i}^{2}\right\}{t}_{c,\tau }^{k\mathrm{-}1}\mathit{d\Gamma }\end{array}\)

where \(k\mathrm{-}1\) represents the index of the previous Newton iteration.

To express \({t}_{c}\) in the global database, you can use \(\mathrm{[}Q\mathrm{]}\).

5.6.8. Matrix writing of the problem with mixed cohesive law#

The matrix system as solved in Newton’s iteration \(k+1\) can be put in the following matrix form:

\(\left[\begin{array}{c}{{\rm K}}_{\text{méca}}+{{\rm A}}_{u}\\ {\rm A}\end{array}\begin{array}{c}{{\rm A}}^{T}\\ C\end{array}\right](\begin{array}{c}\mathit{\delta u}\\ \text{δλ}\end{array})\mathrm{=}(\begin{array}{c}{L}_{\text{méca}}^{1}+{L}_{\text{coh}}^{1}\\ {L}_{\text{coh}}^{2}\end{array})\)

5.6.8.1. Expression of elementary cohesion matrices:#

Let two directions of the fixed base \(I\) and \(J\) have unit vectors \({e}_{I}\) and \({e}_{J}\). As before, we introduce the tangent matrix of the cohesive law into the fixed base \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\) of coefficients \({\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}}_{\mathit{IJ}}\mathrm{=}{e}_{I}\mathrm{\cdot }\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p}\mathrm{\cdot }{e}_{J}\). Through the law of cohesive behavior, we have the tangent matrix \(\mathrm{[}{K}_{\mathit{loc}}\mathrm{]}\) in the local base (see doc [R7.02.11]). We then get \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\) by \(\mathrm{[}{K}_{\mathit{gl}}\mathrm{]}\mathrm{=}{\mathrm{[}Q\mathrm{]}}^{T}\mathrm{[}{K}_{\mathit{loc}}\mathrm{]}\mathrm{[}Q\mathrm{]}\), where \(\mathrm{[}Q\mathrm{]}\) is an orthonormal transition matrix defined by:

Having introduced these notations, we have:

\({\mathrm{\{}{u}^{\text{*}}\mathrm{\}}}_{i}{\mathrm{[}{A}_{u}\mathrm{]}}_{\mathit{ij}}{\mathrm{\{}\delta u\mathrm{\}}}_{j}\mathrm{=}{\mathrm{\int }}_{\Gamma }2{\phi }_{i}{b}_{i}^{\text{*}}r(\mathrm{[}\mathit{Id}\mathrm{]}\mathrm{-}r\mathrm{[}{K}_{\mathit{gl}}\mathrm{]})2{\phi }_{j}{b}_{j}d\Gamma\)

In addition, we choose to discretize \(\lambda\) on a local basis. The coefficient \((1,J)\) of the matrix \(\mathrm{\{}{\lambda }_{n}^{\text{*}}\mathrm{\}}{\mathrm{[}{A}_{u}\mathrm{]}}_{\mathit{ij}}{\mathrm{\{}\delta u\mathrm{\}}}_{j}\) is then \({\mathrm{\int }}_{\Gamma }2{\phi }_{i}{\lambda }_{i}^{\text{*}}n\mathrm{\cdot }(\mathit{Id}\mathrm{-}r\frac{\mathrm{\partial }\delta }{\mathrm{\partial }p})\mathrm{\cdot }{e}_{j}2{\phi }_{j}{b}_{j}d\Gamma\). For the coefficients \((2,J)\) and \((3,J)\), the formula is the same by replacing \(n\) with \({\tau }_{1}\) and \({\tau }_{2}\), respectively. By exploiting the notations introduced in the preceding paragraph, we deduce:

\({\{{\lambda }^{\text{*}}\}}_{i}{[{A}_{u}]}_{\mathit{ij}}{\{\delta u\}}_{j}={\int }_{\Gamma }{\varphi }_{i}{\lambda }_{i}^{\text{*}}\left([\mathit{Id}]-r[{K}_{\mathit{loc}}]\right)\cdot [Q]2{\varphi }_{j}{b}_{j}d\Gamma\)

As for matrix \(\mathrm{[}C\mathrm{]}\), it is simply written:

\({\mathrm{\{}{\lambda }^{\text{*}}\mathrm{\}}}_{i}{\mathrm{[}C\mathrm{]}}_{\mathit{ij}}{\mathrm{\{}\delta \lambda \mathrm{\}}}_{j}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\Gamma }{\phi }_{i}{\lambda }_{i}^{\text{*}}\mathrm{[}{K}_{\mathit{loc}}\mathrm{]}{\phi }_{j}{\lambda }_{j}d\Gamma\)

5.6.8.2. Expression of the elementary cohesion vectors:#

The expression for the coefficient \(I\) of the \({\mathrm{\{}u\mathrm{\}}}_{i}^{\text{*}}{({L}_{\mathit{coh}}^{1})}_{i}\) vector is \({\mathrm{\int }}_{\Gamma }2{b}_{i}^{\text{*}}{\phi }_{i}(\mathrm{-}\lambda \mathrm{\cdot }{e}_{I}+r(\mathrm{〚}u\mathrm{〛}+\delta )\mathrm{\cdot }{e}_{I})d\Gamma\). With the notations introduced in the previous part, we deduce:

\({\mathrm{\{}{u}_{i}\mathrm{\}}}^{\text{*}}{({L}_{\mathit{coh}}^{1})}_{i}\mathrm{=}{\mathrm{\int }}_{\Gamma }2{b}_{i}^{\text{*}}{\phi }_{I}(\mathrm{-}{\mathrm{[}Q\mathrm{]}}^{T}\mathrm{\cdot }\mathrm{\{}\lambda \mathrm{\}}+r(\mathrm{\{}\mathrm{〚}u\mathrm{〛}\mathrm{\}}+{\mathrm{[}Q\mathrm{]}}^{T}\mathrm{\cdot }\mathrm{\{}\delta \mathrm{\}}))d\Gamma\)

where \(\mathrm{\{}\delta \mathrm{\}}\) and \(\mathrm{\{}\lambda \mathrm{\}}\) are given on a local basis, and \(\mathrm{\{}\mathrm{〚}u\mathrm{〛}\mathrm{\}}\) on a fixed basis.

The coefficient \(1\) of the vector \({\mathrm{\{}\lambda \mathrm{\}}}_{i}^{\text{*}}{({L}_{\mathit{coh}}^{2})}_{i}\) is written as \({\mathrm{\int }}_{\Gamma }{\lambda }_{i}^{\text{*}}(\mathrm{〚}u\mathrm{〛}\mathrm{\cdot }n+\delta \mathrm{\cdot }n)d\Gamma\). For the coefficients \(2\) and \(3\), the formula is the same by replacing \(n\) with \({\tau }_{1}\) and \({\tau }_{2}\), respectively. From this we deduce:

\({\mathrm{\{}\lambda \mathrm{\}}}_{i}^{\text{*}}{({L}_{\mathit{coh}}^{2})}_{i}\mathrm{=}{\mathrm{\int }}_{\Gamma }{\lambda }_{i}^{\text{*}}(\mathrm{[}Q\mathrm{]}\mathrm{\{}\mathrm{〚}u\mathrm{〛}\mathrm{\}}+\mathrm{\{}\delta \mathrm{\}})d\Gamma\)

where \(\mathrm{\{}\delta \mathrm{\}}\) is given on a local basis, and \(\mathrm{\{}\mathrm{〚}u\mathrm{〛}\mathrm{\}}\) on a fixed basis.