Tomáš Ligurský

Jaroslav Haslinger

Radek Kučera

Approximation and numerical realization of 3D contact problems with Coulomb friction and a solution-dependent coefficient of friction

Keywords: contact problems, Coulomb friction, solution-dependent coefficient of friction

HAL is

INTRODUCTION

Contact problems describing behaviour of a system of loaded deformable bodies in mutual contact have been of permanent interest in a few last decades. Their mathematical formulation has to include not only unilateral but also frictional conditions on the contact zones. Although in many situations the coefficient of friction F depends on the velocity of the tangential contact displacements (in quasistatic and dynamic problems), one usually assumes that F is a constant or it depends on the spatial variable. The present work concerns contact problems with 1 a solution-dependent coefficient of friction. Since we deal with a static case, the coefficient F will depend on the magnitude of the tangential contact displacements. In addition, we shall consider the Signorini type problem, i.e. the contact problem between an elastic body and a rigid foundation.

There are different approaches for numerical realization of contact problems with Coulomb friction which depend on the way how unilateral and friction conditions are treated. The overview and the comparison of the most frequently used methods can be found in [START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers[END_REF], e.g. Penalty and Lagrange multiplier techniques are the simplest ones. To overcome their drawbacks, augmented Lagrangian methods have been developed. The application of these methods in contact mechanics is described in [START_REF] Laursen | Computational and Impact Mechanics: Fundamentals of Modeling Interfaces Phenomena in Nonlinear Finite Element Analysis[END_REF]. Survey of algorithms of constrained optimization which are used in contact computational mechanics can be also found in [START_REF] Wriggers | Computational contact mechanics[END_REF]. Another efficient realization of contact problems uses Newton's type methods. The problem is written in the form of a system of non-smooth equations involving projection operators onto convex sets by means of which the unilateral and friction conditions are expressed. This system is solved by non-smooth variants of the Newton method (see [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF] for a general description of such methods in finite dimension and [START_REF] Hüber | A primal-dual active set algorithm for three dimensional contact problems with Coulomb friction[END_REF] for numerical realization of 3D contact problems with Coulomb friction). Here we present and apply the different approach which directly follows from the mathematical analysis in the theoretical part of the paper, and which is based on a fixed-point formulation solved by the method of successive approximations. Each iterative step leads to a contact problem with the Tresca model of friction in our case. There are several reasons why we decided just for this way of solving. First of all, it is very easy to implement as soon as one has the solver for Tresca friction at his disposal. Indeed, local Coulomb friction with a coefficient which eventually may depend on the solution can be easily integrated into the code simply by adding an extra iterative loop. Further, it turns out from our numerical experiments that the number of the fixed-point iterations is small (from 8 to 10) and practically it does not depend on the mesh size h. Thus the overall computational efficiency depends on how efficiently each iterative step can be realized. This is done by the newly developed algorithm of convex programming ( [START_REF] Dostál | An optimal algorithm for minimization of quadratic functions with bounded spectrum subject to separable convex inequality and linear equality constraints[END_REF]) which has the convergence rate determined by the condition number of the Hessian. With this result at hand one can prove the scalability of this algorithm when is used together with the TFETI domain decomposition technique. Last but not least, the algorithm treats the original quadratic constraints for friction forces so that no polygonal approximation of Coulomb friction law is necessary.

This paper extends results from [START_REF] Haslinger | Approximation and numerical realization of 3D contact problems with given friction and a coefficient of friction depending on the solution[END_REF], where a 3D Signorini problem with given friction and a coefficient depending on the norm of the tangential contact displacement, i.e. F := F( u t ), is studied. Instead of the model with given friction, the more realistic model with the local Coulomb law of friction will be considered now. On the other hand, only discrete problems are analysed and no convergence of discrete solutions is established.

The same problem has already been studied in [START_REF] Haslinger | Approximation and numerical realization of 2D contact problems with Coulomb friction and a solution-dependent coefficient of friction[END_REF] in 2D. However, its extension to 3D is not straightforward. The fact that the tangential component of displacements in 3D is a vector instead of a scalar in a plane case complicates not only the theoretical analysis but mainly numerical realization since one has to solve minimization problems involving both simple bounds and quadratic constraints instead of simple bounds solely. There is yet another difference between [START_REF] Haslinger | Approximation and numerical realization of 2D contact problems with Coulomb friction and a solution-dependent coefficient of friction[END_REF] and the present paper. Here we use both, the dualization of the unilateral conditions and the regularization of the frictional term in order to define a discrete form of our problem. On the other hand, in [START_REF] Haslinger | Approximation and numerical realization of 2D contact problems with Coulomb friction and a solution-dependent coefficient of friction[END_REF] only the dualization of the nonpenetration condition is used while keeping the original non-differentiable frictional term.

The paper is organized as follows: Section 2 introduces a classical and a weak formulation of the problem. Weak solutions are characterized as fixed-points of a mapping Φ which is defined by a contact problem with given friction and a coefficient of friction independent of the solution. In Section 3 a mixed variational formulation of this auxiliary problem is discretized by the finite element method: displacements, contact stresses are approximated by continuous piecewiselinear and piecewise-constant functions, respectively. This enables us to define an appropriate discretization Φ hH of the mapping Φ. Fixed-points of Φ hH represent discrete solutions of the contact problem with Coulomb friction and a solution-dependent coefficient of friction F. Existence of at least one discrete solution is shown for any continuous, non-negative and bounded function F. Its uniqueness is proved provided that F is small enough and Lipschitz continuous with a sufficiently small modulus of Lipschitz continuity. In addition, we derive the mesh-dependent bounds for F and the Lipschitz constant guaranteeing uniqueness of the solution. The method of successive approximations is proposed for numerical realization. Each iterative step leads to a contact problem with given friction and a fixed coefficient of friction. In Section 4 we introduce its algebraic dual formulation and recall an active-set algorithm for its solving. Finally, the efficiency of this numerical approach is tested on simple model examples.

Throughout the paper we shall denote the Euclidean norm of a vector by . . By H k (D), k ≥ 0 integer (H 0 (D) := L 2 (D)), we denote the standard Sobolev space of functions on D ⊂ R d , d = 2, 3, whose norm will be denoted by . k,D , in what follows. Further, the summation convention is adopted.

SETTING OF THE PROBLEM

Let us consider an elastic body which is represented by a bounded domain Ω ⊂ R 3 with the Lipschitz boundary ∂Ω. The boundary ∂Ω consists of three non-empty, non-overlapping parts Γ u , Γ p and Γ c : ∂Ω = Γ u ∪ Γ p ∪ Γ c . The body is fixed on Γ u , surface tractions of density p = (p 1 , p 2 , p 3 ) ∈ (L 2 (Γ p )) 3 act on Γ p while a rigid foundation S unilaterally supports the body along Γ c . For the sake of simplicity of our presentation we shall suppose that S is the half space

R 2 × R 1 -= {(x 1 , x 2 , x 3 ) ∈ R 3 | x 3 ≤ 0}
and there is no gap between S and Ω, i.e. Γ c is a part of the (x 1 , x 2 )-plane (see Figure 1). In addition, we shall take into account effects of friction on Γ c represented by the model with local Coulomb friction and a coefficient of friction F depending on a solution. Finally, the body is subject to volume forces of density f = (f 1 , f 2 , f 3 ) ∈ (L 2 (Ω)) 3 . We seek an equilibrium state of Ω.

The classical formulation of this problem consists in finding a displacement vector u = (u 1 , u 2 , u 3 ) which satisfies the equilibrium equations with the following boundary conditions:

(equilibrium equations) ∂τ ij ∂x j (u) + f i = 0 in Ω , i = 1, 2 , 3; (2.1)
(kinematical boundary conditions)

u i = 0 on Γ u , i = 1, 2, 3; (2.2)
(static boundary conditions)

T i (u) = p i on Γ p , i = 1, 2, 3; (2.3) p Γ p f Ω Γ c Γ u S x 1 p Figure 1. Geometry of the problem (unilateral conditions) u ν ≤ 0, T ν (u) ≤ 0, u ν T ν (u) = 0 on Γ c ; (2.4) (Coulomb's law of friction) u t = 0 =⇒ T t (u) ≤ -F(0)T ν (u) on Γ c ; u t = 0 =⇒ T t (u) = F( u t )T ν (u) u t u t on Γ c .    (2.5)
A symmetric stress tensor τ (u) = (τ ij (u)) 3 i,j=1 is related to a linearized strain tensor ε(u) = (ε ij (u)) 3 i,j=1 by means of linear Hooke's law:

τ ij (u) = c ijkl ε kl (u) , i, j = 1, 2, 3, ε kl (u) = 1 2 ∂u k ∂x l + ∂u l ∂x k , k, l = 1, 2, 3,
with linear elasticity coefficients c ijkl ∈ L ∞ (Ω), i, j, k, l = 1, 2, 3, satisfying the symmetry and ellipticity conditions:

c ijkl = c jikl = c klij a.e. in Ω, i, j, k, l = 1, 2, 3; ∃ c ell = const. > 0 : c ijkl ξ ij ξ kl ≥ c ell ξ ij ξ ij a.e. in Ω ∀ ξ ij = ξ ji ∈ R 1 . (2.6) Further, u ν = u i ν i , u t = u -u ν ν, where ν = (ν 1 , ν 2 , ν 3
) is the unit outward normal vector to ∂Ω, stand for the normal and tangential component of a displacement vector u on Γ c , respectively. The symbol

T (u) = (T 1 (u), T 2 (u), T 3 (u)), T i (u) = τ ij (u)ν j , i = 1, 2, 3, denotes a stress vector and T ν (u) = T i (u)ν i , T t (u) = T (u) -T ν (u)
ν is its normal and tangential component, respectively. Finally, the coefficient of friction F depending on the spatial variable x as well as the solution u will be given by a continuous, non-negative and bounded function:

F ∈ C Γ c × R 1 + , 0 ≤ F(x, ξ) ≤ F max ∀ (x, ξ) ∈ Γ c × R 1 + . (2.7) 4
T. LIGURSK Ý, J. HASLINGER, R. KU ČERA Let us notice that according to the special geometry of Ω and S we have:

u ν = -u 3 , u t = (u 1 , u 2 , 0) ; T ν (u) = -T 3 (u) , T t (u) = (T 1 (u), T 2 (u), 0) on Γ c .
Remark 2.1. Suppose that t 1 (x), t 2 (x), x ∈ Γ c , are two vectors such that the triplet {ν(x), t 1 (x), t 2 (x)} forms a local orthonormal basis in R 3 with the origin at x ∈ Γ c . Then any vector field v = (v 1 , v 2 , v 3 ) : Γ c → R 3 can be represented in the local coordinate system {ν(x), t 1 (x), t 2 (x)} as

v(x) = (v ν (x), v t (x)) ∈ R 1 × R 2 , x ∈ Γ c , with v ν (x) := v i (x)ν i (x), v t (x) := (v t1 (x), v t2 (x)), v tj (x) := v i (x)t j,i (x), j = 1, 2.
Hence, the tangential component of displacements on Γ c can be treated as a two-dimensional vector.

To present a weak formulation of our problem we introduce the following sets:

V = v ∈ H 1 (Ω) | v = 0 a.e. on Γ u , V = (V ) 3 , K = {v ∈ V | v ν ≤ 0 a.e. on Γ c } .
Further, let

H 1/2 (Γ c ) = ϕ ∈ L 2 (Γ c ) ∃ v ∈ V : v = ϕ a.e. on Γ c , H 1/2 (Γ c ) d = ϕ ∈ L 2 (Γ c ) d ∃ v ∈ (V ) d : v = ϕ a.e. on Γ c , d = 2, 3,
be the trace spaces on Γ c of all functions from V and (V ) d , d = 2, 3, equipped with the norms

ϕ 1/2,Γc = inf v∈V v=ϕ a.e. on Γc |v| 1,Ω , ϕ ∈ H 1/2 (Γ c ) , ϕ 1/2,Γc = inf v∈(V ) d v=ϕ a.e. on Γc |v| 1,Ω , ϕ ∈ H 1/2 (Γ c ) d , d = 2, 3,
respectively. The topological dual space of H 1/2 (Γ c ) will be denoted by H -1/2 (Γ c ) † , its norm and the corresponding duality pairing by . -1/2,Γc , ., . 1/2,Γc , respectively. Finally,

H 1/2 + (Γ c ) = ϕ ∈ H 1/2 (Γ c ) ϕ ≥ 0 a.e. on Γ c , H -1/2 + (Γ c ) = µ ∈ H -1/2 (Γ c ) µ, ϕ 1/2,Γc ≥ 0 ∀ ϕ ∈ H 1/2 + (Γ c )
will denote the non-negative cones in the respective spaces.

The weak formulation of (2.1)-(2.5) is given by the following implicit variational inequality of elliptic type:

Find u ∈ K such that a(u, v -u) -F( u t )T ν (u), v t -u t 1/2,Γc ≥ F (v -u) ∀ v ∈ K, (P)
where

a(v, w) := Ω τ ij (v)ε ij (w) dx , v, w ∈ V , F (v) := Ω f i v i dx + Γp p i v i dS , v ∈ V .
Remark 2.2. If u is smooth enough then applying Green's formula to (P) we recover (2.1)-(2.5). Let us notice, however, that to make sense to the duality term in (P) one needs an additional smoothness of u and F ensuring that [START_REF] Eck | Unilateral contact problems. Variational methods and existence theorems[END_REF]). To overcome this difficulty we shall suppose that

F( u t ) v t ∈ H 1/2 (Γ c ) ∀ v ∈ V (see
T ν (u) ∈ L 2 (Γ c ) in what follows.
If it is so then the duality pairing ., . 1/2,Γc can be replaced by the L 2 (Γ c )-scalar product and (2.7) is sufficient.

Due to (2.6) and Korn's inequality, a is a symmetric bilinear form which is (H 1 (Ω)) 3 -elliptic and continuous on V :

∃ α > 0 : a(v, v) ≥ α v 2 1,Ω ∀ v ∈ V , (2.8)
∃ M > 0 : a(v, w) ≤ M v 1,Ω w 1,Ω ∀ v, w ∈ V .
(2.9)

In the sequel we shall also need the following auxiliary results.

Lemma 2.1. It holds:

(i) if ϕ ∈ (H 1 (Γ c )) 2 then its Euclidean norm ϕ ∈ H 1 (Γ c ) and ϕ 1,Γc ≤ ϕ 1,Γc ; (ii) if ϕ ∈ (H 1/2 (Γ c )) 3 then ϕ t ∈ (H 1/2 (Γ c )) 2 and ϕ t 1/2,Γc ≤ ϕ 1/2,Γc .
For the proofs we refer to [START_REF] Ligurský | Approximation and numerical realization of 3D contact problems with given friction and a coefficient of friction depending on the solution[END_REF].

Below we introduce a fixed-point formulation of (P) which will be used for the discretization of our problem. With any (ϕ, g) ∈ H

1/2 + (Γ c ) × L 2 + (Γ c
) we associate the following auxiliary problem:

Find u := u(ϕ, g) ∈ K such that a(u, v -u) + (F(ϕ)g, v t -u t ) 0,Γc ≥ F (v -u) ∀ v ∈ K, (P(ϕ, g))
where (., .) 0,Γc stands for the scalar product in L 2 (Γ c ). Problem (P(ϕ, g)) is a weak formulation of a contact problem with given friction and the fixed coefficient F(ϕ). It has a unique solution for any (ϕ, g) ∈ H

1/2 + (Γ c ) × L 2 + (Γ c
) as follows from [START_REF] Glowinski | Numerical methods for nonlinear variational problems[END_REF]. This enables us to define the mapping Φ :

H 1/2 + (Γ c ) × L 2 + (Γ c ) → H 1/2 + (Γ c ) × H -1/2 + (Γ c ) by Φ(ϕ, g) = ( u t , -T ν (u)) , (ϕ, g) ∈ H 1/2 + (Γ c ) × L 2 + (Γ c ) , 6
T. LIGURSK Ý, J. HASLINGER, R. KU ČERA

Comparing the definitions (P) and (P(ϕ, g)) it is readily seen that if the couple

( u t , -T ν (u)) is a fixed-point of Φ in H 1/2 + (Γ c ) × L 2 + (Γ c ), i.e.: Φ( u t , -T ν (u)) = ( u t , -T ν (u)) ,
then u is a solution to (P).

Let (ϕ, g) ∈ H

1/2 + (Γ c ) × L 2 + (Γ c ) be arbitrary but fixed. Denote J(ϕ, g, v) = 1 2 a(v, v) + j(ϕ, g, v) -F (v) , v ∈ V ,
the total potential energy functional, where

j(ϕ, g, v) = (F(ϕ)g, v t ) 0,Γc , v ∈ V . (2.10)
It is well-known that (P(ϕ, g)) is equivalent to the following constrained minimization problem:

Find u := u(ϕ, g) ∈ K such that J(ϕ, g, u) ≤ J(ϕ, g, v) ∀ v ∈ K.
Further, let

Λ ν = H -1/2 + (Γ c ) , Λ t (ϕ, g) = µ t ∈ L 2 (Γ c ) 2 µ t ≤ F(ϕ)g a.e. on Γ c
be the Lagrange multiplier sets. It can be proved that

v ∈ K ⇐⇒ v ∈ V & µ ν , v ν 1/2,Γc ≤ 0 ∀ µ ν ∈ Λ ν (2.11) and j(ϕ, g, v) = sup µ t ∈Λt(ϕ,g) (µ t , v t ) 0,Γc . (2.12)
Therefore a solution u of (P(ϕ, g)) can be also characterized as follows:

J(ϕ, g, u) = min v∈K J(ϕ, g, v) = min v∈V sup µν ∈Λν µ t ∈Λt(ϕ,g) L(v, µ ν , µ t ) , where L(v, µ ν , µ t ) = 1 2 a(v, v)-F (v)+ µ ν , v ν 1/2,Γc +(µ t , v t ) 0,Γc , (v, µ ν , µ t ) ∈ V ×Λ ν ×Λ t (ϕ, g) ,
is the Lagrangian.

By a mixed variational formulation of (P(ϕ, g)), (ϕ, g) ∈ H

1/2 + (Γ c ) × L 2 + (Γ c ), we call a problem of finding a saddle-point of L on V × Λ ν × Λ t (ϕ, g) or equivalently: Find (u, λ ν , λ t ) ∈ V × Λ ν × Λ t (ϕ, g) such that a(u, v) = F (v) -λ ν , v ν 1/2,Γc -(λ t , v t ) 0,Γc ∀ v ∈ V , µ ν -λ ν , u ν 1/2,Γc + (µ t -λ t , u t ) 0,Γc ≤ 0 ∀ (µ ν , µ t ) ∈ Λ ν × Λ t (ϕ, g) .      (M(ϕ, g))
It is known (see [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF]) that there exists a unique solution (u, λ ν , λ t ) of (M(ϕ, g)) for any (ϕ, g) ∈ H

1/2 + (Γ c ) × L 2 + (Γ c ).
In addition, u solves (P(ϕ, g)), λ ν = -T ν (u) and λ t = -T t (u). From this we obtain the alternative form of the mapping Φ, namely:

Φ(ϕ, g) = ( u t , λ ν ) , (ϕ, g) ∈ H 1/2 + (Γ c ) × L 2 + (Γ c ) , (2.13)
where u := u(ϕ, g), λ ν := λ ν (ϕ, g) is the first and the second component of the solution to (M(ϕ, g)), respectively.

FINITE ELEMENT DISCRETIZATION

This section deals with a discretization of problem (P). It will be based on the fixedpoint formulation for appropriate discretizations of the mapping Φ and the sets of Lagrange multipliers. The existence as well as the uniqueness of the solutions will be established.

For the sake of simplicity of our presentation we shall suppose that Ω is a polyhedron. Let T h be a partition of Ω into tetrahedra whose diameters do not exceed h which is compatible with the decomposition of ∂Ω into Γ u , Γ p and Γ c . With T h the following sets of continuous, piecewise-linear functions will be associated:

V h = v h ∈ C Ω v h| T ∈ P 1 (T ) ∀ T ∈ T h , v h = 0 on Γ u , V h = (V h ) 3 , W h = ϕ h ∈ C Γ c ∃ v h ∈ V h : v h = ϕ h on Γ c , W + h = {ϕ h ∈ W h | ϕ h ≥ 0 on Γ c } . Obviously, V h ⊂ V and W + h ⊂ H 1/2 + (Γ c
). Further, let T H be a partition of Γ c into rectangles whose diameters do not exceed H. With T H we associate the space L H of piecewise-constant functions:

L H = µ H ∈ L 2 (Γ c ) µ H |R ∈ P 0 (R) ∀ R ∈ T H
and the following discretizations of the sets of Lagrange multipliers:

Λ νH = {µ νH ∈ L H | µ νH ≥ 0 on Γ c } , Λ tH (ϕ h , g H ) = µ tH ∈ (L H ) 2 µ tH | R ≤ 1 meas 2 (R) R F(ϕ h )g H dS ∀ R ∈ T H ,
where meas 2 (R) is the area of R and (ϕ h , g H ) ∈ W + h × Λ νH . For any (ϕ h , g H ) ∈ W + h ×Λ νH we define the mixed finite element approximation of a contact problem with given friction g H and a coefficient F(ϕ h ) as follows:

Find (u h , λ νH , λ tH ) ∈ V h × Λ νH × Λ tH (ϕ h , g H ) such that a(u h , v h ) = F (v h ) -(λ νH , v hν ) 0,Γc -(λ tH , v ht ) 0,Γc ∀ v h ∈ V h , (µ νH -λ νH , u hν ) 0,Γc + (µ tH -λ tH , u ht ) 0,Γc ≤ 0 ∀ (µ νH , µ tH ) ∈ Λ νH × Λ tH (ϕ h , g H ) .          (M hH (ϕ h , g H ))
Let (ϕ h , g H ) ∈ W + h × Λ νH be given. We shall give the interpretation of the first component of the solution (u h , λ νH , λ tH ) to (M hH (ϕ h , g H )). Denote by

K hH = v h ∈ V h R v hν dS ≤ 0 ∀ R ∈ T H , j H (ϕ h , g H , v h ) = sup µ tH ∈ΛtH (ϕ h ,gH ) (µ tH , v ht ) 0,Γc , v h ∈ V h ,
the approximations of K and j, respectively (cf. (2.11) and (2.12)). From

(µ νH -λ νH , u hν ) 0,Γc ≤ 0 ∀ µ νH ∈ Λ νH (3.1)
it easily follows that u h ∈ K hH . Further,

(µ tH , u ht ) 0,Γc ≤ (λ tH , u ht ) 0,Γc ∀ µ tH ∈ Λ tH (ϕ h , g H ) so that j H (ϕ h , g H , u h ) = sup µ tH ∈ΛtH (ϕ h ,gH ) (µ tH , u ht ) 0,Γc = (λ tH , u ht ) 0,Γc . (3.2) Substitution v h := v h -u h , v h ∈ K hH , into (M hH (ϕ h , g H )) 2 leads to a(u h , v h -u h ) = F (v h -u h ) -(λ νH , v hν -u hν ) 0,Γc -(λ tH , v ht -u ht ) 0,Γc ≥ F (v h -u h ) -j H (ϕ h , g H , v h ) + j H (ϕ h , g H , u h )
in virtue of (3.1), (3.2) and the definitions of K hH , Λ νH and j H , respectively. Therefore u h solves the following variational inequality of the second kind:

Find u h := u h (ϕ h , g H ) ∈ K hH such that a(u h , v h -u h ) + j H (ϕ h , g H , v h ) -j H (ϕ h , g H , u h ) ≥ F (v h -u h ) ∀ v h ∈ K hH .      (P hH (ϕ h , g H ))
Remark 3.1. The set K hH is an external approximation of K because the non-penetration condition v hν ≤ 0 on Γ c is satisfied in a weak (integral) sense. Moreover, it holds that

j H (ϕ h , g H , v h ) = (F(ϕ h )g H , π H v ht ) 0,Γc , v h ∈ V h , (3.3)
(cf. (2.10)), where π H : (L 2 (Γ c )) 2 → (L H ) 2 is the (L 2 (Γ c )) 2 -orthogonal projection onto (L H ) 2 : (π H ϕ, µ H ) 0,Γc = (ϕ, µ H ) 0,Γc ∀ µ H ∈ (L H ) 2 , ϕ ∈ L 2 (Γ c ) 2 .
Indeed, let v h ∈ V h be arbitrarily chosen. On the one hand one has:

j H (ϕ h , g H , v h ) = sup µ tH ∈ΛtH (ϕ h ,gH ) R∈TH (µ tH , v ht ) 0,R = sup µ tH ∈ΛtH (ϕ h ,gH ) R∈TH µ tH,i |R R v hti dS ≤ sup µ tH ∈ΛtH (ϕ h ,gH ) R∈TH µ tH |R R v ht dS ≤ R∈TH 1 meas 2 (R) R F(ϕ h )g H dS R v ht dS
using the definition of Λ tH (ϕ h , g H ). On the other hand, taking µ tH ∈ Λ tH (ϕ h , g H ) such that

µ tH |R := 1 meas2(R) R F(ϕ h )g H dS 1 R v ht dS R v ht dS if R v ht dS = 0 ; 0 otherwise
one can verify the opposite inequality. Thus

j H (ϕ h , g H , v h ) = R∈T H R F(ϕ h )g H 1 meas 2 (R) R v ht dS dS = R∈TH R F(ϕ h )g H π H v ht dS = (F(ϕ h )g H , π H v ht ) 0,Γc .
In what follows we shall assume that the following stability condition is satisfied for any

(ϕ h , g H ) ∈ W + h × Λ νH : If (µ νH , µ tH ) ∈ Λ νH × Λ tH (ϕ h , g H ) is such that (µ νH , v hν ) 0,Γc + (µ tH , v ht ) 0,Γc = 0 ∀ v h ∈ V h then (µ νH , µ tH ) = (0, 0).      (3.4)
This condition guarantees (see [START_REF] Hlaváček | Solution of variational inequalities in mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF]) that (M hH (ϕ h , g H )) has a unique solution for any (ϕ h , g H ) ∈ W + h ×Λ νH . It is fulfilled e.g. when the partition T H is coarser than the triangulation T h| Γc (for more details we refer to [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF]).

Let r h : H 1/2 (Γ c ) → W h be a linear interpolation operator possessing the monotonicity property:

ϕ

∈ H 1/2 (Γ c ), ϕ ≥ 0 a.e. on Γ c =⇒ r h ϕ ∈ W + h (3.5)
and the following approximation property:

∃ c r > 0 : ϕ -r h ϕ 0,Γc ≤ c r h Γc ϕ 1,Γc ∀ ϕ ∈ H 1 (Γ c ) , (3.6)
where the constant c r is independent of h Γc := max F ∈T h| Γc diam(F ) provided that T h| Γc belongs to a family T h| Γc , h → 0+, of regular triangulations of Γ c . As an example of such r h we refer to [START_REF] Clement | Approximation by finite element functions using local regularization[END_REF]. With such r h at hand we define the mapping Φ hH :

W + h × Λ νH → W + h × Λ νH by Φ hH (ϕ h , g H ) = (r h u ht , λ νH ) , (ϕ h , g H ) ∈ W + h × Λ νH
, where u h := u h (ϕ h , g H ) and λ νH := λ νH (ϕ h , g H ) are the first two components of the solution to (M hH (ϕ h , g H )). Taking into account that u h , λ νH are approximations of u, λ ν , respectively, the mapping Φ hH can be viewed to be a discretization of Φ (cf. (2.13)).

Definition 3.1. The triplet (u h , λ νH , λ tH ) is called a solution of a discrete contact problem with Coulomb friction and a solution-dependent coefficient of friction if it solves (M hH (r h u ht , λ νH )), i.e. (r h u ht , λ νH ) is a fixed-point of Φ hH in W + h × Λ νH : Φ hH (r h u ht , λ νH ) = (r h u ht , λ νH ) .

Next we prove the existence of at least one fixed-point of Φ hH using the Brouwer fixedpoint theorem. To this end, the spaces L H , (L H ) 2 will be endowed with the following norms, respectively:

µ H -1/2,h = sup v h ∈V h v h =0 (µ H , v h ) 0,Γc v h 1,Ω , µ H ∈ L H , (3.7
)

µ H -1/2,h = sup v h ∈(V h ) 2 v h =0 (µ H , v h ) 0,Γc v h 1,Ω , µ H ∈ (L H ) 2 . (3.8)
The norm in W h × L H will be introduced as follows:

(ϕ h , g H ) W h ×LH := ϕ h 0,Γc + µ H -1/2,h , (ϕ h , g H ) ∈ W h × L H .
Let us notice that (3.7) and (3.8) are mesh-dependent dual norms, making use of (3.4).

Lemma 3.1. Assume that F satisfies (2.7). Then there exist constants R 1 , R 2 > 0 such that

Φ hH maps W + h × Λ νH into (W + h × Λ νH ) ∩ B, where B = {(ϕ h , g H ) ∈ W h × L H | ϕ h 0,Γc ≤ R 1 & µ H -1/2,h ≤ R 2 } .
In addition, if T h| Γc belongs to a system T h| Γc , h → 0+, of strongly regular triangulations of Γ c then both constants R 1 , R 2 are independent of h Γc .

Proof. Let (ϕ h , g H ) ∈ W + h × Λ νH be fixed. Inserting v h := 0 ∈ K hH into (P hH (ϕ h , g H )) we obtain: α u h 2 1,Ω ≤ a(u h , u h ) + j H (ϕ h , g H , u h ) ≤ F (u h ) ≤ F * u h 1,Ω (3.9)
according to (2.8) and the non-negativeness of the functional j H . Here . * stands for the norm in the dual space to (H 1 (Ω)) 3 . From (3.6), (i ) of Lemma 2.1 and the inverse inequality between the (L 2 (Γ c )) 2 and (H 1 (Γ c )) 2 -norms for functions from (W h ) 2 it follows that

r h u ht 0,Γc ≤ r h u ht -u ht 0,Γc + u ht 0,Γc ≤ c r h Γc u ht 1,Γc + u ht 0,Γc ≤ c r h Γc u ht 1,Γc + u ht 0,Γc ≤ (c inv c r + 1) u ht 0,Γc . (3.10)
Further,

u ht 0,Γc ≤ u h 0,Γc ≤ c t u h 1,Ω , (3.11)
where c t is the norm of the trace mapping from V into (L 2 (Γ c )) 3 . This together with (3.9) and (3.10) give:

r h u ht 0,Γc ≤ R 1 := c t (c inv c r + 1) α F * .
From (M hH (ϕ h , g H )) 2 and (2.9) we obtain:

(λ νH , v hν ) 0,Γc + (λ tH , v ht ) 0,Γc = F (v h ) -a(u h , v h ) ≤ F * v h 1,Ω + M u h 1,Ω v h 1,Ω ∀ v h ∈ V .
This and (3.9) imply:

sup v h ∈V h v h =0 (λ νH , v hν ) 0,Γc + (λ tH , v ht ) 0,Γc v h 1,Ω ≤ F * + M u h 1,Ω ≤ 1 + M α F * . Consequently λ νH -1/2,h = sup v h ∈V h v h =0 (λ νH , v h ) 0,Γc v h 1,Ω ≤ sup v h ∈V h v h =0 (λ νH , v hν ) 0,Γc + (λ tH , v ht ) 0,Γc v h 1,Ω ≤ R 2 := 1 + M α F * , (3.12)
taking into account the special geometry of Ω.

If the system T h| Γc , h → 0+, is strongly regular then the constant c inv does not depend on h Γc (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]) and neither does c r . The other constants are clearly independent of h Γc anyway.

⊓ ⊔ Lemma 3.2. The mapping Φ hH is continuous in W + h × Λ νH provided that (2.7) is satisfied. Proof. Let ϕ (k) h , g (k) H → (ϕ h , g H ) in W h × L H , k → ∞ ,
where ϕ

(k) h , g (k)
H , (ϕ h , g H ) ∈ W + h × Λ νH , k ∈ N, and let u

(k) h , λ (k) νH , λ (k) tH ∈ V h × Λ νH × Λ tH (ϕ (k) h , g (k)
H ) be the respective solutions to M hH ϕ

(k) h , g (k) H
, k ∈ N. From (3.9), (3.12) one can see that the sequences u

(k) h , λ (k) νH
are bounded in V h , L H , respectively. In a similar way as in (3.12) it can be shown that λ

(k) tH is bounded in (L H ) 2 . Since V h , L H , (L H ) 2 are finite-dimensional spaces, there exist: a subsequence u (l) h , λ (l) νH , λ (l) tH ⊆ u (k) h , λ (k) νH , λ (k) tH and a triplet (u h , λ νH , λ tH ) ∈ V h × Λ νH × Λ tH (ϕ h , g H ) such that u (l) h → u h in V h , λ (l) νH → λ νH in L H , λ (l) tH → λ tH in (L H ) 2 , l → ∞ .
We shall prove that (u h , λ νH , λ tH ) solves (M hH (ϕ h , g H )). Let (v h , µ νH , µ tH ) ∈ V h × Λ νH × Λ tH (ϕ h , g H ) be arbitrary but fixed. Then one can find a sequence µ

(l) tH , µ (l) tH ∈ Λ tH ϕ (l) h , g (l) H ∀ l, such that µ (l) tH → µ tH in (L H ) 2 , l → ∞ . Inserting (v h , µ νH , µ tH ) := v h , µ νH , µ (l) tH ∈ V h ×Λ νH ×Λ tH ϕ (l) h , g (l) H into M hH ϕ (l) h , g (l) H
and letting l → ∞ there we get:

a(u h , v h ) = F (v h ) -(λ νH , v hν ) 0,Γc -(λ tH , v ht ) 0,Γc , (µ νH -λ νH , u hν ) 0,Γc + (µ tH -λ tH , u ht ) 0,Γc ≤ 0 , i.e. (u h , λ νH , λ tH ) solves (M hH (ϕ h , g H ))
. Since this problem possesses a unique solution, the original sequences u 2 , respectively, and from the trace theorem it follows that

(k) h , λ (k) νH , λ (k) tH tend to u h , λ νH , λ tH in V h , L H , (L H )
u (k) ht → u ht in L 2 (Γ c ) 2 , k → ∞ . (3.13)
By the monotonicity property (3.5) of r h it is easy to verify that

r h u (k) ht -u ht ≤ r h u (k) ht -u ht on Γ c , ∀ k ∈ N . Therefore r h u (k) ht -r h u ht 0,Γc ≤ r h u (k) ht -u ht 0,Γc ≤ c r h Γc u (k) ht -u ht 1,Γc + u (k) ht -u ht 0,Γc ≤ (c inv c r + 1) u (k) ht -u ht 0,Γc , (3.14)
arguing as in (3.10). This combined with (3.13) yield:

r h u (k) ht → r h u ht in W h , k → ∞ .
⊓ ⊔ We arrived at the following result. Theorem 3.1. Let (2.7) be satisfied. Then there exists at least one solution to the discrete contact problem with Coulomb friction and a solution-dependent coefficient of friction.

To obtain the uniqueness of the solution we shall need the stronger assumption on F, namely:

∃ c L > 0 : F(x, ξ) -F x, ξ ≤ c L ξ -ξ ∀ x ∈ Γ c and all ξ, ξ ∈ R 1 + . ( 3
Φ hH (ϕ h , g H ) -Φ hH ( φh , ḡH ) W h ×LH ≤ max{C 1 F max , C 2 c L } (ϕ h -φh , g H -ḡH ) W h ×LH holds for any (ϕ h , g H ), ( φh , ḡH ) ∈ (W + h × Λ νH ) ∩ B,
where B is the same as in Lemma 3.1. Proof. Let (ϕ h , g H ), ( φh , ḡH ) ∈ (W + h × Λ νH ) ∩ B be given and let (u h , λ νH , λ tH ), ūh , λνH , λtH be the solutions of (M hH (ϕ h , g H )), (M hH ( φh , ḡH )), respectively. Then u h , ūh solves (P hH (ϕ h , g H )) and (P hH ( φh , ḡH )), respectively, which making use of (3.3) become:

a(u h , v h -u h ) + (F(ϕ h )g H , π H v ht -π H u ht ) 0,Γc ≥ F (v h -u h ) ∀ v h ∈ K hH , a(ū h , v h -ūh ) + (F( φh )ḡ H , π H v ht -π H ūht ) 0,Γc ≥ F (v h -ūh ) ∀ v h ∈ K hH .
Substituting v h := ūh ∈ K hH into the first and v h := u h ∈ K hH into the second inequality and summing both inequalities we obtain:

α u h -ūh 2 1,Ω ≤ a(u h -ūh , u h -ūh ) ≤ (F(ϕ h )g H -F( φh )ḡ H , π H ūht -π H u ht ) 0,Γc ≤ F(ϕ h )g H -F( φh )ḡ H 0,Γc π H ūht -π H u ht 0,Γc . (3.16)
Using the approximation property of π H , the inverse inequality between the (H 1/2 (Γ c )) 2 and (H 1 (Γ c )) 2 -norms for functions from (W h ) 2 , (ii ) of Lemma 2.1 and the trace theorem we have:

π H u ht -π H ūht 0,Γc ≤ π H u ht -π H ūht 0,Γc ≤ π H (u ht -ūht ) -(u ht -ūht ) 0,Γc + u ht -ūht 0,Γc ≤ c π H u ht -ūht 1,Γc + u ht -ūht 0,Γc ≤ c ′ inv c π H h Γc u ht -ūht 1/2,Γc + u ht -ūht 0,Γc ≤ c ′ inv c π H h Γc u h -ūh 1/2,Γc + u h -ūh 0,Γc ≤ c ′ inv c π H h Γc + c t u h -ūh 1,Ω .
From this and (3.16) we see that

u h -ūh 1,Ω ≤ 1 α c ′ inv c π H h Γc + c t F(ϕ h )g H -F( φh )ḡ H 0,Γc . (3.17)
Further, (2.7) and (3.15) yield:

F(ϕ h )g H -F( φh )ḡ H 0,Γc ≤ F(ϕ h )(g H -ḡH ) 0,Γc + (F(ϕ h ) -F( φh ))ḡ H 0,Γc ≤ F max g H -ḡH 0,Γc + c L ϕ h -φh ∞,Γc ḡH 0,Γc .
Since all norms in finite-dimensional spaces W h , L H are equivalent, we have:

∃ c ′′ inv > 0 : ϕ h ∞,Γc ≤ c ′′ inv h Γc ϕ h 0,Γc ∀ ϕ h ∈ W h , (3.18)
∃ c ′′′ inv > 0 : µ H 0,Γc ≤ c ′′′ inv √ H µ H -1/2,h ∀ µ H ∈ L H so that F(ϕ h )g H -F( φh )ḡ H 0,Γc ≤ c ′′′ inv √ H F max g H -ḡH -1/2,h + c ′′ inv R 2 h Γc c L ϕ h -φh 0,Γc , (3.19)
using that ḡH -1/2,h ≤ R 2 . Let us mention that the constants c ′′ inv , c ′′′ inv may depend on the mesh sizes h and H, respectively. In the same way as in (3.14) one can prove that

r h u ht -r h ūht 0,Γc ≤ (c inv c r + 1) u ht -ūht 0,Γc ≤ c t (c inv c r + 1) u h -ūh 1,Ω . (3.20)
From (3.17), (3.19) and (3.20) we obtain:

r h u ht -r h ūht 0,Γc ≤ c ′′′ inv c t α √ H (c inv c r + 1) c ′ inv c π H h Γc + c t • • F max g H -ḡH -1/2,h + c ′′ inv R 2 h Γc c L ϕ h -φh 0,Γc . (3.21)
Subtracting (M hH ( φh , ḡH )) 2 from (M hH (ϕ h , g H )) 2 one has: λνH -λ νH , v hν 0,Γc + λtH -λ tH , v ht 0,Γc = a(u h -ūh , v h ) ≤ M u h -ūh 1,Ω v h 1,Ω ∀ v h ∈ V h .
Arguing as in (3.12) one can show that

λ νH -λνH -1/2,h ≤ sup v h ∈V h v h =0 λνH -λ νH , v hν 0,Γc + λtH -λ tH , v ht 0,Γc v h 1,Ω ≤ M u h -ūh 1,Ω .
This combined with (3.17) and (3.19) lead to

λ νH -λνH -1/2,h ≤ c ′′′ inv M α √ H c ′ inv c π H h Γc + c t • • F max g H -ḡH -1/2,h + c ′′ inv R 2 h Γc c L ϕ h -φh 0,Γc . (3.22)
Summing (3.21) and (3.22), the assertion of the theorem follows with h × Λ νH ) ∩ B and the discrete contact problem with Coulomb friction and a solution-dependent coefficient of friction has a unique solution. In addition, this solution is the limit of the sequence generated by the method of successive approximations:

C 1 := c ′′′ inv α √ H c ′ inv c π H h Γc + c t [c t (c inv c r + 1) + M ] , (3.23)
C 2 := c ′′ inv c ′′′ inv R 2 αh Γc √ H c ′ inv c π H h Γc + c t [c t (c inv c r + 1) + M ] . ( 3
let ϕ (0) h , g (0) H ∈ W + h × Λ νH be given; for k = 0, 1, . . . , ϕ (k) h , g (k) H ∈ W + h × Λ νH known, set: ϕ (k+1) h , g (k+1) H := Φ hH ϕ (k) h , g (k) H ;        for any choice of ϕ (0) h , g (0) H ∈ W + h × Λ νH .
Suppose that the Babuška-Brezzi condition is satisfied, i.e. there exists a positive constant β independent of h Γc , H such that sup

v h ∈V h v h =0 (µ νH , v hν ) 0,Γc + (µ tH , v ht ) 0,Γc v h 1,Ω ≥ β[ µ νH -1/2,Γc + µ tH -1/2,Γc ] (3.25) holds for any (µ νH , µ tH ) ∈ L H × (L H ) 2 . Setting µ tH := 0 ∈ (L H ) 2 in (3.25
) and taking into account the special geometry of Ω one can easily verify that

β µ νH -1/2,Γc ≤ µ νH -1/2,h ≤ µ νH -1/2,Γc ∀ µ νH ∈ L H .
Thus the mesh-dependent norm . Remark 3.2. Let (3.25), (3.26) be satisfied and let T h| Γc , h → 0+, and {T H }, H → 0+, be such that the ratio H h Γc is bounded from above for h, H → 0+. The sufficient condition guaranteeing the contractivity of Φ hH and hence the uniqueness of the solution is that the parameters F max , c L decay as √ H, h Γc √ H, h, H → 0+, respectively. In two-dimensional problems, the condition on the decay of c L is weaker, namely c L ∼ h Γc H. This is due to the fact that the inverse inequality (3.18) holds with 1 h Γc in 2D. If the coefficient F does not depend on u, i.e. c L ≡ 0, then we arrive at the result from [START_REF] Haslinger | Approximation of the Signorini problem with friction, obeying the Coulomb law[END_REF].

NUMERICAL EXPERIMENTS

Our computations are based on the method of successive approximations mentioned in Corollary 3.1. In order to simplify notation we will use the same symbols for algebraic counterparts of continuous functions from the previous text.

Let us note that the projected gradient enables us to write the optimality criterion characterizing the solution λ ∈ R m + × Λ t (ϕ, g) of (4.2) in the form r(λ) = 0. Our algorithm is based on the fact that the non-zero components of r(µ) at µ = λ determine the descent directions changing appropriately the active-set. To this end, we introduce components of r(µ) called the projected free gradient φ := φ(µ) and the projected boundary gradient β := β(µ), respectively, as follows:

φ A = 0, φ M\A = r M\A , β A = r A , β M\A = 0.
We combine three steps to generate a sequence {µ (l) } that approximates the solution to (4.2):

• the expansion step: µ (l+1) = µ (l) -α φ(µ (l) ),

• the proportioning step: µ (l+1) = µ (l) -α β(µ (l) ),

• the conjugate gradient step: µ (l+1) = µ (l) -α (l) cg p (l) , where α (l) cg and the conjugate gradient directions p (l) are computed recurrently [START_REF] Golub | Matrix computations[END_REF]; the recurrence starts from µ (s) generated by the last expansion or the proportioning step and satisfies A(µ (l+1) ) = A(µ (s) ).

The expansion step may add indices while the proportioning step may release indices to/from the current active-set. The conjugate gradient steps are used to carry out efficiently the minimization of the objective f in the interior of the face W (µ

(s) ) = {µ ∈ R m + × Λ t (ϕ, g)| µ A := µ (s)
A , A = A(µ (s) )}. Moreover, the algorithm exploits a given constant γ > 0 and the releasing criterion

β(µ (l) ) ⊤ r(µ (l) ) ≤ γ φ(µ (l) ) ⊤ r(µ (l) ) (4.3)
to decide which of the steps will be performed.

Algorithm 4.2 Let µ (0) ∈ R m + ×Λ t (ϕ, g), γ > 0, α ∈ (0, S -1 ] and ǫ λ ≥ 0 be given. For µ (l) , µ (s) known, 0 ≤ s ≤ l, where µ (s) is computed by the last step expansion or proportioning, choose µ (l+1) by the following rules:

(i) If r(µ (l) ) ≤ ǫ λ , return λ := µ (l) . (ii) If µ (l) fulfils (4.3), try to generate µ (l+1) by the conjugate gradient step. If µ (l+1) ∈ Int W (µ (s) ), accept it, otherwise generate µ (l+1) by the expansion step. (iii) If µ (l) does not fulfil (4.3), generate µ (l+1) by the proportioning step.

Contrary to simple bound problems analyzed in [START_REF] Dostál | Minimizing quadratic functions subject to bound constraints with the rate of convergence and finite termination[END_REF], the algorithm does not exhibit the finite terminating property while the same convergence rate is achieved. The upper bounds defining the set Λ tH (ϕ h , g H ) are computed by the quadrature formula which is exact for continuous piecewise-linear functions over T h| R i (cf. Figure 4

f (µ (l+1) ) -f (λ) ≤ η f (µ (l) ) -f (λ) , where η = 1 - αα min 2 + 2 γ < 1. T. LIGURSK Ý, J. HASLINGER, R. KU ČERA p S Γ 1 p Γ 2 p Γ c p Ω Γ u x 1 x 3 x 2 Γ 3 p Figure 2.
R 1 R 2 R 3 R i-1 R i R i+1 x 1 x 2 (a) Partition T H of Γc a i1 a i2 a i3 a i4 a i5 a i6 a i7 a i8 a i9 R i (b) Numbering of nodes
(b)): Ri F(ϕ h )g H dS ≈ 1 3 8F(ϕ h (a i1 )) + 2 9 j=2 F ϕ h a ij g H |R i 1 8 meas 2 (R i )
for any R i ∈ T H . Consequently, F i (ϕ) in (4.1) is defined as

F i (ϕ) = 1 12 4F( (ϕ 2i1-1 , ϕ 2i1 ) ) + 9 j=2 F ϕ 2ij -1 , ϕ 2ij , i = 1, . . . , m.
The stability condition (3.4) is obviously satisfied. The initial approximation and the terminating tolerance for the method of successive approximations were chosen to be ϕ (0) = g (0) = 0 and ε = 1e-5, respectively. Both the fixed precision control with r tol = 0.01 and the adaptive precision control with r tol = 0.01 and c f act = 0.9 of the inner loop were tested.

Table I shows how our numerical method behaves for different meshes with F given by param = 6e4, while Table II compares computations for different F on the mesh constructed with n div = 16. Recall that 3n, 3m stands for the total number of the primal and the dual variables, respectively. Further iter denotes the total number of the fixed-point iterations and n A stands for the number of solving the linear system Ax = b, where A is the stiffness matrix. Since this step is the most expensive part of Algorithm 4.2, n A expresses the total cost of computations. The first integer in the respective columns characterizes the fixed precision control while the second integer characterizes the adaptive one. Figure 5 depicts the convergence history of the method of successive approximations (i.e. the dependence of err (k) on k) for the mesh with n div = 16 and the fixed precision control.

The figures illustrate the typical behaviour of the resulting displacement u h and of the Lagrange multipliers λ νH , λ tH (param = 6e4, n div = 16). The deformed body (the deformation is 500× enlarged) is shown in Figure 6. Figures 7 and8 depict the behaviour of -u hν and λ νH on Γ c while Figures 9 and 10 present the graphs of the Euclidean norm u ht and λ tH on Γ c , respectively. The distribution of F( u ht ), of the upper bound F i (u ht )(λ νH ) i , i = 1, . . . , m, on Γ c is shown in Figures 11 and12, respectively. Finally, Figure 13 visualizes in more details λ tH . The radii of the circles situated in the individual elements R i ∈ T H correspond to the bounds F i (u ht )(λ νH ) i while the segments emanating from the centres of these circles represent the components of the vector λ tH in R i .

CONCLUSIONS

The paper is devoted to the discretization and numerical realization of 3D contact problems with Coulomb friction and a coefficient of friction F depending on the solution. The discrete solutions are defined as fixed-points of the mapping Φ hH . Its definition is based on a finite element discretization of a mixed variational formulation of contact problems with given friction and a coefficient of friction which does not depend on the solution. We proved the existence of at least one fixed-point for any F given by a continuous, non-negative and bounded function. Conditions guaranteeing its uniqueness were established. The method of successive approximations was proposed for finding the fixed-points. Each iterative step corresponds to a quadratic programming problem with simple bounds and separable quadratic constraints which was solved by an active-set algorithm. Two strategies for determining the precision control were proposed. The computations guided by the fixed precision control indicate that only a small number of the fixed-point iterations suffices to get a solution with the required accuracy. Moreover, the number of the iterations depends only slightly on both the mesh size and the Lipschitz constant of F. On the other hand the strategy based on the adaptive precision control turned out to be more efficient. In fact, the computational cost of the respective computations, which is characterized by n A , needed for finding the solution with the same accuracy was lower. Although no preconditioning was used in our computations, the values on n A did not increase dramatically for finer meshes. The finite element spaces used in our examples were built on structured meshes. But the proposed numerical approach turned out to be efficient also in contact problems for several elastic bodies

  John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 0:0-0 Prepared using nmeauth.cls Corollary 3.1. If F max , c L from (2.7) and (3.15), respectively, are small enough then Φ hH is contractive in (W +

  -1/2,h can be replaced by the dual norm . -1/2,Γc in the previous estimates. Furthermore assume that (i ) T h| Γc belongs to a system T h| Γc , h → 0 + , of strongly regular triangulations of Γ c . (ii ) T H belongs to a family {T H }, H → 0 + , of strongly regular partitions of Γ c . c inv , c ′ inv , c ′′ inv , c ′′′ inv , c r and c π in (3.23) and (3.24) do not depend on h Γc , H (in particular, for the inverse inequality (3.18) see [1, Th. 3.2.6]).

Theorem 4 . 1 .

 41 Let λ ∈ R m + × Λ t (ϕ, g) be the solution to (4.2), α min denote the smallest eigenvalue of S and γ = max{γ, γ -1 }. Let {µ (l) } be the sequence generated by Algorithm 4.2 with ǫ λ = 0. Then

  Figure 2. Geometry

Figure 3 .

 3 Figure 3. Coefficients of friction

Figure 4 .

 4 Figure 4.

Figure 10 .Figure 11 .Figure 12 . 2 Figure 13 .

 101112213 Figure 5. Convergence history

  .15) Theorem 3.2. Let F satisfy (2.7) and (3.15). Then there exist positive numbers C 1 , C 2 such that

Table I .

 I Different meshes

	n div	3n	3m	iter	n A	Table II. Different coefficients F
	4 6	900 2646	36 81	8/16 2195/ 994 9/18 5765/ 840	param	iter	n A
	8	5832 144 10/21 8864/2215	2e4	9/37 11255/4714
	10 10890 225 10/20 7969/2015	6e4	11/22	9870/1995
	12 18252 324 10/21 9171/2385	3e5	13/43 15379/5608
	14 28350 441 10/21 9686/2113
	16 41616 576 11/22 9870/1995

Copyright c 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 0:0-0 Prepared using nmeauth.cls

22

T. LIGURSK Ý, J. HASLINGER, R. KU ČERA discretized by unstructured meshes.

Contract/grant sponsor: Grant Agency of the Czech Republic; contract/grant number: 201/07/0294 (all authors) Contract/grant sponsor: Ministry of Education, Youth and Sports of the Czech Republic; contract/grant number: MSM0021620839 (the first two authors) Contract/grant sponsor: Charles University Grant Agency; contract/grant number: 18008 (the first author)

Denote n = dim V h , m = dim L H and p the number of the contact nodes of T h . With any ϕ ∈ R 2p and g ∈ R m + we associate the set

The way how to define the upper bounds in (4.1), in particular F i (ϕ), will be shown in the forthcoming parts of this section. The algebraic form of (M hH (ϕ h , g H )) reads as follows:

where (•, •) q stands for the inner product in R q , A ∈ R 3n×3n is the symmetric, positive definite stiffness matrix, F ∈ R 3n is the load vector and N ∈ R m×3n , T ∈ R 2m×3n are the full rowrank matrices coupling u with the dual variables λ ν , λ t . Eliminating u, we obtain the so-called reciprocal variational formulation of the problem:

where

and S := BA -1 B ⊤ , h := BA -1 F , B := (N ⊤ , T ⊤ ) T . We arrive at the following implementation of the method of successive approximations:

ν and go to step (i).

The overall efficiency of our numerical approach depends on a method used in step (i). As (4.2) is a strictly convex problem with the quadratic objective f subject to separable constraints (the simple bounds and the quadratic inequalities), we can solve it by an active-set strategy proposed in [START_REF] Kučera | Minimizing quadratic functions with separable quadratic constraints[END_REF] and analyzed in [START_REF] Kučera | Convergence rate of an optimization algorithm for minimizing quadratic functions with separable convex constraints[END_REF]. Below we shall give a brief description of the algorithm.

Let M = {1, . . . , 3m}. The gradient of f at a point µ ∈ R m + ×Λ t (ϕ, g) is r := r(µ) = Sµ-h and the active-set A ⊆ M at µ := (µ ν , µ t ) is defined by

Using the projection P : R 3m → R m + × Λ t (ϕ, g), we define the projected gradient for a fixed α > 0 as The error in the S-energy norm is bounded by

Proof. See [START_REF] Kučera | Convergence rate of an optimization algorithm for minimizing quadratic functions with separable convex constraints[END_REF]. ⊓ ⊔ Theorem 4.1 yields the optimal value of η for γ = γ = 1 and α = S -1 , when

where κ(S) is the spectral condition number of S.

Remark 4.1. To increase the efficiency of Algorithm 4.1, we initialize Algorithm 4.2 in the kth step by the result obtained in the previous step, i.e., µ (0

t ). Moreover, we choose the inner precisions ǫ λ := ǫ (k) λ sufficiently accurate in order to achieve the terminating tolerance ε for the method of successive approximations. We use two strategies: (a) the fixed precision control ǫ ) h with 0 < r tol < 1, 0 < c fact < 1, err (-1) = 1 and ǫ (-1) λ = r tol /c fact (typically r tol = 0.01 and c fact = 0.9). While the strategy (a) makes it possible to obtain the solution in a small number of outer iterations, the strategy (b) leads to considerably more efficient procedure with a small number of matrixvector multiplications.

In the numerical experiments we consider an elastic, isotropic and homogeneous material characterized by Young's modulus E = 21.19e10 [Pa] and Poisson's ratio σ = 0.277 (steel). The initial configuration is Ω = (0, 3) × (0, 1) × (0, 1) (in m) with Γ u = {0} × (0, 1) × (0, 1), Γ c = (0, 3) × (0, 1) × {0} and Γ p = Γ 1 p ∪ Γ 2 p ∪ Γ 3 p , where Γ 1 p = {3} × (0, 1) × (0, 1), Γ 2 p = (0, 3) × (0, 1) × {1} and Γ 3 p = (0, 3) × {0, 1} × (0, 1). The density of surface tractions is prescribed as follows:

x , 0, p 1 z ) on Γ 1 p , p = (0, 0, p 2 z ) on Γ 2 p , p = (0, 0, 0) on Γ 3 p , where p 1 x = 1e7 [Pa], p 1 z = 2e7 [Pa] and p 2 z = -3e7 [Pa] (see Figure 2). The volume forces are neglected, i.e. f = 0 in Ω. The function F representing the coefficient of friction is independent of the spatial variable, i.e. F(x, ξ) := F(ξ), where

