1 Introduction
###### Abstract

The quasistatic normal-compliance contact problem of isotropic homogeneous linear visco-elastic bodies with Coulomb friction at small strains in Kelvin-Voigt rheology is considered. The discretization is made by a semi-implicit formula in time and the Symmetric Galerkin Boundary Element Method (SGBEM) in space, assuming that the ratio of the viscosity and elasticity moduli is a given relaxation-time coefficient. The obtained recursive minimization problem, formulated only on the contact boundary, has a nonsmooth cost function. If the normal compliance responds linearly and the 2D problems are considered, then the cost function is piecewise-quadratic, which after a certain transformation gets the quadratic programming (QP) structure. However, it would lead to second-order cone programming in 3D problems. Finally, several computational tests are presented and analysed, with additional discussion on numerical stability and convergence of the involved approximated Poincaré-Steklov operators.

Keywords: contact mechanics, tribology, evolution variational inequalities, numerical approximation, boundary element method (BEM), mathematical programming, computational simulations.

AMS Classification: 35Q90, 49N10, 65K15, 65M38, 74A55, 74S15, 90C20.

Quasistatic normal-compliance contact problem of

visco-elastic bodies with Coulomb friction

implemented by QP and SGBEM.

Roman Vodička111Technical University of Košice, Civil Engineering Faculty, Vysokoškolská 4, 042 00 Košice, Slovakia, roman.vodicka@tuke.sk, Vladislav Mantič222University of Seville, School of Engineering, Camino de los Descubrimientos s/n, 41092 Seville, Spain, mantic@etsi.us.es, Tomáš Roubíček333Charles University, Mathematical Institute, Sokolovská 83, CZ-186 75 Praha 8, Czech Republic, roubicek@karlin.mff.cuni.cz444Institute of Thermomechanics of the Czech Acad. Sci., Dolejškova 5, CZ–182 08 Praha 8, Czech Rep.

## 1 Introduction

Mathematically justified modelling and efficient numerical solution of unilateral contact problems with friction is very challenging. It is now well recognized that the combination of the Coulomb friction with unilateral Signorini contact is very difficult and essentially there are no satisfactory mathematical results (i.e. not only for small coefficient of friction or small data or short time) available in the literature. On the other hand, in engineering a so-called normal compliance is a well-accepted compromise (and sometimes considered even as a more realistic variant) to the Signorini problem. This compliance concept is based on an idea of a rough contact surface which is intimately related with the concept of friction resulting microscopically from a certain asperity of the surface and which, when pressed by a normal force, exhibits certain elastic response, macroscopically demonstrated as a certain interpenetration, cf. [25, Sect. 11.4.2], [34, Sect. 4.10], [58, Chap. 5] or [52]. This normal-compliance concept is also well amenable for existence analysis, cf. [3, 26, 43], and for analysis of convergence of numerical approximation, although its algorithmic implementation might be more difficult because (in contrast to the Signorini problem) does not lead directly to a quadratic-programming or a 2nd-order-cone-programming problem.

When nonlinearities are restricted to the boundaries, small-strain concept is adopted, and the material in the bulk is linear and homogeneous (so that the fundamental solutions of the underlying operators are at disposal), the Boundary Element Method (BEM) can be an efficient and very accurate alternative to the Finite Element Method (FEM) usually employed in engineering computations, e.g. [58, 31, 25, 29, 13, 21, 22, 27]. In this paper, we confine ourselves to the mentioned situation in the bulk, so that all nonlinearities arising from the normal-compliance model and the Coulomb friction will be indeed only on the boundary or interfaces between the bodies. It should be noted that Coulomb friction can be reduced to the solution of so-called Tresca friction, see [12, 21] for which error estimates are available, see [17]. Also, it is worthy of mention that we neglect inertia.

There exist several quite successful approaches for the solution of contact problems by BEM, e.g. [5, 15, 8, 16, 33, 38] and references therein. In our present approach we develop a new model, which uses normal compliance contact with visco-elastic behaviour of the bulk but implemented so that the elasto-static boundary element (BE) approach is used in combination with the quadratic programming (QP).

The paper is organized as follows. After formulation of the problem in Section 2, we made the following steps towards devising an approximation and an efficient algorithm in Section 3:

• Discretization made by a semi-implicit formula in time (Sect. 3.1).

• Implementation of BEM, assuming the ratio of the viscosity and elasticity moduli is a given relaxation-time coefficient (Sect. 3.2).

