Boundary Integral Formulation of Frictionless Contact Problems Based on an Energetic Approach and Its Numerical Implementation by the Collocation BEM

ORIGINAL RESEARCH published: 27 November 2018 doi: 10.3389/fbuil.2018.00056 Frontiers in Built Environment | www.frontiersin.org 1November 2018 | Volume 4 | Article 56 Edited by: Vagelis Plevris, OsloMet - Oslo Metropolitan University, Norway Reviewed by: Ivan Argatov, Technische Universität Berlin, Germany Andreas Kampitsis, Faculty of Engineering, Imperial College London, United Kingdom Cyril Fischer, Institute of Theoretical and Applied Mechanics (ASCR), Czechia *Correspondence: Vladislav Manti ˇ c [email protected] Specialty section: This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment Received: 03 July 2018 Accepted: 28 September 2018 Published: 27 November 2018 Citation: Panagiotopoulos CG, Manti ˇ c V, García IG and Graciani E (2018) Boundary Integral Formulation of Frictionless Contact Problems Based on an Energetic Approach and Its Numerical Implementation by the Collocation BEM. Front. Built Environ. 4:56. doi: 10.3389/fbuil.2018.00056 Boundary Integral Formulation of Frictionless Contact Problems Based on an Energetic Approach and Its Numerical Implementation by the Collocation BEM Christos G. Panagiotopoulos1, Vladislav Manti ˇ c2*, Israel G. García2and Enrique Graciani2 1Institute of Applied & Computational Mathematics, Foundation for Research and Technology-Hellas, Heraklion, Greece, 2Grupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Sevilla, Spain A unified methodology to solve problems of frictionless unilateral contact as well as adhesive contact between linear elastic solids is presented. This methodology is based on energetic principles and is casted to a minimization problem of the total potential energy. Appropriate boundary integral forms of the energy are defined and the quadratic problem form of the contact problem is proposed. The problem is solved by the collocation boundary element method (BEM). To solve the quadratic problem two algorithms are developed, both being variants of the well-known conjugate gradient algorithm. The difference between them is given by the explicit construction or not of the quadratic-problem matrix. This matrix has the same physical meaning as the stiffness matrix commonly used in the context of the finite element method (FEM). Both symmetric and non-symmetric formulations of this matrix are presented and discussed, showing that the non-symmetric one provides more accurate results. The present procedure, in addition to its interest by itself, can also be extended to problems where dissipative phenomena take place such as friction, damage and plasticity. Essential steps of the numerical implementation are briefly presented and the numerical solutions of some standard problems of frictionless contact are given and compared to those obtained by other well-known BEM and FEM procedures for contact problems. Keywords: unilateral frictionless contact, adhesive contact, linear elasticity, boundary element method, stiffness matrix in BEM, total potential energy minimization 1. INTRODUCTION Contact problems (Johnson, 1985), often present in engineering applications, represent one of the most important and interesting topics of mechanics. Contact between deformable bodies is a complex and inherently non-linear problem. It is essentially a boundary phenomenon which has a strong effect on stress and displacement fields in the vicinity of the boundary of the contacting bodies. Phenomenologically one of the simplest models of contact is given by the Signorini (1959) contact condition, that imposes the strict non-penetrability of the bodies in contact. After early studies by Hertz (1896), many contact problems have been solved analytically. However, Panagiotopoulos et al. BEM Energetic Approach to Contact complicated geometries or loads, as well as, more complex interfacial conditions (e.g., adhesion, cohesion, friction, etc.), require the use of numerical methods (Wriggers, 2006; Yastrebov, 2013). Unilateral frictionless contact problems for deformable bodies are usually numerically solved by the finite element method (FEM) or boundary element method (BEM). For crack onset and growth modeling (see Panagiotopoulos et al., 2013; Roubíˇ cek et al., 2013), the first two authors with coworkers developed an energy based procedure implemented by BEM. In these applications, contact between crack faces should be considered, hence it is important to present in detail the theoretical background as well as details of the numerical procedure for problems of unilateral (i.e., Signorini) and adhesive frictionless contact. Hereinafter, this framework is referred to as Energetic approach for the solution of elastic Contact problems by BEM (EC-BEM). Unilateral frictionless contact stands for the case where there is no other material in the interfaces (=contact zones) between elastic bodies. For this case the (Signorini, 1959) condition models the exact non-penetrability of the bodies in contact, by the following relations, Eck et al. (2005): •The relative normal displacement at the interface cannot be larger than the initial distance between the bodies in contact. •Stresses can be transmitted only if contact takes place. •Only compressive normal tractions can be transmitted through the contact zone. Adhesive unilateral contact is represented by bodies connected across their common interfaces by a continuous distribution of springs, similar to the Winkler spring model, with (possibly) distinct normal and tangential elastic stiffnesses whose values range from zero to infinity. Several approaches for solving contact problems by BEM were developed in the past. In engineering literature an approach based on incremental-iterative procedures was developed and widely used from the very beginning (see París and Garrido, 1989; Katsikadelis and Kokkinos, 1993; París et al., 1995), among others. In this approach successive increments of load are applied and the Signorini contact condition is verified in each step by an iterative procedure. A BEM methodology for frictionless contact problems, based on the strong formulation of the problem, where the position of contact zones is defined through geometrical parameters is presented in Méchain-Renaud and Cimetière (2000). In contrast to these approaches, the present approach for solving frictionless contact problems is based on general energetic considerations. Usually this approach requires the minimization of the potential energy under unilateral constraints on displacements (Gurtin, 1972; Lazaridis and Panagiotopoulos, 1987; Panagiotopoulos and Lazaridis, 1987), although different approaches may also be found (see Kalker and Randen, 1972), where instead of the potential energy the minimum enthalpy principle is employed, formulated as a boundary integral equations method, valid for frictionless contact problems in the absence of dissipative processes. In the present work the contact problem is modeled as a quasi-static and rate-independent process, (Wriggers and Panagiotopoulos, 1972) (neglecting inertial and viscous forces) for linear elastic solids. In section 2 the energetic approach employed is presented. In section 3, the energetic principles formulated in boundary integral forms amenable to numerical implementation are introduced. A particular BEM implementation is presented in section 4. A procedure for computation of matrices in an explicit form that may be used in quadrature programming algorithms is described in section 5. Some details for analysis of a possible contact between elastic bodies, through interface elements, are presented in section 6, while specific features of the minimization algorithms are described in sections 7 through 10. Finally, in section 11, results of two-dimensional multi-region simulations are presented and compared to solutions obtained by other BEM and FEM techniques available in the literature, showing that the present framework is suitable to confront problems of contact between elastic bodies. 2. THEORETICAL BACKGROUND Consider a finite number Nof linear elastic bodies possibly in frictionless contact. Let these bodies be represented by mutually disjoint subdomains i⊂R2(i=1, ..., N) with Lipschitz boundaries ∂i=Ŵi(see Figure 1) and the unit outward normal and tangential vectors νiand τi, with τi=(−(νi)2, (νi)1). Let C(i)denote the elastic stiffness tensor of subdomain i,ui(x,t) the time dependent displacements on i,ei(ui)= 1 2(▽ui)T+ ▽uithe small strain tensor, σi(ui)=C(i)ei(ui) the stress tensor, pi(ui)=σi(ui)νithe boundary tractions, σi n= pi(ui)νiand σi t=pi(ui)τithe normal and tangential component of boundary tractions, respectively. Let Ŵi j⊂Ŵidenote the part of the boundary Ŵiwhich possibly can contact with j(j=0, ..., Nand j6= i), referred to as potential contact zone of iwith j. We also consider a possible contact with a rigid obstacle on a part of Ŵidenoted as Ŵi 0. We assume that the rest of Ŵiis partitioned into the Dirichlet and Neumann parts of the boundary, Ŵi Dand Ŵi N, respectively. Time-dependent boundary displacement and traction vectors, ui D(x,t) and pi N(x,t), are prescribed on Ŵi Dand Ŵi N, respectively. Thus, e.g., any admissible displacement ui(x,t) on ihas to be equal to the prescribed displacement ui D(x,t) on Ŵi D. Let ŴD=[ 1≤i≤N Ŵi Dand ŴN=[ 1≤i≤N Ŵi N. (1) In the simplest case of conforming or receding contact with a zero initial gap, Ŵi jcan be equal to Ŵj i, but in general, e.g., in the case of advancing contact or a positive gap, they do not coincide. The initial contact zone (in the undeformed configuration) between subdomains iand j, and also the active contact zone (the set of points on each subdomain which will enter in contact in the deformed configuration) both should be included in Ŵi jand Ŵj i. Specially, Ŵi jand Ŵj imay be larger than ∂i∩∂jfor modeling an advancing contact. Frontiers in Built Environment | www.frontiersin.org 2November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 1 | Schematic illustration of the geometry and of the notation for the 2D case of N=3 subdomains in contact including a rigid obstacle. We assume that an intermediate contact surface, denoted as Ŵij =Ŵji, can be defined for each couple of potential contact zones, with one-to-one and onto (bijective) and continuous mappings (projections) between Ŵi j,Ŵj iand Ŵij. A discretized version of such mappings is briefly introduced in section 6. A recent comprehensive revision of different projection techniques for a general contact configuration can be found in Yastrebov (2013). Additionally a normal direction is defined at each point of Ŵij by a unit normal vector νi j(= −νj i) (outward with respect to i). Recall that all Ŵi j,Ŵj iand Ŵij are defined in the undeformed configuration. In the simplest case of conforming or receding contact with a zero initial gap, typically, Ŵi j=Ŵj i=Ŵij. The contact zone ŴCin the present problem is defined by the union of all intermediate contact surfaces ŴC=[ 0≤i<j≤N Ŵij. (2) The initial normal distance (gap) between some points xi∈Ŵi j and xj∈Ŵj icorresponding each other by the above defined mappings is given by the scalar gap function gij =(xj−xi)νi j defined on the intermediate contact surface Ŵij. A positive value of the gap function means an initial separation whereas a negative value means an initial overlapping between the bodies in contact. This gap function is denoted as gnon ŴC. Similarly, we also define [u]ij =uj(xj)−ui(xi) on Ŵij as the signed difference of displacements (displacement jump) between some points xi∈Ŵi jand xj∈Ŵj icorresponding each other by the above defined mappings. The normal and tangential components [u]nand [u]ton ŴC, respectively, are defined as [u]n=[u]ijνi j and [u]t=[u]ijτi j. According to this definition, [u]n, referred FIGURE 2 | Constitutive relations for the unilateral frictionless contact considered. to as opening, decreases when subdomains iand jare getting closer, and [u]n+gn≥0 should hold. [u]tis referred to as sliding. Therefore, [u]ij =[u]nνi j+[u]tτi j. In fact, the model does not depend on the chosen orientation, and indices iand jare omitted when displacement jump is defined on ŴC. The relations between the components of the tractions and displacement jump on a unilateral contact surface are shown in Figure 2, which indicates that σn≤0 if [u]n+gn=0 (contact), no overlapping being allowed, whereas σn=0 if [u]n+gn>0 (separation), and σt=0 in both cases. The elastic strain energy Estored in the volume = S1≤i≤Niis E(u)= N X i=1 Ei(ui)= N X i=1Zi 1 2C(i)ei(ui):ei(ui) d. (3) Then, the total potential energy functional, also referred to as the stored energy functional, is defined as 5(t,u)=   E(u)−W(t,u) if u=uDon ŴD, [u]n≥ −gnon ŴC, ∞elsewhere, (4) where, W(t,u)=PN i=1RŴi Nui(x,t)·pi N(x,t) dSgives the work of external forces in , the body forces being neglected for the sake of simplicity. According to Fichera (1964, see also Kalker and Randen, 1972; Panagiotopoulos and Lazaridis, 1987), the minimizer u(t) of the total potential energy functional 5(t,u) at each time tis the solution of the unilateral frictionless contact considered. Remark 1 (Adhesive unilateral frictionless contact). We refer to adhesive unilateral frictionless contact as the case where a continuous distribution of springs exists on an interface, i.e., a contact zone. The normal and tangential elastic stiffnesses of an interface, κn≥0 and κt≥0, respectively, may have different values (Figure 3). The interface response is considered frictionless, since no frictional dissipation of energy is considered, the tangent stiffness describing an adhesive elastic behavior. The elastic energy stored in such an adhesive interface (layer) is EC(u)=ECn(u)+ECt(u), (5) Frontiers in Built Environment | www.frontiersin.org 3November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 3 | Constitutive relations for the adhesive unilateral contact considered. where the components corresponding to the normal and tangential directions, respectively, are ECn(u)=(RŴC κn 2[u]n2dSif [u]n≥ −gnon ŴC, ∞elsewhere, (6) and ECt(u)=ZŴC κt 2ut 2dS. (7) The total potential energy is augmented by the elastic energy stored in the adhesive interface layer giving 5(t,u)=E(u)+EC(u)−W(t,u). (8) The relations between the components of the tractions and displacement jump on an adhesive interface layer are shown in Figure 3, which indicates that σn≤0 if [u]n+gn=0 (contact), whereas σn=κn[u]nif [u]n+gn>0 (separation). Additionally, σt=κt[u]tin shear. In this case, gnrepresents the width of the adhesive layer. If [u]n>0 then the layer is under tensile stresses and we do not allow to break it, whereas if [u]n<0 then the layer is under compressive stresses and we assume a linear behavior until the two domains enter in contact, no overlapping being allowed. It is easy to realize, that for vanishing interface stiffnesses, κn→0 and κt→0, the above described Signorini unilateral contact model is recovered, cf. Figures 2,3. 3. ENERGY PRINCIPLES IN BOUNDARY INTEGRAL FORMS For given body forces Fiin the elastic domain i, the following equilibrium condition holds: −divσi(ui)=Fiin i, which multiplied by some virtual displacement field ˜uiand integrated by parts over igives Zi σi(ui):ǫi(˜ui)d−ZŴi σi(ui)·νi˜uid x =ZiFi˜uid. (9) This expression may be seen as the principle of virtual work, for virtual displacements ˜ui. The dependence of stresses and subsequently tractions on displacements is explicitly indicated in the notation used. By assuming zero body forces (Fi=0) (in the present work) and choosing ˜ui=ui, the volume integral in (3) is replaced by the boundary integral Ei(ui)=1 2ZŴi uipi(ui) dS. (10) Furthermore, the total potential energy for iis 5i(t,ui)=Ei(ui)−ZŴNiui(x,t)pi N(x,t) dS. (11) In order to formulate the problem as a quadratic programming problem involving boundary values only, we further proceed manipulating (11). Similar formulations may be found in Antes and Panagiotopoulos (1992) (see chapter 8), where also a proof of existence and uniqueness of solutions is provided, and in Aour et al. (2007), however without the extent and the analysis introduced in the present work. In the following, the BVP for a sub-domain iis considered in a similar way as in Panagiotopoulos et al. (2013). Hereinafter in this section we will omit superindex i, for the sake of simplicity. Let uηand pη, respectively, denote the displacement and traction solutions of this BVP restricted to Ŵη,η=C, D and N, e.g. uC=u|ŴCand pD=p|ŴD. We study here a mixed-type operator Mwhich formally assigns (pC,pD,uN) to the known boundary data (uC,uD,pN) of the original BVP POshown in Figure 4, and may be partitioned using the following block structure as pC pD uN!= MCC MCD MCN MDC MDD MDN MNC MND MNN ! uC uD pN!. (12) The columns of the aforementioned block operator Mare associated to the sub-problems Pηdefined in Figure 4. The displacement solution of a subproblem Pηis denoted as uη. From the principle of superposition the displacement solution of the full problem POmay be reconstructed by the sum: u=uC+uD+uN. (13) In view of (10) and (11), the total potential energy for the mixed type BVP POmay be written in an expanded form as 5(t,u)=1 2ZŴC uC(u)pC(u)dS+1 2ZŴD uDpD(u)dS−1 2ZŴN uN(u)pNdS, (14) where we tacitly assumed that ufulfills the Dirichlet and Neumann boundary conditions. By substituting the unknown data for problem POfrom (12), the total potential energy can be written as 5(t,uC)=1 2ZŴC uCMCCuCdS+ZŴC uCMCDuDdS+ZŴC uCMCNpNdS Frontiers in Built Environment | www.frontiersin.org 4November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 4 | Solution of a mixed BVP POfor a single elastic domain given as a superposition of the solutions of the three sub-problems PC,PDand PN. +ZŴD uDMDCuCdS+ZŴD uDMDDuDdS+ZŴD uDMDNpNdS −ZŴN pNMNCuCdS−ZŴN pNMNDuDdS−ZŴN pNMNNpNdS.(15) From (15) it is clear, that since, at time t,uDand pNare a priori known prescribed values, the total potential energy is actually a function of only the contact displacements uC. We further modify (15) in order to keep the unknown contact displacements uCin the integrals on ŴCpart only, by utilizing the following reciprocity relations between the elastic solutions of PCand PN −ZŴN pNMNCuCdS=ZŴC uCMCNpNdS(16) as well as between the solutions of PCand PD, ZŴD uDMDCuCdS=ZŴC uCMCDuDdS. (17) A similar reciprocal relation holds true for the solutions of PN and PD, −ZŴN pNMNDuDdS=ZŴD uDMDNpNdS. (18) Notice a somewhat surprising presence of the negative sign in (16) and (18). It should be mentioned that (16), (17) and (18) hold only approximately in the case of numerical solution of the pertinent problems. Substituting the relations (16) and (17) into the expression (15) gives 5(t,uC)=ZŴC uC1 2MCCuC+MCDuD+MCNpNdS+1 2ZŴD uDMDDuDdS +1 2ZŴD uDMDNpNdS−1 2ZŴN pNMNDuDdS−1 2ZŴN pNMNNpNdS(19) The present problem can be defined now in an explicit form as a quadratic programming problem. By observing here that the last four integrals in (19) are constant with respect to uC, we may omit them in the minimization procedure since they do not have influence on the minimizer of the total potential energy. This leads to the following simplified expression including only the variable part of the total potential energy functional (for the sake of simplicity the notation of this functional is kept the same): 5(t,uC)=1 2ZŴC uCMCCuCdS |{z } A(uC) +ZŴC uCMCDuD(t)+MCNpN(t)dS |{z } -b(t,uC) (20) After a suitable discretization of uCby means of a linear combination of standard vector shape-functions uC(x)= MC X i=1 φi(x)ξi(21) Frontiers in Built Environment | www.frontiersin.org 5November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact where MCis the number of displacement degrees-of-freedom (DOFs), typically nodal values of displacement components, on the contact zone ŴC. The expression in (20) may be rewritten as an algebraic representation of the corresponding quadratic function 5(t,ξ)=1 2ξTAξ−ξTb(t), (22) where ξis a vector of the nodal displacements uCon ŴC, and Ais a kind of the boundary stiffness matrix associated to ŴC. Recall that (20), (21) and (22) refer to a single domain. Then, the total potential energy for the assemblage of Nbodies is defined as follows 5= N X i=1 5i(t,ui C). (23) This leads to the fact that the algebraic representation of the total potential energy 5for the assemblage in of Nbodies in (23) adopts the same form as in (22), however, with the matrix Adefined as a block diagonal matrix formed by the individual matrices Aifor subdomains i, A=block-diag A1,...,AN, (24) and similarly the vector bdefined as a block vector by formed by the individual vectors bifor i bT=n(b1)T,..., (bN)To. (25) 4. BOUNDARY ELEMENT METHOD As can be seen, the possibilities for a numerical implementation of the above framework are quite wide, the key issue is the approximation of the operators Mi CC,Mi CD and Mi CN for each subdomain i. Several numerical methods, such as the finite element method (FEM) or boundary element method (BEM), may be utilized. Since the adjacent bodies represented by i (i=, 1, ...,N), possibly in contact, are linearly elastic, collocation BEM or symmetric Galerkin BEM provide an intrinsic advantage herein. In the present implementation, the collocation BEM has been adopted. The BEM is closely related to the map between the prescribed boundary conditions in displacements or tractions and the unknown boundary displacements or tractions. In pure Dirichlet and Neumann boundary-value problems (BVPs), these maps are called Steklov-Poincaré and Poincaré-Steklov maps, respectively, and BEM can be considered as an approach to discretize these maps. In the present computational procedure, the role of the BEM analysis, applied to each subdomain iseparately (which, in fact, makes this problem very suitable for parallel computers), is to solve the corresponding BVPs on each i. For this goal, we numerically solve the Somigliana displacement identity (see París and Cañas, 1997; Aliabadi, 2002), ci ml(y)ui m(y)+ZŴi −ui m(x)Ti ml(x,y) dSx=ZŴipi m(x)Ui ml(x,y) dSx, (26) where y∈Ŵiand ui m(x) and pi m(x) denote the m-component of displacement and traction vectors, respectively. The weakly singular integral kernel Ui ml(x,y), two-point tensor field, given by the Kelvin fundamental solution (free-space Green’s function) represents the displacement at xin the m-direction originated by a unit point force at yin the l-direction in the unbounded linear-elastic medium defined by C(i). The strongly singular integral kernel Ti ml(x,y), two-point tensor field, represents the corresponding tractions at xin the m-direction. The coefficienttensor ci ml(y) of the free-term is a function of the local geometry of the boundary Ŵiat y, and may be evaluated by a closed analytical formula for isotropic elastic solids (see Mantiˇ c, 1993). The symbol R−in (26) stands for the Cauchy principal value of an integral. Consider a discretization of the boundary Ŵiby a boundary element mesh, which is also used to define a suitable discretization of boundary displacements ui(x) and tractions pi(x) by interpolations of their nodal values along Ŵidenoted as uiand pi, respectively. By imposing (collocating) the Somigliana identity (26) at all boundary nodes (called collocation points) we set the BEM system of linear equations for Ŵi, typically written as Hiui=Gipi(see París and Cañas, 1997; Aliabadi, 2002). The solution of this system provides the unknown nodal values of boundary displacements and tractions. In the present computational implementation of BEM, straight elements with piecewise linear continuous interpolation for displacements and piecewise linear (possibly discontinuous) interpolation for tractions are adopted, respectively. On the ŴC part of the boundaries, usually standard continuous interpolation for tractions is used. In the context of BEM, operators Mi IJ, of (12), with I,J=C,D and Ndo not need to be explicitly constructed, since what we actually need are the resultant tractions or displacements for each sub-problem, as also shown in Figure 4. Thus, we need to solve these three sub-problems for each subdomain iby BEM. 5. EXPLICIT COMPUTATION OF THE MATRIX A AND VECTOR B The matrix Ain (22) has the physical meaning of a boundary stiffness matrix for a subdomain. Several established techniques may be utilized for its computation, a survey of such techniques based on the BEM may be found in Hartmann (1981) and Tullberg and Bolteus (1982). Hereinafter a new methodology for computing Ais introduced. This methodology is consistent with the present energetic approach and leads directly to a symmetric form of A. From (20), the quadratic form Aand linear form b write as A(ξ)=1 2ξTAξ=1 2 MC X i,j=1 aijξiξj, (27) b(ξ)=ξTb= MC X i=1 biξi, (28) with MCbeing the number of displacement DOFs on the contact zone ŴC. Frontiers in Built Environment | www.frontiersin.org 6November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact Choosing ξ(i)defined as ξ(i) i=1, ξ(i) j=0, ∀j6= i, (29) then (27) simplifies to A(ξ(i))=1 2aii →aii = MC X k,l=1 ξ(i) kξ(i) lZŴC φkMCCφldS. (30) Solving the above problem PCby BEM sequentially for all i= 1, ...,MC, and by utilizing the first integral of (20), we may compute all the diagonal components aii of matrix A. Furthermore, in order to compute the non-diagonal components aij we choose ξ(i,j)defined as ξ(i,j) i=ξ(i,j) j=1, ξ(i,j) m=0, ∀m6= i,j, (31) then, similarly to (30), A(ξ(i,j))=1 2aii +ajj +aij +aji=1 2aii +1 2ajj +aij (32) =1 2 MC X k,l=1 ξ(i,j) kξ(i,j) lZŴC φkMCCφldS. (33) Having computed the diagonal elements aii we may compute from the last equation also the non-diagonal elements, considering aij =aji, by aij =1 2  MC X k,l=1 ξ(i,j) kξ(i,j) lZŴC φkMCCφldS−aii −ajj . (34) Similarly, solving the above problems PDand PNby BEM, the computation of vector bis carried out by using the first set of vectors ξ(i)and computed by the second integral in (20), and its components are given as b(ξ(i))=bi→bi= − MC X k=1 ξ(i) kZŴC φkMCDuD+MCNpNdS. A more direct relation of the matrix Ato a respective stiffness matrix could be understood for the case of a single domain where the whole boundary Ŵcoincides with ŴC(see also section 10). Although by the procedure described above a direct symmetric stiffness matrix may be constructed, this matrix was found to lack the properties of equilibrium and rigid motion representation, which are satisfied only approximately and although after a mesh refinement their fulfillment is improved nevertheless they are not fully satisfied. A further detailed study of the properties of this matrix is given in section 10. Additionally in section 9, some alternative modified procedure in order to construct a nonsymmetric matrix Atogether with the reasoning of its necessity is proposed and discussed. 6. THE CONTACT DESCRIPTION VIA INTERFACE ELEMENTS Usually, a unilateral frictionless contact law is defined by a geometric condition of non-penetration (also referred to as compatibility of displacements), a static condition of no tension and no friction and equilibrium of forces (see Renaud and Feng, 2003). In the context of small displacements deformation, we may write the non-penetration condition in the normal direction associated to the surfaces under contact (see below). The exact portion of the boundaries iand jwhich will really contact, the so-called active contact zone, is a priori unknown, for this reason initially we consider Ŵi jand Ŵj ias potential contact boundaries, which should include the points of the active contact zone (see Figure 1). FIGURE 5 | (A) Interface elements constructed on an intermediate surface between two domains under possible contact, and (B) main and secondary elastic domains, corresponding nodes and interface elements. Frontiers in Built Environment | www.frontiersin.org 7November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact The establishment of relations that impose the nonpenetration also referred to as Signorini condition, is one of the most crucial parts of the present formulation. Some of the BEM formulations for contact problems enforce compatibility of displacements on one of the domains, while on the other the equilibrium is enforced (cf. Blázquez et al., 1998b; Graciani et al., 2005). These equilibrium conditions might be given in weak sense instead of their strong form (see Blázquez et al., 1998a). In the present methodology only compatibility conditions are enforced between the bodies in potential contact, by establishing the non-penetration condition which include the evaluation of the displacement jump [u] and gap gnalong ŴC, while the equilibrium of contact tractions arises through the enforcement of the principle of minimum of the total potential energy. The interconnection of the subdomains as well as the consideration of Signorini kinematical conditions is succeeded by the intermediate surface discretized by interface elements, as shown in Figure 5A. Let us consider an intermediate surface between two bodies. The gap gn, at this surface may be considered as finite or zero valued. For each node of the interface elements we define a local reference system and each relative displacement may have an opening (normal) and sliding (tangential) component. In the case of contact problems of two deformable bodies where only small changes in the geometry are assumed and also the mesh of each elastic domain is a conforming one, it is possible to incorporate the contact constraints on a purely nodal basis. However, for the general case of nodes being arbitrarily distributed along the potential contact zone between two bodies, which can occur, e.g., when automatic meshing is used for these bodies, further consideration must be taken into account in the definition of Signorini contact conditions. Such techniques have been developed in the case of BEM in Blázquez et al. (1998b),Graciani et al. (2005),Blázquez et al. (2006), Graciani et al. (2009), and Vodiˇ cka (2000). Nevertheless, the necessity of non-conforming algorithms may also occur in small displacements contact problems and conforming discretizations as shown in Blázquez and París (2009). The non-penetration conditions for the case of non-conforming mesh are considered in the normal direction of the intermediate surface at interface nodes. Each interface node, related to a pair of the adjacent BEM node and boundary point, the latter may coincide with a BEM node or just be located on a BEM element as depicted in Figure 5B. For the latter case, the displacements of the boundary point are computed by the respective displacements of the nodes of this BEM element. The Signorini conditions, that also take into account a possible initial gap (or a finite thickness of the adhesive layer) are imposed on the normal displacements of those pairs. We will write here the conditions taken into account for two different cases referring to Figure 5B. In the first case, the BEM node is denoted as im(main) and the adjacent boundary point coincides with another BEM node denoted as is(secondary), and the following condition is taken into account: um(im)−us(is)n≥ −gn(i). (35) In the next case, the BEM node is denoted as lm(main) and the adjacent boundary point, denoted as ls(secondary), does not coincide with any BEM node (see Figure 5B), then the following condition is taken into account: um(lm)−us(ls)n≥ −gn(l)⇒ um(lm)−Lln Lkn us(ks)−Lkl Lkn us(ns)n≥ −gn(l). (36) The entire set of such equations may be given in a more compact algebraic form in terms of new variables fi=cijuj, (37) where fiis given by the ith constraint expression, cij is the socalled constraint matrix and ujare the components of normal displacements on the potential contact zone, which are the only variables included in the minimization procedure for the frictionless contact as will be shown in the following section. The inclusion of the tangential components of displacements on the contact zone in the case of adhesive frictionless contact is straightforward and will be omitted here for the sake of brevity. 7. MINIMIZATION OF TOTAL POTENTIAL ENERGY As shown in section 2, the contact problem can be formulated as a minimization problem of the total potential energy functional defined in (4). In section 3 appropriate boundary integral forms of this functional were derived in order to generate a pure boundary value problem, as in (10) or (11), which was further manipulated in order to derive a quadratic functional in (20). In section 1 it was assumed, that the problem is a quasi-static one (neglecting inertial and viscous effects), introducing by this way a (pseudo)time variable tinto the stored energy functional in (4). Finally, in section 6 it was shown that in order to define the minimization problem as a simple bound constrained one it is advantageous to work with variables fidefined by (37). Making a time discretization by adopting, for simplicity only and without loss of generality, an equidistant partition of [0, T] with a fixed time-step τ > 0, assuming T/τ ∈N, the minimization of the total potential energy functional defined in (4) leads to the following (formally written) recursive minimization problem: minimize Fk(f)=5(kτ,u(f)) subject to f≥ − g(38) to be solved successively for k=1, ..., T/τ. Here, fis the vector defined by Equation (37) and gthe vector that contains the gap values. For the numerical resolution of the above minimization problem there are two possible approaches, always taking into account the transformation in (37), the first is to directly minimize the objective function of type given in (11). In this case we may use a general algorithm for bound-constraint optimization problems, such as the L-BFGS-B by Zhou et al. Frontiers in Built Environment | www.frontiersin.org 8November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact (1997). However, it is usually advantageous to solve the problem as a quadratic one and for this reason we utilize (23) or its algebraic form as in (22). A variety of algorithms for quadratic programming may be found in Dostál (2009). Notice that, using the algebraic form of (22), the problem may be stated as a Linear Complementarity Problem (LCP), cf. (Gakwaya et al., 1992; Stavroulakis and Antes, 1997), and may be solved, for example, by Lemke’s algorithm. For the present EC-BEM computer implementation, we use a modified conjugate gradient method which does not need an explicit form of matrix Aand vector bin (22) and directly minimizes (23). Its pseudocode is given in Table 1. Nevertheless, also the conventional conjugate gradient method is reproduced in Table 2 following Dostál (2009). A further discussion on differences between these two variants of the conjugate gradient algorithm is given in section 9. Regardless of the minimization algorithm we use, except if it would be a derivative free algorithm, we always need to calculate TABLE 1 | (M1) Conjugate gradient method (CG) adapted for EC–BEM. Given a quadratic form A(x) and a linear form b(x) with x∈Rn Step 0. {Initialization.} Choose x0∈Rn, set g0=A′ x(x0)−b′ x(x0), p1=g0,k=1 Step 1. {Conjugate gradient loop.} while kgk−1k> ǫ Solve BEM for displacements pkon ŴC αk=k gk−1k2/A(pk) xk=xk−1−αkpk gk=gk−1−αkA′ x(pk) βk=k gkk2/kgk−1k2 pk=pk−1−βkpk k=k+1 end while Step 2. {Return the solution.} ˆ x=xk TABLE 2 | (M2) Conjugate gradient method (CG). Given a matrix A and a vector bwith x∈Rn Step 0. {Initialization.} Choose x0∈Rn, set g0=Ax0−b,p1=g0,k=1 Step 1. {Conjugate gradient loop.} while kgk−1k> ǫ αk=k gk−1k2/1 2(pk)TA pk xk=xk−1−αkpk gk=gk−1−αkA pk βk=k gkk2/kgk−1k2 pk=pk−1−βkpk k=k+1 end while Step 2. {Return the solution.} ˆ x=xk the partial derivatives of the potential energy 5. An alternative to compute these derivatives would be to use finite differences but such an approach was proved to be extremely time consuming. For this reason, we set up an analytical computation of derivatives of 5. By applying the chain rule of differentiation we obtain ∂5 ∂fi =∂5 ∂uj ∂uj ∂fi , (39) and by inverting (37) ui=rijfj, (40) where rij =(c−1)ij, which substituted into (39) leads to ∂5 ∂fi =∂5 ∂uj rji. (41) Having established the quadratic formula for (23), the derivatives in (41) can be easily defined. In section 8 we give a more general treatment on how to define the derivatives of the total potential energy (11) with respect to a displacement component. 8. DERIVATIVES OF THE TOTAL POTENTIAL ENERGY WITH RESPECT TO DISPLACEMENTS When the minimization problem (38) is formulated as a quadratic programming problem, derivatives of 5, for a single domain, with respect to the degrees-of-freedom (DOFs) of displacements uCon the potential contact zone can be easily computed. However, we might also treat the minimization problem as a general one. In that case we might need the derivatives of 5as given in (11). The total potential energy functional 5 from (11), for a single domain, is a function of udefined on Ŵ, but it also contains the tractions p, which are functions of displacements, i.e., p(u), and for this reason the derivatives with respect to ucannot be computed in a straightforward manner. The displacement and traction variables, according to the BEM formulation, are discretized for x∈Ŵas u(x)=X i φi(x)ui, (42a) p(x)=X i ψi(x)pi, (42b) where φi(x) and ψi(x) are the vector shape-functions associated to DOFs idefining the boundary values of displacements and tractions, respectively. The partial derivative of the total potential energy with respect to a DOF uiis the functional derivative with the test function being the vector shape-function φiassociated to this DOF ∂5(u) ∂ui =lim δui→0 5u+φi(x)δui−5(u) δui , (43) Frontiers in Built Environment | www.frontiersin.org 9November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 14 | A half elastic cylinder on an elastic rectangle. TABLE 5 | Geometry, elastic properties and imposed displacements. Model A B uD(mm) 37.5 37.5 E1=E2(GPa) 70 70 ν1=ν20.35 0.35 L1(m) 1.5 1.5 L2(m) h2sin φh2sin φ h1(m) 0.75 0.75 h2(m) 0.75 0.75 φ(◦) 18 45 TABLE 6 | Values for the variables nidefining the different meshes. Mesh 1 2 3 n113 26 94 n26 12 100 nh17 13 47 nh27 13 47 nφ6 12 100 n90◦−φ8 16 53 nL1−L211 22 61 is considered on the potential contact zone. Due to the fact that this is an advancing contact problem, we are interested in tracking the evolution of the load, which is a non-linear function of the prescribed displacements. In the numerical solution the vertical displacement is introduced by increasing 10 times this displacement in a monotonic way leading to the value uDshown in Table 5 at the step 10. The problem is solved using the EC-BEM code described above with the mesh 1 (see Table 6), where n1,n2,nh1and nφ correspond to the number of elements on the boundary part defined by L1,L2,h1 and φrespectively. n90◦−φis the number of FIGURE 15 | Deformed shape predicted by EC-BEM for Mesh 1 of Table 6. elements along the arc of the half cylinder which is not a potential contact zone and nL1−L2is the number of elements at the two parts of the upper edge of the rectangle which are not potential contact zones. The deformed shape predicted by this code is represented in Figure 15. Consequently, as expected, the vertical displacement imposed on the half circle results in its indentation onto the rectangle. A part of the potential contact zone is in active contact. In a similar way to the example described in section 11.1.1 and for comparison purposes, the problem is solved by the BEM code implementing a classic algorithm of contact by Blázquez et al. (1998a, 2006). The results of the deformed shape around the contact zone for both meshes in Table 6 are represented in Figure 16. Results obtained by the codes using different meshes (see Table 6) and distinct lengths of the potential contact zone (see Table 5). The aim of this figure is to compare both algorithms, and, for that reason, results obtained with the same mesh are grouped. Macroscopically both solutions for EC-BEM and PM-BEM show a very good agreement. However, it is interesting to see in the detailed plot of solutions for mesh 1, plotted in Figure 17, that results for EC-BEM, even for this coarse mesh, show a very good agreement with the Reference solution. A comparison of tractions along the potential contact zone predicted by EC-BEM, PM-BEM and the reference solution is represented in Figure 18. In this figure, unlike those described previously, results are grouped by the method, in order to evaluate the influence of the mesh on the results for each method. Notice that both codes present more accurate results for finer meshes, thus confirming the expected solution convergence with mesh refinement. The resultant force along the upper part ŴDof the half cylinder can be estimated from the results of tractions predicted by both methods. In particular, we focus on the vertical component Fyof this force because it is the only one which produces work since horizontal displacements are zero by the Frontiers in Built Environment | www.frontiersin.org 16 November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 16 | Comparison of the deformed boundaries around the contact zone predicted by EC-BEM, for two different meshes (Mesh 1 and 2 of Table 6 in the top and bottom plot, respectively) with the Reference solution by PM-BEM (using Mesh 3 of Table 6). boundary conditions. This value is calculated by the integration of stresses over ŴD, Fy=ZŴD σyydx(48) where σyy represent the normal tractions at the nodes. The integral is calculated by a standard quadrature. The values of the resultant force Fy, calculated for the 10 steps by both codes, are represented in Figure 19 as a function of the vertical displacement imposed. Both codes predict that stiffness is increasing with the vertical displacement, as typical for the Hertz contact problems. The reason is that the contact zone between both solids is enlarging, with a consequence that the global system becomes stiffer. FIGURE 17 | Detailed comparison of the deformed boundaries around the contact zone predicted by EC-BEM and PM-BEM for coarse meshes (Meshes 1 and 2 of Table 6) together with the Reference solution by PM-BEM (using Mesh 3 of Table 6). 11.2. Numerical Examples to Test the Adhesive Contact Implementation This section aims to compare the results obtained by the energetic approach and commercial FEM code (ANSYS, 2010) in the case of the elastic contact described in section 2. To this effect, the problem of receding contact studied above is solved again, this time with the contact conditions of the adhesive elastic contact. The problem represented in Figure 9 is studied with the properties of model A, see Table 3, and the mesh 1, see Table 4. The stiffness parameters, assumed for this example, are κn=150 GPa/m and κt=0. Results from the EC-BEM computations for displacements and tractions along the boundary between both solids are presented in Figure 20. As expected, zones with tensions correspond to positive relative displacements (external zones for example) and zones with compressions correspond to vanishing relative displacements (zones adjacent to the previous ones). In this case, a strongly complex zone is predicted at the center of the boundary, which will be studied below. For comparison purposes, the problem is solved using the FEM code (ANSYS, 2010) setting the contact properties in the following way: •Contact behavior: No separation (always) •Contact: Surface to surface •Contact algorithm: Penalty method •Contact detection: At nodal points (normal from contact nodes) •Opening contact stiffness (KFOP): κn Results for tractions in the solids obtained by ANSYS and the EC-BEM code are compared in Figure 21. In the upper figure we can observe that the agreement along the four zones Frontiers in Built Environment | www.frontiersin.org 17 November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 18 | Comparison of the normal tractions computed by EC-BEM and PM-BEM along the interface for different meshes. described above is accurate. It is interesting to focus on the zone in the middle of the boundary where tractions looks like vanishing at the scale of the upper figure. Lower figure shows the complexity of this zone, which is composed by two lateral regions where both solids are separated and a central region where both solids are under contact (compression). The agreement between both results is reasonable in spite of the very small scale compared to the tractions obtained on the whole boundary. A subsequent analysis is carried out by multiplying the displacement uDby 100 without modifying any other parameters in both codes. In this case, the results from the energetic approach presented here remain qualitatively similar, whereas the results from ANSYS present some problems to match the complex central zone. To the knowledge of authors, this could be solved by modifying some of the values of the keyoptions which configure the contact algorithm of ANSYS. This shows the robustness of the energetic approach compared to other contact algorithms, where convergence depends on the value of certain non-physical parameters of the numerical model. FIGURE 19 | Force Fypredicted by both methods as a function of the imposed vertical displacement uD. FIGURE 20 | Normal tractions and normal relative displacements along the potential contact zone in the example of Figure 9 with the adhesive elastic contact computed by EC-BEM, and with uD=0.035 mm. 12. CONCLUDING REMARKS A novel approach for contact problems of elastic bodies was presented. This method is based on the energy principles expressed by boundary integrals. The solution of contact problems is obtained by the minimization of the total potential energy. Appropriate algorithms to solve the quadratic problems were presented. The cases of conforming as well as non-conforming meshes for contact zone were addressed. Numerical examples, showing the suitability of the present framework to solve typical contact problems were examined, Frontiers in Built Environment | www.frontiersin.org 18 November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact FIGURE 21 | Comparison of the normal tractions computed by EC-BEM and ANSYS for the adhesive elastic contact, with uD=0.035 mm. and also successfully compared to other well established BEM methodologies as well as a commercial FEM code. The proposed framework can be easily incorporated to existing BEM or FEM codes. Possible extension to the case of elastodynamics, where instead of the principle of virtual work we may use the principle of virtual power in order to compute the total energy (kinetic plus elastic), using only boundary values, of tractions and velocities, is possible. For this it might be convenient to utilize a formulation such as that introduced in Panagiotopoulos and Manolis (2010, 2011). Finally, the present framework is found to be very useful in the case where also dissipative mechanisms exist on the boundaries or common interfaces of elastic bodies. In such cases we may apply energetic approaches to solve the corresponding non-linear problems as was shown in previous publications by some of the present authors (e.g., Panagiotopoulos et al., 2013). AUTHOR CONTRIBUTIONS All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication. FUNDING The authors acknowledge the support by the Junta de Andalucía and European Social Fund (Proyecto de Excelencia TEP-4051), the Spanish Ministry of Science and Innovation (Project MAT2009-14022), the Spanish Ministry of Economy and Competitiveness and European Regional Development Fund (Projects MAT2012-37387 and MAT2015-71036-P), and Spanish Ministry of Education (FPU grant AP2009-3968). Furthermore, CP acknowledges the financial support of the Stavros Niarchos Foundation within the framework of the project ARCHERS (Advancing Young Researchers’ Human Capital in Cutting Edge Technologies in the Preservation of Cultural Heritage and the Tackling of Societal Challenges). ACKNOWLEDGMENTS The authors would like to express their gratitude to Prof. Tomáš Roubíˇ cek of the Charles University of Prague, as well as to Dr. Roman Vodiˇ cka from the Technical University of Košice, for fruitful discussions. REFERENCES Aliabadi, M. H. (2002). The Boundary Element Method: Applications in Solids and Structures. Chichester, UK: Wiley. ANSYS (2010). ANSYS Mechanical APDL Contact Technology Guide Release 13.0. Canonsburg, PA: Ansys, Inc. Antes, H., and Panagiotopoulos, P. D. (1992). The Boundary Integral Approach to Static and Dynamics Contact Problems. Equality and Inequality Methods, Vol. 108 of International Series of Numerical Mathematics. Basel: Birkhauser Verlag. Aour, B., Rahmani, O., and Nait-Abdelaziz, M. (2007). A coupled FEM/BEM approach and its accuracy for solving crack problems in fracture mechanics. Int. J. Solids Struct. 44, 2523–2539. doi: 10.1016/j.ijsolstr.2006. 08.001 Axelsson, O. (1980). Conjugate gradient type methods for unsymmetric and inconsistent systems of linear equations. Linear Algebra Applic. 29, 1–16. doi: 10.1016/0024-3795(80)90226-8 Blázquez, A., Mantiˇ c, V., and París, F. (2006). Application of BEM to generalized plane problems for anisotropic elastic materials in presence of contact. Eng. Anal. Bound. Elements 30, 489–502. doi: 10.1016/j.enganabound.2005.07.006 Blázquez, A., and París, F. (2009). On the necessity of non-conforming algorithms for “small displacement” contact problems and conforming discretizations by BEM. Eng. Anal. Bound. Elements 33, 184–190. doi: 10.1016/j.enganabound.2008.05.003 Blázquez, A., París, F., and Cañas, J. (1998a). Interpretation of the problems found in applying contact conditions in node-to-point schemes with boundary element non-conforming discretizations. Eng. Anal. Bound. Elements 21, 361– 375. doi: 10.1016/S0955-7997(98)00024-1 Frontiers in Built Environment | www.frontiersin.org 19 November 2018 | Volume 4 | Article 56 Panagiotopoulos et al. BEM Energetic Approach to Contact Blázquez, A., París, F., and Mantiˇ c, V. (1998b). BEM solution of twodimensional contact problems by weak application of contact conditions with non-conforming discretizations. Int. J. Solids Struct. 35, 3259–3278. doi: 10.1016/S0020-7683(98)00016-X de Paula, F. A., and Telles, J. (1989). A comparison between point collocation and Galerkin for stiffness matrices obtained by boundary elements. Eng. Anal. Bound. Elements 6, 123–128. doi: 10.1016/0955-7997(89)90025-8 Dostál, Z. (2009). Optimal Quadratic Programming Algorithms: With Applications to Variational Inequalities, Vol. 23 of Springer Optimization and Its Applications. New York, NY: Springer Science+Business Media. Eck, C., Jarušek, J., and Krbec, M. (2005). Unilateral Contact Problems. Variational Methods and Existence Theorems. Boca Raton, FL: Chapman & Hall; CRC; Taylor & Francis Group. Fichera, G. (1964). Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambigue condizioni al contorno. Mem. Accad. Naz. Lincei VIII, 91–140. Gakwaya, A., Lambert, D., and Cardou, A. (1992). A boundary element and mathematical programming approach for frictional contact problems. Comput. Struct. 42, 341–353. doi: 10.1016/0045-7949(92)90030-4 Graciani, E., Mantiˇ c, V., París, F., and Blázquez, A. (2005). Weak formulation of axi-symmetric frictionless contact problems with boundary elements. application to interface cracks. Comput. Struct. 83, 836–855. doi: 10.1016/j.compstruc.2004.09.011 Graciani, E., Mantiˇ c, V., París, F., and Varna, J. (2009). Numerical analysis of debond propagation in the single fibre fragmentation test. Compos. Sci. Technol. 69, 2514–2520. doi: 10.1016/j.compscitech.2009.07.006 Gurtin, M. E. (1972). The Linear Theory of Elasticity, Vol. VIa/2 of Handbuch der Physik. Berling: Springer-Verlag. Hartmann, F. (1981). The derivation of stiffness matrices from integral equations. Appl. Math. Model. 5, 355–365. doi: 10.1016/S0307-904X(81)80058-3 Hertz, H. (1896). Miscellaneous Papers on the Contact of Elastic Solids. Trans. D. E. Jones (London: McMillan). Hestenes, M. R., and Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand. 49, 409–436. doi: 10.6028/jres.049.044 Johnson, K. L. (1985). Contact Mechanics. Cambridge: Cambridge University Press. Kalker, J., and Randen, Y. V. (1972). A minimum principle for frictionless elastic contact with application to non-Hertzian half-space contact problems. J. Eng. Math. 6, 193–206. doi: 10.1007/BF01535102 Katsikadelis, J., and Kokkinos, F. (1993). Analysis of composite shear walls with interface separation, friction and slip using bem. Int. J. Solids Struct. 30, 1825–1848. doi: 10.1016/0020-7683(93)90236-Z Lazaridis, P. P., and Panagiotopoulos, P. D. (1987). Boundary variational “principles” for inequality structural analysis problems and numerical applications. Comput. Struct. 25, 35–49. doi: 10.1016/0045-7949(87)90216-1 Mantiˇ c, V. (1993). A new formula for the C-matrix in the Somigliana identity. J. Elastic. 33, 191–201. doi: 10.1007/BF00043247 Méchain-Renaud, C., and Cimetière, A. (2000). Influence of a curvature discontinuity on the pressure distribution in two-dimensional frictionless contact problems. Int. J. Eng. Sci. 38, 197–205. doi: 10.1016/S0020-7225(99)00018-X Panagiotopoulos, C. G. (2017). Symplegma: Java Components for Computational Mechanics and Numerical Modelling. Available online at: http://symplegma. org/ Panagiotopoulos, C. G., and Manolis, G. D. (2010). Velocity-based reciprocal theorems in elastodynamics and BIEM implementation issues. Arch. Appl. Mechan. 80, 1429–1447. doi: 10.1007/s00419-009-0376-0 Panagiotopoulos, C. G., and Manolis, G. D. (2011). Three-dimensional BEM for transient elastodynamics based on the velocity reciprocal theorem. Eng. Anal. Bound. Elements 35, 507–516. doi: 10.1016/j.enganabound.2010.09.002 Panagiotopoulos, C. G., Mantiˇ c, V., and Roubíˇ cek, T. (2013). BEM solution of delamination problems using an interface damage and plasticity model. Computat. Mechan. 51, 505–521. doi: 10.1007/s00466-012-0 826-3 Panagiotopoulos, P. D., and Lazaridis, P. P. (1987). Boundary minimum principles for the unilateral contact problems. Int. J. Solids Struct. 23, 1465–1484. doi: 10.1016/0020-7683(87)90064-3 París, F., and Blázquez, A. (1994). On the displacement and load scaling techniques in contact problems using bem. Bound. Elements Commun. 5, 15–17. París, F., Blázquez, A., and Cañas, J. (1995). Contact problems with nonconforming discretizations using boundary element method. Comput. Struct. 57, 829–839. doi: 10.1016/0045-7949(95)92007-5 París, F. and Cañas, J. (1997). Boundary Element Method, Fundamentals and Applications. Oxford: Oxford University Press. París, F., and Garrido, J. A. (1989). An incremental procedure for friction contact problems with the boundary element method. Eng. Anal. Bound. Elements 6, 202–213. doi: 10.1016/0955-7997(89)90019-2 Renaud, C., and Feng, Z.-Q. (2003). BEM and FEM analysis of Signorini contact problems with friction. Comput. Mechan. 31, 390–399. doi: 10.1007/s00466-003-0441-4 Roubíˇ cek, T., Mantiˇ c, V., and Panagiotopoulos, C. G. (2013). A quasistatic mixed-mode delamination model. Disc. Cont. Dynam. Syst. S 6, 591–610. doi: 10.3934/dcdss.2013.6.591 Saad, Y. and Schultz, M. H. (1985). Conjugate gradient-like algorithms for solving nonsymmetric linear systems. Math. Comput. 44, 417–424. doi: 10.1090/S0025-5718-1985-0777273-9 Signorini, A. (1959). Questioni di elasticità non linearizzata e semilinearizzata. Rendiconti di Matematica e delle sue Applicazioni 18, 95–139. Stavroulakis, G. E., and Antes, H. (1997). Nondestructive elastostatic identification of unilateral cracks through BEM and networks. Comput. Mechan. 20, 439–451. doi: 10.1007/s004660050264 Tullberg, O., and Bolteus, L. (1982). “A critical study of different boundary element stiffness matrices,” in Boundary Element Methods in Engineering, ed C. A. Brebbia (Berling: Springer–Verlag), 621–635. Vodiˇ cka, R. (2000). The first and the second–kind boundary integral equation systems for solution of frictionless contact problems. Eng. Anal. Bound. Elements 24, 407–426. doi: 10.1016/S0955-7997(00)00 012-6 Wriggers, P. (2006). Computational Contact Mechanics, 2nd Edn. Berling: Springer. Wriggers, P., and Panagiotopoulos, P. D. (1972). New Developments in Contact Problems, Vol. 384 of CISM courses and lectures. NewYork, NY: Springer Wien. Yastrebov, V. A. (2013). Numerical Methods in Contact Mechanics. London, UK: ISTE Ltd and Wiley. Zhou, C., Byrd, R. H., Lu, P., and Nocedal, J. (1997). Algorithm 778: L-BFGS- B: Fortran subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Soft. 23, 550–560. doi: 10.1145/279232.279236 Conflict of Interest Statement: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Copyright © 2018 Panagiotopoulos, Mantiˇ c, García and Graciani. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. Frontiers in Built Environment | www.frontiersin.org 20 November 2018 | Volume 4 | Article 56