Three different values of param were considered, namely param = 2e4, 6e4 and 3e5 (cf.

Figure 3). The partition T h is constructed in two steps: Firstly, Ω is cut into 3n div × n div × n div cubes, n div even. Secondly, each of these cubes is divided into five tetrahedra. To the partition T h we associate the dual partition T H as shown in Figure 4(a): The fine lines and the black dots represent the triangulation T h| Γc and its nodes, respectively, while the "chessboard" with panes R i comprising eight triangles of T h| Γc constitutes the dual partition T H .


Conditions de vol

La réalisation de 4 vols successifs de 12 minutes chacun a permis de couvrir une surface d'environ 30 ha à une hauteur constante de 40m par rapport au point de décollage et à une vitesse moyenne de 5m/s, seul le vol YS-20161116-085822 a été réalisé à 35 m de hauteur avec un espacement de 15 m entre chaque bande de vol et à une vitesse de 4m/s.

Traitement des nuages de points

Correction GPS et assemblage

Extraction des points sols

Les points sols ont été extraits en utilisant le plugin CSF 2

Pour la zone Dorgnon-Breuil : Cette zone a été traitée en deux parties. Le processus de recalage a été réalisé dans le logiciel CloudCompare grâce aux outils d'alignement de deux nuages de point par recherche d'équivalence de points. Nous avons ainsi chercher des points homologues sur les lignes de vol 1 et 2, recaler les deux nuages entre eux. Puis, rajouté la ligne de vol 3. Chercher de nouveaux points homologues. La progression c'est faite ainsi en rajoutant les lignes de vol au fur et à mesure.

