Numerical modeling of contact problems with the finite element method utilizing a B-Spline surface for contact surface smoothing

Abstract

The basic aim of this work is to present a method to treat structural mechanics problems dealing with tridimensional contact and large elastic deformations. The formulation presented in this article offers a B-Spline based discretization for the contact surface, which solves the continuity problems presented by the classic Lagrangian contact element formulation. A B-Spline surface is utilized to discretize the contact surface, and the B-Spline basis functions are used to distribute the contact forces between the contact surfaces nodes. The finite element utilized is an 8-node hexahedron, and the formulation is continuum mechanics based, written in the current configuration, utilizing a neo-Hookean material model. The contact restrictions are enforced by the augmented Lagrangian method.

Key-words
Contact mechanics; B-Spline; Surface smoothing; Augmented Lagrangian method; Friction

1 INTRODUCTION

Finite element contact simulations exhibiting large deformations has been growing as a subject of study in the last few decades. Many works have been produced about it, and one can cite WRIGGERS in 1989 and 1996, BENSON (1990) and LAURSEN AND SIMO in 1993 as authors giving important contributions to the field, as well as CURNIER in 1999, and BANDEIRA in 2001, who offered further research on the subject.

The classical way to discretize a contact surface on a tridimensional problem is through the Lagrangian contact element. The Lagrangian contact element for a hexahedron finite element mesh is a 4-node, linear element, which results in a flat surface. When discretizing a curved surface, the Lagrangian contact element, being a linear element, leads to a facetization of the contact surface, i. e., the curved surface is represented by several flat surfaces, side by side, with an angle between them. This facetization leads to a geometric discontinuity of the contact surface, and therefore, generates discontinuity of the normal vector. This in turn, leads to convergence problems, especially when large sliding occurs.

To tackle this issue, many approaches have been proposed. One may consider, between many options, the Mortar method, described in (PUSO, LAURSEN and SOLBERG, 2008), the Hermitian element (WRIGGERS, 2006) and (ALAIN BATAILLY, 2013), the smoothing of the contact surface using B-Splines or NURBS surfaces (CALLUM and CORBETT, 2014), a combination of NURBS and Mortar presented by (TEMIZER, 2011), or the isogeometric analysis as proposed by (M.E. MATZEN, 2012) and (DE LORENZIS, WRIGGERS, ZAVARISE, 2012), which does away with the classical Lagrangian element, and utilizes a new contact element based on B-Spline or NURBS basis function.

This paper describes an isogeometric approach, where a B-Spline patch is used to discretize a contact surface over the master body. The difference between this paper and other works which utilize the isogeometric approach is that no smoothing of the slave surface is done. The smoothing is done only on the master surfaces. Each node of the slave surface is tested individually, and the slave surface conforms to the master surface. The authors believe that when contact happens, the slave surface will assume the master surface’s shape, and thus will also be smoothed.

In the first section of this article, a brief description of the contact problem, the virtual work applied to the balance of linear momentum equation and the Neo-Hookean material constitutive equation are presented. Then, the equation for contact forces is shown. Afterwards, the B-Spline concept is briefly described, and the method for discretizing the contact surface using the B-Spline is detailed. The relevant equations for the contact pressures, their linearization and finally the relevant vectors for the discretization of the contact surfaces are presented.

In the last section numerical examples are shown, utilizing the formulation presented in this paper. The results obtained are compared to the commercial software ANSYS and afterwards the results are briefly discussed.

2 PROBLEM DESCRIPTION

It is shown, in Figure 1, a tridimensional contact problem between two deformable solids as described in (SIMO, LAURSEN, 1992). The group of material points, that belong to the bodies in the reference position, are called ${\text{B}}_1$ and ${\text{B}}_2$. Additional configurations can be obtained by the transformations ${\text{φ}}^1:{\text{B}}^1\times \text{I}\to {\text{R}}^3$ and ${\text{φ}}^2:{\text{B}}^2\times \text{I}\to {\text{R}}^3$ in the interval t = [0, T]. The bodies ${\text{B}}_1$ and ${\text{B}}_2$, subject to contact, are called slave and master. In each of these bodies, there is a boundary where contact can be expected. These boundaries are called ${\text{Г}}^1$ and ${\text{Г}}^2$. ${\text{Г}}^1$ is the slave boundary, belonging to the body ${\text{B}}_1$, and ${\text{Г}}^2$ the master boundary, belonging to the body ${\text{B}}_2$. These boundaries are in the reference configuration, and are subject to the same transformations ${\text{φ}}^1$ and ${\text{φ}}^2$ that the bodies ${\text{B}}_1$ and ${\text{B}}_2$ are subject. A point belonging to body ${\text{B}}_1$ is denoted by ${\text{X}}^1$ in the reference configuration and ${\text{x}}^1$ on the current configuration, and the same applies for body ${\text{B}}_2$, exchanging the index.

Figure 1
Notation for the contact problema with friction

The displacement u is defined as,

$u^{\text{α}}=x^{\text{α}}-X^{\text{α}}=\text{φ}\left(X^{\text{α}},\text{t}\right)-X^{\text{α}},$ (2.1)

with $\text{α}$ = 1,2.

To solve the contact problem, it is necessary to calculate the gap between the bodies ${\text{B}}_1$ and ${\text{B}}_2$ at the points where contact might occur. This is done by searching in ${\text{Г}}^2$ for a point ${\overline{\text{x}}}^2$, closest to a point ${\text{x}}^1$ in ${\text{Г}}^1$, by the means of the function

$\left|\text{g}\left(X,\text{t}\right)\right|={\text{min}}_{\text{Y}\in {\text{Г}}^2}\left|\left|{\text{φ}}^1\left(X^1,\text{t}\right)-{\text{φ}}^2\left(X^2,\text{t}\right)\right|\right|.$ (2.2)

The minimum for this function can be calculated by finding the roots of the first derivate.

3 VIRTUAL WORK PRINCIPLE

The equilibrium equation, derived from the balance of linear momentum is defined by,

$\text{Div}{P}_{\left(\text{uh}\right)}+{\text{ρ}}_0\overline{b}-{\text{ρ}}_0\frac{\text{d}v}{\text{dt}}=\mathrm{0,}$ (3.1)

where $P$ is the first Piola Kirchoff tensor, ${\text{ρ}}_0$ is the body’s density and $\overline{b}$ the body forces for unit of volume.

After applying virtual displacements, and bringing the tensors from the reference configuration to the current configuration, one obtains

$\text{w}\left(\text{φ},\text{η}\right)=\underset{\text{φ}\left(\text{B}\right)}{\overset{}{{\displaystyle \int }}}\sigma \cdot {\nabla }^{\text{S}}\eta {\text{dV}}_{\text{t}}-\underset{\text{φ}\left(\text{B}\right)}{\overset{}{{\displaystyle \int }}}{\text{ρ}}_{\text{t}}\left(b-\frac{\text{d}v}{\text{dt}}\right)\cdot \eta {\text{dV}}_{\text{t}}-\underset{\text{φ}\left(\partial \text{Bσ}\right)}{\overset{}{{\displaystyle \int }}}t\cdot \eta {\text{dA}}_{\text{t}}=0.$ (3.2)

