An interface-enriched generalized finite element formulation for locking-free coupling of non-conforming discretizations and contact

We propose an enriched finite element formulation to address the computational modeling of contact problems and the coupling of non-conforming discretizations in the small deformation setting. The displacement field is augmented by enriched terms that are associated with generalized degrees of freedom collocated along non-conforming interfaces or contact surfaces. The enrichment strategy effectively produces an enriched node-to-node discretization that can be used with any constraint enforcement criterion; this is demonstrated with both multi-point constraints and Lagrange multipliers, the latter in a generalized Newton implementation where both primal and Lagrange multiplier fields are updated simultaneously. We show that the node-to-node enrichment ensures continuity of the displacement field—without locking—in mesh coupling problems, and that tractions are transferred accurately at contact interfaces without the need for stabilization. We also show the formulation is stable with respect to the condition number of the stiffness matrix by using a simple Jacobi-like diagonal preconditioner.


Introduction
The computational modeling of problems in contact mechanics requires careful considerations in order to prevent interpenetration between contacting bodies and ensure a proper transfer of contact tractions.A related problem arises in the coupling of non-conforming finite element discretizations, where the presence of hanging nodes, if not handled properly, yields a discontinuous displacement field.In this work these problems are addressed by means of an enriched finite element method that naturally leads to a simple computer implementation and inherently avoids the emergence of locking due to an over-constrained interface.
The numerical analysis of many engineering problems requires the coupling of meshes of different components.These meshes are often non-conforming, i.e., the locations of their nodes do not coincide along the coupling interface, resulting in hanging nodes that call for special treatment.More importantly, the coupling procedure has to avoid over-constraining the coupling interface-to prevent locking-and ensure a proper transfer of tractions between subdomains.Similar issues are shared by engineering applications that involve contact between bodies as, for instance, forging processes, gear systems, and impact problems [1].While coupling non-matching meshes is enforced by means of equality constraints, contact problems use inequality constraints, which makes modeling contact notoriously challenging because of its intrinsic highly nonlinear nature.In general, numerical methods used to resolve contact problems are either incapable of accurately transferring tractions between contacting bodies or require a very intricate computer implementation.
Coupling of non-conforming interfaces and contact is numerically handled by enforcing constraints, to guarantee continuity for the former and avoid interpenetration while allowing both bodies to slide or detach for the latter.Although constraints can be enforced in several manners, their application often leads to locking due to an over-constrained interface.In the multi-point constraint (MPC) method, also known as master/slave approach, non-conformity is handled by constraining slave nodes on one side of the interface to those on the master side.This method, however, has been shown to be sensitive to the choice of master/slave surfaces insofar as guaranteeing continuity.A two-pass MPC approach, where both sides are taken as master and slave of each other [2,3], ensures continuity but suffers from locking due to an over-constrained interface.Dual space methods, such as finite element tearing and interconnecting (FETI) [4] and the mortar method [5,6], use a Lagrange multiplier field to enforce compatibility at the interface.In both methods, the Lagrange multiplier field must be selected carefully to avoid locking by satisfying the inf-sup condition, also known as the Ladyzhenskaya-Babuška-Brezzi (LBB) condition [7,8].A detailed account of several methodologies proposed in the literature to enforce contact constraints is given § 2.
Enriched finite element methods, whereby the primal field is enhanced or enriched by means of appropriate enrichment functions, provide an elegant approach to dealing with non-conforming meshes and contact.For instance, Haikal and Hjelmstad [3] proposed an enriched stabilized discontinuous Galerkin formulation for coupling non-conforming meshes that recovers accurate tractions and is devoid of locking.This was accomplished by modified finite element (FE) shape functions and the addition of enriched nodes, effectively recovering an NTN type of constraint enforcement.
The eXtended/Generalized Finite Element Method has also been applied to solve contact and non-conforming interface problems [9,10,11,12,13].Duarte et al. [9] modified the FE partition of unity by means of clustering and used X/GFEM to deal with non-matching interfaces.For contact problems, Khoei and Nikbakht [13] and Dolbow et al. [10] used X/GFEM to model frictional contact, where contact surfaces were treated as embedded strong discontinuities.Hirmand et al. [11] proposed an augmented Lagrangian-based stabilization technique to model frictional contact and thus obtain smooth contact tractions.Similarly, Akula et al. [12] used the augmented Lagrange method (ALM) within the mortar method to model complex contact surfaces as embedded discontinuities in a background mesh using X/GFEM; a special stabilization technique was required when contacting bodies have high contrast in stiffness and/or mesh density.Nevertheless, X/GFEM is in general unreliable insofar as obtaining accurate interface tractions or requires stabilization techniques.Aragón and Simone [14] reported the oscillatory nature of the traction field recovered by X/GFEM on a notched beam where a cohesive formulation was used to model perfect bonding.It has also been shown that X/GFEM cannot be used as an immerse boundary (or fictitious domain) method since recovered tractions in Dirichlet boundaries oscillate even when using the Barbosa-Hugues stabilization [15].The issue is caused by an over-constrained interface (or boundary), and thus the inf-sup condition is not satisfied [16].Putting aside the issue of oscillatory tractions, X/GFEM may be unstable with respect to the condition number of system matrices, a problem that has prompted recent efforts in pursuit of a stable GFEM (SGFEM) [17,18].Finally, X/GFEM may have poor accuracy in blending elements (elements containing both standard and enriched nodes), prescribing non-homogenous Dirichlet boundary conditions requires especial techniques, and the computer implementation is far from trivial [14].
Inspired by GFEM, another family of enriched FE formulations seeks to solve problems with discontinuities by placing enrichments to nodes created along the discontinuities-as opposed to X/GFEM, where enrichments are added to the nodes of the original mesh.The Interface-enriched Generalized Finite Element Method (IGFEM) [19] was first proposed to resolve weak discontinuities, i.e., situations where the gradient of the primal field is discontinuous, as in interface problems.The method later developed into the Hierarchical Interfaceenriched Finite Element Method (HIFEM) [20], whereby multiple interfaces within a single element are resolved via a hierarchical implementation of enrichment functions.The Discontinuity-Enriched Finite Element Method (DE-FEM) [14] later proposed a generalization of IGFEM/HIFEM for the treatment of both weak and strong discontinuities-e.g., fracture-in a unified formulation.Because enrichments are placed along discontinuities, this family of methods is devoid of many of the disadvantages of X/GFEM: the implementation is straightforward, the method is stable, i.e., the condition number grows with the problem size as in standard FEM [21], and there are no blending elements since enrichment functions are localized to cut elements and vanish at all original mesh nodes.The latter therefore keep the Kronecker-δ property, which significantly simplifies the enforcement of non-homogeneous essential boundary conditions [16,14].Contrary to most non-standard FEMs in use today, in IGFEM/HIFEM/DE-FEM Dirichlet boundary conditions can be enforced strongly as in standard FEM.This is also true even along non-matching Dirichlet boundaries, where HIFEM/DE-FEM is used as an immerse boundary procedure with smooth recovered traction fields [22,23].This is particularly interesting in the context of coupling non-conforming discretizations and contact and has largely inspired the present study.
In this paper we adopt the IGFEM paradigm for coupling non-conforming discretizations and solving contact problems.To this end, we place enriched nodes along non-conforming interfaces and contact surfaces, respectively, collocated at the locations of the non-conforming mesh nodes on the opposite surface.Their associated enriched functions are C 0 -continuous and local to the enriched elements on one side of the non-conforming interface/contact surface.Notably, since the shape functions of the original mesh are kept intact, the partition of unity property is therefore lost in enriched elements, i.e., at any point in the element the sum of shape and enrichment functions does not add to unity.The formulation yields an enriched NTN contact discretization for which is straightforward to enforce constraints by any procedure; here we explore the enriched procedure with both MPC and ALM.The proposed procedure is therefore a two-pass method and therefore shows no bias on the choice of master/slave surfaces.The procedure is applied to several examples with linearized kinematics.Results show that the method passes Taylor and Papodopoulos' contact patch test [24] and avoids locking due to an over-constrained interface; furthermore, traction transfer is fairly accurate without the need of stabilization techniques.This is advantageous over other methods such as the two-pass mortar method which, although passes the patch test and shows no locking, requires interpolation of the pressure field that results in a complex formulation and corresponding computer implementation.We also show the stability of the method insofar as the condition number of the stiffness matrix is concerned.Indeed, the condition number grows with mesh size h as O h −2 -i.e., at the same rate as that of the standard FEM-after the application of a simple diagonal preconditioner.