Il est très difficile de calculer la marge d'erreur de ce type d'alignement manuel. En effet, a chaque recalage d'une ligne de vol, il y a ajout d'une erreur. Chaque erreur de recalage n'est pas strictement cumulable. Il faudra faire une évaluation de cette marge d'erreur, même de façon approximative.

Génération des MNT ombrés

Les points sols ont été importés dans le logiciel ArcGis.

Figure 2 -

 2 Figure 2-Tusson (16) -Mission du 16/11/16 -Localisation des vols sur fond d'image satellite (Google)

Figure 3 -Figure 6 -

 36 Figure 3-Tusson (16) -Zones survolées -Mission du 16/11/2016 -Nuage de points total (~450M) -Vue vers le sud

Figure 7 -Figure

 7 Figure 7-Tusson (16) -Zones survolées -Mission du 16/11/2016 -Nuage de points sol (~ 450M). Vue vers le S

Figure 11 -

 11 Figure 11-Tusson (16)-Zone Dognon-Breuil -du 16/11/2016 -Coupe est/ouest permettant de distinguer les décalages entre les lignes de vol.

Figure 12 -

 12 Figure 12-Tusson (16)-Zone Dognon-Breuil -Mission du 16/11/2016 -Coupe nord/sud permettant de distinguer les décalages entre les lignes de vol.

