Varying-order NURBS discretization: An accurate and efficient method for isogeometric analysis of large deformation contact problems

TL;DR

This paper introduces a varying order NURBS discretization method that improves the accuracy and efficiency of isogeometric analysis in large deformation contact problems by using higher-order NURBS selectively at contact boundaries.

Contribution

The paper proposes a novel varying order NURBS discretization approach that enhances isogeometric contact analysis by combining different NURBS orders for geometry and contact evaluation.

Findings

Achieves accurate solutions with coarse meshes.

Reduces computational cost compared to standard methods.

Maintains robustness and stability in 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 higher-order NURBS for contact integral evaluations. Lower-orders 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 higher-order 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 varying-order NURBS discretization is first examined through the contact patch test. The capabilities and benefits of the proposed methodology are then demonstrated in detail using two-dimensional 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.

Figures

Tables

Table 1. Table 1: Number of elements in different meshes used for ironing problem.

Mesh				Elements
				Die			Slab
m1				$6\times 2$			$8\times 3$
m2				$12\times 4$			$16\times 6$
m3				$24\times 8$			$32\times 12$
m4				$48\times 16$			$64\times 24$

Table 2. Table 2: Degrees of freedom (DOF) density data for the die and slab for different VO and standard NURBS based discretizations at mesh m 1 subscript m 1 \textrm{m}_{1} .

Element	DOF for die			DOF for slab			Total
	Interface	Bulk	Total	Interface	Bulk	Total	DOF
N2	18	36	54	20	60	80	134
N4	26	52	78	24	72	96	174
N6	34	68	102	28	84	112	214
N${}_2{}^{}-$N4	26	36	62	24	60	84	146
N${}_2{}^{}-$N6	24	36	70	28	60	88	158
N${}_2{}^{}-$N2⋅1	30	36	66	36	60	96	162
N${}_2{}^{}-$N2⋅2	42	36	78	52	60	112	190
N${}_2{}^{}-$N2⋅3	54	36	90	68	60	128	218

Table 3. Table 3: Reduction in the oscillation amplitude of vertical and horizontal contact forces for different standard and VO based NURBS discretizations at mesh m 1 subscript m 1 \textrm{m}_{1} . The oscillation amplitude for N 2 is used as a reference.

Element		$\Delta {\text{P}}_y(\%)$		$\Delta {\text{P}}_x(\%)$
N2		100		100
N4		49.88		23.31
N6		24.57		11.81
N${}_2{}^{}-$N4		48.11		27.02
N${}_2{}^{}-$N6		23.83		12.09
N${}_2{}^{}-$N2⋅1		5.98		30.84
N${}_2{}^{}-$N2⋅2		4.31		7.51
N${}_2{}^{}-$N2⋅3		1.89		4.77

Table 4. Table 4: Degrees of freedom density data for both the discretization approaches with different mesh arrangements.

Mesh	Discretizations	DOFs
		Interface		Bulk			Total
	L1	20		960			980
m1	N2	22		1056			1078
	N${}_2{}^{}-$N2⋅1	40		1056			1096
	N${}_2{}^{}-$N2⋅2	58		1056			1114
	L1	38		1824			1862
m2	N2	40		1920			1960
	N${}_2{}^{}-$N2⋅1	76		1920			1996
	N${}_2{}^{}-$N2⋅2	112		1920			2032
	L1	74		3552			3626
m3	N2	76		3648			3724
	N${}_2{}^{}-$N2⋅1	148		3648			3796
	N${}_2{}^{}-$N2⋅2	220		3648			3868
	L1	146		7008			7154
m4	N2	148		7104			7252
	N${}_2{}^{}-$N2⋅1	292		7104			7396
	N${}_2{}^{}-$N2⋅2	436		7104			7540
	L1	290		13920			14210
m5	N2	292		14016			14308
	N${}_2{}^{}-$N2⋅1	580		14016			14596
	N${}_2{}^{}-$N2⋅2	868		14016			14884
m6	N2	580		27840			28420

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} -norm of the error in p N and p T for different discretizations and meshes.

Discretization	${\Vert {\text{p}}_{N_{ref}}-{\text{p}}_N\Vert }_2$			${\Vert {\text{p}}_{T_{ref}}-{\text{p}}_T\Vert }_2$
	m1	m2		m1	m2
N2	$24.38\times {10}^{-2}$	$34.32\times {10}^{-3}$		$24.38\times {10}^{-3}$	$34.33\times {10}^{-4}$
N4	$17.45\times {10}^{-2}$	$25.78\times {10}^{-3}$		$17.45\times {10}^{-3}$	$25.99\times {10}^{-4}$
N${}_2{}^{}-$N4	$17.18\times {10}^{-2}$	$26.80\times {10}^{-3}$		$17.18\times {10}^{-3}$	$27.01\times {10}^{-4}$
N${}_2{}^{}-$N2⋅1	$7.095\times {10}^{-2}$	$3.433\times {10}^{-3}$		$7.097\times {10}^{-3}$	$4.429\times {10}^{-4}$

Table 6. Table 6: Overall analysis time taken by various discretizations with different meshes. The time with N 2 using mesh m 1 is used for normalization.

Discretization		Time
		m1	m2
N2		$1.000$	$1.357$
N4		$1.451$	$2.296$
N${}_2{}^{}-$N4		$1.216$	$1.724$
N${}_2{}^{}-$N2⋅1		$1.058$	$1.475$

Equations

$$\text{for pi=0, }{N}_{l,0}({\xi }^i)$$

$$\text{for pi>0, }{N}_{l,p_i}({\xi }^i)$$

$$R_l^{p_i}({\xi }^i)=\frac{w_l^i}{{\sum }_{A=1}^{n_i}{w}_A^i{N}_{A,p_i}({\xi }^i)}{N}_{l,p_i}({\xi }^i),$$

$$\mathbf{C}({\xi }^i)=\sum _{l=1}^{n_i}{R}_l^{p_i}({\xi }^i){\mathbf{X}}_l^i.$$

$$\mathbf{S}({\xi }^1,{\xi }^2)=\sum _{l=1}^{n_1}\sum _{m=1}^{n_2}{R}_{lm}^{p_1,p_2}({\xi }^1,{\xi }^2){\mathbf{X}}_{lm},$$

$$R_{lm}^{p_1,p_2}({\xi }^1,{\xi }^2)=\frac{w_{lm}}{W({\xi }^1,{\xi }^2)}{N}_{l,p_1}({\xi }^1)N_{m,p_2}({\xi }^2),$$

$${\text{uptau}}_1({\xi }^{\mathrm{m}})=\frac{\partial {\mathbf{x}}^{\text{m}}}{\partial {\xi }^{\text{m}}}={\mathbf{x}}_{,{\xi }^{\text{m}}}^{\text{m}}.$$

$$\delta ({\Pi }_{\text{int}}+{\Pi }_{\text{c}}-{\Pi }_{\text{ext}})=0\iff \delta {\Pi }_{\text{int}}+\delta {\Pi }_{\text{c}}-\delta {\Pi }_{\text{ext}}=0,\forall \delta \mathbf{u}\in \mathcal{V}.$$

$$\delta {\Pi }_{\text{int}}$$

$$\delta {\Pi }_{\text{c}}$$

$$\frac{d}{d{\xi }^{\text{m}}}\hat{d}({\xi }^{\text{m}})=[{\mathbf{x}}^{\text{s}}-{\mathbf{x}}^{\text{m}}({\xi }^{\text{m}})]\cdot {\text{uptau}}_1({\xi }^{\text{m}}).$$

$$\frac{d^2}{d({\xi }^{\text{m}}{)}^2}\hat{d}({\xi }^{\text{m}})=[{\mathbf{x}}^{\text{s}}-{\mathbf{x}}^{\text{m}}({\xi }^{\text{m}})]\cdot {\mathbf{x}}_{,{\xi }^{\text{m}}{\xi }^{\text{m}}}^{\text{m}}({\xi }^{\text{m}})-{\text{uptau}}_1\cdot {\text{uptau}}_1.$$

$${\text{g}}_{\text{N}}=({\mathbf{x}}^{\text{s}}-{\bar{\mathbf{x}}}^{\text{m}})\cdot \bar{\mathbf{n}},$$

$${\text{g}}_{\text{N}}\ge 0,{\text{t}}_{\text{N}}\le 0,{\text{g}}_{\text{N}}{\mathrm{t}}_{\text{N}}=0,\text{on}{\Gamma }_{\mathrm{c}}^s.$$

$$\bm{t}_{\text{N}}=\left\{\begin{array}[]{lr}-\epsilon_{\text{N}}\,\textrm{g}_{\textrm{N}}{\bm{\bar{n}}},&\textrm{if}~{}\textrm{g}_{\textrm{N}}<0\\ \bm{0},&\textrm{otherwise}\,,\end{array}\right.$$

$$\mathbf{t}={\text{t}}_{\text{N}}\bar{\mathbf{n}}-{\text{t}}_{\text{T}}{\overline{\text{uptau}}}_1$$

$${\mathbf{t}}_{\text{T}}=-{\mu }_f{\text{t}}_{\text{N}}\frac{{\dot{\text{g}}}_{\text{T}}}{||{\dot{\text{g}}}_{\text{T}}||}$$

$${\mathbf{g}}_{\text{T}}={\mathbf{g}}_{\text{T}e}+{\mathbf{g}}_{\text{T}sl}.$$

$${\dot{\mathbf{g}}}_{\text{T}sl}=\gamma \frac{{\mathbf{t}}_{\text{T}}}{||{\mathbf{t}}_{\text{T}}||}.$$

$$\gamma \ge 0,\Phi \le 0,\gamma \Phi =0.$$

$${\mathbf{t}}_{\text{T}n+1}^{\text{trial}}$$

$${\Phi }_{n+1}^{\text{trial}}$$

$$\bm{t}_{\textrm{T}}=\left\{\begin{array}[]{lr}\bm{t}_{\text{T}n+1}^{\textrm{trial}},{}{}&\textrm{if}~{}\Phi_{n+1}^{\textrm{trial}}\leq 0\,,\\ \mu_{f}\,\textrm{t}_{\textrm{N}n+1}\bm{n}_{\textrm{T}n+1},{}{}&\textrm{otherwise}\,.\end{array}\right.$$

$${\mathbf{u}}^{ke}={\mathbf{R}}^k({\mathbf{\xi }}^k){\mathbf{u}}^{ke},\delta {\mathbf{u}}^{ke}={\mathbf{R}}^k({\mathbf{\xi }}^k)\delta {\mathbf{u}}^{ke},{\mathbf{x}}^{ke}={\mathbf{R}}^k({\mathbf{\xi }}^k){\mathbf{x}}^{ke}\forall {\mathbf{\xi }}^k\in {\Omega }^{ke}.$$

$$\delta {\mathbf{u}}^{\text{T}}[{\mathbf{f}}_{\text{int}}+{\mathbf{f}}_{\text{c}}]=0,\forall \delta \mathbf{u}\in \mathcal{V}$$

$${\mathbf{f}}_{\text{int}}^{ke}={\int }_{{\Omega }^{ke}}[{\mathbf{B}}^k{]}^{\text{T}}{\mathbf{\sigma }}^k\text{d}v,$$

$$\mathbf{\sigma }=\frac{\lambda }{J}(\text{ln}J)\mathbf{I}+\frac{\mu }{J}({\mathbf{FF}}^{\text{T}}-\mathbf{I}),$$

$${\mathbf{f}}_{\text{c}}^{\text{s}e}=-{\int }_{{\Gamma }_c^{\text{s}e}}[{\mathbf{N}}^{\text{s}}({\xi }^{\text{s}}){]}^{\text{T}}{\mathbf{t}}^{\text{s}}\text{d}\Gamma ,{\mathbf{f}}_{\text{c}}^{\text{m}e}={\int }_{{\Gamma }_c^{\text{s}e}}[{\mathbf{N}}^{\text{m}}({\bar{\xi }}^{\text{m}}){]}^{\text{T}}{\mathbf{t}}^{\text{s}}\text{d}\Gamma ,$$

$$R_1^{p_c,p_2}({\xi }^1,{\xi }^2)$$

$$⋮$$

$$R_{p_c+1}^{p_c,p_2}({\xi }^1,{\xi }^2)$$

$$R_{p_c+2}^{p_1,p_2}({\xi }^1,{\xi }^2)$$

$$⋮$$

$$R_{n_{cp}^e}^{p_1,p_2}({\xi }^1,{\xi }^2)$$

Full text

The capabilities and benefits of the proposed method are demonstrated using the two-dimensional frictionless and frictional contact problems, considering both small and large deformations. The results with the existing fixed-order based NURBS discretizations are used for comparisons. Numerical examples show that with the proposed method, a much higher accuracy can be achieved even with a coarse mesh as compared to the existing NURBS discretization approach. It exhibits a major gain in the numerical efficiency without the loss of stability, robustness, and the intrinsic features of NURBS-based IGA technique for a similar accuracy level. The simplicity of the proposed 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; Higher-order contact boundary; Frictional contact; Non-linear continuum mechanics

1 Introduction

In the last two decades, non-uniform rational B-splines (NURBS) based isogeometric analysis [1] approach has been established as an advantageous computational technology for various classes of problems, especially to those where the geometric approximation and inaccurate capture of the driving forces at the interface adversely influence the accuracy of the result. For example, in cohesive zone modelling problems [2, 3, 4], estimation of the fatigue life of interfacial cracks in the functionally graded material plates [5, 6, 7], evaluation of through-thickness cracks in homogeneous and isotropic plates [8], in three-dimensional problems having straight and curved crack fronts [9], and for the analysis of cracks in the functionally graded piezoelectric materials [10], etc. 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 convex-hull properties, tailorable inter-element continuity, and non-negativeness [1, 11]. Contact modeling belongs to one of these classes that have particularly been benefited from the aforementioned features of the IGA technology. Higher-order smoothness of the NURBS discretized geometry directly provides the continuous normal vector field across the boundary of the contact elements. Thus, it eliminates the need for additional contact surface smoothening approaches [12, 13, 14, 15, 16, 17] often utilized in the context of the finite element method.

The first investigations on the treatment of contact problems using isogeometric analysis are conducted by Temizer et al. [18, 19], Lu [20], and De Lorenzis et al. [21, 22].

Temizer et al. [18] applied the isogeometric analysis to three-dimensional thermomechanical frictionless contact using the mortar contact formulation. For the regularization of the contact constraints, the penalty method is adopted. They demonstrated that NURBS based discretization achieves superior results in terms of quality and robustness over its counterpart Lagrange-polynomial based discretization. However, a very fine mesh is still required to get the result that closely matches with the exact solution. Moreover, the non-mortar contact formulation is found to be over-constrained in nature and leads to inaccurate results. Lu [20] combined the isogeometric analysis with the segment-to-segment based contact formulation presented by Papadopoulos and Taylor [23]. He showed that intrinsically smooth NURBS discretized geometry alleviates the non-physical oscillations of the contact forces and is capable of accurately describing the intricate mechanics of smooth materials such as fabrics. De Lorenzis et al. [21] introduced a mortar-based contact formulation in the context of IGA for the investigation of two-dimensional large deformation contact with Coulomb friction. They demonstrated that the magnitude of the non-physical oscillations of the normal and tangential contact reaction forces reduce on increasing the interpolation order of NURBS discretization. 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-polynomial based discretizations are found to be highly sensitive to the interpolation order. Later, De Lorenzis et al. [22] applied IGA to three-dimensional large deformation frictionless contact with the mortar method. The contact constraints are enforced using the augmented-Lagrangian approach [24]. Temizer et al. [19] introduced a three-dimensional mortar-based frictional contact formulation as an extension of works in [18, 21]. Again, it is shown that on increasing the interpolation order of the NURBS discretization the smoothness of contact forces improves across the contact interface. A point-to-segment contact formulation for two-dimensional frictionless contact is presented by Matzen et al. [25], and recently extended to a weighted point-based contact approach for three-dimensional contact [26], both with the IGA. Many researchers have extensively studied the application of mortar method in the context of IGA [27, 28, 29, 30, 31, 32], e.g. Brivadis et al. [30] investigated the isogeometric mortar method from the theoretical as well as the numerical point of view. Seitz et al. [31] developed a dual mortar method for isogeometric contact analysis and analyzed the spatial rate of convergence for the mesh tying and contact problems. Recently, Duong et al. [32] introduced a segmentation-free isogeometric mortar contact formulation. A number of works have explored the application of collocation methods to different contact problems using IGA, e.g. [33, 34, 35]. De Lorenzis et al. [36] presented a comprehensive overview on the growth of various isogeometric based treatment procedures and their advantages for contact mechanics as compared to traditional finite element analysis.

From the above-reviewed literature on contact modeling using IGA, it is evident that NURBS-based discretization delivers significantly superior performance in terms of accuracy, stability, and robustness over Lagrange-polynomial based discretization. This is due to the aforementioned distinguished intrinsic features of the NURBS over the Lagrange polynomials. However, the application of the existing NURBS-based discretization approach to isogeometric contact analysis is computationally expensive since due to the rigid tensor nature of the NURBS structures, the underlying mesh can be refined only 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 body can only be elevated uniformly. 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, T-splines [37, 38, 39, 40, 3, 4], NURBS-based hierarchical refinement [41, 42], and locally refined (LR) NURBS [43] approaches have been adopted. Dimitri et al. [40, 3, 4] employed the T-splines based discretizations for analyzing the frictionless contact, cohesive/contact interface modeling, and mixed-mode debonding problems. Temizer and Hesch [42], and Hesch et al. [41] used the hierarchical NURBS for the treatment of frictionless and frictional contact, respectively. Zimmermann and Sauer [43] introduced LR NURBS for the analysis of contact computations of solids and membranes. However, the idea to refine the geometry through a controllable order elevation strategy remains unexplored, as already noted by Temizer et al. [19].

In the context of finite element contact analysis, a considerable research progress has been made that integrates the intrinsic features of NURBS and higher-order Lagrange polynomials with the standard FE discretization [44, 45, 46, 47, 48, 49]. Corbett and Sauer [44, 45] introduced a NURBS-enriched contact element formulation that 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 “NURBS-enriched contact finite elements.” Their work is based on the contact element enrichment strategy of Sauer [50, 51], which used higher-order Lagrange and Hermite functions on the contact surface. Konyukhov and Schweizerhof [52, 53, 54] introduced the hierarchical enrichment of the Lagrange functions space, which allows the anisotropic refinement of the contact boundary layer with the higher-order Lagrange polynomial, while the linear order polynomials remain fixed for the construction of the interior domain of a contact body. Later, Maleki-Jebeli et al. [47] proposed a hybrid isogeometric-finite element discretization method where the advantages of NURBS are exploited in two-dimensional 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”, initially presented in [44]. Otto et al. [48, 55] proposed a coupled FE-NURBS discretization approach where an auxiliary NURBS layer is placed between the contact zone of higher-order 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. [49] presented a higher-order mortar-based contact element, where the contact interface and bulk are discretized using the $hp-$FEM. In contrast to the enrichment approach of Corbett and Sauer [44, 45], they used higher-order Lagrange functions as a basis for the discretization of the contact interface, domain, and approximation of the solution field. The 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 that directly allows the controllable order elevation of NURBS discretized structures and accordingly complements the performance of IGA while fully retaining its intrinsic key features has been devoted.

In the present work, we thus propose a novel varying-order (VO) based NURBS discretization method to improve the performance of NURBS-based IGA technique in the framework of computational contact mechanics. In the proposed method, the user-defined higher-order NURBS functions are employed for the contact computations. The minimum-order of NURBS capable of representing the complex geometries exactly are used for the description of overall bulk domain that does not come into contact. To achieve the VO based NURBS discretization, a new higher-order NURBS layer is used as the contact boundary layer of an initially NURBS discretized geometry. This avoids the application of higher-order 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 across the contact interface and matches with the bulk parameterization. The proposed varying-order based NURBS discretization method is accordingly denoted by N${}_{p}-$N${}_{p_{c}}$, where Np is the order of the NURBS functions utilized for the description of the bulk domain and N${}_{p_{c}},(p_{c}>p)$ is the order of NURBS used for the contact boundary layer. Due to its simplicity, the proposed approach can be conveniently embedded into an existing IGA contact code after only a few modifications.

In the this work, the proposed methodology is applied to two-dimensional frictionless and frictional contact problems, considering both small and large deformations. A classical full-pass version of the contact algorithm introduced by Sauer and De Lorenzis [56, 57] is combined with the VO based NURBS discretization method in this work as a first step towards developing a simple yet computationally efficient technique for isogeometric contact analysis. Such a contact algorithm has previously introduced for Lagrange-polynomial based discretizations [58, 49], fixed-order based NURBS discretizations [18, 21, 22, 40, 3, 4] and recently for T-spline based discretizations [40, 3, 4]. The penalty method is employed for the regularization of impenetrability and sticking contact constraints.

The remainder of the paper is structured as follows. In Section 2, the existing geometric modeling procedure using NURBS and the strategies used for the refinement of a NURBS discretized geometry are briefly discussed. Section 3 describes the mathematical formulation for a large deformation contact between two deformable bodies considering friction. The varying-order based NURBS discretization technique and its implementation into the existing IGA contact code are fully described in Section 4. In Section 5, the performance of the proposed method is demonstrated w.r.t. fixed-order based NURBS discretizations using various numerical examples. Finally, Section 6 concludes the paper with future directions.

2 Preliminaries

In this section, we briefly discuss the existing NURBS based discretization procedure used for the continuum and the contact surface in the context of IGA. Next, an overview of the different strategies used for the refinement of a NURBS geometry is presented. The associated issues are also highlighted. For the detailed description we refer to the monographs by Pigel and Tiller [59] and Cottrell et al. [11].

2.1 Geometric modeling using NURBS

The IGA has emerged as a successful integration of computer-aided-design (CAD) and finite element analysis (FEA) technique [11]. It utilizes the CAD polynomials as a basis for the modeling of complex shaped geometries exactly and approximation of unknown solution fields. In current work, NURBS functions are used to describe the continuum as well as contact surface of a geometry.

It is known that NURBS are built from B-splines [59]. Thus, we begin with a brief introduction to B-splines. For a given knot vector $\Xi^{i}$ along the $\xi^{i}~{}(i=1,2)$ parametric direction, a B-spline function of interpolation order333In this work, the linear, quadratic, cubic, etc. piecewise polynomials are denoted with the functions of order $p=1,2,3,etc.,$ respectively. This is as per the notation introduced by Hughes et al. [1]. The term “order” is usually referred to as “degree” in the literature based on computational geometry. $p_{i}$ is defined using the following Cox-de Boor recursive relation [59]:

[TABLE]

where $N_{l,p_{i}}(\xi^{i})\geq 0$ and $\xi_{l}^{i}\in\mathbb{R}$ is the $l^{\textrm{\scriptsize{th}}}$ knot of the knot vector $\Xi^{i}=\{\xi_{1}^{i},\xi_{2}^{i},\dots,{\xi_{n_{i}+p_{i}+1}^{i}}\}$. A knot vector is a set of non-decreasing values of parametric coordinates, i.e. $\xi_{l}^{i}\leq\xi_{l+1}^{i}$. Here, $n_{i}$ denotes the total number of control points for $p_{i}$ order of B-spline defined along the $\xi^{i}$ parametric direction. A $p^{\textrm{\scriptsize{th}}}_{i}$ order of B-spline function offers $C^{p_{i}-m}$ continuity across each knot $\xi^{i}$, where $m$ is the knot multiplicity. In an open knot vector, which are commonly adopted in IGA, the first and last knot entries are repeated by $p_{i}+1$ times that make the NURBS functions interpolatory at the end knots.

NURBS are the projective transformation of the B-splines [59]. Thus, NURBS additionally utilize the weight values in their construction as compared to B-splines. For a given knot vector $\Xi^{i}$, a $p^{\textrm{\scriptsize{th}}}_{i}$ order of univariate NURBS polynomial along the $\xi^{i}$ parametric direction is defined as [59]

[TABLE]

where $R_{l}^{p_{i}}(\xi^{i})\geq 0~{}~{}\forall\,\xi^{i}\in\Xi^{i}\,,$ and $w_{l}^{i}>0$ is the weight value. For a specified control points vector $\mathbf{X}^{i}=\{\mathbf{X}_{l}^{i}\}_{l=1}^{n_{i}}\,,$ a $p_{i}^{\scriptsize{\textrm{th}}}$ order of NURBS curve defined along the $\xi^{i}$ parametric direction can be constructed using the linear combination of univariate NURBS functions and the control points coordinates as

[TABLE]

A bivariate continuum patch is constructed by the tensor product of the two univariate NURBS curves defined in the $\xi^{1}$ and $\xi^{2}$ parametric directions as

[TABLE]

where $\{\mathbf{X}_{lm}\}_{l,m=1}^{n_{1},n_{2}}$ is the coordinate vector of the control points net $\mathbf{X}$. Moreover, $R_{lm}(\xi^{1},\xi^{2})>0$ is a bivariate NURBS function that is given by the tensor product of two univariate B-spline functions defined along the $\xi^{1}$ and $\xi^{2}$ parametric directions as

[TABLE]

where $W(\xi^{1},\xi^{2})={\sum_{l=1}^{n_{1}}\sum_{m=1}^{n_{2}}w_{lm}N_{l,p_{1}}(\xi^{1})N_{m,p_{2}}(\xi^{2})}$ is a normalized weight function that is defined in terms of weight point $w_{lm}$ and the B-spline functions. A composition of knot vector with associated control points accompanied by weights constitutes a patch. In the existing NURBS based discretization procedure, the parametrization for the contact surface is directly inherited from the bulk parametrization, but only in any one of the two parametric directions.

2.2 Refinement strategies

In the context of IGA, a NURBS described geometry can be refined by means of knot insertion ($h-$refinement), order-elevation ($p-$refinement), and $k-$refinement based strategies [1, 11]. In the knot insertion strategy, an additional knot in the knot vector is inserted. If the inserted knot value is unique, an additional knot-span, consequently an element, without changing the inter-element continuity of the NURBS, is introduced. On the other hand, the repetition of knots reduces the smoothness of the NURBS functions across the element boundary.

In the order elevation based strategy, the geometry is refined by means of raising the order of the NURBS interpolations. For the elevation of order from $p_{i}$ to $p_{i}+t$, each knot value is repeated by $t$ times. As a result, the continuity of the NURBS remains unchanged during the order-elevation ($p-$refinement).

The particular order of application of order elevation and knot insertion strategies yields $k-$refinement [1]. During $k-$refinement, first, the inter-element 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 strategy, the inter-element continuity of the NURBS remains unchanged and knot values are repeated due to order-elevation. This, as a result, yields a large number of control points as compared to $k-$refinement based strategy.

Within the existing NURBS-based discretization approach, the interpolation order of the NURBS constructed geometry is uniformly elevated to refine its contact boundary layer. However, this approach is not computationally favorable since the higher-order NURBS used for the evaluation of the contact integrals additionally have to be employed for the computation of the large majority of the region that is away from the contact interface. It, therefore, becomes desirable to develop an improved NURBS-based discretization method that provides a possibility to perform controllable order-elevation for a NURBS described geometry, as also suggested by Temizer et al. [19]. This method is introduced in this paper and is described in detail in Section 4.

3 Frictional contact formulation

In this section, we present the mathematical formulation for two-dimensional large deformation contact between two deformable bodies, considering Coulomb friction. The presented formulation is based on [60, 61, 56, 57].

3.1 Problem description and weak form

It is assumed that two elastic bodies come into contact and undergo large deformations, see Fig. 1 for schematic illustration. According to the convention in [61], one of them is denoted as the slave body $\mathcal{B}^{\textrm{s}}$ and another as master body $\mathcal{B}^{\textrm{m}}$ occupying a bounded domain $\Omega^{k}$ in $\mathbb{R}^{2}$. In the following, the superscript $k=\{\textrm{s},\textrm{m}\}$ is used to specify the slave and master bodies, respectively. The boundary ${\Gamma}^{k}$ of a body $\mathcal{B}^{k}$ can consists following three distinct parts: $\Gamma^{k}_{\sigma}$ where boundary tractions are prescribed, $\Gamma^{k}_{u}$ where displacements are prescribed, and $\Gamma^{k}_{\textrm{c}}$ as the interface where contact will occur. It is assumed that these boundaries satisfy: $\Gamma^{k}_{\sigma}\,\cap\,\Gamma^{k}_{u}\,\cap\Gamma^{k}_{\textrm{c}}=\emptyset$.

The deformed (or current) configuration of each body $\mathcal{B}^{k}$ is given by in terms of their generic point coordinates $\bm{x}^{k}$ that are expressed as $\bm{x}^{k}=\bm{X}^{k}+\bm{u}^{k}\,,$ where $\bm{X}^{k}$ denotes the coordinates of a generic points $\bm{x}^{k}$ in the initial (or reference) configuration and $\bm{u}^{k}$ is the displacement field of the same points of body $\mathcal{B}^{k}$. The master contact surface $\Gamma_{\mathrm{c}}^{\mathrm{m}}$ is parametrized using the convective coordinate ${\xi}^{\textrm{m}}$ and the covariant tangent vector ${{\bm{\uptau}}}_{1}$ at ${\xi}^{\textrm{m}}$ is defined as

[TABLE]

The contravariant vector follows from ${\bm{{{\uptau}}}}^{1}:=m^{11}\bm{\uptau}_{1}$, where the inverse metric is given by $m^{11}:=1/(\bm{\uptau}_{1}\cdot\bm{\uptau}_{1})$.

We write the variational form for a system involving contact using the principal of the stationary total potential energy as [61]

[TABLE]

where $\Pi_{\textrm{int}}$, $\Pi_{\textrm{c}}$ and $\Pi_{\textrm{ext}}$ are the total internal potential, contact potential, and potential stemming from the external loads, respectively. The function $\mathcal{V}$ denotes the space of kinematically admissible variation of the displacement field $\delta\bm{u}$. For the examples analyzed in this work $\Pi_{\textrm{ext}}=0$ holds. According to [61, 57], the total internal potential and contact potential are defined as

[TABLE]

where $\bm{\sigma}^{k}$ denotes the Cauchy stress tensor at a point $\bm{x}^{k}\in\Omega^{k}$ and $\bm{t}^{s}$ is the contact traction at a point on the contact interface of the slave body.

In the following sections, we present the contact search procedure for the determination of contact surface and the kinematical relations for the normal and frictional contact that are necessary for the evaluation of contact potential $\delta\Pi_{\textrm{c}}$.

3.2 Normal contact

It is known that the contact surface in which the one body comes into contact to another is not known a priori. Thus, for the determination of the contact surface in the current configuration, a distance function $\hat{d}({\xi}^{\textrm{m}}):=||\bm{x}^{\textrm{s}}-{\bm{x}}^{\textrm{m}}(\xi^{\textrm{m}})||$ is introduced. Here, ${\bm{x}}^{\textrm{m}}(\xi^{\textrm{m}})$ denotes the physical coordinates of an arbitrary point computed at the parametric point $\xi^{\mathrm{m}}$ on $\Gamma_{\textrm{c}}^{\mathrm{m}}$ in the current configuration. With this definition, the distance between a given (fixed) slave point $\bm{x}^{\textrm{s}}\in\Gamma^{\textrm{s}}_{\textrm{c}}$ on the slave contact surface and an arbitrary point ${\bm{x}}^{\textrm{m}}(\xi^{\textrm{m}})\in\Gamma^{\textrm{m}}_{\text{c}}$ on the master surface is calculated. For the determination of closest contact point $\bar{x}^{\mathrm{m}}$ onto the master surface corresponding to a slave point, the minimum of distance function which can be expressed as [61]

[TABLE]

The above derivative vanishes if the line describing the distance between the slave point $\bm{x}^{\textrm{s}}$ and master point ${\bm{x}}^{\textrm{m}}(\xi^{\textrm{m}})$ is orthogonal to the direction of the covariant tangent vector $\bm{\uptau}_{1}$ computed at ${\xi}^{\textrm{m}}$, see Eq. (9). The corresponding parametric point $\xi^{\textrm{m}}$ at which $\frac{d}{d\xi^{\textrm{m}}}\,\hat{d}({\bar{\xi}}^{\textrm{m}})=0$ is referred as closest projection point $\bar{\xi}^{\textrm{m}}$, which is denoted by a bar over $\xi^{\mathrm{m}}$ in the subsequent description. Since the Eq. (13) is nonlinear, an iterative solution method is required to compute the closest point $\bar{\xi}^{\textrm{m}}$. For this, the local Newton-Raphson method, starting from an initial guess ${\xi}^{\textrm{m}}=\tilde{\xi}^{\textrm{m}}$ is used. It is computed by taking the second derivative of the distance function is given by

[TABLE]

The physical coordinates of the closest point $\bar{\xi}^{\textrm{m}}$ on $\Gamma^{\textrm{m}}_{\text{c}}$ corresponding to the slave point $\bm{x}^{\textrm{s}}$ is denoted by $\bar{\bm{x}}^{\textrm{m}}={\bm{x}}^{\textrm{m}}(\bar{\xi}^{\textrm{m}})$. Further, the covariant tangent vector computed at projection point $\bar{\xi}^{\textrm{m}}$ using Eq. (9) is denoted by $\bar{\bm{\uptau}}_{1}$.

Once the physical projection point $\bar{\bm{x}}^{\textrm{m}}$ corresponding to a slave point $\bar{\bm{x}}^{\textrm{s}}$ is evaluated, the normal gap $\textrm{g}_{\textrm{N}}$ between these points is determined by

[TABLE]

where $\bar{\bm{n}}$ is the outward unit normal vector at $\bar{\bm{x}}^{\textrm{m}}$, see Fig. 1. It is used to distinguish whether the contact between the two bodies is active or non-active. In case of penetration $(\textrm{g}_{\textrm{N}}<0)$, the normal contact tractions $\text{t}_{\text{N}}=\bm{t}\cdot\bar{\mathbf{n}}$ are activated that avoid the penetration of the contacting regions. In order to enforce the impenetrability in the normal direction, following Karush-Kuhn-Tucker (KKT) conditions are considered

[TABLE]

The above-described contact constraints cannot be directly incorporated into the contact variational form as they lead to the non-smooth relationship between the normal gap and the contact pressure, resulting in a non-smooth normal contact constitutive law. In the present work, for the regularization of the contact constraints defined in Eq. (16) the penalty method is adopted. According to [61, 57], the penalty regularized normal contact constraint is defined as

[TABLE]

where $\epsilon_{\text{N}}>0$ denotes the normal penalty parameter.

3.3 Frictional contact

For frictional contact, the contact traction vector $\bm{t}$ is described in terms of the normal and tangential components as [57]

[TABLE]

where $\textrm{t}_{\textrm{N}}$ is defined in Eq. (17), and $\textrm{t}_{\text{T}}=\bm{t}\cdot\bar{\bm{\uptau}}_{1}$ is the tangential traction. Also, $\bar{\bm{n}}$ and $\bar{\bm{\uptau}}_{1}$ are the unit normal and tangent vectors computed at $\bar{\xi}^{\textrm{m}}$. The tangential traction $\bm{t}_{\textrm{T}}=\textrm{t}_{\textrm{T}}\bar{\bm{\uptau}}_{1}$ is determined based on the stick-slip behavior that occurs during frictional contact. The distinction between the stick and slip states is drawn by the slip function $\Phi=||\bm{t}_{\text{T}}||-\mu_{f}\,\text{t}_{\text{N}}\,,$ where $\mu_{f}$ denotes the friction coefficient. It states sticking for $\Phi<0$ and sliding if $\Phi=0$. The function $\Phi=0$ corresponds to a surface in the traction space $\{\text{t}_{\text{N}},\text{t}_{\text{T}}\}$. The tangential traction during the frictional sliding is defined, e.g. for Coulomb’s law, as

[TABLE]

where $\mu_{f}$ denotes the Coulomb friction coefficient and $\dot{\textrm{g}}_{\textrm{T}}$ is the relative tangential sliding velocity.

In the context of computational contact mechanics, the constitutive equations for friction are often formulated in the framework of elasto-plasticity [61]. According to this, the total tangential slip $\mathbf{g}_{\textrm{T}}$ can be decomposed into an elastic (or stick) part $\mathbf{g}_{\textrm{T}e}$ and a plastic (or slip) part $\mathbf{g}_{\textrm{T}sl}$ as

[TABLE]

The tangential stick $\mathbf{g}_{\textrm{T}e}$ is characterized by no relative tangential motion between the surfaces during the frictional contact. Such a physical behavior enforces a geometrical constraint in the relative motion of the contact surface. For its regularization the penalty method is used herein. The tangential slip $\mathbf{g}_{\textrm{T}sl}$ is characterized by the sliding of a slave point $\bm{x}^{\textrm{s}}$ relative to the master surface $\Gamma_{\textrm{c}}^{\textrm{m}}$. It is a dissipative process. The constitutive evolution equation for tangential slip $\mathbf{g}_{\textrm{T}sl}$ can be derived by considering the dissipation during friction. We write the evolution law for the tangential slip as [61]

[TABLE]

where $\gamma$ is the slip parameter that is determined using the following conditions

[TABLE]

The above set of equations is known as KKT conditions for Coulomb friction [61].

The tangential traction is updated based on a classical corrector-predictor approach. In this, first, the (penalty regularized) elastic trial step is computed, followed by the trial slip function, as in [57], with

[TABLE]

Here, all the quantities defined in the current configuration at time step $t=t_{n+1}$ are denoted with subscript $n+1$. Similarly, all the quantities at previous time step $t=t_{n}$ with subscript $n$. The tangential penalty parameter is denoted by $\epsilon_{\textrm{T}}>0$. Moreover, $\bm{x}^{\textrm{m}}_{n+1}(\bar{\xi}^{\textrm{m}}_{n+1})$ and $\bm{x}^{\textrm{m}}_{n+1}({\xi}^{\textrm{m}}_{{sl}\,n})$ are the physical coordinate corresponding to the closest projection point $\bar{\xi}^{\textrm{m}}_{n+1}$ and tangential slip point ${\xi}^{\textrm{m}}_{{sl}\,n}$, respectively, see Fig. 1. According to [61, 57], the tangential traction is determined by based on the changes in stick-slip status with

[TABLE]

where $\bm{n}_{\textrm{T}n+1}={\mathbf{t}_{\text{T}n+1}^{\textrm{trial}}}/{||\mathbf{t}_{\text{T}n+1}^{\textrm{trial}}||}$.

3.4 NURBS discretized weak form

Having discretized the body into a set of elements $\Omega^{k}=\cup_{e}\Omega^{ke}$, the same NURBS functions, as employed in Eq. (7), are used for the approximation of the unknown displacement field, its variation and current coordinates within each element as

[TABLE]

Here, $\mathbf{R}^{k}=[R_{1}^{k}\bm{I},R_{2}^{k}\bm{I},\,\,.\,.\,.\,,R_{n_{cp}^{ke}}^{k}\bm{I}]$ is an array of size $2\times 2n_{cp}^{ke}$ with identity matrix $\bm{I}$ in ${\rm I\!R}^{2}$. The $n_{cp}^{ke}=(p_{1}^{k}+1)\times(p_{2}^{k}+1)$ denotes the total number of control points whose bivariate NURBS basis functions $(R^{k}_{a}(\bm{\xi}^{k}))$ have local support in an element $\Omega^{ke}$. Also, $\mathbf{u}^{ke}$, $\delta\mathbf{u}^{ke}$ and $\bm{x}^{ke}$ denotes the arrays containing the displacement field, its variation and the current coordinates of $n_{cp}^{ke}$ number of control points in an element, respectively. The modification in the approximation of $\bm{u}^{ke}$, $\delta\bm{u}^{ke}$ and $\bm{x}^{ke}$ for VO based NURBS discretization is explained in Sec. 4.

Using the definition of the displacement field given in Eq. (26), the NURBS discretized, Eq. (10), takes the form

[TABLE]

where $\mathbf{f}_{\textrm{int}}$ and $\mathbf{f}_{\textrm{c}}$ are the global control point force vectors corresponding to the $\Pi_{\textrm{int}}$, defined in Eq. (11), and $\Pi_{\textrm{c}}$, in Eq. (12), respectively.

Since, the variation of displacement $\delta\mathbf{u}$ is arbitrary, the Eq. (27) can be written as $\mathbf{f}_{\textrm{int}}+\mathbf{f}_{\textrm{c}}=\mathbf{0}$. The global forces $\mathbf{f}_{\textrm{int}}$ and $\mathbf{f}_{\textrm{c}}$ are obtained based on the assembly of their elemental contributions $\mathbf{f}_{\textrm{int}}^{ke}$ and $\mathbf{f}_{\textrm{c}}^{ke}$ using the control point connectivity arrays for each body $\mathcal{B}^{k}$, see [11, 62].

The elemental internal force vector for a volume element $\Omega^{ke}$ is given by

[TABLE]

where $\mathbf{B}^{k}$ array contains the $n_{cp}^{ke}$ number of the derivatives of bivariate NURBS functions having local support in an element $\Omega^{ke}$. In the examples solved herein, the Neo-Hookean material model is used for the description of bodies. The corresponding constitutive relation is given by [63]

[TABLE]

where $\bm{F}$ is the deformation gradient, $J$ is the determinant of $\bm{F}$, and $\bm{I}$ is the identity tensor. The Lamé constants, i.e. the bulk modulus and shear modulus, are defined via $\lambda=2\mu\nu/(1-2\nu)$ and $\mu=E/2(1+\nu)$, respectively.

3.5 Contact algorithm

In this work, we adopt the classical full-pass version of the contact algorithm introduced by Sauer and De Lorenzis [57] to integrating the proposed VO based NURBS discretization methodology with the isogeometric contact analysis. The full-pass contact algorithm was first introduced by Fischer and Wriggers [58] for the Lagrange-polynomial based FE analysis of the contact problem. It was later extended by Temizer et al. [18], De Lorenzis et al. [21, 22] and Dimitri et al. [40, 3, 4] for the NURBS based isogeometric contact analysis. Recently, Dimitri et al. [40, 3, 4] further extended this algorithm for T-spline based isometric analysis of contact, cohesive zone modeling, and mixed-mode debonding problems. In these literature, a number of terminologies have been adopted to denote the classical full-pass approach such as “mortar” in [58], “Knot-to-Surface” in [18], “non-mortar” in [21, 22], and “Gauss-point-to-surface (GPTS)” in [40, 3, 4]. In the present work, as in [40, 3, 4], we also denote this approach as GPTS, since within this the contact constraints, i.e. impenetrability in Eq. (16) and sticking in Eq. (22), are enforced independently at each quadrature point that enables the evaluation of the contact integrals. This formulation passes the contact patch test within the quadrature error, which ensures the convergence of the solution upon mesh refinement [58, 64]. The only drawback is that the GPTS based contact algorithm leads to over-constrained formulation, as demonstrated by [18, 21, 22]. The over-constraining stems from directly enforcing the contact constraints at too many numbers of locations that, as a result, leads to numerical instability, see [65] for the contact cases. Further, with GPTS approach, oscillatory responses of contact forces are obtained if a large value of penalty parameter is employed. However, in our examples, we use the value of penalty parameter, for which the loss of the numerical stability issues becomes significant with the GPTS approach, lay beyond those required to obtain a satisfactory quality solution from the engineering perspective. The values lies within a range where the ill-conditioning of the global system matrix related issues still arises. The post-processing smoothening scheme of Sauer [51] can also be used as in [40, 3] to alleviate the oscillatory responses of contact forces effectively. The incorporation of the mathematically more sophisticated mortar-based isogeometric contact algorithm by Temizer et al. [18, 19], De Lorenzis et al. [21, 22], Seitz et al. [31], Popp et al. [66, 67, 68] and Duong et al. [32] is subject of future research work.

The elemental contact force vector $\mathbf{f}_{c}^{ke}$ on the contact surfaces $\Gamma_{c}^{ke}$ of a contact element $\Omega^{ke}_{c}$ for a body $\mathcal{B}^{k}$ is given as

[TABLE]

where $\mathbf{N}^{k}$ array consists $n_{cps}^{ke}=p_{1}^{k}+1$ number of univariate NURBS basis functions having support on the contact surface $\Gamma_{c}^{ke}$. In this, all the quantities are computed within the current configuration of the bodies. The contact traction $\bm{t}^{\textrm{s}}$ is computed at each quadrature point using Eqs. (17) and (25) for frictionless and frictional contact, respectively. An active set strategy that incorporates the contribution of the quadrature points in case of penetration is utilized for the evaluation of contact integral. The changes in the evaluation of $\mathbf{f}_{\textrm{c}}^{ke}$ with VO based NURBS discretization is detailed in the next section.

The linearization of the internal force vector $\mathbf{f}_{\textrm{int}}^{ke}$ and contact force vector $\mathbf{f}_{\textrm{c}}^{ke}$ is necessary for the solution of resulting nonlinear system of equations using the Newton-Raphson method. The linearization of the contact force vector yields the contact tangent matrix, see [56, 57], and their corresponding expressions are provided in Appendix A for the sake of completeness.

4 Varying-order based NURBS discretization

In this section, we introduce the theory for the varying-order based NURBS discretization method, followed by the description on its integration into an existing isogeometric contact code.

4.1 Varying-order NURBS

The basic concept of the proposed varying-order (VO) based NURBS discretization for isogeometric contact analysis is illustrated in Fig. 2. Consider a body $\mathcal{B}^{k}$ having contact boundary layer $\Gamma_{\textrm{c}}^{k}$. Let $p_{1}$ and $p_{2}$ be the minimum orders of NURBS functions that are capable of representing the given geometry in an exact manner. The coarse mesh for the geometry is given by the product of open knot vectors $\Xi^{1}~{}\times~{}\Xi^{2}$ defined along the $\xi^{1}$ and $\xi^{2}$ parametric directions, as shown in Fig. 2a. Next, in order to make use of higher-order NURBS functions for contact computations, the originally $p_{1}$ discretized NURBS contact boundary $\Gamma_{\textrm{c}}^{k}$ is replaced with a higher-order $p_{c}>p_{1}$ of NURBS layer as shown in Fig. 2b. The resultant discretization is accordingly denoted by N${}_{p}-$N${}_{p_{c}}$, where Np $(p=~{}\mathrm{max}(p_{1},p_{2}))$ is the order of NURBS used for the description of bulk domain and N${}_{p_{c}}$ $(p_{c}>p)$ for the contact boundary layer.

The new NURBS layer is constructed either using the $k-$refinement or through a combination of $k-$refinement and order-elevation strategies in such a manner that it matches the bulk parametrization, as shown in Fig. 2b. The application of the $k-$refinement strategy to the NURBS contact layer increases the order as well as the inter-element continuity of the NURBS functions, see Fig. 3. On the other hand, application of one or two additional steps of order-elevation to $k-$refined NURBS layer introduces a large number of additional control points along the contact boundary layer while the inter-element continuity, i.e. $C^{p_{c}-1}$, remains unchanged, see Fig. 4. With this approach, the resultant VO based NURBS discretization is denoted by N${}_{p}-$N${}_{p_{c}\cdot p_{s}}$, where $p_{s}$ is the step number of order elevation strategy that is additionally applied to N${}_{p_{c}}$ NURBS layer.

In Figs. 3 and 4, one VO NURBS discretized contact element is highlighted to illustrate the difference between the above two strategies. As shown in Fig. 3, for employing the higher-continuous NURBS for contact computations, the contact boundary layer of an original N2 discretized example geometry is replaced with a user-defined $C^{p_{c}-1}$ continuous N${}_{p_{c}}$ NURBS layer, where $(p_{c}=3,\mathrm{and}~{}4)$. For a N${}_{p_{c}}$ discretization, $(p_{c}+1)$ number of basis functions have local support in a contact element, see right column in Fig. 3. However, as shown in Sec. 5, very smooth basis functions are incapable of accurately capture the sharp changes in the distribution of contact pressure result. This is related to the multi knot-span support of the smooth functions, which widens by an additional knot-span on linearly increasing their smoothness, e.g. from $C^{2}$ to $C^{3}$ as shown in Fig. 3. Thus, the much higher-continuous NURBS does not considerably improve the accuracy of the contact responses results.

On the contrary, as the basis functions plots show in Fig. 4, the application of additional one step of order elevation to the sufficiently $C^{1}$ smooth contact layer leads to the shrinking of the multi-knot span support of the basis functions. Further, the knot-span support of the basis functions remains unvaried on applying an additional or second step of order elevation to the contact layer. Such a sufficiently smooth and higher-order NURBS functions with a moderate knot-span support can adequately capture the local changes in the contact pressure results as shown in Sec. 5. It is also highlighted that although each type of N${}_{2}-$N${}_{p_{c}}$ $(p_{c}=3,~{}\mathrm{and}~{}4)$ discretization, and its equivalent version N${}_{2}-$N${}_{p_{c}\cdot p_{s}}$, $(p_{c}\cdot p_{s}=2\cdot 1,~{}\mathrm{and}~{}2\cdot 2)$, respectively, has same $(p_{c}+1)$ or $(p_{c}+p_{s}+1)$ number of non-zero basis functions across the contact boundary of an element. But, the moderate knot-span support of a N${}_{2}-$N${}_{p_{c}\cdot p_{s}}$ is accompanied by additional large number of overall control points across the contact interface as compared to N${}_{2}-$N${}_{p_{c}}$ discretization, respectively, see Figs. 3 and 4.

This way with the VO based NURBS discretization of a contact geometry: $(i)$ the higher-order NURBS basis functions are utilized only for the evaluation of contact contributions in a fully NURBS discretized geometry, $(ii)$ a large number of additional degrees of freedom are introduced across the contact surface without changing the mesh size, and $(iii)$ the minimum order of NURBS are employed for the description of the bulk domain that does not come into contact. The resulting VO based NURBS discretized contact element is characterized by $n_{cp}^{e}=(p_{c}+1)+(p_{1}+1)\times p_{2}$ number of control points, where $p_{c}+1$ are present on the contact layer and $(p_{1}+1)\times p_{2}$ in its remaining part. The bivariate NURBS basis functions for such an element are defined as

[TABLE]

where the normalizing weight function is given by

[TABLE]

The basis functions defined in Eq. (31) exhibit the non-negativity property

[TABLE]

and they also satisfy the partition of unity property:

[TABLE]

Using the isoparametric concept, the polynomials used for the VO based NURBS discretization, as defined in Eq. (31), are employed for the approximation of the unknown displacement field $\mathbf{u}^{e}$, its variation $\delta\mathbf{u}^{e}$, and the current coordinates $\bm{x}^{e}$ within each contact element as

[TABLE]

Here, the new basis function matrix $\mathbf{R}(\bm{\xi})$ contains a total $n_{cp}^{e}=(p_{c}+1)+(p_{1}+1)\times p_{2}$ number of bivariate NURBS functions $R^{p_{c},p_{2}}_{a}(\bm{\xi})$ having local support in a VO discretized contact element $\Omega^{e}_{c}$, see Eq. (31). The array $\mathbf{B}^{ke}$ in Eq. (28) containing the derivatives of NURBS functions is also modified according to Eq. (35). Moreover, the matrix $\mathbf{N}^{ke}$ in Eq. (30) makes use of higher-order of univariate NURBS functions $R^{p_{c}}(\xi_{1})$ for the evaluation of contact integrals.

The proposed VO based NURBS discretization approach bears some similarities with the NURBS-enriched contact finite element strategy by Corbett and Sauer [44] and hybrid isogeometric-finite element based discretization technique by Maleki-Jebeli et al. [47]. The common focus of these strategies is to improve the performance of the Lagrange-polynomial based finite element method in the context of contact mechanics. For this, the intrinsic features of the NURBS polynomials are coupled with the FE based discretizations. The Bézier extraction operator [69] is required to enable the incorporation of NURBS into the FE structure. In contrast, the present work aims to improve the performance of the NURBS-based IGA technique in the context of contact mechanics while fully retaining its original key purpose, i.e. the unified treatment of design and analysis processes by employing the same functions for the analysis that are used for the construction of a geometry exactly. For this purpose, the VO based NURBS discretization method is introduced. This unexplored idea, as initially noted by Temizer et al. [19] presents a possibility to perform the controllable order-elevation based refinements of the NURBS discretized structures. It is noted that the proposed discretization method can also be converted to the enrichment strategy of Corbett and Sauer [44] if the linear order of NURBS444As they are equivalent to linear order of Lagrange polynomials [11] are chosen for the bulk description. The proposed method, however, unlike the enrichment strategy of [44], does not require the Bézier extraction operator.

4.2 Implementation into existing code

For integrating the VO NURBS discretization strategy into the existing isogeometric contact code, only a few minor modifications are required. First of all, for a given mesh resolution, a $p_{c}>p_{1}$ order of NURBS curve representing the contact boundary layer of the initial NURBS described geometry is constructed. After that, the parametrization for an originally $p_{1}$ NURBS discretized contact layer is replaced with that of the newly constructed $p_{c}$ order of NURBS curve. For this, a certain number of conditions are need to be fulfilled. The total number of control points defining the geometry must be updated in such a manner that it allows the incorporation of the newly constructed $p_{c}$ order of contact layer. This means that the connectivity array for contact elements must be adapted in a way that it contains the underlying control points of the VO based discretized geometry. The derived contact element connectivity arrays can have a different length than the bulk element connectivity arrays. The bivariate NURBS basis functions defined in Eq. (31) are used for the evaluation of elemental quantities for VO based NURBS described contact elements. The univariate $p_{c}$ order of NURBS functions are utilized for the evaluation of contact integrals. With the exception of these modifications, no other changes are need to be made in an existing isogeometric contact code. The local quantities, e.g. elemental 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 [62] for a detailed description on the implementation of IGA in a simplified manner. In this work, a default $(p_{c}+1)\times(p_{2}+1)$ number of Gauss-Legendre quadrature points are employed for the evaluation of the contact integrals unless stated otherwise. Optimal quadrature rules [70, 71, 72], which are well-suited for IGA, can also be opted for the reduced numerical evaluation.

5 Numerical results

In this section, we show the performance of the proposed VO NURBS based discretization method using three, two-dimensional small and large deformations contact problems. In the first example, a frictional ironing problem is analyzed to demonstrate the superior performance of the proposed method in terms of accuracy, robustness, and efficiency over the standard fixed-order based discretizations. The performance is assessed through by varying the inter-element continuity as well as the interpolation order of the NURBS discretizations using the $k-$ and $p-$refinements. In the second and third examples, Hertzian contact and two elastic rings problems are simulated, respectively. In these examples, we demonstrate the coarse-mesh accuracy and cost-efficacy of the proposed method for the evaluation of contact pressure distributions in comparison to standard NURBS discretizations. As compared to the second example, the last one incorporates large deformations and friction.

5.1 Frictional ironing problem

5.1.1 Problem setup

In the first example, we consider the frictional contact between two deformable bodies using the setup similar to the one used by De Lorenzis et al. [21]. The geometric model along with the material details and the boundary conditions are shown in Fig. 5. This example is used to carefully analyze the performance of the proposed VO based NURBS discretization method over the standard fixed-order based NURBS discretization for large deformation and frictional sliding contact problem. The performance is measured in terms of the accuracy of the results, robustness, and the computational efficiency of the proposed approach.

In the considered example, the die is first pressed onto an elastic slab and then moved horizontally relatively, as shown in Fig. 5. A total vertical displacement $U_{y}=-0.23$ mm is applied to the top surface of the die in a uniform manner in $46$ load steps. Thereafter the die is dragged across the slab by applying $U_{x}=1.5$ mm in $250$ load steps while keeping the vertical displacement fixed. An isotropic Neo-Hookean hyperelastic material model is considered to describe the material behaviour of the die and slab. The material parameters are shown in Fig. 5. The penalty parameters are taken as $\epsilon_{\textrm{N}}=\epsilon_{\textrm{T}}=100$ N/mm. The Coulomb’s friction coefficient $\mu_{\textrm{f}}=0.2$ is considered.

First, a very coarse mesh is chosen to amplify the possible difference between the results with VO and standard NURBS discretizations, as discussed in the following. The considered coarsest mesh is denoted by $\textrm{m}_{1}$ and is illustrated in Fig. 6. Three other nested meshes that are obtained through the uniform knot insertion and are used for the convergence study are listed in Table 1. For a Np order of NURBS contact layer, $3\times(p+1)$ number of Gauss-quadrature points are used for the evaluation of contact integrals per contact element.

5.1.2 Prediction of vertical and horizontal contact forces

In this section, we analyze the improvement in the accuracy of the results on increasing the inter-element continuity and interpolation order of the NURBS with the proposed discretization method. The results with standard NURBS based discretizations are used for comparisons.

5.1.2.1   On increasing the inter-element continuity

First, the quality of the results with different VO based discretizations: N${}_{2}-$N${}_{p_{c}}$ ($p_{c}=4,$ and $6$) is compared with that of corresponding standard N${}_{p}~{}(p=2,4$ and $6)$ order of NURBS discretizations on increasing the inter-element continuity of the NURBS at the coarsest mesh $\textrm{m}_{1}$. In case of VO, the quadratic order of NURBS, which are sufficient to describe the setup exactly, are kept fixed for the discretization of the bulk region, while the interpolation order of the bottom surface of the die and upper surface of the slab are elevated to $4$ and $6$. The results with N3 (and N5) are omitted from this investigation, as they are similar to those obtained with N4 (and N6) cases.

The deformed configurations of the setup during downward and sliding motions are shown in Fig. 7 with N2 using mesh $\textrm{m}_{2}$. During the indentation phase, the sum of contact forces on the top surface of the die should increase gradually, and during the sliding phase, forces should be approximately constant. Figure 8a shows the evolution of net vertical contact force Py and the horizontal contact force Px, computed at the top surface of die, as a function of load step $t$ for both Np and N${}_{2}-$N${}_{p_{c}}$ based discretizations. The enlarged views for Py and Px are shown in Figs. 8b and 8c, respectively. Two major observations are made. First, the results verify the validity of the proposed VO based NURBS discretization method, as the contact force curves of $P_{y}$ and $P_{x}$ are nearly indistinguishable for VO based N${}_{2}-$N${}_{p_{c}}~{}(p_{c}=4,\textrm{ and }6)$ and standard N${}_{p}~{}(p=4,\textrm{ and }6)$ discretizations. This is due to employing the same $C^{p-1}-$continuous NURBS for the evaluation of contact integrals and identical number of DOFs present across the contact boundary in both the cases of discretizations.

The DOF density data for N${}_{2}-$N${}_{p_{c}}$ and Np is provided in Table 2 for mesh m1. Second, during sliding, large-amplitude periodic oscillations of vertical and horizontal contact forces are present, especially with standard N2 based discretization. The source of such large-amplitude non-physical oscillations is the insufficient discretization of the NURBS contact layers with the coarse mesh $\textrm{m}_{1}$. This, as a result, limits the conforming ability of the master body contact layer to the finite deformations introduced by the slave body. For a fixed mesh, conforming ability of contact layers increases on increasing the inter-element continuity with NURBS discretization. Consequently, oscillation error in Py and Px reduces as shown in Figs. 8b and 8c, respectively. With N${}_{2}-$N4 and N${}_{2}-$N6, the oscillation amplitude of the horizontal contact force $\Delta\textrm{P}_{x}$ reduce to approximately $27\%$ and $12\%$, respectively of that observed with N2 discretization. A quantitative analysis of the reduction in the oscillation error for N${}_{2}-$N${}_{p_{c}}~{}(p_{c}=4,\textrm{ and }6)$ and N${}_{p}~{}(p=2,~{}3,\textrm{ and }4)$ for both the force components is provided in Table 3. The oscillation amplitude of the reaction forces is computed using $\Delta P_{j}:=\text{max}(P_{j})-\text{min}(P_{j}),~{}\text{where }j=x,\,y$. Here, the amplitude observed with N2 is used as a reference. It is noted that the obtained results and the reduction in the oscillation amplitude on increasing the smoothness of NURBS are in the agreement to those reported in the literature [21, 22, 19, 73, 26, 42, 4]. Moreover, the convergence behaviour and the reduction in the total analysis time with N${}_{2}-$N${}_{p_{c}}$ as compared to Np discretizations are reported in Sec. 5.1.3.

5.1.2.2   On elevating the interpolation order

Within the context of standard IGA, apart from increasing the inter-element continuity of NURBS discretization, the degree of conformity of the contact layer can be increased by increasing the mesh resolution. However, with the usage of a fine mesh, the computational cost increases considerably, which is not desirable. This dismisses one of the key advantages of IGA which is its ability to represent the geometry to high-accuracy even with a coarse mesh. Therefore, the focus of the proposed VO NURBS discretization method is to improve the solution quality while keeping the mesh $\textrm{m}_{1}$ fixed. For this purpose, we carry out the additional order-elevation based refinement to N2 discretized contact layers of the die and slab at mesh m1. With this, although the inter-element continuity remains constant, a large number of additional DOFs are introduced across the contact layers at a fixed mesh attributed to the repetitions of knots with the order elevation process. Performing one and two steps of additional order-elevation to N${}_{2}-$N2 results in N${}_{2}-$N2⋅1 and N${}_{2}-$N2⋅2 discretizations, respectively. The total number of DOFs present in the contact interface and bulk domain corresponding to N${}_{2}-$N${}_{2\cdot p_{s}}~{}(p_{s}=1\textrm{ and }2)$ are listed in Table 2.

Figures 9a, and 9b show the contact force curve of Py and Px, respectively, with VO based N${}_{2}-$N${}_{2\cdot p_{s}}~{}(p_{s}=1,\textrm{ and }2)$ and standard N2 based discretizations. As can be seen, the solution quality improves significantly for the new discretizations, especially with N${}_{2}-$N2⋅2. This is due to the moderate multi knot-span of basis functions as discussed in Sec. 4.1 and a large number of DOFs present across the contact layer with N${}_{2}-$N${}_{2\cdot p_{s}}~{}(p_{s}=1,\textrm{ and }2)$ as compared to N2, see Table 2 for DOFs details. Moreover, the other higher order-elevation based discretization, i.e. N${}_{2}-$N2⋅3, also reduces the oscillation error, as shown in Figs. 9a and 9b. However, only a slight improvement is attained as compared to N${}_{2}-$N2⋅2, see Table 2. This shows that the number of DOFs more than with the N${}_{2}-$N2⋅2 across the contact layer improves the results only slightly. The quantitative details on the reduction of oscillation amplitude for each discretizations for both the force components are shown in Table 3.

Finally, we investigate the influence of order-elevation based refinement to the higher continuous NURBS contact layer, i.e. N4, again while keeping the original mesh m1 fixed. The contact force curve Px corresponding to N${}_{2}-$N${}_{4\cdot p_{s}}~{}(p_{s}=1,~{}2\textrm{ and }3)$ discretizations are shown in Fig. 10. The results with N${}_{2}-$N${}_{2\cdot p_{s}}~{}(p_{s}=1,~{}2\textrm{ and }3)$ are used for the purpose of comparison. The major observation is that on increasing the inter-element continuity of the N${}_{2\cdot p_{s}}$ discretized contact layer to N${}_{2}-$N${}_{4\cdot p_{s}}$, the oscillation error reduces marginally. This is expected as N${}_{2}-$N${}_{4\cdot p_{s}}$ has total only 12 additional DOFs across the contact interface as compared to N${}_{2}-$N${}_{2\cdot p_{s}}$, see Table 2 for DOFs. Furthermore, it can be observed that more smoother response of contact forces are obtained on increasing the continuity of NURBS discretized contact layer.

With the obtained results, it is evident that additional order-elevation based VO NURBS discretizations, particularly N${}_{2}-$N2⋅2, is advantageous over the standard N${}_{p}~{}(p=2,~{}4,\textrm{ and }6)$ as well as N${}_{2}-$Np based VO NURBS discretizations, see Table 3. The N${}_{2}-$N2⋅2 reduces the oscillation amplitude of vertical and horizontal contact forces to $4.31\%$ and $7.51\%$, respectively, as compared to N2 at a fixed mesh $m_{1}$.

5.1.3 Convergence behaviour and analysis time

Next, we investigate the convergence of the vertical contact forces oscillation amplitude $\Delta\textrm{P}_{y}$ and horizontal contact force oscillation amplitude $\Delta\textrm{P}_{x}$ with both the discretization methods upon mesh refinement. Four nested meshes, as described in Table 1, are used. In case of VO, the higher-continuous and additional order-elevation based discretizations are considered. Moreover, the total analysis time taken by VO based discretizations is compared to that with the standard NURBS discretizations for each mesh level.

Figure 11 shows the reduction in the oscillation amplitude of horizontal contact force curve Px for different discretizations. The convergence plots for both $\Delta\textrm{P}_{y}$ and $\Delta\textrm{P}_{x}$ with N${}_{2}-$N${}_{p_{c}}~{}(p_{c}=4,\,6,\,2\cdot 1,\,2\cdot 2,\textrm{ and }2\cdot 3)$ and standard N${}_{p}~{}(p=2,\,4,\,\textrm{6})$ based discretizations for all the four meshes are shown in Fig. 12. The total processing time taken by different discretizations over the DOFs are shown in Fig. 13. For each discretization, the time in $\%$ is calculated using the following expression

[TABLE]

where the maximum total analysis time is with the standard N6 order of NURBS using the finest mesh m4 (which is the most expensive).

From Figs 12a, 12c, and 13a it can be observed that to attain the accuracy equivalent to N4, N${}_{2}-$N4 takes approximately $33.74\%$ lower computational cost even at the coarsest mesh m1. The reduction in the computational cost with N${}_{2}-$N4 is due to fact that it employs the lower order of NURBS, i.e. N2, for the bulk computations as compared to N4 based uniform discretization. The gain in the computational efficiency further improves on increasing the mesh resolution, e.g. N${}_{2}-$N4 takes approximately $52.44\%$ lower cost as compared to N4 at the finest mesh m4. A similar observation is made on comparing the results with more higher-continuous NURBS discretizations, e.g. to deliver the result equivalent to N6, N${}_{2}-$N6 takes approximately $59.58\%$ lower analysis time at mesh m4.

Further, from Figs. 12 and 13 it can be observed that for a fixed mesh, a significant improvement in the accuracy is achieved with additional order-elevation based VO discretizations, i.e. with N${}_{2}-$N2⋅1 and N${}_{2}-$N2⋅2, among all tested cases. The N${}_{2}-$N2⋅2 delivers the most accurate results at a cost slightly more than with the N2. Taking a closer look at Figs. 12b and 12d reveals that the accuracy attained with N${}_{2}-$N2⋅2 at mesh m1 and m2 is comparable to N2 using mesh m3 and m4, respectively. The corresponding total analysis time taken by N${}_{2}-$N2⋅2 is approximately $56.62\%$ and $72.78\%$ lower than with the N2 for the similar accuracy level, see Fig. 13b. With other higher order based VO discretizations, i.e. with N${}_{2}-$N2⋅3, a slight improvement in the accuracy is achieved as compared to N${}_{2}-$N2⋅2 for the same mesh level.

In summary, this example demonstrates the advantages of the proposed VO based discretizations in terms of the accuracy of results and computational efficiency over the standard fixed-order based NURBS discretizations. To attain the accuracy equivalent to Np based discretization, N${}_{2}-$Np takes much lower computational cost. Further, with the additional order-elevation based VO discretizations, a significant improvement in the accuracy is achieved even at the coarse mesh. It attains a major gain in the numerical efficiency to provide the accuracy as similar to N2 based discretization.

5.1.4 Performance at large indentation depth

Finally, in this section, we examine the influence of the interpolation order of the underlying bulk elements on the solution quality for the ironing problem with large indentation depth. It is known that numerical simulation of this problem can be computationally challenging if the stiffness ratio of the bodies is high and the harder body is dragged relative to the softer one after the significant indentation [74, 73]. In such a case, the strong interaction between the bulk and contact surface discretizations takes place. Therefore, it becomes important to assess the accuracy of the solution with VO discretization to that of corresponding standard NURBS discretization for large deformation contact. The problem setup used for this analysis is same as described in Sec 5.1.1 except the modification that higher stiffness for the die $E_{\textrm{die}}=10^{4}$ N/mm2 and vertical displacement $U_{y}=-0.5$ mm are considered. The Coulomb’s friction coefficient is $\mu_{\textrm{f}}=0.1$. Four nested meshes that are listed in Table 1 are used for the analysis. The deformed configurations of the setup at the end of compression and sliding processes are shown in Fig. 14 with N2 using mesh m2.

Figure 15 shows the comparison of the oscillation amplitude of the vertical contact force $\Delta\textrm{P}_{y}$ and horizontal contact force $\Delta$Px with both the discretization methods for meshes m1 to m4. Two major observations are made. First, although N2 based discretization is $C^{1}-$continuous across the contact layer, it fails to converge at the coarsest mesh m1. Clearly, N2 with the mesh m1 is not sufficient to analyze this problem. On the other hand, the higher-order based NURBS discretization is capable of providing the solution even at the coarsest mesh m1. Second, the results obtained with the VO based discretization verify the application of the proposed methodology to large deformation contact case. For the same mesh level, the accuracy of the result achieved with VO based N${}_{2}-$N${}_{p_{c}}~{}(p_{c}=4\textrm{ and }6)$ discretizations is equivalent to that of fixed-order based N${}_{p}~{}(p=4\textrm{ and }6)$ discretizations, respectively.

5.2 Hertzian contact problem

In the context of isogeometric contact analysis, Hertzian-type contact problem has been widely used to asses the contact pressure distribution across the contact interface on increasing the mesh resolution and interpolation order of the NURBS with the different type of contact formulations, see e.g. [18, 21, 19, 22, 27, 25, 26, 40, 31]. Therefore, as a second example, we consider it to demonstrate the coarse mesh accuracy of the VO NURBS discretization method as compared to standard NURBS discretizations.

The example considers the contact between a deformable infinitely long cylinder, having an outer radius $R_{o}=1$, and a rigid surface under the plane strain setting [75]. The setup of the problem along with the boundary condition, which is taken from [18], is shown in Fig. 16. Due to symmetry, only a quarter of the geometry is utilized. The natural NURBS based description of geometry consists a non-zero internal radius $R_{i}=0.1$. The numerical investigation reveals that it exerts no influence on the quality of the solution. For its modeling, a linearly elastic material with Young’s modulus $E=1$ and Poisson’s ration $\nu=0.3$ is used. Here, the top surface of the cylinder is subjected to prescribed vertical displacement $\bar{v}$ and a penalty parameter $\epsilon_{\text{N}}=2\times 10^{3}$ is chosen as a default value. Six uniformly refined meshes that are driven by $9n\times 48$ number of elements along each parametric directions $(\textrm{where }n=1,\,2,\,4,\,8,\,16\textrm{ and }32)$ and are denoted by m1 to m6 are used for this analysis. For numerical efficiency, approximately $80\%$ of the total elements in each parametric direction are relocated in such a manner that they lie within the $10\%$ surface length of the geometry, as shown in Fig. 17. It has been tested that $48$ number of redistributed elements along the radial direction is adequate for the considered example. The undeformed and deformed configurations of the setup with N2 using mesh m3 are shown in Figs. 17a and 17b, respectively.

Figure18 shows the comparison between the local solution quality of contact pressure distribution for both the standard N2, and VO based N2⋅1 and N2⋅2 discretizations with that of reference contact pressure distribution for meshes m3 to m5. Here, the dimensionless normal contact pressure $\bar{p}_{N}=p_{N}/p_{N_{ref}}$ is plotted versus the dimensionless contact coordinate $\bar{X}=x/a_{ref}$, where $p_{N}$ is the normal contact pressure evaluated at an active integration point and $x$ is the distance of this points from the first point of contact. The result obtained with N2 at the finest mesh level m6 is used as a reference. The corresponding maximum contact pressure $p_{N_{ref}}$ and the contact area $a_{ref}$ are used for the normalization purpose. It is noted that the results obtained with N2 based discretization for different meshes are in agreement to those reported in literature [18, 21, 19, 22, 27, 40, 26, 31] for Hertzian-type contact problem. In case of VO based discretizations, the quadratic order of NURBS that are sufficient to represent the considered CAD geometry exactly are kept fixed for the bulk description. For contact layer, one and two steps of additional order elevation refinements that yields N${}_{2}-$N2⋅1 and N${}_{2}-$N2⋅2 discretizations, respectively, are performed. The total number of DOFs present in the contact interface and in the remaining bulk region with both the standard and VO based NURBS discretizations for each mesh are listed in Table 4.

From Fig 18, it can be observed that for the same mesh resolution, VO based discretizations provide more accurate results than with the standard N2 based discretization. With N2, the oscillation of the contact pressure near the contact interface boundary is clearly visible at the coarse mesh m3, see Fig 18b. As the mesh is refined, the magnitude and extent of the oscillation reduces, as can be seen from Figs. 18a, 18d, and 18g. With a very fine mesh m5, although the contact pressure curve seems to match exactly with the reference solution, a slight kink is still present near the contact interface boundary. It is noted that in case of Hertzian-type contact, the quality of the solution is often affected by the sudden change of contact status from active to non-active within an element that lies across the edge of contact region. To this date, a number of solution approaches have been introduced to fix this issue, see e.g [76, 27]. Moreover, T-splines based discretization that allows adaptive mesh refinement delivers higher accurate solution than standard NURBS for a fixed number of DOFs, see [40]. However, in the present work, we particularly focus on alleviating such an error with simpler NURBS based discretizations that exclude the remeshing strategies and is computationally inexpensive.

From Figs. 18a-18c it can be observed that with VO based discretizations, especially with N${}_{2}-$N2⋅2, superior quality result is obtained even at the coarse mesh m3 as compared to N2. It is attributed to the additional number of DOFs present across the contact interface with N2⋅2 over N2. This, improves the resolution of edge region at the fixed mesh. At the intermediate mesh level m4, the contact pressure distribution with N${}_{2}-$N2⋅2 seems to match very well with the reference solution, see Fig. 18f. This is impressive as N${}_{2}-$N2⋅2 takes only half of the elements to deliver the solution comparable to N2 at mesh m5, see Fig. 18g. Further, from Figs. 18e-18f and 18h-18i it can be seen that the quality of results remain unchanged on further increasing the mesh resolution.

Next, in Fig. 19, we quantitatively show the comparison between the error in the normalized contact pressure distribution for both the standard and VO based discretizations for each mesh level. Moreover, the error decay over these meshes with the classical linear order of Lagrange discretization, i.e. L1, is also included for comparison purpose. The four noded element mesh describing the circular geometry in an approximate form is generated through the conversion of NURBS elements for all mesh level, i.e. m1 to m5. For each mesh, the error is calculated using the $L_{2}$-norm

[TABLE]

where $\bar{p}_{N_{ref}}$ is the normalized reference normal contact pressure distribution that is computed with N2 at mesh m6. In the current work, the spatial rate of convergence is not monitored as the number of elements along the radial direction is fixed during the mesh refinement (driven by $9n\times 48$). Here, we are rather interested in the difference between the accuracy of the result with VO and standard NURBS based discretizations on increasing the mesh resolution along the contact boundary layer of the geometry. It is highlighted that for the unilateral contact problems, the spatial convergence rate is generally limited by the reduced regularity of the solution under uniform mesh refinement, see e.g. [31] for the detailed description on the convergence rate with a different order of standard NURBS discretizations using the dual mortar based contact algorithm.

Figure 19 clearly shows the advantage of using VO as compared to standard N2 and L1 based discretizations for the same mesh level. The error curve with N${}_{2}-$N2⋅2 lies below the N${}_{2}-$N2⋅1, N2, and L1 curves. This is attributed to additional number of DOFs present across contact interface with N${}_{2}-$N2⋅2, see Table 4 for DOFs density data. A closer look reveals that accuracy in the result with N${}_{2}-$N2⋅2 at mesh m3 is comparable to N2 using mesh m4. Further, at mesh m4, N${}_{2}-$N2⋅2 provides a much more accurate result than with N2 using a finer mesh m5. As compared to L1, a considerable gain in accuracy is obtained with N${}_{2}-$N2⋅2 for the same mesh level. The result with N${}_{2}-$N2⋅2 using intermediate mesh m3 are of similar accuracy to that obtained with L1 at a very fine mesh m5.

Besides, a comparison with the higher than $C^{1}-$continuous based NURBS discretizations, i.e. with N${}_{2}-$N4, its counterpart standard N4, and N${}_{2}-$N3⋅2, is also carried out and the obtained results are shown in Fig. 20. It can be seen that although on increasing the continuity from $C^{0}$ (with L1) to $C^{1}$ (with N2), a more accurate result can be obtained at a fixed mesh, i.e. at intermediate m4 and fine m5 mesh level. But, with the more than $C^{1}-$continuous NURBS discretization, i.e. with N${}_{2}-$N4 or N4, the absolute error values are nearly equivalent to those obtained with N2 for the same mesh. This observation also holds for the higher-continuous version, i.e. N${}_{2}-$N3⋅2, of the N${}_{2}-$N2⋅2 discretization. Consequently, it shows that NURBS functions smoother than $C^{1}$ are not advantageous to analyze the Hertzian contact problem. This behaviour is due to the inability of higher smooth NURBS functions to adequately capture the local changes in the distribution of the contact pressure result.

Finally, we compare the associated computational efforts in terms of the overall analysis time with both the discretization approaches for similar accuracy. Figure 21 shows the plots for $L_{2}$-norm of the error in the contact pressure with the total analysis time with meshes m1 to m5 for L1, N2, N${}_{2}-$N2⋅1 and N${}_{2}-$N2⋅2 discretizations. Here, the time in $\%$ is computed using Eq. (36) and the maximum analysis time, that is used for normalizing the time values, is with N2 at mesh m6. From Fig. 21 it is evident that to attain the accuracy similar to N2, N${}_{2}-$N2⋅1 and N${}_{2}-$N2⋅2 require much lower computational efforts. With a closer look, it turns out that N${}_{2}-$N2⋅2 at mesh m3 and m4 take approximately $42.48\%$ and $43.74\%$ lower analysis time than that with N2 at m4 and m5, respectively, to attain the similar accuracy. Subsequent comparison with L1 discretization shows that to attain the comparable accuracy, N${}_{2}-$N2⋅2 at mesh m1 and m2 takes a slightly more analysis time as compared to L1 at m2 and m3, respectively. Further, N${}_{2}-$N2⋅2 at mesh m3 requires approximately $34.86\%$ lower computational efforts as compared to L1 at m5 for the similar accuracy level. The additional usage of a fine mesh, e.g. m4 or m5, with N${}_{2}-$N2⋅2 leads to a significant gain in the accuracy, which, however, comes at the expense of additional computational efforts as compared to L1 at m5.

In summary, this example shows the much higher accuracy and the computational efficiency of the VO based NURBS discretizations as compared to standard N2 based discretization. For the fixed mesh, VO NURBS yields much more accurate results at a cost slightly more than with the standard NURBS based discretization. To deliver a similar accuracy, it takes much lower computational efforts as compared to N2 discretization.

5.3 Two elastic rings contact problem

In the final example, the frictional contact between the two deformable rings undergoing the large deformation is simulated. The setup of the problem, taken from [49], is shown in Fig. 22. The outer and inner radii of the upper and lower rings are $(10,8)$ and $(12,10)$ m, respectively. The bottom surface of the lower ring is held stationary by fixing the displacement in both directions. A Neo-Hookean hyper-elastic material model, according to Eq. (29), is considered. The material parameters are: Young’s modulus $E_{\textrm{upper}}=1\times 10^{2}$ N$/$m2 and $E_{\textrm{lower}}=3\times 10^{2}$ N$/$m2, and Poisson’s ratio $\nu=0.3$ for both the rings. The plane strain assumption is used. Due to the symmetry with respect to the $y-$axis, only half of the setup is utilized. The bottom surface of the upper ring is used as a slave surface, while the outer top surface of the lower one as the master surface. A vertical displacement $\bar{u}_{y}=4.0$ m is applied in the downward direction to the top surface of the upper ring in $40$ uniform load steps. The penalty parameters are taken as $\epsilon_{\text{N}}=100$ and $\epsilon_{\text{T}}=10$. The coefficient of friction $\mu_{f}=0.1$ is used. The objective is to predict the normal and tangential reaction responses at the contact interface of the two rings using the VO based NURBS discretizations. The results obtained with N2 based discretization are used for comparisons.

A coarse mesh having $12\times 3$ number of elements in the angular and radial directions of the upper geometry, and $20\times 3$ of the lower one is shown in Fig. 23a. Overall, three different meshes are used to carry out this simulation. These meshes have $(12,20),(24,40)$, and $(48,80)$ number of elements across the contact boundary layer of the upper and lower geometries, and are denoted by m1, m2, and m3, respectively. Numerical investigation reveals that N2 based discretization using fine mesh m3 delivers the converged normal and tangential contact pressure profiles for this example. Thus, the obtained result are used as the reference solutions. The deformed shapes of the setup with the vertical displacement field $u_{y}$ at four different values of load steps with N2 using mesh m1 are shown in Figs. 23b-23e.

The normal contact pressure pN and tangential contact pressure pT with standard N2 and N4, and VO based N${}_{2}-$N4 and N${}_{2}-$N2⋅1 discretizations using meshes m1 and m2 are shown in Fig. 24. The corresponding $L_{2}$-norm of the errors in pN and pT, calculated using Eq. (37), are summarized in Table 5. The overall computational time taken by different discretizations using both the meshes are provided in Table 6. From Fig. 24 two observations are made. First, the contact pressure curves of pN and pT with higher-continuous VO based N${}_{2}-$N4 are nearly identical to those obtained with standard N4 discretization for meshes m1 and m2. However, only a marginal improvement in the accuracy of results, attributed to higher smooth NURBS functions, is obtained as compared to N2 discretization. Second, concerning the distribution of contact pressure across the contact interface at a fixed mesh, N${}_{2}-$N2⋅1 is able to capture it more accurately than N2 and N${}_{2}-$N4 (or N4) based discretizations. As can be seen from Tables 5 and 6, it improves the accuracy of the normal contact pressure distribution by a factor of $3.4$ and $9.9$, and the tangential contact pressure distribution by $3.4$ and $7.7$ at a slightly more cost, i.e. $5.8\%$ and $7.3\%$, than with the N2 based discretization with meshes m1 and m2. In summary, this example demonstrates the benefit of using VO based NURBS discretization as compared to standard NURBS discretization. It provides a major gain in the accuracy at a cost only slightly more than with the standard NURBS discretization.

6 Conclusion

In this work, a novel varying-order based NURBS discretization methodology is introduced for isogeometric contact analysis. With the proposed method, a much higher numerical accuracy is achieved in the contact computations even at a coarse mesh as compared to existing fixed-order based NURBS discretization. The improvement in the accuracy is obtained for both the small and large deformation contact problems, considering with or without friction. In other words, the proposed method exhibits a major gain in computational efficiency over the standard NURBS based discretizations for the similar accuracy level. The gain in the accuracy is attributed to the application of additional order-elevation based refinement strategy to the contact boundary layer of the NURBS discretized geometry. With this, a large number of additional DOFs are introduced across the contact interface than with the standard NURBS discretizations at a fixed mesh. The accuracy further improves on further elevating the interpolation-order of the contact boundary layer. As shown in the ironing example, the quality of the result improves marginally on additionally increasing the inter-element continuity of order-elevation discretized contact layer. Refining the geometry with the proposed method is quite robust, as illustrated by the numerical examples. The simplicity of the current approach allows itself to be conveniently embedded into the existing isogeometric contact codes. The extension of the developed methodology to three-dimensional contact problems in the framework of isogeometric analysis is currently underway.

Acknowledgments

The authors gratefully acknowledge the support from SERB, DST, under projects SB/FTP/

ETA-0008/2014 and IMP/2019/000276. The authors would like to express their gratitude towards Prof. Roger A. Sauer, AICES RWTH Aachen University, and the anonymous reviewers for their helpful comments and suggestions. The first author also gratefully acknowledges the support from IIT Guwahati and Ministry of Human Resource and Development (MHRD), Government of India, for the financial assistantship.

Appendix A Contact tangent matrices

In this appendix, we provide the contact tangent matrices for the full-pass contact formulation used for the isogeometric analysis of the above numerical examples. In the following, the subscript $n+1$, denoting the contact quantities in the current configuration, is dropped for convenience. The presented expressions for the elemental contact tangent matrices $\mathbf{k}_{\textrm{c}}^{kke}$ (where $k=\{s,m\}$) are obtained from the linearization of the contact force vector $\mathbf{f}_{\textrm{c}}^{ke}$. We refer to the works by Sauer and De Lorenzis [56, 57] for the detailed description of the formulation, including the linearization procedure.

[TABLE]

where $\mathbf{N}^{\textrm{s}}=\mathbf{N}(\xi^{\textrm{s}})$ and $\mathbf{N}^{\textrm{m}}=\mathbf{N}(\bar{\xi}^{\textrm{m}})$. The normal and tangential components of the contact traction vector $\bm{t}^{\textrm{s}}:=\bm{t}^{\textrm{s}}_{\textrm{N}}-\bm{t}^{\textrm{s}}_{\textrm{T}}$ are defined in Eqs. (17) and (25) for frictionless and frictional contact, respectively. The derivatives of $\bm{t}^{\textrm{s}}$ are given as

[TABLE]

The derivative of normal component of contact traction vector for each contact surface is given by

[TABLE]

with $c^{11}=1/[m_{11}-\textrm{g}_{\textrm{N}}\bar{\bm{n}}\cdot\bar{\bm{x}}_{,11}]$. The derivative of tangential component of contact traction vector for a stick step are given by

[TABLE]

and for a slip step are

[TABLE]

Moreover, the derivative of closest projection point are defined as

[TABLE]