Previous work on contact discretizations
We limit the scope of this survey to contact between deformable bodies, i.e., deformable-deformable contact.Many contact discretization procedures have been proposed through the years.Node-to-node (NTN) [25] and the node-to-segment (NTS) [26,27,28] discretizations are widely used, where constraints are enforced between a node pair or a node-segment pair, respectively.Their applicability is, however, hindered by some intrinsic limitations: Even though NTN can transfer tractions accurately, it can only be used with conforming discretizations along master/slave surfaces and only for infinitesimal sliding.The NTS approach overcomes this restriction but a one-pass approach-where nodes of one contacting surface are constrained to the segments of the other surface-is biased insofar as the choice of master/slave surfaces.A one-pass NTS approach not only fails to prevent interpenetration at times, it has been shown not to pass the contact patch test [24], i.e., correct contact tractions cannot be recovered.Several works have attempted to address this issue.Papadopoulos and Taylor [2] proposed to enforce constraints in an average sense and integrate the contact pressure using Simpson's rule, thereby passing the patch test.Later, Zavarise and De Lorenzis [29] improved the one-pass NTS contact formulation to pass the contact patch test at the expense of losing symmetry in the contact contribution to the stiffness matrix.
A two-pass NTS approach could also be used to avoid interpenetration of contacting bodies, but this approach leads to an over-constrained contact interface and ultimately to locking.As a result, the traditional two-pass NTS contact discretization does not fulfill the LBB condition [30,2] and therefore no convergence can be attained with mesh refinement.Papadopoulos et al. [31] and Papadopoulos and Solberg [32] proposed a method that enforces continuity of tractions weakly, whereby nodes are divided into groups in which the gap is constrained and pressure continuity is enforced.Their approach avoids locking and passes the patch test when the same interpolation for geometry and traction field is used.The method was later extended to 3D by Jones and Papadopoulos [33].
The idea of enforcing traction continuity or constraints in a weak sense can be understood as a segmentto-segment (STS) method, first proposed by Simo et al. [34], where displacement field is interpolated with linear shape functions, the pressure field is defined as piecewise constant over contact segment, so that the local variables (gaps) are evaluated in an "average" sense.Another STS-based procedure by Zavarise and Wriggers [35] also enforces contact constraints in a weak sense, whereby local variables are evaluated at integration points and then contact contributions to the stiffness matrix and force vector are calculated based on the integration of these local variables.El-Abbasi and Bathe [36] proposed yet another STS method in which integration points are projected onto the opposite contact surface.Although their formulation is capable of handling both linear and quadratic elements in contact, the LBB condition is satisfied only when the pressure is interpolated with linear continuous functions.In addition, a "composite" integration rule that combines Gaussian and Newton-Cotes rules is required to pass the patch test.All these formulations and their corresponding computer implementations are more complex than those of NTS contact.
The mortar method was first used to handle domain decomposition problems [5] and then contact [37,38,39].The method also has single-pass and dual-pass versions.The traction field is defined on a mortar surface, which is one of the two contacting surfaces or an intermediate surface that is introduced, and contact constraints are enforced in a weak sense.The single pass mortar method has been shown to pass the contact patch test [40].This method was later extended to 3D and large deformations [40,41,42,43,44,45,46,47,48].For large deformation kinematics, gap functions interpolated by shape functions in the mortar surface and nodally-averaged normal vectors are used [44].To avoid non-smooth normals, special techniques are also employed to obtain smooth surfaces, for instance by means of Hermite functions [41] or by other averaging techniques [49].Results of this method are generally much more accurate than those of traditional NTS procedures-but at the same the computer implementation is more involved and requires more computational resources [50].Although some mortar-based methods satisfy the LBB condition, they are biased insofar as the choice of contact surfaces is concerned.Still, an appropriate choice of dual (Lagrange multiplier) shape functions should be made.Furthermore, the computer implementation is much more complex than standard NTS, which is also less demanding in terms of computational costs, especially in 3D.To avoid the contact surface bias issue, Solberg et al. [51] proposed a two-pass Mortar method with contact constraints enforced by means of Lagrange multipliers, thereby nodes are separated into active and inactive sets; however, a penalty-based stabilization is needed to avoid traction oscillations, and this method is not guaranteed to work for arbitrary geometries, particularly in 3D.Recently, this method was also extended to solve plasticity and self-contact by Puso and Solberg [52].Park et al. [53] and Rebel et al. [54] proposed a method, different from the mortar method, where nodes from two contact surfaces are also projected into an intermediate surface, and an independent variable is defined to describe the motion of this surface.By enforcing the traction continuity weakly, constant stress can be transferred exactly for contact patch test.This method was also extended to handle frictional contact by Gonzalez [55].On the basis of the mortar method, a dual Mortar method was proposed by Flemisch et al. [6] to solve contact problems with curved interfaces in 2D, where discontinuous shape functions for dual variables are used because they are more stable and ensure more accurate solutions.This method was then extended to 3D contact problems with large deformation by Popp et al. [56] and Popp and Wall [57], and extended to frictional contact by Gitterle et al. [58].In general, the dual Mortar method offers the possibility of condensing the multipliers out of the system matrix (another advantage of using dual shape functions), resulting in a symmetric positive definite matrix.In their overview, Popp and Wall [57] highlight several advantages of the dual mortar method for contact over the traditional NTS approach.They show that, although smooth interpolations of the contact surfaces are needed (for instance using NURBS), especially for large deformation contact, solving problems with dual mortar can be more efficient than with the traditional single-pass mortar procedure.
Because the smoothness of contact surfaces influences the transfer of tractions greatly, several surface smoothing techniques have been proposed over the years: Belytschko et al. [59] suggested to compute gap functions based on a least square fit of the original non-smooth surface, obtaining a smooth signed distance function which reduces traction oscillations; Padmanabhan and Laursen [60] and Sauer [61] proposed to use Hermite functions to interpolate contact surface and solution fields; El-Abbasi et al. [62] suggested to use cubic splines; Krstulovic-Opara et al. [63] proposed a Bézier formulation; Sauer [64] used high-order Lagrange shape functions; and Stadler et al. [65] used a NURBS-based formulation.These smoothing techniques aim at reducing traction oscillations and to improve convergence.
Contact constraints have also been enforced recently in the context of isogeometric analysis (IGA), whereby CAD features are used directly in the numerical analysis.Lu [66] and De Lorenzis et al. [67] proposed IGA formulations for frictionless and frictional contact, respectively; a similar approach to NTS was established in IGA, the so-called knot-to-surface, where the main difference with respect to the traditional NTS lies in that the knot is used instead of mesh nodes and NURBS are used to describe contact surfaces.This method was later extended to 3D [68,69,70], and a mortar-based IGA procedure was found to be more accurate than the standard knot-to-surface method.Corbett and Sauer [71,72] proposed to use a NURBS interpolation along contact surface, while the rest is discretized using linear elements.This modification makes the approach more efficient.For contact in large deformation, accurate results are still generally hard to obtain without mesh refinement.Dimitri et al. [73,74] proposed to use T-splines to refine the discretization.Similar to the dual mortar method with FEM, a dual mortar method in IGA was also proposed by Seitz et al. [75] and offers better efficiency than the mortar-based IGA contact formulation [76].An isogeometric collocation method was also used to solve contact problems [77,78], whereby contact forces are regarded as Neumann boundary conditions, and contact constraints are enforced as Dirichlet conditions.Since these conditions are fulfilled strongly, quadrature of Neumann BCs is eliminated and efficiency is improved greatly.Recently, Duong and Sauer [79] and Duong et al. [80] proposed an IGA contact formulation based on the surface potential method [81] and equipped with the refined boundary quadrature method [82].In this procedure the "two-half-pass" method [83] is used, i.e., only tractions on the slave surface are considered in each pass.XFEM is also used in this method to capture weak/strong discontinuities around contact boundaries [80].With a smoothing technique for post-processing, results are much more accurate than those obtained with traditional two-pass methods.However, IGA-based methods still suffer from locking, with the exception of the two-half-pass method, and the procedure requires a completely different discretization based on NURBS-which may not be straightforward to implement in general displacement-based FEM codes.
For most methods above, it is still not trivial to obtain accurate contact tractions for complex problems due to numerical artifacts [29].Even for mortar methods traction oscillations are generated for curved contact interfaces-which are relatively small and concentrated near the contact edge with the correct choice of mortar surface [12].These numerical errors do not necessarily reduce with mesh refinement.Some STS and mortar based methods [32,33,51,53,84,52] are quite accurate in the definition of contact tractions, pass the contact patch test, and avoid locking.Their formulations and implementations are, however, much more complicated than those of traditional NTN or NTS procedures.
There are also procedures that aim at transforming the contact problem into an equivalent NTN contact constraint enforcement.For instance, the virtual element method (VEM) has recently been explored to solve problems in contact mechanics [85].In VEM non-conforming contacting nodes are projected onto the opposite contact surface, and new nodes are generated at those locations, transforming the problem into a VEM-conforming mesh for which NTN contact constraints can be enforced straightforwardly.Besides, adding contact constraints to a general VEM code is straightforward.Results show that VEM passes the patch test.However, a penalty-based stabilization is needed to avoid traction oscillations.This method was later developed to solve frictional contact in large deformation kinematics by Wriggers and Rust [86], and solve contact with curved contact surfaces by Aldakheel et al. [87].The coupling between two domains Ω 2 and Ω 3 (perfect bonding) is also illustrated.Traction t1 and t3 are applied on boundaries Γ t 1 and Γ t 3 , and boundaries Γ u 1 and Γ u 3 are subject to a prescribed displacement.Outer surface normals are denoted n i , i = 1 . . .3. The insets show non-matching meshes at the contacting interface, with original standard FEM nodes shown in black, and enriched nodes, on the opposite surface, shown as red circles.
Enriched formulations have also been used to convert non-conforming contact discretizations into an equivalent NTN enforcement.For instance, the enriched discontinuous Galerkin formulation for coupling non-conforming discretizations by Haikal and Hjelmstad [88] was later extended to solve frictionless contact with finite deformation kinematics [3].Masud et al. [89] further combined the ideas put forward by Haikal et al. [88,3] with Nitsche's method in a variational multiscale framework that could be used not only for coupling non-conforming meshes but also for solving frictional contact problems.The methodology proposed in this paper follows the enrichment strategy proposed by Haikal and Hjelmstad [3] in which, instead of modifying the shape functions of elements in contact, we keep the standard basis intact.Both strategies could be understood as a form of h-hierarchical refinement along contact surfaces and thus share some similarities with the non-hierarchical h-refinement approach used to handle non-conforming meshes in 3D proposed by Jiao and Heath [90,91].Their approach is based on the concept of common refinement of two meshes which is defined as "the intersections of the elements of the input meshes".The new discretization contains therefore nodes that can be found in both non-matching meshes.Our approach follows a similar strategy in that all standard nodes from one contact surface are projected to the other and vice versa.
3 Problem description and formulation

