Varyingorder NURBS discretization: An accurate and efficient method for isogeometric analysis of large deformation contact problems
Abstract
In this paper, a novel varying order NURBS discretization method is proposed to enhance the performance of isogeometric analysis within the framework of computational contact mechanics. The method makes use of higherorder NURBS for contact integral evaluations. Lowerorders NURBS capable of modelling complex geometries exactly are utilized for the bulk discretization. This unexplored idea provides the possibility to refine the geometry through controllable order elevation strategy for isogeometric analysis. To achieve this, a higherorder NURBS layer is used as the contact boundary layer of the bodies. The NURBS layer is constructed using the surface refinement strategies such that it is accompanied by a large number of additional degrees of freedom and matches with the bulk parametrization.
The validity of the presented isogeometric mortar contact formulation with varyingorder NURBS discretization is first examined through the contact patch test. The capabilities and benefits of the proposed methodology are then demonstrated in detail using twodimensional frictionless and frictional contact problems, considering both small and large deformations. It is shown that using the proposed method, accurate solutions can be achieved even with a coarse mesh. It is also shown that the current method requires a considerably lower computational cost compared to standard NURBS discretization while retaining robustness and stability. The simplicity of the method lends itself to be conveniently embedded in an existing isogeometric contact code after only a few minor modifications.
Keywords: Computational contact mechanics; Isogeometric analysis; NURBS; Mortar method; Patch test; Frictional contact.
isbn
\clearfieldissn
\clearfieldmonth
\renewbibmacrovolume+number+eid\printfieldvolume\setunit \printfieldnumber\setunit \printfieldeid
\DeclareFieldFormat[article]number\mkbibparens#1
\DeclareFieldFormat[article]volume\mkbibbold#1
^{†}^{†}affiliationtext: Department of Mechanical Engineering,
Indian Institute of Technology Guwahati, Guwahati, Assam, India, 781039
1 Introduction
In the last two decades, NURBSbased isogeometric analysis (IGA) [hughes2005] has been established as an advantageous computational technology for various classes of problems, especially to those where the geometric approximation adversely influence the accuracy of the solution. This is attributed to the distinguished intrinsic features of its underlying basis function, viz. the ability to represent complex geometry exactly even with a coarse mesh, variational diminishing and convexhull properties, tailorable interelement continuity, and nonnegativeness [hughes2005, Cottrell2009]. Contact modelling belongs to one of these classes that has particularly been benefited from the aforementioned features of the IGA technology. Higherorder smoothness of the NURBS discretized geometry directly provides the continuous normal vector field across the boundary of the contact elements. Thus, eliminates the need of additional contact surface smoothening approaches [Padmanabhan2001, ElAbbasi2001, Wriggers2001, Krstulovi2002, AlDojayli2002, Stadler2003] often utilized in the context of finite element method.
The first investigations on the treatment of contact problems using isogeometric analysis are conducted by Temizer et al. [Temizer2011, Temizer2012], Lu [Lu2011], and De Lorenzis et al. [DeLorenzis2011, DeLorenzis2012].
Temizer et al. [Temizer2011] applied the isogeometric analysis to threedimensional thermomechnical frictionless contact using the mortar contact formulation. They demonstrated that NURBS based discretization achieves superior results in terms of quality and robustness over its counterpart Lagrange discretization. Their investigation of the classical Hertz contact problem showed that the stress oscillation error near the boundary of the contact zone reduces with increasing the order of NURBS interpolations. However, a very fine mesh is still required to get the result which closely matches with the exact solution. Moreover, the nonmortar contact formulation is found to be overconstrained in nature and leads to inaccurate results. Lu [Lu2011] combined the isogeometric analysis with the segmenttosegment contact formulation presented by Papadopoulos and Taylor [Papadopoulos1992]. He showed that intrinsically smooth NURBS discretized geometry alleviates the nonphysical oscillations of the contact forces and is capable of accurately describing the intricate mechanics of smooth materials such as fabrics. De Lorenzis et al. [DeLorenzis2011] introduced a mortarbased contact formulation in the context of IGA for the investigation of twodimensional large deformation contact with Coulomb friction. They demonstrated that the magnitude of the nonphysical oscillations of the normal and tangential contact reaction forces reduce on increasing the order of NURBS interpolations. They also showed that the oscillation error in the Hertzian contact stress at the boundary of the contact region is reduced if a very fine mesh is used. Moreover, the distributions of the contact responses obtained with the Lagrange based discretizations are found to be highly sensitive to the interpolation order. Later, De Lorenzis et al. [DeLorenzis2012] applied IGA to threedimensional large deformation frictionless contact with the mortar method. The contact constraints are enforced using the augmentedLagrangian approach [ALART1991]. Temizer et al. [Temizer2012] introduced a threedimensional mortarbased frictional contact formulation as an extension of works in [Temizer2011, DeLorenzis2011]. Again, it is shown that increasing the order of the NURBS interpolations improves the smoothness of contact reaction forces. A pointtosegment contact formulation for twodimensional frictionless contact is presented by Matzen et al. [Matzen2013], and recently extended to a weighted pointbased contact approach for threedimensional contact [Matzen2016], both with the IGA. Many researchers have extensively studied the application of mortar method in the context of IGA [Kim2012, Temzier2013Multiscale, Dittmann2014, BRIVADIS2015, SEITZ2016, Duong2018], e.g. Brivadis et al. [BRIVADIS2015] investigated the isogeometric mortar method from the theoretical as well as the numerical point of view. Seitz et al. [SEITZ2016] developed a dual mortar method for isogeometric contact analysis. Recently, Duong et al. [Duong2018] introduced a segmentationfree isogeometric mortar contact formulation. A number of works have explored the application of collocation methods to different contact problem using IGA, e.g. [KRUSE2015, DELORENZIS2015, Weeger2018]. De Lorenzis et al. [DeLorenzis2014] presented a detailed review of various isogeometric based treatment procedures for contact problems.
From the abovereviewed literature on contact modelling using IGA, it is evident that NURBSbased discretization delivers significantly superior performance in terms of accuracy, stability, and robustness compared to Lagrange discretization. However, the application of existing NURBSbased discretization approach to isogeometric contact analysis is computationally expensive since due to the rigid tensor nature of the NURBS structures the mesh can only be refined in a global manner. Moreover, the interpolation order of NURBS functions employed for the discretization of the contact boundary layer and for the remaining bulk domain of the geometry has to be uniformly elevated. From the analysis point of view, this may not be desirable since the accuracy of the contact solution is primarily governed by the description of the contact interface. According to literature, to enable the local mesh refinement in the context of isogeometric contact, Tsplines [BAZILEVS2010Tspline, DORFEL2010, Scott2011, Dimitri2014Tspline, Dimitri2014, Dimitri2017], NURBSbased hierarchical refinement [Temizer20161, Temizer20162], and locally refined (LR) NURBS [Sauer2017] approaches have been adopted. Tsplines have been utilized for studying cohesive/contact interface debonding [Dimitri2014Tspline, Dimitri2014, Dimitri2017]. Temizer and Hesch [Temizer20162], and Hesch et al. [Temizer20161] used the hierarchical NURBS for the treatment of frictionless and frictional contact, respectively. Zimmermann and Sauer [Sauer2017] introduced LR NURBS for the analysis of contact computations of solids and membranes. However, the idea to refine the geometry through controllable order elevation strategy remains unexplored as already noted by Temizer et al. [Temizer2012].
In the context of finite element contact analysis, to avoid the usage of higherorder NURBS for the bulk description, considerable research progress has been made that combines the intrinsic features of NURBS with the FE discretization [Corbett2014, Corbett2015, RASOOL2016182, MALEKIJEBELI2018, Otto2018, DIAS2019]. Corbett and Sauer [Corbett2014, Corbett2015] introduced a NURBSenriched contact element formulation which combines the geometric smoothness of NURBS with the efficiency characteristic of the FE discretization. The potential contact layer of each finite element is locally replaced by a NURBS layer, resulting in “NURBSenriched contact finite elements.” Their work is based on the contact element enrichment strategy of Sauer [Roger2011_enrichment, Sauer2013], which used higherorder Lagrange and Hermite functions on the contact surface. Later, MalekiJebeli et al. [MALEKIJEBELI2018] proposed a hybrid isogeometricfinite element discretization method where the advantages of NURBS are exploited in twodimensional cohesive interface contact/debonding. The bulk domain is given by FE discretization. The transition from NURBS to FE discretizations is carried out using the so called “transition elements”, originally presented in [Corbett2014]. Otto et al. [Otto2018] proposed a coupled FENURBS discretization approach where an auxiliary NURBS layer is placed between the contact zone of higherorder FE discretized contacting bodies. In order to tie the NURBS layer with FE discretization the pointwise and mortar mesh tying approaches are used. Recently, Dias et al. [DIAS2019] presented a higherorder mortarbased contact element, where the contact interface and bulk are discretized using the FEM. In contrast to enrichment approach of Corbett and Sauer [Corbett2014, Corbett2015], they used higherorder Lagrange functions as a basis for the discretization of the contact interface, domain, and approximation of the solution field. NURBS are only used to map the position of FE contact layer nodes.
On the other hand, in the context of isogeometric contact analysis, no such effort has been devoted that directly allows the controllable order elevation [Temizer2012] and accordingly complements the performance of IGA while fully retaining its intrinsic key features.
In the present work, a novel varyingorder NURBS discretization method is thus proposed where the userdefined higherorder NURBS interpolations are used for the contact computations. The minimumorder NURBS capable of representing the complex geometries exactly are used for the description of the bulk domain. To achieve this, a higherorder NURBS layer is used as the contact boundary layer of the body. This, as a result, avoids the application of higherorder NURBS in the region away from the contact interface. The layer is constructed using different surface refinement strategies in such a manner that it is accompanied by a large number of additional degrees of freedom and matches the bulk parametrization. The proposed varying order NURBS discretization method is accordingly denoted by NN, where N is the order of the NURBS basis utilized for the description of the bulk domain and N is the order of NURBS for the contact boundary layer. In the current work, this methodology is applied to twodimensional frictionless and frictional contact problems, considering both small and large deformations. Due to its simplicity, the proposed approach can be conveniently embedded into an existing IGA contact code after only a few modifications.
The mortar method, originally introduced as a domain decomposition technique [mortar1990], ensures the stability of the contact solution and provides optimal convergence rates, see e.g. Fisher and Wriggers [Fischer2005], Puso and Laursen [PUSO2004], Tur et al. [TUR2009], Hesch and Betsch [Hesch2009], and Popp et al. [POPP2013]. In this work, the mortar isogeometric contact formulation presented by De Lorenzis et al. [DeLorenzis2011] is extended for varyingorder NURBS discretization. The impenetrability and sticking constraints are enforced in a weak sense at the active control points. For the regularization of the impenetrability and tangential sticking constraints the penalty method is adopted.
The remainder of the paper is structured as follows. In Section 2, the mortar based contact formulation for the large deformation frictional problem in the continuum setting along with the standard NURBS discretization procedure are presented. The varyingorder NURBS discretization method and its implementation into the existing IGA contact code are described in Section 3. In Section 4, three, twodimensional numerical examples are presented. First example examines the validity of the mortar isogeometric contact formulation with varyingorder NURBS discretization using the contact patch test. In the next two examples, the performance of the proposed discretization method is demonstrated using the classical Hertz and large deformation ironing contact problems. Finally, Section 5 concludes the paper with possible future directions.
2 Contact problem description and NURBS discretization
In this section, first, we present the computational formulation for twodimensional large deformation frictional contact between two deformable bodies [laursen2003, Wriggers2006]. Next, the existing NURBSbased discretization approach used for the continuum and the contact boundary layer is described. The associated issues are also highlighted.
2.1 Computational contact formulation
It is assumed that two hyperelastic bodies, where one of them is denoted as slave, , and other as master, , come into contact and undergo finite deformations.
The current configuration of body , is given by , where and represent the coordinates of a generic point in the reference and current configuration, respectively, and is the displacement field.
The master surface is parametrized using the convective coordinate that defines covariant vector , where the contravariant vector and .
In the current configuration, the contact surface of body is described as (assuming perfect contact). The closest contact point corresponding to fixed slave point is computed via intersecting the master surface with a line that passes through in the direction of unit normal vector with closest point projection procedure [KONYUKHOV2008]. In the present formulation, the contact region is pulled back to the reference configuration of the slave contact surface , henceforth, all the contact integrals will be computed on . For the unbiased treatment of the contacting bodies, we refer to work by Sauer and De Lorenzis [unbiased2013, unbiased2015].
The minimum normal gap is defined as
(1) 
Here, is a point on the slave surface , is the corresponding physical coordinate of the closest contact point on the master surface and is the unit normal vector at .
The tangential relative gap is defined in incremental form as
(2) 
which is in agreement to the timediscretized backward Euler formulation. Here, as well in the following, all the contact quantities will be expressed by default to the current time step , while the subscript is used to define quantities at previous step .
The contact traction vector in the current configuration is defined in terms of its normal and tangential components as
(3) 
In case of , the normal contact tractions are activated that avoid the penetration of the contact regions. In order to enforce the impenetrability in the normal direction, following KarushKuhnTucker (KKT) conditions must be satisfied [Wriggers2006]
(4) 
and the KKT conditions for Coulomb friction law in the tangential direction are [Wriggers2006]
(5) 
where and are the Coulomb friction coefficient and tangential slip velocity, respectively.
It is noted that the abovedescribed contact constraints cannot be directly incorporated into the variational form as they lead to the nonsmooth relationship between the normal gap and the contact pressure, resulting in a nonsmooth normal contact constitutive law. In the present work, the penalty method is adopted for the regularization of the contact constraints defined in Eqs. (4) and (5). The penalty regularized normal contact constraint is defined as
(6) 
where denotes the normal penalty parameter. The regularized frictional contact constraints in the tangential direction yields:
(7) 
where is the tangential penalty parameter. The frictional contact traction is determined based on the classical returnmapping algorithm. The trial tangential traction and contact state are first computed with the assumption that
(8) 
and the status: slip or stick is identified using
(9) 
Based on the penalty solution method, the contact contribution to the weak form (contact virtual work) is given by
(10) 
For the evaluation of above contact integral an active set strategy is utilized. The linearization of the contact virtual work , which is necessary for the NewtonRaphson iterative solution method, leads to
(11) 
Further details, together with the variation and linearization of all quantities such as , , etc., can be found in [DeLorenzis2011, Wriggers2006].
2.2 NURBSbased discretization
In this subsection, the existing NURBSbased discretizing procedure used for the continuum and the contact surface in the context of IGA is briefly discussed.
The IGA has emerged as a successful integration of computeraideddesign (CAD) and finite element analysis (FEA) [Cottrell2009]. It utilizes the CAD polynomials as a basis for the modelling of complex geometries exactly and the approximation of unknown solution fields. In the present work, NURBS are used for the discretization of the continuum including the contact surface.
NURBS are built from Bsplines. For a given knot vector along the parametric direction, the Bspline interpolations of order are defined using the following Coxde Boor recursive relation [nurbsbook]:
(13)  
(14)  
where and is the knot of the knot vector which is a set of nondecreasing value of parametric coordinates, i.e. . Also, denotes the total number of control points for order of Bspline functions along the parametric direction. A order Bspline function offers continuity across each knot , where is the knot multiplicity. In an open knot vector, the first and last knot entries are repeated by times that make the functions interpolatory at the end knots.
NURBS are the projective transformation of Bsplines. Thus, NURBS additionally utilize the weight values in their construction. For a given knot vector , a order univariate NURBS function along the parametric direction is defined as [nurbsbook]
(15) 
where and is the weight value. For a specified control points vector, a order NURBS curve along the parametric direction can be constructed using the linear combination of univariate NURBS interpolations and the control point’s coordinate vector as
(16) 
A bivariate continuum patch is defined by the tensor product of the two univariate NURBS curves defined along the and parametric directions as
(17) 
where is the coordinate vector of the control points net . Here, is the bivariate NURBS basis function defined as the tensor product of two univariate Bspline basis functions along the and parametric directions as
(18) 
with as a normalized weight that is given in terms of weight point and Bspline basis functions. A composition of knot vector with associated control points accompanied by weights constitutes a patch. For in detailed description, we refer to monograph by Pigel and Tiller [nurbsbook] and Cottrell et al. [Cottrell2009].
In the context of IGA, a NURBS discretized geometry can be refined by means of knot insertion (refinement), order elevation (refinement), and refinement strategies [hughes2005]. Within the knot insertion scheme, an additional knot in the knot vector is inserted. If the inserted knot value is unique, an additional knotspan, consequently an element, without changing the interelement continuity of the NURBS, is introduced. On the other hand, repetition of knots reduce the smoothness of the NURBS basis functions.
In case of order elevation scheme, the geometry is refined by means of raising the order of the NURBS interpolations. For the elevation of order from to , each knot value is repeated by times. As a result, the continuity of the NURBS remains unchanged during refinement.
The particular order of application of order elevation and knot insertion schemes yields refinement [hughes2005]. During refinement, first, the interelement continuity of the NURBS interpolations is increased and then, the additional elements within the given knot vector are introduced. On the other hand, if their application order is reversed, i.e. the knot insertion is performed first before the order elevation scheme [Corbett2014], the interelement continuity of the NURBS remains unchanged and knot values are repeated due to orderelevation. This, as a result, yields a large number of control points compared to refinement strategy.
Within the existing standard NURBSbased discretization approach, to improve the accuracy of the contact solution with a fixed mesh, the interpolation order of the NURBS discretized structure is uniformly elevated. However, this approach is not computationally favorable since the higherorder NURBS used for the approximation of contact integrals have to be employed for the computation of the vast majority of the bulk region away from the contact surface. It therefore becomes desirable to develop an improved NURBSbased discretization method that provides a possibility to perform controllable orderelevation for a NURBS discretized geometry as already suggested in the paper by Temizer et al. [Temizer2012]. Such a method is proposed in this paper and is described in detail in the next Section.
3 Varyingorder NURBSbased discretization
In this section, we introduce the theory for the varyingorder NURBS discretization method, followed by the description on its integration into an existing isogeometric contact code.
3.1 Varyingorder NURBS
The basic concept of the varyingorder (VO) NURBS discretization for isogeometric contact analysis is illustrated in Fig. 1. Consider a body having contact boundary layer . Let and be the minimum orders of NURBS basis functions capable of representing the geometry in an exact manner and a coarse mesh is given by the product of open knot vectors defined along the and parametric directions as shown in Fig. (a)a. Next, in order to make use of higherorder NURBS interpolation functions for contact computations, the discretized contact boundary is replaced with a NURBS layer of order as shown in Fig. (b)b.
The NURBS layer is constructed using refinement or through a combination of refinement and orderelevation [Corbett2014] strategies in such a manner that it matches the bulk parametrization. The application of the refinement strategy to the NURBS contact layer increases the order as well as the interelement continuity of the NURBS basis functions. On the other hand, application of one or two additional steps of orderelevation to refined NURBS layer introduces a large number of additional control points along the contact boundary layer while the interelement continuity, i.e. , remains unchanged. This way, the higherorder NURBS basis functions are utilized for the evaluation of contact contributions, a large number of additional degrees of freedom are introduced across the contact surface with a fixed mesh, and the minimum order NURBS are employed for the description of the bulk domain that does not come into contact. The resulting isogeometric contact element is characterized by number of control points, where on the contact layer and on the remaining three faces of each element. The bivariate NURBS basis functions for such an element are defined as
(19)  
where the normalizing weight is given by
(20) 
The basis functions in Eq. (3.1) exhibit the nonnegativity property
(21) 
and they also satisfy the partition of unity property:
(22) 
Using the isoparametric concept, the basis functions defined in Eq. (3.1) are employed for the approximation of the unknown displacement field, its variation , and the current coordinates as
(23) 
where, , and are the displacement, its variation and the current coordinate vector corresponding to the control point, respectively.
The proposed VO NURBS based discretization approach generalizes the enrichment strategy of Corbett and Sauer [Corbett2014], which can be constructed using the current methodology by discretizing the region other than the contact boundary layer with the linearorder NURBS interpolations. Moreover, unlike the enrichment strategy of [Corbett2014], the current discretization approach does not require Bézier extraction operator [Bezier_extraction_NURBS], which is required to enable the incorporation of NURBS basis into the finite element code structure.
3.2 Implementation into existing code
For integrating the VO NURBS discretization strategy into the existing isogeometric contact code, only minor modifications are required. First of all, for a given mesh resolution, a order of NURBS curve representing the contact boundary layer of the body is constructed. After that, the parametrization for the initially described order NURBS contact layer is replaced with that of the newly constructed order curve. For this, a certain number of conditions are needed to be fulfilled. The total number of control points defining the body must be updated to allow the incorporation of the order contact curve. This means that the connectivity array for contact elements must be adapted in a manner that it contains all the underlying control points. The contact element connectivity arrays can have a different length than the bulk element connectivity arrays. The bivariate NURBS basis functions for contact elements are evaluated at the quadrature points according to Eq. (3.1). Within each contact element, number of GaussLegendre quadrature points are employed for the evaluation of the weak form obtained through the NURBS interpolation. This ensures the quadratic convergencerate during the NewtonRaphson iterations. Optimal quadrature rules [Hughes2010, AURICCHIO2012, FAHRENDORF2018], which are wellsuited for IGA, can also be opted for the reduced numerical evaluation. With the exception of these modifications, no other changes are needed to be made in an existing code. The local quantities, e.g. stiffness matrices and force vector, are assembled to their global part in the same way as with the standard procedure. The reader is referred to Agrawal and Gautam [Agrawal2018] for detailed description on the implementation of IGA in a simplified manner.
4 Numerical Examples
In this section, we show the capabilities and performance of the proposed VO NURBS discretization method over the existing standard NURBS discretization approach using three, twodimensional small and large deformation contact problems. In the first example, the contact patch test is considered to test the validity of the mortar contact formulation for VO NURBS discretizations. Next, the numerical results and the analytical solution of the Hertz contact problem are compared on the basis of accuracy and overall processing time for VO and standard NURBS discretizations. In the final example, a large deformation frictional ironing problem is solved to demonstrate the superiority of the proposed discretization method over the standard approach in terms of accuracy, efficiency, and robustness.
4.1 Contact patch test
The first example considers a simplified form of the contact patch test originally proposed by Taylor and Papadopoulos [Taylor1991]. The purpose of this example is to verify whether the assemblage of displacement based isogeometric contact elements for VO NURBS discretization is complete. The contact algorithm that passes the contact patch test ensures the convergence of contact solution upon the mesh refinement. Otherwise, solution errors at the contact interface do not necessarily reduce with decreasing the mesh size [Bathe2001].
The setup of the patch test, taken from Matzen and Bischoff [Matzen2016], is illustrated in Fig. 2. In this, two blocks, which are modelled using isotropic linear elastic material with Young’s modulus N/mm and Poisson’s ratio under the plane stress condition, are positioned on top of each other. The top line of the upper block is subjected to prescribed vertical displacement mm. The bottom line of the lower block is fixed against vertical displacement. For the considered displacementbased loading condition, the contact line is expected to deform uniformly by mm.
To deal with the worst case scenario, the nonconforming mesh at the contact interface is considered as shown in Fig. 2. For the bulk, second order NURBS are used and the contact boundary is given by , and order of the NURBS layer. For the purpose of comparison, quadratic order of NURBS discretization, i.e. N, is employed.
Figure 3 shows the contour plot of vertical displacement for NN discretization. It can be observed that a uniform gradient of develops between the top and the bottom lines of the setup that ranges from to mm. The gradient of remains constant in the horizontal direction. Identical displacement contour plots have been obtained for NN and N discretizations.
Figure 4 illustrates the deviation of numerically evaluated vertical displacement along the contact interface from the exact solution, i.e. mm, for N, NN, and NN discretizations. It can be observed that the deviation of vertical displacement from analytical solution reduces on increasing the order of NURBS interpolations for contact boundary from to , and . This is due to employing higher order (, and ) NURBS for the evaluation of contact integrals, which in turn improves the accuracy of the result compared to N discretization with a fixed mesh. It is noted that with standard N and N discretizations identical to NN and NN results are achieved since same order of NURBS functions are used for the contact computations in both the discretization approaches. This supports the proposition that the accuracy of the solution improves on increasing the order of NURBS interpolations for the discretization of the contact surface.
Next, to investigate the convergence behaviour of the mortar contact formulation with VO NURBS discretization, the relative vertical displacement error using the norm with successively refining the coarsest mesh shown in Fig. 2 is reported. Figure 5 shows the convergence plots for N, NN, and NN discretizations. It can be observed that the error norm of decreases and approaches toward the exact solution as the mesh is refined. The rate of convergence for each case of discretization is one. For a fixed mesh resolution, NN provides the most accurate result followed by NN compared to N discretization. Quantitatively, the value of relative errors even with the coarsest mesh is below for NN and NN discretizations, which are probably within the acceptable margins. With the obtained results, it is remarked that the isogeometric mortar contact formulation with VO NURBS is numerically consistent as the obtained result converges to the exact solution with refining the mesh.
4.2 Hertz contact problem
In this example, the frictionless contact of an infinitely long cylinder having outer radius and inner radius with a rigid surface is considered. The setup of the problem, adopted from Temizer et al. [Temizer2011], is illustrated in Fig. 6. Due to symmetry, only a quarter of the geometry is considered. The cylinder is subjected to the overall force applied as a distributed vertical load over the top surface of the cylinder. A linearly elastic material with Young’s modulus and Poisson’s ration under the plane strain condition is used for the analysis. The penalty parameter is chosen as a default value to prevent the penetration of the cylinder into the rigid surface.
Four different mesh arrangements having , , , and number of elements along each parametric directions are considered. In order to capture the variation of the contact force along the contact interface, approximately of the total elements in each parametric direction are redistributed in such a manner that they lie within the surface length of the geometry [Temizer2011]. First, we compare the results for NN discretizations, where the quadratic order NURBS are used to discretize the bulk domain and , and order NURBS are used for the contact boundary layer, with that of the standard N order of NURBS discretizations. Figure 6 illustrates the obtained results for both the discretization approaches with different mesh arrangements. The dimensionless contact pressure is plotted versus the dimensionless contact coordinate , where is the normal contact pressure evaluated at the active control points and is the distance of these points from the first point of contact. The analytical solution for the current setup, i.e. the maximum contact pressure and the contact area radius are given by the Hertz theory [johnson1987].








