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

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.


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 (1963Fichera ( , 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 where the relative error for the derivative value is very stable and close to the mantissa for increments of 9 1 10   to 300 1 10   .
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.

FORMULATION OF THE CONTACT PROBLEM
Let the problem of two elastic bodies 1 The possible contact boundary have conditions which depends on the contact state, defined by the gap 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.

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 2  ) are positioned as closely as possible to the master (region 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  The gap g for any pair of points in this problem is obtained through the following relation where 1 x and 2 x 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., where n is the outward normal vector at 1 x and 1 2 , t t are the tangential vectors. Therefore the gap can be decomposed in its tangential and normal components: 2.2 Unilateral and frictional contact conditions The unilateral and frictional laws for a boundary point are grouped as suggested in , resulting in the well-known contact conditions which depends on the normal gap n g , normal ( ) . ; Contact -Slip case, 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 t The displacement boundary integral equation for a point in the boundary is where the boundary is always smooth due to the use of discontinuous elements, i.e., 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 where N is the number of nodes of the element and n  is a vector containing the shape functions evaluated at  .
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 : 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 where 1 2 ( , ) n n   correspond to the nodal coordinates of the n-th geometrical node, 1 2 ( , )   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: which reduces to the continuous functions n  when a =0. The interpolation functions for the 8-node element are also written in compact form, Collocation for node i , generates the discretized form of Eq. (6): . e e e e i i ij ij Performing the collocation process over all boundary nodes results in the algebraic system of equations that leads to the BEM solution, i.e.,

CONTACT PROJECTION FUNCTIONS
The formulation used for the imposition of contact restrictions was implemented according to . The handling of inequalities, typical of contact conditions, was avoided using suitable projection functions.

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 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 projecting the variable v onto   , the admissible region of the contact normal tractions, i.e.
A mixed variable is now introduced, the augmented normal traction Similarly, the augmented tangential traction vector is defined as Considering a plane defined by the tangential tractions,  represents admissible for the tangential tractions under a given normal traction n t . The tangential contact restrictions are finally written as

Normal-Tangential operator
In order to combine the normal and tangential restrictions in a single statement, the normal-tangential projection function where the region f  is the augmented friction cone, with radius

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, where the compatibility conditions 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.

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 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 where  ji  C 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 j k in Eq. (26) and = ( ) m I  , we have: 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 where k is the gap vector accounting for all the contact pairs. The vector go k is the sum of the initial gap and rigid body displacements vectors. Substituting Eq. (26) and (27) where  are the BEM nodal tractions rotated from the global to the master nodes local coordinate system.

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 Substituting *   rk λ λ with r being the penalization parameters vector, the previous equation becomes For the normal direction one writes the restriction as  (20)), the tangential restriction results which, using  As for the slip case, • Slip case:

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  -differentiable functions are described as non Fréchet-differentiable due mainly to the absence of linearity on the  -derivative. Pang x h x f x , 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 ( ) n , the GNMls algorithm can be summarized in the following steps: z be an arbitrary initial vector, and , the direction ( ) n z is obtained by solving the equation: 3. Find the first integer = 1,2, m  which fulfills the decreasing error condition: 4. Updated solution vector: , stop. Otherwise = 1 n n  , and return to step 2.

Search direction: Directional derivatives
The system of equations (46) can be split in its linear and nonlinear parts, is the linear part of the Jacobian matrix, is the nonlinear part of the Jacobian. A linear approximation for this derivative can be used (Strömberg, 1997): , and (Rodríguez-Tembleque and Abascal, 2010). To find the iteration direction ( ) n Δ z , this simplification leads to the following equation to be solved: which is the same as solving the following SLE, 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 x , one has to take the complex part of the solution and divide it by the perturbation h, as in 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 x is a real function.

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 z u were applied to the upper face of the half sphere, while symmetry conditions were used over the corresponding planes. An offset of 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 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 = 0.1  . Table 1 also lists the maximum contact pressure 0 p and the contact area radius predicted for the prescribed displacements z u , from Johnson (1987).   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 x z . Both numerical and analytical displacements were normalized with the prescribed displacements z u . 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 n  normalized with the maximum analytic pressure 0 p . 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 = 0 x , 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.

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 1 E .
Due to the position of the collocation points in discontinuous elements, there is no node exactly at = 0 x , 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: where o p is the maximum pressure and r is the distance from the center of the sphere. o p is related to the material and geometrical properties by the equation where 0 u is the prescribed displacement, * 1 1 2 and * E is the equivalent elastic modulus of the pair: The pressure sensitivity is therefore Eq. (63) results, for 4 = 6.3 10 r m   , a sensitivity of 3 7.56 10   . 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.  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. 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.

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 2 1 10   . 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.