Sensitivity analysis of 3D frictional contact with BEM using complex-step differentiation

Abstract

This paper presents a study of the complex step differentiation method applied to a parameter sensitivity analysis for 3D elastic contact problem. The analysis is performed with the Boundary Element Method (BEM) using discontinuous elements and the Generalized Newton Method with line search (GNMls). A standard BEM implementation is used and the contact restrictions are fulfilled through the augmented Lagrangian method. This methodology in conjunction with the BEM avoids the calculation of the nonlinear derivatives during the solution process, allowing a fast and reliable solution procedure. The parameter sensitivity is evaluated using complex-step differentiation. This well-known method approximates the derivative of a function analogously to the standard finite differences method, with the advantages of being numerically exact and nearly insensitive to the step-size. As an example, a Hertz-type problem is solved and the sensitivity of the contact pressures with respect to the Young Modulus variation is evaluated. The obtained results are compared with analytical and numerical solutions found in the literature.

Keywords
boundary element method; frictional contact; complex step method; sensitivity analysis

1 INTRODUCTION

Contact problems are often found in engineering applications. While some cases can be simplified or even assumed to be irrelevant, there are others where the contact itself is the reason of the engineering problem. Wear, tear, fatigue, friction, among others, are all problems which arise as a consequence of the contact occurrence. Advanced and high performance materials, which are widely used in engineering, increase the need to predict the contact conditions as well as the mechanical response, both locally and globally.

Contact problems are often found in engineering applications. While some cases can be simplified or even assumed to be irrelevant, there are others where the contact itself is the reason of the engineering problem. Wear, tear, fatigue, friction, among others, are all problems which arise as a consequence of the contact occurrence. Advanced and high performance materials, which are widely used in engineering, increase the need to predict the contact conditions as well as the mechanical response, both locally and globally.

Unilateral contact imposes an ambiguous boundary condition, where the contact area and hence the equilibrium configuration, may have a nonlinear dependency on the loads acting on the structure. Since the works published by Fichera (1963, 1973), the contact between two solids is also known as the Signorini problem. The author obtained analytical solutions for an elastic sphere lying on a rigid plane, with and without friction, considering gravitational loads. Alart and Curnier (1991) state that the frictional contact problem is governed by a multi-valued tribological law which does not derive from a natural potential and therefore cannot be formulated as standard optimization problems using inequality constraints.

The Boundary Element Method (BEM) is well-known for its ability to solve contact problems, since its formulation intrinsically treats the displacements and tractions with same order of approximation. This enables the direct application of the contact constraints without the need of penalty parameters or Lagrange multipliers. Since the pioneer work of Andersson (1981), which used this property to develop an incremental loading technique to solve 2D contact problems, few other works are also found on the literature using the same principles, such as Paris and Garrido (1989), Garrido et al. (1991), Man et al. (1993a), Man et al. (1993b), Paris et al. (1995) and Garrido and Lorenzana (1998). An iterative procedure similar a predictor-corrector algorithm was used by Ghaderi-Panah and Fenner (1998). Other works using the Lagrange multiplier or the penalty parameter methods, which are mandatory for the Finite Element Method (FEM), are found in Gakwaya et al. (1992) and Yamazaki et al. (1994).

González et al. (2008) and Rodríguez-Tembleque and Abascal (2010) treated FEM-BEM coupled problems using the Augmented Lagrangian formulation to avoid some drawbacks of Lagrange multiplier and penalty methods. These approaches were based on the work of Alart and Curnier (1991). The contact restrictions are imposed in the form of projection functions, resulting in a very robust framework to both FEM and BEM frictional contact analysis. The nonlinear system of equations is then solved by the Generalized Newton Method with line search (GNMls) (Pang, 1990): a generalization of the standard Newton method to B-differentiable functions whose convergence is independent of the penalization parameter used. With an unconstrained optimization between each step of the Newton method it is also possible to accelerate the convergence. The method was also used by Rodriguez-Tembleque et al. (2011) to study 3D frictional contact on anisotropic media using the BEM. Among the existing contact treatment techniques, the Augmented Lagrangian is the most adequate in conjunction with BEM because it eliminates the need of trial and error contact state estimations (Rodríguez-Tembleque et al., 2008). This method imposes the contact restrictions exactly and treats the stick case like a standard two-region BEM formulation.

The idea of numerical differentiation by the complex step method (CSM) is due to Lyness and Moler (1967). Later, Squire and Trapp (1998) proved that if the increment is performed on the imaginary part, the resulting variation on the complex part of the function exactly approximates its derivative at that point. This avoids the evaluation of the function at two points as in finite differences. Another feature related to the use of complex variable is that the derivative is nearly independent of the increment size, which avoids common instabilities regarding the selection of this parameter. In Martins et al. (2003) an example is presented for the function $f\mathrm{(}x\mathrm{)=}{e}^x\mathrm{/(}\sqrt{{{\displaystyle \sin }}^3\mathrm{(}x\mathrm{)}+{{\displaystyle \cos }}^3\mathrm{(}x\mathrm{)}}\mathrm{)}$ where the relative error for the derivative value is very stable and close to the mantissa for increments of $1\times 10\text{}{\text{}}^{-9}$ to $1\times 10\text{}{\text{}}^{-300}$.

To the best of the author’s knowledge, the first application of the complex step method with BEM found on the literature is in Gao et al. (2002). It was used to obtain the stresses at internal points, avoiding the use of integral equations for stresses. Among other works, the methodology was used by Mundstock and Marczak (2009) to evaluate sensitivities in shape optimization problems, obtaining excellent results in comparison to analytical solutions.

This work describes a study of sensitivity in contact problems, applying the Complex Step Method with BEM. The paper is organized as follows. First, the contact problem is described along with the restrictions imposed. The boundary integral formulation leading to BEM equations is then presented. Next, the contact variables on the discrete form, the projection functions and the augmented variables are shown. The nonlinear system of equations is presented in a way similar to problems containing multiple regions. The GNMls and the Complex Step Method are also detailed. Results for a classical Hertz problem considering friction is solved and the results obtained are compared with its analytical solution and a previous work from the literature. The sensitivity of the contact pressure is analyzed as an application example and compared to the associated analytical derivative.

2 FORMULATION OF THE CONTACT PROBLEM

Let the problem of two elastic bodies ${\Omega }_1$ and ${\Omega }_2$, with boundary $\Gamma \mathrm{=}{\Gamma }_{\overline{p}}\cup {\Gamma }_{\overline{u}}\cup {\Gamma }_c$(disjoint set), where ${\Gamma }_{\overline{p}}$, ${\Gamma }_{\overline{u}}$, ${\Gamma }_c$, correspond to the portion of the boundary prescribed with tractions, displacements and possible contact interface, respectively, as illustrated in Fig. 1. As a boundary value problem, or the prescribed tractions, $\overline{p}$, or the prescribed displacements, $\overline{u}$, are specified along the boundaries, and the complementary variables, $u$ and $p$, respectively, are the unknowns.

Figure 1
Solids under contact and boundary conditions.

