# Nonlinear Cauchy problem and identification in contact mechanics: a solving method based on Bregman-gap

Stéphane Andrieux, Thouraya Baranger

► To cite this version:

Stéphane Andrieux, Thouraya Baranger. Nonlinear Cauchy problem and identification in contact mechanics: a solving method based on Bregman-gap. Inverse Problems, 2020, 36 (11), 10.1088/1361-6420/abbc76. hal-03084876

HAL Id: hal-03084876

https://hal.science/hal-03084876

Submitted on 21 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
NONLINEAR CAUCHY PROBLEM AND IDENTIFICATION IN CONTACT
MECHANICS : A SOLVING METHOD BASED ON BREGMAN-GAP

STÉPHANE ANDRIEUX

ONERA, University Paris-Saclay
Chemin de la Hunière, BP 80100, 92123 Palaiseau, France

THOURAYA NOURI BARANGER

Université de Lyon, Université Claude Bernard Lyon 1
LMC2 EA7427, 69622 Villeurbanne, France

**ABSTRACT.** This paper proposes a solution method for identification problems in the context of contact mechanics when overabundant data are available on a part $\Gamma_m$ of the domain boundary while data are missing from another part of this boundary. The first step is then to find a solution to a Cauchy problem. The method used by the authors for solving Cauchy problems consists of expanding the displacement field known on $\Gamma_m$ towards the inside of the solid via the minimization of a function that measures the gap between solutions of two well-posed problems, each one exploiting only one of the superabundant data. The key question is then to build an appropriate gap functional in strongly nonlinear contexts. The proposed approach exploits a generalization of the Bregman divergence, using the thermodynamic potentials as generating functions within the framework of Generalized Standard Materials, but also Implicit Generalized Standard Materials in order to address Coulomb friction. The robustness and efficiency of the proposed method are demonstrated by a numerical bi-dimensional application dealing with a cracked elastic solid with unilateral contact and friction effects on the crack's lips.

**Keywords.** Cauchy problem, Bregman divergence, full-field displacements, Contact Mechanics, identification, Inverse problems, bi-potentials.

CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Cauchy problems in linear elasticity with contact</td><td>3</td></tr><tr><td>2.1. Formulation of direct problems in contact mechanics</td><td>3</td></tr><tr><td>2.2. Formulation of Cauchy or data completion problems</td><td>4</td></tr><tr><td>3. A general variational method for the solution of Cauchy problems</td><td>5</td></tr><tr><td>4. The symmetrized Bregman Gap and gap functionals</td><td>6</td></tr><tr><td>5. Two solution algorithms for the Cauchy problem</td><td>8</td></tr><tr><td>6. Bregman Gaps in contact mechanics</td><td>9</td></tr><tr><td>6.1. Linear elasticity with unilateral contact and Tresca friction law</td><td>9</td></tr><tr><td>6.2. Linear elasticity with unilateral contact and Coulomb friction</td><td>10</td></tr><tr><td>7. Application to some identification problems in contact mechanics</td><td>11</td></tr><tr><td>7.1. Identification of missing boundary conditions</td><td>11</td></tr><tr><td>7.2. Identification of the friction coefficient</td><td>16</td></tr><tr><td>8. Conclusion</td><td>21</td></tr><tr><td>References</td><td>22</td></tr></table>

E-mail addresses: Stephane.Andrieux@onera.fr, Thouraya.Baranger@univ-lyon1.fr.
Date: January 1, 2020.
# 1. INTRODUCTION

Contact mechanics deals with the special situation where quasi-static loaded deformable bodies are in mutual contact. Because of non-penetration and possibly friction conditions, these problems turn out to be very nonlinear. Nevertheless, the nature of the nonlinearity is very different as no dissipation of energy occurs in the first case, whereas the dissipative feature of the second one requires treating it as an evolution problem [1, 2]. Furthermore the Tresca friction law (constant slip bound) can be recast into the convex formalism of Standard Generalized Materials [3], but this is no more the case for the Amontons-Coulomb's friction law or derived laws where slip threshold depends on the normal traction [4, 5]. Adhesion phenomena can be also present and can be modeled for example by using an extra variable like in the formulation of [6, 7], but although being straightforward to include in the approach proposed in the present paper, adhesion is not addressed here.

Contact mechanics is involved in numerous inverse or identification problems, some of them being specific such as the identification of contact zone or friction parameters. We are here interested in problems where the available data consist in full mechanical fields ($U_m, F_m$), i.e. displacement and stress vector fields, known on some parts of the boundary of the solid bodies (not including the possible contact areas). Practical situations include in particular displacement measurements obtained by digital image correlation (DIC) techniques on a stress-free part $\Gamma_m$ of the external boundary of the solids [8, 9]. In this paper, advantage is taken from this large amount of data and exploitation is made of the superabundant mechanical boundary conditions given by the pair ($U_m, F_m = 0$) by solving a Cauchy problem [10] delivering the whole mechanical fields inside the solids. Then identification procedures are designed or identification results arise directly from the knowledge of the mechanical field itself.

One of the key-points of this approach is of course the solution method for the Cauchy problem for the operators encountered in contact mechanics, that is usually the Lamé system if linear elasticity is used for description of the behavior of the bulk material within the bodies, and the Signorini's unilateral boundary condition [11] combined with the friction law on the contact areas. Both of the last conditions are nonlinear. Solving the Cauchy problem attracted the attention of numerous researchers with a wide range of methods, a large amount of them address nevertheless only linear elliptic operators and 2D situations [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 14]. Very few papers can be found for nonlinear elliptic Cauchy problems [13, 23, 24].

The authors developed in a series of papers solution algorithms for the resolution of the Cauchy problem in linear and nonlinear mechanics [25, 26, 27, 28, 29, 30, 31, 32], including nonlinear elasticity and elastoplasticity, stationary and heat equation [33, 34, 35, 36, 37, 38, 39, 40, 41]. The method relies on a gap functional which is minimized in the proposed approach, provided it is positive, and zero when the gap vanishes. These properties of the gap functional can be automatically fulfilled if this functional is a Bregman divergence. The reason lies in the convexity of the function generating the Bregman divergence, and in thermomechanics the thermodynamic potentials and dissipation functions are generally convex functions so they are natural candidates for generating Bregman divergences well suited to the problem in consideration. However, this has to be adapted in the case with friction because it cannot generally be described by a convex dissipation function. We turn then to the implicit Generalized Standard Material framework which enables to recover convexity. Equipped with the adapted solution method for Cauchy problems in contact mechanics, some inverse or identification problems are tackled with.

The paper is organized as follows. Firstly, the Cauchy problem in the context of linear elasticity and contact condition is recalled, the general variational method proposed for solving Cauchy problem is described. It involves the definition of a gap between the solutions of two forward problems. The adaptation of the method for contact mechanics is described in part 3 where specific gap functions are derived; it is suggested to ground the building of the function on the concept of Bregman divergence by introducing the Bregman Gap. A distinction has to be made for unilateral contact conditions with (or without) Tresca like slip threshold, and unilateral contact with Coulomb's friction condition. The part 5 is devoted to various inverse or identification applications. Then a brief final section offers some concluding observations.
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Notation and geometry of the contacting solids.</figcaption></figure>

## 2. CAUCHY PROBLEMS IN LINEAR ELASTICITY WITH CONTACT

### 2.1. Formulation of direct problems in contact mechanics.
Consider two three-dimensional elastic solids bodies $\Omega_\alpha$, $\alpha = 1, 2$ the boundaries of which are separated in three non-overlapping parts:

$$ (1) \qquad \partial\Omega_\alpha = \bar{\Gamma}_D^\alpha \cup \bar{\Gamma}_N^\alpha \cup \bar{\Gamma}_C^\alpha, \quad \Gamma_D^\alpha \cap \Gamma_N^\alpha = \Gamma_N^\alpha \cap \Gamma_C^\alpha = \Gamma_D^\alpha \cap \Gamma_C^\alpha = \emptyset $$

where the subscripts $D, N, C$ stand respectively for Dirichlet boundary, the Neumann boundary and the contact boundary. In the context of small deformation, the common zone of possible contact is denoted by $\Gamma_C (= \Gamma_C^1 = \Gamma_C^2)$, and the unit normal $\mathbf{n}$ on $\Gamma_C$ is chosen such that: $\mathbf{n} = \mathbf{n}_1 = -\mathbf{n}_2$ where $\mathbf{n}_\alpha$ is the external unit normal to the domain $\Omega_\alpha$.

The relative displacement between the two solids on $\Gamma_C$ is defined by $[\mathbf{u}] = \mathbf{u}^2 - \mathbf{u}^1$, where $\mathbf{u}^\alpha$ is the displacement field in the domain $\Omega_\alpha$. The relative displacement is split into the normal (scalar) and tangential part as follows:

$$ (2) \qquad [\mathbf{u}] = u_N \mathbf{n} + \mathbf{u}_T \quad \text{where } u_N = \mathbf{u} \cdot \mathbf{n} $$

so that the non-interpenetration condition reads: $u_N \ge 0$. If the solid occupying the domain $\Omega_2$ is rigid ($\mathbf{u}^2 = \mathbf{0}$), then the usual Signorini condition $u_N \le 0$ is recovered. Similarly, the normal stress vector $\sigma \cdot \mathbf{n} = \sigma^1 \cdot \mathbf{n} = \sigma^2 \cdot \mathbf{n}$ is decomposed into the normal force density $R_N$ and the vector tangential force density $R_T$.

$$ (3) \qquad \sigma \cdot \mathbf{n} = R_N \mathbf{n} + \mathbf{R}_T \quad \text{with } R_N = \sigma \cdot \mathbf{n} \cdot \mathbf{n} $$

Equipped with these definitions, it is possible to set the boundary conditions on the contact surface, namely:

(1) Signorini's condition or unilateral contact without friction

$$ (4) \qquad u_N \ge 0, \quad R_N \le 0, \quad u_N R_N = 0, \quad \mathbf{R}_T = \mathbf{0} $$

(2) Unilateral contact with dry friction

