Adaptive local surface refinement based on LR NURBS and its application to contact

TL;DR

This paper introduces an adaptive local surface refinement method using LR NURBS, enabling precise modeling and efficient contact analysis of complex 3D geometries with reduced computational effort.

Contribution

It develops a new LR NURBS-based refinement technique with an automatic remeshing algorithm for improved contact simulations in 3D solids and membranes.

Findings

LR NURBS achieve high accuracy with less computation than uniform refinement.

The method effectively models complex geometries in contact problems.

Numerical examples demonstrate the efficiency and accuracy of the approach.

Abstract

A novel adaptive local surface refinement technique based on Locally Refined Non-Uniform Rational B-Splines (LR NURBS) is presented. LR NURBS can model complex geometries exactly and are the rational extension of LR B-splines. The local representation of the parameter space overcomes the drawback of non-existent local refinement in standard NURBS-based isogeometric analysis. For a convenient embedding into general finite element code, the B\'ezier extraction operator for LR NURBS is formulated. An automatic remeshing technique is presented that allows adaptive local refinement and coarsening of LR NURBS. In this work, LR NURBS are applied to contact computations of 3D solids and membranes. For solids, LR NURBS-enriched finite elements are used to discretize the contact surfaces with LR NURBS finite elements, while the rest of the body is discretized by linear Lagrange finite elements. For membranes, the entire surface is discretized by LR NURBS. Various numerical examples are shown, and they demonstrate the benefit of using LR NURBS: Compared to uniform refinement, LR NURBS can achieve high accuracy at lower computational cost.

Tables

Table 1. Table 5.1: Rigid sphere in contact with a square membrane sheet: Comparison of the mesh parameters and the computational error between the LR mesh and the uniformly refined mesh. Cubic discretization is considered. Corresponding refinement depths have almost equal error but hugely different number of dofs.

LR mesh	dofs	$♯$ el.	rel. error
depth 0	$96$	$16$	$1.95\cdot {10}^{-1}$
depth 1	$280$	$64$	$1.72\cdot {10}^{-2}$
depth 2	$416$	$112$	$1.45\cdot {10}^{-3}$
depth 3	$701$	$211$	$1.09\cdot {10}^{-4}$
depth 4	$1340$	$430$	$1.63\cdot {10}^{-5}$

Equations

$$\sum _{k=\mathrm{I}}^{\mathrm{II}}\left[\delta {\Pi }_{\mathrm{int}}^k-\delta {\Pi }_{\mathrm{c}}^k-\delta {\Pi }_{\mathrm{ext}}^k\right]=0,$$

$$\begin{array}[]{lll}\delta\Pi^{k}_{\mathrm{int}}=\displaystyle\begin{cases}\displaystyle\int_{\mathcal{S}_{k}}\delta\boldsymbol{\varphi}^{k}_{,\alpha}\,\cdot\,{\sigma_{k}^{\alpha\beta}}\,\boldsymbol{a}^{k}_{\beta}\,\text{d}a_{k}~{},\quad\forall\,\delta\boldsymbol{\varphi}_{k}\in\mathcal{V}_{k}~{},&\text{for}\quad\mathcal{S}_{k}~{},\\[11.38109pt] \displaystyle\int_{\mathcal{B}_{k}}\text{grad}(\delta\boldsymbol{\varphi}_{k})\,:\,\boldsymbol{\sigma}_{k}\,\text{d}v_{k}~{},\quad\forall\,\delta\boldsymbol{\varphi}_{k}\in\mathcal{V}_{k}~{},\displaystyle&\text{for}\quad\mathcal{B}_{k}~{}.\end{cases}\end{array}$$

$$\begin{array}[]{lll}\delta\Pi^{k}_{\mathrm{c}}=\displaystyle\int_{\Gamma^{k}_{\mathrm{c}}}\delta\boldsymbol{\varphi}_{k}\cdot\boldsymbol{t}^{k}_{\mathrm{c}}\,\text{d}a_{k}~{},\quad\forall\,\delta\boldsymbol{\varphi}_{k}\in\mathcal{V}_{k}~{},\end{array}$$

$$\mathbf{x}\approx {\mathbf{x}}^h=\mathbf{N}{\mathbf{x}}_e.$$

$${\mathbf{v}}^{\mathbf{T}}\left[{\mathbf{f}}_{\mathrm{int}}-{\mathbf{f}}_{\mathrm{c}}\right]=0,\forall \mathbf{v}\in {\mathcal{V}}^h,$$

$$\begin{array}[]{lll}\mathbf{f}_{\mathrm{int}}=\displaystyle\begin{cases}\displaystyle\int_{\mathcal{S}_{k}}\mathbf{N}^{\mathrm{T}}_{,\alpha}\,\sigma^{\alpha\beta}\,\mathbf{N}_{,\beta}\,\text{d}a_{k}~{},&\text{for}\quad\mathcal{S}_{k}~{},\\[11.38109pt] \displaystyle\int_{\mathcal{B}_{k}}\mathbf{N}_{,\boldsymbol{x}}\,\mbox{\boldmath$\sigma$}\,\text{d}v_{k}~{},\displaystyle&\text{for}\quad\mathcal{B}_{k}~{},\end{cases}\end{array}$$

$${\mathbf{f}}_{\mathrm{c}}={\int }_{{\Gamma }_{\mathrm{c}}^k}{\mathbf{N}}^{\mathrm{T}}{\mathbf{t}}_{\mathrm{c}}\text{d}{a}_k.$$

$$\mathbf{g}={\mathbf{x}}_{\mathrm{I}}-{\mathbf{x}}_{\mathrm{II}}.$$

