Continuum framework for finite element modelling of finite wear

Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Continuum framework for finite element modelling of finite wear Jakub Lengiewicz, Stanisław Stupkiewicz ⇑ ´ skiego 5B, 02-106 Warsaw, Poland Institute of Fundamental Technological Research (IPPT), Pawin

a r t i c l e

i n f o

Article history: Received 6 May 2010 Received in revised form 23 November 2010 Accepted 21 December 2010 Available online 30 December 2010 Keywords: Frictional contact Wear Sensitivity analysis

a b s t r a c t A finite deformation contact problem with friction and wear is studied in which the shape changes due to wear are finite. Accordingly, in addition to the initial configuration and the current configuration, an intermediate time-dependent configuration is introduced that corresponds to the undeformed body of the shape changed due to wear. Two time scales are also introduced in order to distinguish the fast time of the actual deformation (contact) problem from the slow time of the wear process (shape evolution problem). Separation of these time scales allows us to partially decouple the deformation problem and the shape evolution problem. Shape parameterization is introduced and the corresponding shape update scheme is formulated as a minimization problem. In particular, a second-order scheme is developed which exploits shape sensitivities of the deformation problem. Numerical examples are provided to illustrate the performance and accuracy of the proposed numerical schemes. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Wear is a process in which material is removed from a contact surface as a result of damage processes induced at the contact interface during frictional contact interaction. There are many mechanisms responsible for wear, and their activity and importance heavily depend on the materials and surface properties of the contact pair, as well as on the actual contact conditions, e.g. [1]. Accordingly, formulation of a general predictive wear law, and even selection of the relevant influential variables, is difficult, if not impossible [2]. The related aspects of constitutive modelling of wear are out of scope of this work, and the classical wear law of Archard [3] is adopted. In this work we address another aspect of modelling of wear, namely how to consistently formulate and efficiently solve frictional contact problems with wear. Specifically, we focus on the finite changes of shape due to wear, and a possibly general contact problem with wear is considered in which both the deformation and the shape changes due to wear are finite. The formulation covers thus contact and wear of materials, such as polymers or rubberlike materials, but also more usual metallic or ceramic contacts, as a special case. In general, two classes of contact problems with wear can be distinguished. In the problems belonging to the first class, wear is considered as a postprocessing quantity that is evaluated after the solution of the contact problem is obtained, e.g. [4]. The effect

⇑ Corresponding author. E-mail address: [email protected] (S. Stupkiewicz). 0045-7825/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2010.12.020

of the accumulated wear (and of the related change of shape) on the solution of the contact problem is thus neglected. In the second class of problems, the coupling of deformation and wear is fully accounted for, i.e. the contact interaction is affected by the time-dependent change of shape of contacting bodies due to the accumulated wear. In practical terms, the problems of the second type are usually solved in an incremental manner, where the problem of the first type is solved at each increment, and the shape is updated using the wear profile accumulated over that increment (typically with a time speedup as the wear processes are usually very slow). This general incremental scheme is at the basis of the many available wear simulation approaches, e.g. [5–10], and it is also adopted in this work. In the context of the finite element method, the incremental shape update is associated with remeshing or remapping of the finite element mesh according to the computed incremental wear profile. This can be avoided if no configuration changes are assumed and linear kinematics is adopted. Then the accumulated wear can be added to the initial normal gap, and wear increments can be computed at each time increment of the frictional contact problem, e.g. [11,12]. The latter formulation is also beneficial because the implicit backward-Euler scheme can be quite easily applied for the time integration of wear [12]. An implicit time integration scheme has also been applied in [13] in combination with the non-smooth contact dynamics method. Otherwise, the remeshing-based techniques are usually combined with the explicit forward-Euler scheme which is conditionally stable, hence problems related to the instability of the time integration scheme are necessarily encountered, e.g. [11,5–7]. Let us also mention here alternative approaches in which the steady-state or asymptotic

179

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

wear states are determined directly rather than incrementally, e.g. [14,15]. In this work, we focus on the finite changes of shape (or configuration) due to wear, and the corresponding continuum framework is proposed along with its finite element treatment. Apparently, the first trial to develop a continuum formulation including finite wear is due to Jourdan [16]. However, that formulation is questionable because an additive split of the transformation mapping into elastic and wear parts is assumed by analogy to the additive split postulated in [17] in the small-strain format. Accordingly, an additive split of the deformation gradient is obtained in [16] and the elastic deformation gradient is evaluated with respect to the initial (unworn) configuration. In the present work, we consider three configurations, and the shape evolution mapping is introduced which transforms the initial configuration to the undeformed configuration, as well as the deformation mapping which transforms the undeformed configuration to the current configuration. Furthermore, two time scales are introduced in order to distinguish the fast time of the contact problem from the slow time of the wear process. These notions constitute the basis of the continuum formulation presented in Section 2. Next, a shape parameterization is introduced in Section 3, which is independent from the spatial discretization used for numerical solution of the contact problem. First- and second-order shape update schemes are also proposed, the latter employing shape sensitivity analysis. The finite element treatment is discussed in Section 4 and illustrative numerical examples are provided in Section 5, where the performance and the accuracy of the proposed numerical schemes are compared. The theoretical framework, the time-dependent shape parameterization and the second-order shape update scheme constitute the major new contribution of this work.

The time evolution of the undeformed configuration Xi results _ i P 0 the nominal wear rate, i.e., the volfrom wear. Denote by W ume removed from the surface due to wear per unit area and unit time, where both the volume and the area refer to the undeformed configuration Xi, see also Section 2.4. Accordingly, the following relationship, the shape evolution law, exists between the wear rate _ i and the time derivative of the shape evolution mapping Wi, W