From Fig. 6, it clear that the contact pressure distributions for NN are in excellent agreement with their corresponding standard N order of NURBS discretizations for the same mesh level. This is due to employing the same order of NURBS interpolations for the evaluation of contact integrals in both the discretization approaches. Further, the accuracy of the contact solution increases monotonically on increasing the mesh resolution and nearly exact result is obtained with a very fine mesh arrangement, see Figs. (g)g and (h)h. These results closely match those reported by Temizer et al. [Temizer2011] and Lorenzis et al. [DeLorenzis2011] for the classical Hertz contact problem.
Besides comparing the accuracy of the obtained results for VO NURBS to standard NURBS discretization approach, the associated computational cost in terms of the overall analysis time is also investigated. Figure 7 shows the required time for pre and postprocessing (i/o), contact and contact elements computations, and the bulk element computations for both the discretization approaches. The analysis time associated with each discretization is calculated using
(24) 
where the maximum total analysis time is for the standard N order of NURBS discretization at the very fine mesh level which is the most expensive.
From Figs. (a)a and (b)b, it can be observed that NN utilizes at least lower computational cost than that of corresponding standard discretization with the fixed mesh resolution. It is impressive as NN take at most half the computational cost than N to deliver identical results. The major difference between the overall analysis time for the two discretization methods lies within their bulk element computations. The VO NURBS discretization method employs lowerorder NURBS for the bulk description. This, as a result, uses a fewer number of degrees of freedom per element during the bulk element computations compared to corresponding standard N order of NURBS discretization. The degree of freedom density data for both the discretization approaches with different mesh arrangements is given in Table 1.


Mesh  Discretizations  DOFs  

Interface  Bulk  Total  
N  150  7200  7350  
NN  150  7104  7254  
N  152  7296  7448  
NN  152  7104  7256  
N  222  10656  10787  
NN  222  10560  10782  
N  224  10752  10976  
NN  224  10560  10784  
N  294  14112  14406  
NN  294  14016  14310  
N  296  14208  14504  
NN  296  14016  14312  
N  366  14568  17934  
NN  366  17472  17838  
N  368  17664  18032  
NN  368  17572  17840  
NN  292  7104  7396  
NN  436  7104  7540 


From Figs. (a)a to (f)f, it can be further observed that the distribution of contact pressure curve matches exactly with the analytical solution except at the boundary of contact region. It is affected by the element that lies between the contact and noncontacting region. This is a well known issue for which a number of solution approaches have been introduced, e.g. see [Franke2010, Duong2018]. Within the context of IGA, such oscillation error is alleviated if a very fine mesh is used as shown in Figs. (g)g and (h)h. The present analysis particularly focuses on the use of higherorder NURBS to improve the approximation of the contact pressure at a very coarse mesh. For this, an additional orderelevation refinement scheme over the NURBS discretized contact layer is applied. Performing one and two steps of additional order elevation on NN results in NN and NN discretizations, respectively. With this, a large number of additional degrees of freedom across the contact interface are introduced due to the repetition of knots during order elevation. The total numbers of degrees of freedom across the contact interface and in the bulk domain for both the discretizations are also listed in Table 1.
Figure 8 illustrates the results for NN and NN with mesh arrangement. The associated analysis time for these two discretizations are shown in Fig. 9.
From Fig. (b)b, it can be observed that NN discretization with mesh achieves the accuracy comparable to standard N order of NURBS discretization based results with a very fine mesh , see Figs. (g)g and (h)h. This is due to a large number of degrees of freedom present across the contact interface with NN which majorly improves the approximation of contact responses with a very coarse mesh. Further, from Fig. 9, it is evident that NN takes approximately lower analysis time than that with standard N discretization, which is remarkable. The other VO NURBS discretization, e.g. NN and NN, deliver very similar results to . However, they utilize slightly higher computational cost as compared to NN.
4.3 Frictional ironing problem
The third example considers the ironing problem with two deformable bodies [DeLorenzis2011]. This example is used to demonstrate the performance of VO NURBS discretization in terms of accuracy, robustness, and cost of solution with respect to standard NURBS discretization approach for large deformation and large relative tangential motion with Coulomb friction. The geometric model along with the material details and the boundary condition are shown in Fig. 10. Three nested meshes, which are shown in Fig. 11 and are denoted as , , and , respectively, are used for the die and the slab.
In this example, the die is first pressed downwards into the slab by uniform vertical displacement mm in load steps, and then moved horizontally across the slab by applying a displacement mm in load steps. An isotropic, hyperelastic NeoHookean material model is considered for both the die and the slab for large deformations. The corresponding constitutive relation is given by [bonet1997]
(25) 
where is the deformation gradient, , and is the identity tensor. In this work, the values for the Lamé constants are defined as and . The penalty parameters are taken as and the coefficient of sliding friction is . The mortar contact integrals are evaluated using equidistant quadrature points per contact element.



4.3.1 Prediction of contact reaction forces
Figure 12 illustrates the horizontal contact reaction force and the vertical contact reaction force as a function of load step for VO NN and the standard N order of NURBS discretizations with mesh . The results are in agreement with those reported by De Lorenzis et al. [DeLorenzis2011] for the sliding case.



The enlarged views for and are shown in Figs. (b)b and (c)c, respectively. It can be observed that the nonphysical oscillations of vertical and horizontal contact reaction forces are present for all the cases during sliding. The oscillation error reduces with increasing the order of NURBS discretization, i.e. from N to N. The curves of and are indistinguishable for NN and standard N order of NURBS discretizations as expected. This is due to employing the same order of NURBS for the evaluation of contact integrals and identical number of degrees of freedom across the contact boundary in both the cases of discretizations. Table 2 lists degrees of freedom density data for both NN and N order of NURBS discretizations. With NN and NN, the oscillation amplitude of horizontal contact reaction force reduce to approximately and , respectively, of that observed with N. A quantitative analysis of the reduction in the oscillation error for NN and N is given in Table 3 for both the force components. The oscillation amplitude of reaction forces is computed using . The associated total computational time for the two discretization approaches is reported in Sec. 4.3.2.
We now present the results for NN and NN discretizations, where one and two steps of additional order elevation scheme are applied to the N discretized contact boundary layer. Figure 13 shows the variation of the horizontal and the vertical contact reaction forces for NN and NN with respect to most accurate result with the standard N order of NURBS. It can be noticed that the oscillation error reduces significantly for the new discretizations. Among the two tested cases, NN delivers the best result. It reduces the oscillation amplitude of vertical and horizontal reaction forces to and , compared to N. This is attributed to a large number of additional degrees of freedom present across the contact interface when compared to other cases (see Table 2 for degree of freedom details). The other higherorder NURBS layer based discretizations, e.g. NN and NN, also reduce the oscillation error more than N, although equivalent to NN results are obtained as shown in Figs. (b)b and (a)a for NN arrangement.
It should be noted that with the most accurate enrichment strategy of Corbett and Sauer [Corbett2014], i.e. QN, the maximum reductions that are achieved in the oscillation amplitude of the horizontal and the vertical contact reaction forces are and , respectively, of that observed with quadratic NURBSenriched finite element QN for the similar sliding case. On comparing the results for NN with that of QN, it is evident that the VO NURBS based discretization NN reduces the oscillation amplitude of horizontal and the vertical contact reaction forces by a factor of approximately and , respectively, of that achieved with QN.
Element  DOF for die  DOF for slab  Total  

Interface  Bulk  Total  Interface  Bulk  Total  DOF  
N  30  120  150  28  168  196  346 
N  34  136  170  30  180  210  380 
N  38  152  190  32  192  224  414 
NN  34  120  154  30  168  198  352 
NN  38  120  158  32  168  200  358 
NN  54  120  174  52  168  220  394 
NN  78  120  198  76  168  244  442 
NN  82  120  202  78  168  246  448 
Element  

N  100  100 
N  60  44.12 
N  32  23.53 
NN  62  44.12 
NN  32  26.47 
NN  7.1  14.71 
NN  4.4  8.82 
NN  4.4  8.82 
Such a difference between the performance of NN and QN is due to the fact that the accuracy of contact reaction forces does not only depends on the discretization of the contact surface, but also on the bulk deformations as shown in [Corbett2014, Corbett2015]. In contrast to quadraticorder NURBS N based mesh, bilinear Lagrange finite elements Q based mesh suffers from locking. Thus, NN leads to more accurate results than QN.
With the obtained results it is clear that the VO NURBS discretization, particularly NN, is a more effective strategy than the standard NURBS based discretization approach in the context of IGA and the existing enrichment strategy of [Corbett2014] for this case.
Figure 14 shows the contour plot of vertical displacement field in the deformed configurations of the setup at different load steps for N order of NURBS discretization with mesh arrangement.