In this equation, $\sigma$ is Cauchy’s tensor and ${\text{ρ}}_{\text{t}}$ is density for current unit of volume. This equation represents equilibrium between the internal work and the external work defined at the current configuration.

It’s necessary to include the material behavior in the finite element simulation, thus an equation for a material model is necessary. The Neo-Hookean material model is adopted on this article. This material model is adequate for simulating large tridimensional displacements. Therefore, Cauchy’s tensor on the equation (3.2) is obtained through the deformation energy function for a Neo-Hookean material, given by

$\sigma =\frac{\text{Λ}}{2\text{J}}\left({\text{J}}^2-1\right)I+\frac{{\text{μ}}_{\text{d}}}{\text{J}}\left(b-I\right).$ (3.3)

Equation (3.2) can be divided into two parts, internal and external forces. The first term represents the internal forces, and the second and third terms represent the external forces,

${\text{П}}_{\text{Int}}=\underset{\text{φ}\left(\text{B}\right)}{\overset{}{{\displaystyle \int }}}\sigma \cdot {\nabla }^{\text{S}}\eta {\text{dV}}_{\text{t}}$ (3.4)

and

${\text{П}}_{\text{Ext}}=\underset{\text{φ}\left(\text{B}\right)}{\overset{}{{\displaystyle \int }}}{\text{ρ}}_{\text{t}}\left(b-\frac{\text{d}v}{\text{dt}}\right)\cdot \eta {\text{dV}}_{\text{t}}-\underset{\text{φ}\left(\partial \text{Bσ}\right)}{\overset{}{{\displaystyle \int }}}t\cdot \eta {\text{dA}}_{\text{t}}.$ (3.5)

Moreover, since the problem being simulated is a contact problem, the equations need to depict the behavior of both bodies, , as exhibited on (LAURSEN, 2002). Equation (3.2) is then expanded to include both bodies,

$\text{W}\left(\text{φ},\text{η}\right)={\displaystyle \sum }_{\text{α}=1}^2\left[\underset{{\text{φ}}^{\text{α}}{\left(\text{B}\right)}^{\text{α}}}{\overset{}{{\displaystyle \int }}}{\left(\sigma \cdot {\nabla }^{\text{S}}\eta {\text{dV}}_{\text{t}}\right)}^{\alpha }-\underset{{\text{φ}}^{\text{α}}{\left(\text{B}\right)}^{\text{α}}}{\overset{}{{\displaystyle \int }}}{\text{ρ}}_{\text{t}}\left(b-\frac{\text{d}v}{\text{dt}}\right)\cdot \eta {\text{dV}}_{\text{t}}{}^{\text{α}}-\underset{{\text{φ}}^{\text{α}}{\left(\partial \text{Bσ}\right)}^{\text{α}}}{\overset{}{{\displaystyle \int }}}t\cdot \eta {\text{dA}}_{\text{t}}{}^{\text{α}}\right]=0.$ (3.6)

To satisfy the impenetrability condition and to represent the contact forces, a contact energy term is added to equation (3.6).

${\text{П}}_{\text{Int}}+{\text{П}}_{\text{Ext}}+{\text{П}}_{\text{C}}=0$ (3.7)

The definition of ${\text{П}}_{\text{C}}$ depends upon the restriction method utilized. In this paper, the augmented Lagrangian method was utilized, as described in (BANDEIRA, et al., 2010), (BERTSEKAS, 1995) and (LUENBERGER, 1984). The contact equation is defined, as shown in (BANDEIRA, 2001), (WRIGGERS, 2008) and (LAURSEN and SIMO, 1993a), by

${\text{П}}_{\text{C}}=\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}{\text{t}}_{\text{N}}{\text{δg}}_{\text{N}}\text{dA}+\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}{t}_{\text{T}}\text{δ}{g}_{\text{T}}\text{dA}.$ (3.8)

The terms ${\text{t}}_{\text{N}}$ and $t_{\text{T}}$ represent the normal and tangencial contact pressures, respectively. ${\text{g}}_{\text{N}}$ represents the gap function, which notates the penetration between the bodies, and $g_{\text{T}}$ is the tangencial displacement function, which represents how much one of the bodies slides over the other.

The first integral on equation (2.8) accounts for the work done by normal forces, while the second integral accounts for the work done by the tangential forces. If the problem doesn’t include friction, the second integral is null.

Equation (3.7) then can be discretized and solved by the finite element method. In this paper, the isoparametric hexahedral element was used. More details on the finite element discretization can be found on (BATHE, 1996), (ZIENKIEWICZ, TAYLOR, ZHU, 2005), and (COOK, 2002).

4 NODE TO SURFACE FORMULATION

In this section, the node-to-surface contact formulation is expressed, for more details, refer to (WRIGGERS, 2008) and (BANDEIRA, 2001). The contact formulation is independent of the finite element utilized for its discretization, hence it can be used for a B-Spline based finite element.

The equation for the contact forces as found in BANDEIRA is defined in (3.8) and is repeated here for convenience,

${\text{C}}_{\text{c}}=\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}{\text{t}}_{\text{N}}{\text{δg}}_{\text{N}}\text{dA}+\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}{t}_{\text{T}}\text{δ}{g}_{\text{T}}\text{dA}.$ (4.1)

The first term accounts for the normal forces and the second one to the tangential forces. When the augmented Lagrangian method is used to enforce the contact restriction, the term ${\text{t}}_{\text{N}}$ is defined by

${\text{t}}_{\text{N}}=\langle {\text{λ}}_{\text{N}}+{\text{ξ}}_{\text{N}}{\text{g}}_{\text{N}}\rangle ,$ (4.2)

where ${\text{λ}}_{\text{N}}$ is the Lagrangian multiplier, ${\text{ξ}}_{\text{N}}$ is the penalty factor and ${\text{g}}_{\text{N}}$ is the gap function. In turn, ${\text{g}}_{\text{N}}$ which is defined by