1-

  Création d'une géodatabase 2-Création d'un jeu de classe d'entité 3-Conversion des .las en Multipoint avec les outils de l'extension 3D Analyst 4-Création d'un jeu "Terrain" dans la géodatabase 5-Conversion du multipoint en MNT avec les outils de l'extension 3D Analyst 6-Conversion du MNT en raster avec les outils de l'extension 3D Analyst 7-Création d'un ombrage

FigureFigure 15 -

 15 Figure 13-Tusson (16)-Zone Dognon-Breuil -Mission du 16/11/2016 -MNT ombré

Figure 16 -

 16 Figure 16-Tusson (16) -Zone Dognon-Breuil -Mission du 16/11/2016 -MNT ombré (Zoom)

  Les 4 vols ont été corrigés des imprécisions GPS selon la procédure PPP (Precise Point Positioning) offerte par le site Natural Resources Canada 1

		puis assemblés dans CloudCompare (v2.8).
	Les nuages de points comportes 28.593.346 points avant filtrage des points sols, répartit de la façon
	suivante :
	Zone Dognon-Breuil	=> 17.554.124 points
	Zone Justice	=> 3.772.880 points
	Zone Le Magnou	=> 7.266.341 points

https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php?locale=en

Pour la zone Le Magnou, le relevé Lidar par drone a permis la délimitation du site. Les contours de ce dolmen apparaissent. Il est de forme légèrement circulaire ; une dépression marque son centre.


