Abstract
The quasistatic normalcompliance contact problem of isotropic homogeneous linear viscoelastic bodies with Coulomb friction at small strains in KelvinVoigt rheology is considered. The discretization is made by a semiimplicit 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 relaxationtime 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 piecewisequadratic, which after a certain transformation gets the quadratic programming (QP) structure. However, it would lead to secondorder 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 normalcompliance contact problem of
viscoelastic bodies with Coulomb friction
implemented by QP and SGBEM.
Roman Vodička^{1}^{1}1Technical University of Košice, Civil Engineering Faculty, Vysokoškolská 4, 042 00 Košice, Slovakia, roman.vodicka@tuke.sk, Vladislav Mantič^{2}^{2}2University of Seville, School of Engineering, Camino de los Descubrimientos s/n, 41092 Seville, Spain, mantic@etsi.us.es, Tomáš Roubíček^{3}^{3}3Charles University, Mathematical Institute, Sokolovská 83, CZ186 75 Praha 8, Czech Republic, roubicek@karlin.mff.cuni.cz^{4}^{4}4Institute 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 socalled normal compliance is a wellaccepted 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 normalcompliance 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 quadraticprogramming or a 2ndorderconeprogramming problem.
When nonlinearities are restricted to the boundaries, smallstrain 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 normalcompliance 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 socalled 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 viscoelastic behaviour of the bulk but implemented so that the elastostatic 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 semiimplicit formula in time (Sect. 3.1).

Implementation of BEM, assuming the ratio of the viscosity and elasticity moduli is a given relaxationtime 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 linearlyresponding 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 (apriori 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 KelvinVoigt 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 parabolictype equation (e.g.Jeffreys’ or the socalled 4parameter 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 wellposedness of the model. Sometimes the former type of rheology is referred to as solidtype and the latter as fluidtype, 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 displacementtraction 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 anticlockwise 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 boundaryvalue problem we will consider is the following:
(2.2a)  
(2.2b)  
(2.2c)  
(2.2d)  
with (resp. ) the positivedefinite symmetric 4thorder tensors of elastic (resp. viscous) moduli on domain, a coefficient of friction and a decreasing smooth and convex function describing the compression energy of the normalcompliance contact.  
The normal and tangential components of the traction stress defined in (2.2c) are denoted by and in (2.2d), respectively. The evolution boundaryvalue problem (2.2) is to hold for a.a. time instances with a fixed time horizon. We will further consider an initialboundaryvalue problem by prescribing the initial conditions at time for the displacement:  
(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 normalcompliance model, for some purposes of this paper, we can assume that, for some for , otherwise , satisfying the following growth/continuity conditions:
(2.3a)  
(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 boundaryvalue problem (2.2a) with (2.2d), when written in a weak formulation, takes the form of a doublynonlinear evolution governed by a convex stored energy and a nonsmooth convex potential of dissipative forces . In the socalled Biotequation 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:
(2.5a)  
(2.5b)  
(2.5c)  
(2.5d)  
(2.5e) 
Note that in (2.5a) and (2.5b) we split into the bulk part and the contactboundary 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.2ad) in terms of some shifted solution satisfying the homogeneous Dirichlet conditions on because this (otherwise usual) procedure would give rise a nonhomogeneous righthand 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
(2.6) 
For a general case , like in [9], one should use the formulation as a system:
(2.7a)  
(2.7b)  
(2.7c) 
However, one of the advantages of BIE and BEM is that the Dirichlet conditions are ‘hidden’ into the storedenergy functional so that timevarying 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
(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,
(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 viscoelastic materials to homogeneous solids to allow for BEMtreatment, and further of the normalcompliance 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 Semiimplicit discretization in time
For the sake of simplicity, the timestepping 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 semiimplicit recursive scheme
(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 semiimplicit formula rather than fullyimplicit one, brings a variational structure to the incremental recursive problem (3.1) in the sense that (3.1) represents the firstorder optimality condition for the strictly convex nonsmooth functional
(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:
(3.3) 
It provides the following energy estimate, returning back to the original :
(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:
(3.5) 
The residuum in (3.4) can serve for varying adaptively the time step to keep this residuum under a preselected 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
(3.6) 
In most materials (like metals) the relaxation time is small comparing to outerloading 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 nonpotential 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 quasistatic 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
(3.7) 
Realizing that
(3.8) 
and that, thanks to (3.6),
(3.9a)  
(3.9b) 
we can rewrite (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:
(3.11) 
where . Like before, (3.10) represents the 1storder optimality condition for to be a minimizer of subject to the boundary condition for .
The BIEmethod 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
(3.12a)  
(3.12b) 
Moreover, we define the linear mapping with being the weak solution to the Transmission Boundary Value Problem (TBVP):
(3.13a)  
(3.13b)  
(3.13c)  
(3.13d) 
Later in Section 3.33.4, we will solve (3.13) approximately by SGBEM, see Appendix A, using the formulation for multidomain 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:
(3.14a)  
(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
(3.15a)  
(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
(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 nonprecise but more conventional notation in engineering throughout this article.
Let us further realize that holds if solves the boundaryvalue 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:
(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 boundaryvalue 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 BEMtreatment. 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.
with nodal variables gathered in respective vectors , and .
The SGBEMapproximation of is the discretized energy from the equation (3.16) with approximations (A.5), let us denote it by with referring to the meshsize of the BE discretization of , is then
,  
(3.18a)  
with from (A.5), i.e. ,  (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 BEMversion of (3.17) seeks (denoted here for notational simplicity without ) such that
(3.19a)  
where and is defined by from (3.18). Of course, again we use the recursion (3.14b), which is written in the terms of the spatial discretizations as  
(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 QPminimization algorithm in 2dimensional 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 linearlyresponding compliance on , i.e. piecewise quadratic :
(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 quadraticprogramming (QP) algorithms.
First, let us reconsider the nonquadratic terms in . A rather classical trick of handling such unpleasant terms by displacing them into a position of constraints by additional variables; Moscotype transformation, although the original work [36] was designed rather for the optimality conditions as variational inequalities, which works even for nonpotential problems. For the rateindependent 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 linearconstraint 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 :
(3.21a)  
(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
(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 righthand 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
(3.23a)  
(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 :