Governing equations
and their outer normals n i .Domains Ω 1 and Ω 2 are in frictionless contact along the surface Γ 12 = Γ 1 ∩ Γ 2 ≡ Γ c , and domains Ω 2 and Ω 3 are perfectly bonded along the interface Γ 23 = Γ 2 ∩ Γ 3 ≡ Γ g .Such interface could be physical (e.g., the interface between two different materials), numerical (e.g., a non-conforming discretization with hanging nodes), or a combination thereof.Boundaries Γ u i and Γ t i denote regions with prescribed Dirichlet and Neumann boundary conditions, respectively.These regions are disjoint, i.e., Γ u i ∩ Γ t i = ∅.The displacement (primal) field is denoted by u and the traction (dual) field by t.These are composed by considering the fields within all domains Ω i .We therefore denote u i the restriction of u to the ith domain, i.e., u i = u| Ωi .The same notation applies to the traction field t i and to other subscripted quantities.
We are interested in solving the elastostatics frictionless contact boundary value problem whose strong form is expressed as: Given body force b i : Ω i → R d , prescribed displacement ūi : Γ u i → R d , and traction ti : with boundary conditions interface conditions for perfect bonding and contact conditions In Equation ( 1), σ i = D i : ε i is the Cauchy stress tensor which obeys Hooke's law, D i is a fourth-order constitutive tensor, and ε i = 1 2 (∇u i + ∇u i ) the linearized strain tensor (small deformation theory).For the perfect bonding between domains Ω 2 and Ω 3 , equal displacements and tractions are enforced through interface Γ g .Contact between domains Ω 1 and Ω 2 is enforced through the classical Hertz-Signorini-Moreau conditions, also known as Karush-Kuhn-Tucker conditions in the theory of optimization [1]; these consist of a contact pressure p n , which can only be in compression, and a gap function g n , which ensures that the contact surfaces cannot penetrate one other.Finally, Equation (8) ensures that pressure vanishes when the gap between the two contact surfaces is open, and the gap vanishes when they are in contact.
Equations ( 1)-( 8) can be seen as the solution to a constrained optimization problem on the functional given by where the first term represents the potential energy in all domains.For the ith domain, the potential energy is The other two terms in (9) are associated with constraints at the perfectly bonded interface Γ g and at the contact interface Γ c , enforced with Lagrange multipliers λ g and λ c , respectively; their form varies depending on how the constraints are enforced.In this work we use the standard Lagrange multiplier method for the perfectly bonded interface, and an augmented Lagrangian type of enforcement for contact (which has been shown to outperform both Lagrange and penalty methods).Their forms are therefore respectively written as where • denotes the Macaulay brackets, λ n = λ c • n the normal component of the Lagrange multiplier (where we could use either n 1 or n 2 , as at this stage the choice is arbitrary), and λn the augmented Lagrange multiplier given by λn = λ n + n g n = p n , with n the augmentation (penalty) parameter in the normal direction.