• The changes in energy formulation related to implementation by the Symmetric Galerkin Boundary Element Method (SGBEM) (Sect. 3.3).

• A transformation of the nonsmooth incremental problems to get the quadratic programming (QP) structure in 2D contact problems with linearly-responding normal compliance (Sect. 3.4).

Let us emphasize that the method is robust in the sense that it avoids any nonconvex minimization, and even the minimizers at each particular time level can be found in a finite number of steps, and has a guaranteed numerical stability (a-priori estimates) for arbitrary coefficient of friction and large data and long times, as shown in Appendix C. Most importantly, involving a viscous rheology and the normal compliance ensures that the continuous friction model has a unique solution, which was proved by Han et al. [19] (see also [20, 49, 7]) for the Kelvin-Voigt rheology, used also here. It is worth mentioning that uniqueness results for linear elastic solids are only available when assuming some bounds for contact problem parameters [24, 26, 43]. In fact, any ‘enough dissipative’ rheology leading to a parabolic-type equation (e.g.Jeffreys’ or the so-called 4-parameter solid) will serve similarly but the less dissipative rheologies leading to a rather hyperbolic equation (like Maxwell’s, Boltzmann’s or Burgers’), would not guarantee well-posedness of the model. Sometimes the former type of rheology is referred to as solid-type and the latter as fluid-type, cf. [40].

Eventually, in Section 4, the algorithm exploiting the above constructed incremental QP (solved here by a conjugate gradient based method) is tested in numerical simulations for various types of contacts, including receding contact and conforming contact problems.

The paper also includes tree appendices. In Appendix A, we give some details of SGBEM implementation for spatial discretization of a suitable system of Boundary Integral Equations (BIE), in Appendix B, we present the block structure of the displacement-traction map for interface variables, in Appendix C, we sketch some theoretical aspects of the model concerning numerical stability and convergence of the numerical implementation.

## 2 Contact problem with friction at small strains

For the notational simplicity, only two bodies in contact will be considered. Let these bodies occupy domains with bounded Lipschitz boundaries , Figure 1. Later, we will consider only but some parts allow also .

Let denote the unit outward normal vector defined a.e. at . For , denotes the unit tangential (surface) vector such that it defines anti-clockwise orientation of .

The contact zone is defined as the common part of and , i.e. . Notice, that this definition, adopted for the sake of formulation simplicity, covers just receding and conforming contact types. Nevertheless, the present computational implementation covers also the case of advancing contact. The Dirichlet and Neumann boundary conditions are defined on the outer boundary parts, respectively, prescribing displacements as on and tractions as on at a current time , cf. (2.2b,c) below. Of course, we assume that and are disjoint and relatively open parts of , disjoint also with , and . We also assume that the Dirichlet part is far from the contact boundary, i.e. . The difference, i.e. a gap, on of (traces of) functions defined on and will be denoted by . In particular, the gap of displacements on the contact boundary means . Further, we will use the gap of the normal and tangential displacements

 (2.1)

We also use the convention that the ‘dot’ will stand for the partial time derivative .

The classical formulation of the evolution boundary-value problem we will consider is the following:

 divση=0,  ση=Dηe(\LARGE.uη)+Cηe(uη),  e(u):=12(∇u)⊤+12∇u on Ωη, η=A,B, (2.2a) uη=gη\tiny\rm D(t) on Γη\tiny\rm D, (2.2b) pη=fη\tiny\rm N(t)    with the % traction  pη:=ση→nη on Γη\tiny\rm N, (2.2c) on Γ\tiny\rm C; (2.2d) with Cη (resp. Dη) the positive-definite symmetric 4th-order tensors of elastic (resp. viscous) moduli on η-domain, μ>0 a coefficient of friction and γ:R→R+ a decreasing smooth and convex function describing the compression energy of the normal-compliance contact. The normal and tangential components of the traction stress defined in (2.2c) are denoted by σn and σt in (2.2d), respectively. The evolution boundary-value problem (2.2) is to hold for a.a. time instances t∈[0,T] with T a fixed time horizon. We will further consider an initial-boundary-value problem by prescribing the initial conditions at time t=0 for the displacement: uη(0,⋅)=uη0 on  Ωη,  η=A,B. (2.2e)