4.3.2 Computational time
Next, we investigate the reduction in the oscillation amplitude of and in terms of associated overall computational time upon refining the mesh for the die and slab shown in Fig (a)a. Figure 15 illustrates the results for different VO and standard NURBS discretizations with three meshes , , and . The associated total analysis time for each discretization is calculated using Eq. (24), where the time for N order of NURBS at a fine mesh is used as a reference.
From Figs. (a)a and (c)c it is evident that in order to attain the accuracy comparable to N discretization, NN takes at least lower computational time as compared to standard N discretization for the same mesh level. Further, from Figs. (b)b and (d)d it is clear that a significant improvement in the accuracy is achieved with NN and NN for the fixed mesh resolution. It turns out that NN with the mesh and yield results comparable to N discretizations with the mesh level and , respectively. The analysis time for NN is still considerably lower, i.e. , than that for N discretization. It can also be observed that a good convergence behaviour with VO NURBS is achieved.




5 Conclusion
In the present work, a novel varyingorder NURBS discretization method is proposed that provides a possibility to achieve controllable order elevation in the framework of isogeometric analysis. The method makes use of higherorder NURBS polynomials for the evaluation of the contact integrals. Lowerorder of NURBS capable of representing the geometry exactly are used for the bulk computations. To achieve this, a higherorder NURBS layer, which accompanies a large number of degrees of freedom and matches with the bulk parametrization, is utilized as the contact boundary layer of the geometry. The isogeometric contact formulation with varying order NURBS discretization is presented in the nonlinear kinematic setting using the mortar method. Nonpenetrability and tangential sticking constraints are regularized using the penalty method.
The consistency of the presented isogeometric mortar contact formulation with VO NURBS discretization is verified through the contact patch test. The performance and capabilities of the proposed discretization method are then demonstrated through twodimensional frictionless and frictional contact problems, considering small and large deformations. Based on the obtained results, it can be concluded that the newly proposed varyingorder NURBS discretization method displays a superior performance in terms of accuracy and efficiency compared to that of the standard NURBS discretization approach.
For the same mesh resolution, the proposed method requires a considerably lower computational cost to obtain identical results as that of standard NURBS discretization. The gain in the efficiency stems from the lowerorder NURBS polynomials employed for the computation of the vast majority of the bulk domain that does not come into contact. This gain further improves as the mesh is refined.
It is next shown that the accuracy of the solution significantly improves on applying the additional steps of order elevation to NURBS contact layer. This is attributed to the presence of a large number of additional degrees of freedom across the contact interface which consequently improves the accuracy of contact computations even with a coarse mesh. It is also shown that the associated computational cost still remains considerably lower as compared to the standard NURBS discretization.
Refining the geometry using varying order NURBS is observed to be quite robust for the considered numerical examples. Additionally, a good convergence behaviour is obtained with the proposed discretization method. Simplicity of the current approach allows itself to be conveniently embedded into existing isogeometric contact codes.
Several extensions are considered in the near future. The extension of the proposed VO NURBS discretization method to threedimensional contact problems within the context of isogeometric analysis is planned. An efficient quadrature technique [Duong2015] that adaptively refines the integration domain along the contact boundary of the geometries to further improve the integration accuracy of the contact forces and their tangent contributions with a very coarse mesh can also be incorporated in the current framework.
6 Acknowledgment
The authors gratefully acknowledge the support from SERB, DST under project SR/FTP/ETA0008/2014. The authors would like to express their gratitude towards Prof. Roger A. Sauer, AICES RWTH Aachen University for his valuable comments.