The possible contact boundary ${\Gamma }_c$ have conditions which depends on the contact state, defined by the gap $g$ between the two regions. When the gap positive, the surfaces are free and the tractions are null. Compatibility conditions must be forced when the gap aperture is zero and the tractions are not null, and they depend on the existence of friction or not.

2.1 Kinematic conditions

In this work, the BEM formulation assumes small strains and infinitesimal displacements and the contact regions are treated as a node-to-node contact. The nodes are positioned in a conforming scheme, as done in Rodríguez-Tembleque et al. (2008), forming pairs, i.e., the slave nodes (the ones which belong to region ${\Omega }_2$) are positioned as closely as possible to the master (region ${\Omega }_1$) nodes, or at least matching the displacement path performed by the contact node pair. The contact variables in the discrete form are then related to the possible contact node pairs, as in Fig. 2.

Figure 2
Master ($P_1$) and slave ($P_2$) nodes along a contact interface.

The gap $g$ for any pair of points in this problem is obtained through the following relation

$g=B^T(x_2-x_1)+B^T(u_2-u_1),$ (1)

where $x_1$ and $x_2$are coordinates on boundaries and respectively, while and are the displacements of these points, both referring to a global coordinate system. is the rotation matrix responsible for changing the variables to the local coordinate system with origin at, i.e.,

$B=\left[\begin{array}{ccc}{t}_1& t_2& n\end{array}\right],$ (2)

where $n$ is the outward normal vector at $x_1$ and $t_1\mathrm{,}{t}_2$ are the tangential vectors. Therefore the gap can be decomposed in its tangential and normal components:

$g=\left[\begin{array}{ccc}{\text{g}}_{t_1}& {\text{g}}_{t_2}& {\text{g}}_n\end{array}\right].$ (3)

2.2 Unilateral and frictional contact conditions

The unilateral and frictional laws for a boundary point are grouped as suggested in González et. al. (2008), resulting in the well-known contact conditions which depends on the normal gap $g_n$, normal $(t_n)$ and tangential $(t_t)$ tractions,