It should be emphasized that the classical formulation (2.2) is rather formal in some aspects and gain a sense only in a suitable weak formulation. In particular, the jump of stress vector is well defined only in the space , and is imposed to equal zero, in contrast to which has a good meaning in the space if is qualified as (2.3) below, being equal to . It should be noted that the jump of the stress tensor is quite difficult to define and it is actually not needed in the formulation, additionally it does not have to vanish, while the jump of the stress vector does.

As for the normal-compliance model, for some purposes of this paper, we can assume that, for some for , otherwise , satisfying the following growth/continuity conditions:

 ∃C∈R (2.3a) ∀u,˜u∈R:∣∣γ′(u)−γ′(˜u)∣∣≤C(|u|q−2+|˜u|q−2)∣∣u−˜u∣∣. (2.3b)

Typically, for unilateral contact and, e.g., satisfies (2.3) provided . If , may have an arbitrary polynomial growth, which may imitate a finite interpenetration. Let us remark that the finite interpenetration up the depth, say, , means that and for . This is a quite realistic replacement of the Signorini contact, although any rigorous analysis for the evolutionary model does not seem to be available yet. Anyhow, especially for an efficient implementation by QP, we will later confine ourselves to a piecewise quadratic case, i.e. in particular .

For lucidity, it is noteworthy that the evolution boundary-value problem (2.2a) with (2.2d), when written in a weak formulation, takes the form of a doubly-nonlinear evolution governed by a convex stored energy and a nonsmooth convex potential of dissipative forces . In the so-called Biot-equation form, it reads as the differential inclusion

 (2.4)

for and , where the symbol refers to partial subdifferentials relying on the convexity of and . It involves the functionals , , and in the specific form:

 E(g\tiny\rm D;u):=E% B(g\tiny\rm D;u)+E\tiny\rm C(u)  % with   E\tiny\rm C(u):=∫Γ\tiny\rm Cγ([[u]]n)dΓ, (2.5a) (2.5b) (2.5c) (2.5d) ⟨F(t),v⟩:=∑η=A,B∫Γη\tiny\rm Nfη\tiny\rm N(t)vηdΓ∀v=(v\tiny A,v\tiny B), vη∈H1(Ωη;Rd), η=A,B. (2.5e)

Note that in (2.5a) and (2.5b) we split into the bulk part and the contact-boundary part to be used later especially in Section 3.2, while in (2.5c) and (2.5d) we denoted the split of the functional according to the degree of homogeneity of its terms with respect to the displacement rate (the subscript index). The dissipation rate reflects these different homogeneity degrees, namely , as occurs in (2.10) below.

We intentionally do not want to write (2.2a-d) in terms of some shifted solution satisfying the homogeneous Dirichlet conditions on because this (otherwise usual) procedure would give rise a non-homogeneous right-hand side in the transformed equilibrium equation (2.2a) and, thus, it would not be well suited with BEM for spatial discretization we want to use later.

The weak formulation of (2.4) is simple if , namely

 ∀a.a.t∈[0,T]  ∀% \small˜\normalsizeu∈H1(Ω\tiny A∪Ω\tiny B;Rd),  \small˜\normalsizeu|Γ\tiny\rm D% =0: (2.6)

For a general case , like in [9], one should use the formulation as a system:

 ∀a.a.t∈[0,T]:    ξ(t)+ζ(t)=F(t),   and (2.7a) (2.7b) ∀\small˜\normalsizeu% ∈H1(Ω\tiny A∪Ω\tiny B;Rd), \small˜\normalsizeu|Γ\tiny\rm D=0:  ⟨ζ(t),% \small˜\normalsizeu−u(t)⟩+E(g\tiny\rm D(t);\small˜\normalsizeu)≥E(g\tiny\rm D(t);u(t)). (2.7c)

However, one of the advantages of BIE and BEM is that the Dirichlet conditions are ‘hidden’ into the stored-energy functional so that time-varying Dirichlet condition does not represent any problem if is far from , cf. (3.17) and (3.19).

Here, when we want to write either the energy balance or (2.7) in a simpler form (2.6), we use such shift only for testing (2.4). Let define a (suitably smooth) prolongation of the boundary condition . Since is assumed, we can always consider on so that this transformation does not influence (2.2d) and then also the -term in (2.8) below. Then, has zero traces on and, in view of (2.4), satisfies

 ∂% \large.uR1(uz;\LARGE.uz)+R′2(\LARGE.uz+\LARGE.u\tiny\rm D)+∂uE(g\tiny\rm D(t);uz+u% \tiny\rm D)∋F(t), (2.8)

where stands for Gâteaux derivative. Now, testing (2.8) by having zero traces on becomes legal and renders

 (2.9)

and by integrating it over the time interval with the initial conditions (2.2e) provides the following energy balance, by substituting according to its definition,

 E(g\tiny\rm D(T);u(T))+∫T0(R1(u;\LARGE.u)+2R2(\LARGE.u))%dt=E(g\tiny\rm D(0);u0)+∫T0(⟨F(t),\LARGE.u−\LARGE.u\tiny% \rm D⟩+∑η=A,B∫Ωη(Cηe(uη)+Dηe(\LARGE.uη)):e(% \LARGE.uη\tiny\rm D)dΩ)dt. (2.10)

This expresses energy conservation in terms of the original displacement : the sum of the energy stored in the system at time and the energy dissipated during the time interval due to friction and viscous processes equals the energy stored at the system at time plus the work made by the external boundary loadings and during the same time interval.

## 3 Numerical approximation and computer implementation

The numerical procedure devised for solving the above problem considers time and spatial discretizations separately, as usual. The procedure is formulated in terms of the boundary data only, with the spatial discretization carried out by SGBEM. We restrict generality of the visco-elastic materials to homogeneous solids to allow for BEM-treatment, and further of the normal-compliance response by assuming a quadratic one. Eventually we also confine ourselves to 2D problems (i.e. ) in order to facilitate an efficient implementation by SGBEM combined with quadratic programming (QP); note that in 3D problems, the structure of the linear constraints (3.21) and (3.23) below cannot be obtained, although efficient tools do exist for this case, too, cf. Remark 3.2 below.

### 3.1 Semi-implicit discretization in time

For the sake of simplicity, the time-stepping scheme is initially defined by a fixed time step size such that is integer. Yet, it should be remarked that in specific calculations mixing stick/slip/jump regimes with typically very different time scale, a certain adaptivity is desirable and here very simple (as the inertial term is neglected). The velocity is approximated by the finite difference , where denotes the approximate solution at the discrete time . We discretise the differential inclusion (2.4) by the semi-implicit recursive scheme

 ∂% \large.uR(uk−1;uk−uk−1τ)+∂uE(g\tiny\rm D(kτ);uk)∋F(kτ),     k=1,...,T/τ,     u0=u0. (3.1)

We consider fixed thorough the whole section and thus, for notational simplicity, we will not explicitly emphasize it when writing . Occurrence of rather than , i.e. the semi-implicit formula rather than fully-implicit one, brings a variational structure to the incremental recursive problem (3.1) in the sense that (3.1) represents the first-order optimality condition for the strictly convex nonsmooth functional

 Hk(u):=τR(uk−1;u−uk−1τ)+E(g\tiny\rm D(kτ);u)−⟨F(kτ),u⟩. (3.2)

Thus, we can take the unique minimizer of as the solution to (3.1).

An analogy to (2.10) can be obtained by making a test of (3.1) with the transformation by shifting the solution by into the form like (2.8). Like the test of (2.8) by , we now make a test by at the time instant , where . In view of (2.5), it yields:

 ∫Γ\tiny\rm C−μγ′([[uk−1]]n)⋅∣∣[[uk−uk−1]]t∣∣dΓ +∑η=A,B1τ∫Ωηe(uη,kz−uη,k−1z+uη,k\tiny\rm D−uη,k−1\tiny\rm D):Dη:e(uη,kz−uη,k−1z)dΩ +∑η=A,B∫Ωηe(uη,kz+uη,k\tiny\rm D):Cη:e(uη,kz−uη,k−1z)dΩ +∫Γ\tiny\rm Cγ′([[uk]]n)⋅([[uk−uk−1]]n)dΓ=∑η=A,B∫Γη\tiny\rm Nfη\tiny\rm N⋅(ukz−uk−1z)dΓ. (3.3)

It provides the following energy estimate, returning back to the original :

 R1(uk−1;uk−uk−1)+2τR2(uk−uk−1)+E(g\tiny\rm D(kτ),uk)≤E(g\tiny\rm D((k−1)τ),uk−1) +∑η=A,B∫Ωη(Cηe(uη,k−1)+Dηe(uη,k−uη,k−1τ)):e(uη,k\tiny\rm D−uη,k−1\tiny\rm D)dΩ +∑η=A,B12∫Ωηe(uη,k\tiny\rm D−uη,k−1\tiny\rm D):Cη:e(uη,k\tiny\rm D−uη,k−1\tiny\rm D)dΩ+⟨F(kτ),uk−uk−1⟩. (3.4)