${\text{g}}_{\text{N}}=\{\begin{array}{c}\left(x^{\text{s}}-{\overline{x}}^{\text{m}}\right)\cdot {\overline{n}}_{\text{c}}if\left(x^{\text{s}}-{\overline{x}}^{\text{m}}\right)\le 0\\ 0if\left(x^{\text{s}}-{\overline{x}}^{\text{m}}\right)>0\end{array},$ (4.3)

where $x^{\text{s}}$ is the point of the slave surface currently analyzed, or rather, a node, and ${\overline{x}}^{\text{m}}$ is the projection of the slave point on the master surface.

The term ${\text{δg}}_{\text{N}}$ is given by

${\text{δg}}_{\text{N}}=(v^{\text{s}}-{\overline{v}}^{\text{m}})\cdot {\overline{n}}_{\text{c}}.$ (4.4)

The second term of equation (4.1) accounts for the friction forces, parallel to the contact surface. On that equation, $t_{\text{T}}$ depends on the friction model utilized. In this paper, the Coulomb friction model was deemed sufficient. In the Coulomb friction model, the term $t_{\text{T}}$ depends on whether the friction is on the stick case, where the two bodies don’t move relative to each other, or the slip case, where one body slides over the other. Depending on which is the current case, $t_{\text{T}}$ is calculated differently. The following functions are used to calculate it.

${\text{t}}_{{\text{T}}_{\text{α}}}=t_{\text{T}}\cdot {\overline{a}}_{\text{α}}$ (4.5)

and

${\text{t}}_{{\text{Tα}}_{\text{n}+1}}=\{\begin{array}{c}{\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}se{\text{Φ}}_{\text{n}+1}^{\text{trial}}\le 0\\ \mu t_{{\text{N}}_{\text{n}+1}}\frac{{\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}}{{\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}}se{\text{Φ}}_{\text{n}+1}^{\text{trial}}>0\end{array}$ (4.6)

The value of ${\text{Φ}}_{\text{n}+1}^{\text{trial}}$ defines whether the current node is on a stick or slip case, and then ${\text{t}}_{{\text{Tα}}_{\text{n}+1}}$ is defined accordingly. ${\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}$ is defined by

${\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}={\text{t}}_{{\text{Tα}}_{\text{n}}}^{}+\Delta {\text{λ}}_{\text{Tα}}+{\text{ξ}}_{\text{T}}{\text{u}}_{\text{Tα}},$ (4.7)

and ${\text{Φ}}_{\text{n}+1}^{\text{trial}}$ is given by

${\varphi }_{\text{n}+1_{}}^{\text{trial}}={\left[{\text{t}}_{{\text{Tα}}_{\text{n}+1}}^{\text{trial}}{\text{m}}_{\text{αβ}}^{\left(\text{t}\right)}{\text{t}}_{{\text{T}}_{\text{n}+1}}^{{\text{β}}^{\text{trial}}}\right]}^{1/2}-{\text{μ t}}_{{\text{N}}_{\text{n}+1}}.$ (4.8)

For more details regarding the node-to-surface formulation utilizing the augmented Lagrangian method to enforce contact restrictions, one can refer to (BANDEIRA, 2001).

Going back to (4.1), $\text{δ}{g}_{\text{T}}$ is defined as

$\mathrm{\delta }{\mathbf{g}}_{\mathrm{T}}={\mathbf{v}}^{\mathrm{s}}-{\stackrel{-}{\mathbf{v}}}^{\mathrm{m}}={\dot{\stackrel{-}{\mathrm{\xi }}}}^{\mathrm{\alpha }}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }},$ (4.9)

with ${\dot{\overline{\text{ξ}}}}^{\text{α}}$ defined by

${\dot{\overline{\text{ξ}}}}^{\text{α}}=\frac{1}{{\overline{\text{a}}}_{\text{αβ}}-{\text{g}}_{\text{N}}{\overline{\text{b}}}_{\text{αβ}}}\left\{\left[v^{\text{s}}-{\overline{v}}^{\text{m}}\right]\cdot {\overline{a}}_{\text{α}}+{\text{g}}_{\text{N}}{\overline{n}}_{\text{c}}{\overline{v}}_{,\text{α}}^{\text{m}}\right\}$ (4.10)

and the following terms defined as follows,

${\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}={\mathbf{x}}_{,\mathrm{\xi }}^{\mathrm{m}}\left(\stackrel{-}{\mathrm{\xi }}\right)\bullet {\mathbf{x}}_{,\mathrm{\xi }}^{\mathrm{m}}\left(\stackrel{-}{\mathrm{\xi }}\right)={\stackrel{-}{\mathrm{a}}}_{\mathrm{\beta }}\bullet {\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }},$ (4.11)

${\overline{\text{b}}}_{\text{αβ}}=x_{,\text{ξξ}}^{\text{m}}\left(\overline{\text{ξ}}\right)\cdot {\overline{n}}_{\text{c}}={\widehat{x}}_{,\text{αβ}}^{\text{m}}\cdot {\overline{n}}_{\text{c}},$ (4.12)

and

${\overline{a}}_{\text{α}}=\frac{\partial x^{\text{m}}\left({\text{ξ}}^{\text{α}},\text{t}\right)}{\partial {\text{ξ}}^{\text{α}}}=x_{,\text{α}}^{\text{m}}\left({\text{ξ}}^{\text{α}},\text{t}\right).$ (4.13)

4.1 Linearization

The solution for the equation (4.1) is discovered using Newton-Raphson’s method, and as such, the derivatives of equation (4.1) are necessary. The derivative is given by

$\text{K}=\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}\left[\Delta {\text{t}}_{\text{N}}{\text{δg}}_{\text{N}}+{\text{t}}_{\text{N}}\Delta \left({\text{δg}}_{\text{N}}\right)\right]\text{dA}+\underset{{\text{dB}}_{\text{OC}}^{\text{S}}}{\overset{}{{\displaystyle \int }}}\left[\Delta {\text{t}}_{{\text{T}}_{\text{α}}}\text{δ}{\overline{\text{ξ}}}^{\text{α}}+{\text{t}}_{{\text{T}}_{\text{α}}}\Delta \left(\text{δ}{\overline{\text{ξ}}}^{\text{α}}\right)\right]\text{dA}.$ (4.14)

For the first integral, the term $\Delta {\text{t}}_{\text{N}}$ is given by

$\Delta {\text{t}}_{\text{N}}={\text{ε}}_{\text{N}}\Delta {\text{g}}_{\text{N}},$ (4.15)

with $\Delta {\text{g}}_{\text{N}}$ defined as

$∆{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}^{}=\left\{\begin{array}{c}\left(∆{\mathrm{u}}^{\mathrm{s}}-∆{\stackrel{-}{\mathrm{u}}}^{\mathrm{m}}\right)\cdot {\stackrel{-}{\mathrm{n}}}_{\mathrm{c}},se{\mathrm{g}}_{\mathrm{N}}<0\\ 0,se{\mathrm{g}}_{\mathrm{N}}>0\end{array},\right.$ (4.16)

and ${\text{t}}_{\text{N}}$ and ${\text{δg}}_{\text{N}}$ as shown previously. The term $\Delta \left({\text{δg}}_{\text{N}}\right)$ is given by

$∆\left(\mathrm{\delta }{\mathrm{g}}_{\mathrm{N}}\right)=-{\stackrel{-}{\mathbf{v}}}_{,\mathrm{\alpha }}^{\mathrm{m}}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\alpha }}\cdot {\stackrel{-}{\mathbf{n}}}_{\mathrm{C}}-\left({\mathbf{v}}^{\mathrm{s}}-{\stackrel{-}{\mathbf{v}}}^{\mathrm{m}}\right)\cdot {\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}{\stackrel{-}{\mathrm{a}}}^{\mathrm{\alpha }\mathrm{\beta }}{\stackrel{-}{\mathbf{n}}}_{\mathrm{C}}\cdot \left(∆{\mathbf{u}}_{,\mathrm{\beta }}^{\mathrm{m}}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\beta },\mathrm{\gamma }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\gamma }}\right),$ (4.17)

where $\Delta {\overline{\text{ξ}}}^{\text{α}}$ is defined by

$∆{\mathrm{\xi }}^{\mathrm{\beta }}=\frac{1}{{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}-{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\stackrel{-}{\mathrm{b}}}_{\mathrm{\alpha }\mathrm{\beta }}}\left\{\left[{∆\mathbf{u}}^{\mathrm{s}}-{∆\stackrel{-}{\mathbf{u}}}^{\mathrm{m}}\right]\cdot {\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}+{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\mathrm{}\stackrel{-}{\mathbf{n}}}_{\mathrm{C}}\cdot {\stackrel{-}{\mathbf{u}}}_{,\mathrm{\alpha }}^{\mathrm{m}}\right\}.$ (4.18)