( _ i  Ni ¼ W

_i W

on Cic ;

0

on Ci n Cic ;

where Ci is the boundary of Xi , Cic  Ci denotes the potential contact surface, and Ni is the unit outer normal defined in the undeformed configuration Xi. Note that Eq. (2) is written in the _ i denotes here undeformed configuration Xi, i.e. for Xi 2 Ci, hence W ^ i ðX i Þ; tÞ, in view of the one-to-one correspondence be_i¼W _ i ðX W ^ i imposed by the shape evolution mapping Wi. In ortween Xi and X der to simplify the notation, this convention is used whenever convenient. 2.2. Two time scales In a typical situation, the wear process is very slow at the time scale of the deformation problem, so that the rate of change of configuration Xi due to wear is negligible at this scale (e.g., compared to the sliding velocity). At the same time, it is the accumulation of wear over a long period of time that results in a significant variation of the configuration and contact conditions. Following this line, we introduce two time scales, namely the sscale of the deformation problem and the t-scale of the shape changes due to wear. The shape and deformation mappings (1) are rewritten accordingly,

^ i ; tÞ; X i ¼ Wi ðX

2. Continuum formulation

ð2Þ

xi ¼ uit ðX i ; sÞ;

t 2 ½0; T;

s 2 ½t; t þ Ds;

ð3Þ

i

2.1. Three configurations Let us consider two hyperelastic bodies Bi ; i ¼ 1; 2, in frictional contact. The bodies undergo finite deformation and finite shape changes due to wear at the contact interface. Thus, we introduce three configurations of each of the bodies Bi : the initial configura^ i , the time-dependent undeformed configuration Xi, and the tion X current (deformed) configuration xi, cf. Fig. 1. The shape change due to wear is described by the shape evolution mapping Wi, and the deformation is described by the deformation mapping ui,

^ i ; tÞ; X i ¼ Wi ð X

xi ¼ ui ðX i ; tÞ;

t 2 ½0; T;

ð1Þ

^i 2 X ^ i ; X i 2 Xi ; xi 2 xi . If the initial configuration X ^ i is a where X stress-free (natural) configuration then the undeformed configuration Xi is also stress-free. The deformation mappings ui are governed by the equilibrium equations, constitutive relations, and boundary and contact conditions (the respective formulation is discussed later).

^ i ; tÞ, and Ds is a characteristic or represenwhere X 2 Xit ; Xit ¼ Wi ðX tative time of the deformation problem, for instance, one cycle of a cyclic loading program. Now, the separation of the two time scales is assumed, so that the deformation problem can be analyzed at fixed t, and it becomes a standard frictional contact problem. The corresponding formulation is briefly introduced in Appendix A. Solution of the deformation problem yields the deformation mappings u1t ðX 1 ; sÞ and u2t ðX 2 ; sÞ, as well as the contact variables such as the sliding velocity vT and the contact tractions tN and tT. _ i ðX ^ i ; sÞ is determined for all points Furthermore, the wear rate W t i i ^i X ¼ W ðX ; tÞ on the potential contact surface Cic using a suitable wear law, for instance, the Archard law discussed in Section 2.5. _ i , the wear increment DW i ðX ^ i Þ accumulated durBy integrating W t t ing the time period [t, t + Ds] is then obtained, which is used to de_ i ðX ^ i Þ, the average wear rate at the t-scale, viz. fine W t i

^ _ i ðX ^ i Þ ¼ DW t ð X Þ ; W t Ds i

^ iÞ ¼ DW it ðX

Fig. 1. The three configurations at two time instants (shown for one body only).

Z t

tþDs

_ i ðX ^ i ; sÞds; W t

^i 2 C ^i : X c ð4Þ

180

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

The assumption of separation of the two time scales requires that the shape change associated with the wear increment DW it accumulated over one deformation cycle is sufficiently small so that its effect on the solution of the deformation problem is negligible. As a result, the wear rate at the t-scale can be determined as a postprocessing quantity at the s-scale, and the problem of shape evolution due to wear and the deformation problem are partially decoupled. Such decoupling is, in fact, implicitly assumed in the simulation approaches developed, e.g., in [5–10].

nodes are determined. Alternative shape update schemes are discussed in Section 3.

Remark 1. Specification of the deformation problem at the sscale may be nontrivial in general. This is because friction is a path-dependent phenomenon. Suitable boundary and initial conditions must thus be applied in order to properly describe the complex evolution of stick and slip zones while the shape evolves due to wear. However, in some specific situations, the related problems are easily overcome, for instance, when each loading cycle starts with an open contact (as assumed, e.g., in [9,10]) or when gross sliding occurs during the loading cycle, so that the contact memory is erased (as in the numerical example discussed in Section 5.2).

_ i ds ¼ dV ¼ 1 dm ; W i .i0 dSi dS

2.4. Nominal and spacial wear rate _ i , as introduced in Section 2.1, refers to The nominal wear rate W the undeformed configuration Xi and is expressed in terms of the volume of the material removed due to wear. It can be equivalently expressed in terms of the mass, namely i

i

where .i0 is the mass density in the undeformed configuration. By _ i can be defined, analogy, the spacial wear rate w

_ i ds ¼ w

dv

i

i

ds

i

¼

2.3. Time integration of the shape evolution problem In the time-discretized setting, the shape transformation Wi is only evaluated at discrete time instants. This will be denoted by ^ i Þ ¼ Wi ðX ^ i ; tÞ. Application of the forward-Euler time integraWit ðX tion scheme to the shape evolution law (2) yields

h

i

_ i ðX ^ i Þ  Wi ðX ^ i Þ  N i ðX ^ i Þ ¼ DtW ^ i Þ; WitþDt ðX t t t

X it 2 Cic;t ;

. ds i

i

;

ð8Þ

i i _ i; _ ¼ Ji W jw

ð9Þ i

i

where j = ds /dS is the area transformation factor (which follows from the Nanson’s formula, n ds = JFTN dS), and Ji = dvi/dVi = det Fi is the determinant of the deformation gradient Fi = @xi/@Xi. 2.5. Archard wear law The classical wear law of Archard [3] is adopted in this work, however, other wear laws can equally be used. The Archard law states that the wear rate is proportional to the normal pressure and sliding distance. Assuming that the Coulomb friction law holds, the normal contact pressure is proportional to the sliding friction stress, hence the Archard law can be equivalently expressed in terms of the latter, or, more generally, in terms of the frictional dissipation rate density D_ i , cf. [18–20], namely

_ i ¼ K i D_ i ; W

i i D_ i ¼ j d_ ¼ j t T  v T ;

ð10Þ

i

ð5Þ

_ i depend on the position X ^ i in the initial configurawhere N it and W t tion through the shape evolution mapping Wit . For the computational efficiency reasons, the time increment Dt can be adopted much larger than Ds (but limited by the desired accuracy of the time integration scheme). As the forward-Euler scheme is only conditionally stable, the maximum time increment Dt is constrained by the corresponding stability limit. ^ iÞ The incremental form (5) constitutes a constraint on WitþDt ðX as it prescribes its normal component on the boundary. In order ^ i Þ, additional assumptions must be to fully determine WitþDt ðX adopted concerning the tangential part and the shape transformation in the interior. For instance, one can assume that the shape transformation is such that the points on the contact boundary Cic are transformed along the normal direction. This yields the following shape update scheme

X itþDt ¼ X it  M DW it ðX it ÞN it ;

1 dm

where .i is the mass density in the current configuration xi, and we have

i

Remark 2. A class of quasi-steady-state processes can be defined in which the deformation problem is a steady-state process independent of the time s. An example is the pin-on-disk tribological test, in which the reference frame is attached to the pin, and the disk is analyzed in an Eulerian frame. In the case of a quasisteady-state process, the average wear rate at the t-scale is simply _ i ðX ^ iÞ ¼ W _ i ðX ^ i ; sÞ. equal to the actual wear rate at the s-scale, W t t

ð7Þ

ð6Þ

and X itþDt ¼ X it for X it 2 Cit n Cic;t . Here, M P 1 is a time-integration _ i ¼ M DW i , speedup parameter, Dt = MDs, so that we have DtW t t i and DW t is defined by Eq. (4). The shape update scheme (6) is, in fact, a commonly used scheme that is applied at the contact nodes of the finite element mesh, e.g. [5–10]. The shape update (6) of the boundary nodes is then followed by remeshing in which the positions of the interior

where K is the wear coefficient, tT is the spatial friction stress, and vT is the sliding velocity. Note that the wear law (10) defines the nominal wear rate, so that, consistently, the dissipation rate density D_ i is referred to the area in the undeformed configuration Xi, hence the area transformation factor in Eq. (10)2. At the same time, the spatial density of frictional dissipation rate d_ ¼ t T  v T is referred to the area in the current configuration xi. Formally, one could state the wear law by postulating proportionality between the spacial quantities, namely

_ _ i ¼ ki d: w

ð11Þ

In view of Eq. (9), the corresponding nominal form would then be

_ i ¼ 1 ki D_ i ; W Ji

ð12Þ

which is not equivalent to Eq. (10) as it involves Ji. Eqs. (10) and (12) define a transformation relationship Ki = ki/Ji between the wear coefficients Ki and ki, so that if one of them is a constant then the other depends on the deformation in the subsurface layer. Obviously, the difference between the two forms vanishes when the materials are incompressible or when the deformation is small. We are not aware of any experimental observations that would justify the choice of one or the other form of the wear law. In the following, we choose the nominal form (10) of the wear law because, in practice, the wear profile would rather be measured in the undeformed configuration, hence the nominal wear rate can be considered an observable quantity. Secondly, the computational

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

treatment of the spatial form (12) would be more difficult because _ i , used in the computational framework, the nominal wear rate W would then depend not only on the contact quantities, but also on the deformation in the bulk material (through Ji). Note that there is an alternative description in which wear criterion is formulated in terms of the thermodynamic driving force for the propagation of a damage interface within the contact subsurface layer, cf. [21,14,22]. In this case, evaluation of the wear criterion would also involve the stresses or strains in the bulk material. 3. Shape parameterization and shape update schemes 3.1. Shape parameterization Let us introduce an approximation of the shape evolution mapping Wi expressed in terms of time-dependent shape parameters assembled in the vector /(t), namely

  ^ i ; /ðtÞ : ^ i ; tÞ  W ei X X i ¼ Wi ð X

ð13Þ

The corresponding deformation mapping also depends on the shape parameters /,

~ i ðX i ; s; /ðtÞÞ: xi ¼ uit ðX i ; sÞ  u

ð14Þ

This dependence is, however, implicit. The governing equations of the deformation problem are now enforced on the domain Xit , which is defined by the shape parameterization (13), i.e. ^ i ; /ðtÞÞ. e i ðX Xit ¼ W Note that the shape parameterization (13) is defined separately for each contacting body with the corresponding vector of shape parameters /i. However, the shape parameters /1 and /2 of both ~ 1 and u ~ 2, bodies do affect each of the deformation mappings u hence the global vector of shape parameters / = {/1, /2} is introduced in Eqs. (13) and (14). Note also that the shape parameterization (13) is in general independent of the discretization (e.g. finite element discretization) used to solve the deformation problem. The shape evolution law (2) cannot be now fulfilled locally because of the finite number of shape parameters in the shape parameterization (13). Accordingly, the shape evolution law (2) is enforced approximately through the following minimization problem:

!2

min D/

181

2 Z  h i 2 X ^ i ; / þ D/Þ  W ^ i ; / Þ þ DW i ^ji d^Si ; e i ðX e i ðX N it  W t t t t i¼1

^i C

ð17Þ where

_ i ¼ M DW i ; DW it ¼ Dt W t t

ð18Þ

DW it

is computed from Eq. (4) after the deformation problem corresponding to / = /t is solved, and Dt = MDs, as in Eq. (6). Note that N it ; DW it and ^jit are here evaluated at the beginning of the time increment, i.e. for / = /t. We also note that the minimization in Eq. (17) can be performed separately for each body. Upon solving the minimization problem (17) for the unknown increment D/, the shape parameters are updated according to

/tþDt ¼ /t þ D/;

ð19Þ

and the solution proceeds to the next time step. Remark 3. The nodal shape parameterization is a special case of the general shape parameterization introduced above, in which the number of shape parameters is equal to the number of nodes on the discretized contact surface, and the positions of the surface nodes (e.g., along the normal direction) constitute the shape parameters. Then the first-order shape update scheme (17) and (18) is essentially equivalent to the usual forward-Euler wear integration scheme, as used in [5–10], except that the shape parameterization (13) directly defines the mesh update (reme^ shing) within the whole domain Xi , see also Eq. (6) and the following discussion. 3.3. Second-order approximation of wear increment DW it Consider now the following second-order approximation of the wear increment,

_i € i ¼ Dt W _ i þ 1 Dt 2 W _ i þ 1 Dt 2 @ W t /_ DW it ¼ Dt W t t t 2 2 @/ ! i 1 @ DW t ¼ M DW it þ D/ : 2 @/

ð20Þ

3.2. First-order forward-Euler time integration scheme

Here, the second term of the expansion accounts for the variation of the wear rate caused by the shape evolution during the time increment Dt. In order to evaluate this term, the shape sensitivity @ DW it =@/ of the wear increment DW it must be computed at the fast time scale. Details of shape sensitivity analysis are not provided here. The reader is referred to [23,24] for a general overview of analytical sensitivity analysis methods and to [25–27] for the details of the present implementation and automation of the direct differentiation method in the context of the frictional contact problems. The shape update scheme based on the minimization problem (17) with DW it given by Eq. (20) will be referred to as the second-order scheme. The performance of this method, compared to the first-order scheme using DW it given by Eq. (18), is studied in Section 5. Note that the variation of shape parameters of one body affects the wear increment of the other body. Thus the sensitivity @ DW it =@/ introduces coupling of shape evolution of both bodies so that the minimization problem (17) must now be solved simultaneously for both contacting bodies.

Application of the forward-Euler time integration scheme to the minimization problem (16) yields the following shape update scheme

Remark 4. The framework introduced in this work relies on the assumption of hyperelastic material behaviour, so that the changes due to wear are fully defined by the shape evolution mappings Wi.

min /_

2 Z X i¼1

Cit

N it

ei @W _i  /_ þ W t @/

i

dS ;

ð15Þ

_ i ¼ ð@ W e i =@/Þ/_ and, in order to simplify the notation, we dewhere W _ i ¼ 0 on Ci n Ci . The integral in Eq. (15) can be transformed fine W t c;t t ^ i , and the equivalent minimization to the initial configuration X problem reads

min /_

2 Z X i¼1

^i C

N it

ei @W _i  /_ þ W t @/

!2 ^ji d^Si ; t

ð16Þ

i where ^jit is the area transformation factor, dS ¼ ^jit d^Si , which depends on the shape transformation mapping Wi through the Nanson’s formula. The minimization problem (16) uniquely defines the shape evolution due to wear in the whole domain. This is in contrast to the general shape evolution law (2) which enforces only a constraint on the shape evolution mapping Wi.

182

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

Direct generalization to elastoplasticity or other path-dependent material behaviour is possible if no time-integration speedup is employed, i.e. for M = 1. Then the shape update schemes (6) or (17) with the corresponding remeshing of the bodies should be accompanied by an appropriate mapping of the internal variables (e.g., plastic strains) from the old mesh to the new one. Otherwise, for M > 1, development of an adequate description of the evolution of internal variables on the slow time scale t is not straightforward.

C1-continuity of the interpolation at x22 and x23 is enforced, cf. [30,29]. The Bézier curve is thus given by 3 X

x2 ðnÞ ¼

Bk3 ðnÞdk ;

Bkn ¼

k¼0

1 2n

  n ðn þ 1Þk ðn  1Þnk ; 1 6 n 6 1; k ð22Þ

Bkn

where is the Bernstein polynomial of degree n and the control points are defined by

d0 ¼ x22 ;

d1 ¼ x22 þ aðx23  x21 Þ;

d2 ¼ x23 þ aðx22  x24 Þ;

ð23Þ

4. Finite element discretization and implementation 4.1. Node-to-segment contact element In the node-to-segment discretization strategy, the contact conditions, i.e. the unilateral contact condition (A.3) and the Coulomb friction law (A.4), are enforced at the nodes of the slave surface, while the master surface is represented by segments (see Appendix A for the continuum formulation of frictional contact). Each slave node forms a contact element, and the global contact search procedure is used to find the master segment associated with the slave node at the current time step. Upon the finite element discretization, the contact contribution Gc to the virtual work, Eq. (A.7), is approximated by the summation over the contact elements,

Gc ðu; duÞ ¼

Z C1c

ðT N dg N þ T T a dna ÞdS 

X

A1e ðT N;e dg N;e þ T T a;e dnae Þ;

e2S act c

ð21Þ S act c

d3 ¼ x23 ;

A1e

where denotes the set of active contact elements and is the tributary area associated with the slave node. An active contact element is thus formed by the slave node X 1e , the master nodes X 2k;e defining the master segment, and auxiliary nodes X 1;a k;e on the slave surface that are used to determine the tributary area A1e (here, the nodes are denoted by their positions in the respective undeformed configurations). In the following, we shall omit the contact element index e. More details of the finite element implementation are provided below for the case of the two-dimensional smooth contact element shown in Fig. 2 (the nodes are shown in the current configuration in which smoothing and local contact search are performed). Here we assume that both bodies are discretized into quadrilateral fournode elements, so that the discretized boundaries consist of linear segments. Smoothing of the master surface is introduced in order to improve the convergence behaviour and the accuracy of the finite element solution [28,29]. In the present formulation, cubic splines are used for that purpose. The control points dk defining the Bézier curve joining the master nodes x22 and x23 are defined in terms of the additional two neighboring nodes x21 and x24 so that

and a = 1/6. Finally, the tributary area associated with the slave node X1 is determined in the undeformed configuration and is given by, cf. [31,32],

A1 ¼

 1 1 1 1;a kX  X 1;a 1 k þ kX  X 2 k : 2

ð24Þ

4.2. Wear The wear increments evaluated at the slave node during one time increment are obtained by backward-Euler integration of the rate equations (A.9) and (A.10),

DW 1 ¼ K 1 DD1 ;

DW 2 ¼ K 2 DD1 ;

DD1 ¼ T T a Dna ;

ð25Þ

DW 2

where denotes the wear volume of the master surface per unit area of the undeformed slave surface, cf. Eq. (A.10). In order to consistently transfer the wear increment DW 2 to the nodes of the master surface, we first compute the corresponding wear volume DV2,

DV 2 ¼ A1 DW 2 ;

ð26Þ (X 22

which is next split between the two nodes and master segment associated with the slave node,

1 DV 22 ¼ ð1  nÞDV 2 ; 2

1 DV 23 ¼ ð1 þ nÞDV 2 ; 2

X 23 )

defining the

ð27Þ

where 1 6  n 6 1 is the local coordinate of the projection point 2 ¼ x2 ð x nÞ, so that DV 2 ¼ DV 22 þ DV 23 . Finally, the wear increments at the master nodes are obtained by dividing the wear volume DV 2k by the corresponding tributary area A2k , namely

DW 2k ¼

DV 2k A2k

¼

1 A1 ð1  nÞ 2 DW 2 ; 2 Ak

ð28Þ

where the tributary area A2k corresponds to the undeformed configuration of the master body and is determined analogously to the tributary area A1 in Eq. (24),

A2k ¼

 1 2 kX k  X 2k1 k þ kX 2k  X 2kþ1 k : 2

ð29Þ

Note that DW 2k is the contribution from one contact element, thus the total wear increment at a master node is a sum of all contributions resulting from all slave nodes associated with the master segments formed by this master node. Once the deformation problem is solved at the current instant of the slow time t, the wear increments DW it can be determined according to Eq. (18) or (20) with DW it given by Eq. (4). In order to apply the shape update scheme (17), the normal vectors N it must be known at all the nodes of the slave and master surfaces. The normal at the typical node X ik is determined as the average of the normals of the two linear segments adjacent to the node, according to

N ik;t ¼ Fig. 2. Two-dimensional smooth contact element (in the current configuration).

X ikþ1  X ik1 kX ikþ1  X ik1 k

e3 ;

ð30Þ

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

where X ik1 and X ikþ1 are the positions of the two neighboring nodes, e3 is the unit vector normal to the plane of analysis, and denotes the vector product. 5. Numerical examples 5.1. Steady-state pin-on-flat problem As the first example, a highly simplified two-dimensional pinon-flat problem is studied, which roughly corresponds to the pinon-disk tribological test. A hyperelastic pin is pressed into an hyperelastic block with a constant force P, cf. Fig. 3(a). The problem is analyzed in two dimensions, and plane strain conditions are assumed, as if both the contacting bodies were long in the e3-direction, perpendicular to the plane of analysis. Relative sliding in the e3-direction is assumed, and the effect of the associated friction stresses on the in-plane deformation is neglected. As a result, the deformation problem resulting from the assumption of separation of time scales (Section 2.2) is a steady-state frictionless contact problem, and the wear rate is assumed proportional to the (nominal) normal traction TN and to the sliding velocity v3, so that we _ i ¼ ðji =j1 ÞK i T N v 3 ; i ¼ 1; 2, cf. Eqs. (A.9) and (A.10). have W The geometry of the problem is specified by the pin radius R = 3 mm, and the height and length of the block, respectively, H = 10 mm and L = 15 mm. The hyperelastic neo-Hookean material model is adopted for both bodies with identical elastic properties: Young’s modulus E = 10 MPa and Poisson’s ratio m = 0.4 (quadrilateral F-bar elements [33] are used to avoid volumetric locking effects). The loading force (per unit length in the out-of-plane direction) is equal to P = 20 N/mm. The wear coefficients at both surfaces are K1 = K2 = 105 mm2/N, the sliding velocity is v3 = 10 mm/s, and the wear process is simulated for t 2 [0, T], where T = 9600 s. The computations have been carried out in the AceGen/AceFEM environment [34]. The adopted shape parameterization is sketched in Fig. 3(b). The shape of the boundary is described by a B-spline curve, and the positions of its control points constitute the shape parameters. In Fig. 3(b), the open circles denote the positions of the control points in the initial configuration. The modified shape is defined by the control points marked by the full circles, which are obtained by scaling the initial positions (along the dotted lines) by the values of the corresponding shape parameters. For a smooth parame-

183

terization, quadratic B-splines are used, while a piecewise linear function (linear B-spline) is used to define the nodal parameterization, cf. Remark 3. In both cases, the parameterization preserves the symmetry with respect to the vertical axis. The mesh update scheme based on the above shape parameterization is analogous. The positions of the nodes of the block are scaled along the vertical lines, while the pin nodes are scaled along the family of lines passing through point A, cf. the dotted lines in Fig. 3(b). Fig. 4 shows the undeformed and deformed finite element mesh at three time instants (at the t-scale). The finite configuration changes due to both deformation and shape evolution are clearly visible. A more detailed view of the evolution of the shape of the undeformed contact surfaces is presented in Fig. 5 for the case of nodal shape parameterization (20 and 15 shape parameters for the pin and block, respectively). The instability of the time integration scheme is observed when the time step is too large, cf. Fig. 5(b), which is a well-known consequence of using the forward-Euler integration scheme [11,5–7]. Fig. 6 illustrates the application of a smooth shape parameterization independent of the finite element discretization. Now, the instability is not observed even when large time steps are used (the solution is obtained in 8 increments). However, the exact shape cannot be properly represented when the number of shape parameters is too small. This is clearly seen in the case of 5/6 (pin/block) shape parameters in Fig. 6(b). The effect of the time-step length on the accuracy of the shape update scheme (17) employing either first- or second-order approximation of wear increments is collectively shown in Fig. 7. The time step used for integration of the shape evolution problem is varied between Dt = 150 s (64 time increments) and Dt = 2400 s (4 time increments). The final shape obtained in each of the analyzed cases is referred to the final shape obtained for the nodal shape parameterization using the second-order scheme with 128 increments (Dt = 75 s), and the root-mean-square of the difference of nodal positions is reported in Fig. 7. Note that the number of time increments can be used as a measure of the computational cost (time) of the analysis. Care is only needed when comparing the first- and second-order shape update schemes. The latter are more expensive because the shape sensitivity analysis is carried out in addition to the direct analysis. The additional cost of sensitivity analysis depends on several factors, including the number of shape parameters, and it is proportional

Fig. 3. Pin-on-flat problem: (a) geometry and boundary conditions, (b) sketch of the shape parameterization (accounting for the symmetry, the pin and the block are parameterized here by, respectively, 5 and 6 independent shape parameters).

184

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

Fig. 4. Steady-state pin-on-flat problem: undeformed (top) and deformed (bottom) finite element mesh corresponding to (a) t = 0, (b) t = T/2 and (c) t = T (nodal shape parameterization, 128 increments, first-order scheme).

Fig. 5. Shape evolution due to wear obtained in (a) 128 and (b) 8 time increments (first-order scheme) using nodal shape parameterization. The initial, final and three intermediate shapes are shown; the final shape is plotted with a solid line.

to this number in the case of the direct differentiation method (DDM) employed in this work. However, in the present example, this additional cost is rather small (about 20% for 20/15 shape parameters). This is because the actual contact problem is pathindependent (frictionless contact), so the sensitivity problem is solved only once, after solving the direct problem. The cost of solving the minimization problem (17) itself is negligible. As expected, the second-order scheme is more accurate than the first-order scheme, although the difference between the two schemes decreases with decreasing number of shape parameters. Also, the error increases with decreasing number of shape parameters (except for large time increments). When smooth shape parameterization is applied, the error quickly stabilizes with increasing number of time increments (at 16–32 increments),

and further reduction of the time step does not give noticeable reduction of the error. 5.2. Reciprocating pin-on-flat problem The second example is a simplified two-dimensional arrangement corresponding to the reciprocating pin-on-flat test. The geometry is similar to that in the previous example except that the length of the block is now L = 25 mm. The pin is slid against the block at the constant normal force P, and the complete cycle consists of sliding the pin to the right and to the left with the amplitude of ±2R starting from the center of the block. The friction coefficient is equal to l = 0.1, the wear coefficients in the wear law (10) are now K1 = K2 = 5 105 mm2/N, and 9600 cycles are

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

185

root mean square error

Fig. 6. Shape evolution due to wear obtained in 8 time increments (first-order scheme) using smooth shape parameterization with (a) 10/12 and (b) 5/6 shape parameters. The initial, final and three intermediate shapes are shown; the final shape is plotted with a solid line.

10

1

10

2

10

3

1st order approximation of wear 2nd order approximation of wear nodal shape parametrization 20 15 parameters smooth shape parametrization 10 12 parameters 10

smooth shape parametrization 5 6 parameters

4

4

8

16

32

64

number of time increments Fig. 7. Steady-state pin-on-flat problem: root-mean-square error of the final shape as a function of the number of time increments used for the solution of the shape evolution problem.

Fig. 8. Reciprocating pin-on-flat problem: initial (top) and final (bottom) undeformed and deformed configurations at the beginning of the initial stage (a), at the beginning of the actual loading cycle (b), and in the right-most position after one quarter of the loading cycle (c).

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

root mean square error

186

10

1

10

2

10

3

10

1st order approximation of wear 2nd order approximation of wear nodal shape parametrization 20 25 parameters smooth shape parametrization 10 20 parameters smooth shape parametrization 5 10 parameters

4

4

8

16

32

64

number of time increments Fig. 9. Reciprocating pin-on-flat problem: root-mean-square error of the final shape as a function of the number of time increments used for the solution of the shape evolution problem.

simulated. All the remaining parameters of the problem are identical to those in the previous example. The deformation problem that is solved at each time instant at the t-scale includes an initial stage, in which the pin is brought into contact with the block in the left-most position, cf. Fig. 8(a), and then it is slid to the position at the center of the block, cf. Fig. 8(b). At this time instant, the actual loading cycle begins, which ends in the same position after the pin is slid the distance of 2R to the right, cf. Fig. 8(c), then the distance of 4R to the left, and again the distance of 2R to the right. Of course, accumulation of wear is suppressed during the initial stage. The above procedure assures that the path-dependent frictional response is treated properly, cf. Remark 1. In qualitative terms, the detailed results of the present cyclic loading problem are similar to those presented above for the case of the steady-state problem, and thus they are omitted here. In Fig. 9, we only summarize the results of the study of accuracy of the various shape update schemes. As in Fig. 7, the final shape obtained for the nodal parameterization using the second-order scheme with 128 increments (M = 75) is used as a reference solution in all cases. All the comments to Fig. 7 equally apply in the present case. As the numerical efficiency of the first- and second-order schemes is concerned, the additional computational cost of sensitivity analysis is here higher than in the previous example. In our implementation, the computational cost increases by a factor of 1.8, 2.7 and 3.4 for, respectively, 5/10, 10/20 and 20/25 shape parameters. However, considering the accuracy and the associated cost, the second-order scheme is still competitive, particularly in the case of the nodal parameterization (20/25 parameters). 6. Conclusion A continuum formulation has been developed for a general class of finite-wear contact problems. The concept of three configurations has been introduced to properly treat finite deformations due to contact interactions and finite shape changes due to wear. Furthermore, two time scales have been distinguished, so that the deformation problem and the wear problem can be partially decoupled. Considering the computational strategies, the space of shape evolution mappings has been approximated by a prescribed, finite-dimensional counterpart, parameterized by a time-dependent vector of shape parameters. The problem of shape evolution due to

wear has been transformed to an approximate problem of shape fitting with respect to shape parameters. The shape parameterization allows to uniquely define the shape evolution mapping for each point in the domain, not only the normal component on the boundary. It can thus be directly applied to map the finite element mesh after the shape update. Moreover, the associated smoothing of the wear profile appears to have a stabilizing effect on the shape-evolution integration schemes. Two shape update schemes have been considered in this work. The first-order scheme results from the direct application of the forward-Euler scheme to the rate shape-evolution minimization problem. As an extension, the second-order approximation of the wear increment has also been introduced, which requires shape sensitivities of the wear increments. Several combinations of the proposed shape update schemes and shape parameterizations have been considered, and their accuracy has been studied for two representative numerical examples. Improved accuracy of the second-order scheme has been confirmed, and this scheme combined with the nodal shape parameterization proved to be the most accurate. However, if shape sensitivities are not available and large time speedups need to be applied, then the use of shape smoothing techniques should be considered. In particular, the beneficial effect of smoothing of the wear profile, as introduced by smooth shape parameterizations, on the stability of the shape update schemes has been illustrated. Acknowledgements This work has been partially supported by the KomCerMet project (contract No. POIG. 01.03.01-14-013/08-00 with the Polish Ministry of Science and Higher Education) in the framework of the Operational Programme Innovative Economy 2007–2013. Appendix A. Frictional contact problem The deformation problem, as introduced in Section 2.2, is a standard frictional contact problem in which the deformation mappings uit ðX i ; sÞ; i ¼ 1; 2, are to be found, that satisfy the equilibrium equation, boundary and initial conditions, and contact constraints. The corresponding formulation is briefly outlined in this Appendix. For details the reader is referred to the monographs [28,29] and references cited therein. The subscript t denoting parameterization by the slow time t is omitted below.

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

The formulation of the large deformation contact problem is based on the master–slave approach and on the notion of the closest-point projection. The contact pair is defined by the projection of a point x1 of the deformed slave surface c1c ¼ u1 ðC1c ; sÞ on the deformed master surface c2c ¼ u2 ðC2c ; sÞ. The projection point is de2 . Introducing parameterization of the master surface c2c noted by x 2 ¼ x2 ð by convective coordinates n = {n1, n2}, we have x nÞ, where  n ¼ f n1 ;  n2 g are the coordinates of the projection point. The normal gap gN and the sliding velocity vT are the basic kinematic contact variables defined as follows:

2

t N ¼ t  n; 2

t T ¼ tT a sa ;

ðA:2Þ a

2

where t = r n, r is the Cauchy stress, and s is the cobasis, sa  sb ¼ dab . In view of the action–reaction principle we have t2 = t1. The kinematic contact variables and the contact tractions are related by the contact constraints, namely the unilateral contact condition,

g N P 0;

t N 6 0;

g N tN ¼ 0;

ðA:3Þ

and the Coulomb friction law,

kv T kt T ¼ v T kt T k;

kt T k þ lt N 6 0;

kv T kðkt T k þ ltN Þ ¼ 0:

ðA:4Þ

The virtual work principle, which is the basis for the finite element treatment, is written in the following form:

Gðu; duÞ ¼ G1 ðu1 ; du1 Þ þ G2 ðu2 ; du2 Þ þ Gc ðu; duÞ ¼ 0; 1

ðA:5Þ

2

where u = {u , u }. Here, Gi denotes the virtual work of internal and external forces, excluding the contact interactions, defined for each body separately,

Gi ðui ; dui Þ ¼

Z

i

Xi

P i  rdui dV 

Z Cit

i

T i  dui dS ;

ðA:6Þ

where Pi is the first Piola–Kirchhoff stress, and T⁄i is the surface traction prescribed on the boundary Cit . The virtual work Gc developed by the contact forces takes the form

Gc ðu; duÞ ¼

Z

1

C1c

ðT N dg N þ T T a dna ÞdS :

ðA:7Þ

The contact contribution is integrated over the undeformed contact surface C1c of the slave body, hence the nominal contact tractions are introduced, 1

T N ¼ j tN ;

1

T T ¼ j tT ;

ðA:8Þ

1

where j is the area transformation factor of the slave surface (ds1 = j1 dS1), so that the nominal contact tractions TN and TT refer to the unit area in the undeformed configuration X1 of the slave body. The virtual work principle (A.5) must be combined with a suitable regularization technique in order to enforce the contact conditions (A.3) and (A.4). In this work, the augmented Lagrangian method is used in the version proposed in [35,36]. Using the Archard wear law (10), the nominal wear rate of the slave surface is directly determined in terms of the nominal friction traction TT, namely

_ 1 ¼ K 1 D_ 1 ; W

D_ 1 ¼ T T a n_ a :

_ 2 ¼ K 2 D_ 1 ; W 

ðA:10Þ

_ 2 ds can be interpreted as the wear volume of the master where W  surface referred to the unit area of the undeformed slave surface. _ 1 and W _ 2 correspond to We remind here that the wear rates W 2 Þ, which is defined by the closest-point prothe contact pair ðx1 ; x jection at the current time instant.

ðA:1Þ

where the unit outer normal n = n2 is adopted as the normal of the contact pair and sa = @x2/@na, a = 1, 2, is the tangent basis. The normal contact traction tN and the tangential traction tT are defined in terms of the spatial (Cauchy) traction vector t = t2,

t ¼ tN n þ tT ;

2

_ 2 ¼ j K 2 D_ 1 ¼ j W _ 2; W  1 1 j j

References

v T ¼ n_ a sa ;

2 Þ  n; g N ¼ ðx1  x

2

187

ðA:9Þ

At the same time, the nominal wear rate of the master surface is given by

[1] S.C. Lim, Recent developments in wear-mechanism maps, Tribol. Int. 31 (1998) 87–97. [2] H.C. Meng, K.C. Ludema, Wear models and predictive equations: their form and content, Wear 181–183 (1995) 443–457. [3] J.F. Archard, Contact and rubbing of flat surfaces, J. Appl. Phys. 24 (1953) 981– 988. [4] C. Agelet de Saracibar, M. Chiumenti, On the numerical modeling of frictional wear phenomena, Comput. Methods Appl. Mech. Engrg. 177 (1999) 401–426. [5] P. Podra, S. Andersson, Simulating sliding wear with finite element method, Tribol. Int. 32 (1999) 71–81. [6] M. Oqvist, Numerical simulations of mild wear using updated geometry with different step size approaches, Wear 249 (2001) 6–11. [7] I.R. McColl, J. Ding, S.B. Leen, Finite element simulation and experimental validation of fretting wear, Wear 256 (2004) 1114–1127. [8] V. Hegadekatte, N. Huber, O. Kraft, Modeling and simulation of wear in pin on disc tribometer, Tribol. Lett. 24 (2006) 51–60. [9] C. Paulin, S. Fouvry, C. Meunier, Finite element modelling of fretting wear surface evolution: application to a Ti–6Al–4V contact, Wear 264 (2008) 26– 36. [10] L. Gallego, B. Fulleringer, S. Deyber, D. Nelias, Multiscale computation of fretting wear at the blade/disk interface, Tribol. Int. 43 (2010) 708– 718. [11] L. Johansson, Numerical simulation of contact pressure evolution in fretting, Trans. ASME J. Tribol. 116 (1994) 247–254. [12] N. Stromberg, An augmented Lagrangian method for fretting problems, Eur. J. Mech. A/Solids 16 (1997) 573–593. [13] F. Jourdan, A. Samida, An implicit numerical method for wear modeling applied to a hip joint prosthesis problem, Comput. Methods Appl. Mech. Engrg. 198 (2009) 2209–2217. [14] M. Peigney, Simulating wear under cyclic loading by a minimization approach, Int. J. Solids Struct. 41 (2004) 6783–6799. [15] I. Paczelt, Z. Mróz, On optimal contact shapes generated by wear, Int. J. Numer. Methods Engrg. 63 (2005) 1250–1287. [16] F. Jourdan, Numerical wear modelling in dynamics and large strains: application to knee joint prostheses, Wear 261 (2006) 283–292. [17] N. Stromberg, L. Johansson, A. Klarbring, Derivation and analysis of a generalized standard model for contact, friction and wear, Int. J. Solids Struct. 33 (1996) 1817–1836. [18] Z. Mróz, S. Stupkiewicz, An anisotropic friction and wear model, Int. J. Solids Struct. 31 (1994) 1113–1131. [19] S. Fouvry, P. Kapsa, L. Vincent, Quantification of fretting damage, Wear 200 (1996) 186–205. [20] A. Ramalho, J.C. Miranda, The relationship between wear and dissipated energy in sliding systems, Wear 260 (2006) 361–367. [21] M. Dragon-Louiset, On a predictive macroscopic contact-sliding wear model based on micromechanical considerations, Int. J. Solids Struct. 38 (2001) 1625– 1639. [22] C. Stolz, A thermodynamical approach to contact wear as application of moving discontinuities, Arch. Appl. Mech. 77 (2007) 165–175. [23] P. Michaleris, D.A. Tortorelli, C.A. Vidal, Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity, Int. J. Numer. Methods Engrg. 37 (1994) 2471–2499. [24] M. Kleiber, H. Antunez, T.D. Hien, P. Kowalczyk, Parameter Sensitivity in Nonlinear Mechanics, Wiley, Chichester, 1997. [25] S. Stupkiewicz, J. Korelc, M. Dutko, T. Rodicˇ, Shape sensitivity analysis of large deformation frictional contact problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3555–3581. [26] S. Stupkiewicz, J. Lengiewicz, J. Korelc, Sensitivity analysis for frictional contact problems in the augmented Lagrangian formulation, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2165–2176. [27] J. Lengiewicz, J. Korelc, S. Stupkiewicz, Automation of finite element formulations for large deformation contact problems, Int. J. Numer. Methods Engrg., in press. doi:10.1002/nme.3009. [28] T.A. Laursen, Computational Contact and Impact Mechanics, Springer-Verlag, Berlin, 2002. [29] P. Wriggers, Computational Contact Mechanics, Wiley, Chichester, 2002. [30] G. Pietrzak, Continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems, Ph.D. Thesis, Lausanne, EPFL, 1997.

188

J. Lengiewicz, S. Stupkiewicz / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 178–188

[31] G. Zavarise, P. Wriggers, E. Stein, B.A. Schrefler, Real contact mechanisms and finite element formulation—a coupled thermomechanical approach, Int. J. Numer. Methods Engrg. 35 (1992) 767–785. [32] S. Stupkiewicz, Extension of the node-to-segment contact element for surfaceexpansion-dependent contact laws, Int. J. Numer. Methods Engrg. 50 (2001) 739–759. [33] E.A. de Souza Neto, D. Peric´, M. Dutko, D.R.J. Owen, Design of simple low order finite elements for large strain analysis of nearly incompressible solids, Int. J. Solids Struct. 33 (1996) 3277–3296.

[34] J. Korelc, AceGen and AceFEM user manuals, 2009. Available at: . [35] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92 (1991) 353–375. [36] G. Pietrzak, A. Curnier, Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment, Comput. Methods Appl. Mech. Engrg. 177 (1999) 351–381.