Carine Calastrenc

Nicolas Poirier

François Baleux

Lagarde (65) -Castet-Crabe

teaching and research institutions in France or abroad, or from public or private research centers.

Objectif et contexte

La mission d'acquisition de données Lidar par drone du 06/12/2016 a été réalisée sur la commune de Lagarde (65), sur l'emprise d'une forêt.

Conditions de vol

La réalisation d'un vol de 12 minutes a permis de couvrir une surface d'environ 14 ha à une hauteur constante de 80m par rapport au point de décollage (situé à 30 mètres en contrebas de la forêt) et à une vitesse moyenne de 5m/s.

Traitement des nuages de points

Correction GPS et assemblage

La presque totalité de la surface du site (hormis la pointe Sud et la face Est) a été couverte pas un vol. La dernière ligne de vol qui devait documenter la limite Est du site n'a pu être réalisé en raison d'une coupure de communication entre la station sol et le drone.

Les données issues de ce vol ont été corrigés des imprécisions GPS selon la procédure PPP (Precise Point Positioning) offerte par le site Natural Resources Canada 1 puis assemblés dans CloudCompare (v2.8).

Le nuage de points global comporte 8 332 248 points, soit environ 64 pts/m² avant filtrage des points sols.

Extraction des points sols

Les points sols ont été extraits en utilisant le plugin Lasgroud 2 implémenté dans QGis avec le paramètre wilderness. Les points sols obtenus sont au nombre de 1 282 691, soit environ 10 pts/m².