In the second term of equation (4.14), $\Delta {\text{t}}_{{\text{T}}_{\text{α}}}$ is defined as

$\Delta {\text{t}}_{{\text{T}}_{\text{α}}}=\Delta {\text{t}}_{{\text{T}}_{\text{α}}}^{\text{trial}}={ϵ}_{\text{T}}\Delta {\overline{\text{ξ}}}^{\text{β}}{\overline{\text{a}}}_{\text{αβ}}+{ϵ}_{\text{T}}\left({\overline{\text{ξ}}}_{\text{i}+1}^{\text{β}}-{\overline{\text{ξ}}}_{\text{i}}^{\text{β}}\right)\Delta {\overline{\text{a}}}_{\text{αβ}}$ (4.19)

for the stick case, with $\Delta {\overline{\text{a}}}_{\text{αβ}}$ given by

$\Delta {\overline{\text{a}}}_{\text{αβ}}=\Delta {\overline{a}}_{\text{α}}\cdot {\overline{a}}_{\text{β}}+{\overline{a}}_{\text{α}}\cdot \Delta {\overline{a}}_{\text{β}}.$ (4.20)

For the slip case, $\Delta {\text{t}}_{{\text{T}}_{\text{α}}}$ is defined as

$\Delta {\text{t}}_{{\text{T}}_{\text{α}}}=\text{μ}{ϵ}_{\text{N}}\Delta {\text{g}}_{\text{N}}\frac{{\text{t}}_{{\text{T}}_{\text{α}}}^{\text{trial}}}{t_{\text{T}}^{\text{trial}}}+\text{μ}\frac{{\text{t}}_{\text{N}}^{}}{t_{\text{T}}^{\text{trial}}}\Delta {\text{t}}_{{\text{T}}_{\text{α}}}^{\text{trial}}+\text{μ}\frac{{\text{t}}_{\text{N}}^{}}{t_{\text{T}}^{\text{trial}}{}^3}{\text{t}}_{{\text{T}}_{\text{α}}}^{\text{trial}}{\text{t}}_{{\text{T}}_{}}^{\text{trial β}}(t_{\text{T}}^{\text{trial}}\Delta a_{\text{β}}-\Delta {\text{t}}_{{\text{T}}_{\text{β}})}^{\text{trial}},$ (4.21)

with $\Delta {\text{t}}_{{\text{T}}_{\text{α}}}^{\text{trial}}$ being given by

$∆{\mathrm{t}}_{{\mathrm{T}}_{\mathrm{\alpha }}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}={\mathrm{ϵ}}_{\mathrm{T}}\mathrm{}{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}\mathrm{}∆{\mathrm{\xi }}^{\mathrm{\beta }}+{\mathrm{ϵ}}_{\mathrm{T}}\left({\stackrel{-}{\mathrm{\xi }}}_{\mathrm{i}+1}^{\mathrm{\beta }}-{\stackrel{-}{\mathrm{\xi }}}_{\mathrm{i}}^{\mathrm{\beta }}\right)\left\{(∆{\mathbf{u}}_{,\mathrm{\alpha }}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\theta }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\theta }}\right)\bullet {\stackrel{-}{\mathbf{a}}}_{\mathrm{\beta }}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\bullet (∆{\mathbf{u}}_{,\mathrm{\beta }}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\beta },\mathrm{\theta }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\theta }})\}.$ (4.22)

Following that, $\text{δ}{\overline{\text{ξ}}}^{\text{α}}$ is defined as

$\mathrm{\delta }{\mathrm{\xi }}^{\mathrm{\beta }}=\frac{1}{{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}-{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\stackrel{-}{\mathrm{b}}}_{\mathrm{\alpha }\mathrm{\beta }}}\left\{\left[{\mathbf{v}}^{\mathrm{s}}-{\stackrel{-}{\mathbf{v}}}^{\mathrm{m}}\right]\cdot {\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}+{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\mathrm{}\stackrel{-}{\mathbf{n}}}_{\mathrm{C}}\cdot {\stackrel{-}{\mathbf{v}}}_{,\mathrm{\alpha }}^{\mathrm{m}}\right\},$ (4.23)

and $\Delta \left(\text{δ}{\overline{\text{ξ}}}^{\text{α}}\right)$ given by

$∆\left(\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}\right)=\frac{1}{{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}-{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\mathrm{b}}_{\mathrm{\alpha }\mathrm{\beta }}}\{-{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\cdot \left({\stackrel{-}{\mathrm{v}}}_{,\mathrm{\beta }}^{\mathrm{m}}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}+∆{\stackrel{-}{\mathrm{u}}}_{,\mathrm{\beta }}^{\mathrm{m}}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\beta },\mathrm{\gamma }}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\gamma }}\right)+\left(∆{\stackrel{-}{\mathrm{u}}}_{,\mathrm{\alpha }}^{\mathrm{m}}+{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}\right)\cdot \left({\mathbf{v}}^{\mathrm{s}}-{\stackrel{-}{\mathbf{v}}}^{\mathrm{m}}-{\stackrel{-}{\mathbf{a}}}_{\mathrm{\gamma }}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\gamma }}\right)+\left({\stackrel{-}{\mathrm{v}}}_{,\mathrm{\alpha }}^{\mathrm{m}}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}^{}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}\right)\cdot \left(∆{\mathbf{u}}^{\mathrm{s}}-∆{\stackrel{-}{\mathbf{u}}}^{\mathrm{m}}-{\stackrel{-}{\mathbf{a}}}_{\mathrm{\gamma }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\gamma }}\right)+{\mathrm{g}}_{\mathrm{N}}{\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\mathrm{}\cdot \left({\stackrel{-}{\mathrm{v}}}_{,\mathrm{\alpha }\mathrm{\beta }}^{\mathrm{m}}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}+∆{\stackrel{-}{\mathrm{u}}}_{,\mathrm{\alpha }\mathrm{\beta }}^{\mathrm{m}}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}+{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }\mathrm{\gamma }}\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\beta }}\mathrm{}∆{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\gamma }}\right).$ (4.24)

5 CONTACT ELEMENT DISCRETIZATION

The solids bodies involved in the contact problem must be discretized into a finite element mesh. In this paper, the bodies are discretized by 8-node hexahedron elements, and the contact surface is discretized using a B-Spline based finite element.

Using a B-Spline based finite element for the discretization of the contact surfaces gives a few advantages over the Lagrangian 4-node contact element, such as a smooth surface and a higher degree of continuity between elements, avoiding the numerical problems associated with the node-surface formulation. The definition for B-Spline surfaces is briefly shown in the next session, and it’s followed by the procedure to discretize the contact surface using a B-Spline based finite element.

5.1 B-Spline surface

A B-Spline surface is constructed based on the interpolation of a number of spatial points. These points are called control points, and their spatial coordinates are interpolated using a number of functions called base functions to construct a tridimensional surface.