Variational formulation
The solution of the boundary value problem makes the optimization problem (9) stationary with respect to variations of all arguments.These three equations are obtained by taking the directional derivative along the directions of the variations: for all admissible variations δu, δλ g , and δλ c .For the variations of the displacement, we define the vector-valued function space where H 1 (Ω i ) denotes the first-order Sobolev space on Ω i .The primal field is chosen from the set Lagrange multipliers λ and their variations δλ are taken from a fractional Sobolev space, i.e., λ a ∈ Λ (Γ a ) ≡ H −1/2 (Γ a ) d , with a = g, c.
Multiple-point constraints (MPCs) are also investigated in this work to enforce both perfect bonding and contact.In such case, the enforcement works by writing simple constraint equations that describe the relationship between slave degrees of freedom (DOFs) as a function of the master DOFs.As there is no need for Lagrange multipliers, the rightmost two terms in ( 9) disappear.
Next we discuss the discretization of Equations ( 13) through (15).We start with the enriched formulation used in elements where constraints are to be enforced.

The finite-dimensional interface-enriched generalized finite element formulation
It is not straightforward to couple meshes or simulate contact when the discretizations of the domains are non-matching, i.e., when there are hanging nodes in the former or contacting nodes do not occupy the same location in space in the latter (no node-to-node contact).To alleviate the burden associated with the lack of mesh conformity, we employ an enrichment scheme inspired by the Interface-enriched Generalized Finite Element Method (IGFEM) [19,21], whereby enriched nodes are created along material interfaces to resolve weak discontinuities (those in the field gradient).Consequently, all domains are discretized into non-overlapping finite elements such that Ω h i = ∪ i e i , e i ∩ e j = ∅ (∀i = j), with the entire computational domain given by In other to solve numerically the finite-dimensional form of the equations in §3.2, the trial solution and the weight function are chosen from the interface-enriched generalized finite element space In (18) the first term represents the standard finite element component and the second term the enrichment.In the former, ι h is the index set of all standard nodes in the mesh, N i is the ith standard Lagrange interpolation function associated with degrees of freedom u i (which physically represent the nodal displacement at standard node x i ).In the enrichment term, ι e is the index set of enriched nodes, ψ i is the enrichment function associated with enriched DOFs α i .Finally, s i is a scaling factor that is used to improve the condition number of the system matrix after assembly.This factor, which was studied thoroughly in a recent publication [21], is required for a robust implementation that handles interfaces getting arbitrarily close to standard mesh nodes.The enrichment function ψ i is constructed with the aid of Lagrange shape functions in integration elements (see enrichment ψ i corresponding to enriched node x ⊥ i in Figure 2), which not only form the support ψ i , but, as the name implies, are also used for the numerical quadrature of the local stiffness and force arrays.Contrary to the original IGFEM formulation, where the support of an enrichment comprises integration elements at both sides of a material interface, here the support comprises only integration elements on one side of the non-conforming interface or contact surface.In our implementation, a computational geometric engine creates enriched nodes along the non-conforming interfaces or along contact surfaces, so that each enriched node corresponds to a standard mesh node on the surface of the opposite domain.As shown in Figure 1, the location of mismatching mesh nodes can be directly determined along the interface Γ g in mesh coupling problems, whereas for contact problems the closest node projection [1,92] to an element edge is used to determine the location of enriched nodes at either side of contact surfaces.

Stiffness matrix contributions
In this section we derive the discrete expressions for the stiffness contributions given by the first term in (13).The assembly of the local stiffness matrix k i and force vector f i for finite elements that do not contain enriched nodes follows standard procedures.For enriched integration elements, and with reference to [19,14], the local arrays can be expressed as Figure 2: Enrichment function ψ i associated with enriched node x ⊥ i , which coincides with a non-conforming standard node x i .The support of ψ i comprises integration elements only in Ω 1 ; represents an enriched node and a standard mesh node. with where D is the constitutive matrix, and B u and B α are the strain-displacement matrices of shape and enrichment functions, respectively [22,93].Finally, the global stiffness matrix K and global force vector F are obtained by standard assembly of their local counterparts.Denoting by A the standard finite element assembly operator, global arrays are therefore which, by making explicit the contributions of standard and enriched components, can be expressed as

Constraint enforcement
The enriched space alone does not ensure continuity across non-conforming interfaces and does not avoid interpenetration of contact surfaces.To properly resolve the field at mesh coupling and contact interfaces, constraints are imposed between enriched and standard DOFs.Although the use of discrete MPCs and the augmented Lagrange method (ALM) is well established, their application in conjunction with this enriched framework is different and both methods are therefore thoroughly discussed in this section.

Multiple-point constraint method
Conventionally, multiple-point constraints (MPCs) are used in a "master and slave" situation to enforce some sort of compatibility between nodes, e.g., for enforcing continuity along an interface or to enforce periodicity at opposite sides of a domain; an MPC can be simply regarded as an approach to create a "tie" among nodes.This method can be used in a single-pass or in a two-pass approach, which means that either side is selected as the master surface, or that either side serves as the master to its opposite side, respectively.In both approaches, the constraint is expressed as a linear combination of displacement vectors as in where u s is the displacement of the slave node x s , ι m is the index set of the slave's master nodes, and g is the gap vector.In mesh coupling problems, as the objective is to ensure continuity across an interface, the homogeneous MPC is used, i.e., g = 0.This constraint method is simple and straightforward to implement.However, as already discussed in the introduction, continuity cannot be ensured in a single-pass approach, and the system Figure 3: Schematic for the enriched discretization with ALM.The locations of enriched nodes are detected using the projections of mismatching mesh nodes with the normal vector (e.g., n + i for enriched node x + i ).The symbols represent enriched nodes and the standard mesh nodes.may become over-constrained in a two-pass approach [3].In the following, we provide an alternative formulation of MPCs combined with the enriched approach outlined in the previous section for handling mesh coupling and contact problems without slipping and separation.In the context of IGFEM, MPCs have been employed in immersed domain problems [22] and in the enforcement of Bloch-Floquet periodic boundary conditions [94].Note, however, that for general contact problems it is preferred to adopt ALM (this will be explained in detail in the next section).For the cases where MPCs work, the objective is to ensure C 0 -continuity across non-conforming interfaces or contact surfaces in a given direction.
For a standard hanging node x i , an enriched node x ⊥ i is created at the same location although placed in the adjacent element (see Figure 2).The constraint equation enforces that the displacement of standard and enriched nodes, u i and u ⊥ i , respectively, have to be equal.Referring back to (18), this condition is expressed as where ι ⊥ h ⊂ ι h is the subset of mesh nodes with standard shape functions whose supports intersect the support of the enrichment function.Mathematically, if the support of a standard shape function is defined as ω i = { x| N i (x) = 0}, and similarly for an enriched function (26) we use ψ i (x ⊥ i ) = 1.A similar constraint is enforced for every hanging standard node and enriched node pair at both sides of the non-conforming interface, making the proposed procedure actually correspond to a two-pass approach.However, when compared to just using Equation (25), enriched DOFs render the interface response less stiff.These constraints can be written in matrix form for all enriched nodes in the system as

Note that in Equation
where U = u 1 . . .α 1 . . . is the vector containing all standard and enriched DOFs, Ū is the vector of independent (unconstrained) DOFs, and T is a transformation matrix storing coefficients N i and s i in Equation (26).Following standard procedures for MPCs [22], the unconstrained system KU = F , where K and F were given in (24), is transformed into K Ū = F , where K = T KT , F = T F , and the transformation matrix T ties both standard and enriched nodes.Notice that only a few additional enriched DOFs need to be added, and applying the MPC method between master and slave DOFs is straightforward.Numerical examples for coupling non-conforming discretizations will be shown in § 5 to illustrate the accuracy and robustness of the approach.

Augmented Lagrange method
For general contact problems where relative slip and separation occur, it is not convenient to use MPCs to enforce constraints because the contact status needs to be detected during the nonlinear iterative contact step.The augmented Lagrange method can be regarded as a combination of penalty and Lagrange multiplier methods, whereby the inequality-constrained minimization contact problem is transformed into an unconstrained saddle point problem.Compared with pure penalty and Lagrange multiplier methods, this approach decreases the conditioning number of the stiffness matrix and is able to satisfy constraints exactly with finite penalty parameters [95].In this section we obtain the discrete form corresponding to the weak form ( 13)- (15).When considering Lagrange multipliers, we seek a solution to the system K Û = F with where Û = u α λ is the vector of unknowns in terms of standard DOFs, enriched DOFs, and Lagrange multipliers.The solution is found incrementally by making use of a generalized Newton loop [96,97], i.e., K∆ Û = ∆ F , where K is mostly linear (we assume linear kinematics) and the nonlinear components arise due to contact.In the above expressions, an explicitly differentiates this system from that used in the MPC method.
Notice that a hat is also used for submatrix components, implying that the original matrices are modified by applying coupling terms.In other words, the system arrays ( 24) are augmented with the contributions of the Lagrange multipliers, which are obtained by assembling the contributions of n λ constraints: The explicit expressions of the Lagrange multiplier contributions k c i and f c i are given next.

Contact contribution
Figure 3 shows three standard mesh nodes x i , x j , and x k .Without loss of generality and with reference to mesh node x i , we denote by x ⊥ i the position of the corresponding enriched node on the opposite surface.The position of x ⊥ i is determined by projecting x i using the unit normal vector n i .The gap function g n,i along n i and between x i and x ⊥ i is expressed as where g 0,i = X i − X ⊥ i • n i is the initial gap, X i and X ⊥ i represent the initial positions of standard and enriched nodes, respectively, and u i = x i − X i and u ⊥ i = x ⊥ i − X ⊥ i their corresponding displacements.Referring back to Figure 3, the standard FEM shape functions attached to nodes x j and x k are the only ones that contribute to the displacement field at the location of enriched node x ⊥ i .By using (18) to express its displacement, the gap reads where ⊗ denotes the Kronecker product, and I is the d × d identity matrix.Since the initial gap is constant during the analysis, the variation of ( 31) is given by The weak form of the contact contribution given by ( 15) is approximated by a summation over the n λ active contact nodes.Substituting g n,i and δg n,i with ( 31) and (32), respectively, the total contact contribution reads Noteworthy, in the discretized right hand side of (33) the Lagrange multipliers represent the force acting on the enriched nodes.By introducing δ Û i = δU ⊥ i δλ n,i , (33) can be expressed as: with Equation ( 34) can be solved iteratively, and in this work we use the generalized Newton method [96,97], where Lagrange multipliers are updated together with the primal field in the same iterative loop.In this method, which converges faster than Uzawa's algorithm, linearization of the contact contribution is required, leading to where Here k c i and f c i refer to the contact contribution of mesh node i to the stiffness matrix and the load vector in a generalized Newton loop.
It is worth noting that this formulation and its solution procedure are similar to those of NTN contact (we refer the reader to [1]).The only thing we needed to modify is the C i vector, which is based on the vector N i of enriched node x ⊥ i .

Inactive constraints
In an inactive state of contact, i.e., when λn > 0, Equation ( 15) is approximated by where n i is the number of inactive constraints, and the incremental equation is therefore expressed as For the inactive constraints, the contributions are considered in (29), which practically implies the update of the Lagrage multiplier as described in Algorithm 1. Similar to the active case previously described, the solution increments are solved using the generalized Newton method.It is worth noticing that, compared to the standard way of applying constraints, only the enriched part described in § 3.4 needs to be added.

Implementation
In this section we discuss the implementation of the method in a displacement-based FEM framework.Since the calculation of the element local arrays considering enrichment functions has been detailed elsewhere [14,19], here we only focus on the implementation of the constraints to enforce non-conforming mesh coupling and contact.Because the use of MPCs for coupling non-conforming meshes is standard [98, pp. 325-340], we only provide the detailed pseudo-code for handling contact problems (where the coupling using LMs is also included).

Coupling of non-conforming meshes
For each non-conforming node, an enriched node with the same coordinates is created on the other side of the non-conforming interface as shown in Figure 2. The element that contains this enriched node is then split into integration elements [93].An ordered tree data structure is recommended to store the associations among integration elements and their mesh parent elements.Integration elements are used to perform the numerical quadrature of elements' local stiffness and force arrays given by ( 19)- (22).The assembly of these contributions into the corresponding global counterparts follows standard procedures.MPCs or LMs are then used to enforce continuity constraints among non-conforming mesh nodes and their enriched slaves; while the application of MPCs is done according to standard procedures, the use of LMs can be regarded as a simplified version of the contact problem without a convergence test (more on this below).

Contact
We use a generalized Newton method to solve the nonlinear contact load increment following the procedure outlined in Algorithm 1 which is now described.At each load increment we create enriched nodes using the closest projection method [1,92] to determined their location (refer to Figure 3).Integration elements are created afterwards, similarly to the procedure just discussed for the coupling of non-conforming meshes.Here, as the relative displacement between contacting bodies can be neglected, we only determine the locations of enriched nodes at the beginning of each contact step, in analogy with the active set strategy in NTS contact [1].For simplicity, we first detect enriched nodes without checking their gap in the normal direction; we then compute the gap and, if it is larger than zero, these nodes make no contribution to the stiffness or internal force calculation based on the inactive constraint formulation described in § 3.5.This approach helps the convergence within a step in the case where nodes are not in contact at the start of a step, but come in contact during the iteration.
Once the number of enriched nodes is determined, the number of required Lagrange multipliers is also determined (the cardinality of set N is denoted |N | in Algorithm 1).The total number of DOFs, which includes the DOFs for displacement and Lagrange multipliers in the generalized Newton method, is used to initialize the global arrays.Notice that in the displacement and Lagrange multiplier vectors, only the elements that correspond to active constraints are kept (and zeros are added for newly added enriched DOFs and multipliers).The unconstrained global arrays are then assembled considering also the enriched contributions [14].
The following loop over the index set of enriched nodes adds the contribution of contact constraints to the global arrays.Augmented Lagrange multipliers are calculated based on the gap and penalty parameter.Note that contributions are added only if there is contact, i.e., λn,i ≤ 0.
After the solution of the system of equations and the update of the primal (displacement) and dual (Lagrange multiplier) fields, a convergence check is made.Whenever the norm of both vectors is lower than a user-specified tolerance, the analysis continues to the next contact step.As mentioned above, LMs can be used for the mesh coupling problem, which can be regarded as a contact analysis where the enriched positions (mismatching nodes) are known at the beginning of the analysis; the solution to this problem can then be obtained directly in one iteration.

Algorithm 1 Generalized Newton pseudo-code for solving a contact step
Input: Solution U and Lagrange multiplier λ vectors from previous step, sets of standard finite nodes N and elements E, ordered tree H, contact surfaces C, boundary conditions B, material properties P, and penalty parameter -get total DOFs, d ≡ DOFs per node {U , λ} ← copy(U old, λold) -copy solution vector of previous step -initialize global arrays {K, F } ← assemble(U , N ∪ Ne, E ∪ H, B, P) -assemble global arrays [14] for each i = 1, . . ., |Ne| do -loop over index set of enriched nodes -update augmented Lagrange multipler -update global residual vector -update Lagrange multipliers if ∆U < tol and ∆λ < tol then break return U , λ

Numerical examples
The accuracy and robustness of the proposed method is now demonstrated by means of numerical examples.Homogeneous linearly elastic materials and plane strain conditions are chosen.For convenience, no units are adopted so results are valid for any consistent unit system.Constant strain triangular elements are used with one point quadrature rule for both standard and integration elements.

Contact patch test
The commonly used contact patch test of Taylor and Papodopoulos [24] is used to study our method's ability to correctly transfer contact tractions.The problem consists of an elastic substrate with a rectangular punch as illustrated in Figure 4.Both the punch and the substrate are subject to a uniformly distributed vertical unit traction t along their top edges.The substrate and the punch have the same Young's modulus E = 10 and Poisson's ratio ν = 0.3.The problem is setup so that the punch comes into contact with the substrate because of the applied pressure.Since the same material is used for the punch and substrate, and the pressure is uniformly distributed on the top surfaces, the substrate should experience a constant state of strain and stress.The problem is then solved using the following discretization methods: (a) Standard FEM using a single-pass MPC; (b) Standard FEM with a two-pass MPC; (c) Our enriched method using a two-pass MPC; and (d) Our enriched method using ALM.For the single-pass MPC example, we defined the top surface of the substrate as master and the lower surface of the top block as slave (switching master and slave surfaces leads to an unconstrained system).Finally, the same strategy is used to integrate the applied pressure in all cases.
Figure 5 shows the stress field on the deformed configuration for all methods.The single-pass MPC method (panel (a)) is unable to ensure continuity and results in a non-constant stress field; notice also the interpenetration between the substrate and the punch.The two-pass method (panel (b)) ensures continuity along the contact boundary and passes the patch test.The results obtained with our enriched formulation (two-pass MPC in panel (c) and ALM in panel (d)) show that the method ensures C 0 -continuity and passes the patch test.In addition, exact integration of the applied pressure is readily possible because of the presence of integration elements with an enriched node at the location where the applied pressure is discontinuous.