1 https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php?locale=en

Génération des MNT ombrés

Les points sols ont été importés dans le logiciel Saga Gis 2.1.2.

La génération du MNT ombré a suivi la procédure décrite par V. Wichmann 3 à savoir :

Interprétation

Une première interprétation des images produites permet d'identifier un site fortifié de hauteur. Il s'agit d'un site a double système de fossé. Le premier, le plus large, délimite un espace d'environ 98 hectares. Ce fossé mesure 24 m de largeur en moyenne et a une profondeur maximale de 10 m. Le fossé intérieur encercle un espace de 2,5 hectares. Il se développe sur 10 m de large et 6 mètres de hauteur. Dans la partie Sud, on observe une déformation du sol qui divise l'espace en deux "terrasses". Dans l'angle nord-est de l'enceinte intérieur on remarque un linéament. Ce premier traitement de l'ensemble du nuage de point laisse voir, dans la partie centrale du site, des anomalies qui pourraient être identifiées comme de possibles structures archéologiques.

Les points documentation l'intérieur de ce site ont été isolés avec Lasground (utilisation de la même procédure et même critère que pour la nuage de point complet). Une fois la classification faites, le MNT a été réalisé avec ArcGis.

  L'objectif était d'évaluer le potentiel de la technologie Lidar embarquée sur drone pour la détection de sites fortifiés de hauteur. Le contexte topographique et végétal est celui d'une petite forêt située dans une zone assez vallonnée couverte d'une végétation mixte composée d'arbres, arbustes et de milieu ouvert (prairie).