$$ (5) \qquad \left\{ \begin{array}{l} u_N \ge 0, \quad R_N \le 0, \quad u_N R_N = 0 \\ \|\mathbf{R}_T\| \le g(R_N), \quad \dot{\mathbf{u}}_T = \lambda \mathbf{R}_T, \quad \text{with} \\ \qquad \left| \begin{array}{ll} \lambda = 0 & \text{if } \|\mathbf{R}_T\| < g(R_N) \\ \lambda > 0 & \text{if } \|\mathbf{R}_T\| = g(R_N) \end{array} \right. \end{array} \right. $$

with the slip threshold $g(R_N) = \tau$ for Tresca's friction law, and $g(R_N) = \rho_0|R_N|$ for Coulomb's law. The dissipation has the following expression:

$$ D = R_N \dot{u}_N + \mathbf{R}_T \cdot \dot{\mathbf{u}}_T \equiv \mathbf{R}_T \cdot \dot{\mathbf{u}}_T $$
The direct simulation of friction and unilateral contact mechanics necessitates a time discretization. For preserving convexity and achieving accuracy on the geometry of the contact, sliding and adhesion zones, the backward Euler method (or implicit Euler method) should be used. The equation system is then the Lamé system completed with the boundary conditions on the contact zone $\Gamma_C$ given either by (4) or (5) with their implicit incremental form:

$$ (6) \qquad \left\{ \begin{array}{ll} \operatorname{div}(\sigma^{\alpha} + \Delta\sigma^{\alpha}) = 0 & \text{in } \Omega_{\alpha} \\ \sigma^{\alpha} + \Delta\sigma^{\alpha} = \mathbb{A}^{\alpha} : \epsilon(\mathbf{u}^{\alpha} + \Delta\mathbf{u}^{\alpha}) & \text{in } \Omega_{\alpha} \\ \varepsilon(\mathbf{u}^{\alpha} + \Delta\mathbf{u}^{\alpha}) = \nabla^{s}(\mathbf{u}^{\alpha} + \Delta\mathbf{u}^{\alpha}) & \text{in } \Omega_{\alpha} \\ \mathbf{u}^{\alpha} + \Delta\mathbf{u}^{\alpha} = U^{\alpha} + \Delta U^{\alpha} & \text{on } \Gamma_{D}^{\alpha} \\ (\sigma^{\alpha} + \Delta\sigma^{\alpha}) \cdot \mathbf{n} = F^{\alpha} + \Delta F^{\alpha} & \text{on } \Gamma_{N}^{\alpha} \end{array} \right. $$

A variational formulation of this equation system has been established [7, 42] and reads:

$$ \left\{ \begin{array}{l} \text{for given } \mathbf{u} = (\mathbf{u}^1, \mathbf{u}^2), \\ \text{find } \Delta \mathbf{u} = (\Delta \mathbf{u}^1, \Delta \mathbf{u}^2) \text{ with } \mathbf{u} + \Delta \mathbf{u} \in K \text{ such that :} \\ a(\mathbf{u}+\Delta\mathbf{u}, \mathbf{v}-\mathbf{u}-\Delta\mathbf{u})+J(\mathbf{u}+\Delta\mathbf{u}, \mathbf{v}-\mathbf{u})-J(\mathbf{u}+\Delta\mathbf{u}, \Delta\mathbf{u}) \geq (\mathbf{F}+\Delta\mathbf{F}, \mathbf{v}-\mathbf{u}-\Delta\mathbf{u}) \\ \forall \mathbf{v} \in K \end{array} \right. $$

with the following definition of the space of admissible displacement fields $V$, and the convex cone $K$ of displacement fields fulfilling additionally the non-interpenetration condition:

$$ (7) \qquad \begin{aligned} V &= \left\{ (\mathbf{u}^1, \mathbf{u}^2) \in [H^1(\Omega_1)]^3 \times [H^1(\Omega_2)]^3, \quad \mathbf{u}^\alpha = U^\alpha \text{ on } \Gamma_D^\alpha, \alpha = 1, 2 \right\} \\ K &= \left\{ (\mathbf{u}^1, \mathbf{u}^2) \in V, \quad u_N \ge 0, \text{ on } \Gamma_C \right\} \end{aligned} $$

and the following definition of the bi-linear form $a(\bullet, \bullet)$ and the convex function $J$:

$$ (8) \qquad \begin{aligned} a(\mathbf{u}, \mathbf{v}) &= \sum_{\alpha=1,2} \int_{\Omega_{\alpha}} \mathbb{A}^{\alpha} : \epsilon(\mathbf{u}^{\alpha}) : \epsilon(\mathbf{v}^{\alpha}) d\Omega \\ J(\mathbf{u}) &= \int_{\Gamma_C} g(R_N \mathbf{P}(\mathbf{u}^1)) \|\mathbf{R}_T \mathbf{P}(\mathbf{u}^1)\| d\Gamma \end{aligned} $$

Here $\mathbf{P}$ is the projection on the space $U_s^1$ of elastic statically admissible fields in the domain $\Omega_1$:

$$ (9) \qquad U_s^1 = \left\{ \mathbf{v} \in [H^1(\Omega_1)]^n, \int_{\Omega_1^1} \mathbb{A}^1 : \epsilon(\mathbf{v}) : \epsilon(\mathbf{w}) d\Omega = \int_{\Gamma_N^1} F^1 \cdot \mathbf{w} d\Gamma, \forall \mathbf{w} \in [H^1(\Omega_1)]^n, \mathbf{w} = 0 \text{ on } \Gamma_D^1 \cup \Gamma_C^1 \right\} $$

With the usual properties of the Hooke tensors $\mathbb{A}$ (continuity and coerciveness), the variational problem (7) admits a unique solution, provided the friction coefficient $\rho_0$ is small enough as a supplementary condition in the case of Coulomb's friction law [42].

## 2.2. Formulation of Cauchy or data completion problems.
In order to establish the Cauchy or data completion problem, which exploits superabundant boundary data on a part on the boundary of the contacting solids, another partition of the boundaries $\partial\Omega_i$ has to be defined, namely by introducing:

* a part denoted $\Gamma_m$ where the superabundant data pair $(\mathbf{U}_m, \mathbf{F}_m)$ is available;

* a part denoted by $\Gamma_u$ where the data pair $(\mathbf{U}_u, \mathbf{F}_u)$ is unknown;

* a last part denoted by $\Gamma_b$ where usual combination of Dirichlet and/or Neumann conditions are prescribed (this part can possibly include contact and friction boundary conditions).

Some applications necessitate to redefine the contacting bodies and to use interior surfaces as parts of the boundary $\Gamma_u$. The Cauchy problem and the data completion problem can be defined as follows:

**Cauchy problem.** To find the displacement field $\mathbf{u}$ satisfying the equilibrium condition and the constitutive equation within $\Omega_1 \cup \Omega_2$, the boundary condition on $\Gamma_u$ and $\Gamma_b$ and meeting the superabundant data on $\Gamma_m : \mathbf{u} = \mathbf{U}_m, \sigma(\mathbf{u}) \cdot \mathbf{n} = \mathbf{F}_m$

**Data completion problem.** To find the displacement field $\mathbf{U}_u$ on $\Gamma_u$ such that there exists a displacement field $\mathbf{u}$ satisfying the equilibrium condition and the constitutive equation within
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Notation and geometry of the contacting solids in the case of the Cauchy problem.</figcaption></figure>

$\Omega_1 \cup \Omega_2$, the boundary condition on $\Gamma_u$ and $\Gamma_b$ and meeting the superabundant data on $\Gamma_m$:
$$ \mathbf{u} = \mathbf{U}_m, \sigma(\mathbf{u}) \cdot \mathbf{n} = \mathbf{F}_m $$

The Cauchy problem and the data completion problem are very similar; they differ only by the unknowns of the problem. In the sequel, we shall use the missing boundary Dirichlet data $\mathbf{U}_u$ as an auxiliary variable in order to determine the full displacement field $\mathbf{u}$ (and other mechanical quantities) in the solids, via the solution of the elastic equilibrium problem. The two problems will be then solved equivalently and the terms Cauchy problem or Data completion problem will be used as well.

In the case of quasi-static evolution problems, the Euler implicit scheme is used so the incremental Cauchy or Data completion problems have to be addressed. Nevertheless, the adaptation of the above definition is straightforward.

### 3. A GENERAL VARIATIONAL METHOD FOR THE SOLUTION OF CAUCHY PROBLEMS

As mentioned in the introduction, various methods have been proposed to solve the Cauchy or data completion problem. We describe here a method belonging to the variational approaches, that is based on the minimization of an *ad hoc* functional, the arguments of which are the missing data on the boundary $\Gamma_u$. For that purpose, two auxiliary well-posed direct problems ($\mathcal{P}_1, \mathcal{P}_2$), based on the equations (6), (4) or (5), and parametrized by the Dirichlet data $\Delta\eta$ on $\Gamma_u$, are defined as follows:

$$ (10) \qquad \mathcal{P}_{\alpha}: \left\{ \begin{array}{ll} \operatorname{div}(\sigma_{i}^{\alpha} + \Delta\sigma_{i}^{\alpha}) = 0 & \text{in } \Omega_{\alpha} \\ \sigma_{i}^{\alpha} + \Delta\sigma_{i}^{\alpha} = \mathbb{A}^{\alpha}: \varepsilon(\mathbf{u}_{i}^{\alpha} + \Delta\mathbf{u}_{i}^{\alpha}) & \text{in } \Omega_{\alpha} \\ \epsilon(\mathbf{u}_{i}^{\alpha} + \Delta\mathbf{u}_{i}^{\alpha}) = \nabla^{s}(\mathbf{u}_{i}^{\alpha} + \Delta\mathbf{u}_{i}^{\alpha}) & \text{in } \Omega_{\alpha} \\ (\mathbf{u}_{i}^{\alpha} + \Delta\mathbf{u}_{i}^{\alpha})_{n} \geq 0, (R_{i}^{\alpha} + \Delta R_{i}^{\alpha})_{n} \leq 0, & \text{on } \Gamma_{c} \\ (\mathbf{u}_{i}^{\alpha} + \Delta\mathbf{u}_{i}^{\alpha})_{n} (R_{i}^{\alpha} + \Delta R_{i}^{\alpha})_{n} = 0 & \text{on } \Gamma_{c} \\ (\Delta R_{i}^{\alpha})_{T} = 0 & \text{on } \Gamma_{c} \\ \Delta\mathbf{u}_{i}^{\alpha} = \Delta\eta & \text{on } \Gamma_{u} \end{array} \right. \quad \alpha = 1, 2 \quad i = 1, 2 $$

$$ (11) \qquad \text{with} \quad \left\{ \begin{array}{ll} \Delta\sigma_1^1 \cdot \mathbf{n} = \Delta\mathbf{F}_m & \text{on } \Gamma_m \text{ for } \mathcal{P}_1 \\ \Delta\mathbf{u}_2^1 = \Delta\mathbf{U}_m & \text{on } \Gamma_m \text{ for } \mathcal{P}_2 \end{array} \right. $$

The boundary part $\Gamma_b$ is free of load and only Signorini's condition is used. Each problem uses only one of the superabundant data ($\Delta\mathbf{U}_m, \Delta\mathbf{F}_m$) on the boundary $\Gamma_m$. Now it is clear that if $\Delta\eta$ is such that the two solution fields of problems ($\mathcal{P}_1, \mathcal{P}_2$) are equal, $\Delta\mathbf{u}_1 = \Delta\mathbf{u}_2 = \Delta\mathbf{u}$,
<figure><img src="image_3.png" /><figcaption>FIGURE 3. Geometrical interpretation of the Bregman divergence for $n = 2$.</figcaption></figure>

hence the Cauchy problem is solved with solution $\Delta u$. If a gap functional $E(v, w)$ can be built between the fields $\Delta u_1$ and $\Delta u_2$ with properties:

$$ (1) \quad E(v, w) \geq 0 \quad \forall v, w $$

$$ (2) \quad E(v, w) = 0 \Leftrightarrow v = w $$

then the variational method for solving the Data Completion Problem is simply:

$$ (12) \qquad \Delta \mathbf{U}_u = \operatorname{Arg} \min_{\Delta \eta} f_{\mathrm{obj}}(\Delta \eta) = E(\Delta \mathbf{u}_1(\Delta \eta), \Delta \mathbf{u}_2(\Delta \eta)) $$

and $\Delta \mathbf{u} = \Delta \mathbf{u}_2(\Delta \mathbf{u}_u)$ is the solution of the Cauchy problem.

### 4. THE SYMMETRIZED BREGMAN GAP AND GAP FUNCTIONALS

As just seen, the solution of the Cauchy problem in the proposed approach relies on the design of appropriate gaps between field solutions of equations systems. In the convex realm, advantage can be taken of existence of convex thermodynamic or dissipation potentials. For that purpose, a Bregman Gap, inspired by the notion of Bregman divergence is derived.

**Definition (Bregman Divergence, [43, 44]).** Let $\mathcal{J}$ be a proper convex differentiable function, $\mathcal{J}(\mathbf{e}) : \mathbb{R}^n \to \mathbb{R}$, the Bregman divergence, generated by $\mathcal{J}$, between two points belonging to dom($\mathcal{J}$) is the non-negative scalar :

$$ (13) \qquad D_{\mathcal{J}}(\mathbf{e}_1, \mathbf{e}_2) = \mathcal{J}(\mathbf{e}_1) - \mathcal{J}(\mathbf{e}_2) - \langle \nabla_{\mathcal{J}}(\mathbf{e}_2), \mathbf{e}_1 - \mathbf{e}_2 \rangle $$

Bregman divergences are generally not distances as they are neither symmetric nor satisfying the triangle inequality. Nevertheless, interesting properties have been proved [45] and have been extended to the infinite dimensional context [46], particularly in the context of learning applications, although the Bregman divergence appeared first in optimization algorithms design [43, 44]. Among these properties, the one of additivity is of particular interest in applications because it allows combining two or more generating functions. The geometric interpretation of the Bregman divergence is illustrated on the figure 3 for $n = 2$.

**Property. Additive property.**

$$ D_{\lambda\mathcal{J}+\mu\mathcal{F}}(\mathbf{e}_1, \mathbf{e}_2) = \lambda D_{\mathcal{J}}(\mathbf{e}_1, \mathbf{e}_2) + \mu D_{\mathcal{F}}(\mathbf{e}_1, \mathbf{e}_2), \text{ for } \mathcal{J}, \mathcal{F} \text{ convex and } (\lambda, \mu) \text{ positive.} $$

Nevertheless, the original notion of Bregman divergence cannot be used for building a gap functional between the solutions of the two auxiliary problems, mainly because of two limitations. The first one is the necessary differentiability of the generating function $\mathcal{J}$ in the definition, which is disqualifying for using non differentiable dissipation potentials as generating functions. The second one is that the Bregman divergence addresses distance between points in vector
spaces whereas it is of interest to deal rather with couples of dual quantities for application in mechanics. Actually the pertinent fields involved in the auxiliary problems are pairs such as $(\sigma, \epsilon)$ and $(\mathbf{A}, \alpha)$, where $\mathbf{A}$ is a thermodynamic force and $\alpha$ an internal variable, each element of the couple being related to the other by a relation involving a potential and incorporating derivatives or subdifferentials, for example:

$$ (14) \qquad \sigma = \frac{\partial \varphi}{\partial \epsilon}(\epsilon, \alpha) \quad \dot{\alpha} \in \partial \psi(\mathbf{A}) $$

$\varphi$ and $\psi$ are free energy and potential of dissipation, respectively. We then define the Bregman Gap as follows:

**Definition (Bregman Gap (BG)).** Let $\mathcal{J}$ be a proper convex function, $\mathcal{J}(\mathbf{e}) : \mathbb{R}^n \to \mathbb{R}$, the Bregman Gap, generated by $\mathcal{J}$, between a point $\mathbf{e}_1$ and the couple $(\mathbf{e}_2, \mathbf{p}_2)$, is the non-negative scalar:

$$ (15) \qquad \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_1, [\mathbf{e}_2, \mathbf{p}_2]) = \mathcal{J}(\mathbf{e}_1) - \mathcal{J}(\mathbf{e}_2) - \langle \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle $$

where $\mathbf{p}_2$ is a sub-gradient of $\mathcal{J}$ at point $\mathbf{e}_2$, and $\partial \mathcal{J}(\mathbf{e})$ is the subdifferential of $\mathcal{J}$ at point $\mathbf{e}$:

$$ \partial \mathcal{J}(\mathbf{e}) = \{\mathbf{p}, \mathcal{J}(\mathbf{d}) \geq \mathcal{J}(\mathbf{e}) + \langle \mathbf{p}, \mathbf{d} - \mathbf{e} \rangle \quad \forall \mathbf{d} \in \mathrm{dom}(\mathcal{J})\}. $$

Because of the lack of symmetry of the Bregman Gap, it is preferable to symmetrize it in order to build a gap between two pair $(\mathbf{e}_i, \mathbf{p}_i)$ and to get a similar role of the two solutions within the functional $E$.

**Definition (Symmetrized Bregman Gap).** Let $\mathcal{J}$ be a proper convex function, $\mathcal{J}(\mathbf{e}) : \mathbb{R}^n \to \mathbb{R}$, the symmetrized Bregman gap, generated by $\mathcal{J}$, between two couples $(\mathbf{e}_1, \mathbf{p}_1)$ and $(\mathbf{e}_2, \mathbf{p}_2)$, with $\mathbf{p}_1 \in \partial \mathcal{J}(\mathbf{e}_1)$ and $\mathbf{p}_2 \in \partial \mathcal{J}(\mathbf{e}_2)$ is:

$$ (16) \qquad \mathrm{BG}_{\mathcal{J}}^{s}([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) = \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_1, [\mathbf{e}_2, \mathbf{p}_2]) + \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_2, [\mathbf{e}_1, \mathbf{p}_1]) $$

It is worth noticing that if the generating function $\mathcal{J}$ is differentiable, then the symmetrized Bregman divergence and the symmetrized Bregman gap will coincide,

$$ (17) \qquad \mathcal{J} \text{ differentiable} \Rightarrow \mathrm{BG}_{\mathcal{J}}^{s}([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) = D_{\mathcal{J}}(\mathbf{e}_1, \mathbf{e}_2) + D_{\mathcal{J}}(\mathbf{e}_2, \mathbf{e}_1) $$

Furthermore for twice differentiable generating functions, the Bregman Gap takes the form:

$$ D_{\mathcal{J}}^{s}(\mathbf{e}_1, \mathbf{e}_2) = (\mathbf{e}_1 - \mathbf{e}_2)^{t} \nabla^2 \mathcal{J}(\mathbf{e}) (\mathbf{e}_1 - \mathbf{e}_2), \quad \mathbf{e} \in [\mathbf{e}_1, \mathbf{e}_2] $$

More important for application to the Cauchy problem are the following properties.

$$ \bullet \mathrm{BG}_{\mathcal{J}}^{s}([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) = \langle \mathbf{p}_1 - \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle $$

$$ \bullet \text{ If } \mathcal{J} \text{ is a quadratic function then: } \mathrm{BG}_{\mathcal{J}}^{s}([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) = 2\mathcal{J}(\mathbf{e}_1 - \mathbf{e}_2) $$

Then provided a generating function $\mathcal{J}$ can be identified as a function of the state variables associated to the solution fields $\mathbf{u}_1$ and $\mathbf{u}_2$ of the two auxiliary problems, we can generally define the gap functional as:

$$ (18) \qquad \begin{aligned} E(\mathbf{u}_1, \mathbf{u}_2) &= \int_{\Omega} \mathrm{BG}_{\mathcal{J}}^{s}([\mathbf{e}(\mathbf{u}_1), \mathbf{p}(\mathbf{u}_1)], [\mathbf{e}(\mathbf{u}_2), \mathbf{p}(\mathbf{u}_2)]) d\Omega \\ &= \int_{\Omega} \langle \mathbf{p}(\mathbf{u}_1) - \mathbf{p}(\mathbf{u}_2), \mathbf{e}(\mathbf{u}_1) - \mathbf{e}(\mathbf{u}_2) \rangle d\Omega \end{aligned} $$

or equivalently in the case of quadratic function $\mathcal{J}$:

$$ (19) \qquad E(\mathbf{u}_1, \mathbf{u}_2) = \int_{\Omega} \mathcal{J}(\mathbf{e}(\mathbf{u}_1) - \mathbf{e}(\mathbf{u}_2)) d\Omega $$

Note that the generating function $\mathcal{J}$ does not appear anymore in eq. (18), but its existence is mandatory in order ensure that the gap $\langle \mathbf{p}_1 - \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle$ is always positive. In contact mechanics, natural generating functions can be found with thermodynamics potentials: namely the free (Helmholtz) energy associated with the linear elastic behavior of the solids, and the pseudo-dissipation potential associated with the dissipative contact conditions. The Generalized Standard Material (GSM) formulation [3] can be used to identify the state variables and the corresponding potentials. Nevertheless, as the Coulomb law of friction cannot be recast into the GSM formulation, this case is addressed separately in the sequel.
## 5. TWO SOLUTION ALGORITHMS FOR THE CAUCHY PROBLEM

The Cauchy problem can now be solved by defining the two auxiliary problems ($P_1, P_2$) parametrized by the Dirichlet data $\eta$ on $\Gamma_u$, and by minimizing the gap obtained, thanks to a suitable symmetrized Bregman gap between the solutions of these two problems, as a function of $\eta$:

$$ (20) \qquad \min_{\eta} f_{obj}(\eta), \quad f_{obj}(\eta) = E(\mathbf{u}_1(\eta), \mathbf{u}_2(\eta)) $$

As described before [25, 26], a way to minimize $f_{obj}$ is to use trust region method with conjugate gradient process and a computation of the gradient by an adjoint method. Then each iteration of this algorithm necessitates the computation of the two solutions $\mathbf{u}_1$ and $\mathbf{u}_2$ and the computation of two adjoint problems with solutions $\mathbf{u}_1^*$ and $\mathbf{u}_2^*$. This algorithm has been used previously with some efficiency in various 2D and 3D Cauchy problems and this general approach compares very favorably with existing method as the KMF algorithm [12], by demanding a number of iterations to convergence of several orders of magnitude lesser. The KMF algorithm, initially presented as a fixed point algorithm, has been reinterpreted in the framework of the present method as an alternating descent algorithm, which explains its poor performance [33, 40].

Nevertheless, to address more nonlinear problems and also problems where singularities can appear in the unknown data $\eta$, as it is the case where the boundary $\Gamma_u$ contains cracks, a new formulation has been designed. For the sake of clarity, it will be described below on the model problem of elasticity. The two auxiliary problems read then:

Problem $P_1$

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_1) = 0 & \text{in } \Omega \\ \sigma_1 = \mathbb{A} : \epsilon(\mathbf{u}_1) & \text{in } \Omega \\ \epsilon(\mathbf{u}_1) = \nabla^s(\mathbf{u}_1) & \text{in } \Omega \\ \sigma_1 \cdot \mathbf{n} = \mathbf{F}_m & \text{on } \Gamma_m \\ \mathbf{u}_1 = \eta & \text{on } \Gamma_u \end{array} \right. $$

<figure><img src="image_4.png" /></figure>

Problem $P_2$

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_2) = 0 & \text{in } \Omega \\ \sigma_2 = \mathbb{A} : \epsilon(\mathbf{u}_2) & \text{in } \Omega \\ \epsilon(\mathbf{u}_2) = \nabla^s(\mathbf{u}_2) & \text{in } \Omega \\ \mathbf{u}_2 = \mathbf{U}_m & \text{on } \Gamma_m \\ \mathbf{u}_2 = \eta & \text{on } \Gamma_u \end{array} \right. $$

<figure><img src="image_5.png" /></figure>

The gap function associated with the symmetrized Bregman Gap generated by the function $J(\nabla \mathbf{u}) = \frac{1}{2}(\sigma(\mathbf{u}) : \epsilon(\mathbf{u}))$ is:

$$ (21) \qquad f_{obj}(\eta) = E(\mathbf{u}_1(\eta), \mathbf{u}_2(\eta)) = \int_{\Omega} (\sigma_1(\eta) - \sigma_2(\eta)) : (\epsilon_1(\eta) - \epsilon_2(\eta)) \, dV $$

and the solution of the Cauchy problem is obtained by minimizing $f_{obj}$. To derive an alternative formulation, the following auxiliary problem $P_3$ is defined, depending on the solution $\mathbf{u}_1$ of problem $P_1$:

Problem $P_3$

$$ \left\{ \begin{array}{ll} \text{div}(\sigma_3) = 0 & \text{in } \Omega \\ \sigma_3 = \mathbb{A} : \epsilon(\mathbf{u}_3) & \text{in } \Omega \\ \epsilon(\mathbf{u}_3) = \nabla^s(\mathbf{u}_3) & \text{in } \Omega \\ \mathbf{u}_3 = \mathbf{U}_m & \text{on } \Gamma_m \\ \sigma_3 \cdot \mathbf{n} = \sigma_1 \cdot \mathbf{n} & \text{on } \Gamma_u \end{array} \right. $$

<figure><img src="image_6.png" /></figure>

where the Dirichlet condition on $\Gamma_u$ in problem $P_2$ is replaced with a Neumann condition involving the solution of the problem $P_1$ which depends on $\eta$, so that $\mathbf{u}_3$ is a function of $\eta$ as well.
A new functional of the Dirichlet data $\eta$ on $\Gamma_u$ can be defined as:

$$ (22) \qquad g_{obj}(\eta) = E(\mathbf{u}_1(\eta), \mathbf{u}_3(\eta)) = \int_{\Omega} (\sigma_1(\eta) - \sigma_3(\eta)) : (\epsilon_1(\eta) - \epsilon_3(\eta)) \, dV $$

Some interesting properties can be derived making the minimization of $g_{obj}$ an equivalent solution method of the Cauchy problem. Firstly, the field $\mathbf{u}_3$ solution of $\mathcal{P}_3$ can be characterized as a minimum of the function $E$ on $V_m$:

$$ (23) \qquad V_m = \{ \mathbf{v} \in H^1(\Omega), \mathbf{v} = \mathbf{U}_m \text{ on } \Gamma_m \} $$

To prove this assertion, the stationary point $\mathbf{w}$ of the convex function $E(\mathbf{u}_1(\eta), \mathbf{v})$ is studied on $V_m$:

$$ (24) \qquad D_{\mathbf{v}} E(\mathbf{u}_1(\eta), \mathbf{w}) \cdot \delta\mathbf{v} = 0 \quad \forall \delta\mathbf{v} \in V_0 = \{ \mathbf{v} \in H^1(\Omega), \mathbf{v} = \mathbf{0} \text{ on } \Gamma_m \} $$

this leads to:

$$ \int_{\Omega} (\sigma(\mathbf{u}_1) - \sigma(\mathbf{w})) \cdot \epsilon(\delta\mathbf{v}) dV = 0 \quad \forall \delta\mathbf{v} \in V_0 $$

Notice that $\sigma(\mathbf{w}) = \mathbb{A} : \epsilon(\mathbf{w})$, $\sigma(\mathbf{u}_1) = \mathbb{A} : \epsilon(\mathbf{u}_1)$ and $\sigma(\delta\mathbf{v}) = \mathbb{A} : \epsilon(\delta\mathbf{v})$ and using the fact that $div(\sigma_1) = 0$ in $\mathcal{P}_1$ one obtains:

$$ \int_{\Omega} \mathrm{div}(\sigma(\mathbf{w}) \delta\mathbf{v}) dV - \int_{\Gamma_u} (\sigma(\mathbf{w}) - \sigma(\mathbf{u}_1)) \cdot \mathbf{n} \delta\mathbf{v} dS = 0 \quad \forall \delta\mathbf{v} \in V_0 $$

hence it can be deduced that $\mathrm{div}(\sigma(\mathbf{w})) = 0$ and $\sigma(\mathbf{w}) \cdot \mathbf{n} = \sigma(\mathbf{u}_1) \cdot \mathbf{n}$ on $\Gamma_u$. These two last equalities and the condition $\mathbf{w} \in V_m$ just as $\mathbf{u}_3$ show that $\mathbf{w} = \mathbf{u}_3$. We have then:

$$ (25) \qquad \mathbf{u}_3 = \arg \min_{\mathbf{v} \in V_m} E(\mathbf{u}_1(\eta), \mathbf{v}) $$

Thence, it is now straightforward to show following relations:

$$ i. \quad g_{obj}(\eta) \le f_{obj}(\eta) \quad \forall \eta $$

$$ ii. \quad \frac{\partial g_{obj}}{\partial \eta}(\eta) \cdot \delta \eta = D_{\mathbf{u}_1} E(\mathbf{u}_1(\eta), \mathbf{u}_3(\eta)) \cdot (D_{\eta} \mathbf{u}_1(\eta) \cdot \delta \eta) $$

the inequality (i.) can be deduced from the minimum property of the field $\mathbf{u}_3$, see (24). Indeed, as $\mathbf{u}_2 \in V_m$:

$$ g_{obj}(\eta) = E(\mathbf{u}_1(\eta), \mathbf{u}_3(\eta)) \le E(\mathbf{u}_1(\eta), \mathbf{u}_2(\eta)) = f_{obj}(\eta) $$

The relation (ii.) is derived by:

$$ \frac{\partial g_{obj}}{\partial \eta}(\eta) \cdot \delta \eta = D_{\mathbf{u}_1} E(\mathbf{u}_1(\eta), \mathbf{u}_3(\eta)) \cdot (\mathcal{D}_{\eta} \mathbf{u}_1(\eta) \cdot \delta \eta) + D_{\mathbf{u}_3} E(\mathbf{u}_1(\eta), \mathbf{u}_3(\eta)) \cdot \underbrace{(\mathcal{D}_{\eta} \mathbf{u}_3(\eta) \cdot \delta \eta)}_{\delta \eta \in V_0} $$

Remark that the second term vanishes because of the stationarity definition of $E(\mathbf{u}_1(\eta), \mathbf{v}(\eta))$ for $\mathbf{v} = \mathbf{u}_3$ as shown in (24). Then, the term $(\mathcal{D}_{\eta} \mathbf{u}_3(\eta) \cdot \delta \eta)$ vanishes on $\Gamma_m$ for all $\delta \eta$, hence $(\mathcal{D}_{\eta} \mathbf{u}_3(\eta) \cdot \delta \eta) \in V_0$.

Consequently, as $g_{obj}(\eta^{opt}) = 0 = f_{obj}(\eta^{opt})$, one can use the function $g_{obj}$ instead of $f_{obj}$ in order to solve the Cauchy problem by minimization. Secondly the second property (ii.) makes less expensive the computational cost of gradient of $g_{obj}$. In fact, only the adjoin problem associated to the field problem $\mathcal{P}_1$ has to be computed.

# 6. BREGMAN GAPS IN CONTACT MECHANICS

## 6.1. Linear elasticity with unilateral contact and Tresca friction law.
When the solids under scrutiny are linearly elastic (in the context of small transformations) and when the friction law can be cast into the formalism of Standard Generalized Materials (GSM), then the bulk Helmholtz free energy and the surface pseudo-potential of dissipation can naturally be used as generating functions of the Bregman Gap appearing in the general form (16) of the gap
functional $E(\mathbf{u}_1, \mathbf{u}_2)$. As an example, the Tresca friction law is associated with the following pseudo-potential of dissipation:

$$ (26) \qquad \mathcal{D}(\dot{\mathbf{u}}_T) = \tau \|\dot{\mathbf{u}}_T\| $$

the dual variable of $\dot{\mathbf{u}}_T$ being $\mathbf{R}_T$ and the friction constitutive equation reads: $\mathbf{R}_T \in \partial\mathcal{D}(\dot{\mathbf{u}}_T)$. In the incremental formulation, a linear combination of the Helmholtz free energy, reduced here to the elastic energy density, $\varphi(\epsilon(\Delta\mathbf{u}))$ and of the pseudo-potential of dissipation $\mathcal{D}(\Delta\mathbf{u}_T)$ leads to the gap functional:

$$ (27) \qquad E(\Delta\mathbf{u}_1, \Delta\mathbf{u}_2) = \chi \int_{\Omega} (\Delta\sigma_1 - \Delta\sigma_2) : (\epsilon(\Delta\mathbf{u}_1) - \epsilon(\Delta\mathbf{u}_2)) \, d\Omega \\ + (1-\chi) \int_{\Gamma_c} (\mathbf{R}_{T1} - \mathbf{R}_{T2}) \cdot (\Delta\mathbf{u}_{T1} - \Delta\mathbf{u}_{T2}) \, d\Gamma $$

with $0 < \chi \le 1$, the zero value of $\chi$ has to be excluded because of the non-strict convexity of $\mathcal{D}$. The first (bulk) part of the functional $E$ can be transformed into a surface integral, leading then to the final expression:

$$ (28) \qquad f_{obj}(\eta) = E(\Delta\mathbf{u}_1, \Delta\mathbf{u}_2) = \chi \int_{\Gamma_m} (\Delta\mathbf{F}_m - \Delta\sigma_2 \cdot \mathbf{n}) \cdot (\Delta\mathbf{u}_1 - \Delta\mathbf{U}_m) \, d\Gamma \\ + \int_{\Gamma_c} (\mathbf{R}_{T1} - \mathbf{R}_{T2}) \cdot (\Delta\mathbf{u}_{T1} - \Delta\mathbf{u}_{T2}) \, d\Gamma $$

where $\Delta\mathbf{u}_1$ and $\Delta\mathbf{u}_2$ depend on $\eta$ via the two problems $(\mathcal{P}_1, \mathcal{P}_2)$.

## 6.2. Linear elasticity with unilateral contact and Coulomb friction.
As the Coulomb friction criterion does not fall into the framework of Generalized Standard Materials, a generalization of the Bregman Gap definition is necessary. An extension of the theory of GSM has been proposed by De Saxcé [47], encompassing internal variables evolution laws which do not enjoy the normality condition like Armstrong-Frederick hardening plasticity, Drücker-Prager non associated plasticity, Clam-Clay soil constitutive relations, as well as the Coulomb friction law [5]. The notion of Implicit Standard Generalized Materials is introduced using what is named bi-potentials. More precisely [48, 49], if $\mathbf{x}$ and $\mathbf{y}$ are two elements of real vectors spaces $X$ and $Y$ respectively, the duality product between them denoted by $\langle \bullet, \bullet \rangle$ having an energetic interpretation, a real function $b$ on $X \times Y$ is called a bi-potential if and only if :

i. $b$ is separately convex and lower semi-continuous for $\mathbf{x}$ and $\mathbf{y}$;

ii. $\inf_{\mathbf{y} \in \mathcal{Y}} b(\mathbf{x}, \mathbf{y}) \in \mathbb{R}^+ \cup \{+\infty\} \forall \mathbf{x} \in X$, $\inf_{\mathbf{x} \in X} b(\mathbf{x}, \mathbf{y}) \in \mathbb{R}^+ \cup \{+\infty\} \forall \mathbf{y} \in Y$,

iii. the function $b$ fulfills the generalized Legendre-Fenchel property:

$$ \text{a. } b(\mathbf{x}, \mathbf{y}) \ge \langle \mathbf{x}, \mathbf{y} \rangle \quad \forall (\mathbf{x}, \mathbf{y}) \in X \times Y $$

$$ \text{b. } b(\mathbf{x}, \mathbf{y}) = \langle \mathbf{x}, \mathbf{y} \rangle \Leftrightarrow \mathbf{x} \in \partial b_{\mathbf{y}}(\mathbf{x}, \mathbf{y}) \Leftrightarrow \mathbf{y} \in \partial b_{\mathbf{x}}(\mathbf{x}, \mathbf{y}) $$

Obviously, the GSM framework is recovered as soon as the bi-potential is a separated variables function. The generalized Bregman Gap and its symmetrized counterpart can now be defined as:

**Definition.** The generalized Bregman Gap generated by the bi-potential $b(\mathbf{x}, \mathbf{y})$, between a point $\mathbf{e}_1$ and the couple $(\mathbf{e}_2, \mathbf{p}_2)$ with $\mathbf{p}_2 \in \partial\mathcal{J}(\mathbf{e}_2)$, is the non-negative scalar:

$$ \begin{aligned} \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_1, [\mathbf{e}_2, \mathbf{p}_2]) &= b(\mathbf{e}_1, \mathbf{p}_2) - b(\mathbf{e}_2, \mathbf{p}_2) - \langle \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle \\ &\equiv b(\mathbf{e}_1, \mathbf{p}_2) - \langle \mathbf{p}_2, \mathbf{e}_1 \rangle \quad \mathbf{p}_2 \in \partial_{\mathbf{y}} b(\mathbf{e}_2, \mathbf{p}_2) \end{aligned} $$

The last expression shows that, thanks to the definition of the bi-potentials, the generalized Bregman Gap is always positive. The symmetrized generalized Bregman Gap is:

$$ (29) \qquad \begin{aligned} \mathrm{BG}_{\mathcal{J}}^s([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) &= \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_1, [\mathbf{e}_2, \mathbf{p}_2]) + \mathrm{BG}_{\mathcal{J}}(\mathbf{e}_2, [\mathbf{e}_1, \mathbf{p}_1]) \\ &\equiv b(\mathbf{e}_1, \mathbf{p}_2) + b(\mathbf{e}_2, \mathbf{p}_1) - \langle \mathbf{p}_2, \mathbf{e}_1 \rangle - \langle \mathbf{p}_1, \mathbf{e}_2 \rangle \end{aligned} $$

but this expression cannot be reduced, as it is the case for SGM, to the cross-product $\langle \mathbf{p}_1 - \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle$, because the positivity of this cross-product is equivalent to the normality rule, which is precisely violated for non-Standard Generalized Materials.

$$ (30) \qquad \mathrm{BG}_{\mathcal{J}}^s([\mathbf{e}_1, \mathbf{p}_1], [\mathbf{e}_2, \mathbf{p}_2]) = \langle \mathbf{p}_1 - \mathbf{p}_2, \mathbf{e}_1 - \mathbf{e}_2 \rangle + b(\mathbf{e}_1, \mathbf{p}_2) + b(\mathbf{e}_2, \mathbf{p}_1) \\ - \langle \mathbf{p}_1, \mathbf{e}_1 \rangle - \langle \mathbf{p}_2, \mathbf{e}_2 \rangle $$
It can be shown that the Coulomb friction law can be defined using the following bi-potential [48, 49], the dissipation being given by $\mathbf{R}_T \dot{\mathbf{u}}_T = \lambda \|\mathbf{R}_T\|^2$ when friction occurs:

$$ (31) \qquad b(R_N, \mathbf{R}_T; \dot{u}_N, \dot{\mathbf{u}}_T) = I_{\{u_N\}}^{+} + I_{\{-R_N\}}^{+} + I_{\{-\rho R_N - \|\mathbf{R}_T\|\}}^{+} + \rho_0 R_N \|\dot{\mathbf{u}}_T\| $$

where $I^+$ is the indicator function of the set $[0, +\infty[$. For the incremental formulation, the bi-potential reads:

$$ (32) \qquad b(\Delta R_N, \Delta \mathbf{R}_T; \Delta u_N, \Delta \mathbf{u}_T) = I_{\{\Delta u_N\}}^{+} + I_{\{-R_N - \Delta R_N\}}^{+} + I_{\{-\rho(R_N + \Delta R_N) - \|\mathbf{R}_T + \Delta \mathbf{R}_T\|\}}^{+} + \rho_0 (R_N + \Delta R_N) \|\Delta \mathbf{u}_T\| $$

Because the two solutions of the auxiliary problems satisfy the contact condition and the Coulomb's friction yield function, the gap functional in elastostatics with unilateral contact and Coulomb's friction is simply:

$$ (33) \qquad E(\Delta \mathbf{u}_1, \Delta \mathbf{u}_2) = \chi \int_{\Omega} (\Delta \sigma_1 - \Delta \sigma_2) : (\epsilon(\Delta \mathbf{u}_1) - \epsilon(\Delta \mathbf{u}_2)) d\Omega \\ -(1-\chi) \int_{\Gamma_c} \rho R_N \cdot (\|\Delta \mathbf{u}_{T1}\| + \|\Delta \mathbf{u}_{T2}\|) d\Gamma \\ -(1-\chi) \int_{\Gamma_c} [\Delta R_{N1}(\rho \|\Delta \mathbf{u}_{T2}\| + \Delta \mathbf{u}_{T1}) + \Delta R_{N2}(\rho \|\Delta \mathbf{u}_{T1}\| + \Delta \mathbf{u}_{T2})] d\Gamma $$

with $0 < \chi \le 1$, the zero value of $\chi$ has to be excluded because of the non-strict convexity of $b$. The first (bulk) part of the functional $E$ can be transformed into a surface integral, leading then to the final expression:

$$ (34) \qquad f_{obj}(\eta) = \chi \int_{\Gamma_m} (\Delta \mathbf{F}_m - \Delta \sigma_2 \cdot \mathbf{n}) \cdot (\Delta \mathbf{u}_1 - \Delta \mathbf{U}_m) \, d\Gamma \\ -(1-\chi) \int_{\Gamma_c} \rho R_N \cdot (\|\Delta \mathbf{u}_{T1}\| + \|\Delta \mathbf{u}_{T2}\|) \, d\Gamma \\ -(1-\chi) \int_{\Gamma_c} [\Delta R_{N1}(\rho \|\Delta \mathbf{u}_{T2}\| + \Delta \mathbf{u}_{T1}) + \Delta R_{N2}(\rho \|\Delta \mathbf{u}_{T1}\| + \Delta \mathbf{u}_{T2})] \, d\Gamma $$

where $\Delta \mathbf{u}_1$ and $\Delta \mathbf{u}_2$ depend on $\eta$ via the two problems $(\mathcal{P}_1, \mathcal{P}_2)$.

## 7. APPLICATION TO SOME IDENTIFICATION PROBLEMS IN CONTACT MECHANICS

As can be seen, the functional $E$ in both cases is not differentiable at points where the tangential displacement increment vanishes, because the norm of the increment appears in its expression. However, thanks to the convexity of $E$, the subdifferential at these points is well defined although not reduced to a singleton. Because the set of these points, corresponding to the limit of the sliding zone, is a set of zero measure, we can pick arbitrarily an element in the subdifferential (namely 0) in the process of derivation of the gradient of $E$ without degrading significantly the overall quality of it.

### 7.1. Identification of missing boundary conditions.
In order to show the efficiency and robustness of the presented method, we consider the example of two-dimensional cracked elastic solid submitted to bending load, as shown on figure 4. The crack buried into the solid brings additional difficulties because of the singularity of the solution in its vicinity for linear elastic materials, but also presents special feature when contact and friction are taken into account such as canceling some of the singularities. This makes this illustration a non-trivial one. The domain is rectangular ($10 \, m \times 20 \, m$) with the assumption of plane strains. The material has an elastic behavior with Young modulus $E = 10^6 \, Pa$ and Poisson coefficient $\nu = 0.3$. The crack is at the middle of the domain $[-a, a]$ on the $x$ axis with $a = 2.5 \, m$. The solid is submitted to a two kind of boundary conditions :

(a) Distributed linear loading $f_y = 1000 \, Pa$ and $f_x = 0$ on the top boundary $\Gamma_{top}$ and $f_y = -1000 \, Pa$ and $f_x = 0$ on the bottom boundary $\Gamma_{bottom}$. To ensure the symmetry of the response conditions ($u_y = 0$) and $u_x = 0$ are assumed on the point $B$, then $u_y = 0$ on the point $A$. This loading exerts a compression on the right half of domain, a traction on other half. Thanks to this loading, it results: the closing of the crack, in the vicinity crack tip $(a, 0)$, contained in half of the domain in compression; the opening of the crack, in the vicinity crack tip $(-a, 0)$, contained in half of the tensile domain. This case is built in order to simulate a crack in an infinite domain, for which a closed form solution exists.
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Geometry of cracked solid and boundary conditions of forward problems.</figcaption></figure>

(b) Distributed linear loading $f_y = 1000 \times Pa$ and $f_x = 100(1-x^2/25)Pa$ on the top boundary $\Gamma_{top}$ and conditions ($u_x = 0, u_y = 0$) on the boundary $\Gamma_{bottom}$. This loading exerts a compression on the right half of the domain, a traction on other half associated to shearing effect. Thanks to this loading, it results: the closing of the crack, in the vicinity crack tip $(a, 0)$, contained in half of the domain in compression with stick-slip in this part of the crack; the opening of the crack, in the vicinity crack tip $(-a, 0)$, contained in half of the tensile domain.

The software used to compute forward and inverse problems is CODE-ASTER© [50] which is an open source distribution. CODE-ASTER is coupled with the numerical optimization toolbox of SciPy [51] and the Optimization Toolbox of Matlab [52].

The above forward problems are used to generate synthetic data ($\mathbf{U}_m, \mathbf{F}_m$), that will be used for the superabundant boundary conditions to solve inverse problems. Hereafter, we consider two inverse problems :

(A) the first one involving friction-less contact and linear forward problems defined on the sub-domain shown on figure 6a. The Cauchy data ($\mathbf{F}_m, \mathbf{U}_m$) are extracted from the forward problem defined in figure 4a.

(B) the second one involving contact with Coulomb friction $\mu = 0.3$ and nonlinear forward problems defined on whole domain shown on figure 6b. The Cauchy data ($\mathbf{F}_m, \mathbf{U}_m$) are extracted from the forward problem defined in figure 4b.

In order to avoid the "Inverse Problem Crime" as defined in [53], two different meshes are used for producing the reference solution and computing the inverse problems and fracture parameters. The difference between the two interpolations of the "data" in the two meshes, introduce also a form of small noise into the problem.

Figure 5a shows the distribution of Von-Mises stress around the crack resulting from the forward problem (a). Figure 5b shows the distribution of Von-Mises stress around the crack
<figure><img src="image_8.png" /><figcaption>FIGURE 5. Von-Mises stress distribution obtained for the two forward problems.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 6. Geometry and boundaries used in inverse problems.</figcaption></figure>

resulting from the forward problem (b), and where we can easily observe the sliding phenomenon on the crack.

7.1.1. *First example: linear inverse problem (A).* In this inverse problem we consider only the top half of the domain as shown on figure 6a. Using definitions given in section 2.2, we denote by $\Gamma_u = \Gamma \cup \Gamma_{top}$ the boundary where data are missing; by $\Gamma_m = \Gamma_r \cup \Gamma_l$ the boundary where data are overspecified. Thus, we assumed that:
* the interface $\Gamma$ where the crack is located is known, whereas the crack's position and its length on this interface are unknown and have to be identified.

* boundary conditions on $\Gamma_u$ are missing and have to be recovered.

* on $\Gamma_m$ overspecified data are given: displacement components are given or measured $\mathbf{u} = \mathbf{U}_m$, and as this boundary is free of loading Neumann boundary condition is also known: $\mathbf{F}_m = \sigma(\mathbf{u}) \cdot \mathbf{n} = \mathbf{0}$.

The components of $\mathbf{U}_m$ are extracted from the forward problems described in the item (a) of section 7.2. Recall that in these cases the Neumann condition $\mathbf{F}_m = \mathbf{0}$ is exactly known because $\Gamma_m$ is free of load. While the Dirichlet condition $\mathbf{U}_m$ extracted on the nodes of forward problem mesh is noisy. The noise rate depends on the differences between the mesh used in the direct problem and the one used in the inverse problem. Figures 8 show the identified data compared to reference ones. These results can be exploited as follows:

i. One can observe that displacement and surface traction components on the top boundary are both close to the reference values, see figures 8a and 8b.

ii. On the bottom boundary, we observe that the identified data present important errors in the vicinity of the crack tips. However the opened part of the crack can be easily identified by the displacement gap of $U_y^{Id}$ as shown on figure 8c. It is straightforward to estimate that a tip crack is located in the vicinity of $x = 2.5 m$.

iii. Using these identified data and above information, one can compute the stress intensity factor $K_I$ for Mode-I loading by means of contour integral for energy release rate $G$ (J-Integral). Notice that this computation can be hold without knowing precisely the location of the crack tip. It is only necessary to know the area where it is located. Following, J-Integral formula takes into account the symmetry of this case.

$$ (35) \qquad G = \frac{K_I^2(1 - \nu^2)}{E} = 2 \int_C W \mathbf{n}_x - \mathbf{T}_i \frac{\partial U_i}{\partial x}, \quad \mathbf{T} = \sigma \cdot \mathbf{n}. $$

As J-Integral is known to be path-independent, $C$ can be any boundary defining region surrounding the supposed crack tip, $\mathbf{n}$ is the outward normal at the contour $C$ and $\mathbf{T}_i$ is surface traction component ($i = x, y$) on $C$. Nevertheless, because of the use of approximated reconstruction of the field, the choice of $C$ location was made in a manner to avoid the noisy area surrounding the crack tip. Indeed, even if the path-independent property remains true, the stability property is lost due to the noise affecting the identified field.

<figure><img src="image_10.png" /><figcaption>FIGURE 7. Distribution of Von Mises stress obtained from identified data of the linear inverse problem: case A.</figcaption></figure>

In this case, $C$ has rectangular shape, it begins at $x = -1.1m$ and finishes at $x = -4m$ with height of 3.9 m as shown on figure 7 with white line. The reference value of the stress intensity factor is $K_I \approx 3967 Pa\sqrt{m}$, it arises from a closed form solution in a semi-infinite medium, see [54] for more details. The identified value is $K_I^{Id} \approx 3960 Pa\sqrt{m}$, the relative error is then less than 1%.
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Identified data of Linear inverse problem: case A. Exponents "R" and "Id" denote reference and identified values.</figcaption></figure>

7.1.2. *Second example: nonlinear inverse problem (B).* In this inverse problem, the whole domain is considered, as shown on figure 6b. Using definitions given in section 2.2, we denote by $\Gamma_u = \Gamma_{top}$ the boundary where data are missing, by $\Gamma_m = \Gamma_r$ the boundary where data are overspecified and by $\Gamma_b = \Gamma_l \cup \Gamma_{bottom}$ the boundary where Neumann and Dirichlet conditions are prescribed. Thus, we assumed that:

* the crack's position and the Coulomb coefficient are known.

* boundary conditions on $\Gamma_u$ are missing and have to be recovered.

* on $\Gamma_m$ overspecified data are given: displacement components are given or measured $\mathbf{u} = \mathbf{U}_m$, and as this boundary is free of loading Neumann boundary condition is also known: $\mathbf{F}_m = \sigma(\mathbf{u}) \cdot \mathbf{n} = \mathbf{0}$.

* on $\Gamma_b$ the left boundary is free of load and on the bottom boundary Dirichlet condition is prescribed as shown on figure 6b.

The components of $\mathbf{U}_m$ are extracted from the forward problems described in the item (b) of section 7.2. Recall that in these cases the Neumann condition $\mathbf{F}_m = \mathbf{0}$ is exactly known because $\Gamma_m$ is free of load. While the Dirichlet condition $\mathbf{U}_m$ extracted on the nodes of the forward problem mesh is noisy. The noise rate depends on the differences between the mesh used in the direct problem and the one used in the inverse problem. Figures 9 and 10 show the identified data compared to the reference ones. These results can be exploited as follows:

i. The reader can observe that on figure 9a identified and reference values of displacement components on $\Gamma_u$ are quite close. However $U_x^{Id}$ is little bit more noisy than $U_y^{Id}$.
ii. Figure 9b shows identified and reference values of stress vector components on $\Gamma_u$. Notice here that both identified components are noisy. However tangential component $\sigma_{nt}^{Id}$ is much more noisy than the normal component $\sigma_{nn}^{Id}$.

In order to improve these data it is possible to use post-regularization using Total Variation denoising method [55], which was carried out by minimizing the following functions.

$$J_{TV} = \frac{1}{2} \int_{\Gamma_{top}} \left[ (\mathbf{X} - \mathbf{Y})^2 + \frac{\alpha}{2} \sqrt{(\nabla \mathbf{X})^2 + \varepsilon^2} \right]$$

where $\alpha$ is regularization parameter, $\varepsilon$ is a small parameters to avoid derivative singularity, $\mathbf{X}$ is the denoised field and $\mathbf{Y}$ is the noisy one. Figures 9c and 9d show reference and regularized-identified displacement and stress vector components. We observe that post-regularized data are in good agreement the reference ones.

However, the reader can notice that the noise in the identified stress vector components is huger than that in the displacements components because it is not only due to the noise present in the superabundant data but also to the stress computation process. Indeed, stress field is computed by derivative operation of the Finite Element Method (FEM). This is done on the boundary where displacement components are identified and are noisy. Consequently, the noise is amplified by the derivative procedure in the post-treatment step. In order to avoid this effect, one can use a fictive boundary $\Gamma_{top}^f$ inside the domain and parallel to $\Gamma_u$, where data are less noisy. So, instead of computing the missing data on $\Gamma_u$ they were computed on the fictive boundary $\Gamma_{top}^f$ as shown on figures 9e and 9f. In this case $\Gamma_{top}^f$ is chosen at $y = 9 m$ enough close to $\Gamma_u$, we observe that identified data are in good agreement with reference ones and present better precision than regularized ones. This trick can be used to minimize the sensitivity to the perturbation on $\Gamma_u$.

iii. Figures 10a to 10d show the reference and identified data on the top lip of the crack. Figure 10a shows the displacement component, figure 10b shows Coulomb criterion and contact indicator, figure 10c shows gap and slip values and figure 10c shows stress vector components and Coulomb criterion. The major remark here is that identified and reference values are identical due to the fact that this boundary is far from the one where data were missing and have to been identified. Consequently they are much less sensitive to noise on $\Gamma_u$.

From above observations, it can be concluded that by defining problem $P_1$ and $P_3$ using a fictitious boundary $\Gamma_u^f$ instead of $\Gamma_u$, one can significantly improve the quality and precision of the identified data on $\Gamma_u$. The choice of the fictitious boundary must be made judiciously in the vicinity of the real boundary $\Gamma_u$ taking into account the physical phenomena of the problem.

In order to show the influence of noise on the identification results, we introduced a mixed Gaussian-impulse noise in the Dirichlet condition $\mathbf{U}_m$ as follows :

$$\mathbf{U}_m^\delta = \mathbf{U}_m + \mathbf{U}_m(a \mathbf{G} + b \mathbf{I})$$

where: $a = 0.05$ is a percentage of the amplitude of white Gaussian noise vector represented by $\mathbf{G}$; $b = 0.05$ is a percentage of amplitude of random-valued impulse noise represented by the vector $\mathbf{I}$. The impulse noise is computed by assuming that 5% of data are corrupted, which are randomly chosen. The Neumann Condition $\mathbf{F}_m = 0$ remains exactly known because the surface is free of load. A stopping rule based on the estimate derived in [34] was employed in order to avoid numerical explosion during the optimization process.

Figure 11 shows the exact and identified data with noisy Cauchy data. We notice that within 5% of noise the identified components of displacement field and loading are in good agreement with the exact ones as shown in figures 11a and 11b. TV post-regularization was also carried out on these identified data. It improves the quality of the identified components as shown on figures 11c and 11d. However, on figures 11e and 11f we observe that identified parameters on the crack interface are slightly perturbed, they remain insensitive to the effect of this noise level.

7.2. **Identification of the friction coefficient.** This section is dedicated to the identification of the friction coefficient of the Coulomb's law, provided Cauchy data are given on the part $\Gamma_m$
<figure><img src="image_12.png" /><figcaption>(A) Displacement components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(C) Regularized displacement components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(E) Displacement components on fictive boundary $\Gamma_{top}^f$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(B) Stress vector components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(D) Regularized Stress vector components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(F) Stress vector components on fictive boundary $\Gamma_{top}^f$</figcaption></figure>

<caption>FIGURE 9. Identified data of nonlinear inverse problem on $\Gamma_u$: case B. Exponents "R" and "Id" denote reference and identified values.</caption>

of the boundary of the solid and data are missing on the boundary $\Gamma_u$ where contact occurs. As indicated before, the identification is performed once the Cauchy problem has been solved and when the mechanical fields on $\Gamma_c$ have been determined. The contact surface $\Gamma_c$ decomposes into three non-overlapping parts:

• $\Gamma_{nc}$ the non-contact or separation part: $R_N(x) = R_T(x) = 0$,

• $\Gamma_s$ the sliding part: $\rho_0 |R_N(x)| = \|R_T(x)\|$, $R_N(x) < 0$,

• $\Gamma_a$ the (friction) adhesion part: $\rho_0 |R_N(x)| > \|R_T(x)\|$, $R_N(x) < 0$,

$\rho_0$ being the actual (unknown) Coulomb's friction coefficient. A necessary condition for the identification problem to be solvable is of course that friction actually occurs on the contact
<figure><img src="image_18.png" /><figcaption>(A) Displacement components on $\Gamma_{Crack}$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(B) Coulomb criterion and Contact indicator $\Gamma_{Crack}$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(C) Gap and slip on $\Gamma_{Crack}$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(D) Stress vector components and Coulomb criterion on $\Gamma_{Crack}$</figcaption></figure>

<caption>FIGURE 10. Identified data of nonlinear inverse problem on $\Gamma_{Crack}$: case B. Exponents "R" and "Id" denote reference and identified values.</caption>

surface, that is $meas \Gamma_s \stackrel{\Delta}{=} |\Gamma_s| > 0$. Let us remark first that $\rho_0$ is the smallest scalar such that the Coulomb criterion is satisfied over the whole contact surface $\Gamma_c$:

$$ \rho_0 = \inf \{ \rho, \rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| \ge 0 \quad \forall \mathbf{x} \in \Gamma_c \} $$

For any positive scalar $\rho$, one has, on the partition of the contact surface $\Gamma_c$:

* on $\Gamma_{nc}$ (non-contact or separation part) $\rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| = 0 \forall \rho$

* on $\Gamma_s$ (sliding part) $$ \begin{cases} \rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| < 0 & \text{if } \rho < \rho_0 \\ \rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| = 0 & \text{if } \rho = \rho_0 \\ \rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| > 0 & \text{if } \rho > \rho_0 \end{cases} $$

* on $\Gamma_a$ (adhesion part) $\rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\| > 0 \forall \rho \ge \rho_0$

We define then $J(\rho)$ as the function:

$$ J(\rho) = \int_{\Gamma_c} H(\rho|R_N(x)| - \|\mathbf{R}_T(\mathbf{x})\|) \, d\Gamma, $$

where $H$ is the Heaviside function: $H(x) = 1$ if $x > 0$, and $H(x) = 0$ if $x \le 0$. It is straightforward to show the following properties for $J$.
<figure><img src="image_22.png" /><figcaption>(A) Displacement components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(B) Stress vector components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(C) Regularized displacement components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(D) Regularized Stress vector components on $\Gamma_{top}$</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(E) Gap and slip on $\Gamma_{Crack}$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(F) Stress vector components and Coulomb criterion on $\Gamma_{Crack}$</figcaption></figure>

<caption>FIGURE 11. Identified data of nonlinear inverse problem on $\Gamma_u :$ case B with noisy Cauchy data (5%). Exponents "R" and "Id" denote reference and identified values.</caption>

(1) $J$ is a positive monotonic function, and $J(0) = 0$

(2) For $\rho > \rho_0$, $J$ is constant equal to $|\Gamma_s| + |\Gamma_a|$

(3) $J$ is discontinuous at $\rho = \rho_0$ and $[J(\rho_0)] = |\Gamma_s|$

<figure><img src="image_28.png" /></figure>

Finding $\rho_0$ is then simply looking for the largest scalar such that $J$ is strictly lesser than $|\Gamma_s| + |\Gamma_a|$.
For that purpose a dichotomy algorithm can be used starting with a segment $[0, \rho_{max}]$, where $J(\rho_{max}) = |\Gamma_s| + |\Gamma_a|$. With such algorithm it is possible to define an a priori precision or equivalently a number of iterations, as function of the estimated noise on the data $\sigma \cdot \mathbf{n}$:

$$ pres = \frac{\Delta \rho}{\rho_{max}} \quad \text{or} \quad N_{max} = -\frac{\log(pres)}{\log(2)} \text{ with the estimation : } \Delta \rho \approx \frac{\rho_g n_{\sigma} + n_{\tau}}{|R_N(x)|_{max}} $$

where $\rho_g$ is a first guess of the friction coefficient and $(n_{\sigma}, n_{\tau})$ the noises on $|R_N(x)|$ and $||\mathbf{R}_T(x)||$ respectively.

**Remark.**

(1) Any measure $d\Gamma$ can be used, for example if the data is available at $M$ points $x_i$ on $\Gamma_c$:

$$ \int_{\Gamma_c} f(x)d\Gamma = \sum_{i=1}^{M} f(x_i) \quad \text{then} \quad J(\rho) = \sum_{i=1}^{M} H(\rho|R_N(x_i)| - ||\mathbf{R}_T(x_i)||) $$

(2) Once $\rho$ has been determined, the three zones in $\Gamma_c$ are also identified by using the value of the friction coefficient criterion or the simultaneous zero (or negligible) values of $|R_N(x_i)|$ and $||\mathbf{R}_T(x_i)||$.

(3) Being global, this approach should be less sensitive to the noise in the data than a direct estimation based on local computations of the ratio $||\mathbf{R}_T(x)|| / |R_N(x)|$.

(4) An alternative way for determining the friction coefficient, especially if the expected noise levels on the fields $|R_N(x_i)|$ and $||\mathbf{R}_T(x_i)||$ are high, is to use the total variation of the function $J$ between two succeeding points $i\Delta\rho$ and $(i+1)\Delta\rho$ and to select the friction coefficient by:

$$ \rho_{opt} = \frac{2i_{opt} + 1}{2} \Delta\rho \quad i_{opt} = \text{ArgMax} \left| J[(i+1)\Delta\rho] - J[i\Delta\rho] \right| $$

**Illustration.** In an axisymmetric situation, consider the following expressions of $\sigma = -R_n(x) = -\sigma_{zz}$ and $\tau = ||R_T(x)|| = |\sigma_{rz}|$ on the disc with radius $R$ and normal $e_z$ representing the surface of contact $\Gamma$:

$$ (36) \qquad \begin{cases} \sigma(r) = (1 - r^2)\chi_{[0, \varepsilon R]} \\ \tau(r) = \beta\rho_0 r\chi_{[0, \alpha(\beta)R]} + \rho_0(1 - r^2)\chi_{[\alpha(\beta)R, \varepsilon R]} \end{cases} \quad \text{with} \quad \alpha(\beta) = \frac{-\beta + \sqrt{\beta^2 + 4}}{2R} $$

which correspond to : $\Gamma_a = [0, \alpha R]$, $\Gamma_s = [\alpha R, \varepsilon R]$, $\Gamma_{nc} = [\varepsilon R, R]$, $\Gamma_C = [0, R]$. $\chi_\gamma$ is the characteristic function of the set $\gamma$.

Consider noise-free stress components $\sigma$ and $\tau$ as displayed on figure 12a. The function $J(\rho)$ and its variation $dJ(\rho)$ are displayed on the figure 12b, and exhibit as expected, a jump of amplitude $(\varepsilon - \alpha)R = 0.5858$ at the point $\rho = 0.5$. This jump is well highlighted by the variation of $J(\rho)$ denoted $dJ(\rho)$. If some noise is added to stress components as shown on figure

<figure><img src="image_29.png" /><figcaption>FIGURE 12. Exact and noisy stress components and function $J$ for the case where $R = 2$, $\beta = 2$, $\rho_0 = 0.5$, $\alpha = 0.2071$, $\varepsilon = 0.5$, with 10% of uniform white noise.</figcaption></figure>
12a, the noisy function $J^\delta(\rho)$ is by far more irregular and no more strictly monotonic. Consider for example, a partition of the segment $[0, R]$ consisting of $m = 41$ points regularly spaced, the function $J^\delta(\rho)$ is plotted on figure 12b for $p = 101$ points regularly spaced on the segment $[0, 1]$. Notice that the jump remains well highlighted by the variation of $J^\delta(\rho)$ denoted $dJ^\delta(\rho)$. The value of the friction coefficient identified by the jump of $J^\delta(\rho)$ is $\rho^{Id} = 0.51$.

**Application.** In this section we consider the linear inverse problem (A) with loading defined in (b) section. This case of loading generates traction-compression associated to shearing effect. The crack is then closed in the vicinity of the crack tip $(a, 0)$, a compression with stick-slip acting in this part of the crack while in the vicinity of the second crack tip $(-a, 0)$ the crack is open. The identified stress vector components are extracted on $\Gamma$ as shown on figure 6a where the crack is located. Figure 13a shows the normal and tangential stress components. The reference value are of course plotted only on the crack boundary $\Gamma_{crack}$ however the identified components are plotted on $\Gamma$. Figure 13b shows the reference function $J(\rho)^R$ computed on $\Gamma_{crack}$ which is compared to the function $J(\rho)^{Id}$ computed by using identified stress component on $\Gamma$. One can notice that the jump of $J(\rho)$ is located in the vicinity of the exact value of friction coefficient $\mu = 0.3$ which was used in the computational procedure of the forward problem. The identified value of the friction coefficient is $\rho^{Id} \in [0.29, 0.31]$.

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>FIGURE 13. Identified stress vector components and friction coefficient on $\Gamma$: case A load (b). Exponents "R" and "Id" denote reference and identified values.</figcaption></figure>

## 8. CONCLUSION

Cauchy problems in the context of contact and friction mechanics have been addressed in this paper. The approach used here is an extension of previous work in convex nonlinear mechanics where the minimization of an energy-like gap between the solutions of two auxiliary well-posed forward problems is the driving force behind the method. The extension to contact and friction necessitated to introduce a new concept of gap, namely the Bregman gap, derived from the Bregman divergence. Also, in order to deal with Coulomb friction, a further extension of the Bregman gap concept including implicit standard generalized material potentials generating functions has been proposed. Furthermore an improvement of the minimization algorithm has been described. Illustrations of the performance of the approach are provided in the case of cracked solids with unilateral contact and friction. Since Cauchy problem is severely (exponentially [56]) ill-posed, the greatest attention must be paid to the noise inevitably present in the data. In our problem, the most important source of noise comes from images correlation which provides the surface displacement fields that are used as Cauchy data. In order to avoid as much as possible that this type of noise is exponentially reflected in the results, pre-processing is mandatory. Another source of noise comes, when using synthetic data, from the interpolation between the mesh that produced the data and the one used to solve the Cauchy problem. Since this noise is directly reflected in the results, post-TV regularization is used on the results themselves. For the same reason, a technique of immersed boundary has also been used. Very good identification results
are then obtained for recovering missing boundary conditions and identifying parameters such as stress intensity factors. A new algorithm for the identification of the coulomb coefficient using the stress vector field on the contact area was also given.

## REFERENCES

[1] J. Haslinger, I. Hlaváček, and J. Necas. Numerical methods for unilateral problems in solid mechanics. In *Finite Element Methods (Part 2)*, Numerical Methods for Solids (Part 2), volume 4 of *Handbook of Numerical Analysis*, pages 313–485. Elsevier, 1996.

[2] P. Wriggers. *Computational contact mechanics*. Springer-Verlag Berlin Heidelberg, 2006.

[3] B. Halphen and Q. S. Nguyen. Sur les matériaux standards généralisés. *Journal de Mécanique*, 14:39-63, 1975.

[4] G. Amontons. De la résistance causée dans les machines. *Mémoires de l'Académie Royale A*, pages 257-282, 1699.

[5] C. A. Coulomb. *Théorie des machines simples*. Bachelier Paris, 1821.

[6] M. Fremond. *Contact with Adhesion*, pages 177-221. Springer Vienna, Vienna, 1988.

[7] M. Raous, L. Cangémi, and M. Cocu. A consistent model coupling adhesion, friction, and unilateral contact. *Computer Methods in Applied Mechanics and Engineering*, 177(3):383–399, 1999.

[8] A. Sutton, J.J. Orteu, and Schreier H. *Image Correlation for Shape, Motion and Deformation Measurements Basic Concepts, Theory and Applications*. Springer, 2009.

[9] S. Avril, M. Bonnet, A. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Latourte, D. Lemosses, S. Pagano, E. Pagnacco, and F. Pierron. Overview of identification methods of mechanical parameters based on full-fields measurements. *Experimental Mechanics*, 48:381-402, 2008.

[10] J. Hadamard. *Lectures on Cauchy's problem in linear partial differential equations*. Dover, New York, 1923.

[11] G. Fichera. Sul problema elastostatico di Signorini con ambigue condizioni al contorno. *Atti Accad. Naz. Lincei, VIII. Ser., Rend., Cl. Sci. Fis. Mat. Nat.*, 34:138-142, 1963.

[12] V. A. Kozlov, V. G. Maz'ya, and A. F. Fomin. An iterative method for solving the Cauchy problem for elliptic equations. *Comput. Math. Math. Phys.*, 31:45-52, January 1992.

[13] J. Baumeister and A. Leitão. On iterative methods for solving ill-posed problems modeled by partial differential equations. *Journal of Inverse and Ill-Posed Problems*, 9(1):13-29, 2001.

[14] L. Marin and D. Lesnic. Boundary element solution for the cauchy problem in linear elasticity using singular value decomposition. *Computer Methods in Applied Mechanics and Engineering*, 191(29-30):3257-3270, 2002.

[15] L. Marin and D. Lesnic. The method of fundamental solutions for the cauchy problem in two-dimensional linear elasticity. *International Journal of Solids and Structures*, 41(13):3425-3438, 2004.

[16] F. Ben Belgacem and H. El Fekih. On Cauchy's problem: I. a variational steklov-poincaré theory. *Inverse Problems*, 21(6):1915, 2005.

[17] M. Azaiez, F. Ben Belgacem, and H. El Fekih. On Cauchy's problem : Ii. completion, regularization and approximation. *Inverse Problems*, 22(4):1307-1336, 2006.

[18] D. L. Young, C. C. Tsai, C. W. Chen, and C. M. Fan. The method of fundamental solutions and condition number analysis for inverse problems of laplace equation. *Computers & Mathematics with Applications*, 55(6):1189-1200, 2008.

[19] Y. C. Hon and T. Wei. Backus-gilbert algorithm for the Cauchy problem of the laplace equation. *Inverse Problems*, 17(2):261, 2001.

[20] A. Cimetière, F. Delvare, M. Jaoua, and F. Pons. Solution of the Cauchy problem using iterated tikhonov regularization. *Inverse Problems*, 17(3):553-570, 2001.

[21] L. Bourgeois. A mixed formulation of quasi-reversibility to solve the Cauchy problem for laplace's equation. *Inverse Problems*, 21(3):1087-1104, 2005.

[22] L. Bourgeois. Convergence rates for the quasi-reversibility method to solve the cauchy problem for laplace's equation. *Inverse Problems*, 22(2):413, 2006.
[23] H. Egger and A. Leitão. Nonlinear regularization methods for ill-posed problems with piecewise constant or strongly varying solutions. *Inverse Problems*, 25(11):115014, 2009.

[24] P. Klüger and A. Leitão. Mean value iterations for nonlinear elliptic cauchy problems. *Numer. Math.*, 96(2):269-293, 2003.

[25] S. Andrieux and T. N. Baranger. An energy error-based method for the resolution of the Cauchy problem in 3D linear elasticity. *Computer Methods in Applied Mechanics and Engineering*, 197:902-920, 2008.

[26] T. N. Baranger and S. Andrieux. An optimization approach for the Cauchy problem in linear elasticity. *Structural and Multidisciplinary Optimization*, 35(2):141-152, 2008.

[27] T. N. Baranger and S. Andrieux. Data completion for linear symmetric operators as a Cauchy problem: an efficient method via energy like error minimization. *Vietnam Journal of Mechanics*, 31:247-261, 2009.

[28] S. Andrieux and T. N. Baranger. On the determination of missing boundary data for solids with nonlinear material behaviors, using displacement fields measured on a part of their boundaries. *Journal of the Mechanics and Physics of Solids*, 97:140-155, 2016.

[29] S. Andrieux and T. N. Baranger. Solution of nonlinear cauchy problem for hyperelastic solids. *Inverse Problems*, 31(11):115003, 2015.

[30] T. N. Baranger, S. Andrieux, and T. B. T. Dang. The incremental cauchy problem in elastoplasticity: General solution method and semi-analytic formulae for the pressurised hollow sphere. *Comptes Rendus Mécanique*, 343(5-6):331-343, 2015.

[31] S. Andrieux and T. N. Baranger. Three-dimensional recovery of stress intensity factors and energy release rates from surface full-field displacements. *International Journal of Solids and Structures*, 50(10):1523-1537, 2013.

[32] S. Andrieux and T. N. Baranger. Emerging crack front identification from tangential surface displacements. *Comptes Rendus Mécanique*, 340(8):565-574, 2012.

[33] T. N. Baranger and S. Andrieux. Constitutive law gap functionals for solving the Cauchy problem for linear elliptic PDE. *Applied Mathematics and Computation*, 218(5):1970-1989, 2011.

[34] R. Rischette, T. N. Baranger, and N. Debit. Numerical analysis of an energy-like minimization method to solve Cauchy problem with noisy data. *Journal of Computational and Applied Mathematics*, 235:3257-3269, 2011.

[35] R. Rischette, T. N. Baranger, and S. Andrieux. Regularization of the noisy cauchy problem solution approximated by an energy-like method. *International Journal for Numerical Methods in Engineering*, 95(4):271-287, 2013.

[36] R. Rischette, T. N. Baranger, and N. Debit. Numerical analysis of an energy-like minimization method to solve a parabolic cauchy problem with noisy data. *Journal of Computational and Applied Mathematics*, 271:206-222, 2014.

[37] T. N. Baranger, S. Andrieux, and R. Rischette. Combined energy method and regularization to solve the Cauchy problem for the heat equation. *Inverse Problems in Science and Engineering*, 0(0):1-14, 2013.

[38] T. N. Baranger, B. T. Johansson, and R. Rischette. On the alternating method for Cauchy problems and its finite element discretisation. In Larisa Beilina, editor, *Applied Inverse Problems*, volume 48 of *Springer Proceedings in Mathematics & Statistics*, pages 183-197. Springer New York, 2013.

[39] X. Escriva, T. C. N. Baranger, and N. H. Tlatli. Leak identification in porous media by solving the Cauchy problem. *Comptes Rendus Mécanique*, 335(7):401-406, JUL 2007.

[40] S. Andrieux, T. N. Baranger, and A. Ben Abda. Solving Cauchy problems by minimizing an energy-like functional. *Inverse Problems*, 22(1):115-133, 2006.

[41] S. Andrieux, A. Ben Abda, and T. N. Baranger. Data completion via an energy error functional. *Comptes Rendus Mécanique*, 333(2):171-177, 2005.

[42] M. Cocu. Existence of solutions of signorini problems with friction. *International Journal of Engineering Science*, 22(5):567-575, 1984.

[43] L. M. Bregman. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. {USSR} Computational
*Mathematics and Mathematical Physics*, 7(3):200 – 217, 1967.

[44] Y. Censor and A. Lent. An iterative row-action method for interval convex programming. *Journal of Optimization Theory and Applications*, 34(3), Jul 1981.

[45] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh. Clustering with bregman divergences. *J. Mach. Learn. Res.*, 6:1705–1749, December 2005.

[46] Béla A. Frigyik, Santosh Srivastava, and Maya R. Gupta. Functional bregman divergence and bayesian estimation of distributions. *Information Theory, IEEE Transactions on*, 54(11):5130–5139, November 2008.

[47] G. De Saxce and Z. Q. Feng. New inequality and functional for contact with friction: The implicit standard material approach. *Mechanics of Structures and Machines*, 19(3):301–325, 1991.

[48] P. Laborde and Y. Renard. Fixed point strategies for elastostatic frictional contact problems. *Mathematical Methods in the Applied Sciences*, 31(4):415–441, 2007.

[49] H. Boumediène Khenous, J. Pommier, and Y. Renard. Hybrid discretization of the signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers. *Applied Numerical Mathematics*, 56(2):163 – 192, 2006.

[50] EDF Electricité de France. Finite element *code_aster*, analysis of structures and thermomechanics for studies and research. Open source on www.code-aster.org, 1989–2018.

[51] E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001-2018.

[52] Mathwork. *Matlab version R2018a*, The Mathworks Inc, 2018.

[53] D. Colton and R. Kress. *Inverse Acoustic and Electromagnetic Scattering Theory*. Springer, 1998.

[54] H. D. Bui. *Fracture Mechanics: inverse problems and solutions*. Springer, 2006.

[55] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. *Physica D: Nonlinear Phenomena*, 60(1):259 – 268, 1992.

[56] F. Ben Belgacem. Why is the Cauchy problem severely ill-posed ? *Inverse Problems*, 23(2):823–836, 2007.