The above inequality is useful when we permit the variations of the time step parameter between subsequent time steps providing a form of adaptivity in the algorithm: the inequality is required to be satisfied as an approximation of equality with a given tolerance – if the difference is greater than , the time step is diminished, if the difference is much more less than , the time step is enlarged, otherwise it is kept the same.

In the derivation of the inequality (3.4), the convexity of the functional from (2.5a) and (2.5b) has been used, which provides, e.g., the following relation for the pertinent boundary term:

 ∫Γ\tiny\rm Cγ′([[uk]]n)⋅[[uk−uk−1]]ndΓ≥∫Γ\tiny\rm Cγ([[uk]]n)dΓ−∫Γ\tiny\rm Cγ([[uk−1]]n)dΓ. (3.5)

The residuum in (3.4) can serve for varying adaptively the time step to keep this residuum under a pre-selected tollerance. This, however, needs the energy conservation in the limit problem, cf. Remark C.1 in the Appendix C.

### 3.2 Implementation of viscosity towards BIE

Now we make a slight restriction of generality by assuming that

 Dη=χCη,        χ>0 given. (3.6)

In most materials (like metals) the relaxation time is small comparing to outer-loading time scale but we assume that it is rather large to affect the tangential contact with friction. We assume to be the same for both subdomains; note that the potential (3.11) cannot be formulated if would be different in each subdomain; anyhow, some considerations (and computations) can be performed for non-potential problems, too, as suggested in [40, Sect.3.3]. This allows for using the simple trick from [45, Remark 6.2], which can be extended even for more complex rheologies as shown in [40], which facilitates usage of standard static BIEs even for quasi-static evolution problems.

The simple viscosity ansatz (3.6) is chosen in order to exploit the reformulation of the viscoelastic problem in the bulk in terms of an recursive static elastic problem in the bulk, similarly as in [40], which is solved by the conventional elastostatic SGBEM here. For this aim, let us introduce a new variable (a ‘fictitious’ displacement) at the time level , as

 vk:=uk+χuk−uk−1τ. (3.7)

Realizing that

 uk=ττ+χvk+χτ+χuk−1     and also     uk−uk−1τ=vk−uk−1τ+χ (3.8)

and that, thanks to (3.6),

 (3.9a) with  % \small˜\normalsizeg\tiny\rm D(kτ):=g\tiny\rm D(kτ)+χτ(g\tiny\rm D(kτ)−g\tiny% \rm D((k−1)τ)), (3.9b)

we can re-write (3.1) in terms of this new variable as:

 (3.10)

here, we have also used that is (positive) homogeneous of degree 1 so that is homogeneous of degree 0 and the factor from (3.8) disappears from the -term in (3.10). The recursive scheme then involves also the first equality in (3.8) and starts with as in (3.1). It is important that all bulk integrals are now contained in the ‘static’ part . The inclusion (3.10) has again the variational structure employing the convex, nonsmooth functional:

 Kk(v) =∑η=A,B∫Ωη12e(vη):Cη:e(vη)dΩ−∑η=A,B∫Γη\tiny\rm Nfη\tiny\rm N⋅vηdΓ (3.11)

where . Like before, (3.10) represents the 1st-order optimality condition for to be a minimizer of subject to the boundary condition for .

The BIE-method used here is based on a formulation of the problem only on the contact boundary . To this goal, we abbreviate the displacement and the fictitious displacement gaps by

 z:=[[u]], (3.12a) w:=[[v]], (3.12b)

Moreover, we define the linear mapping with being the weak solution to the Transmission Boundary Value Problem (TBVP):

 div(Cηe(vη))=0 in Ωη,\ \ η=A,B, (3.13a) v=gη\tiny\rm D on Γη\tiny\rm D,\ \ η=A,B, (3.13b) pη=fη\tiny\rm N    with   pη=(Cηe(v))⋅→nη on Γη\tiny\rm N,\ \ η=A,B, (3.13c) [[v]]=wp\tiny A+p\tiny B=0} on Γ\tiny\rm C. (3.13d)