The B-Spline parametric surface is constructed using the base functions. In turn, the base functions are built based on a parameterization of space given by the knot vectors, i. e.

$\text{Ξ}=\left[{\text{ξ}}_1,{\text{ξ}}_2,{\text{ξ}}_3\dots {\text{ξ}}_{\text{n}+\text{k}}\right].$ (5.1)

${\xi }_n$ represents the knots. The knots define the limits of the base functions domain, by dividing the parametric space into discrete elements. The knot vector has $n$ elements, $n$ being the sum of the number of control points plus the order k of the base function. The order of the base function represents the number of intervals between knots that composes the function’s domain. An open knot vector has its first and last entries repeated k times. The open knot vector makes the curve or surface to behave the same when close to its borders as it behaves in the middle.

A point in the B-Spline surface is given by

$Q\left(\text{t},\text{u}\right)={\displaystyle \sum }_{\text{i}=1}^{\text{n}+1}{\displaystyle \sum }_{\text{j}=1}^{\text{m}+1}{B}_{\text{i},\text{j}}{\text{N}}_{\text{i},\text{k}}\left(\text{ξ}\right){\text{M}}_{\text{j},\text{l}}\left(\text{η}\right),$ (5.2)

where $B$ is a control point, $\text{N}$ and $\text{M}$ are the base functions, one increasing in a different direction to the other. $\text{k}$ and $\text{l}$ are the order of the base functions. The B-Spline base function is given by