Figure 1 -

 1 Figure 1-Lagarde (65) -Mission du 06/12/2016 -Localisation du site

Figure 2 -

 2 Figure 2-Lagarde (65) -Mission du 06/12/2016 -Nuage de points total (~30M) -Vue vers le Nord

Figure 3 -

 3 Figure 3-Lagarde (65) -Mission du 06/12/2016 -Nuage de points sol (~ 15M). Vue vers le Nord

Figure 4 -

 4 Figure 4-Lagarde (65) -Mission du 06/12/2016 --MNT ombré

Figure 5 -

 5 Figure 5-Lagarde -mission du 06/12/2016 -Modèle finalisé

Figure 7 -

 7 Figure 7-Lagarde (62) --mission du 06/12/2016 -Interprétation des résultats dans la zone centrale du site

Figure 8 -

 8 Figure 8-Lagarde -mission du 06/12/2016 -Zoom intérieur du site après traitement


Mei Zhu

Cun Li

Xiufeng Wang

Licun Li

Jean-Pascal Sutter

Thermal Magnetic Hysteresis in a Copper-Gadolinium-Radical Chain Compound

Magnetic bistability spanning over a temperature domain of 40 K can result from a small structural deformation of the gadolinium aminoxyl coordination. This is illustrated for a nitronyl nitroxide 3d-4f chain, [Ln(hfac) 3 Cu(hfac) 2 (NIT-Pyrim) 2 ] (Ln III = Gd, Dy), which is the first example of a bistable lanthanide-based complex.

Temperature-triggered bistable molecular materials have received great attention because of their potential applications in switching, sensing, and information storage. 1 The most remarkable examples of such materials are the iron(II) spin-crossover (SCO) compounds taking advantage of conversion of the electronic configuration between high and low spins. 2 Related magnetic bistability has also been observed in valence tautomeric complexes, 3 exchange-coupled systems such as organic radicals, 4 and radical-transition metal complexes, 5 in which a structural phase transition leads to significant modifications of the exchange interactions. Spectacular examples include a dithiadiazolyl radical (i.e., 1,3,5-trithia-2,4,6-triazapentalenyl) found to exhibit bistability over a large temperature domain around room temperature as a result of a dimerizationdissociation process between radical units 4a and copper nitroxide complexes for which the bistability is associated with an axial-toequatorial coordination isomerism for the radical ligand. 5a The temperature domain for which the bistability exists and the width of thermal hysteresis are two characteristics of special importance in these materials. One of the approaches followed to achieve large hysteresis for the bistability domains is to strengthen the cooperation, i.e., the elastic interactions, between magnetic molecules. A classical way is to increase the intermolecular interactions by hydrogen-bonding, π-π, etc., 6 interactions or to extend the dimensions of the molecular system. 7 It has been reported recently that magnetic bistability with a rather wide thermal hysteresis loop may also result from a structural phase transition associated with a disorder-order phenomenon for a molecular fragment. 8 To date, SCO-like behaviors relying on a disruption of the metal-radical exchange interaction mainly concern copper(II) derivatives. 5a,9 Herein we report a case of spin-transition-like behavior with wide hysteresis resulting from variation of the lanthanide-radical coordination within 2p-3d-4f mixed-spin compounds. While the amplitude of the magnetic changes is less spectacular than that for the copper-radical complexes, the magnetic transitions are sharp and the domain of bistability spans over about 40 K.

