Among the wide spectrum of engineering applications, the class of problems involving contact interactions between solids are complex both with regard to their mathematical description and numerical treatment. The transfer of mechanical forces, thermal or electrical conduction, tectonic plate movements are merely a few examples of the ubiquity of contact in nature and engineering applications. The vast majority of phenomena occurring at the contact interface, such as wear, adhesion, lubrication and fretting, as well as mass and energy transfer, determine to the greater extent the service life of engineering components. This lays a strong emphasis on the importance to gain a fine understanding of these mechanisms for the accurate analysis and timely prediction of failure. In the last decade, a class of surface-to-surface contact discretizations, such as mortar method, coupled with appropriate treatment of inequality contact constraints has been well established and proved its ability to efficiently treat contact problems puso_mortar_2004 ; fischer_mortar_2006 ; popp_mortar_2012 . However, the numerical treatment of the contact and related phenomena remains particularly challenging when it involves the evolution of complex surfaces. These complexities can involve cumbersome remeshing algorithms and field transfer procedures, as for instance those encountered in the context of wear. In addition, construction of adequate finite element meshes near contact interfaces and stability of the contact formulations are necessary ingredients to ensure the overall efficiency, accuracy and robustness of the numerical procedures.
In this work, first we present a coherent framework for mortar frictional contact based on the monolithic augmented Lagrangian method for both normal and tangential contact components alart_mixed_1991 ; cavalieri2015numerical . We then extend this framework for a two dimensional unified framework, called MorteX, capable to treat frictional contact problems between a real and a virtual (intra-mesh) surfaces. It combines the features of Generalized/Extended Finite element methods (GFEM/X-FEM) dolbow_extended_1999 ; belytschko_review_2009 and mortar methods for solving the contact problems belgacem_mortar_1998 ; mcdevitt_mortar-finite_2000 ; puso_mortar_2004 ; popp_finite_2009 . This unified framework was elaborated in akula_tying_paper for mesh tying problems and is extended here to handle frictional contact problems. Within the proposed framework, in a discretized setting (see Fig. 1) the contact constraints are imposed between a “real” surface (explicitly represented by edges of finite elements) and an embedded virtual surfaces (internal surface geometrically non-conformal with the underlying discretization).
The X-FEM is an enrichment method based on the partition of unity (PUM) for discontinuous fields melenk_partition_1996 ; babuska_partition_1997 . In this framework, discontinuities are not explicitly modeled, and as a result, the mesh is not required to conform to these discontinuities. In X-FEM, enrichment functions are added to the finite element interpolation using the framework of PUM, to account for non-smooth behavior without compromising on the optimal convergence ferte_interface_2014 . The X-FEM methods are extensively used in applications such as fracture mechanics, shock wave front and oxidation front propagation, and other applications involving discontinuities both strong and weak sukumar_arbitrary_2000 ; sukumar_modeling_2001 ; diez_stable_2013 ; gross_extended_2007 ; ji_hybrid_2002 , and in particular for the contact between crack lips dolbow_extended_2001 ; khoei_contact_2006 ; ribeaucourt_new_2007 ; elguedj_mixed_2007 ; liu_contact_2008 ; gravouil_stabilized_2011 ; mueller-hoeppe_crack_2012 . The ability of the X-FEM to describe geometric features (e.g. inclusions and voids) independently of the underlying mesh, will be exploited in the proposed unified framework to describe complex and evolving surfaces.
As opposed to the node-to-node discretisation which is limited to small sliding and the node-to-surface method which is inaccurate, the surface-to-surface based mortar methods are stable and very accurate in treating contact problems under finite sliding conditions el_abbasi_stability_2001 . The mortar methods provide us with a comprehensive framework to address the limitations of incompatible interface discretizations. They are a subclass of domain decomposition methods (DDM), that are tailored for non-conformal spatial interface discretizations bernardi_new_1994 , and were originally introduced for spectral elements belgacem_spectral_1994 ; bernardi_coupling_1990
. The coupling and tying of different physical models, discretization schemes, and/or non-matching discretizations along interfaces between the domains can be ensured by mortar methods. The mathematical optimality and applicability of the mortar methods in spectral and finite element frameworks were studied extensively for elliptic problems inbernardi_coupling_1990 ; belgacem_spectral_1994 ; wohlmuth_discretization_2001 . The mortar methods were brought into the contact realm with the initial contributions from Belgacem et al., Hild, Laursen et al. belgacem_mortar_1998 ; hild_numerical_2000 ; mcdevitt_mortar-finite_2000 . Subsequent works of Fischer and Wriggers fischer_frictionless_2005 ; fischer_mortar_2006 , Puso and Laursen puso_mortar_2004 ; puso_mortar_2004-1 have extended the application of mortar methods to the class of non-linear problems. The work of Popp, Gietterle, et al gitterle_finite_2010 ; popp_mortar_2012 dealt with the class of dual Lagrangian resolution schemes wohlmuth_mortar_2000 .
In the present methodology, we use the monolithic augmented Lagrange multiplier resolution scheme to treat the frictional contact problems. It represents a so-called mixed formulation brezzi_mixed_2012 . The choice of Lagrange multipliers’ functional space strongly affects the convergence rate and can lead to loss of accuracy in the interfacial tractions. These difficulties arise from a locking type phenomena reported for mixed variational formulations as a result of non-satisfaction of Ladyzhenskaya-Babuška-Brezzi (LBB) babuska_finite_1973 ; brezzi_mixed_2012 . This locking, present in standard FEM with Lagrange multipliers on the boundary barbosa_finite_1991 , but also in the mixed X-FEM framework when Dirichlet boundary conditions are applied along embedded surfaces, was extensively studied in fernandez-mendez_imposing_2004 ; moes_imposing_2006 ; bechet_stable_2009 ; burman_fictitious_2010 ; hautefeuille_robust_2012 ; ramos_new_2015 . However, to the best of our knowledge, the manifestation of mesh locking effects when imposing contact constraints using Lagrange multipliers, was never reported. The manifestation of the locking effect will be illustrated in the present work for contact problems both in the context of standard mortar and MorteX methods. In addition, stabilization strategies initially proposed for the mesh tying problems in the context of MorteX method akula_tying_paper , are adapted here for contact applications.
The paper is organized as follows. In Section 2 and 3, we present the strong and weak continuum problem description of two-body frictional contact. The classical spatial mortar interface discretization scheme for the frictional contact problem is derived in Section 4. In Section 6 this framework is extended to account for embedded contact surfaces, and referred as MorteX framework. It will explicitly detail the features of the X-FEM and adaptation of the mortar which are needed to construct a coherent computational scheme. In Section 8, few examples including the contact patch test are considered to illustrate the methods performance for both frictionless and frictional contact problem settings. In Section 9, we give the concluding remarks and present the prospective works.
2 Problem definition
The case of two deformable bodies in contact is considered, for simplicity the derivations are made in infinitesimal deformation formalism, but is readily extandable to the finite strain framework. Fig. 2 shows two statements of this contact problem: (a) is a classical formulation where the solids come in contact over their outer surfaces. In (b) an equivalent, from continuum point of view, setting is shown, whereby contact boundary is still an outer boundary of , while the homologue deformable body is now and its contact boundary is embeded in the domain , whose part is discarded from the computation. In the current and next two sections, we will state the problem and suggest the methods to handle the configuration shown in Fig. 2(a), whereas in Section 6, we will extend these methods to deal with the problem shown in Fig. 2(b). The domains , , have non-intersecting Dirichlet boundaries , Neumann boundaries and contact boundaries , i.e. . The standard boundary value problem (BVP) is described by the balance of linear momentum along with the imposed boundary conditions:
is the Cauchy stress tensor. Under the small deformation assumption the strain-displacement and constitutive relations are:
where and are the Cauchy strain tensor and the fourth-order elasticity tensor describing Hooke’s law, respectively. represents the density of body forces, is the unit outward normal to and , are the prescribed tractions and displacements on and , respectively. The contact between the bodies introduces constraints which can be treated as additional configuration-dependent traction (), which have to be only compressive (in absence of adhesion) and ensure non-penetration of solids (4).
The gap function defines the signed distance between the two surfaces and is the fundamental measure that determines the status of contact. Here, we define the gap vector as a vector from a pointon the surface to its closest fprojection along the normal onto the surface (the choice of these surfaces is arbitrary):
The normal gap function is thus obtained as a dot product of the gap vector (6) and the outward unit normal
The MorteX method is, at this stage, elaborated for two-dimensional problems only, which allows some simplifications in case of frictional contact, which is also presented in 2D. The tangential relative sliding velocity is given as in puso_mortar_2004 :
where is the unit tangential vector to at point . In order to preserve the consistency of units, hereinafter the incremental slip will be used instead of the slip velocity , where is an arbitrary time increment. This replacement does not change the physical sense of all following equations. With these notations at hand, the contact traction can be decomposed into normal and tangential contributions:
Along the contact interface, the classical Hertz-Signorini-Moreau conditions kikuchi_contact_1988 also known as KKT conditions in the theory of optimization hestenes_multiplier_1969 ; powell_algorithms_1978-1 have to be satisfied. These conditions are formulated using the normal gap and the contact pressure as:
The tangential friction is described using the Coulomb friction law. The frictional constraints are formulated using the tangential slip increment , the tangential traction and the contact pressure . The incremental slip vanishes when the tangential traction is below the frictional threshold , non-zero incremental slip is possible when . These conditions can be also formulated as KKT conditions:
where is the Coulomb’s coefficient of friction.
3 Weak form
The Eqs. (1)-(3) denote the strong form of the standard solid mechanics BVP, and Eq. (4), along with the contact conditions (10)-(11) represents the contact part of the BVP. The weak form of the contact problem reduces to a variation inequality and the associated displacements should be selected in a constrained functional space, which has to include the non-penetration condition fichera_boundary_1973 ; kikuchi_contact_1988 . It is a constrained minimization problem. However, using the augmented Lagrangian method permits us to convert this problem into a fully unconstrained one, with displacements and virtual displacements (test functions) selected in unconstrained functional spaces: and the test function space definitions:
where denotes the first order Sobolev space. The virtual work corresponding to the structural part is given by:
where are the change in internal energy and the virtual work of forces, respectively.
Regarding the contribution of the virtual work associated with the contact problem, we adapt the monolithic augmented Lagrangian method (ALM) framework of Alart and Curnier alart_mixed_1991 to the mortar framework. The ALM is a mixed dual formulation that introduces Lagrange multipliers as dual variables along with the primal displacement variables. The fields of Lagrange multipliers and are respectively equivalent to the contact pressure and the tangential friction shear introduced in Eq. (9). They and their virtual variations are selected from the appropriate functional space , namely fractional Sobolev space of order . The augmented Lagrangian functionals and for the normal and frictional contact respectively are given by the following expressions pietrzak_continuum_1997 ; pietrzak_large_1999 :
where are the augmented normal and tangential Lagrange multipliers, and the augmented pressure, respectively:
where and are the normal and tangential augmentation parameters, respectively. The augmented pressure determining the frictional threshold is replaced by the corresponding augmented Lagrange multiplier , however is not subjected to variation, therefore a different notation is used to highlight this subtle difference. Integrating the functionals (15) and (16) over the contact surface integrates the contact constraints in the Lagrangian:
where is the potential energy of the system.
The problem of contrained optimization transforms into the min-max or saddle point problem for the Lagrangian, whose solution is equivalent to the solution of the variational inequality optimization problem on a non-convex domain and its boundary hestenes_multiplier_1969 ; powell_algorithms_1978-1 ; pietrzak_continuum_1997 . The min-max problem is solved by searching the stationary point minimizing the Lagrangian with respect to primal (displacement) field and maximizing it with respect to the dual field (Lagrange multipliers). The variation of the structural part, which contains only primal variables is given in (14), whereas the variation of the contact contribution is given by:
The contact virtual work (19) is expressed as a summation of the contact contribution to the virtual work resulting from the variations of the positions (20), the weak contribution of normal contact constraints (21) from the variation of the normal Lagrange multiplier and the weak contribution of tangential contact constraints (22) from the variation of the tangential Lagrange multiplier. Note that the variations and derivatives of piece-wise smooth potentials and can be seen as sub-derivatives pietrzak_continuum_1997 . We recall that the variation of is not required. Based on the three possible contact statuses (i.e. stick, slip or no contact), the contact surface is divided into three distinct sub-surfaces:
The frictional contact contribution to the virtual work can be split as follows:
4 Mortar interface discretization
Within the mortar discretized framework, the contacting surfaces are typically classified into mortar and non-mortar sides. In the following, the superscript “” refers to the mortar side of the interface and “” to the non-mortar side. The choice of mortar and non-mortar is rather arbitrary: however, in general the choice of the finer meshed surface as mortar side is known to result in higher accuracy. All the contact related integral evaluations will be carried out on the mortar side of the interface. In the considered two dimensional case, the edges of potentially contacting surfaces are parametrized by . In the current configuration the interpolations of mortar and non-mortar surface edges are given by:
and are the interpolation functions of mortar and non-mortar sides, respectively; whereas and is the number of nodes per segment at mortar and non-mortar sides, respectively. Note that we use isoparametric linear elements, which implies that the normal remains constant along the edge. Moreover, that in the current framework, for the sake of simplicity, we omitted averaging of normals at nodes and their interpolation, which was elaborated in yang_two_2005 ; popp_mortar_2012 . Hereinafter, Einstein summation over repetitive index is used. We assign scalar Lagrange multipliers and for the normal and tangential directions to the nodes of the mortar side. The Lagrange multiplier vector is interpolated functions :
where is the number of nodes used for interpolation over every edge, where L can be less than or equal to M.
Every mortar contact element is created between parts of mortar and non-mortar segments which are projected one on the other (see Fig. 3). The segment-to-segment projection is described in detail in wriggers_computational_2006 ; yang_two_2005 ; popp_finite_2009 and will not be reproduced here. Each contact element consists of nodes, each with primal DoFs in 2D, and of or dual DoFs associated with Lagrange multipliers in case of frictionless or frictional contact, respectively.
4.1 Discrete integral kinematic quantities
As a first step towards the elaboration of the contact part of the residual vector, the discretisations (24) and (25) are inserted in the continuous weak form (23). First, we focus on the derivation of a single term of this integral, which enables us to introduce coincise notations and simplifications to be used for the other terms. It is obtained for every active contact element by integrating over a part or an entire mortar edge :
where is the unit outward normal of the mortar edge. Note that this equality holds if the parametrization ensures collinearity between the vector in brackets and the normal vector. Therefore, depends on positions of all the nodes both mortar and non-mortar forming the contact element. This is not explicitly reflected in the equation but should be kept in mind. The hat-notations, which denoted projections in (6,8), are therefore omitted hereinafter.
where and are the so-called mortar
integrals popp_finite_2009 , which are defined as
is the integral normal gap. The variation of the integral normal gap is given by:
By analogy, the contribution of the tangential contact term can be introduced as:
where is the integral slip increment given by puso_mortar_2004 :
where the current time step and the previous time step come in an interplay. The variation of the incremental slip, again following puso_mortar_2004 is given by:
4.2 Evaluation of mortar integrals
The evaluation of the mortar side integrals is straightforward, as it only involves the product of shape functions from the mortar side (28). This is in contrast with the evaluation of the mortar integral which combines shape functions from both the mortar and non-mortar sides. Its evaluation requires a mapping between them (Fig. 3). For this purpose, non-mortar quantities are projected onto the mortar side of the interface using the mortar segment normal vector. The non-mortar nodes are projected onto the mortar side segment along the mortar segment normal . The local coordinates of the projections are found by solving
The extremities of the integration domain are defined either by mortar or non-mortar edge nodes, details of this procedure are provides in yang_two_2005 and not repeated here. The segment of the mortar edge over which the evaluation of mortar integrals is carried out is parametrized by . The projection coordinates , which determine the limits of the integration domain, are mapped on the segment parameterization:
Fig. 3 illustrates how the integration domain is determined between two edges. To evaluate the integrals using Gauss quadratures, the mortar-side Gauss points are projected along mortar segment normal onto the non-mortar side and the corresponding local coordinates are determined by:
The mortar matrices, evaluated with Gauss quadratures, take the following form:
where, as previously, , , and , where is number of Gauss integration points, is the normalized Jacobian
The factor reflects the fact that the integral is evaluated only over a part of the mortar edge.
4.3 Residual vector
Within the augmented Lagrangian formulation, every created contact element contributes to the virtual work of the system irrespectively of its contact status (active or inactive). This translates into a smoother energy potential in comparison to the standard method of Lagrange multipliers. For each contact element, the discretized form of the virtual work can be obtained for the three possible contact statuses, namely stick, slip and no contact as follows:
5 Extended finite element methods
In this section, we present some elements of the extended finite element method (X-FEM) dolbow_extended_1999 ; belytschko_review_2009 . In complement to the mortar method, X-FEM is a key ingredient required to construct a unified MorteX framework to treat contact problems along along surfaces embedded in the bulk mesh. By analogy with the MorteX method for tying akula_tying_paper , the mesh hosting the embedded surface will be referred to as the host mesh. Such an embedded (or virtual) surface splits the domain of the host solid into two parts. One part is discarded (i.e. considered as empty space) and the remaining effective part represents the solid, see Fig. 4. The virtual line runs through elements, which will be called blending elements. As such, these elements are also made of discarded and effective parts. The finite elements are then essentially partitioned into three distinct categories, namely the discarded, blending and standard elements. Here, X-FEM will be used to account for the fact that only the effective volume contributes to the internal virtual work.
The virtual surface of the host domain is treated as an internal discontinuity, the outer part of the volume is disregarded using the X-FEM framework. The X-FEM relies on enhancement of the FEM shape functions used to interpolate the displacement fields. The enrichment functions describing the field behavior are incorporated locally into the finite element approximation. This feature allows the resulting displacement to capture discontinuities. The subdivision of the host mesh is defined by indicator function (where is the spatial position vector in the reference configuration in domain ) sethian_level_1999 . The indicator function is non-zero only in the efficient part of the host domain :
The discontinuity surface can be seen as a level-set defined as follows:
In the considered framework, for the sake of simplicity, we will assume that the embedded surface is represented by a piece-wise linear line with possible kink points occurring inside host elements.
In practice, the enrichment of shape functions in case of void/inclusion problem can be simply replaced by a selective integration scheme sukumar_modeling_2001 . For the standard elements, there is no change in volume of integration and the discarded elements are simply excluded from the volume integration procedure. In order to obtain the effective volume of integration for each blending element, we perform the clipping of the blending elements by the embedded surface [Fig. 5] . The clipped element is virtually111Virtually stands here in the sense that the interpolation order and connectivity of nodes do not change. triangulated to perform the integration of the internal virtual work using a Gauss quadrature rule. In this framework, we use three newly initialized Gauss points per triangle which ensures a rather accurate integration in view of linear or bi-linear shape functions for triangular and quadrilateral elements, respectively. The clipping of a single element by an arbitrary embedded surface results in one or several various polygons222Hereinafter, we assume that all elements use first order interpolation, therefore all edges of elements are straight. It enables us to assume that an intersection or difference of elements can be always represented as one or several polygons. either convex or non-convex or even disjoint, which represent the effective volumes of integration. The tilde () notations are adopted for the quantities related to the effective part of the host mesh ().
6 MorteX framework
Here, the MorteX framework, which is presented in a separate paper akula_tying_paper in the context of mesh tying problems, is extended to treat frictional contact problems. Similar to the tying problem, it consists of two distinct procedures. The first procedure invokes the X-FEM method to account for the fact that only a portion of the body is contributing to the internal virtual work; it is called the effective volume, see on Fig. 4. The second procedure deals with the enforcement of the contact constraints between a “real” surface (explicitly represented by edges) and an embedded virtual surface , which is geometrically non-conformal with the finite element discretization and passing through the elements. The contact treatment is adjusted through few modifications of the classical mortar contact schemes introduced in Section 4.
6.1 MorteX interface discretization
Similar to the tying problem in the MorteX framework, the rationale behind the choice of mortar and non-mortar sides, holds for the contact problems as well akula_tying_paper : the real surface is always selected as the mortar surface. The interpolations of the nodal positions in the current configuration and the Lagrange multipliers along the mortar side remain the same as for the classical schemes, and are given by:
where as previously and is the number of nodes per mortar segment and those carrying Lagrange multipliers respectively. Vector represents the nodal positions of the mortar nodes forming the segment, and is the parametric coordinate of the mortar side. The parametrization on the non-mortar side requires more consideration. As shown in Fig. 6, the embedded surface is divided into straight segments, whose vertices can either be a clip intersection (i.e. the intersection between and the edge of a bulk element of the non-mortar side) or a kink point (a vertex of the discrete virtual surface that lies inside a blending element). Hereinafter, we will refer to these vertices as “points of interest”, without distinction between kink and clip intersection points; their notation will be equipped with a tilde . The shape of the embedded surface in the current configuration parametrized by between two points of interest, can be defined using two-dimensional shape functions inherited from the host element [see Fig. 6]:
Irrespective of the fact that the embedded surface is assumed piece-wise linear in the reference configuration, it can become piece-wise non-linear in the current configuration, therefore the displacement along every linear segment cannot be parametrized by linear shape functions. Note that here, in the MorteX contact set-up, the kink points represent vertices of the discretized surface , whereas they represented imprints of mortar side nodes in case of MorteX tying set-up akula_tying_paper .
6.1.1 MorteX contact element
A MorteX contact element is formed between a single mortar segment and a blending element, and it consists of nodes, from the mortar segment, and from blending element. In addition the contact element stores or Lagrange multipliers for the frictionless and frictional cases, respectively. As was shown in akula_tying_paper , to avoid mesh-locking (or over-constraining of the interface), which results in spurious oscillations, the number of Lagrange multipliers should be balanced by the number of blending elements, i.e. often it is reasonable to use one set per blending element, a detailed discussion can be found in Section 7.
6.1.2 Discrete integral kinematic quantities
The aforementioned selective integration scheme, used to accommodate for the presence of an embedded virtual surface, leads to changes in the discrete contact integral quantities: the integral gap vector (27), the integral normal gap (31) and the incremental slip (34). For a MorteX contact element, these nodal quantities are now evaluated along the interface formed by a pair of real and embedded segments (tilde-notations are used):
where , , , and are the mortar integrals in the MorteX framework, which will be defined in (58). Note that within the MorteX framework, the definitions of purely mortar side quantities, like the mortar segment normal and mortar side integral remain the same as in the classical mortar framework, whereas now presents a convolution of volumetric and surface shape functions of primal and dual quantities, respectively.
6.1.3 Discrete contact virtual work
The MorteX residual vector is larger than the standard mortar one as it involves the displacement DoFs from the bulk of the non-mortar side of the interface. The residuals for stick, slip and non-contact statuses are: