Isogeometric dual mortar methods for computational contact mechanics

,


Introduction
Robust and accurate contact discretizations for nonlinear finite element analysis have been an active field of research in the past decade and a new class of formulations emerged with the introduction of isogeometric analysis (IGA) [1].IGA is intended to bridge the gap between computer aided design (CAD) and finite element analysis (FEA) by using the smooth non-uniform rational B-splines (NURBS) or T-splines common in CAD also as a basis for the numerical analysis.The use of such smooth basis functions has some advantages over classical Lagrange polynomials for FEA such as a possibly higher accuracy per degree of freedom [2,3] and, more importantly, a higher inter-element continuity.While finite elements based on Lagrange polynomials are limited to C 0 inter-element continuity independent of the polynomial order p, NURBS can be constructed with a maximum of C p−1 continuity.This high continuity results, amongst others, in a smooth surface representation which makes the application to computational contact mechanics particularly appealing, which has already been anticipated in the original proposition of IGA by Hughes et al. [1].
As a consequence, in the past five years various discretization techniques have been developed for IGA or transferred from finite element based contact mechanics to IGA, such as node-to-segment [4], Gauss-point-to-segment [5][6][7][8][9][10] and mortar methods [5,6,[11][12][13][14].We refer to the recent review in [15] for a comprehensive discussion of such methods, comparisons to their finite element counterparts and further references.In addition to the mentioned methods based on an isogeometric Galerkin approximation, the higher inter-element continuity of NURBS basis functions allows for the use of collocation methods, see [16] for a general introduction and [17,18] for an application to computational contact mechanics.Besides the discretization technique, computational contact algorithms can be distinguished with respect to the underlying solution procedure.While Gauss-point-to-segment approaches are, due to their lack of inf-sup stability (see e.g.[5,7] for numerical investigations), usually combined with a penalty approach, see [7][8][9][10], node-to-segment and mortar formulations can be combined with penalty methods [5,6], Uzawa-type algorithms [11], Lagrange multiplier methods [13,14] or augmented Lagrange methods [12].In contrast to penalty methods, the other mentioned methods fulfill the contact constraints in a discrete sense exactly.
Proposing an alternative approximation of the Lagrange multiplier, so-called dual mortar methods were originally introduced in the context of domain decomposition [19] and later extended to small deformation frictionless contact in 2D [20] and 3D [21], small deformation frictional contact [22], finite deformation frictionless contact [23,24] and finite deformation frictional contact [25].Since NURBS are naturally associated with higher order approximations, the extension of the dual mortar method to second order Lagrange finite elements in [26,27] including optimal a priori error estimates is also worth mentioning.A comprehensive review on dual mortar methods for contact mechanics can be found in [28,29].In the context of domain decomposition in IGA, optimality and stability of standard mortar methods have only very recently been investigated in [30][31][32][33], where also the construction of dual B-spline basis functions has been outlined theoretically [33].
In this contribution, we present a mortar contact formulation for IGA using a dual Lagrange multiplier method where the inequality constraints are reformulated using nonlinear complementarity (NCP) functions.In contrast to standard mortar methods, the use of a dual Lagrange multiplier basis in mortar contact algorithms yields a localization of the contact constraints due to the biorthogonality property of the dual basis functions.Thus, they can be combined with efficient nonlinear solution schemes such as semi-smooth Newton methods as an active set strategy, and the additional Lagrange multiplier degrees of freedom can easily be condensed from the global system of equations.The aim of this work is to transfer the concept of dual basis functions from finite elements to IGA, thus combining the efficiency of dual mortar methods with the beneficial IGA inherent concept of a smooth geometrical representation.The newly developed isogeometric dual mortar method is applied to both domain decomposition and (unilateral) frictional contact problems with finite deformations.In both applications, spatial convergence orders for the discretization error using uniform mesh refinement are studied numerically.To the best of our knowledge, this contribution presents the first realization of a dual NURBS approach in both domain decomposition and contact mechanics, and probably also a first proper investigation on spatial convergence orders of higher order NURBS contact algorithms in general.
The remainder of this paper is organized as follows.Section 2 reviews the three-dimensional continuum mechanical description of a two-body frictional contact problem with finite deformations.Next, NURBS basis functions are briefly introduced in Section 3 as a general tool for isogeometric analysis.These NURBS basis functions are then used in Section 4 to discretize the contact problem.Most importantly, biorthogonal basis functions with the same support as their primal counterparts are constructed and used for the discretization of the Lagrange multiplier.Numerical examples for tied contact as well as unilateral frictionless and frictional contact are presented in Section 5 and some concluding remarks are given in Section 6.

Problem definition of finite deformation frictional contact
We consider a two-body contact problem in two or three dimensions with finite deformation kinematics as illustrated in Fig. 1.
Therein, the open sets where F = ∂x /∂X and S denote the material deformation gradient and second Piola-Kirchhoff stress tensor, respectively.Dirichlet boundary conditions are prescribed by given displacement û on Γ u and Neumann boundary conditions by a given first Piola-Kirchhoff traction vector t0 using the unit outward normal N on Γ σ .In this setting, inertia forces have been neglected for the sake of simplicity.Constitutive nonlinearities are exemplarily incorporated through a hyperelastic strain-energy function ψ and the second Piola-Kirchhoff stress tensor is given as the derivative of the strain-energy function with respect to the right Cauchy-Green tensor C = F T F: Of course, other linear or nonlinear constitutive relations, such as plasticity are possible.A combination of the later presented contact algorithm with small strain plasticity [34] and finite strain anisotropic plasticity [35] based on NCP functions have already been presented.
In the following, we will refer to Ω (1) as the slave body (with the slave contact surface γ (1) c ) and to Ω (2) as the master body (with the master contact surface γ (2) c ).On the potential contact surface, we define the fundamental kinematic quantities of the contact problem with the normal gap function g n and the relative tangential velocity v τ at a point x (1) on the slave surface as with the projection P t : γ 2) mapping a point on the slave surface along its normal n onto the master surface.Such a projection along slave normals can be found e.g. in [36].Alternatively, a closest point projection is used in some cases, which, however, may not be uniquely defined everywhere; see e.g.[37] for a discussion on the solvability of closest point projections.Analogously to the kinematic quantities above, the slave side contact traction is split into a normal and tangential component: The master side contact traction follows directly from Newton's third law of motion as the negative slave side traction t (2)  = −t (1) .At the contact interface we enforce the well known Hertz-Signorini-Moreau conditions for normal contact and Coulomb's law of frictional sliding: where f is the so-called slip function for Coulomb's law bounding the Euclidean norm of the tangential contact traction by the product of the friction coefficient F and the normal contact pressure p n .Eqs. ( 1), ( 2), ( 7) and ( 8) constitute the strong form of the two body contact problem.
Next, the boundary value problem (1) and the contact constraints ( 7) and ( 8) are transformed into a variational form.Therefore, we first define the solution space U (i) and the testing space V (i) for the displacements as Moreover, a vector-valued Lagrange multiplier λ = −t c is introduced as the negative slave side contact traction to enforce the contact inequality constraints (7) and (8).Just like the contact traction in (6), the Lagrange multiplier can be decomposed into a normal and tangential part.The Lagrange multiplier is chosen from the convex set Here, W represents the trace space of U on γ c .With these function spaces, the strong form of the contact problem can be transformed into an equivalent mixed variational form by applying the method of weighted residuals and subsequent integration by parts, see e.g.[20] for a detailed derivation.As a consequence of the balance of linear momentum at the contact interface, the contact terms can be formulated as slave side integrals only and one obtains: where the first line represents the usual structural internal and external virtual work, and the second line states the contact virtual work δW c .The contact constraints are reformulated in form of a variational inequality, viz.
In Section 4, this weak form of the contact problem will be discretized using NURBS based finite elements.Special emphasis will be put on the choice of the Lagrange multiplier ansatz functions.The inequalities will be transformed into a set of non-smooth nonlinear complementarity (NCP) functions, which make them compatible with semi-smooth Newton based solvers.

Basics of isogeometric analysis
Before we go into the details of the mortar contact formulation, some basics of isogeometric analysis are recalled.This introduction is by no means all-embracing but rather meant to serve as a starting point for the following section.For a more thorough introduction to the topic, the reader is exemplarily referred to the monographs [38] for NURBS functions and [39] for isogeometric analysis.In this work, non-uniform rational B-splines (NURBS) are employed for the structural discretization in two or three dimensions, thus resulting, in a one-dimensional NURBS curve or a two-dimensional NURBS surface at the contact interface, respectively.Higher flexibility in mesh construction can be achieved using hierarchical NURBS or T-splines, see e.g.[7,8] for the application of T-splines to contact mechanics.
We start from a tensor product B-spline function space with the open non-uniform knot vector of non-descending knot values associated with the ith parametric direction of the patch.The first and last entries in the knot vector are repeated p i + 1 times where p i is the polynomial order of the B-spline basis function in the ith direction.An element is defined by two consecutive knots in each parametric direction, not counting the zero-sized elements stemming from repeated knots in the knot vector.The inter-element continuity of the B-spline function is determined by C p−m i where m i is the multiplicity of the knot common to the two elements.Throughout this work, we will consider second and third order B-splines/NURBS, the presented methods are however readily applicable to higher orders as well.A three dimensional volume V is then parameterized using control point coordinates X i, j,k as where R i, j,k (ξ, η, ζ ) represent the B-spline, or more general NURBS, basis functions of each control point.The NURBS basis functions are constructed from B-splines by introducing a weighting function: where B r are the B-spline basis functions in each parametric direction and w i, j,k is the weight associated with a control point (i, j, k).NURBS share some properties with B-splines, such as the inter-element continuity, non-negativity and partition of unity, but have the additional capability to exactly represent many conic geometries (e.g.circles).If all weights within a patch are equal, (rational) NURBS functions reduce to (polynomial) B-splines, such that the latter can be seen as a special case of the former.Since we have used open knot vectors in each parametric direction, the surface representation is solely determined by the control points and knot vectors on the respective surface of the parametric domain.Consequently, as for Lagrange finite elements, the contact formulation can later be derived with the knowledge of the boundary discretization only.Hence, the trace space restriction of the discretization of a three dimensional NURBS mesh is given by a two dimensional NURBS manifold S, viz.
For brevity of notation, we will drop the dependency on the parametric coordinates in the following.Consequently, surface and volume discretization may not be explicitly distinguishable in the nomenclature, however, the given context should provide the necessary information, e.g. when speaking of displacement shape functions at the contact interface of a three dimensional problem we are referring to (17) rather than (15).Since IGA is usually employed as an isoparametric method, said NURBS basis functions are used not only for the interpolation of discrete control point coordinates to describe the geometry, but also for the interpolation of discrete displacements d a at the control points and their variation.Assuming, for simplicity, a clamped Dirichlet boundary one gets where the dependency of the parametric direction has been dropped for improved readability, and a sum over all n cp = l • m • n control points has been introduced instead.The interpolation (18) serves as a spatial approximation of the displacements in the continuous weak form (12) with the basis functions R a .For the discretization of the Lagrange multiplier, two different approaches will be presented in the following section, namely the use of the same NURBS functions as well as a set of so-called dual NURBS basis functions.
As aforementioned, there also exists an "element-wise" point of view to B-spline or NURBS basis functions, with an element being the (non-zero) interval of two consecutive entries in the knot vector.As in classical finite elements, the parameterization of geometry and solution variables can also be represented by a set of element-local interpolations For each element n cpe =  n d i=1 ( p i + 1) NURBS functions have a non-zero value, where p i is the polynomial order in direction i and n d is the number of parametric directions, i.e. n d = 3 for volume elements and n d = 2 for (contact) surface elements.