${\text{N}}_{\text{i}\mathrm{,1}}\left(\text{ξ}\right)=\{\begin{array}{c}1se{\text{x}}_{\text{i}}\le \xi <{\text{x}}_{\text{i}+1}\\ 0casocontrário\end{array}{\text{N}}_{\text{i},\text{k}}\left(\text{ξ}\right)=\frac{\left(\text{ξ}-{\text{x}}_{\text{i}}\right){\text{N}}_{\text{i},\text{k}-1}\left(\text{ξ}\right)}{{\text{x}}_{\text{i}+\text{k}-1}-{\text{x}}_{\text{i}}}+\frac{\left({\text{x}}_{\text{i}+\text{k}}-\text{t}\right){\text{N}}_{\text{i}+\mathrm{1,}\text{k}-1}\left(\text{ξ}\right)}{{\text{x}}_{\text{i}+\text{k}}-{\text{x}}_{\text{i}+1}}.$ (5.3)

The base functions are built recursively, with higher orders needing more calculations and thus processing power on a numerical analysis. An order of 1 reduces the base function to a Lagrangian base function.

B-Spline surfaces enjoy $C^{k-m_i}$ continuity, where k is the order of the base function, and m is the multiplicity of an inner knot $\text{i}$. In this paper, other than an open knot vector, no knot vector manipulation was used, so the multiplicity of any inner knot is 1.

More detailed information about B-Spline and NURBS curves can be found on (LES PIEGL, 1997) and (ROGERS, 2001).

5.2 B-Spline based finite elements

When the classic Lagrangian 4-node contact element is used, multiple contact elements are utilized to form a contact surface. When there’s contact between a slave node and the master surface, the resulting contact forces are distributed between the slave node, and the four nodes belonging to the single contact element where contact was detected.

On this paper, B-Spline based isogeometric analysis is utilized. A single B-Spline surface, or patch, is utilized to represent the whole contact surface on the master body. As such, a B-Spline surface is constructed using the nodes of the master’s surface contact area.

More details about isogeometric analysis can be found on (AUSTIN COTTRELL, 2009) and a bi-dimensional isogeometric analysis can be seen on (DE LORENZIS, 2011).

In case contact between a slave node and the master surface is found, the resulting contact forces are distributed among the slave node and contact surface nodes, using the B-Spline base functions as criteria for the distribution. As a general rule, the further away from the contact point the node is, smaller is the contact force acting on that node.

For the numerical examples executed, the methodology found in (BANDEIRA, 2001) was followed, but utilizing a B-Spline based finite element surface in place of a Lagrangian element contact surface mesh. As such, only the master surface is discretized with the B-Spline finite elements, the slave surface is not discretized, but each of its nodes is checked for contact against the master surface.

The knot vectors used to construct the B-Spline surface patch are constructed dynamically, based on the number of nodes taken as control points in the appropriate direction.

6 DISCRETIZATION OF CONTACT FORCES

Discretization of the contact work described in equation (4.1) can be written in vector form as

${\text{W}}_{\text{C}}\left(u,v\right)={\displaystyle \sum }_{\text{i}=1}^{\text{ns}}{\text{A}}_{\text{i}}{R}_{\text{contact i}}\cdot V_{\text{contact i}},$ (6.1)

where ${\text{A}}_{\text{i}}$ represents the contact finite element area, $V_{\text{contact i}}$ is the total number of slave nodes, $V_{\text{contact i}}$ represents the derivative of the spatial position of the element nodes and $R_{\text{contact i}}$ represents the contact stresses.

Expressing this in a vector representation

$R_{\text{contact}}={\text{t}}_{\text{N}}N+{\text{t}}_{{\text{T}}_{\text{α}}}{D}^{\text{α}}={\text{t}}_{\text{N}}N+{\text{t}}_{{\text{T}}_1}{D}^1+{\text{t}}_{{\text{T}}_2}{D}^2$ (6.2)

with the vectors being defined by

${\mathbf{V}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{}\mathrm{i}}=\left[\begin{array}{c}{\mathbf{v}}_{\mathrm{s}}\\ {\mathbf{v}}_{\frac{1}{1}}\\ \begin{array}{c}{\mathbf{v}}_{\frac{2}{1}}\\ {\mathbf{v}}_{\frac{3}{1}}\\ \begin{array}{c}\dots \\ {\mathbf{v}}_{\mathrm{a}-\frac{2}{\mathrm{b}}}\\ \begin{array}{c}{\mathbf{v}}_{\mathrm{a}-\frac{1}{\mathrm{b}}}\\ {\mathbf{v}}_{\frac{\mathrm{a}}{\mathrm{b}}}\end{array}\end{array}\end{array}\end{array}\right]\mathrm{}\mathbf{N}=\left[\begin{array}{c}{\stackrel{-}{\mathbf{n}}}_{\mathbf{c}}\\ -{{\mathrm{Q}}_{\frac{1}{1}}\stackrel{-}{\mathbf{n}}}_{\mathbf{c}}\\ \begin{array}{c}-{\mathrm{Q}}_{\frac{2}{1}}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ -{{\mathrm{Q}}_{\frac{3}{1}}\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ \begin{array}{c}\dots \\ -{\mathrm{Q}}_{\mathrm{a}-\frac{2}{\mathrm{b}}}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ \begin{array}{c}-{\mathrm{Q}}_{\mathrm{a}-\frac{1}{\mathrm{b}}}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ -{\mathrm{Q}}_{\frac{\mathrm{a}}{\mathrm{b}}}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\end{array}\end{array}\end{array}\end{array}\right]{\mathbf{N}}_{\mathbf{\alpha }}=\left[\begin{array}{c}0\\ -{{{\mathrm{Q}}_{\frac{1}{1}}}_{,\mathrm{\alpha }}\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ \begin{array}{c}-{{{\mathrm{Q}}_{\frac{2}{1}}}_{,\mathrm{\alpha }}\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ -{{\mathrm{Q}}_{\frac{3}{1}}}_{,\mathrm{\alpha }}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ \begin{array}{c}\dots \\ {{-\mathrm{Q}}_{\mathrm{a}-\frac{2}{\mathrm{b}}}}_{,\mathrm{\alpha }}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ \begin{array}{c}{{-\mathrm{Q}}_{\mathrm{a}-\frac{1}{\mathrm{b}}}}_{,\mathrm{\alpha }}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\\ -{{\mathrm{Q}}_{\frac{\mathrm{a}}{\mathrm{b}}}}_{,\mathrm{\alpha }}{\stackrel{-}{\mathbf{n}}}_{\mathrm{c}}\end{array}\end{array}\end{array}\end{array}\right]{\mathbf{T}}_{\mathrm{\alpha }}=\left[\begin{array}{c}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ {{-\mathrm{Q}}_{\frac{1}{1}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ \begin{array}{c}{{-\mathrm{Q}}_{\frac{2}{1}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ {{-\mathrm{Q}}_{\frac{3}{1}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ \begin{array}{c}\dots \\ {{-\mathrm{Q}}_{\mathrm{a}-\frac{2}{\mathrm{b}}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ \begin{array}{c}{{-\mathrm{Q}}_{\mathrm{a}-\frac{1}{\mathrm{b}}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ {{-\mathrm{Q}}_{\frac{\mathrm{a}}{\mathrm{b}}}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\end{array}\end{array}\end{array}\end{array}\right]$

${\mathrm{D}}^{\mathrm{\alpha }}=\frac{1}{{\stackrel{-}{\mathrm{a}}}_{\mathrm{\alpha }\mathrm{\beta }}-{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\stackrel{-}{\mathrm{b}}}_{\mathrm{\alpha }\mathrm{\beta }}}\left[{\mathrm{T}}_{\mathrm{\beta }}-{\mathrm{g}}_{\stackrel{-}{\mathrm{N}}}{\mathrm{N}}_{\mathrm{\beta }}\right]\mathrm{}\mathrm{\delta }{\stackrel{-}{\mathrm{\xi }}}^{\mathrm{\alpha }}={\mathrm{D}}^{\mathrm{\alpha }}\bullet {\mathrm{V}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{}\mathrm{i}}$

where ${\text{Q}}_{\text{a}/\text{b}}$ is defined as

${\text{Q}}_{\text{a}/\text{b}}={\text{N}}_{\text{a},\text{k}}\left(\text{t}\right){\text{M}}_{\text{b},\text{l}}\left(\text{u}\right),$ (6.3)

with “a” and “b” being the index number of the control points that construct the surface patch in each direction, and “k” and “l” are the selected base function’s order, omitted for brevity.

The linearization of equation (6.1) is necessary to solve the problem for the displacements, so it’s given by

$\frac{\partial {\text{W}}_{\text{c}}}{\partial u}\Delta u={\displaystyle \sum }_{\text{i}=1}^{\text{ns}}{\text{A}}_{\text{i}}{\left(V_{\text{contact i}}\right)}^{\text{T}}{K}_{\text{contact i}}\Delta U_{\text{contact i}}$ (6.4)

with $K_{\text{contact i}}$ representing the stiffness matrix of the element.

To create the stiffness matrix contribution due to contact forces, a discretization of the equation is necessary. This discretization, in vector form, is given by

$K_{\text{N}}={ϵ}_{\text{N}}N{N}^{\text{T}}+{\text{t}}_{\text{N}}\left[N_{\text{α}}{D}^{\text{αT}}+{\text{a}}^{\text{βα}}{T}_{\text{α}}\left(N_{\text{B}}^{\text{T}}-D^{\text{γT}}\left({\overline{\text{n}}}_{\text{c}}\cdot {\overline{\text{a}}}_{\text{β},\text{γ}}\right)\right)\right]$ (6.5)

and the additional vectors necessary to write the equation are defined by

$∆{\mathbf{U}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{}\mathrm{i}}=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{∆\mathbf{u}}_{\mathbf{S}}\\ {∆\mathbf{u}}_{1/1}\end{array}\\ {∆\mathbf{u}}_{2/1}\\ \begin{array}{c}{∆\mathbf{u}}_{3/1}\\ \dots \\ {∆\mathbf{u}}_{\mathrm{a}-2/\mathrm{b}}\end{array}\end{array}\\ {∆\mathbf{u}}_{\mathrm{a}-1/\mathrm{b}}\\ {∆\mathbf{u}}_{\mathrm{a}/\mathrm{b}}\end{array}\right]\mathrm{}{\mathbf{P}}_{\mathrm{\alpha }}=\left[\begin{array}{c}0\\ -{{\mathrm{Q}}_{1/1}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\\ \begin{array}{c}-{{\mathrm{Q}}_{2/1}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\\ -{{\mathrm{Q}}_{3/1}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\\ \begin{array}{c}\begin{array}{c}\dots \\ {{-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\\ {{-\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\end{array}\\ {{-\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}}_{,\mathrm{\alpha }}{\mathbf{t}}_{\mathrm{t}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\end{array}\end{array}\end{array}\right]\mathrm{}{\mathbf{T}}_{\mathrm{\alpha }\mathrm{\beta }}=\left[\begin{array}{c}0\\ {{-\mathrm{Q}}_{1/1}}_{,\mathrm{\beta }}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ \begin{array}{c}-{{\mathrm{Q}}_{2/1}}_{,\mathrm{\beta }}\mathrm{}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ {{-\mathrm{Q}}_{3/1}}_{,\mathrm{\beta }}{\mathrm{}\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ \begin{array}{c}\dots \\ \begin{array}{c}{{-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}}_{,\mathrm{\beta }}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\\ {{-\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}}_{,\mathrm{\beta }}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\end{array}\\ {{-\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}}_{,\mathrm{\beta }}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha }}\end{array}\end{array}\end{array}\right]\mathrm{}{\widehat{\mathbf{T}}}_{\mathrm{\alpha }\mathrm{\beta }}=\left[\begin{array}{c}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\\ -{\mathrm{Q}}_{1/1}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\\ \begin{array}{c}-{\mathrm{Q}}_{2/1}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\\ {-\mathrm{Q}}_{3/1}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\\ \begin{array}{c}\begin{array}{c}\dots \\ {-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\end{array}\\ {-\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\\ {-\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}{\stackrel{-}{\mathbf{a}}}_{\mathrm{\alpha },\mathrm{\beta }}\end{array}\end{array}\end{array}\right]\mathrm{}{\mathbf{N}}_{\mathrm{\alpha }\mathrm{\beta }}=\left[\begin{array}{c}0\\ -{{{\mathrm{Q}}_{1/1}}_{,\mathrm{\alpha }\mathrm{\beta }}\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\\ \begin{array}{c}-{{{\mathrm{Q}}_{2/1}}_{,\mathrm{\alpha }\mathrm{\beta }}\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\\ -{{\mathrm{Q}}_{3/1}}_{,\mathrm{\alpha }\mathrm{\beta }}{\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\\ \begin{array}{c}\dots \\ {{-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}}_{,\mathrm{\alpha }\mathrm{\beta }}{\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\\ \begin{array}{c}{{-\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}}_{,\mathrm{\alpha }\mathrm{\beta }}{\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\\ -{{\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}}_{,\mathrm{\alpha }\mathrm{\beta }}{\stackrel{-}{\mathrm{n}}}_{\mathrm{c}}\end{array}\end{array}\end{array}\end{array}\right]\mathrm{}\mathbf{E}=\left[\begin{array}{c}1\\ -{\mathrm{Q}}_{1/1}1\\ \begin{array}{c}-{\mathrm{Q}}_{2/1}1\\ -{\mathrm{Q}}_{3/1}1\\ \begin{array}{c}\dots \\ \begin{array}{c}{-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}1\\ -{\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}1\end{array}\\ -{\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}1\end{array}\end{array}\end{array}\right]\mathrm{}{\mathbf{E}}_{\mathrm{\alpha }}=\left[\begin{array}{c}\widehat{0}\\ -{{\mathrm{Q}}_{1/1}}_{,\mathrm{\alpha }}1\\ \begin{array}{c}-{{\mathrm{Q}}_{2/1}}_{,\mathrm{\alpha }}1\\ -{{\mathrm{Q}}_{3/1}}_{,\mathrm{\alpha }}1\\ \begin{array}{c}\dots \\ \begin{array}{c}{{-\mathrm{Q}}_{\mathrm{a}-2/\mathrm{b}}}_{,\mathrm{\alpha }}1\\ {{-\mathrm{Q}}_{\mathrm{a}-1/\mathrm{b}}}_{,\mathrm{\alpha }}1\end{array}\\ {{-\mathrm{Q}}_{\mathrm{a}/\mathrm{b}}}_{,\mathrm{\alpha }}1\end{array}\end{array}\end{array}\right]\mathrm{}$ (6.6)

7 NUMERICAL EXAMPLES

Three numerical examples were considered to test the B-Spline based finite element contact algorithm. The first one was the Parisch model (PARISCH, 1989), which was executed to compare directly with the results obtained by the Lagrangian finite element under large displacements. The second one was an authors created example, two beams subject to a significant bending moment in which both show large displacements. Finally, the third example shows two parallel beams, and it was designed to test the B-Spline finite element when the contact surface is generated after the body is significantly deformed. Lastly, the final example was designed to exhibit an amount of sliding over the contact surface.

7.1 Parisch

In Figure 2, one can see the example created by (PARISCH, 1989) and referenced in (BANDEIRA, 2001). It consists of two blocks, shaped like square based prisms, one supporting the other. The inferior block has all its base nodes restricted in all directions, while none of the superior block’s nodes has any restrictions. The superior block is supported only by contact forces. There are 50 nodes and 13 elements on this example, and more geometry information can be seen on Figure 2.

This example was executed without any tangential forces (no friction). The normal penalty parameter was ${\text{ε}}_{\text{N}}={10}^5$, the blocks elasticity modules were ${10}^9$ Pa for the superior block and ${10}^{10}$ Pa for the inferior block, and it needed only a single load increment to reach convergence.

Figure 2
Parisch example

The resulting deformations can be observed on Figure 3. The results obtained were very close to the values obtained by Parisch. The vertical coordinate of the node where the load was applied is 21.95 on the Parisch article, and the value found on this paper was 21.94, a negligible difference.

Figure 3
Parisch displacements

Figure 4 shows the maximum and minimum main stresses acting on the Parisch example. ANSYS doesn’t run this example, because it considers it not constrained enough.

Figure 4
Parisch maximum and minimum main stresses

7.2 Crossed Beams

The second example can be observed on Figure 5. Created by the authors, the example consists of two beams crossed and superposed, with a 0.2 mm gap between them. On the free end of the superior beam, at the edge, a load totaling 9000N is applied. This example has 288 nodes and 120 elements. More information about the example geometry can be observed on Figure 5.

Figure 5
Crossed Beams example

This example was simulated with friction, the contact penalty parameters were ${\text{ε}}_{\text{N}}={10}^5$ for normal penalty and ${\text{ε}}_{\text{T}}={10}^3$ for the tangential penalty. The superior beam has an elasticity modulus of ${4.10}^{11}$ Pa, while the inferior beam has ${10}^{12}$ Pa. The friction coefficient used was $\text{μ}=\mathrm{0.2,}$ and three load increments were applied.

The resulting displacements can be observed in Figure 6 and Figure 7.

Figure 6
Crossed Beams displacements

Figure 7
Crossed Beams displacements (lateral view)

In Figure 8, one can observe the results obtained on ANSYS for this example, with the same loads, number of nodes and elements.

Figure 8
Crossed Beams, ANSYS displacements

The obtained results for maximum and minimum main stresses for the crossed beams example are shown in figures 9 and 10.

Figure 9
Crossed Beams maximum main stress

Figure 10
Crossed Beams minimum main stress

Figures 11 and 12 show the maximum and minimum stresses obtained by ANSYS.

Figure 11
Crossed Beams ANSYS maximum main stress

Figure 12
Crossed Beams ANSYS minimum main stress

7.3 Parallel Beams

Illustrated in Figure 13, this example consists in two parallel beams, restricted on both ends, with forces applied perpendicularly. There are 384 and 180 elements in this example. The loads were applied in a way to arch the beams before contact happens, so the contact surfaces are curved.

The superior beam has double the width of the inferior one; otherwise they’re the same size. Information about restricted nodes, load positioning and details about the dimensions can be observed on Figure 8. The elasticity modulus of both beams is ${10}^{11}$ Pa. Friction was used and the friction coefficient is $\text{μ}=0.2$. 30 load increments were utilized, and the normal and tangential penalty values are ${\text{ε}}_{\text{N}}={10}^8$ and ${\text{ε}}_{\text{T}}={10}^4$, respectively.

Figure 13
Parallel Beams example

The deformations can be observed on Figure 14 and Figure 15.

Figure 14
Parallel Beams displacements

Figure 15
Parallel Beams displacements (frontal view)

In Figure 16, one can observe the ANSYS results, again with the same loads and discretization.

Figure 16
Parallel Beams ANSYS displacements

Figures 17 and 18 show the maximum and minimum main stresses obtained.

Figure 17
Parallel Beams maximum main stress

Figure 18
Parallel Beams minimum main stress

Figures 19 and 20 show the maximum and minumum main stresses obtained on ANSYS.

Figure 19
Parallel Beams ANSYS maximum main stress

Figure 20
Parallel Beams ANSYS minimum main stress

7.4 Supported Beam

The last example shown, observed in Figure 21, consists of a beam fixed on one end, suspended over another beam which is fixed on both ends. There’s a 0.5 mm gap between them, and there are twelve 150N loads applied on the nodes of the free end of the superior beam. Details about the loads positions and the example’s geometry can be observed on Figure 12. This example consists of 441 nodes and 200 elements.

Figure 21
Supported Beam example

The elasticity modulus for both beams is ${10}^{11}$ Pa, the normal penalty utilized was ${\text{ε}}_{\text{N}}=5.{10}^3$ and 30 load increments were necessary for convergence.

The displacement results can be observed on Figure 22 and Figure 23.

Figure 22
Supported Beam (frontal view)

Figure 23
Supported Beam (Isometric view)

The ANSYS results from the same mesh can be observed on Figure 24.

Figure 24
Supported Beam ANSYS displacements

In figures 25 and 26 one can see the maximum and minimum main stresses for the Supported Beam example.

Figure 25
Supported Beam maximum main stress

Figure 26
Supported Beam minimum main stress

In figures 27 and 28, the ANSYS maximum and minimum main stresses are shown.

Figure 27
Supported Beam ANSYS maximum main stress

Figure 28
Supported Beam ANSYS minimum main stress

8 CONCLUSION

In this paper, problems exhibiting large displacements, curved contact surfaces, as well as large sliding were presented. In all cases presented, convergence was attained, and the displacements obtained were consistent with the results generated in ANSYS, as one can observe on table 1. The Parisch numerical example isn’t constrained enough to run on ANSYS, and was left out of the comparisons.

Table 1
Result comparison

	Max displacement	Max ANSYS displacement	Percentage (result/ANSYS results)
Crossed beams	2.9357	3.0112	97.49%
Parallel beams	0.5365	0.5488	97.76%
Supported beam	1.0972	1.1081	99.02%

The stress values found on this paper diverged from the ones found on ANSYS by a fairly significant amount. The reason for this discrepance is a difference in the Neo-Hookean energy strain equation of the material. The equation used in this paper was developed by Ciarlet-Simo, as described on (WRIGGERS, 2008),

${\text{W}}_{\text{d}}\left({\text{I}}_{\text{c}},\text{J}\right)=\frac{\text{λ}}{2}\left[\frac{1}{2}\left({\text{J}}^2-1\right)-\ln \left(\text{J}\right)\right]+\frac{{\text{μ}}_{\text{d}}}{2}\left[\left({\text{I}}_{\text{c}}-3\right)-2\ln \left(\text{J}\right)\right].$ (8.1)

On equation 8.1, $\text{J}$ is the jacobian, and ${\text{μ}}_{\text{d}},\text{λ}$ are the Lamé parameters. The Neo-Hookean energy strain equation used by ANSYS, as found on its help file, is given by

${\text{W}}_{\text{d}}\left({\text{I}}_{\text{c}},\text{J}\right)=\frac{\text{K}}{2}{\left(\text{J}-1\right)}^2+\frac{1}{2}{\text{μ}}_{\text{d}}\left(b-I\right),$ (8.2)

where b is the left Cauchy Green tensor, and K is given by

$\text{K}=\text{λ}+\frac{2}{3}{\text{μ}}_{\text{d}}.$ (8.3)

The ANSYS stress results can still be very valuable for a qualitative comparison. One can observe that the stress distributions displayed on the results found on this paper and the ones derived from ANSYS are consistently similar.

It was shown that B-Spline based finite element and the discretization technique utilized on this paper are valid forms of smoothing the contact surface and solving the facetization problem. It provides consistent stresses across the contact surfaces and forbids penetration, especially with better discretized meshes.

The implementation of the B-Spline based finite element is relatively simple, and further studies can be done to compare performance and convergence between it and other existing smoothing techniques.

REFERENCES

• ALAIN BATAILLY, B. M. N. C. A comparative study between two smoothing strategies for the simulation of contact with large sliding. Comput Mech, p. 581-601, 2013.
• BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado) - Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.
• BANDEIRA, A. A. et al. Algoritmos de otimização aplicados à solução de sistemas estruturais não-lineares com restrições: uma abordagem utilizando os métodos da Penalidade e do Lagrangiano Aumentado. Exacta (São Paulo. Impresso), 8, 2010. 345-361.
• BATHE, K.-J. Finite Element Procedures. First. ed. New Jersey: Prentice-Hall, Englewood Cliffs, 1996.
• BERTSEKAS, D. P. Nonlinear programming. Belmont: Athena Scientific, 1995.
• CALLUM J. CORBETT, R. A. S. NURBS-enriched contact finite elements. Computer methods in applied mechanics and engineering, 2014.
• COOK, M. P. W. Concepts and Applications of Finite Element Analysis. 4th. ed. [S.l.]: John Wiley & Sons, INC., 2002.
• D.J. BENSON, J. O. H. A single sufrace contact algorithm for the post-buckling analysis of shell structures. Comput. Methods Appl. Mech. Eng 78 (1990) 141-163, 1990.
• DE LORENZIS, L.; WRIGGERS, P.; ZAVARISE, G. A mortar formulation for 3D large deformation contact using NURBS-based isogeometric analysis and the augmented Lagrangian method. Comput Mech, 49, 2012. 1-20.
• I. TEMIZER, P. W. T. J. R. H. Contact treatment in isogeometric analysis with NURBS. Comput. Methods Appl. Mech. Engrg., p. 1100-1112, 2011.
• J. AUSTIN COTTRELL, T. J. R. H. Y. B. Isogeometric Analysis. 1st. ed. [S.l.]: Wiley, 2009.
• L. DE LORENZIS, I. T. P. W. G. Z. A large deformation frictional contact formulation using NURBS-based isogeometric analysis. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, 2011.
• LAURSEN, T. A. Computational Contact and Impact Mechanics. [S.l.]: Springer, 2002.
• LAURSEN, T. A.; SIMO, J. C. A Continuum-Based Finite Element Formulation for the Implicit Solution of Multibody, Large Deformation Frictional Contact Problems. Int. J. Num. Meth. Engng., 36, 1993a. 3451-3485.
• LES PIEGL, W. T. The NURBS Book. 2nd. ed. [S.l.]: Springer, 1997.
• LUENBERGER, D. G. Linear and Nonlinear Programming. 2. ed. Reading. ed. Massachusetts: Addison-Wesley Publishing Company, 1984.
• M.E. MATZEN, T. C. M. B. A point to segment contact formulation for isogeometric, NURBS based finite elements. Elselvier, 2012.
• PARISCH, H. A Consistent Tangent Stiffness Matrix for Three-Dimensional Non-Linear Contact Analysis. International Journal for Numerical Methods in Engineering, 28, 1989. 1803-12.
• PUSO, M. A.; LAURSEN, T. A.; SOLBERG, J. A segment-to-segment mortar contact method for quadratic elements and large deformations. Comput Methods Appl Mech Eng, 197, 2008. 555-566.
• ROGERS, D. F. An Introduction to NURBS. [S.l.]: Morgan Kaufman Publishers, 2001.
• SIMO, J. C.; LAURSEN, T. A. An augmented Lagrangian Treatment of Contact Problems Involving Friction. Computers & Structures, 42, 1992. 97-116.
• WRIGGERS, P. Computational Contact Mechanics. Second Edition. ed. Hannover, Germany: Springer-Verlag Berlin Heidelberg, 2006.
• WRIGGERS, P. Nonlinear Finite Element Methods. Hannover, Germany: Springer-Verlag Berlin Heidelberg, 2008.
• ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The Finite Element Method for Solid and Structural Mechanics. 6. ed. Oxford: Elsevier Butterworth-Heinemann, 2005.

Publication Dates

• Publication in this collection
  2018

History

• Received
  02 Feb 2018
• Reviewed
  16 Apr 2018
• Accepted
  20 Apr 2018