Convergence study
The convergence of the proposed method is investigated by means of the classical problem of a circular hole in an infinite plate, for which the exact solution can be found in References [99,100].As shown in the schematic of Figure 6a, a square computational domain of size L = 20 with a centered whole of radius r = 4 is chosen.The material properties of the plate are E = 10 and ν = 0.3.On the boundary of the square domain we prescribe the exact displacement field corresponding to the uniform far-field traction t = ±σ ∞ e 1 , with σ ∞ = 1.
Figure 6 shows the two discretizations for this problem: a standard conforming FE mesh in panel (a) and a mesh composed of two parts that are non-conforming along the coupling interface in panel (b), where the ratio between the the element sizes in the top and bottom domains is equal to 2 and is kept constant with mesh refinement.Four different analysis approaches are compared: i) Standard FEM using conforming meshes; ii) Two-pass MPC on non-conforming discretizations; iii) Our enriched method using a two-pass MPC; and iv) Our enriched method using LM.Since the single-pass MPC cannot ensure continuity along the coupling interface, as demonstrated in the previous example, it has been discarded in this analysis.Convergence is studied via the standard, L 2 , and energy, E, norms of the error defined as where quantities with and without the superscript h refer to approximate and exact solutions, respectively.The plots in the top row of Figure 7 show the convergence results with the configuration shown in Figure 6 (horizontal interface).While the two-pass MPC method does not converge due to locking originated from an overconstrained interface [3,101], the curves of the two enriched approaches overlap and achieve the same rate of convergence as that of standard FEM with conforming discretizations: about 1 for the L 2 norm and 0.5 for the energy norm.Similar results, reported in the second row of Figure 7, are obtained when the non-conforming mesh is rotated 90 • , resulting in a vertical interface.
The energy norm of the error corresponding to the four methods using the mesh with horizontal and vertical non-conforming interfaces are shown in Figure 8, plotted per element (average value).The results corresponding to the enriched methods and standard FEM are in good agreement, whereas those obtained by the two-pass MPC (panel (b)) show some clear differences along the coupling interface.
This example shows that the performance of the proposed method for mesh coupling is basically identical to that of the standard FEM on conforming meshes.Because the same convergence rates are obtained, it can be concluded that the LBB condition is fulfilled.In contrast to the two-pass MPC method, our method avoids interface locking because enriching the primal field gives more kinematic freedom to the interface.The enrichment therefore enables an accurate representation of the mechanical behavior at the coupling interface.

Stability
Following our previous work on enriched FEM [22,93,21], we investigate the stability of the proposed method.Using the same problem geometry of the contact patch test used in § 5.1, we examine the influence of punch location and mesh size on the condition number of the system matrix.We compute the condition number of the global system matrix K as   where λ max and λ min denote, respectively, the highest and lowest (non-zero) eigenvalues of the system matrix.
No Dirichlet boundary conditions are enforced on the system and therefore we discard the lowest six eigenvalues, which correspond to the rigid body modes of both blocks.We investigate the condition number of the matrices with MPC and ALM.For each approach three variations of the enriched method are compared: ns) The enriched method without scaling enrichment functions, i.e., s i = 1 in Equation (18); os) The enriched method with the optimal scaling proposed in Reference [21], i.e., s i = 2ζ (1 − ζ), where 0 ≤ ζ ≤ 1 denotes the (relative) location of the enriched node in the finite element side that contains it; and pc) The enriched method without scaling, but using a diagonal preconditioner such that K pc = ∆K∆, where ∆ ij = δ ij / K ij is a diagonal matrix with δ ij denoting the Kronecker delta.