$\{\begin{array}{l}\begin{array}{lllll}{\text{t}}_n=0;\hfill & {\text{g}}_n\ge 0;\hfill & t_t=0;\hfill & \hfill & \text{No contact case}\text{,}\hfill \\ {\text{t}}_n\le 0;\hfill & {\text{g}}_n=0;\hfill & \Vert t_t\Vert =\mu \left|t_n\right|;\hfill & {\dot{g}}_t\cdot {\dot{t}}_t=-\Vert {\dot{g}}_t\Vert .\Vert t_t\Vert ;\hfill & \text{Contact - Slip case}\text{,}\hfill \\ {\text{t}}_n\le 0;\hfill & {\text{g}}_n=0;\hfill & {\dot{g}}_t=0;\hfill & \hfill & \text{Contact - Stick case}\text{.}\hfill \end{array}\hfill \end{array}$ (4)

The time rate in equation (4) must be approximated by a finite differences scheme. Rodríguez-Tembleque and Abascal (2010) suggested the following approximation for${\dot{g}}_t$at time ${\tau }_k$as

${\dot{g}}_t\simeq \frac{\Delta g_t}{\Delta \tau }.$ (5)

2 BOUNDARY ELEMENT METHOD FORMULATION

The displacement boundary integral equation for a point in the boundary is

$c_{lk}^i{u}_k^i+{\displaystyle \underset{\Gamma }{\int }}{p}_{lk}^{\mathrm{*}}{u}_{k}d\Gamma \mathrm{=}{\displaystyle \underset{\Gamma }{\int }}{u}_{lk}^{\mathrm{*}}{p}_{k}d\Gamma \mathrm{,}$ (6)

where the boundary is always smooth due to the use of discontinuous elements, i.e., $c_{lk}^i\mathrm{=}\frac{1}{2}$. In order to obtain a numerical solution from Eq. (6), the boundary is discretized into a finite number of elements.

The geometry is interpolated by standard shape functions written in terms of the parametric coordinates $\xi =({\xi }_1,{\xi }_2)$:

$x\mathrm{=}{\displaystyle \sum _{n\mathrm{=1}}^N}{\varphi }_n\mathrm{(}\xi \mathrm{)}{x}_n^j\mathrm{,}$ (7)

where $N$ is the number of nodes of the element and ${\varphi }_n$ is a vector containing the shape functions evaluated at $\xi$. The variables are calculated at collocation points placed at an offset distance from the geometrical nodes. Displacements and tractions are interpolated with discontinuous interpolation functions $\overline{\varphi }$:

$u\mathrm{=}{\displaystyle \sum _{n\mathrm{=1}}^N}{\overline{\varphi }}_n\mathrm{(}\xi \mathrm{)}{u}_n^j\mathrm{,}p\mathrm{=}{\displaystyle \sum _{n\mathrm{=1}}^N}{\overline{\varphi }}_n\mathrm{(}\xi \mathrm{)}{p}_n^j\mathrm{.}$ (8)

In the numerical implementation of this work linear (Q4) and quadratic (Q8) quadrilateral elements were used. The interpolation functions for the linear element are written in compact form as

$\begin{array}{c}{\overline{\varphi }}_n\mathrm{=}\frac{1}{4d^2}\mathrm{(}d+{\xi }_1{\overline{\xi }}_1^n\mathrm{)(}d+{\xi }_2{\overline{\xi }}_2^n\mathrm{),}n\mathrm{=1}\dots \mathrm{4,}\end{array}$ (9)

where $\mathrm{(}{\overline{\xi }}_1^n\mathrm{,}{\overline{\xi }}_2^n\mathrm{)}$ correspond to the nodal coordinates of the n-th geometrical node, $\mathrm{(}{\xi }_1\mathrm{,}{\xi }_2\mathrm{)}$ defines the coordinates where the variable is being calculated and $d$ denotes the distance from the center of the element to the physical nodes (collocation points) in local coordinates. This distance $d$ is related to the offset ($a$) by the relation:

$d\mathrm{=1}-\frac{\text{a}}{100}\mathrm{,}$ (10)

which reduces to the continuous functions ${\varphi }_n$ when $\text{a}\mathrm{=0}$. The interpolation functions for the 8-node element are also written in compact form,

$\begin{array}{lll}{\overline{\varphi }}_n\hfill & \mathrm{=}\frac{\mathrm{(}d+{\xi }_1{\overline{\xi }}_1^n\mathrm{)(}d+{\xi }_2{\overline{\xi }}_2^n\mathrm{)(}{\xi }_1{\overline{\xi }}_1^n+{\xi }_2{\overline{\xi }}_2^n-d\mathrm{)}}{4d^3}\hfill & n\mathrm{=1}\dots \mathrm{4,}\hfill \\ {\overline{\varphi }}_n\hfill & \mathrm{=}\frac{\mathrm{(}d+{\xi }_1\mathrm{)(}d-{\xi }_1\mathrm{)(}d+{\xi }_2\mathrm{)(}d-{\xi }_2\mathrm{)}}{2d^3\mathrm{(}d-{\xi }_1{\overline{\xi }}_1^n-{\xi }_2{\overline{\xi }}_2^n\mathrm{)}}\mathrm{,}\hfill & n\mathrm{=5}\dots 8.\hfill \end{array}$ (11)

Collocation for node $\text{i}$, generates the discretized form of Eq. (6):

$C_i{u}_i+H_{ij}^e{u}^e=G_{ij}^e{p}^e.$ (12)

Performing the collocation process over all boundary nodes results in the algebraic system of equations that leads to the BEM solution, i.e.,

$\begin{array}{l}Hu\mathrm{=}Gp\mathrm{.}\hfill \end{array}$ (13)

4 CONTACT PROJECTION FUNCTIONS

The formulation used for the imposition of contact restrictions was implemented according to González et. al. (2008). The handling of inequalities, typical of contact conditions, was avoided using suitable projection functions.

4.1 Unilateral and frictional contact conditions

To fulfill the non-penetration condition, the complementarity of the normal components of the gap and tractions variables are written in the following manner

$\{\begin{array}{ll}{\text{t}}_n{\text{g}}_n=\mathrm{0,}\hfill & \hfill \\ {\text{t}}_n\le \mathrm{0,}\hfill & {\text{g}}_n\ge \mathrm{0,}\hfill \end{array}$ (14)

which means that not only one of the two variables must be null, whether is the normal tractions or the normal gap, but also that the gap must be no less than zero and the normal traction to be never positive. The above inequalities can be eliminated using the following projection function

${\mathbb{P}}_{{ℝ}_{-}}(v):ℝ\to {ℝ}_{-},$ (15)

projecting the variable $v$ onto ${ℝ}_{-}$, the admissible region of the contact normal tractions, i.e.

${\mathbb{P}}_{{ℝ}_{-}}(v)=\min (v\mathrm{,0}).$ (16)

A mixed variable is now introduced, the augmented normal traction $t_n^{\mathrm{*}}$, defined as

$t_n^{*}=t_n+r_n{\text{g}}_n,$ (17)

where $r_n$ is a positive penalization parameter $(r\in {ℝ}_{+})$. The complementarity condition can then be fulfilled using a single equation,

$t_n-{\mathbb{P}}_{{ℝ}_{-}}\mathrm{(}{t}_n^{\mathrm{*}}\mathrm{)=0.}$ (18)

4.2 Unilateral and frictional contact conditions

Similarly, the augmented tangential traction vector is defined as

$t_t^{\mathrm{*}}\mathrm{=}{t}_t-r_t{\dot{g}}_t\mathrm{,}$ (19)

the gap must be

${\mathbb{P}}_{{ℂ}_g}\mathrm{(}v\mathrm{)=}(\begin{array}{cc}v& \text{if}\Vert v\Vert \mathrm{<}\left|\mu t_n\right|\mathrm{,}\\ \left|\mu t_n\right|e_t& \text{if}\Vert v\Vert ⩾\left|\mu t_n\right|\mathrm{,}\end{array}$ (20)

where $e_t\mathrm{=}v/\Vert v\Vert$. This function projects the variable onto the Coulomb disk, ${ℂ}_g$, of radius $g\mathrm{=}\left|\mu t_n\right|$. Considering a plane defined by the tangential tractions, $t_{t_1}\perp t_{t_2}$, ${ℂ}_g$ represents admissible for the tangential tractions under a given normal traction $t_n$. The tangential contact restrictions are finally written as

$t_t-{\mathbb{P}}_{{ℂ}_g}\mathrm{(}{t}_t^{\mathrm{*}}\mathrm{)=0.}$ (21)

4.3 Normal-Tangential operator

In order to combine the normal and tangential restrictions in a single statement, the normal-tangential projection function ${\mathbb{P}}_{{ℂ}_f}$ is defined as

${\mathbb{P}}_{{ℂ}_f}\mathrm{(}{t}^{*}\mathrm{)=}\left[\begin{array}{c}{\mathbb{P}}_{{ℂ}_g}\mathrm{(}{t}_t^{\mathrm{*}}\mathrm{)}\\ {\mathbb{P}}_{{ℝ}_{-}}\mathrm{(}{t}_n^{\mathrm{*}}\mathrm{)}\end{array}\right]\mathrm{,}$ (22)

where the region ${ℂ}_f$ is the augmented friction cone, with radius $\left|\mu {\mathbb{P}}_{{ℝ}_{-}}\mathrm{(}{t}^{\mathrm{*}}{}_n\mathrm{)}\right|$ and inclination defined by the friction coefficient. The inequalities in Eq. (4) are then completely fulfilled though the following set of equations

$t-{\mathbb{P}}_{{ℂ}_f}\mathrm{(}{t}^{*}\mathrm{)=0.}$ (23)

5 CONTACT FORMULATION WITH BOUNDARY ELEMENTS

The system of linear equations (SLE) for a contact problem is assembled similarly as in multiple-region problems, resulting the system of equations,

$\begin{array}{c}{H}^{{\Gamma }_1^{nc}}{u}^{{\Gamma }_1^{nc}}-G^{{\Gamma }_1^{nc}}{p}^{{\Gamma }_1^{nc}}+H^{{\Gamma }_1^c}{u}^{{\Gamma }_1^c}-G^{{\Gamma }_1^c}{p}^{{\Gamma }_1^c}\mathrm{=}{\overline{G}}^{{\Gamma }_1^{nc}}{\overline{p}}^{{\Gamma }_1^{nc}}-{\overline{H}}^{{\Gamma }_1^{nc}}{\overline{u}}^{{\Gamma }_1^{nc}}\mathrm{,}\\ H^{{\Gamma }_2^{nc}}{u}^{{\Gamma }_2^{nc}}-G^{{\Gamma }_2^{nc}}{p}^{{\Gamma }_2^{nc}}+H^{{\Gamma }_2^c}{u}^{{\Gamma }_1^c}+G^{{\Gamma }_2^c}{p}^{{\Gamma }_1^c}\mathrm{=}{\overline{G}}^{{\Gamma }_2^{nc}}{\overline{p}}^{{\Gamma }_2^{nc}}-{\overline{H}}^{{\Gamma }_2^{nc}}{\overline{u}}^{{\Gamma }_2^{nc}}\mathrm{.}\end{array}$ (24)

where the compatibility conditions

$\begin{array}{l}{u}^{{\Gamma }_2^c}=u^{{\Gamma }_1^c},\\ p^{{\Gamma }_2^c}=-p^{{\Gamma }_1^c}.\end{array}$ (25)

must be enforced along the contacting interfaces.

Although Eq. (24) suffices to solve the interface problem, for a more general condition such as contact, the compatibility conditions have to be written differently than Eq. (25), taking in to account the current contact state. This leads to the introduction of two auxiliary sets of equations: the kinematic relations and the contact restrictions, introduced in the following sections.

5.1 Discrete contact variables

In order to write the contact tractions with a distinct symbol of the global ones and keeping consistency with the references notation, the local contact tractions for a node pair $i$ are defined as ${\lambda }_i$. In this work, only BEM is used, thus there exist no requirement of using Lagrange multipliers, as one of the well-known advantages of BEM in contact problems is that the tractions are already part of the unknowns. The transformation between the element tractions and contact tractions can be accomplished using the following relations

$p_j^{{\Gamma }_1^c}={\tilde{C}}_{ji}^1{\Lambda }_i,$ (26)

$p_k^{\Gamma 2c}=-{\tilde{C}}_{ki}^2{\Lambda }_i,$ (27)

where ${\tilde{C}}_{ji}^{\alpha }$ is formed with the rotation matrix $B$ assembled for each region .., according to the assembly incidence $j$ of each contact node $i$. That is, for a contact node pair $I\mathrm{=(}j\mathrm{,}k\mathrm{)}$ in Eq. (26) and $m\mathrm{=}I\mathrm{(}\alpha \mathrm{)}$, we have:

${\tilde{C}}_{mI}^{\alpha }=B_I.$ (28)

The discrete form of the gap $g$ is now introduced using the same transformation. The kinematic relations of the possible contact region are written as

$k=k_{go}+{({\tilde{C}}^2)}^T{u}^{\Gamma 2c}-{({\tilde{C}}^1)}^T{u}^{\Gamma 1c},$ (29)

where $k$ is the gap vector accounting for all the contact pairs. The vector $k_{go}$ is the sum of the initial gap and rigid body displacements vectors.

Substituting Eq. (26) and (27) in Eq. (24) and adding Eq. (29) to the SLE results

$\begin{array}{ll}{A}^{{\Gamma }_1^{nc}}{x}^{{\Gamma }_1^{nc}}+H^{{\Gamma }_1^c}{u}^{{\Gamma }_1^c}-{\tilde{C}}^1{G}^{{\Gamma }_1^c}\Lambda \hfill & \mathrm{=}{b}^{{\Gamma }_1^{nc}}\mathrm{,}\hfill \\ A^{{\Gamma }_2^{nc}}{x}^{{\Gamma }_2^{nc}}+H^{{\Gamma }_2^c}{u}^{{\Gamma }_2^c}+{\tilde{C}}^2{G}^{\Gamma 2c}\Lambda \hfill & \mathrm{=}{b}^{{\Gamma }_2^{nc}}\mathrm{,}\hfill \\ {({\tilde{C}}^1)}^T{u}^{\Gamma 1c}-{({\tilde{C}}^2)}^T{u}^{\Gamma 2c}+C_{g}k\hfill & \mathrm{=}{C}_g{k}_{go}\mathrm{,}\hfill \end{array}$ (30)

where $\Lambda$ are the BEM nodal tractions rotated from the global to the master nodes local coordinate system.

5.2 Contact restrictions

The discrete form of the contact restrictions can be incorporated to the SLE through the application of the projection functions on a node pair. Writing equation (23) for the discrete form of the contact tractions results in

$\lambda -{\mathbb{P}}_{{ℂ}_f}({\lambda }^{*})=0.$ (31)

Substituting ${\Lambda }^{*}=\Lambda +rk$ with $r$ being the penalization parameters vector, the previous equation becomes

$\lambda -{\mathbb{P}}_{{ℂ}_f}(\lambda +rk)=0.$ (32)

For the normal direction one writes the restriction as

${\lambda }_n-\min ({\lambda }_n{}^{*}\mathrm{,0})=\mathrm{0,}$ (33)

which, for the no contact case ${\mathrm{(}{\lambda }_n^{\mathrm{*}}\mathrm{)}}_I\ge 0$ and by the Coulomb friction law, $\Vert {\lambda }_t\Vert \le \mu \left|{\lambda }_n\right|$results in

${\lambda }_n\mathrm{=0,}{\lambda }_t\mathrm{=0.}$ (34)

Writing the restrictions for the contact case for the tangential direction, one needs to write the restrictions separately for stick and slip. Along the normal direction ${\mathrm{(}{\lambda }_n^{\mathrm{*}}\mathrm{)}}_I\mathrm{<0}$, so the restriction equation of a node pair $I$, for both stick and slip cases, is written as

${\lambda }_n-\min \mathrm{(}{\lambda }_n-r_n{k}_n\mathrm{,0)=0,}$ (35)

which leads to

$r_n{k}_n\mathrm{=0.}$ (36)

For the stick case ${({\lambda }_t^{\mathrm{*}}\mathrm{)}}^I\mathrm{<}\mu \left|{\lambda }_n\right|$ (the first case of equation (20)), the tangential restriction results

${\lambda }_t-{\mathbb{P}}_{{ℂ}_g}({\lambda }_t+r_t{k}_t)=\mathrm{0,}$ (37)

which, using ${\mathbb{P}}_{{ℂ}_g}\mathrm{(}v\mathrm{)=}v$, reads

$r_t{k}_t=0.$ (38)

As for the slip case, ${\mathrm{(}{\lambda }_t^{\mathrm{*}}\mathrm{)}}_I\ge \left|\mu {\lambda }_n\right|$ (second case of equation (20)), ${\mathbb{P}}_{{ℂ}_f}({\lambda }_t^{*})=\left|\mu t_n\right|(e_t^{*})$, is used, so the tangential restriction is expressed as

$\begin{array}{l}{\lambda }_{t_1}+\mu e_{t_1}^{*}{\lambda }_n=\mathrm{0,}\\ {\lambda }_{t_2}+\mu e_{t_2}^{*}{\lambda }_n=\mathrm{0,}\end{array}$ (39)

where

$e_t^{*}=\frac{{\lambda }_t^{*}}{\Vert {\lambda }_t^{*}\Vert }.$ (40)

Grouping Eqs. (34), (36), (38) and (39), in a single system for all contact restrictions, leads to

$P_{\lambda }\Lambda +P_{\text{g}}k=0\mathrm{,}$ (41)

where $P_{\lambda }$ and $P_{\text{g}}$ are assembled according to the contact state of each pair $I$ and depend on the contact pair state (Rodríguez-Tembleque and Abascal, 2010):

• No contact case: ${\left({\lambda }_n^{\mathrm{*}}\right)}_I⩾0$

${\mathrm{(}{P}_{\lambda }\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\mathrm{,(}{P}_{\text{g}}{\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\mathrm{.}$ (42)

• Stick case: ${\left({\lambda }_n^{\mathrm{*}}\right)}_I\mathrm{<0}$ and $\Vert {\lambda }_t^{\mathrm{*}}\Vert \mathrm{<}\mu \left|{\lambda }_{nI}^{\mathrm{*}}\right|$

${\mathrm{(}{P}_{\lambda }\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\mathrm{,(}{P}_{\text{g}}{\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}{r}_t& 0& 0\\ 0& r_t& 0\\ 0& 0& r_n\end{array}\right]\mathrm{.}$ (43)

• Slip case: ${\lambda }_{nI}^{\mathrm{*}}\mathrm{<0}$ and $\Vert {\lambda }_{tI}^{\mathrm{*}}\Vert \mathrm{<}\mu \left|{\lambda }_{nI}^{\mathrm{*}}\right|$

${\mathrm{(}{P}_{\lambda }\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}1& 0& \mu e_{t_1}^{\mathrm{*}}\\ 0& 1& \mu e_{t_2}^{\mathrm{*}}\\ 0& 0& 0\end{array}\right]\mathrm{,(}{P}_{\text{g}}{\mathrm{)}}_I\mathrm{=}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -r_n\end{array}\right]\mathrm{.}$ (44)

5.3 Contact restrictions

After incorporating the contact restrictions Eq. (41), Eq. (30) can be recast as a system of the type $Rz\mathrm{=}f$:

${\left[\begin{array}{cccccc}{A}^{{\Gamma }_1^{nc}}& H^{{\Gamma }_1^c}& 0& 0& {\tilde{C}}^1{G}^{{\Gamma }_1^c}& 0\\ 0& 0& A^{{\Gamma }_2^{nc}}& H^{{\Gamma }_2^c}& -{\tilde{C}}^2{G}^{{\Gamma }_2^c}& 0\\ 0& {({\tilde{C}}^1)}^T& 0& -{({\tilde{C}}^2)}^T& 0& I\\ 0& 0& 0& 0& P_{\lambda }& P_{\text{g}}\end{array}\right]}_{}^{}\left\{\begin{array}{c}{x}^{{\Gamma }_1^{nc}}\\ u^{{\Gamma }_1^c}\\ x^{{\Gamma }_2^{nc}}\\ u^{{\Gamma }_2^c}\\ \Lambda \\ k\end{array}\right\}=\left\{\begin{array}{c}{b}^{{\Gamma }_1^{nc}}\\ b^{{\Gamma }_2^{nc}}\\ k_{\text{go}}\\ 0\end{array}\right\}$ (45)

where $A^{{\Gamma }_{\alpha }^{nc}}$ are the columns from $H^{{\Gamma }_{\alpha }}$ or $-G^{{\Gamma }_{\alpha }}$ (Eq. (24)) corresponding to the independent unknowns $x^{{\Gamma }_{\alpha }^{nc}}$ for the $\alpha$ region. $H^{{\Gamma }_{\alpha }^c}$ are the matrices relative to the unknown tractions on the possible contact regions. The contact tractions ($\Lambda$), as well as the gap ($k$), are assumed in their local coordinate system, by the incorporation of the rotation matrices on the system of equations, i.e., Eq. (2), which are assembled for each contact pair on the main rotation matrix presented as ${\tilde{C}}^{\alpha }$ in Eq. (45). $I$ is the identity matrix and the vector $k_{\text{go}}$ accumulates the initial gap and rigid body displacements.

6 NONLINEAR SYSTEM SOLUTION: GENERALIZED NEWTON METHOD WITH LINE SEARCH

As suggested in Rodríguez-Tembleque and Abascal (2010), the nonlinear system of equations is solved using the Generalized Newton Method with line search, first proposed by Pang (1990). The $\mathcal{B}$-differentiable functions are described as non Fréchet-differentiable due mainly to the absence of linearity on the $\mathcal{B}$-derivative. Pang (1990) presents an example of complementarity function $H\mathrm{:}{ℝ}^n\to {ℝ}^n$, $H\mathrm{(}x\mathrm{)=}\min \mathrm{(}h\mathrm{(}x\mathrm{),}f\mathrm{(}x\mathrm{))}$, which is non Fréchet-differentiable, the same case as the normal projection operator. Contrary to the Uzawa method, the GNMls does not require damping or stabilization parameters and its convergence is independent of the penalization factor used in the augmented Lagrangian (Alart and Curnier, 1991). Denoting the current iteration by the superscript $\mathrm{(}n\mathrm{)}$, the GNMls algorithm can be summarized in the following steps:

1. Let $z^{\mathrm{(0)}}$ be an arbitrary initial vector, and

$\Theta (z^{(n)})=R^{(n)}{z}^{(n)}-f$ (46)

2. Given $z^{\mathrm{(}n\mathrm{)}}$ with $\text{Q}\mathrm{(}{z}^{\mathrm{(}n\mathrm{)}}\mathrm{)}\ne 0$, the direction $\Delta z^{\mathrm{(}n\mathrm{)}}$ is obtained by solving the equation:

$\Theta (z^{(n)})+\mathcal{B}\Theta (z^{(n)},\Delta z^{(n)})=0$ (47)

3. Find the first integer $m\mathrm{=1,2,}\dots$ which fulfills the decreasing error condition:

$\Psi {(z^{(n)}+{\alpha }^{(n)}\Delta z^{(n)})}^{}{⩽}_{}{(1-2\sigma {\alpha }^{(n)})}^{}\Psi (z^{(n)}),$

with ${\alpha }^{\mathrm{(}n\mathrm{)}}\mathrm{=}q\mathrm{.}{\beta }^m$, $q\mathrm{>0}$, $\beta \in \mathrm{(0,1)}$, $\sigma \in \mathrm{(0,1/2)}$, and

$\Psi (z^{(n)})=\frac{1}{2}{\Vert \Theta (z^{(n)})\Vert }^2.$ (48)

4. Updated solution vector: $z^{\mathrm{(}n+\mathrm{1)}}\mathrm{=}{z}^{\mathrm{(}n\mathrm{)}}+{\alpha }^{\mathrm{(}n\mathrm{)}}\Delta z^{\mathrm{(}n\mathrm{)}}$

5. If $\Psi (z^{(n+1)})\le {\epsilon }_1$, stop. Otherwise $n\mathrm{=}n+1$, and return to step 2.

6.1 Search direction: Directional derivatives

The system of equations (46) can be split in its linear and nonlinear parts,

$\Theta (z^{(n)})={\Theta }_{LD}(z^{(n)})+{\Theta }_{NLD}(z^{(n)}),$ (49)

where ${\Theta }_{LD}(z^{(n)})$ is the part which has the linear directional derivative at all points, and ${\Theta }_{NLD}(z^{(n)})$ has the nonlinear derivative at some points. The $\mathcal{B}$-derivative $\mathcal{B}\Theta (z^{(n)},\Delta z^{(n)})$ can be separated in the same manner:

$\mathcal{B}\Theta (z^{(n)},\Delta z^{(n)})=\left(\nabla {\Theta }_{LD}(z^{(n)})+\partial {\Theta }_{NLD}(z^{(n)})\right)\Delta z^{(n)},$ (50)

where $\nabla {\Theta }_{LD}(z^{(n)})$ is the linear part of the Jacobian matrix,

$\nabla {\Theta }_{LD}(z^{(n)})={\left[\begin{array}{cccccc}{A}^{{\Gamma }_1^{nc}}& H^{{\Gamma }_1^c}& 0& 0& {\tilde{C}}^1{G}^{{\Gamma }_1^c}& 0\\ 0& 0& A^{{\Gamma }_2^{nc}}& H^{{\Gamma }_2^c}& -{\tilde{C}}^2{G}^{{\Gamma }_2^c}& 0\\ 0& {({\tilde{C}}^1)}^T& 0& -{({\tilde{C}}^2)}^T& 0& I\\ 0& 0& 0& 0& 0& 0\end{array}\right]}_{}^{},$ (51)

and $\partial {\Theta }_{NLD}(z^{(n)})$ is the nonlinear part of the Jacobian. A linear approximation for this derivative can be used (Strömberg, 1997):

$\partial {\Theta }_{NLD}(z^{(n)})={\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& J_{\lambda }& J_g\end{array}\right]}_{}^{(n)},$ (52)

where $J_{\lambda }^{(n)}$ and $J_{\lambda }^{(n)}$ are written for each contact pair, as are the matrices $P_{\lambda }$ and $P_g$:

• No contact case: ${\lambda }_{nI}^{\mathrm{*}}⩾0$

${\mathrm{(}{J}_{\lambda }\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]}_I\mathrm{,}{\mathrm{(}{J}_{\text{g}}\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]}_I\mathrm{.}$ (53)

• Stick case: ${\lambda }_{nI}^{\mathrm{*}}\mathrm{<0}$ and $\Vert {\lambda }_{tI}^{\mathrm{*}}\Vert \mathrm{<}\mu \left|{\lambda }_{nI}^{\mathrm{*}}\right|$

${\mathrm{(}{J}_{\lambda }\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]}_I\mathrm{,}{\mathrm{(}{J}_{\text{g}}\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}{r}_t& 0& 0\\ 0& r_t& 0\\ 0& 0& r_n\end{array}\right]}_I\mathrm{.}$ (54)

• Slip case: ${\lambda }_{nI}^{\mathrm{*}}\mathrm{<0}$ and $\Vert {\lambda }_{tI}^{\mathrm{*}}\Vert \mathrm{>}\mu \left|{\lambda }_{nI}^{\mathrm{*}}\right|$

${\mathrm{(}{J}_{\lambda }^{\mathrm{(}n\mathrm{)}}\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}{\Psi }_{11}& {\Psi }_{12}& \mu {\omega }_{t_1}^{\mathrm{*}}\\ {\Psi }_{21}& {\Psi }_{22}& \mu {\omega }_{t_2}^{\mathrm{*}}\\ 0& 0& 0\end{array}\right]}_I^{\mathrm{(}n\mathrm{)}}\mathrm{,}{\mathrm{(}{J}_{\text{g}}^{\mathrm{(}n\mathrm{)}}\mathrm{)}}_I\mathrm{=}{\left[\begin{array}{ccc}-r_t{\tilde{\Psi }}_{11}& -r_t{\tilde{\Psi }}_{12}& 0\\ -r_t{\tilde{\Psi }}_{21}& -r_t{\tilde{\Psi }}_{22}& 0\\ 0& 0& -r_n\end{array}\right]}_I^{\mathrm{(}n\mathrm{)}}\mathrm{,}$ (55)

where $\Psi =(1+\alpha )I-\beta R$,$\tilde{\Psi }=\Psi -I$,$R={({\Lambda }_t^{*(n)})}_I\otimes {({\Lambda }_t^{*(n)})}_I$, and

$\alpha =\frac{\mu {({\Lambda }_n^{*(n)})}_I}{\Vert {({\Lambda }_t^{*(n)})}_I\Vert },\beta =\frac{\mu {({\Lambda }_n^{*(n)})}_I}{{\Vert {({\Lambda }_t^{*(n)})}_I\Vert }^3}$

6.2 Linearized Derivatives

It is possible to simplify the contact nonlinear system even further by assuming that the nonlinear part of the directional derivatives ${\mathrm{(}{J}_{\lambda }^{\mathrm{(}n\mathrm{)}}\mathrm{)}}_I$ and ${\mathrm{(}{J}_g^{\mathrm{(}n\mathrm{)}}\mathrm{)}}_I$ can be approximated by the projection functions ${\mathrm{(}{P}_{\lambda }\mathrm{)}}_I$ and ${\mathrm{(}{P}_g\mathrm{)}}_I$ (Rodríguez-Tembleque and Abascal, 2010). To find the iteration direction $\Delta z^{(n)}$, this simplification leads to the following equation to be solved:

$\Theta (z^{(n)})+R^{(n)}\Delta z^{(n)}=\mathrm{0,}$ (56)

which is the same as solving the following SLE,

$R^{(n)}(z^{(n)}+\Delta z^{(n)})-f=0.$ (57)

7 COMPLEX STEP METHOD

Mundstock and Marczak (2009) obtained good results calculating shape derivatives for 2D elasticity problems using the complex step method with the BEM. In order to calculate the derivative of a function $f\mathrm{(}x\mathrm{)}$, one has to take the complex part of the solution and divide it by the perturbation $h$, as in

$\frac{\partial f\mathrm{(}x\mathrm{)}}{\partial x}\mathrm{=}\frac{\text{Im}\mathrm{[}f\mathrm{(}x+i\mathrm{.}h\mathrm{)]}}{h}+O\mathrm{(}{h}^2\mathrm{).}$ (58)

The procedure can be extended to evaluate any kind of derivative in existing numerical codes by adding the complex increment to the desired variable, as long as $f\mathrm{(}x\mathrm{)}$ is a real function.

8 RESULTS

The problem analyzed in this work consists of two spheres in contact (Fig. 3a), allowing the direct comparison with the classical Hertz solution. Both solids were considered elastic which is not a case commonly seen in the literature - generally, one of the bodies is considered rigid and the other a half-space BEM region. The mesh was simplified by taking advantage of the symmetries, as shown in Fig. 3b. The discretization used allows the evaluation of the algorithm when a small region of the solid is transferring the load, resulting in significant stress concentration at a few nodes.

The mesh used for each sphere has 264 quadrilateral elements, and is illustrated in Fig. 3b for quadratic elements. The mesh has a large element size ratio in order to reduce the overall number of degrees of freedom (DOF) without compromising the solution accuracy. The linear and quadratic elements resulted in 6936 and 14400 DOF respectively. The kinematic relations and the contact restrictions accounted for 9 percent of these DOF.

Prescribed displacements ${\overline{u}}_z$ were applied to the upper face of the half sphere, while symmetry conditions were used over the corresponding planes.

Figure 3
(a) Problem description and (b) mesh used on the Hertzian contact example.

An offset of $\mathrm{12.5\%}$ was used in all the elements, and all material and geometrical properties considered are listed in Table 1. To the best of authors knowledge there is no closed form solution for this problem when $\mu \ne 0$. In order to provide a reference solution for comparison of the results, the geometrical and material properties selected were the same as in Rodríguez-Tembleque and Abascal (2010). The friction coefficient used was $\mu \mathrm{=0.1}$. Table 1 also lists the maximum contact pressure $p_0$ and the contact area radius predicted for the prescribed displacements ${\overline{u}}_z$, from Johnson (1987).

Table 1
Sphere on Sphere contact problem parameters and results from the hertz solution.

BEM model parameters
Young's Modulus (both spheres)	$E^1=E^2$	1e4	Pa
Poisson Ratio (sphere 1)	${\nu }_1$	0.3	-
Poisson Ratio (sphere 2)	${\nu }_2$	0.4	-
Coulomb friction coefficient	$\mu$	0.1	-
Radius (both spheres)	$R^1=R^2$	1.0	m
Initial Separation	$k_0$	1.0E-1	m
final stress (MPa)	${\overline{u}}_z-k_0$	8.0E-4	m
Hertz solution
Maximum contact presure (Eq.(60)	$p_0$	145.5161	Pa
Radius of the contact area	$a$	2.0E-2	m

A single load step was used to set the boundary conditions, and the initial parameters for the GNMls were $q\mathrm{=1.0}$, $\beta \mathrm{=0.90}$, $r_t\mathrm{=}{r}_n\mathrm{=0.70}$. The residue ${\Psi }^{(n)}$ and the scaling parameter ${\alpha }^{\mathrm{(}n\mathrm{)}}$ during the GNMls iterations are shown in Fig. 4. The residue during the first iterations does not change, i.e., once the algorithm finds an optimal direction, it converges with a logarithmic rate.

Figure 4
(a) GNMls convergence: Residue (${\Psi }^{(n)}$) and scaling parameter (${\alpha }^{\mathrm{(}n\mathrm{)}}$) obtained by the line search procedure for each iteration (n).

The displacements along the contact area are plotted in Fig. 5, as a function of the coordinates $x$. These values were extracted from a line of nodes close to plane of symmetry $xz$. Both numerical and analytical displacements were normalized with the prescribed displacements ${\overline{u}}_z$. In this particular example, the quadratic element (Q8) showed a larger difference to the Hertz solution for the displacements than the linear element (Q4). The tangential slip obtained for both elements is also depicted in the figure, normalized with respect to the maximum slip of the quadratic element. In this case it is possible to verify the transition from stick to slip states.

Figure 6 shows the normal contact tractions ${\lambda }_n$ normalized with the maximum analytic pressure $p_0$. The normal displacements and tractions agree well with the reference solution. The norm of tangential tractions divided by the friction coefficient is also plotted.

At approximately 75% of the normalized contact distance the tangential traction reach the maximum value predicted by the Coulomb law, precisely where the nodes initiate the slip state. Although the maximum normal traction does not present a pronounced difference between the linear and quadratic elements, at the tangential direction it is possible to see a larger difference between the two. The quadratic elements also predict a larger slip region. At this transition region the normal tractions are in much better agreement with the reference solution, what indicates a better approximation by the quadratic element. This effect is also noted in the results of Rodríguez-Tembleque and Abascal (2010), as the number of nodes is even more reduced in their mesh. The tangential tractions should have vanished at the symmetry line $x\mathrm{=0}$, but they are producing a residual value. In numerical experiments we found out that if entire meshes are used without the aid of symmetry, correct values are obtained.

Figure 5
Normalized displacements and tangential slip along possible contact region.

Figure 6
Normalized tractions on the contact region. Reference tangential tractions from Rodríguez-Tembleque and Abascal (2010).

8.1 Sensitivity analysis

The initial analysis for any finite difference scheme should be a verification of step size. However, in order to show the independence of complex step method with respect to the step size, a variation ranging from 1e-15 to 1e-300 was tested.

The complex step was added to the value of the elastic modulus used in the top sphere, and the maximum contact pressure sensitivity was analyzed. The comparison was carried out against the analytical solution, obtained by deriving the Hertz equation for pressure with respect to $E_1$.

Due to the position of the collocation points in discontinuous elements, there is no node exactly at $x\mathrm{=0}$, so the analytical solution was evaluated at the position of the closest node.

According to Johnson (1987), the normal pressure distribution along the contact region is given by:

$p\mathrm{(}r\mathrm{)=}{p}_0\left(1-\frac{r^2}{a^2}\right)\mathrm{,}$ (59)

where $p_o$ is the maximum pressure and $r$ is the distance from the center of the sphere. $p_o$ is related to the material and geometrical properties by the equation

$p_0\mathrm{=}\frac{2}{\pi }{E}^{\mathrm{*}}\sqrt{\mathrm{(}{u}_0\mathrm{/}{R}^{\mathrm{*}}\mathrm{)}}\mathrm{,}$ (60)

where $u_0$ is the prescribed displacement, $R^{\mathrm{*}}\mathrm{=(1/}{R}_1+\mathrm{1/}{R}_2{\mathrm{)}}^{-1}$ and $E^{\mathrm{*}}$ is the equivalent elastic modulus of the pair:

$E^{\mathrm{*}}\mathrm{=((1}-{\nu }_1^2\mathrm{)/}{E}_1+\mathrm{(1}-{\nu }_2^2\mathrm{)/}{E}_2{\mathrm{)}}^{-1}\mathrm{.}$ (61)

Considering ${\nu }_1\mathrm{=}{\nu }_2\mathrm{=}\nu$, the result for the derivative of the maximum pressure is

$\frac{\partial p_o}{\partial E_1}\mathrm{=}-\frac{\text{(2/}\pi \text{)}\sqrt{{\text{(u}}_{\text{o}}\text{/R)}}\left({\nu }_{\text{1}}{}^2-1\right)}{{\left(E_1\right)}^2{\left(\frac{{\nu }_{\text{1}}{}^2-1}{E_1}+\frac{{\nu }_{\text{2}}{}^2-1}{E_2}\right)}^2}\mathrm{.}$ (62)

The pressure sensitivity is therefore

$\frac{\partial p}{\partial E_1}\mathrm{=}\frac{\partial p_o}{\partial E_1}\left(1-\frac{r^2}{a^2}\right)$ (63)

Eq. (63) results, for $r\mathrm{=6.3}\times {10}^{-4}m$, a sensitivity of $7.56\times {10}^{-3}$. Table 2 shows the difference obtained by the present approach to the analytical solution, for different values of complex increments. As can be seen, the actual difference is larger than the variation of the results, revealing stability of the method for any step size chosen.

Table 2
Difference of the numerical solution to the analytic values of Eq. (63) - element Q8.

Log(h)	Difference to analytical solution ($\times {10}^{-4}$)	relative difference ($\times {10}^{-2}$)
-15	1.736 641 581 748 989 83	1.6
-25	1.736 641 581 752 129 68	1.6
-35	1.736 641 581 754 003 18	1.6
-50	1.736 641 581 749 284 74	1.6
-75	1.736 641 581 750 672 51	1.6
-100	1.736 641 581 748 105 12	1.6
-150	1.736 641 581 753 465 42	1.6
-200	1.736 641 581 746 006 11	1.6
-250	1.736 641 581 746 734 69	1.6
-300	1.736 641 580 409 436 36	1.6

Figure 7 illustrates the normal and tangential traction sensitivities over the contact area. It is worth to note traction shows the higher sensitivity, differently from Fig. 6, where they were divided by the friction coefficient.

Figure 7
Sensitivity of normal and tangent contact tractions to the variation of ${\text{E}}_1$

8.1 Computational Aspects

The code was developed under Matlab environment. All computations were performed on a Desktop PC running Linux OS. The maximum memory usage was approximately 5 Gb during assembly procedures, reducing to 1.7 Gb during the solution process, for quadratic elements. The system size was small enough to fit on the system's RAM, allowing the use of UMFPACK direct solver from Davis (2004). The total run time for the cases using complex algebra was 600 seconds, while the non complex time was around 300 seconds.

Preliminary tests suggests avoiding iterative solvers whenever possible due to their need of efficient preconditioning schemes. The iterative solvers tested (GMRES, BICG and other variants) needed a considerable number of iterations, related to the poor matrix conditioning, resulting in a longer solution time. Even with the lower operation complexity of iterative solvers, the poor floating point per second performance of the sparse matrix-vector operations, when compared with the BLAS3 routines used by the UMFPACK direct solver, resulted in a longer solution time. Investigations on iterative solvers with BEM found on Barra et al. (1992), Valente and Pina (2006), Xiao and Chen (2007), also show a large number of iterations.

9 CONCLUSIONS

This work evaluated the use of complex step method to obtain sensitivities on BEM-BEM contact problems with friction. For the particular case evaluated, the nonlinear solution obtained by the GNMls presented a fast convergence.

Despite all the simplifications considered, the results obtained for the normal direction of displacements and tractions on the contact regions have shown a relative difference to the analytical solution of $1\times {10}^{-2}$. Moreover, the tangential tractions appear to be more affected by such simplifications, as can be seen in the non vanishing tangential tractions at center of the contact. This result is also affected by the symmetry boundary condition in conjunction with the use of discontinuous elements.

A direct comparison between the results from linear and quadratic elements show that the latter was in fact accurate for contact problems, mostly due to the superior number of nodes in the contact region. The quadratic elements delivered results in much better agreement for normal tractions. The quadratic elements also predicted a larger slip area, due to the higher DOF count on the contact region.

The numerical sensitivities evaluated with present scheme also agreed well with the reference solution, making the present approach a potential candidate to be implemented in structural optimization cases containing contact.

The run times indicate that the complex variables do not affect the integration, assembly, or the solution performance more than the additional time needed to handle double data sizem, if compared to a standard BEM code.

It is also important to point out that the complex increment stands as an optimum scheme to be used in other contact problems since it does not affect the contact state when it is added to the nodal positions, and has proven to be very stable and reliable, regardless the step size. This robustness more than compensates its two-fold memory requirements and solution time.

10 ACKNOWLEDGEMENTS

The authors wish to express his thanks to CNPq for the financial support.

References

• Alart P., Curnier A. (1991) A mixed formulation for frictional contact problems prone to newton like solution methods, Comput. Methods Appl. Mech. Eng. 92(3): 353-375.
• Andersson T. (1981) The boundary element method applied to two-dimensional contact problems with friction, in ‘Boundary element methods’, Springer, pp. 239-258.
• Barra L., Coutinho A., Mansur W., Telles J. (1992) Iterative solution of bem equations by gmres algorithm, Computers & Structures 44(6): 1249 - 1253.
• Davis T. A. (2004) Algorithm 832: Umfpack v4.3-an unsymmetric-pattern multifrontal method, ACM Trans. Math. Softw. 30(2): 196-199.
• Fichera G. (1963) Sul problema elastostatico di Signorini con ambigue condizioni al contorno., Atti Accad. Naz. Lincei, VIII. Ser., Rend., Cl. Sci. Fis. Mat. Nat. 34: 138-142.
• Fichera G. (1973) Boundary value problems of elasticity with unilateral constraints, in ‘Linear Theories of Elasticity and Thermoelasticity’, Springer, pp. 391-424.
• Gakwaya A., Lambert D., Cardou A. (1992) A boundary element and mathematical programming approach for frictional contact problems, Computers & Structures 42(3): 341 - 353.
• Gao X., Liu D., Chen P. (2002) Internal stresses in inelastic BEM using complex-variable differentiation, Computational Mechanics 28(1): 40-46.
• Garrido J., Foces A., Paris F. (1991) BEM applied to receding contact problems with friction, Mathematical and Computer Modelling 15(3): 143-153.
• Garrido J., Lorenzana A. (1998) Receding contact problem involving large displacements using the BEM, Engineering Analysis with Boundary Elements 21(4): 295 - 303. Contact Mechanics.
• Ghaderi-Panah A., Fenner R. (1998) A general boundary element method approach to the solution of three-dimensional frictionless contact problems, Engineering Analysis with Boundary Elements 21(4): 305 - 316.
• González J. A., Park K., Felippa C. A., Abascal R. (2008) A formulation based on localized lagrange multipliers for bem-fem coupling in contact problems, Computer Methods in Applied Mechanics and Engineering 197(6):623-640.
• Johnson K. L. (1987) Contact mechanics, Cambridge university press.
• Lyness J. N., Moler C. B. (1967) Numerical differentiation of analytic functions, SIAM Journal on Numerical Analysis 4(2):202-210.
• Man K., Aliabadi M., Rooke D. (1993a) BEM frictional contact analysis: load incremental technique, Computers & structures 47(6):893-905.
• Man K., Aliabadi M., Rooke D. (1993b) BEM frictional contact analysis: modeling considerations, Engineering analysis with boundary elements 11(1):77-85.
• Martins J. R., Sturdza P., Alonso J. J. (2003) The complex-step derivative approximation, ACM Transactions on Mathematical Software (TOMS) 29(3):245-262.
• Mundstock D. C., Marczak R. J. (2009) Boundary element sensitivity evaluation for elasticity problems using complex variable method, Structural and Multidisciplinary Optimization 38(4):423-428.
• Pang J.-S. (1990) Newton’s method for b-differentiable equations, Mathematics of Operations Research 15(2):311-341.
• Paris F., Blazquez A., Canas J. (1995) Contact problems with nonconforming discretizations using boundary element method, Computers & structures 57(5): 829-839.
• Paris F., Garrido J. (1989) An incremental procedure for friction contact problems with the boundary element method, Engineering Analysis with Boundary Elements 6(4): 202 - 213.
• Rodríguez-Tembleque L., Abascal R. (2010) A FEM-BEM fast methodology for 3d frictional contact problems, Computers and Structures. 88(15-16): 924-937.
• Rodriguez-Tembleque L., Buroni F., Abascal R., Sáez A. (2011) 3d frictional contact of anisotropic solids using bem, European Journal of Mechanics-A/Solids 30(2): 95-104.
• Rodríguez-Tembleque L., González J. A., Abascal R. (2008) A formulation based on the localized Lagrange multipliers for solving 3D frictional contact problems using the BEM, Numerical Modeling of Coupled Phenomena in Science and Engineering: Practical Use and Examples p. 359.
• Squire W., Trapp G. (1998) Using complex variables to estimate derivatives of real functions, SIAM Review 40(1):110-112.
• Strömberg N. (1997) An augmented Lagrangian method for fretting problems, European Journal of Mechanics-A/Solids 16, 573-593.
• Valente F., Pina H. (2006) Conjugate gradient methods for three-dimensional BEM systems of equations, Engineering Analysis with Boundary Elements 30(6):441 - 449.
• Xiao H., Chen Z. (2007) Numerical experiments of preconditioned krylov subspace methods solving the dense non-symmetric systems arising from bem, Engineering Analysis with Boundary Elements 31(12): 1013 - 1023.
• Yamazaki K., Sakamoto J., Takumi S. (1994) Penalty method for three-dimensional elastic contact problems by boundary element method, Computers & structures 52(5): 895-903.

• Available Online: May 28, 2018

Publication Dates

• Publication in this collection
  2018

History

• Received
  01 Aug 2017
• Reviewed
  24 May 2018
• Accepted
  24 May 2018