$${\mathbf{t}}_{\mathrm{n}}({\mathbf{x}}_k)=\{\begin{array}{ll}-{\epsilon }_{\mathrm{n}}{g}_{\mathrm{n}}{\mathbf{n}}_{\mathrm{p}},& g_{\mathrm{n}}<0,\\ 0,& g_{\mathrm{n}}\ge 0,\end{array}$$

$${\mathbf{t}}_{\mathrm{t}}=-\mu p\frac{{\dot{\mathbf{g}}}_{\mathrm{t}}}{{||\dot{\mathbf{g}}}_{\mathrm{t}}||}.$$

$${\mathbf{g}}_{\mathrm{t}}=\Delta {\mathbf{g}}_{\mathrm{e}}+{\mathbf{g}}_{\mathrm{s}},$$

$$f_{\mathrm{s}}=||{\mathbf{t}}_{\mathrm{t}}||-\mu p,$$

$$N_{i}^{0}(\xi)=\begin{cases}1~{},&$if $\xi_{i}\leq\xi<\xi_{i+1}~{},\\ 0~{},&$otherwise~{},$\end{cases}$$

$$N_i^p(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_i^{p-1}({\xi }_i,...,{\xi }_{p+i})+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1}^{p-1}({\xi }_{i+1},...,{\xi }_{p+i+1}).$$

$$B_{ij}^{pq}(\xi ,\eta )=N_i^p(\xi )M_j^q(\eta ).$$

$$\mathbf{x}(\xi ,\eta )=\sum _{i=1}^n\sum _{j=1}^m{B}_{ij}^{pq}(\xi ,\eta ){\mathbf{x}}_{ij},$$

$$\mathbf{x}(\xi ,\eta )=\sum _{i=1}^n\sum _{j=1}^m{B}_{ij}^{pq}(\xi ,\eta ){\mathbf{x}}_{ij}{\gamma }_{ij}.$$

$$\begin{array}[]{lll}\Xi_{1}&=\left[\xi_{1},...,\xi_{i-1},{\hat{\xi}},\xi_{i},...,\xi_{p+1}\right]~{},\\[5.69054pt] \Xi_{2}&=\left[\xi_{2},...,\xi_{i-1},{\hat{\xi}},\xi_{i},...,\xi_{p+2}\right]~{}.\end{array}$$

$$\gamma N_{\Xi }^p(\xi )={\gamma }_1{N}_{{\Xi }_1}^p(\xi )+{\gamma }_2{N}_{{\Xi }_2}^p(\xi ).$$

$$\begin{array}[]{lll}\gamma_{1}&=\alpha_{1}\,\gamma~{},\\ \gamma_{2}&=\alpha_{2}\,\gamma~{}.\end{array}$$

$$\begin{array}[]{lll}\alpha_{1}=\begin{cases}1~{},&\xi_{p+1}\leq\hat{\xi}<\xi_{p+2}~{},\\ \displaystyle\frac{\hat{\xi}-\xi_{1}}{\xi_{p+1}-\xi_{1}}~{},&\xi_{1}<\hat{\xi}<\xi_{p+1}~{},\\[11.38109pt] \end{cases}\\ \alpha_{2}=\begin{cases}\displaystyle\frac{\xi_{p+2}-\hat{\xi}}{\xi_{p+2}-\xi_{2}}~{},&\xi_{2}<\hat{\xi}<\xi_{p+2}~{},\\ 1~{},&\xi_{1}<\hat{\xi}\leq\xi_{2}~{}.\end{cases}\end{array}$$

$$\begin{array}{rl}\gamma B_{\mathbf{\Xi }}^{pq}(\xi ,\eta )& =\gamma N_{\Xi }^p(\xi )M_{\mathcal{H}}^q(\eta ),\\ & =\gamma ({\alpha }_1{N}_{{\Xi }_1}^p(\xi )+{\alpha }_2{N}_{{\Xi }_2}^p(\xi ))M_{\mathcal{H}}^q(\eta ),\\ & ={\gamma }_1{B}_{{\mathbf{\Xi }}_1}^{pq}(\xi ,\eta )+{\gamma }_2{B}_{{\mathbf{\Xi }}_2}^{pq}(\xi ,\eta ),\end{array}$$

$${\mathbf{x}}_j=\frac{{\mathbf{x}}_j{\gamma }_j+\mathbf{x}\gamma {\alpha }_j}{{\gamma }_j+{\alpha }_j\gamma },$$

$${\gamma }_j={\gamma }_j+\gamma {\alpha }_j.$$

$$({\mathbf{x}}_i{)}_k=\frac{({\mathbf{x}}_i^{\mathrm{w}}{)}_k}{{\mathrm{w}}_i},k=1,...,d,\text{where}{\mathrm{w}}_i=({\mathbf{x}}_i^{\mathrm{w}}{)}_{d+1},$$

$$W(\xi ,\eta )=\sum _{i=0}^n\sum _{j=0}^m{B}_{ij}^{pq}(\xi ,\eta ){\mathrm{w}}_{ij}.$$

$$R_{ij}^{pq}(\xi ,\eta )=\frac{B_{ij}^{pq}(\xi ,\eta ){\mathrm{w}}_{ij}}{W(\xi ,\eta )}=\frac{B_{ij}^{pq}(\xi ,\eta ){\mathrm{w}}_{ij}}{{\sum }_{\hat{i}=0}^n{\sum }_{\hat{j}=0}^m{B}_{\hat{i}\hat{j}}^{pq}(\xi ,\eta ){\mathrm{w}}_{\hat{i}\hat{j}}}.$$

$$\mathbf{x}(\xi ,\eta )=\sum _{i=0}^n\sum _{j=0}^m{R}_{ij}^{pq}(\xi ,\eta ){\mathbf{x}}_{ij}{\gamma }_{ij}.$$

$${\mathbf{x}}_{ij}^{\mathrm{w}}=[x_{ij}{\mathrm{w}}_{ij},y_{ij}{\mathrm{w}}_{ij},z_{ij}{\mathrm{w}}_{ij},{\mathrm{w}}_{ij}].$$

$${N}_{a}^{e}({\mbox{\boldmath$\xi$}})={\mathbf{c}}_{a}^{e}\,\mathbf{B}({\mbox{\boldmath$\xi$}})~{}.$$

Full text

**Adaptive local surface refinement based on LR NURBS and its application to contact

**

Christopher Zimmermann, Roger A. Sauer111corresponding author, email: sauer@aices.rwth-aachen.de

*Aachen Institute for Advanced Study in Computational Engineering Science (AICES),

RWTH Aachen University, Templergraben 55, 52062 Aachen, Germany*

Published222This pdf is the personal version of an article whose final publication is available at http:/​/link.springer.com in Computational Mechanics, DOI: 10.1007/s00466-017-1455-7

Submitted on 18. April 2017, Revised on 13. July 2017, Accepted on 21. July 2017

Abstract

A novel adaptive local surface refinement technique based on Locally Refined Non-Uniform Rational B-Splines (LR NURBS) is presented. LR NURBS can model complex geometries exactly and are the rational extension of LR B-splines. The local representation of the parameter space overcomes the drawback of non-existent local refinement in standard NURBS-based isogeometric analysis. For a convenient embedding into general finite element code, the Bézier extraction operator for LR NURBS is formulated. An automatic remeshing technique is presented that allows adaptive local refinement and coarsening of LR NURBS.

In this work, LR NURBS are applied to contact computations of 3D solids and membranes. For solids, LR NURBS-enriched finite elements are used to discretize the contact surfaces with LR NURBS finite elements, while the rest of the body is discretized by linear Lagrange finite elements. For membranes, the entire surface is discretized by LR NURBS. Various numerical examples are shown, and they demonstrate the benefit of using LR NURBS: Compared to uniform refinement, LR NURBS can achieve high accuracy at lower computational cost.

Keywords: Adaptive local refinement, computational contact mechanics, isogeometric analysis, LR B-splines, nonlinear finite element analysis

1 Introduction

A wide range of engineering applications, especially those that are governed by local surface effects, necessitate an accurate surface description. An example are contact problems. Due to their strong nonlinear behavior, those are usually only solvable by numerical methods. For this, isogeometric finite element discretizations are advantageous, since they offer at least $C^{1}$-continuity across element boundaries. Standard NURBS based discretizations lack local refinement, but their ability to discretize complex geometries exactly makes them powerful.

Isogeometric Analysis (IGA) was developed by Hughes et al., (2005) and the work of Cottrell et al., (2009) summarizes the concept of IGA. Borden et al., (2011) introduced the Bézier extraction operator, which allows a suitable embedding of isogeometric analysis into existing finite element codes. The work of Scott et al., (2011) introduces the extension of the Bézier extraction operator to T-splines. This offers a local representation of the Bézier extraction operator. B-splines and NURBS are the most widespread element types in IGA. Those element types are only globally refinable. Hierarchical B-splines were introduced by Forsey and Bartels, (1988) and allow local refinement. An analysis-suitable approach that bases on hierarchichal NURBS discretization can be found in, e.g. Schillinger et al., (2012). The work of Sederberg et al., (2003) introduces T-splines as an approach to discretize surfaces more efficiently than hierarchical B-splines. Scott et al., (2012) developed a local refinement approach based on analysis-suitable T-splines.

Locally Refined B-splines (LR B-splines) are developed as a new approach to allow local refinement in IGA. LR B-splines were introduced by Dokken et al., (2013) and further advanced by Bressan, (2013) and Johannessen et al., (2014). Johannessen et al., 2015b studied the similarities and differences of LR B-spline, classical hierarchical and truncated hierarchical discretizations. The classical hierarchical approach leads to a much denser sparsity pattern of matrices than LR B-spline or truncated hierarchical discretizations. Whether an LR B-spline or a truncated hierarchical discretization is beneficial is depending on the ratio of the locally refined domain w.r.t. the entire mesh. In the IGA community, T-splines are more widespread than LR splines and commercial software for the automatic mesh generation is available, e.g. Rhinoceros (McNeel,, 2012) with a T-splines plug-in (Autodesk,, 2015). In contrast to T-splines, the refinement by LR B-splines is directly performed in the parameter domain, which reflects the piecewise polynomial structure. This is more convenient than the multiple vertex grids of T-splines. With an LR B-spline description new possibilites are achieved on how the computational mesh can be generated. This makes them interesting in many fields of engineering. In the works of Nørtoft and Dokken, (2014) and Johannessen et al., 2015a the LR B-splines are successfully applied to the computation of Navier-Stokes flows. To the best of the authors’ knowledge rational basis functions in the context of LR splines have not been used so far. In this work, LR B-splines are extended to LR NURBS to combine the advantages of local refinement and the ability of modeling complex geometries exactly.

The numerical examples presented here consider 3D solids and membranes using a classical finite element formulation for solids (Wriggers,, 2008) and the membrane finite element formulation of Sauer et al., (2014). This work is based on computational contact mechanics (Laursen,, 2002; Wriggers,, 2006) using isogeometric finite elements (De Lorenzis et al.,, 2011; Lu,, 2011; Temizer et al.,, 2011; De Lorenzis et al.,, 2012) and an unbiased contact algorithm (Sauer and De Lorenzis,, 2013, 2015). In the literature several approaches can be found that account for local refinement and adaptive meshing. Lee and Oden, (1994) present an adaptive method for $h$-$p$ refined meshes in the context of frictional contact and non-isogeometric finite elements. The work of Hager et al., (2012) presents a non-conforming domain decomposition method. The domain decomposition consists of a global coarse mesh and an overlapping fine mesh for the contact computation. In the context of IGA and frictionless contact a recent work based on analysis-suitable T-splines is presented in Dimitri et al., (2014). The work by Dimitri and Zavarise, (2017) is based on analysis-suitable T-splines and investigates frictional contact and mixed mode debonding problems. In the work of Temizer and Hesch, (2016) hierarchical NURBS are used to allow local refinement. In this work, LR NURBS are applied to frictionless and frictional sliding contact.

LR B-splines have been successfully used in 2D, since they are linearly independent in 2D. In 3D, the linear independence for arbitrary, locally refined meshes has not been proven yet and requires further study. To overcome this issue, this work uses 2D LR NURBS and combines them with an enrichment technique (Corbett and Sauer,, 2014, 2015) along the third direction. This combines the high accuracy, which is achieved by isogeometric analysis and the efficiency of standard finite elements. This formulation is adapted for LR NURBS-enrichment to gain the possibility of local refinement. In this work we present a novel technique for adaptive local surface refinement and coarsening that

• •

considers the extension of LR B-spline finite elements to LR NURBS finite elements

• •

uses a local formulation of the Bézier extraction operator

• •

is automatically controlled by a proposed refinement criterion

• •

is extended to LR NURBS-enriched finite elements for 3D solids

• •

is applied to frictionless and frictional sliding contact.

In Sec. 2.1 a brief summary of the nonlinear finite element formulation for membranes and 3D solids is given. A short introduction to computational contact mechanics follows in Sec. 2.2. The fundamentals of LR B-splines are presented in Sec. 2.3. The geometric modeling using LR NURBS and the formulation of the Bézier extraction operator is presented in Sec. 3. The technique of adaptive local refinement and coarsening using LR NURBS discretizations is presented in Sec. 4. Numerical results are shown in Sec. 5. The performance of LR NURBS within the proposed local refinement technique is compared with uniform meshes using NURBS discretizations. A conclusion is given in Sec. 6.

2 Preliminaries

This section presents an overview of the fundamental mathematical formulations that are used within this work. A general finite element formulation for membranes and 3D solids in the context of contact computations is discussed in Sec. 2.1 and Sec. 2.2. The fundamentals of LR B-splines are given in Sec. 2.3.

2.1 Finite element formulation

The finite element formulation for nonlinear membranes and 3D solids is similar. Their governing equations (e.g. equilibrium) are different but obtaining their discretized weak form follows the same approach. The quasi-static weak form of a system of two deformable bodies in contact is described by

[TABLE]

considering two objects $(k=\mathrm{I},\mathrm{II})$333the index $k$ is used as super- and subscript but the position has no special meaning. This describes the equilibrium between the internal virtual work $\delta\Pi^{k}_{\mathrm{int}}$, the virtual contact work $\delta\Pi^{k}_{\mathrm{c}}$ and the external virtual work. In our examples $\delta\Pi^{k}_{\mathrm{ext}}~{}=~{}0$ is considered. The computational formulation for nonlinear membranes in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{3}$ in the framework of curvilinear coordinates is taken from Sauer et al., (2014). A general constitutive setting and finite element formulation for solids in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{3}$ can be found in Wriggers, (2008). The internal virtual work for membranes $\mathcal{S}_{k}$ in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{3}$ and solids $\mathcal{B}_{k}$ in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{3}$ is

[TABLE]

Here, $\sigma^{\alpha\beta}_{k}$ are the contra-variant components of the Cauchy stress tensor $\sigma$ with $\alpha=1,2$ and $\beta=1,2$. The parametric derivative of the virtual displacement field is denoted by $\delta\boldsymbol{\varphi}^{k}_{,\alpha}$ and $\boldsymbol{a}^{k}_{\alpha}$ are the co-variant tangent vectors. The virtual contact work for the contact surfaces $\Gamma^{k}_{\mathrm{c}}$ is expressed as

[TABLE]

where $\boldsymbol{t}^{k}_{\mathrm{c}}$ denotes the contact traction on the surface $\Gamma^{k}_{\mathrm{c}}\subset\mathcal{S}_{k}$. The computation of the contact integrals is performed using the two-half pass algorithm of Sauer and De Lorenzis, (2013). This algorithm evaluates Eq. (2.3) on each body equivalently. System (2.1) is solved for the unknown deformation field $\boldsymbol{\varphi}_{k}\in\mathcal{V}_{k}$, with $\mathcal{V}_{k}$ as a suitable space of kinematically admissible variations. The membrane and the contact surface are discretized into a set of membrane elements $\Gamma^{k}_{e}$ and a set of contact elements $\Gamma^{k}_{\mathrm{c}\,e}$ such that $\mathcal{S}_{k}\approx\mathcal{S}_{k}^{h}=\bigcup_{e}\Gamma^{k}_{e}$ and $\Gamma^{k}_{\mathrm{c}}\approx\Gamma^{k\,h}_{\mathrm{c}}=\bigcup_{e}\Gamma^{k}_{\mathrm{c}\,e}$, respectively. Here, the element domains $\Gamma^{k}_{e}$ and $\Gamma^{k}_{\mathrm{c}\,e}$ are considered equal and have the same finite element discretization. Similarly, the solids are discretized into a set of volume elements $\Omega^{k}_{e}$ and a set of surface contact elements $\Gamma^{k}_{\mathrm{c}\,e}$ such that $\mathcal{B}_{k}\approx\mathcal{B}_{k}^{h}=\bigcup_{e}\Omega^{k}_{e}$ and $\Gamma^{k}_{\mathrm{c}}\approx\Gamma^{k\,h}_{\mathrm{c}}=\bigcup_{e}\Gamma^{k}_{\mathrm{c}\,e}$, respectively. The superscript $h$ denotes the approximated discrete setting. With a set of $n_{e}$ basis functions $\mathbf{N}=\left[N_{1}\boldsymbol{I},...,N_{n_{e}}\boldsymbol{I}\right]$ with identity matrix $\boldsymbol{I}$ in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{3}$ and discrete points $\mathbf{x}_{e}$, a point $\boldsymbol{x}\in\mathcal{B}_{k}$ is interpolated within each element as

[TABLE]

The same holds for the reference configuration $\boldsymbol{X}\approx\boldsymbol{X}^{h}$, the displacement field $\boldsymbol{u}\approx\boldsymbol{u}^{h}$ and the virtual displacement field $\delta\boldsymbol{\varphi}\approx\boldsymbol{v}^{h}$. In here we drop the index $k$ for convenience. In standard finite elements the discrete points $\mathbf{x}$ are known as nodes, while in IGA they are called control points. For standard finite element basis functions Lagrange polynomials are used, while for IGA the basis functions will be discussed in Sec. 2.3. As a consequence of the above definitions, the dicretized weak form can be written as

[TABLE]

with the internal forces

[TABLE]

and the contact forces

[TABLE]

A brief overview of computational contact mechanics is given in the next section.

2.2 Computational contact formulation

This section gives a brief summary of computational contact mechanics. For more detailed information, the monograph by Wriggers, (2006) is recommended. Here, the penalty method is used in order to enforce contact. The gap vector $\boldsymbol{g}$ between two points $\boldsymbol{x}_{\mathrm{I}}$ and $\boldsymbol{x}_{\mathrm{II}}$ on the surfaces $\mathcal{S}_{\mathrm{I}}$ and $\mathcal{S}_{\mathrm{II}}$ is given by

[TABLE]

The closest projection point $\boldsymbol{x}_{\mathrm{p}}$ of $\boldsymbol{x}_{\mathrm{I}}$ is obtained by orthogonal projection of $\boldsymbol{x}_{\mathrm{I}}$ onto $\mathcal{S}_{\mathrm{II}}$. The contact traction is expressed as $\boldsymbol{t}_{\mathrm{c}}=\boldsymbol{t}_{\mathrm{n}}+\boldsymbol{t}_{\mathrm{t}}$, i.e. it is decomposed into its normal and tangential components. According to the penalty formulation, the normal traction

[TABLE]

is applied at each contact point. The traction is proportional to the normal gap $g_{\mathrm{n}}=(\boldsymbol{x}_{\mathrm{I}}-\boldsymbol{x}_{\mathrm{p}})\cdot\boldsymbol{n}_{\mathrm{p}}$, the surface outward normal $\boldsymbol{n}_{\mathrm{p}}$ at $\boldsymbol{x}_{\mathrm{p}}$ and the constant penalty parameter $\varepsilon_{\mathrm{n}}$. $k$ is the index for the body under consideration. The tangential contact traction $\boldsymbol{t}_{\mathrm{t}}$ during frictional sliding is given by Coulomb’s law

[TABLE]

Here, we find the friction coefficient $\mu$, the contact pressure $p=||\boldsymbol{t}_{\mathrm{n}}||$ and the relative tangential sliding velocity $\boldsymbol{\dot{g}}_{\mathrm{t}}$. The tangential contact slip

[TABLE]

consists of the reversible $\Delta\boldsymbol{g}_{\mathrm{e}}$ and irreversible $\boldsymbol{g}_{\mathrm{s}}$ part. The slip criterion

[TABLE]

is used to distinguish between sticking and slipping. Sliding occurs for $f_{\mathrm{s}}=0$ and sticking for $f_{\mathrm{s}}<0$. The examples presented below deal with frictionless and frictional contact. For frictionless contact, the tangential traction $\boldsymbol{t}_{\mathrm{t}}$ vanishes and the contact traction is given by $\boldsymbol{t}_{\mathrm{c}}=\boldsymbol{t}_{\mathrm{n}}$. The implementation and formulation of frictional contact is taken from Sauer and De Lorenzis, (2015).

2.3 LR B-spline fundamentals

The core idea of LR B-splines is to disband the tensor-product mesh structure of classical B-splines in order to obtain a local representation of the parameter domain. As introduced in the B-spline theory, a knot vector $\Xi$ of size $n+p+1$ generates $n$ linearly independent basis functions of order444also known as degree $p$. The local representation of the parameter space and the geometric discretization are obtained by splitting the global knot vector $\Xi=[\xi_{1},...,\xi_{n+p+1}]$ into local knot vectors $\Xi_{i}=[\xi_{i},...,\xi_{i+p+1}]$ each constructing a single basis function. With the local knot vector representation, the domain of the basis function is minimal, i.e. the basis has minimal support. The basis functions are defined by the Cox-de Boor recursion formula, which depends on the entries of the local knot vector $\Xi_{i}$ and desired polynomial order $p$, see Cox, (1971) and De Boor, (1972). For $p=0$

[TABLE]

and for $p>0$

[TABLE]

The LR B-spline basis $B^{p\,q}_{i\,j}(\xi,\eta)$ of order $p$ and $q$ in 2D is defined as a separable function $B:{\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{2}\rightarrow{\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}$

[TABLE]

$N_{i}^{p}(\xi)$ and $M_{j}^{q}(\eta)$ are the shape functions in each parametric direction. The properties of LR B-splines follow directly from standard isogeometric analysis. The basis is non-negative and forms a partition of unity, and the resulting geometry lies within the convex hull of the control points. The continuity across the element boundaries is ${C}^{p-{m_{v}}}$, i.e. it is defined by the polynomial order $p$ and the multiplicity of the knot vector entry $m_{v}$. The LR B-spline surface is then constructed by a set of control points $\mathbf{x}_{i\,j}$ and the local basis functions $B^{p\,q}_{i\,j}(\xi,\eta)$, as

[TABLE]

with $n,m$ as the number of single basis functions in each parametric dimension. To ensure that LR B-splines keep the partition of unity property during local refinement the basis functions are multiplied by a scaling factor $\gamma_{i\,j}\in(0,1]$

[TABLE]

With the above formulation, the global B-spline representation is split into a local representation consisting of a set of locally refinable B-splines. In the following we use the terminology LR B-spline for both locally refinable and locally refined B-splines. It is important to note that the global B-splines do not exist anymore. Global B-splines do not allow local refinement. The ability is only given by the LR B-splines. Locally refined meshes and the local refinement procedure are discussed in the following sections.

2.3.1 Locally refined meshes

The local representation of the parametric domain leads to new possibilities on how the finite element mesh can be constructed. The ability of local refinement is the major advantage in comparison to classical tensor meshes. The LR mesh is the result of a series of meshline insertions into an initial tensor mesh. Fig. 2.1a shows an example of an LR mesh. The meshlines never stop in the center of an element (knot span). As a knot vector of an LR B-spline basis function has size of $p+2$ knot vector entries, the meshlines cross at least $p+2$ knots. A meshline insertion can be either a new meshline, an elongation of an existing one, a joining of two meshlines or the increase of the multiplicity of a meshline. Increasing the multiplicity decreases the continuity of the LR B-spline. In the bi-variate case, horizontal and vertical meshlines can be inserted. Any type of meshline insertion implies that an LR B-spline loses the property of minimal support. An LR B-spline has minimal support if the support domain of the basis functions is not fully crossed by any other meshline. The basis function domain of the LR B-spline in Figure 2.1b, marked in gray, has minimal support because its support domain is not fully crossed by any meshline insertion completely. The LR B-spline basis marked in the Fig. 2.1c does not have minimal support. Its basis function domain is crossed by the vertical meshline, which is spanned by $\eta=[0.2,0.8]$ at $\xi=0.5$. If an LR B-spline loses the property of minimal support refinement is performed as discussed in the next section.

2.3.2 Procedure of local refinement

Local refinement is perfomed if an LR B-spline loses the property of minimal support due to meshline extensions. In general, the procedure of local refinement is realized by single knot insertion. From the classical theory it is known that this enriches the basis while the geometry remains unchanged. The insertion of a single knot $\hat{\xi}$ into the knot span \mbox{\boldmath\xi}_{i}=[\xi_{i-1},\xi_{i}] of a local knot vector $\Xi$ of size $p+2$ leads to a knot vector of size $p+3$. Consequently, two locally refined B-splines are generated by splitting the enlarged knot vector $\Xi=\left[\xi_{1},...,\xi_{i-1},\hat{\xi},\xi_{i},...,\xi_{p+2}\right]$ into the local knot vectors of size $p+2$

[TABLE]

The relation for an LR B-spline basis in one parametric direction on the LR mesh is then given by

[TABLE]

with

[TABLE]

The associated $\alpha_{1}$ and $\alpha_{2}$ are determined by

[TABLE]

Consider the knot $\hat{\xi}=0.375$ is inserted into the knot vector $\Xi=[0,0.25,0.5,0.75]$. The knot vector entries construct the quadratic basis function $N^{2}_{\Xi}$. As illustrated in Fig. 2.2, the original basis function is split into two new functions each describing the LR B-spline basis functions $\gamma\,\alpha_{1}\,N^{2}_{\Xi_{1}}$ and $\gamma\,\alpha_{2}\,N^{2}_{\Xi_{2}}$.

The refinement process for bivariate functions is performed in one parametric direction at a time, i.e. first all horizontal meshlines and then all vertical meshlines are inserted (or vice versa). Consider two local knot vectors $\Xi$ and $\mathcal{H}$, one in each parametric direction. For a weighted, bivariate LR B-spline basis ${\gamma}\,B^{p\,q}_{\boldsymbol{\Xi}}(\xi,\eta)$ with $\mathbf{\Xi}=[\Xi,\mathcal{H}]$ the following relation is given

[TABLE]

for the case that $\Xi$ is split into $\Xi_{1}$ and $\Xi_{2}$. The refinement process is given by two steps. At first all LR B-spline basis functions whose support domain is crossed by meshline extensions will be split. The second step is to check every newly created basis function if it has minimal support. In the case that a newly created LR B-spline does not have minimal support an additional splitting is performed. By splitting $\gamma\,B^{p\,q}_{\boldsymbol{\Xi}}(\xi,\eta)$ into ${\gamma_{j}}\,B^{p\,q}_{\boldsymbol{\Xi}_{j}}(\xi,\eta)$, with $j=1,2$, the following cases can occur:

The LR B-spline does not exist and a new LR B-spline is created. In this case the new control point $\mathbf{x}_{j}$ is a copy of the control point $\mathbf{x}$ of the parent LR B-spline $\mathbf{x}_{j}=\mathbf{x}$. The weight $\gamma_{j}$ is set by Eq. (2.20). 2. 2.

If $B^{p\,q}_{\boldsymbol{\Xi}_{j}}(\xi,\eta)$ already exists, the control point and the weight are set to

[TABLE]

and

[TABLE]

After splitting the former LR B-spline, $\gamma\,B^{p\,q}_{\boldsymbol{\Xi}}(\xi,\eta)$ is deleted in both cases. The algorithm proceeds with the second step and checks if the support domain of the new local basis functions is fully crossed by any existing meshline. If one does not have minimal support the first step is performed again. Note that at every step of the refinement process the partition of unity is maintained and the geometric mapping is left unchanged. For the use of LR B-splines in IGA it is required to ensure that the resulting spline space is linearly independent. In the bivariate case this can be guaranteed by only using primitive meshline extensions555A primitive meshline extension is (a) a meshline spanning $p+1$ elements, (b) elongating a meshline by one element or (c) raising the multiplicity of a meshline (length of $p+1$ elements)., see Mourrain, (2014) and Johannessen et al., (2014). In this work, all meshline extensions are formulated as primitives.

3 Geometric modeling using LR NURBS

So far, non-rational LR B-splines have been discussed. In this section, their extension to Locally Refined Non-Uniform $Rational$ B-Splines (LR NURBS) is introduced. With a NURBS representation of objects one gains the ability to describe many geometries, which cannot be represented by polynomials, exactly. Especially conic sections, ellipsoids, spheres, cylinders, etc. can be constructed by a projective transformation of piecewise quadratic curves exactly, see Farin, (1992). This is one of the defining features of NURBS.

3.1 LR NURBS

The fundamentals of standard NURBS in isogeometric analysis can be found in, e.g. Hughes et al., (2005). Combining NURBS with the LR B-spline theory from Sec. 2.3, the extension of LR B-splines to LR NURBS follows. An LR NURBS object in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{d}$ is constructed by the projective transformation of an LR B-spline entity in ${\rm\kern 2.40005pt\vrule width=0.20004pt,height=6.58746pt,depth=-0.21529pt\kern-3.00003ptR}^{d+1}$. The first $d$ entries of a projective control point $\boldsymbol{\mathbf{x}}_{i}^{\mathrm{w}}$ represent the spatial coordinates and the $d+1$ entry the weight, e.g. $\boldsymbol{\mathbf{x}}_{i}^{\mathrm{w}}=[x_{i},y_{i},z_{i},\mathrm{w}_{i}]$ for $d=3$. The control points of the LR NURBS object result from the projective transformation

[TABLE]

with $(\mathbf{x}_{i})_{k}$ as the $k^{th}$ component of the vector $\mathbf{x}_{i}$ and $\mathrm{w}_{i}$ as the $i^{th}$ weight. To generalize this relation the weighting function for the bivariate case is introduced by

[TABLE]

$B^{p\,q}_{i\,j}(\xi,\eta)$ are the standard LR B-spline basis functions and the transformation is applied by dividing every point of the curve by $W(\xi,\eta)$. Each element of surface ${\mathcal{S}}$ is a polynomial, which is divided by another polynomial of the same order. The resulting LR NURBS basis is then defined as

[TABLE]

With a set of control points $\mathbf{x}_{i\,j}$ and the scaling factors $\gamma_{i\,j}$ the LR NURBS surface is then defined as

[TABLE]

The properties of the LR NURBS basis, such as its continuity and support, follow directly from the knot vectors as before. The basis is still non-negative and it still forms a partition of unity. This leads to the strong convex hull property of the LR NURBS. Note that the weights are separated from any explicit geometric representation. The weights are each associated with a specific control point and a manipulation of them leads to a change of the resulting geometry. In the case that all weights are equal, the surface is again a polynomial and consequently $R^{p\,q}_{i\,j}(\xi,\eta)=B^{p\,q}_{i\,j}(\xi,\eta)$. Therefore, LR B-splines are a special case of LR NURBS. In order to locally refine an LR NURBS object it is necessary to take special care of the weights. Before the refinement procedure starts, the projective control points $\mathbf{x}_{i\,j}^{\mathrm{w}}$ of the LR NURBS are computed by Eq. (3.1). The weights $\mathrm{w}_{i\,j}$ are included in vector $\mathbf{x}_{i\,j}^{\mathrm{w}}$ and treated as the fourth part of the control points. This results in

[TABLE]

After the refinement process, the control points are transformed back by dividing the spatial coordinates by their associated weight. Note that while the refinement is performed, $\mathrm{w}_{i\,j}$ is treated in the same way as the control points. The linear independence of LR NURBS spaces follow directly from the LR B-splines.

3.2 Bézier extraction operator for LR NURBS

This section discusses the Bézier extraction of LR NURBS. The Bézier extraction is an advantageous technique for the decomposition of splines into their Bézier elements. It allows for the embedding of isogeometric analysis into a standard finite element framework. The Bézier extraction of NURBS is introduced by Borden et al., (2011) and advanced to T-splines by Scott et al., (2011). The Bézier extraction of LR NURBS shares similarities with the formulation of the Bézier extraction of T-splines. Both element types have a local parametric representation. The Bézier decomposition of each local basis function, that is nonzero within an element, is performed separately. The basis function in one parametric direction of an LR NURBS element $\Gamma^{e}_{\square}$ within the knot span $\xi$ can be expressed in terms of a set of Bernstein polynomials \mathbf{B}({{\mbox{\boldmath\xi}}}) and a linear operator $\mathbf{c}_{a}^{e}$. Each local basis function {N}_{a}^{e}({\mbox{\boldmath\xi}}) of an element $\Gamma^{e}_{\square}$ can then be expressed by

[TABLE]

Here, $a=1,2,...,n_{e}$ where $n_{e}$ is the number of nonzero basis functions of $\Gamma^{e}_{\square}$. In matrix-vector form this becomes

[TABLE]

$\mathbf{C}^{e}$ is the element-wise Bézier extraction operator. The linear operator ${\mathbf{c}}_{a}^{e}$ is of dimension $[1\times(p+1)]$ and represents a single row of $\mathbf{C}^{e}$. The coefficients of ${\mathbf{c}}_{a}^{e}$ are calculated by an algorithm presented in Scott et al., (2011). The algorithm inserts knots into the knot vector until all entries have multiplicity $p$. The knot insertion performs the spline decomposition into its Bézier elements. With ${\mathbf{c}}_{a}^{e}$ and the Bernstein polynomials the relation in Eq. 3.6 is obtained. Due to the unstructured LR meshes it is required to account for two aspects:

To determine ${\mathbf{c}}_{a}^{e}$ it is required that the knot vectors are open, i.e. the first and last knot vector entry have multiplicity of $p+1$. The knot vectors of LR NURBS are in general not open. For this, the local knot vectors are extended, which is simply done by adding knot vector entries. The extension of the knot vectors does not affect the basis of LR NURBS itself because it is only used for the Bézier decomposition. 2. 2.

The unstructured meshes of LR NURBS or T-splines necessitate a special consideration. An example is shown in Fig. 3.1. Consider the element $\Gamma^{e}_{\square}$ with \tilde{}\mbox{\boldmath\xi}=[0,0.1] marked in yellow in Fig. 3.1 on the left. The support domain of one bivariate local basis function is marked in gray. In $\xi$-direction, the basis function $N_{a}(\xi)$ is spanned by the local knot vector ${\Xi}=[0,0.2,0.4,0.6]$ as shown in Fig. 3.1 on the right. The basis function is non-zero within the domain of $\Gamma^{e}_{\square}$. From $\Xi$ follows that ${\mathbf{c}}_{a}^{e}$ is defined within the knot span {\mbox{\boldmath\xi}}=[0,0.2] and exceeds \tilde{}\mbox{\boldmath\xi}=[0,0.1] (the domain of $\Gamma^{e}_{\square}$). Those non-matching domains must be considered and treated. For T-splines an additional knot is inserted into the knot vector so that ${\mathbf{c}}_{a}^{e}$ can be determined for the correct parametric domain \tilde{}\mbox{\boldmath\xi}=[0,0.1]. In contrast to T-splines, no additional knot insertion is required for the formulation for LR NURBS.

For LR NURBS the linear operator ${\mathbf{c}}_{a}^{e}$ is mapped by using an additional linear operator $\mathbf{T}_{a}^{e}$ to obtain a new $\tilde{\mathbf{c}}_{a}^{e}$ such that

[TABLE]

With $\tilde{\mathbf{c}}_{a}^{e}={\mathbf{c}}_{a}^{e}\,\mathbf{T}_{a}^{e}$

[TABLE]

follows. The matrix $\mathbf{T}_{a}^{e}$ follows from the relation

[TABLE]

After some algebraic manipulations this equation is expressed as

[TABLE]

such that

[TABLE]

If the knot spans coincide, i.e. \tilde{}\mbox{\boldmath\xi}=\mbox{\boldmath\xi}, then

[TABLE]

[TABLE]

and $\tilde{\mathbf{c}}_{a}^{e}={\mathbf{c}}_{a}^{e}$ follows. With this formulation, the basis functions having support in each element $\Gamma^{e}_{\square}$ are expressed by the Bernstein polynomials \mathbf{B}({{\mbox{\boldmath\xi}}}) and the linear operators $\tilde{\mathbf{c}}_{a}^{e}$.

4 Adaptive local refinement and coarsening

This section presents a technique for adaptive local refinement and coarsening for LR NURBS discretizations. The technique is formulated in the context of frictional contact. This requires the mapping of contact variables from one mesh to another.

4.1 Adaptive local refinement

An adaptive local refinement technique necessitates an indicator for refinement. Commonly used indicators use a posteriori error estimation that provide a reliable error distribution, see e.g. Ainsworth and Oden, (1997), Johannessen et al., (2014) and Kumar et al., (2015). In this work, a refinement indicator based on the contact state is used. If new contact is detected in an element, refinement is considered. If contact is lost in an element, coarsening is considered. After the contact domain is determined, the local refinement is performed by meshline extensions in the parameter domain. The refinement is performed until a prescribed refinement depth or smallest element size, is obtained. For sliding contact it is desirable to automatically control the technique of adaptive local refinement and coarsening. A changing contact domain can lead to a high number of newly detected contacts and lost contacts. To reduce the number of refinement and coarsening events, three parameters are used. In Fig. 4.1 four cases are sketched to illustrate the automatic control by using the parameters $d_{\mathrm{ref}}^{d}$, $d_{\mathrm{safe}}^{d}$ and $d_{\mathrm{crs}}^{d}$. The superscript $d$ denotes the current refinement depth that is performed.

• •

The parameter $d_{\mathrm{ref}}^{d}\geq 0$ enlarges the refinement domain by a number of unrefined elements (see Fig. 4.1a). In this example, $d_{\mathrm{ref}}^{d}=3\,d_{\mathrm{e}}^{1-d}$ in both directions. Here,

[TABLE]

is the minimum element length of the previous refinement depth and $d^{0}_{\mathrm{e}}$ is the original element length. The elements within the enlarged contact domain are flagged for refinement. The LR mesh is the result of splitting all flagged elements equally.

• •

The parameter $d_{\mathrm{safe}}^{d}\geq 0$ sets a safety domain at the boundary between coarse and refined elements. If contact is detected within this domain, local refinement is performed. An example with $d_{\mathrm{safe}}^{d}=d_{\mathrm{e}}^{1-d}$ in both directions is shown in Fig. 4.1b. This parameter ensures that at every time-step the contact domain is surrounded by refined elements. $d_{\mathrm{safe}}^{d}$ detects a refinement event and does not specify a domain for refinement.

• •

The parameter $d_{\mathrm{crs}}^{d}\geq 0$ controls the coarsening. An example with $d_{\mathrm{crs}}^{d}=5\,d_{\mathrm{e}}^{1-d}$ in both directions is shown in Fig. 4.1c. Coarsening is performed if refined elements are detected beyond this distance from the contact domain. $d_{\mathrm{crs}}^{d}$ detects a coarsening event and does not specify a domain for coarsening.

Fig. 4.1d shows the three parameters combined. The parameters ensure that the contact domain and the neighboring elements are at every time-step represented in the desired discretization. The parameters $d_{\mathrm{ref}}^{d}$, $d_{\mathrm{safe}}^{d}$ and $d_{\mathrm{crs}}^{d}$ can be set arbitrarily and a good choice is depending on the specific problem setup. A good compromise of the size of the locally refined domain and the number of refinement/coarsening events is desirable. The coarsening is discussed in the next section.

4.2 Adaptive coarsening

Coarsening LR meshes in the context of contact computations is challenging. The reasons are the unstructured parametric representation of LR NURBS and that frictional contact is history dependent. The contact variables need to be preserved during the process of coarsening and refinement. A technique for NURBS and T-spline coarsening is presented in Thomas et al., (2015). This technique has not been advanced to LR NURBS yet. To overcome this here, the coarsening is combined with local refinement. Therefore, an intermediate step is introduced, in which the entire mesh is coarsened. Based on the coarse mesh, the local refinement is then performed to obtain the desired LR mesh. The coarsening and refinement process requires that the contact domain has the same discretization and identical contact variables before and after the process. Due to the intermediate step, the contact domain is temporary coarsened. To preserve the contact variables a mapping is required. The contact variables are stored before the mesh is coarsened, and then mapped to the new mesh after the refinement step. The technique consists of six main steps to perform a consistent adaptive coarsening within the context of contact computations. The technique is sketched in Fig. 4.2 and the steps are:

Coarsening is detected when refined elements are detected beyond the distance $d_{\mathrm{crs}}^{d}$ (Fig. 4.2a). Store the current configuration and the contact variables for each element in contact. 2. 2.

Recover the undeformed, reference mesh (Fig. 4.2b). 3. 3.

Interpolate the control points of the mesh in Fig. 4.2a to obtain the deformed, coarse mesh (Fig. 4.2c). The interpolation leads to a geometric approximation error, which is addressed in step 6. 4. 4.

Perform local refinement within the enlarged contact domain $d_{\mathrm{ref}}^{d}$ (Fig. 4.2d) to obtain the desired LR mesh (Fig. 4.2e). 5. 5.

Preserve the contact variables by mapping them to the current configuration. 6. 6.

Reduce the geometric approximation error by recovering control points from the configuration in Fig. 4.2a: By comparing the meshes in Fig. 4.2a and Fig. 4.2e it turns out that their discretization only differ slightly. The control points of the domains, that remain unchanged in their discretization, are recovered.

As Fig. 4.2 shows, the discretization of the contact domain, at a given time-step, is unaffected by the coarsening and refinement procedure. With the adaptive local refinement and coarsening technique, a refinement/coarsening event is not detected at every computational step. A refinement/coarsening event is performed only when it is detected according to the criteria of Fig. 4.1. The results of this technique can be seen in the following examples, e.g. Fig. 5.8.

5 Numerical examples

The performance of the adaptive local refinement and coarsening technique using LR NURBS elements is illustrated in this section by several numerical examples. The first example is used for validation. In the remaining examples, the performance of LR NURBS for frictionless and frictional contact is investigated. In the last example, LR NURBS-enriched contact elements are considered for 3D friction of two deformable solids. To investigate the benefit of LR NURBS discretization, comparisons to reference models with uniform NURBS discretization are made.

5.1 Inflation of a hemisphere

The first example considers the inflation of a spherical membrane balloon. The inflation can be described by an analytical formula, which is used for validation. The system is modeled as a hemisphere and the boundary conditions are chosen such that the symmetry of the model is maintained during inflation, see the left side of Fig. 5.1. The discretized geometry consists of five patches and the NURBS weights differ for each. This chosen geometry is not perfectly spherical but still an approximation. The material is described by the incompressible Neo-Hookean membrane material model (Sauer et al.,, 2014)

[TABLE]

Here, $\mu$ is the shear modulus, $J$ is the surface area change and $A^{\alpha\beta}$ and $a^{\alpha\beta}$ are the contra-variant components of the metric tensor in the reference and current configuration. The problem setup is similar to the one of Sauer et al., (2014) in which the balloon is modeled as $1/8$th of a sphere with NURBS and Lagrange discretizations. Initially, the balloon has the volume $V_{0}=4\pi R^{3}/3$ with the radius $R$. The volume is increased step-wise until it reaches $V=50\,V_{0}$. The analytical pressure-volume relation is given by

[TABLE]

with the internal pressure $p_{\mathrm{int}}$.

We now investigate the behavior of the proposed local refinement procedure using quadratic LR NURBS elements. Five different uniform meshes with quadratic NURBS discretizations are taken each as the starting point. The number of Gaussian quadrature points is $n_{qp}=3\times 3$ for all elements. As an homogeneous hemisphere has no domain of major interest, the center is chosen to apply local refinement of depth 2, see the right side of Fig. 5.1. The performance of the five meshes is compared to quadratic NURBS and quadratic Lagrange discretizations examining the pressure error, see Fig. 5.2. The isogeometric elements behave as observed in Sauer et al., (2014) and show better convergence behavior than standard Lagrange elements. The LR NURBS meshes behave equal well than the uniform NURBS meshes. Better behavior can not be expected in this example, since the hemisphere has a homogeneous behavior and local refinement does not improve the numerical results significantly. But with this first example we can conclude that the applied local refinement procedure works successful and the computations with LR NURBS are valid. In the following examples, the performance of adaptive local refinement with respect to the accuracy and the computational cost is investigated.

5.2 Rigid sphere in contact with a square membrane sheet

The second example considers a square membrane sheet with dimension $2\,L_{0}\times 2\,L_{0}$ that is initially pre-stretched by $\lambda=1.1$. The pre-stretching is applied to avoid the membrane instability. A rigid sphere, located at $\boldsymbol{X}=[0,0,L_{0}]$ with radius $R=L_{0}$, is pushed gradually downwards until the bottom of the sphere reaches $z=-L_{0}/2$ as shown in Fig. 5.3. The initial mesh consists of $4\times 4$ LR NURBS elements with $5\times 5$ Gaussian quadrature points each. The material behavior is the same as in the previous example. The boundaries at $X=2\,L_{0}$ and $Y=2\,L_{0}$ are fixed in all directions. At $X=0$ and $Y=0$ the deformation is constrained to be zero in $x$- and $y$-direction, respectively. Contact between the rigid sphere and the membrane is treated by the penalty method with the penalty parameter

[TABLE]

depending on the individual element size. Here, $l_{el_{x}}$ and $l_{el_{y}}$ denote the current element lengths, and $l_{0_{x}}$ and $l_{0_{y}}$ the initial element lengths. $p$ is the polynomial order666the order in both parametric directions is set equally, i.e. $q=p$ and $\varepsilon_{\mathrm{n}}^{0}=10\,E_{0}/L_{0}$ is the constant penalty parameter. A mesh and order depending penalty parameter was studied in Roohbakhshan and Sauer, (2016) and Sauer and De Lorenzis, (2015) giving a good balance between accuracy and convergence behavior.

The contact domain is locally refined by the procedure described in Sec. 4. A series of resulting meshes for quadratic LR NURBS elements is illustrated in Fig. 5.4 from refinement depth 1 to 5. A strong local aggregation of elements is obtained in the contact domain. In this example, the parameter $d_{\mathrm{ref}}^{d}$ is taken as zero.

The performance of LR NURBS meshes and uniform meshes using NURBS discretizations is investigated by examining the relative error

[TABLE]

of the normal contact force $f_{\mathrm{n}}$ w.r.t. a reference solution. For the reference solution, a uniform mesh with cubic NURBS discretizations is taken. This mesh has more than $5\cdot 10^{4}$ dofs and $\varepsilon_{\mathrm{n}}^{el}>10^{6}\,E_{0}/L_{0}$, which is highly accurate. The relative contact force error of LR NURBS and homogeneous meshes is shown in Fig. 5.5. The proposed local refinement technique using LR NURBS elements show very good convergence behavior, both for quadratic and cubic elements.

It turns out that for cubic elements the obtained relative error decreases to $\approx 1.63\,\cdot 10^{-5}$, while for quadratic elements $\approx 2.1\,\cdot 10^{-4}$ is achieved. Since the penalty parameter is increased by Eq. (5.3) and since it is depending on the polynomial order, higher accuracy for cubic elements can be expected. The dofs are reduced by a factor of more than $16$ for LR NURBS meshes with quadratic elements at the $5^{\mathrm{th}}$ refinement depth. For LR NURBS meshes with cubic elements at the $4^{\mathrm{th}}$ refinement depth the dofs are reduced by a factor of more than $37$. This is a huge decrease of the computational cost. In Tab. 5.1 the mesh properties and numerical results for cubic elements are listed. This example shows that with the local refinement procedure, a high accuracy is achieved while decreasing the computational cost.

In a second part of this example, the performance of LR NURBS elements for different sizes of the contact domain is investigated. Computations for different radii of the rigid sphere are performed. The initial radius is set to $R=1/10\,L_{0}$ and is increased step-wise to $R=L_{0}$. The sphere is positioned at $X=Y=0$ and lowered gradually to $z=-R/2$. The performance of LR NURBS meshes using cubic elements is investigated. Local refinement of depth 5 is applied. The relative error of the normal contact force is defined w.r.t. uniformly discretized meshes using cubic NURBS elements. The smallest element size of each mesh is identical, leading to the same value for $\varepsilon_{\mathrm{n}}^{el}$. The percentage of the used dofs with respect to the radius $R$ is shown in Fig. 5.6. The relative contact force error for all radii $R$ is in the range of $5\cdot 10^{-8}$ and $3\cdot 10^{-7}$. The LR NURBS meshes capture the reference solutions nicely. A high accuracy is achieved while only using few dofs. Consequently, for $R=1/10L_{0}$ less than $2\%$ of the time is used in comparison to the reference model. The time for the local refinement procedure is negligible. This reduction of the computational effort that still achieves high accuracy demonstrates the benefit of the local refinement technique using LR NURBS.

5.3 Frictionless sliding contact

In the third example the performance of the adaptive local refinement and coarsening technique from Sec. 4 is investigated. For this, frictionless sliding contact is considered. The problem setup consists of a rectangular membrane sheet having dimension $8\,\lambda\,L_{0}\times 2\,\lambda\,L_{0}$ with isotropic pre-stretch $\lambda=1.25$. The boundaries are clamped. A rigid sphere with radius $L_{0}$ is initially located at $\boldsymbol{X}_{0}=\lambda\cdot[L_{0},L_{0},L_{0}]$ and first moved downwards to $\boldsymbol{x}_{0}=\lambda\cdot[L_{0},L_{0},L_{0}/2]$. This is followed by horizontal motion until the sphere reaches $\boldsymbol{x}_{0}=\lambda\cdot[7\,L_{0},L_{0},L_{0}/2]$, see left side of Fig. 5.7. The volume enclosed by the membrane is constrained to be constant so that it behaves like a cushion (Sauer et al.,, 2014). The material and the computational parameters are the same as in the previous example. During frictionless sliding the contact domain changes, and the adaptive local refinement and coarsening procedure from Sec. 4 is applied to obtain highly resolved meshes in the contact domain.

The right side of Fig. 5.7 shows the deformed configuration after the downward displacement of the rigid sphere has been applied. The parameters for the automatic control from Sec. 4 are $d_{\mathrm{ref}}^{d}=3\,d_{\mathrm{e}}^{1-d}$, $d_{\mathrm{safe}}^{d}=2\,d_{\mathrm{e}}^{1-d}$ and $d_{\mathrm{crs}}^{d}=4\,d_{\mathrm{e}}^{0}$. In Fig. 5.8 two meshes resulting from the adaptive local refinement and coarsening technique are illustrated. It can be observed that during sliding only the contact domain has a high aggregation of LR NURBS elements.

In the following, the performance of adaptive local refinement and coarsening is investigated. For this, LR NURBS meshes are compared to uniformly discretized meshes using standard NURBS. The normal and tangential contact forces for the LR NURBS, the uniformly refined and the initial, coarse mesh are illustrated in Fig. 5.9. The normal and tangential contact forces for the LR NURBS and the uniformly refined mesh match each other nicely. The contact forces of the initial, coarse mesh are oscillating. The oscillations arise from the coarse discretization. The normal contact force of the initial, coarse mesh is much smaller than the normal contact forces of the refined meshes. The reason is the mesh dependent penalty parameter, see Eq.(5.3). The net tangential contact force $f_{\mathrm{t}}$ is small in comparison to the net normal contact force $f_{\mathrm{n}}$, as the sliding is considered to be frictionless.

In Fig. 5.10 the relative contact force error of $f_{\mathrm{t}}$ and $f_{\mathrm{n}}$ is illustrated. It can be seen that the error is below $10^{-4}$ for quadratic LR NURBS and below $10^{-6}$ for cubic LR NURBS. The LR NURBS meshes capture the reference solution nicely with less than $30\%$ dofs in both examples. The high relative error of quadratic LR NURBS elements in comparison to cubic elements is caused by the coarse periphery. The evolution of the relative error in tangential direction shows abrupt changes. These result from the geometric approximation error of the coarsening steps. The setting of the parameters for the automatic control leads to an almost constant evolution of the relative normal contact force error. This example shows that the adaptive local refinement and coarsening technique using LR NURBS elements leads to highly accurate results while decreasing the computational cost.

5.4 Frictional contact of two deformable membranes

In the next example, the performance of LR NURBS is investigated for frictional contact considering two deformable membranes. This example addresses two new aspects. First, the preservation of sliding variables during coarsening and refinement. Second, the local refinement and coarsening of two objects. We consider two deformable rectangular solid membranes with dimension $2\,L_{0}^{u}\times 0.5\,L_{0}^{u}$ and $0.5\,L_{0}^{l}\times 2\,L_{0}^{l}$. They are fixed at their boundaries in all parametric directions. An isotropic pre-stretch of $\lambda=1.5$ in longitudinal direction is applied to avoid membrane instabilities (like wrinkling) during frictional sliding. The left side of Fig. 5.11 shows the problem setup. Both surfaces are inflated until they almost touch. This is followed by sliding the upper membrane along the lower one. During sliding, the inflated volume of each membrane is constrained to be constant. The adaptive locally refined surfaces at initial contact detection are illustrated on the right side of Fig. 5.11. The friction coefficient is set to $\mu_{\mathrm{f}}=0.25$ and the enclosed volume of both surfaces is constrained to be constant after inflation. The material behavior is the same as in the previous examples. The number of Gaussian quadrature points is $5\times 5$ for each element. Contact is treated by the penalty method. The penalty parameter is increased by Eq. (5.3) with $\varepsilon_{0}=80\,E_{0}/L_{0}$ and it is taken equal in tangential and normal direction. The two-half-pass algorithm (Sauer and De Lorenzis,, 2015) is used to evaluate the contact forces. Due to frictional sliding the contact domain changes and the presented adaptive local refinement and coarsening procedure is applied on both surfaces. In Fig. 5.12 two meshes resulting from the adaptive local refinement and coarsening technique are illustrated.

The performance of quadratic and cubic LR NURBS is investigated. The LR NURBS meshes are compared to their uniform counterparts using quadratic and cubic NURBS discretizations separately. The normal contact force $f_{\mathrm{n}}$ and the tangential contact force $f_{\mathrm{t}}$ acting on the lower surface are shown in Fig. 5.13. Their behavior is as expected, and both discretizations match each other nicely.

Next, the relative contact force errors

[TABLE]

are investigated. They are illustrated in Fig. 5.14. The evolution of the relative error in normal direction shows a similar behavior as in the previous example. The cubic LR NURBS elements perform better than quadratic elements due to the coarse periphery. The evolution of the relative error in tangential direction is more smooth than in the previous example, see Sec. 5.3. The reason is that the geometric approximation error shows no remarkable influence in this example.

A maximum relative error of $\approx 3\cdot 10^{-4}$ for quadratic LR NURBS discretizations is achieved. The relative error for cubic LR NURBS elements is much smaller and a maximum of $\approx 3\cdot 10^{-6}$ is achieved. The parameters for the automatic control from Sec. 4 are $d_{\mathrm{ref}}^{d}=3\,d_{\mathrm{e}}^{1-d}$, $d_{\mathrm{safe}}^{d}=2\,d_{\mathrm{e}}^{1-d}$ and $d_{\mathrm{crs}}^{d}=4\,d_{\mathrm{e}}^{0}$. With this setting less than $30\%$ of the dofs of the reference meshes are used for both quadratic and cubic LR meshes. This example shows that the technique of adaptive local refinement and coarsening can be successfully applied to frictional contact of membranes. The computational cost is reduced while still achieving high accuracy. In the next example the adaptive local refinement and coarsening technique is applied to LR NURBS-enriched volume elements considering frictional contact.

5.5 Frictional ironing contact between two deformable solids

Local refinement for general 3D isogeometrical models is challenging and the linear independency of LR meshes in 3D has not been developed yet. To overcome this issue, the NURBS-enrichment technique by Corbett and Sauer, (2014) is used to extend LR NURBS surface discretizations to 3D.

5.5.1 LR NURBS-enriched finite elements

The core idea of this technique is to use standard finite Lagrange elements in the bulk domain and enrich the surface by isogeometric finite elements. The local refinement of NURBS-enriched elements is performed in two main steps

Local refinement of the isogeometric surface mesh using LR NURBS. The adpative local refinement technique from Sec. 4 is applied. 2. 2.

Extend the local refinement through all the Lagrange elements in the thickness direction

This work considers the local refinement of NURBS-enriched and Lagrange elements in two dimensions only. The refinement of the thickness direction is not considered here, but it is possible to use standard refinement procedures. The local refinement of the Lagrange elements typically involve the generation of hanging nodes. For the treatment of these hanging nodes the approach of Demkowicz et al., (1989) is used.

5.5.2 Numerical results

The last example considers a deformable hollow hemisphere in frictional contact with a deformable block. Both solids have a quadratic LR NURBS-enriched contact surface. The hemisphere is discretized with three elements in its thickness direction and $7\times 7$ elements on the contact surface. The block is discretized with eight elements in its thickness direction and $4\times 40$ elements on the contact surface. Each contact element has $5\times 5$ Gaussian quadrature points. The bulk domain is discretized by linear Lagrange finite elements. Initially, the block has dimension $L_{0}\times L_{0}\times 10\,L_{0}$ and the hollow hemisphere has the inner radius $R_{0}^{\mathrm{in}}=L_{0}/6$ and the outer radius $R_{0}^{\mathrm{out}}=L_{0}/2$, see Fig. 5.15. The base of the block is fixed in all directions. Periodic boundary conditions are used on the left and right boundary of the block. First, a downward displacement of $2/3R_{0}^{\mathrm{out}}$ is applied to the top of the hollow hemisphere, which is followed by a displacement of $6\,L_{0}$ in tangential direction. The downward displacement is considered to be frictionless while the tangential displacement is considered to be frictional. The solids follow the isotropic, non-linear Neo-Hookean material model (e.g. see Zienkiewicz et al., (2005))

[TABLE]

with the deformation gradient $\boldsymbol{F}$ and the identity tensor $\boldsymbol{I}$. The Lamé constants $\Lambda$ and $\mu$ can be expressed in terms of the Poisson’s ratio $\nu$ and Young’s modulus $E$ by

[TABLE]

Considered here are $\nu_{1}=\nu_{2}=0.3$ for both solids and $E_{1}=E_{0}$ for the block and $E_{2}=10\,E_{0}$ for the hollow hemisphere. Contact is computed by the two-half-pass algorithm using the penalty method with $\varepsilon_{0}=80\,E_{0}/L_{0}$ and increasing $\varepsilon^{el}_{\mathrm{n}}=\varepsilon^{el}_{\mathrm{t}}$ by Eq. (5.3). The friction coefficient is set to $\mu_{\mathrm{f}}=0.25$.

Next, the performance of adaptive local refinement and coarsening using quadratic LR NURBS-enriched finite elements is investigated. For this, LR NURBS meshes are compared to uniformly discretized meshes using quadratic NURBS-enriched finite elements. Local refinement of depth 2 is applied. Fig. 5.16 shows the locally refined and deformed block after the hollow hemisphere has been moved downward. The coloring shows the normalized stress invariant I_{1}=\mathrm{tr}\,\mbox{\boldmath\sigma}/E_{0}. It can be observed that the area of major interest is locally refined while the periphery remains coarse. As the surface of the hemisphere is almost entirely in contact with the block the local refinement leads to a uniformly refined hemisphere, see Fig. 5.17.

The normal and tangential contact forces for LR NURBS and uniform meshes are shown in Fig. 5.18. These are evaluated at the block’s contact surface. The contact forces for the two meshes only differ slightly. By first investigating the normal contact force an oscillatory and periodic behavior can be observed. The oscillations can be reduced by further mesh refinement. The normal contact force of the LR NURBS mesh shows abrupt increases at several points. This is caused by an adaptive local refinement step. The same can be observed for the tangential contact force, see right side of Fig. 5.18. The absolute and relative contact force errors in normal and tangential direction are illustrated in Fig. 5.19. The relative contact force errors are approximately in the range of $4\cdot 10^{-4}$ to $2\cdot 10^{-3}$. The errors show a continuous increase followed by an abrupt decrease. The errors reflect the influence of the refinement and coarsening technique also seen in Fig. 5.18. Comparing the relative contact force error to the results of the previous example, one order of magnitude is lost (in case of quadratic LR NURBS elements). The third dimension of the solids affect the numerical results strongly. For further improvements one could take local error measures into account. It is expected that by using cubic LR NURBS-enrichment the error would decrease but this would not reduce the influence of the solid domain. The parameters for the automatic control from Sec. 4 are $d_{\mathrm{ref}}^{d}=5\,d_{\mathrm{e}}^{1-d}$, $d_{\mathrm{safe}}^{d}=4\,d_{\mathrm{e}}^{1-d}$ and $d_{\mathrm{crs}}^{d}=6\,d_{\mathrm{e}}^{0}$. With this setting the LR NURBS mesh has still $\approx 40\%$ of the dofs of the reference mesh. Reducing the size of the locally refined domain would lead to an increase of the error of the contact forces. Still, the adaptive local refinement and coarsening of LR NURBS-enriched elements show the benefits of this procedure.

6 Conclusion

This work presents a novel concept of adaptive local surface refinement using LR NURBS elements in the framework of computational contact mechanics. The Bézier extraction of LR NURBS elements is obtained by an additional mapping of the Bézier extraction operator. With this a convenient embedding of LR NURBS elements into general finite element codes is achieved.

The adaptive local refinement and coarsening technique is automatically controlled by the presented refinement indicators. The numerical results show that the geometric approximation error arising from the coarsening is negligible. The technique is applied to pure isogeometric surface elements formulation and to isogeometrically-enriched volume elements. The numerical examples show the benefit of using LR NURBS. Due to the local refinement, the computational cost decreases, while still achieving high accuracy. A good convergence behavior for LR NURBS elements is also achieved.

The numerical examples in Sec. 5.1 and Sec. 5.5 consider curved surfaces that truly use NURBS instead of B-splines. So they show that the extension from LR B-splines to LR NURBS is implemented successfully. The automatic control of the adaptive local refinement and coarsening technique works robustly for frictionless and frictional contact. As seen in the examples, quadratic LR NURBS achieve less accuracy than cubic LR NURBS. The reason for this is that the latter capture the solution in the coarse (i.e. unrefined) mesh more accurately. This motivates the incorporation of local error measures that indicate elements for refinement. Especially, the computations performed with quadratic LR NURBS-enriched elements show the demand for such an error measure. The refinement depth is prescribed in the presented examples. By using error measures the refinement depth can be specified automatically to achieve a desired accuracy. This work considers the local refinement of LR NURBS-enriched and Lagrange elements in two dimensions. More complex 3D objects call for an enhanced adaptive local refinement technique, which includes local refinement in the third dimension.

In this work, only bivariate LR NURBS are considered. Their extension to trivariate LR NURBS would be interesting, since they allow local refinement of 3D LR NURBS elements. This is a topic of its own because of the complex structure of trivariate LR NURBS meshes. Ensuring that the resulting LR mesh is linearly independent has not been proven yet and is challenging for arbitrary meshline extensions.

Acknowledgements

The authors are grateful to the German Research Foundation (DFG) for supporting this research through project GSC 111.