Effect of punch location
In the first test, the influence of the punch location on the condition number as it moves on the substrate is evaluated (see Figure 9a).Both substrate and punch are discretized with two triangular elements (see Figure 9b).Their material properties are, respectively, E 1 = 10, E 2 = 10 000 and ν 1 = ν 2 = 0.3.
Figure 10 shows the condition number of the stiffness matrix as a function of the punch location.For MPC, the condition number for the unscaled enriched method (labelled κ( Kns )) rises slightly when the punch approaches the sides of the substrate.However, when using the optimal scaling proposed in Reference [21] (labelled κ( Kos )), the condition number is the same as that of unscaled method, showing that the ineffectiveness of the scaling factor demonstrated in the one-dimensional example in Appendix A holds also in this case.Applying the diagonal preconditioner (labelled κ( Kpc )) improves the condition number.Overall, the condition number remains bounded as enriched nodes are placed arbitrarily close to standard finite element nodes when dealing with MPCs.
For ALM, the condition number of the original stiffness matrix without enrichment scaling (labelled κ( Kns )) also overlaps the one that uses optimal enrichment scaling (labelled κ( Kos )).However, with the diagonal preconditioner (labelled κ( Kpc ), the condition number improves significantly, but is still higher than that of the MPC method.

Effect of mesh size
In this second test, which is illustrated schematically in Figure 11a, we study the influence of mesh refinement with a fixed punch.The material properties are the same as those used in the first test.The results of our enriched approach are compared to those of standard FEM using conforming node-to-node contact discretizations, as shown in Figure 11b.Figure 11c shows a typical finite element discretization used for all other results.The results in Figure 12 show the condition number as a function of mesh size (left) and total number of DOFs (right).The reference curve (labeled κ(K std )) is computed using the conforming mesh shown in Figure 11b.It is well known that the condition number in standard FEM scales as O h −2 with mesh size h and O (n d ) with the total number of DOFs n d .For MPC, the condition number of the original enriched system matrix indeed deteriorates with mesh refinement (curve κ( Kns )).Also in this case, the optimal scaling proposed in Reference [21] (curve κ( Kos )) has no effect on the conditioning.However, applying the simple diagonal preconditioner improves the condition number significantly (curve κ( Kpc )).Noteworthy, the condition number of the enriched system matrices increases at the same rate as that of standard FEM with node-to-node contact.Therefore, the enriched method constrained using MPCs with a simple preconditioner is as stable as standard FEM.
For ALM, the condition numbers are generally worse than those obtained with MPCs.The condition number also deteriorates with mesh refinement, and the magnitude with and without scaling is also the same (curves κ( Kns ) and κ( Kos ), respectively).We find that the condition number of the enriched system matrix with a simple pre-conditioner (curve κ( Kpc )) is close to that of the preconditioned matrix using MPCs and grows at the same rate as FEM with node-to-node contact.

Hertzian contact problem
where t is the magnitude of applied traction, p n denotes the contact pressure, b the contact area, and E * the effective stiffness.The load is applied in twenty load increments of equal magnitude.We used the enriched method with ALM, and convergence (with a tolerance ∆U / U < 10 −5 ) was reached within four iterations per step.Two different meshes are used in this numerical test: for the first one, the substrate discretization is coarser than that of the punch, while in the second one punch and substrate are discretized with elements of similar sizes.In both cases, Figure 14 shows the presence of enriched nodes also far from the contact area.Enriched nodes are actually added also far from the contact area for convenience (in practice, all standard nodes from one contact surface are projected to the other and vice versa).Enriched nodes that do not come into contact with the corresponding standard node will be regarded as inactive in the calculation of stiffness matrix and force vector components.The results are reported in terms of stress distribution σ 22 in Figure 14 and contact pressure profile in Figure 15.The stress field shows a typical Hertzian contact distribution.Figure 15b shows the contact pressure relative error, computed as e = p n − p h n / |p n |, where p h n is the pressure obtained numerically.We find that the steep pressure gradients in elements that transition from contact to no contact are responsible of yielding inaccurate contact tractions.The error in the interior region is within 7% for both discretizations studied.Therefore, in this region the numerical contact pressure profile approximates the analytical solution closely.

Discussion
Compared to traditional contact and coupling formulations, the proposed method has a number of advantages: As the enriched formulation essentially transforms the problem into a node-to-node discretization, it is possible to utilize the most straightforward coupling and contact techniques.This also facilitates the implementation in existing standard displacement-based finite element packages.Furthermore, as the tractions are properly transferred and over-constrained locking is avoided, no contact stabilization techniques are required.This important intrinsic property of the formulation, which was first noticed while comparing DE-FEM with X/GFEM in Reference [14], relates also to the recovery of smooth tractions in immersed analysis using IGFEM and DE-FEM [16,22].Noteworthy, our traction profiles are on par with those obtained by the virtual element method (VEM) [85], whereby a non-conforming problem is transformed into a node-to-node VEM-conforming discretization; in our procedure, however, formulation and computer implementation are much simpler than VEM.It is important to note that the original element's shape functions are kept intact, and that enrichments are only nonzero in the elements along the contact boundary (thus the partition of unity property in these elements is lost).Because enrichment functions are local by construction and vanish at the original mesh nodes, the Kronecker delta property on those nodes is retained, and all standard DOFs preserve their physical interpretation.As a result, post-processing is only required to compute the solution at enriched nodes, and because every enriched node is matched to a standard node, obtaining the latter's DOFs further reduces post-processing.It has been acknowledged in previous works that interface-and discontinuity-enriched formulations may have difficulties in properly reconstructing field gradients (strains and thus stresses).This issue stems from the construction of the enriched finite element space, which may use sliver integration elements that degrade the accuracy of field gradients (as in standard FEM).While the issue is more pronounced for material interfaces [102,103], it has been shown recently that the issue is negligible near Dirichlet boundaries [22] and near traction-free cracks [93].In the context of contact and mesh coupling problems, our numerical experiments indicate that this issue is not present.As enriched nodes are only placed along the contact/coupling interfaces, we conjecture that the presence of sliver integration elements does not adversely affect the gradient field accuracy and thus the method can properly recover strains and stresses.It is worth noting that there is work that aims at improving the accuracy of recovered gradient fields from enriched FEMs (and in fact unfitted FEMs in general) [104].
Although in this work we considered linearized kinematics and frictionless contact, the extension of the current enriched framework to more advanced problems such as frictional contact, contact in 3D, and contact in large deformation is relatively straightforward.The only drawbacks we see at the moment are related to the non-symmetric global stiffness matrix stemming from the frictional contact formulation, the necessity of a more efficient way of contact detection in a three-dimensional implementation, and the possible need of smoothing techniques to achieve better convergence properties in large deformation problems.Finally, the proposed method could also be applied to high-order approximations, albeit high-order enriched functions would be needed to properly describe curved edges (assuming that geometry is described nonlinearly).
A An analytical 1D problem to investigate the role of the optimal scaling factor A simple 1D problem is proposed to investigate the ineffectiveness of the optimal scaling factor in the results reported in § 5.3.Figure 16 shows the problem, consisting of a 1D rod (element 1) that is coupled somewhere along its length at location ξ to another 1D rod (element 2) through an enriched node.
For simplicity, both elements have unit lengths, cross sections, and stiffnesses.The shape functions of both elements are The enrichment function below, endowed with scaling factor s, is defined in element 1: The element matrix for element 1, including the enrichment function contribution, is computed as while that for element 2 reads The assembled stiffness matrix of the entire system is then which clearly shows that standard and enriched elements are decoupled.Coupling is achieved by means of the multiple-point constraint relation from which, since ψ (ξ) = 1, Considering that the transformation matrix T for the multiple-point constraint is therefore from which one obtains the modified stiffness matrix Note that the scaling factor s does not appear in the modified stiffness matrix K.The scaling factor, therefore, does not have any influence on eigenvalues and condition number of the system matrix.It can be concluded that scaling of enrichment functions is ineffective towards the improvement of the condition number of enriched coupling and contact problems using MPCs.

Figure 1 :
Figure 1: Problem schematic showing the contact between two domains Ω 1 and Ω 2 along their contacting surface Γ c12 .The coupling between two domains Ω 2 and Ω 3 (perfect bonding) is also illustrated.Traction t1 and t3 are applied on boundaries Γ t 1 and Γ t 3 , and boundaries Γ u 1 and Γ u 3 are subject to a prescribed displacement.Outer surface normals are denoted n i , i = 1 . . .3. The insets show non-matching meshes at the contacting interface, with original standard FEM nodes shown in black, and enriched nodes, on the opposite surface, shown as red circles.

Figure 4 :
Figure 4: Geometry and boundary conditions of the contact patch test.The top block is placed on the middle of the substrate.

Figure 5 :
Figure 5: Deformed configurations (4× magnification) for the contact patch test showing the element stress in the e 2 direction and the contact tractions plotted at nodes with red arrows: Standard FEM with single-(a) and two-pass (b) MPC; Enriched approach with a two-pass MPC (c) and ALM (d).

Figure 6 :
Figure6: (a) Schematic of an infinite plate with a circular hole subject to a uniform traction at x 1 = ±∞.The square region indicates the computational domain used in the convergence study.The problem is analyzed using the discretizations in panels (b-c) which represent a typical conforming mesh used for standard FEM (b), and a typical mesh with a non-conforming interface (c).Convergence results are shown for the cases of horizontal and vertical interfaces (in the latter case, the mesh is rotated).

Figure 7 :
Figure 7: Convergence results for the horizontal (top row) and vertical (bottom row) non-conforming interface.The figures show the error in the L 2 norm (left column) and the energy norm (right column) as a function of the total number of degrees of freedom n d .The curves for enriched methods with MPC and LM overlap.

Figure 11 :Figure 12 :
Figure 11: (a) Fixed punch schematic; (b) Standard FEM conforming node-to-node contact mesh used used as reference; and (c) Non-conforming mesh used for the enriched methods.

Figure 14 :
Figure 14: Stress component σ 22 for different discretization designs (integration elements are drawn with red edges): (a) substrate mesh much coarser than punch mesh, (b) similar mesh sizes for punch and substrate.

Figure 13
Figure 13 illustrates a Hertzian contact problem, where a semi-circular punch with radius r = 10 and material properties E 2 = 7 × 10 5 and ν 2 = 0.3 is subject to a uniformly distributed load t = −25e 2 and pushed against a substrate with material properties E 1 = 7000 and ν 1 = 0.3.The substrate has length l = 20 and height h = 10 and is simply supported along the bottom edge.The solution to this problem is given in Reference [85]:

Figure 15 :
Figure 15: Traction along the contact interface in the Hertzian contact problem: Comparison of the numerical profiles obtained with the discretizations reported in Figure 14 to the analytical solution (a), and corresponding relative error with respect to the analytical solution (b).