Later in Section 3.3-3.4, we will solve (3.13) approximately by SGBEM, see Appendix A, using the formulation for multi-domain problems in [54]. To this goal, it is desirable to replace all bulk integrals in (3.10). The resulted recursive inclusion will be then formulated only on , being written exclusively in terms of and only:

 ∂\large.wR1(zk−1n;wkt−zk−1t)+[EB]′w(\small˜\normalsizeg\tiny\rm D(kτ);wk)+E′% C(τwkn+χzk−1nτ+χ)∋F′w(kτ), (3.14a) zk=ττ+χwk+χτ+χzk−1,k=1,...,T/τ; (3.14b)

for the recursion (3.14b), cf. (3.8). Of course, for , we start with the initial condition , calculate from (3.14a), then from (3.14b), and continue for etc.. In (3.14a), we have used

 EB(\small˜\normalsizeg\tiny\rm D;w):=EB(% \small˜\normalsizeg\tiny\rm D;S(\small˜\normalsizeg\tiny\rm D,f\tiny\rm N,w)), (3.15a) EC(wn)=E%C([[v]]n):=EC(v), F(t,w):=⟨F(t),S(\small˜\normalsizeg\tiny\rm D,f\tiny\rm N,w)⟩. (3.15b)

Note that and depend only on the values on and thus the mapping is not used in the definition of and . Note also that, in contrast to , the functional is smooth. The functional (3.11) can be written as a function of as

 Kk(w) :=∫Γ\tiny\rm C(τ+χτγ(ττ+χwn+χτ+χzk−1n)−μγ′(zk−1n)∣∣wt−zk−1t∣∣ +∑η=A,B12∫Γη\tiny\rm Dpη(\small˜% \normalsizeg% \tiny\rm D(kτ),f\tiny\rm N(kτ),w)⋅˜gη\tiny\rm D(kτ)dΓ (3.16)

where and denote the displacement and traction, respectively, resulting in TBVP (3.13) through the (generalized) Poincaré-Steklov operator which is described in Appendix B and whose representation is given by the relation (B.6).

In fact, the expressions and are to be understood rather as -dualities than Lebesgue integrals. We dare to keep this non-precise but more conventional notation in engineering throughout this article.

Let us further realize that holds if solves the boundary-value problem (3.13). The optimality conditions for minimization of represent a weak formulation of the discrete problem (3.10). The advantage of (3.14a) in contrast to (3.10) is that the constraints coming from the Dirichlet conditions (which prevent to write the weak formulation of the original problem generally as in (2.6)) are eliminated. To be more specific:

 +⟨pk,\small˜\normalsizew−wk⟩≥R1(zk−1n;wkt−zk−1t)+⟨F′w(kτ),\small˜\normalsizew−wk⟩ (3.17)

with the traction being from (3.16).

### 3.3 Spatial discretization and SGBEM

The role of the SGBEM in the present computational procedure is to provide a complete boundary-value solution from the given boundary data to calculate the elastic strain energy in these domains. Thus, at each time step and at each iteration of the minimization algorithm, the SGBEM code calculates, in a similar way as in [54], unknown tractions along and unknown displacements along , assuming the displacement gap on to be known from the used minimization procedure.

Here, we rely on the assumption that the material is homogeneous in each solid to allow for BEM-treatment. Then, the present implementation of SGBEM, deduced from [54, 55], is briefly described in Appendix A. In fact, the necessary SGBEM calculations are besad on the solution of (A.3). The variables of the (modified) displacement , the displacement gap and the traction appearing in (3.16) are approximated as given by the formula (A.5), i.e.

 vη(x)=∑nNηψn(x)vηn,w(x)=∑n\tiny\rm CN\tiny Bψn\tiny\rm C(x)wn\tiny\rm C,pη(x)=∑mNηφm(x)pηm,

with nodal variables gathered in respective vectors , and .

The SGBEM-approximation of is the discretized energy from the equation (3.16) with approximations (A.5), let us denote it by with referring to the mesh-size of the BE discretization of , is then

 , Kkh(w) :=∫Γ\tiny\rm C(τ+χτγ(ττ+χ[wh]n+χτ+χ[wk−1h]n)−μγ′(zk−1n,n\tiny\rm C)∣∣[wh]t−[zk−1h]n∣∣ −12p\tiny Bh(\small˜\normalsizeg\tiny\rm D(kτ),f\tiny\rm N(kτ),wh)⋅wh)dΓ −∑η=A,B∫Γη\tiny\rm N12fη\tiny\rm N⋅vηh(\small˜\normalsizeg\tiny\rm D(kτ),f\tiny\rm N% (kτ),wh)dΓ +∑η=A,B12∫Γη\tiny\rm Dpηh(\small˜% \normalsizeg\tiny% \rm D(kτ),f\tiny\rm N(kτ),wh)⋅% \small˜\normalsizegη\tiny\rm D(kτ)dΓ, (3.18a) with wh=wh(w) from (A.5), i.e. wh(x)=∑n\tiny\rm CN\tiny Bψn\tiny% \rm C(x)wn\tiny\rm C, (3.18b)

and with and obtained by the approximate (generalized) Poincaré-Steklov operator of Appendix B, also used in solution of contact problems by BEM [33, 32]. Here also, is the master surface for displacements, similarly as in [54], and verifies the recursive formula analogous to (3.14b), i.e. . The BEM-version of (3.17) seeks (denoted here for notational simplicity without ) such that

 (3.19a) where Vh={w∈W1,∞(Γ\tiny\rm C;Rd); w in the form (???)}⊂H1/2(Γ% \tiny\rm C;Rd) and pkh is defined by pηh(\small˜\normalsizeg\tiny\rm D(kτ),f\tiny\rm N(kτ),wh) from (3.18). Of course, again we use the recursion (3.14b), which is written in the terms of the spatial discretizations as zk=ττ+χwk+χτ+χzk−1, (3.19b)

provided the initial condition is assumed so that for sufficiently small (or it is approximated by an element of the space ).

### 3.4 QP-minimization algorithm in 2-dimensional problems

Once all the boundary data are known, the energy of the state given by in (3.16) can be calculated. It is worth seeing how it can be carried out in the present implementation provided we confine ourselves to the most common case of linearly-responding compliance on , i.e. piecewise quadratic :

 γ(g)=kg2min(0,g)2,    kg>0 % given. (3.20)

Further, except Remark 3.2, we confine ourselves to 2D problems where has, after discretization, a polyhedral graph, and similarly also , with a help of (3.20), leads after discretization to polyhedral constraints, see (3.23) below. All this will allow for usage of efficient quadratic-programming (QP) algorithms.

First, let us re-consider the non-quadratic terms in . A rather classical trick of handling such unpleasant terms by displacing them into a position of constraints by additional variables; Mosco-type transformation, although the original work [36] was designed rather for the optimality conditions as variational inequalities, which works even for nonpotential problems. For the rate-independent problems where the unpleasant terms are piecewise affine, it was used in [44, Lemma 4] and then also in [46, Lemma 4]. Here we have also piecewise quadratic terms, but they need only a simple modification to consider their square root to obtain again the linear-constraint structure but then add the square power of these auxiliary variables into the cost function. Let these additional auxiliary variables handling the piecewise linear and the piecewise quadratic terms be denoted as and , respectively, and the following constraints hold for and :

 α−wt α+wt (3.21a) β ≥0, β+wn ≥−χτ[[uk−1]]n. (3.21b)

For the discretization, the approximation formulas for both auxiliary parameters and given by the pertinent BE mesh should be considered. In what follows, the same mesh and approximation as used in (A.5) for displacements on the boundary part is considered. The approximation formulas can be written in the form

 α(x)=∑n\tiny\rm CN\tiny Bψn\tiny% \rm C(x)αn\tiny\rm C,β(x)=∑n\tiny% \rm CN\tiny Bψn\tiny\rm C(x)βn% \tiny\rm C, (3.22)

where , are the nodal unknowns associated to the node .

The restrictions (3.21) written for the discretization can be obtained as relations between the nodal values, realizing that the right-hand side terms can be expressed in terms of applying to them the transformation (3.8). If runs through all nodes in the possible contact zone with respect to the boundary , the nodal relations read

 αn\tiny\rm C−wt,n\tiny\rm C ≥−\small˜\normalsizeuk−1t,n\tiny\rm C, αn\tiny\rm C+wt,n\tiny\rm C (3.23a) βn\tiny\rm C ≥0, βn\tiny\rm C+wn,n\tiny\rm C ≥−χτ[[\small˜\normalsizeuk−1n\tiny\rm C]]n, (3.23b)

where in the nodal values at the node subscript indices and are used to distinguish between tangential and normal components, respectively.

Then, in terms of these new variables, from (3.18) takes the following form denoted by :

 ˆKkh(w,α,β)=∫Γ\tiny\rm C(12τkgτ+χ(∑n\tiny\rm CN% \tiny Bψ,n\tiny\rm C(x)β