The coordination polymers [Ln(hfac) 3 Cu(hfac) 2 (NIT-Pyrim) 2 ][ L n III =G d( 1), Dy (2); NIT-Pyrim = 2-(5pyrimidinyl)-4,4,5,5-tetramethylimidazoline-1-oxyl-3-oxide] were obtained by reacting Ln(hfac) 3 (Ln III = Gd, Dy) with NIT-Pyrim followed by Cu(hfac) 2 in a 1:2:1 ratio. The crystal structures for 1 and 2 (Figure 1) were determined from singlecrystal X-ray diffraction data collected first at room temperature (HT) and subsequently at 113 K (LT), that is, above and below the magnetic anomaly (vide infra). Because the two compounds are isomorphous, only the structure of 1 will be briefly described. At 293 K, 1 crystallizes in the monoclinic P2/c space group. It consists of a one-dimensional (1D) array of Gd(hfac) 3 and Cu(hfac) 2 units bridged by NIT-Pyrim linked to the metal centers respectively through one aminoxyl moiety and the pyrimidine group (Figure 1). The copper(II) ion is surrounded by two nitrogen atoms from two pyrimidine groups coordinated in a cis configuration and by four oxygen atoms from two hfac anions. The average equatorial Cu-O/N bond distance is 2.17 Å, and the apical Cu-O bond length is 1.942( 7) Å. Accordingly, the coordination polyhedron of the copper(II) ion may be represented as a compressed tetragonal geometry. The gadolinium(III) ion is eight-coordinated with two oxygen atoms from two nitronyl nitroxide moieties and six oxygen atoms of three hfac anions. Continuous-shape-measure analysis indicates that the gadolinium(III) ion is located in a distorted dodecahedral environment (D 2d ; Table S4). 10 The Gd-O rad distance with 2.389(8) Å compares well with those of the reported Ln(hfac) 3 -nitronyl nitroxide complexes; 11 the Gd-O-N-C(sp 2 ) torsion angle is 91(1)°. The intrachain Gd•••Cu distance is 8.202 Å, and the shortest interchain metal-metal separation is found between gadolinium and copper with a distance of 10.945 Å. The shortest distance between two aminoxyl moieties belonging to different chains is 4.336 Å (O5••• O5 -x,1-y,-z ). Finally, it can be pointed out that the fluorine atoms for two CF 3 groups in the surroundings of gadolinium are disordered over two positions (Figure 1).

For the structure collected at 113 K, the space group remains unchanged, and no important modifications are found for the cell parameters (Table S1). These are limited to a shortening of the a (the chain growth direction) and c parameters by about 0.5 Å (i.e., 3.5%) and 0.1 Å (0.5%) respectively, while b is increased by 0.1 Å (0.9%). The overall molecular structure is identical with that found at 293 K with similar bond parameters with the noticeable exception of the Gd-O-N-C(sp 2 ) torsion angle, which is reduced from 91(1)°to 86(1)°(Figure 1). Moreover, the CF 3 groups in the Gd(hfac) 3 moieties are no longer found to be disordered. Hence, the structural modifications between the HT and LT phases mainly come down to a disordered-toordered location of the CF 3 groups, accompanied by an alteration of the Gd-O-N-C(sp 2 ) torsion angles. Interestingly, these rather minor structural alterations have a significant effect on the magnetic feature of the complex.

The χ M T versus T behavior for 1 (Figure 2; χ M stands for the molar magnetic susceptibility for a formula unit) is characterized by marked discontinuities that are observed in both the cooling and warming modes (ramp, 1 K min -1 ). Starting from a value of 9.48 cm 3 K mol -1 (the anticipated spin-only value is 9.0 cm 3 K mol -1 ) at 300 K, χ M T ↓ ( ↓ indicates data collected while cooling) gradually increases upon cooling. However, for about 176 K, a slight but clear drop of χ M T ↓ is observed, followed by a more pronounced decrease below 131 K (from 9.69 to 9.47 cm 3 K mol -1 for 115 K). Upon further cooling, χ M T ↓ smoothly increases again, reaching a maximum of 9.87 cm 3 K mol -1 for 16 K before falling steeply for lower T. Upon heating, the sample χ M T ↑ superposes to the behavior observed during cooling until 55 K; afterward, it runs below (Figure 2). However, between 155 and 172 K, a sudden increase is found, with χ M T ↑ coming close to that observed in the cooling mode. A second but small hop at 198 K leads to merging of the two curves for the higher temperature domain. The same behavior is found for a subsequent coolingheating cycle (Figure 2, green plot) apart from the disappearance of the small drop at 176 K seen initially during cooling. It can be noticed that the related small step found while heating is shifted to 185 K. Such an open-S-shaped behavior reveals thermal hysteresis where the magnetic response of 1 depends on its thermal history. Considering only the domain for which the difference in χ M T for the two states remains the largest, the bistability temperature range spans from 123 to 163 K, that is, a width of 40 K. The possibility of freezing the HT phase by fast cooling from 298 to 10 K was also considered; however, the subsequent χ M T ↑ behavior upon warming was characteristic of the LT phase. The field dependence of magnetization for 1 recorded at 2.0 K (Figure S6) shows a rapid increase of M for low fields and reaches 9.82 Nβ at 70 kOe, a behavior in agreement with ferromagnetic interactions among the spin carriers. In addition, alternating-current (ac) magnetic susceptibility measurements confirmed the absence of magnetic ordering above 2 K (Figure S7).

The steps observed in the χ M T versus T behavior are typical for modification of the exchange interactions taking place in 1.On the basis of the HT and LT crystal structures, the chains are well separated and interchain interactions may be excluded to account for the observed behavior. A structural modification more likely to alter an exchange interaction is the change in the gadoliniumaminoxyl coordination. Indeed, experimental observations and theoretical calculation have shown that the strength of the gadolinium-nitroxide interaction depends on the Gd-O-N-C torsion angle: the larger the Gd-O-N-C torsion angle, the stronger should be the Gd-ON magnetic coupling. 12 Analysis of the magnetic behavior for 1 confirms this trend. Bearing in mind that the copper(II) ion is located in a compressed octahedral environment, its magnetic orbital is d z 2 ; hence, the possible magnetic interaction with the nitronyl nitroxide radical through the pyrimidine unit (located in an equatorial position) should be very weak. Thus, from a magnetic point of view, a chain consists of independent Rad-Gd-Rad units plus one noninteracting copper(II) ion. For the Rad-Gd-Rad unit, two kinds of exchange interactions are operative, namely, the direct Gd-ON interaction (J 1 ) and a next-neighbor interaction between two NO groups via the gadolinium(III) ion (J 2 ). The theoretical expression for χ M T derived from the Hamiltonian H = -2J 1 (S Rad1 S Gd1 + S Gd1 S Rad2 ) -2J 2 S rad1 S rad2 (see the Supporting Information) was used to model the behaviors between 300 and 136 K and between 14 and 140 K respectively for the HT and LT phases. 13 The best fit to the experimental behaviors yielded J 1 = 3.52 cm -1 , J 2 = -7.14 cm -1 , g = 2.02, and g Cu = 2.10 for the HT behavior and J 1 = 1.00 cm -1 , J 2 = -6.36 cm -1 , g = 2.02, and g Cu = 2.10 for the LT domain, values consistent with those reported in similar compounds. 14 These J 1 values support the role of the Gd-O-N-C torsion angle in the major transition of the magnetic behavior found in 1. The different χ M T profiles observed for the HT and LT phases can clearly be ascribed to the marked difference in the strength of the ferromagnetic gadolinium-radical interaction (3.5 versus 1.0 cm -1 ); the jumps of χ M T result from the temperature-triggered structural transitions, modifying this exchange parameter. Differential scanning calorimetry (DSC) confirmed the coincidence of the magnetic steps and the phase transition (Figure S5). Interestingly, like the two-step transition found in the magnetic behavior, DSC also shows two signals, which suggests that the structural modification involves two steps close in temperatures. While the LT and HT structures are informative on the phase before and after the structural transition, its stepping remains unknown.

As mentioned above, the width of the bistability domain is correlated with the cooperative effect between the molecules. In the absence of obvious chemical interactions between the chains, the observed wide hysteresis in 1 might be associated with the disorder-to-order transformation of CF 3 on the gadolinium moieties. At room temperature, two CF 3 per lanthanide are disordered, while they are ordered at low temperature. This transformation couples with the alteration of the gadoliniumradical coordination, which results in thermal magnetic bistability. It is tempting to suggest that ordering-disordering of the CF 3 groups and structural distortion take place consecutively and contribute respectively to the small and larger steps seen in the χ M T behaviors.

For the dysprosium derivative 2, the χ M T versus T behaviors in the cooling and warming modes are almost the same (Figure S8) and do not exhibit the steps found for 1 despite a related structural behavior (the Dy-O-N-C torsion angles are 92(1)°a nd 87(1)°at 293 and 113 K, respectively). The anticipated change in χ M T with the phase transition might be less marked because of a smaller variation of the ferromagnetic dysprosiumradical interaction and/or be overwhelmed by the contributions of the crystal-field effect of dysprosium(III). 15 The absence of slow magnetization relaxation for 2 was checked by ac susceptibility measurements (Figure S10).

Interestingly, a rather modest change in the Gd-O-N-C torsion angle is sufficient to signific a n t l ym o d i f yt h e ferromagnetic gadolinium-radical interaction and, in turn, result in marked differences in the magnetic behavior of the HT and LT phases. Moreover, the dimensionality of the system and the disorder-order transformation of CF 3 groups may play crucial roles in the wide thermal hysteresis loop characterizing the compound. Exchange-coupled metal-radical compounds exhibiting spin-transition-like behavior with hysteresis are very scarce, 9,16 and to the best of our knowledge, the reported chain is the first example of a lanthanide-based compound with such a desirable feature. This finding may suggest new opportunities to construct bistable systems.

Figure 1 .

 1 Figure 1. Crystal structure of complex 1 (some fluorine and hydrogen atoms are hidden for clarity): (top) view of the 1D coordination polymer at 293 K developing along the a axis showing the structurally disordered CF 3 groups; (bottom) detail of the arrangement of the gadolinium surroundings at 293 and 113 K.

Figure 2 .

 2 Figure 2. χ M T versus T plots for 1 in the cooling and warming modes for two cycles. The red lines represent the calculated behaviors.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (Grant 21471083) and MOE Innovation Team (IRT13022) of China.

AUTHOR INFORMATION

Corresponding Authors E-mail: llicun@nankai.edu.cn. E-mail: sutter@lcc-toulouse.fr.

Notes

The authors declare no competing financial interest.