Dual basis functions
Mortar methods are based on a variationally consistent discretization of the weak form ( 12) and ( 13), hence a suitable discrete approximation of the dual space M (11) is necessary.Given some still to be defined basis functions Φ as a basis of M h and discrete vector-valued Lagrange multipliers λ j at each control point on the potential contact surface, the Lagrange multiplier field on the slave side is approximated by where the restriction to the admissible convex set of Lagrange multipliers in (11) still has to be imposed.
Remark.In this contribution, we will always assume that every control point of the structural approximation of the slave contact surface also carries a discrete Lagrange multiplier.This is however not mandatory; for contact mechanics with second order Lagrangian finite elements, also Lagrange multiplier interpolations of lower order can be found e.g. in [26,27,40].The same is true for isogeometric thermo-mechanical contact analysis in e.g.[14].In [33], different orders of approximation of the primal (i.e.displacement) and dual (i.e.Lagrange multiplier) variables have been analyzed in the context of isogeometric domain decomposition, and it has been shown that equal order approximation is inf-sup stable and yields optimal convergence, whereas an approximation of one order lower for the Lagrange multiplier is inf-sup unstable and an approximation of two orders lower for the Lagrange multiplier is inf-sup stable but no longer optimal.
Two different types of Lagrange multiplier interpolations can be distinguished and will be termed standard and dual mortar methods in the following.Within a standard mortar approach, the same ansatz functions R a are used for approximation of the primal variable at the interface as well as the Lagrange multiplier.For contact algorithms based on Lagrangian finite elements such methods can be found in e.g.[36,[41][42][43][44] and within IGA in [12,13].Using these methods based on a standard Lagrange multiplier interpolation, a system of increased size containing both displacement and Lagrange multiplier degrees of freedom has to be solved.In this work, we follow a different approach using dual shape functions for the Lagrange multiplier, which were initially introduced in domain decomposition problems [19,45] and extended to contact problems in [20][21][22][23][24][25][26][27] and recently reviewed in [29].While dual mortar methods are meanwhile well-established in finite elements, the present work, to the authors knowledge, is the first application of dual basis functions in the context of IGA for both domain decomposition and finite deformation frictional contact.Dual basis functions are characterized by fulfilling a biorthogonality condition [19] with the Kronecker symbol δ ab .Note, that since the weak form of the contact work (i.e. the second line in ( 12)) is integrated over the current domain, also the biorthogonality is required to hold in the current configuration.Different methods to construct such dual bases exist, and we want to follow the most simple one, where the dual basis functions have the same support as their primal counterparts, fulfill a partition of unity and are constructed via element-wise linear combinations of the primal shape functions [21,[45][46][47].On each element e one readily obtains with the coefficient matrix for each element k de, (1) In the construction of the coefficient matrix, the local integration for every slave element is only performed on that part of the element domain, for which a feasible projection to the master surface is possible.This is crucial for the consistent treatment of partially projecting elements in complex contact scenarios, as has been investigated for Lagrangian finite elements in [48] for two-dimensional mortar formulations and in [49] for the general three-dimensional case.To properly detect the integration domain and reduce the integration error to a minimum, a segmentation process for isogeometric contact analysis will be described later on.For a well-defined construction of dual shape functions according to (23), the primal shape functions are required to have a non-zero integral value on the integration domain.
Higher order Lagrange polynomials do, in general, not meet this requirement which necessitates the use of a local basis transformation of the primal basis to obtain well defined dual shape functions, see [26,27].NURBS, on the other hand, are positive on the entire element, such that the construction ( 22), ( 23) is well defined without any further modifications and for any approximation order.For a two dimensional contact problem, i.e. a one dimensional contact boundary, an exemplary set of primal and dual basis functions of second order is depicted in Fig. 2. It should be pointed out that the dual basis functions generated by ( 22), (23) only guarantee a partition of unity.Consequently, the global approximation order is limited to one in the L 2 -norm, independent of the local approximation.Since the dual NURBS do not possess the optimal reproduction order, optimal convergence rates as proven in [33] cannot be guaranteed.For dual mortar methods based on Lagrange polynomials optimality can be recovered by a transformation of the primal basis [46] or by extending the support of the dual basis functions [50].An extension of the latter approach to B-splines is outlined in [33], but in general still an open question.However, for contact problems the solution is typically in H t (Ω (i) ) 3 with t < 5 /2, such that a priori estimates are already limited by the regularity of the solution.Even this simple construction of dual shape functions meets the requirements in [27] for optimal a priori estimates for the displacements in the H 1 -norm of order O(h 3 /2 ).Namely, we have  where the equality directly follows from (21) in conjunction with the element-wise construction ( 22), ( 23) together with the partition of unity property of the NURBS basis functions, and the inequality results from the non-negativity of NURBS.
• Finally, a non-negative integral of each dual basis function tested with a primal one which follows from (21) and the non-negativity of the NURBS basis functions R b .
These properties, together with a regularity assumption on the shape of the active contact zone gives the optimal a priori estimates [27] for the discrete contact problem.In terms of an application to classical domain decomposition, however, the presented dual shape functions lack the optimal approximation order, such that convergence in the H 1 -norm of order O(h p ) cannot be achieved but only O(h 3 /2 ) can be guaranteed.
Applying the spatial discretization for the displacement and its variation (18) as well as the Lagrange multiplier discretization (20) to the contact virtual work δW c in (12) gives where S and M are the sets of all slave and master side control points on the contact surface with the virtual displacements δd S and δd M , respectively, and P t,h is the discrete counterpart to the projection already introduced in (3).In the last step, we have introduced the discrete mortar coupling matrices k dγ , j = 1, . . ., n (1)  cp , k = 1, . . ., n (1)  cp , with the identity I dim ∈ R n dim ×n dim in n dim spatial dimensions.Here, the use of dual shape functions yields a diagonal mortar matrix D, which allows for an easy condensation of the Lagrange multiplier degrees of freedom from the global system of equations, see e.g.[19,[22][23][24][25].

Treatment of inequality constraints
Just like for finite elements based on Lagrange polynomials, the use of dual shape functions in the inequality constraints (13) results in a control point wise decoupled enforcement of the contact inequality constraints ( 7) and (8), which can then be reformulated as non-smooth nonlinear complementarity (NCP) functions C n, j and C τ, j for each control point j, viz.
with the two complementarity parameters c n , c t > 0 [20,22].Typical for mortar methods, the contact constraints are enforced in a weak sense by introducing the weighted normal gap gn, j and tangential slip ũτ, j : where g n,h is a discrete version of the normal gap function in (3) and u τ,h a tangential slip increment over one time step stemming from a spatial discretization and time integration of (4).Although only quasi-static problems are considered in this work, a time integration is necessary at this point, since friction (similar to plasticity) is path dependent and the (pseudo-) time acts as a path parameter, see [36] for an objective time integration algorithm for the tangential slip.A consistently linearized system of both the discrete balance equation and the complementarity function can be solved with a non-smooth variant of Newton's method [52], where the NCP functions serve as the basis for a primal dual active set strategy [53].Since the algorithm and linearization details are the same as for Lagrangian finite elements, we refer to [23][24][25] and omit the details here.
Remark.If the primal basis functions R j are used to discretize the Lagrange multiplier instead of the dual basis Φ j , the discrete inequality constraints following from (13) do not decouple for the control points.Hence, a strict variational form would require a coupled NCP function containing all n (1) cp control points.For practical reasons, however, the constraints are commonly still enforced at the control points independently [5,6,11,12], which can be interpreted as a lumping technique.This technique will also be applied here, so that for a standard mortar method we replace Φ j by R j in ( 27), ( 28) and (31).

Accurate contact integration for IGA
As aforementioned, the construction of dual shape functions in ( 22) and ( 23) needs to be restricted to the part of the element for which a feasible projection onto the master surface exists [48,49].But even if all slave elements are fully covered by the master surface, an integration based on a segmentation process as illustrated in [36,54,55] is more accurate than an element based integration usually performed on slave elements [41,56], see also [57] for a comparison in terms of accuracy and efficiency.Recently, element based integration for isogeometric mortar methods has been analyzed in [58], where also an alternative non-symmetric integration rule for the mortar integrals is used to reduce the integration error.
In an element-based integration scheme, the quadrature points to evaluate the integrals over the slave surface in ( 27), (28) and (31) are chosen as the usual Gauss integration points for each slave element.This is a rather efficient process, since no information from the master side has to be included when determining the Gauss point locations.However, the integrand in (28) contains shape functions from both the slave and master side.Considering the fact that these shape functions are defined as C 0 continuous piece-wise polynomials in FEA and C p−m continuous piece-wise rational functions in IGA, this reduced continuity is ignored in (28) if the integration points are simply chosen as the Gauss points of a slave element.This integration error becomes even worse, if the slave element is not fully covered by the master surface, because this will result in of C −1 continuity, i.e. jumps.Segment based integration schemes aim Fig. 3. Construction of integration points for one pair of isogeometric contact facets.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)at resolving this issue by splitting the integration domain into segments on which both slave and master side shape functions are C ∞ continuous.Therefore, segment based integration algorithms usually provide more accurate results.This accuracy, however, comes at the expense of a rather intricate determination of the Gauss point locations (and their linearizations), since information about the master side mesh needs to be included.For isogeometric contact analysis, the element based integration is far more popular, only [30] and [14] have already used some sort of segmentation procedure in domain decomposition and contact mechanics, respectively.
Here, we propose a very simple but efficient extension of the segmentation procedure described in [24,36,54,55] for standard first order finite elements to IGA similar to [14], that merely consists of pre-and post-processing our already existing segmentation algorithm as depicted in Fig. 3. Starting from a pair of slave and master elements (solid black lines) given by the control points (dots) and control polygons (dashed lines), we construct an auxiliary integration element s int and m int (solid red lines), respectively, as a four node quadrilateral connecting the corners of the NURBS surface element in the spatial domain (see top left of Fig. 3).These integration elements can then be treated using standard mortar segmentation algorithms, for example the ones from [24,36,54,55], that construct an auxiliary plane at the slave integration element center (Fig. 3, box top left), project both slave and master element into this plane (Fig. 3, box top right), perform a polygon clipping of the projected elements (Fig. 3, box bottom left) and a triangulation of the resulting clip polygon (Fig. 3, box bottom right).On those triangular integration cells, quadrature points are defined.In a final post-processing step, these quadrature points are mapped back onto the original slave element (Fig. 3, bottom right).We emphasize, that the simplified quadrilateral integration elements are only used to determine the Gauss point locations ξ gp , whereas the actual integration as well as the entire actual contact algorithm is still based on the fully isogeometric surface description, such that the smoothness of the isogeometric surface description is not lost.The presented segmentation algorithm can be viewed as a beneficial combination of the two existing isogeometric segmentation procedures in [14,30].Where [30] searches for intersections between NURBS curves, our method only needs intersections of straight lines, since the intersection phase is already based on bilinear integration elements.On the other hand, [14] also constructs a piece-wise linear approximation of the contact surface.However, the introduced linear approximation therein is used not only for determining the Gauss point location but also the interpolation of the Lagrange multiplier.

Numerical results
In the following, we investigate the proposed isogeometric dual mortar methods in three numerical examples and compare the results to standard isogeometric mortar methods as well as contact formulations based on Lagrange polynomials.The first order Lagrange polynomial (and the equivalent first order NURBS) approximation will be denoted as L 1 , second order Lagrange polynomials as L 2 , second order NURBS as N 2 and third order NURBS as N 3 .All of which can be combined with either a standard or dual Lagrange multiplier interpolation.
In a first example, the dual mortar method for domain decomposition is analyzed for linear elasticity.The problem of an infinite plate with a circular hole under tension is simulated to investigate the convergence order in a domain decomposition case and shows predicted reduced convergence order O(h 3 /2 ).The second example is a Hertzian type contact problem and demonstrates the contact stress distribution in a two dimensional frictionless contact setting.Due to the switch-over from domain decomposition to unilateral contact, the reduced regularity of the solution already limits the spatial convergence order.Therefore, both standard and dual mortar methods now converge equally optimal at order O(h 3 /2 ).Finally, a three dimensional rotating ironing example is simulated to compare the isogeometric formulation to approximations based on Lagrange polynomials of first and second order.This last example incorporates both finite deformations and large frictional sliding.All simulations have been performed with our parallel in-house research code BACI [59].

2-dimensional tied contact: convergence study
In a first validation step, we analyze a problem of tied contact, i.e. without the inequality nature of the unilateral contact constraints.Such a tied contact problem results in a classical domain decomposition application of the mortar method which then simply couples two independent discretizations with non-matching meshes at the interface.To be able to compare numerical results with an analytical solution, we simplify the non-linear elasticity problem stated above to linear kinematics, introducing the linearized strain tensor ϵ = 1 /2(∇u + (∇u) T ).Linear elasticity is assumed relating the Cauchy stress σ with the linearized strain tensor via Hooke's law σ = C : ϵ using the fourth order elasticity tensor with the components with Young's modulus E and Poisson's ratio ν.
Similar to [39], we analyze an infinite plate with a circular hole under constant in-plane tension in plane strain condition, for which an analytical solution can be found e.g. in [39].We adopt the geometrical setting given in [39] using a quarter of the plate with appropriate symmetry conditions and exact surface traction boundary conditions, see Fig. 4(a) for the dimensions and material properties.In [39], a repeated control point is used to describe the geometry with only one patch of two elements, which, however, results in a singular point in the inverse mapping at the repeated node.Instead, we use two patches such that the mapping of the isogeometric surface and its inverse is well defined everywhere within the domain.Two different interfaces are considered, a straight one where the NURBS basis functions at the interface reduce to B-splines (i.e.piece-wise polynomials) and a curved interface with general NURBS basis functions.We perform uniform knot insertion for both patches to obtain an element size ratio of 2:3 at the interface, see Fig. 4(b) and (c) for the coarsest levels, and keep this ratio fixed in the following.Uniform refinement via knot insertion is analyzed for second and third order NURBS basis functions.Note, that the two patches do not share any control points at the interface, but the coupling is enforced using the mortar method.
Convergence of the approximated solution u h towards the analytical solution u is assessed in terms of the energy norm For an interpolation of the Lagrange multiplier using equal order NURBS basis functions one can expect convergence of order O(h p ) in the energy norm where p is the polynomial order of the B-spline basis, whereas for the presented dual Lagrange multiplier basis only O(h 3 /2 ) can be guaranteed theoretically for any p ≥ 2. This is due to the reduced approximation order of one of the dual Lagrange multiplier basis.Figs. 5 and 6 compare the convergence behavior for standard and dual Lagrange multiplier shape functions for second and third order NURBS and different choices of the slave side for the straight and curved interface, respectively.As expected, a standard Lagrange multiplier method yields optimal convergence of order O(h p ) for both choices of the slave side.In the isogeometric dual mortar approach, on the other hand, the choice of the slave side has an influence on the convergence orders.If the coarser side is chosen as slave, the convergence order drops to the theoretically predicted O(h 3 /2 ).If, however, the finer side is chosen as slave, order O(h 2 ) convergence can be observed.The gain of order O(h 1 /2 ) may be the result of super convergence effects.However, theoretical investigation of which is beyond the scope of this paper.Next, we compare the results above with approximations based on Lagrange polynomials of first and second order.For the sake of brevity, we restrict ourselves to the more general case of a curved interface, since the results are qualitatively the same.Fig. 7 displays the convergence behavior of first and second order finite elements with different approximations of the Lagrange multiplier.Here, all methods converge with the optimal order of O(h 1 ) and O(h 2 ), respectively, even the theoretically most critical one using second order dual basis functions in the s:m = 2:3 setting.This is due to the fact that the used discretization with second order 9-node quadrilaterals always places the side-mid nodes and the central nodes equally spaced, i.e. at the locations of the Gauss-Lobatto points.Therefore, O(h 2 ) convergence can be expected, see [46].For orders p > 2, equally spaced nodes do no longer correspond to Gauss-Lobatto points and convergence orders for the dual case may deteriorate.Such higher order Lagrange polynomial approximations are, however, beyond the scope of this paper and the reader is referred to [46] for a detailed investigation of this topic.
Finally, we compare the convergence rates of different dual mortar interpolations in the most critical case, i.e. the curved interface with s:m = 2:3 and a dual Lagrange multiplier approximation.Fig. 8 compares the isogeometric results from Figs. 5(b) and 6(b) with the results using classical finite elements in Fig. 7.
Basically, two ways are possible to compare finite elements and IGA, either by the element size (Fig. 8(a)) or by the number of nodes, respectively control points (Fig. 8(b)).Since IGA requires less control points per element, these two approaches yield different results.While the former compares the efficiency in terms of element evaluations, the latter rather focuses on the size of the resulting linear system to be solved.The complexity of solving this linear system is, however, not only determined by its size, but also the bandwidth, which is increasing with the polynomial order.Obviously, the first order approximation yields the largest errors and converges slowly.Comparing the second order methods L 2 and N 2 , one observes that when taking the element size as a reference the finite element version converges faster and yields lower errors from the beginning.However, taking the number of control variables as a reference, the isogeometric discretization N 2 is more accurate than the finite element counterpart L 2 , despite the reduced convergence order.Only if the already very fine mesh is further refined, the L 2 version would gain the advantage due to faster convergence.The third order NURBS version N 3 gives the best results in coarse meshes and behaves similarly to the N 2 case in the limit.
In summary, the presented simple element-wise construction of dual NURBS basis functions may yield a deterioration in the convergence orders for domain decomposition.Thanks to the higher accuracy per degree of freedom, however, the isogeometric dual mortar method is, despite sub-optimal convergence, competitive to classical finite elements using Lagrange polynomials.To obtain optimal convergence, more sophisticated methods of constructing the dual basis are necessary, e.g.dual shape functions with an extended support [33,50].Nonetheless, the presented method provides an easy and efficient way of coupling patches without increasing the global system size (compared to [33]) and gaining a localization of the coupling (compared to [32]).

2-dimensional Hertzian-type contact: convergence study
Although the presented element-wise construction of dual shape functions yields sub-optimal convergence in domain decomposition applications, they may still be interesting for unilateral contact applications.In this case, the spatial convergence is usually limited by the reduced regularity of the solution, such that even the simple element-wise construction gives optimal convergence in finite element analysis [27].Hence, in our second example, we want to investigate the spatial convergence properties of the isogeometric dual mortar contact algorithm in detail.We therefore use a two dimensional Hertzian-type contact of a cylindrical body (radius R) with a rigid planar surface under plane strain conditions.The two horizontal upper boundaries undergo a prescribed vertical displacement.To avoid singularities in the isogeometric mapping, we introduce a small inner radius (radius r ), see Fig. 9 for the geometric setting, the material parameters and the coarsest mesh.Again meshes using second and third order NURBS basis functions are used as depicted in Fig. 9 for the coarsest level, where different Bézier elements are marked with different shading.In this setup half of the elements on the potential contact surface are located within one ninth of the circumferential length and C p−1 continuity is ensured over the entire active contact surface.In the convergence study, uniform mesh refinement via knot insertion is performed on each of the patches resulting in a constant local element aspect ratio.Although only relatively small deformations are to be expected, we use a fully nonlinear description of the continuum using nonlinear kinematics and a Saint-Venant-Kirchhoff material under plane strain assumption as well as the nonlinear contact formulation from Section 4.
In Table 1, we compare different refinement levels and study the convergence behavior in terms of the energy norm.Since no analytical solution is available, we use the finest mesh of level 7 with standard third order NURBS as a numerical reference solution.Table 1 gives the error decay over six refinement levels for both a standard and dual Lagrange multiplier interpolation of second and third order together with the numerical convergence order in each step.In the limit, all methods converge with the expected order of O(h 3 /2 ) in the energy norm and also the absolute error values are quantitatively very similar.Only the N 3 standard case gives a slightly higher order in the last step since the next level of this mesh is chosen as the numerical reference solution.In view of Table 1, the use of dual shape functions for the Lagrange multiplier instead of primal ones does not come at the expense of a reduced accuracy but yields equally accurate results while reducing the total system size to the number of displacement degrees of freedom only.In contrast to the domain decomposition case above, the convergence is now limited by the regularity of the solution, such that both standard and dual interpolations converge with the same order.The use of higher order NURBS, i.e. third order in Table 1 or even higher seems questionable from this viewpoint, since no faster convergence is gained from the higher order interpolation with uniform mesh refinement.
In the context of isogeometric contact formulations, Hertzian-type contact settings are commonly used to assess the smoothness of the pressure distribution using different computational methods, see e.g.[5][6][7]12,13,15].Strictly  speaking, the Lagrange multiplier in a variationally consistent contact formulation is not a pointwise defined function but an element of the dual space of the trace space.Thus, we cannot evaluate the Lagrange multiplier pointwise with respect to the space coordinates, but only as a linear functional.This observation motivates that in the discrete setting we associate with λ h in (20) a NURBS function λh specified in terms of the non-negative basis functions R a and the already computed discrete control point values.More precisely the definition guarantees that λh is strongly non-negative.Moreover, due to (24) we find that Thus, both λh and λ h impose the same mean contact pressure per face e on the contact zone.Following [60], it can be shown that the modified discrete Lagrange multiplier has the same convergence order and, moreover, in contrast to the standard low order finite element approach, a C p−1 continuous contact pressure is found.Figs.10-12 compare the pressure distributions for the meshes 2, 4 and 6 using second order NURBS with the standard Lagrange multiplier approach and the dual approach using either dual or standard shape functions for the pressure post-processing.Therein, red circles indicate the discrete control point values of the Lagrange multiplier.Since the given example does not exactly model a Hertzian contact problem (nonlinear kinematics vs. linear Hertzian theory, thick hollow cylinder vs. full circular disk, displacement control vs. point load), the contact pressure of the level 7 with the post-processed dual N 2 approximation is depicted as a reference with dashed blue lines instead of the solution according to Hertzian theory.As expected, the visualization using dual shape functions results in a strongly oscillatory behavior due to the discontinuous and not even strictly non-negative nature of the dual shape functions.However, the post-processed contact traction for the dual Lagrange multiplier method yield even smoother and less oscillatory results than the standard case on the same mesh.Especially with the medium mesh (Fig. 11) the use of standard NURBS basis functions for the Lagrange multiplier yields oscillations in the contact stress; a similar behavior can also be observed in various IGA contact formulations, e.g. in [5,6,12,15].Those oscillations vanish completely for the post-processed dual mortar method.This effect also transfers to higher order NURBS basis functions: Fig. 13 depicts the contact stress distributions for third order NURBS and mesh 4. Again, the standard isogeometric mortar method shows some oscillations in the contact pressure whereas the dual mortar method yields a perfectly smooth result.Since the oscillations in the standard mortar method usually occur at the boundary of the active contact zone, we suspect them to be a result of the localized active set strategy.As mentioned in Section 4.2, only the use of dual basis functions yields decoupled constraints in a consistent way.For standard mortar methods, the imbalance of coupled virtual work contributions for the control points in (26) and decoupled constraints may be responsible for the oscillations.This being only a suspicion at the current stage, the oscillations in standard mortar methods may require further study in the future.Next, we compare four different dual mortar methods: we therefore use the NURBS approximation of mesh 4 from above and generate an approximation using the same number of elements, but based on first order Lagrange polynomials (denoted as L 1 ).Additionally, a second order Lagrange polynomial approximation is generated by keeping the number of nodes fixed and elevate the order to second order Lagrange polynomials (denoted as L 2 ).(c) and 13(c) for illustrative purposes.One observes the usual behavior for higher order contact formulations: starting from the common ancestor L 1 /N 1 and a piece-wise linear approximation of the contact stress, order elevation for Lagrange polynomials yields oscillations and locally negative contact traction.Order elevation using NURBS basis functions, on the other hand, gives smoother contact pressure distributions.This trend becomes even more pronounced, if the order is further increased, see [15].Comparing the two isogeometric methods N 2 and N 3 in Figs.14(c) and 14(d), it has to be noted that in the higher order case N 3 the contact pressure distribution gets "smeared" over a larger region.This is a consequence of the inability of smooth NURBS basis functions to accurately represent the non-smooth contact pressure solution.Increasing the approximation order in the IGA case while keeping the inter-element continuity at a maximum aggravates this deficiency.Finally, the adequate representation of frictional contact tractions should be investigated using this setup.Therefore, the problem is enhanced with Coulomb friction with a friction coefficient F = 0.75.Fig. 15displays the solution of mesh 4 and the results are again compared to the solution of mesh 7 with the post-processed dual N 2 approximation.Obviously, none of the relatively coarse discretizations are able to fully capture the sharp kink in the frictional contact traction at the stick/slip transition.However, the observations already made for the normal contact pressure appear even more pronounced here.For the standard Lagrange multiplier approximation the oscillations in the normal contact pressure also influence the frictional response, thus yielding oscillatory tangential tractions in the central region.The dual approximation seems to be oscillatory at first sight, but this is not the case since the discrete control point values are smoothly distributed.Quite in contrast, the oscillations in the standard approximation are actually completely removed in the post-processed dual solution where, a smooth distribution of both normal and tangential contact traction is obtained.

3-dimensional rotating ironing: global contact forces
A final example demonstrates the applicability of the presented isogeometric dual mortar method to large frictional sliding and again compares the result to contact formulations based on Lagrange polynomials.Inspired by [15], we therefore choose a rotating ironing example and enhance it with Coulomb friction (friction coefficient F = 0.1).An   any qualitative difference.In summary, the known advantages of isogeometric contact algorithms also transfer to the dual mortar approach: keeping the inter-element continuity at a maximum, higher order NURBS give smoother contact forces, whereas higher order Lagrange polynomial finite elements do not.

Conclusions
In this paper a dual mortar formulation for isogeometric analysis has been derived and investigated numerically.In contrast to standard mortar methods, the use of dual basis functions allows for an easy condensation the Lagrange multiplier from the global system of equations.Here, we used a very simple construction of the dual basis at element level.Those basis functions have the same local support as the primal functions and fulfill a biorthogonality condition, but lack the optimal reproduction order.
The application of the isogeometric dual mortar method is twofold.First, in the classical sense of mortar methods, it can be used for domain decomposition problems such as tied contact in solid mechanics.In IGA, this is an interesting topic since it is often necessary to construct a geometry by multiple patches that need to be coupled.Here, the dual mortar method presents a promising approach in principle, since it provides a sound mathematical basis, does not enlarge the global system of equations and allows for a localized transfer between the slave and master side.However, due to the currently chosen element wise construction of the dual basis functions, the convergence order may reduce to order O(h 3 /2 ) in this application.This may be cured in two ways, either by an appropriate choice of the primal basis [46] or by an enlarged support of the dual basis functions [50].Despite being quite well analyzed in the context of finite elements, the construction of optimal dual basis in IGA is far from trivial.However, since efficient and accurate patch coupling methods are of great interest in IGA, this will be a topic of future research.
The second field of application for dual mortar methods is contact mechanics.Here, large deformation frictional contact has been considered.In contrast to domain decomposition methods, due to the less regular solutions only lower order of convergence can be expected anyhow, even for standard mortar methods.Hence, the dual interpolation exhibits the same order of convergence as the standard one, while enabling an easy condensation of the additional Lagrange multiplier degrees of freedom and resulting in a positive definite system.Compared to finite element analysis, the isogeometric approach yields smoother contact pressure distributions, especially for higher order ansatz functions.The isogeometric dual mortar method is combined with a simple post-processing of the contact pressure and gives even less oscillatory contact pressure distributions as the standard isogeometric approximations.Therefore, the presented isogeometric dual mortar method constitutes an accurate and efficient and thus very promising algorithm for contact treatment in IGA.

Fig. 1 .
Fig. 1.Notation of a two body contact problem.

c
define the reference configuration of the two contacting bodies and Ω (i) t the respective spatial configuration.The boundaries of the two bodies are divided into three disjoint sets Γ denoting the Dirichlet boundary of prescribed displacements, Neumann boundary of prescribed traction and the potential contact boundary, respectively.The spatial configuration of these surfaces is given by γ each body i the boundary value problem to be solved is defined as Div

( 1 )
c , M its corresponding dual space, and ⟨•, •⟩ γ (1) c denotes the scalar or vector-valued duality pairing between W and M on γ cp a=0 Φ a = 1 by construction; proof see[21] • inf-sup stability, see e.g.[51, Remark 2.11], • a positive mean value on each element within the support of a dual basis function  e Φ a de =  e R a de > 0 (24) (a) Standard one dimensional B-spline basis functions of second order.(b) Corresponding dual basis functions of second order.

Fig. 2 .
Fig.2.Primal and dual basis functions for a one dimensional B-spline example.The equations for the shape functions in the central element are given to underline the desired biorthogonality and partition of unity.

Fig. 4 .Fig. 5 .
Fig.4.Infinite plate-Problem setup and coarsest mesh represented by the control points (dots and boxes) and mesh (dashed and dotted lines) for the two patches of second order NURBS.

Fig. 6 .
Fig.6.Infinite plate-Convergence of different slave-master pairings and Lagrange multiplier shape functions using third order NURBS basis functions.
(a) First order finite elements.(b)Second order finite elements.

Fig. 7 .Fig. 8 .
Fig. 7. Infinite plate-Convergence of different slave-master pairings and Lagrange multiplier shape functions for finite elements based on Lagrange polynomials.

Fig. 14 compares
Fig.14compares the obtained contact stress distribution from the four dual mortar methods.Figs.14(c) and 14(d) are reproduced from Figs.11(c) and 13(c) for illustrative purposes.One observes the usual behavior for higher order contact formulations: starting from the common ancestor L 1 /N 1 and a piece-wise linear approximation of the contact stress, order elevation for Lagrange polynomials yields oscillations and locally negative contact traction.Order elevation using NURBS basis functions, on the other hand, gives smoother contact pressure distributions.This trend becomes even more pronounced, if the order is further increased, see[15].Comparing the two isogeometric methods N 2 and N 3 in Figs.14(c) and14(d), it has to be noted that in the higher order case N 3 the contact pressure distribution gets "smeared" over a larger region.This is a consequence of the inability of smooth NURBS basis functions to accurately represent the non-smooth contact pressure solution.Increasing the approximation order in the IGA case while keeping the inter-element continuity at a maximum aggravates this deficiency.Finally, the adequate representation of frictional contact tractions should be investigated using this setup.Therefore, the problem is enhanced with Coulomb friction with a friction coefficient F = 0.75.Fig.15displaysthe solution of mesh 4 and the results are again compared to the solution of mesh 7 with the post-processed dual N 2 approximation.Obviously, none of the relatively coarse discretizations are able to fully capture the sharp kink in the frictional contact traction at the stick/slip transition.However, the observations already made for the normal contact pressure appear even more pronounced here.For the standard Lagrange multiplier approximation the oscillations in the normal contact pressure also influence the frictional response, thus yielding oscillatory tangential tractions in the central region.The dual approximation seems to be oscillatory at first sight, but this is not the case since the discrete control point values are smoothly distributed.Quite in contrast, the oscillations in the standard approximation are actually completely removed in the post-processed dual solution where, a smooth distribution of both normal and tangential contact traction is obtained.

Fig. 17 .
Fig. 17.Rotating ironing-comparison of vertical contact forces for different approximation schemes.

Fig. 18 .
Fig. 18.Rotating ironing-comparison of horizontal contact forces for different approximation schemes.