- Email: [email protected]
Determination of the shakedown limit for large, discrete frictional systems
European Journal of Mechanics A/Solids 49 (2015) 242e250
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
Determination of the shakedown limit for large, discrete frictional systems R.C. Flicek a, *, D.A. Hills a, J.R. Barber b, D. Dini c a
Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, United Kingdom Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109-2125, USA c Department of Mechanical Engineering, Imperial College London, South Kensington Campus, Exhibition Road, London SW7 2AZ, United Kingdom b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 25 March 2014 Accepted 3 August 2014 Available online 13 August 2014
It has been shown that complete frictional contacts subjected to cyclic loads can shake down such that the interface becomes fully stuck after some number of cycles (even if partial-slip occurs at the onset of loading). However, if the amplitude of the cyclic load is greater than a particular value (the shakedown limit), it is impossible for shakedown to occur. In this paper, we examine a numerical approach for determining the shakedown limit for elastic frictional systems subjected to quasi-static loads using a discrete formulation. We then use this technique to determine the shakedown limit for a finite element model of a (coupled) complete contact with ~50,000 total degrees of freedom and ~250 along the contact interface. Finally, we compare the calculated value of the shakedown limit to a series of over 1000 transient simulations and investigate the influence of initial conditions on steady state frictional energy dissipation. The results demonstrate that the dissipative properties of complete contacts can be highly dependent on the initial residual slip displacement state. © 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Shakedown limit Cyclic loading Fretting fatigue
1. Introduction It is extremely common for frictional contacts in engineering systems to be subjected to a combination of a static load (to form the contact) and a time varying cyclic load (e.g. due to vibrations). Some examples include riveted and bolted joints (Abad et al., 2012; Law et al., 2006), dove-tail connections in jet engines (Rajasekaran and Nowell, 2006; Xi et al., 2000), and spline couplings (Limmer et al., 2001; Banerjee and Hills, 2006; Banerjee et al., 2008) to name just a few. Furthermore, it is well known that even when the contact interface appears nominally to be stuck, these loading conditions frequently result in small zones of micro-slip (or partialslip) along the contact interface (where the relative tangential displacement is generally of the order 25100 mm (Szolwinski and Farris, 1996)). The effect of micro-slip on system performance is two-fold: on the one hand, it results in the well known damage process of fretting fatigue, which significantly reduces the service life of components (Farris et al., 2000); but on the other hand, it provides a source of damping, which can be beneficial for system performance (Law et al., 2006).
* Corresponding author. E-mail address: robert.fl[email protected] (R.C. Flicek). http://dx.doi.org/10.1016/j.euromechsol.2014.08.001 0997-7538/© 2014 Elsevier Masson SAS. All rights reserved.
Much of the previous work on micro-slip and fretting fatigue has focused on incomplete (or non-conforming) contacts (e.g. (Nowell and Hills, 1987; Hills et al., 2012)), which are characterized by the edge of contact being smooth such that the bodies do not ‘conform’ in the absence of applied loads. This is probably due to their relevance in engineering practice but also because closed form solutions exist for many relevant loading cases, e.g. (Mindlin €ger, 1998; Ciavarella, 1998). For these and Deresiewicz, 1953; Ja contacts, the contact area is a function of the applied normal load, and the contact pressure distribution smoothly tends to zero as the contact edge is approached. As a direct consequence of this, a zone of micro-slip will almost always arise at the edge of an incomplete contact that is subjected to cyclic loads (Barber et al., JanuaryeFebruary 2008). Complete (or conforming) contacts, such as the contact shown in Fig. 1, are characterized by the two contacting bodies ‘conforming’ in the absence of the applied loads. For this class of contacts, the edge of contact is sharp, and the contact area either does not change in the presence of applied loads or becomes smaller (due to the edges ‘lifting off’).1 Because sharp features result in a large
1 For this reason, this type of contact is sometimes referred to as a receding contact (Dundurs and Stippes, 1970).
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 1. The complete contact formed between a square elastic punch and an elastically similar half-plane showing: the static normal load, P, the time-varying load, s(t), the size of the contact, 2a, and the (x,y) coordinate set.
stress concentration in their vicinity, complete contacts are usually avoided when possible. However, these contacts do arise in some critical applications such as in the spline couplings that join the central shafts in jet engines (Limmer et al., 2001; Banerjee and Hills, 2006). Numerical analysis of complete contacts subjected to cyclic loads sometimes predicts that the contact shakes down (Churchman and Hills, June 2006; Banerjee and Hills, 2006): that is, after the first few loading cycles, a favourable residual displacement state is developed that inhibits slip for all subsequent cycles. This behaviour clearly has some similarities to cyclic plasticity problems, which has prompted the question as to whether a theorem analogous to Melan's plastic shakedown theorem (Melan, 1936) could be developed for frictional contacts. In other words, if there exists a residual displacement distribution that would inhibit slip at all times during load cycle (i.e. a safe displacement distribution), is the existence of this displacement distribution sufficient to guarantee that the contact will shake down? This question has been answered unequivocally by Klarbring et al. (December 2007) for discrete contact problems and by Barber et al. (JanuaryeFebruary 2008) for continuum problems. These authors have proven that the existence of a safe displacement distribution is both a necessary and sufficient condition to guarantee shakedown if the contact is uncoupled: that is, if relative tangential contact displacements have no influence on the distribution of normal contact traction. Thus, for uncoupled contacts, there exists a load factor, i.e. a ratio of cyclic load amplitude to static load (e.g. s/P in Fig. 1), above which shakedown is impossible and below which it is guaranteed; hence the contact's steady state response never depends on initial conditions. On the other hand, if the contact is coupled, the existence of a safe displacement distribution is still a necessary condition for shakedown to be possible, but it is no longer sufficient. In fact, for coupled contacts, there almost always exists a range of load factors for which the steady state response depends on the initial residual displacement distribution (Klarbring et al., December 2007; Barber et al., JanuaryeFebruary 2008). Ahn et al. (December 2008) have demonstrated this behaviour for a two-node system and have shown that for coupled contacts there exists: 1. a load factor below which shakedown is guaranteed to occur, denoted l1 2. a load factor above which shakedown is impossible, denoted l2. Note that these two limits on the applied load may be zero for many contacts, i.e. l1 ¼ l2 ¼ 0, (e.g. for incomplete contacts), implying that shakedown is impossible. Moreover, even if l2 > 0, it
243
is still possible that l1 ¼ 0, implying that although shakedown is sometimes possible, whether it actually occurs is always dependent on initial conditions. Ahn et al. (December 2008) also presented a method for calculating both l1 and l2. Unfortunately, their approach involves solving a system of equations that becomes combinatorially more complex as the number of degrees of freedom in the system increases. For instance, Jang and Barber (2011) were able to apply this approach to determine both l1,l2 for a system of 10 cracks, but applying this analysis to a significantly larger frictional system would be computationally prohibitive. As it is common for the finite element models used in practice to incorporate one or two orders of magnitude more nodes along the frictional interface than are considered by Jang and Barber (2011), it is of practical interest to develop a more computationally efficient approach. € rkman and Klarbring (1987) have presented an efficient Bjo approach for calculating l2, which we will henceforth refer to as the shakedown limit. To do this, these authors formulate the problem in discrete form, and show that the shakedown limit calculation can be posed as a Linear Programming problem for which several efficient solution algorithms exist. These authors then apply this calculation to two example problems, and also compare the results to numerical results obtained using the finite element method. The aim of this paper is: (i) to re-visit the Linear Programming approach for calculating the shakedown limit and (ii) to examine the influence of initial conditions on the level of frictional dissipation present in the steady state. This is of practical importance because energy dissipation is associated with damage processes such as fretting fatigue and wear, which are themselves linked to component life. In particular, we wish to investigate whether energy dissipation is significantly influenced by initial conditions or if their effect can safely be neglected. The structure of the paper is as follows. In x 2, we formulate the discrete contact problem. In x 3, we outline how to determine when the onset of slip or separation will occur for a contact subjected to monotonic loading. In x 4, we describe how the shakedown limit calculation can be posed as a Linear Programming optimization problem for which several solution algorithms already exist. In x 5, these calculations are applied to the example (coupled) complete contact shown in Fig. 1. Over 1000 numerical simulations of the time-evolution of the contact are then performed under a wide range of cyclic loads and initial displacement states, and the influence of initial conditions on steady state frictional energy dissipation is investigated. Finally, the implications of these results are discussed in x 6.
2. Formulation Let us consider some two-dimensional complete contact geometry, such as that shown in Fig. 1, which has been discretized using the finite element method such that there are N nodes in potential frictional contact along the interface. We can write the reaction forces, r, and the relative displacements, u, at the contact nodes as
r ¼ ½q; pT ¼ ½q1 ; …; qN ; p1 ; …; pN T
(1a)
u ¼ ½v; wT ¼ ½v1 ; …; vN ; w1 ; …; wN T ;
(1b)
where qi,pi are the shear, normal reactions and vi,wi are the tangential, normal relative displacements, and where i2{1,…,N}. Here we adopt the convention that pi is positive in compression (the negative y-direction), wi positive for a positive gap (the
244
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
positive y-direction), and qi,vi are positive in the positive x-direction with respect to the x,y coordinate system shown in Fig. 1. We can write an expression for the time-evolution of the contact reactions as
rðtÞ ¼ KuðtÞ þ r w ðtÞ;
(2)
where K is a symmetric 2N 2N contact stiffness matrix and rw accounts for the influence of the external loads on the contact reactions. More specifically, rw represents the reactions, r, that would be generated by the external loads if all the contact displacements were fixed to be zero, i.e. if u ¼ 0. Note that both K and rw can be obtained from a standard finite element model using the static reduction procedure described by Thaitirarot et al. (2014). Also note that the contact stiffness matrix, K, is positive semi-definite (and singular) if either contacting body has rigid body degrees of freedom; otherwise K is positive definite (and non-singular).
qðtÞ pðtÞ
¼
A B
BT C
vðtÞ wðtÞ
þ Fs
qws pws
þ Fc ðtÞ
qwc pwc
;
(5)
where q,p,v,w are defined in equation (1), and K is partitioned into sub-matrices corresponding to the following forceedisplacement relationships: tangentialetangential, A; normaletangential, B; and normalenormal, C. Note that as K is symmetric and positive semidefinite (or positive definite), so too are A and C, whereas B is generally not symmetric. We also define the load factor, l, as
l¼
Fc ðtÞjmax Fc ðtÞjmin Fcmax Fcmin ¼ ; Fs Fs
(6)
where Fcmin , Fcmax are the minimum, maximum values of the cyclic load, respectively. Note that the load factor is only defined for cyclic loading regimes. 3. First violation of the stick condition
2.1. The Coulomb friction law Before defining the friction law, we note that the contacting bodies cannot inter-penetrate and that the normal contact reactions cannot be tensile:
wi 0;
pi 0:
(3a)
Furthermore, if there is a positive gap at node i, its normal and shear reactions must be zero. Conversely, if the normal reaction at node i is compressive, its gap must be zero:
wi > 00pi ¼ qi ¼ 0
(3b)
pi > 00wi ¼ 0:
(3c)
Finally, we assume that the Coulomb friction condition applies at each node, and this requires that
jqi j fpi ;
The load ratio, Fc/Fs, at which slip or separation first occurs under monotonic loading is always dependent on the initial residual slip displacement distribution, v(0). Nevertheless, this information may be helpful when planning or analysing experiments, particularly when the case of a null initial slip distribution, i.e. v(0) ¼ 0, is considered. Moreover, this calculation is trivially simple to carry out. Thus, suppose we have a complete contact with v(0) ¼ 0 subject to a loading regime in which the static load, Fs, is first applied and held constant, and the time-varying load, Fc, is subsequently applied monotonically. We can determine the load ratio at which slip or separation first occurs by checking when the stick condition is first violated. In other words, we search for the lowest value of Fc/Fs that violates the condition pi 0 or jqij < fpi at any node, i. Recall that for complete contacts, w ¼ 0 in the absence of external loads by definition. Thus we can use equation (5) to write these conditions as
pws i þ
Fc wc p 0 Fs i
(7a)
Fc wc Fc wc qi < f pws p : þ i Fs Fs i
(7b)
f >0
(3d)
jqi j < fpi 0v_i ¼ 0
(3e)
qws i þ
0 < jqi j ¼ fpi 0sgnðv_i Þ ¼ sgnðqi Þ;
(3f)
Note that these conditions are only valid for predicting first violations, so once slip, for example, is predicted to initiate from equation (7b), separation is no longer accurately predicted from equation (7a).
where f is the coefficient of friction and a superposed dot denotes the derivative with respect to time, t. Thus equation (3d) defines the admissible range of shear reactions; equation (3e) states that if the magnitude of qi is less than a critical value, node i is stationary; and equation (3f) states that if jqij is equal to a critical value, slip occurs in the direction that opposes the shear reaction at that node.
2.2. Loading regime We consider a loading regime comprised of a static component, Fs, superposed with a quasi-statically applied time-varying (cyclic) component, Fc(t), where Fs,Fc are scalars that account for the strength of the applied load. The effect of these applied loads on contact reactions can be expressed as
rw ðtÞ ¼ Fs r ws þ Fc ðtÞr wc ;
(4)
where rws,rwc are the rw vectors corresponding to a unit application of Fs,Fc, respectively. With this notation, equation (2) can be rewritten as
4. The shakedown limit It is perhaps most natural to think about frictional shakedown as it might arise in the context of an experiment: namely, there is some test configuration that is subjected to a particular combination of loads, and the coefficient of friction (which we will assume remains constant) is determined by the material properties, surface treatment, etc. Thus if we wanted to determine the shakedown limit experimentally, we might apply cyclic loads of various amplitudes, vary the initial displacement state, and use digital image correlation to determine the maximum cyclic load amplitude for which the contact shakes down. In contrast, the approach we describe here for calculating the shakedown limit is, in some sense, the reverse of this hypothetical experimental approach. To pose the calculation of the shakedown limit as a Linear €rkman and Klarbring, 1987), we Programming problem as in (Bjo begin by choosing a particular cyclic loading regime in which the cyclic load ranges from some minimum value, Fcmin , to some
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
maximum value, Fcmax . We then assume that the load factor (see equation (6)) corresponding to this load range is the shakedown limit if the coefficient of friction is some currently unknown value, denoted fmin. Our aim therefore becomes to determine fmin: that is, the coefficient of friction below which it is impossible for the contact to shake down for the specified load range. We can then perform this calculation for a series of load ranges to determine the shakedown limit for a range of friction coefficients. Note that if the load range is sufficiently large, it may be impossible for the contact to shake down irrespective of the coefficient of friction, making our initial assumption false. This can occur if there does not exist any residual slip displacement that inhibits separation throughout the load cycle. In this case, the calculation will simply fail to converge, implying shakedown is impossible for the specified load range.
To determine fmin, first recall that if a contact is to shake down, separation cannot occur at any time in the load cycle, i.e. w(t) ¼ 0. Therefore, for a contact that has shaken down, we can use equation (5) to write the ratio of the (element-wise) absolute value of the shear reactions to the normal reactions, R, as
Rðv; Fc Þ ¼
Note that whereas we have formulated the problem as a €rkman and ‘minimax’ constrained optimization problem, Bjo Klarbring (1987) formulate it as a constrained maximization. However, both approaches are equivalent, and the reader is € rkman and Klarbring, 1987) for a directed to the original paper (Bjo discussion of the mathematical structure of this optimization. 4.2. Improving convergence One additional constraint that is often helpful for improving convergence of the optimization in equation (11) is to eliminate the potential for rigid body motion in the slip displacement distribution. This can be done by removing one node's slip displacement from the optimization and optimizing against a subset of the slip displacements, v*, where
v ¼ ½v1 ; …; vN1 T :
4.1. Calculation
jqðFc Þj jAv þ Fs qws þ Fc qwc j ¼ : pðFc Þ Bv þ Fs pws þ Fc pwc
(8)
Notice that we have written the contact reactions as a function of Fc as opposed to a strict function of time, t. We can do this because every load cycle for a shaken down contact is identical, and the reactions simply vary with the applied loads. Also notice that v is not a function of Fc because, by definition, slip cannot occur at any point during the load cycle for a contact that has shaken down. According to the friction law defined in x 2.1, slip will occur if any element of R becomes equal to the coefficient of friction. As we are considering the case of quasi-static loading, R will take on its maximum value at one (or both) of the extreme points (in time) of the load cycle, i.e. at Fcmin or Fcmax. We can construct a vector, S(v), of length 2N that is comprised of the contact reaction ratios, R, at the extreme points in the load cycle as
245
(12)
The slip displacement vN, which is not optimized against, can then be set to
vN ¼
N1 X
vi
(13)
i¼1
to force the sum of all slip displacements to be zero. Note that it is not important which node's slip displacement is removed from the optimization.2 5. Example application
(11)
We now intend to illustrate how these calculations can be used in practice by applying them to the example problem shown in Fig. 1. This figure shows the complete contact formed between two elastic continuum bodies: a square punch and an elastically similar half-plane. As our approach uses the discrete formulation, the first step in our analysis is to discretize the contact, which we have done using the commercial finite element software ABAQUS/CAE. Of course once we discretize these bodies, we are no longer studying the same problem in a rigorous sense. However, in practice it is usually assumed that a continuum is well represented by a discrete system with a sufficiently refined mesh. The details of the finite element model we have used are in Appendix A, but here it suffices to state that both bodies have elastic modulus E ¼ 200 GPa, Poission's ratio n ¼ 0.3, the model assumes plane strain with unit depth, the contact interface comprises 128 uniformly spaced nodes, and the full model incorporating the punch and the half-plane comprises 50,040 degrees of freedom, where every node in the model is associated with two degrees of freedom. Note that the finite element model we used does not take advantage of the symmetry about the vertical centreline of the contact. This is simply because we generated these results from a set of models we created to study several different loading regimes, and some of these loads do not maintain this horizontal symmetry, e.g. when a uniform shear traction is applied to the top of the punch. To determine the sensitivity of the contact solution to the discretization, we performed a mesh convergence study in which we varied both the number of contact nodes and the total number of nodes in the model. The results of this study suggest that a model with 128 contact nodes provides a converged result at the interior
where j2{1,…,2N}. Several algorithms exist for the solution of such problems, e.g. (Brayton et al., 1979; Powell, 1978). Note that the solution to this problem, vopt, generally is not unique.
2 We tested this by removing nodes at various positions along the contact from the optimization and specifying their slip displacement according to equation (13); the solution was insensitive to which node was chosen.
) ( R v; Fcmin SðvÞ ¼
: R v; Fcmax
(9)
Thus further slip will be inhibited throughout the load cycle if
f SðvÞjmax ;
(10)
i.e. if the coefficient of friction is greater than the maximum element value of S. To determine fmin, we simply need to optimize v to determine an optimal slip distribution, vopt. For our purposes, an optimal slip distribution is one that: (i) prevents separation at all times in the load cycle and (ii) requires the lowest coefficient of friction to prevent slip from occurring. From equation (10), we see that (ii) requires that vopt minimizes the maximum value of S while (i) requires that the condition p 0 is simultaneously satisfied. Therefore, finding vopt is a standard ‘mini-max’ optimization problem of the form
min max Sj ðvÞ v
j
such that p 0;
246
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
of the contact but that the solution close to the corner is not fully converged. However, we found that obtaining a converged result at the corners with a uniform mesh along the contact would be computationally impractical, so a more complex biased mesh would be required. As our primary object is to demonstrate the technique, we feel that this simple mesh is sufficient. However, note that the technique presented in x 4.1 can be applied to larger problems without modification. The next step in this analysis is to apply the static reduction technique described by Thaitirarot et al. (2014), reducing the finite element model to a discrete system of the N contact nodes alone. In doing this, we obtain the contact stiffness matrix, K, and the rw vectors corresponding to each applied load, P,s(t). Note that P corresponds to the mean contact pressure, and s(t) corresponds to the mean bulk load; hence both loads have units of [FL2], where F,L denote force, length, respectively. Thus we can write the applied loads for this example problem in the form given by equation (4) as
rw ðtÞ ¼ PrP þ sðtÞr s ;
(14)
where rP,rs have units of [F/(FL2)] or simply [L2] and correspond to the application of a unit load in the directions shown in Fig. 1. Finally, the load factor is defined as
l¼
sðtÞjmax sðtÞjmin smax smin ¼ ; P P
that is implied for monotonically increasing/decreasing s, which corresponds to a bulk tension/compression, respectively. Note that we adopt the convention that leading edge slip refers to the edges of the punch slipping outwards, whereas trailing edge slip refers to the edges of the punch slipping inwards. The first thing to notice from this figure is that if f < 0.48, the first violation of the stick condition occurs before s is applied; that is, leading edge slip initiates from the contact edge on the application of P alone. It is worth mentioning that an asymptotic analysis of the contact edge in the continuum formulation predicts that edge slip will occur if f < 0.543 (see (Hills and Dini, January 2004)). The discrepancy between these results is due to the mesh being insufficiently fine near the contact edges. In any case, when f > 0.48, the entirety of the contact interface remains stuck on the application of P. If s is then monotonically applied in tension, slip initiates in the trailing edge sense when s/ P~2, depending on the coefficient of friction. Conversely, if s is monotonically applied in compression, slip initiates in the leading edge sense within the range 2.3 < s/P < 0. Note that although slip initiates somewhat inwards from the contact edge when f > 0.48 (in both tension and compression), the point of first slip always remains quite close to the contact edge. Also notice that for f < 1, the first violation of the stick condition always results from slip, not separation.
(15)
where smin and smax are the minimum and maximum values of s, respectively. 5.1. First violation of the stick condition We can now determine the load ratio, s/P, at which slip or separation first occurs for this contact when v(0) ¼ 0. To do this, we assume that P is first raised to some value and then held constant, and that s(t) is subsequently applied monotonically (in tension or compression). We then check for violations of the conditions specified in equation (7), and the results of this calculation are shown in Fig. 2 on a plot of s/P vs. f. The way to interpret this figure is to pick a coefficient of friction, which corresponds to some vertical line, and to view the behaviour
5.2. The shakedown limit We now intend to calculate the shakedown limit, l2, and to compare it with the results of a series of transient simulations. The loading history we consider is shown in Fig. 3: namely, P is first applied and held constant, and s(t) is then varied within the range 0 s(t) smax. In other words, we consider the case of ‘zero-topeak’ loading in cyclic tension; thus the load factor is defined as l ¼ smax/P. Also note that the initial slip displacement condition, i.e. v(0), need not be null; indeed, the transient simulations we carried out explored a wide range of v(0) distributions.3 We calculated l2 for the loading regime described above in MATLAB R2013a using the fminimax function, and the results are shown in Fig. 4 on a plot of l vs. f. Also shown in this figure are the results of 347 ‘marching in time’ transient simulations,4 which we carried out using the Gauss-Seidel algorithm described by Ahn and Barber (2008). Simulations that eventually shook down are marked with a circle, and those that failed to shake down are marked with a dot. The main thing to recognize about this figure is that it demonstrates that the calculated values of l2 are consistent with results of all the transient simulations we performed; that is, every transient simulation for l > l2 failed to shake down. 5.3. Frictional energy dissipation All the transient simulations we performed were subjected to the loading regime shown in Fig. 3 until a steady state response was reached, which we determined using a criterion based on the frictional energy dissipation per cycle, W. Specifically, the contact was determined to have reached a steady state response if: (i) W was similar in size to the smallest meaningful value of W that can be represented with a double precision floating point number or
Fig. 2. The slip-stick behaviour of the contact shown in Fig. 1 when P is first applied to a contact with a null initial slip displacement distribution, and s is then applied monotonically in tension (positive) or compression (negative).
3 Note that we varied v(0) by assuming a particular displacement distribution and then simply scaled this distribution. This approach enables any distribution to be used, but we mostly assumed either a linear or sinusoidal distribution or the optimal distribution obtained from the calculation of l2. 4 It is difficult to identify some of these simulations in this and subsequent figures since we ran multiple simulations with the same l but with different v(0).
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 3. The cyclic loading regime considered for the shakedown limit calculation.
(ii) the change in dissipation from the previous load cycle was less than 0.01%. We calculated the frictional energy dissipation per cycle as
W¼
I N X i¼1
display a value of zero; thus, we plot simulations that shook down (for which WSS ¼ 0) as the arbitrarily small value of 108. Both Figs. 4 and 5 demonstrate that the range of l for which shakedown is conditional (i.e. where some simulations shook down while others did not) can be quite large. For example, Fig. 5 shows that simulations failed to shake down when l was as low as 66%, 50%, and 62% the value of l2 when f was 0.3, 0.6, and 0.9, respectively. In addition, Figs. 4 and 5 suggest that the size of the conditional region (as a percentage of l2) remains quite similar over the range of f we have considered. However, it is almost certain that the conditional region is at least somewhat larger than is suggested by these figures because it is unlikely that any of the v(0) distributions we tested were the absolute ‘worst’ v(0) possible (in terms of encouraging shakedown). 5.4. Optimal initial residual slip displacement
qi ðtÞv_i ðtÞdt;
(16)
cycle
where the time-integration was performed numerically with a piecewise-linear approximation in each load increment. Note that W is a non-negative quantity because equation (3) states that the sign of qi always opposes that of v_i for v_i s0. We performed a convergence study and found a load increment of 0.05 (with respect to s) was sufficient to obtain a converged result; thus this load increment was used for all the transient simulations we performed. As this model is two-dimensional, W is the dissipation per load cycle per unit depth; hence its units are [(FL)/L] or simply [F]. Thus we can write the dissipation per cycle in the steady state in dimensionless form, WSS, which is given by
WSS ¼
247
E
W; a2 P 2 1 n2
(17)
where 2a is the contact-width. To illustrate the effect that initial conditions can have on WSS, we plot WSS against l in Fig. 5 for 873 transient simulations carried out at three coefficients of friction: 0.3, 0.6, and 0.9. Each of these plots also includes a vertical line corresponding to l2 at the relevant coefficient of friction. Note that as the vertical scale on these plots is logarithmic, they cannot
Fig. 4. The (calculated) shakedown limit and the results of 347 transient simulations that were run until a steady state was reached. Transient simulations that shook down are shown with a circle, and those that did not shake down are shown with a dot.
Fig. 5 illustrates that the frictional energy dissipated by a complete contact after several load cycles can be very sensitive to the initial residual slip displacement state, v(0). Hence, it is of practical interest to determine which sorts of initial displacement conditions result in a more/less dissipative contact. Conveniently, we automatically obtain an optimal slip distribution for promoting shakedown when calculating l2. Fig. 6 shows three of these distributions in dimensionless form plotted against normalized position along the contact interface, x/a, for f ¼ 0.3,0.6,0.9, where the dimensionless slip displacement is given by
vdm ¼
E
v: aP 1 n2
(18)
Although the information provided by Fig. 6 is useful, it is clear that pre-loading a contact with distributions such as these might be difficult in practice. However, a linear distribution of v(0) can be achieved by the simple expedient of applying an appropriate bulk stress, s, before the constant normal load, P, is applied. Hence we might ask: what is the best/worst initial condition for promoting shakedown that can be achieved by this procedure? To address this question, we performed a series of simulations with different bulk pre-loads, and the results of these simulations are shown in Fig. 4. In addition, we show the effect of a (dimensionless) bulk pre-load, sdm, on steady state dissipation, WSS, in Fig. 7 for the example case of f ¼ 0.6 and l ¼ 3, where the dimensionless bulk load is given by sdm ¼ s/P. Fig. 7 illustrates that shakedown can be encouraged by applying a tensile bulk pre-load prior to the formation of the contact, or it can be inhibited with a compressive bulk pre-load. In fact, even though l ¼ 3 is 97% the value of l2 for f ¼ 0.6, Fig. 7 demonstrates that shakedown can still be achieved by applying a bulk pre-load of sdm 1.5. Although it is perhaps surprising how large a range of WSS can be excited using a bulk pre-load alone, it should not be surprising that tensile/compressive bulk pre-loads promote/inhibit shakedown. This is because all the optimal slip distributions shown in Fig. 6 consist mainly of trailing edge slip, and trailing/leading edge slip is generated by a tensile/compressive bulk load. Also notice that the effect of the bulk pre-load on WSS ‘saturates’ quite quickly and that the largest effect occurs near sdm~1. This saturation effect is not surprising because if excessively large bulk pre-loads are applied, there is simply a larger amount of slip in the first cycle, but the steady state remains largely unchanged. Also, note that the reason the largest effect on WSS is observed for tensile pre-loads (and not about zero pre-load) is probably because we consider a zero-to-peak loading regime, where the s(t) cycles about the mean bulk load of smax/2, as opposed to a fully reversing loading regime, where s(t) cycles about zero.
248
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 5. The dimensionless frictional energy dissipated per load cycle in the steady state, WSS, plotted against load factor, l, for 873 transient simulations at three coefficients of friction: (a) 0.3, (b) 0.6, and (c) 0.9.
5.5. Comparison of results €rkman and Klarbring (1987) apply the shakeIn their work, Bjo down limit optimization calculation to two example problems. They also perform a series of finite element transient simulations (until a steady state is reached) with a null initial displacement condition and determine the cyclic load above which shakedown does not occur, which can be regarded as an estimate of the shakedown limit.
Fig. 6. The optimal residual slip displacement distributions plotted in dimensionless form against normalized position along the contact, x/a.
These authors then find that this estimate of the shakedown limit gives similar results to the optimization calculation in most cases. This result could imply that l1~l2 for their example problems, so initial conditions are unimportant. Conversely, it may imply that a null initial displacement condition is a near optimal displacement condition for promoting shakedown.
Fig. 7. The effect of a (dimensionless) bulk pre-load, sdm, on the dimensionless energy dissipated per cycle in the steady state, WSS, for the example case of f ¼ 0.6 and l ¼ 3.
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
However, it seems unlikely for a null initial displacement distribution to be an optimal distribution for both examples problems €rkman and Klarbring (1987). Thus, this may be an considered by Bjo indication that l1~l2 for the two examples considered by these authors, and this would be in sharp contrast to the results shown in Figs. 4 and 5, which demonstrate that contacts may fail to shake down significantly below the shakedown limit. However, more transient simulations would need to be performed with different initial conditions to investigate the size of the conditional region for € rkman and Klarbring (1987). the contacts examined by Bjo 6. Discussion The aim of the analysis we have presented is to bring attention to a computationally efficient way to calculate the shakedown limit. This work has practical relevance because frictional contacts in engineering structures are very frequently subjected to a combination of static and cyclic loads resulting, at least transiently, in partial slip conditions. Although it is well known that partial-slip plays a significant role in determining component performance (Farris et al., 2000; Barber, 2011), it is not always clear if a given configuration will result in these conditions after several load cy€ rkman and cles. The Linear Programming approach presented by Bjo Klarbring (1987) provides an efficient way to determine the load range above which partial-slip is guaranteed. The desirability of partial-slip is application specific and generally involves a trade off between the positive effects of damping and the detrimental effects of fretting damage. Still, in many cases, the aim is to eliminate or minimize the extent of partial-slip and the fretting damage with which it is associated. For engineers with this aim, the results presented in Figs. 4 and 5 may be a bit disheartening; these figures illustrate that the dissipative properties of practical complete contact geometries, which invariably exhibit significant levels of coupling even if both bodies are elastically similar, are highly sensitive to initial conditions. Moreover, partial-slip can persist in the steady state even for load factors that are well below the shakedown limit. Although the calculation we have presented does not provide sufficient information to determine if a given load level is ‘safe’ (i.e. that shakedown is guaranteed), it does provide some very valuable information on how best to install components to encourage or discourage shakedown. This information is illustrated in Fig. 6, which shows the optimal residual slip displacement distributions for three coefficients of friction. Notice that this information is automatically obtained from the calculation of l2. Furthermore, the nature of the optimal slip distribution for this particular configuration suggests that shakedown can be encouraged by pre-loading the contact with bulk tension, and this is confirmed in Fig. 7. This suggests that for some configurations there may be a simple procedure for encouraging/inhibiting shakedown in practice. It is not clear how much can be inferred about other configurations from the small set of results we have considered here. This is meant both with regard to changes in contact geometry and changes in the types of loads that are applied, e.g. if we considered a shear load instead of a bulk load. However, one effect we expect will apply generally is that the size of the conditional region (see x 5.3) will increase as the level of coupling increases. This is simply because shakedown is never conditional for the uncoupled case, so one might expect the size of the conditional region to be related to how much coupling is present. This prompts the question as to whether the level of coupling can be quantified and used to predict the size of the conditional region. Some preliminary work has been done on this subject (Brake et al., September 2013), but more work is required to determine the relationship between metrics of frictional coupling and the size of the conditional region.
249
7. Conclusion We have revisited a technique for efficiently calculating the € rkman and Klarbring (1987) that can be shakedown limit due to Bjo applied to coupled, discrete frictional systems with a large number of degrees of freedom. This calculation is then applied to an example complete contact problem between a square punch and an elastically similar half-plane, and the influence of initial conditions on steady state energy dissipation is examined. To do this, we first discretize the problem using the finite element method and then apply a static reduction technique to reduce the problem to a discrete frictional system with 128 contact nodes. We then compare the calculated shakedown limit to a series of transient simulations that were run until a steady state was reached. These results verify that the shakedown limit calculation is correct. They also demonstrate that the dissipative properties of coupled complete contacts can be highly dependent on the initial residual slip displacement state. For instance, some transient simulations failed to shake down for load factors as low as 66%, 50%, and 62% the value of the shakedown limit when coefficient of friction was 0.3, 0.6, and 0.9, respectively.
Acknowledgements R. C. Flicek would like to thank Rolls-Royce plc and the Technology Strategy Board for financial support under the programme SILOET-II. The authors wish to thank Rolls-Royce plc for granting permission to publish this work. Appendix A. The finite element model We modelled the punch and the half-plane using ABAQUS/CAE version 6.11-1 as two-dimensional, plane-strain geometries, which incorporated a combination of CPE4R and CPE3 elements. The square punch was modelled with no external boundary conditions, and each of its edges were of length 2a. The model of the punch was comprised of 5316 nodes, and it had 128 uniformly spaced nodes along the contact interface. Because a true half-plane cannot be modelled using the finite element method, we approximated the half-plane as a large rectangular block of width 20a and height 10a. This block had two boundary conditions applied to it: (i) the displacement in the y-direction was constrained to be zero along the entirety of its bottom edge and (ii) the displacement in the xdirection was constrained to be zero at one node at the centre of the bottom edge of the block as shown in Fig. 1. The model of the block was comprised of 19,704 nodes. Like the punch, the block also had 128 uniformly spaced nodes along the contact interface, which was centred about its top edge. Thus the full model incorporated a total of 50,040 degrees of freedom with 256 degrees of freedom along the contact interface.
References Abad, J., Franco, J.M., Celorrio, R., Lez aun, L., 2012. Design of experiments and energy dissipation analysis for a contact mechanics 3d model of frictional bolted lap joints. Adv. Eng. Softw. 45 (1), 42e53. http://dx.doi.org/10.1016/ j.advengsoft.2011.09.021. Ahn, Y.J., Barber, J.R., October-November 2008. Response of frictional receding contact problems to cyclic loading. Int. J. Mech. Sci. 50 (10e11), 1519e1525. http://dx.doi.org/10.1016/j.ijmecsci.2008.08.003. Ahn, Y.J., Bertocchi, E., Barber, J.R., December 2008. Shakedown of coupled twodimensional discrete frictional systems. J. Mech. Phys. Solids 56 (12), 3433e3440. http://dx.doi.org/10.1016/j.jmps.2008.09.003. Banerjee, N., Hills, D.A., 2006. Analysis of stickeslip and contact-edge behaviour in a simplified fretting fatigue test. J. Strain Anal. Eng. Des. 41 (3), 183e192. http:// dx.doi.org/10.1243/03093247JSA83.
250
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Banerjee, N., Dini, D., Hills, D.A., 2008. Frictional complete contacts subject to shear and bulk tension. P. I. Mech. Eng. C J. Mech. 222 (12), 2301e2309. http:// dx.doi.org/10.1243/09544062JMES1076. Barber, J.R., 2011. Frictional systems subjected to oscillating loads. Ann. Solid Struct. Mech. 2 (2e4), 45e55. http://dx.doi.org/10.1007/s12356-011-0017-5. Barber, J.R., Klarbring, A., Ciavarella, M., January-February 2008. Shakedown in frictional contact problems for the continuum. C. R. Mec. 336 (1e2), 34e41. http://dx.doi.org/10.1016/j.crme.2007.10.013. € rkman, G., Klarbring, A., 1987. Shakedown and residual stresses in frictional Bjo systems. In: Gladwell, G.M.L., Ghonem, H., Kalousek, J. (Eds.), Contact Mechanics and Wear of Rail/Wheel Systems II: Proceedings of the 2nd International Symposium, pp. 27e39. Brake, M.R., Flicek, R.C., Hills, D.A., September 2013. Development of a coupling metric to assess the shakedown limits for a contact interface. In: 5th World Tribology Congress, Torino, Italy. Brayton, R.K., Director, S.W., Hachtel, G.D., Vidigal, L., 1979. A new algorithm for statistical circuit design based on quasi-newton methods and function splitting. IEEE Trans. Circuits Syst. 26 (9), 784e794. http://dx.doi.org/10.1109/ TCS.1979.1084701. Churchman, C.M., Hills, D.A., June 2006. General results for complete contacts subject to oscillatory shear. J. Mech. Phys. Solids 54 (6), 1186e1205. http:// dx.doi.org/10.1016/j.jmps.2005.12.005. Ciavarella, M., 1998. The generalized Cattaneo partial slip plane contact problem. I e Theory. Int. J. Solids Struct. 35 (18), 2349e2362. http://dx.doi.org/10.1016/ S0020-7683(97)00154-6. Dundurs, J., Stippes, M., 1970. Role of elastic constants in certain contact problems. J. Appl. Mech. 37 (4), 965e970. http://dx.doi.org/10.1115/1.3408725. Farris, T.N., Szolwinski, M.P., Harish, G., 2000. Fretting in aerospace structures and materials. In: Hoeppner, D.W., Chandrasekaran, V., Elliott, C.B. (Eds.), Fretting Fatigue: Current Technology and Practice, ASTM STP, vol. 1367. ASTM, West Conshohocken, PA, pp. 523e537. Hills, D.A., Dini, D., January 2004. What level of friction guarantees adhesion in a complete contact? J. Strain Anal. Eng. Des. 39 (5), 549e551. http://dx.doi.org/ 10.1243/0309324041896425. Hills, D.A., Thaitirarot, A., Barber, J.R., Dini, D., 2012. Correlation of fretting fatigue experimental results using an asymptotic approach. Int. J. Fatigue 43, 62e75. http://dx.doi.org/10.1016/j.ijfatigue.2012.02.006.
J€ ager, J., 1998. A new principle in contact mechanics. J. Tribol. 120 (4), 677e684. http://dx.doi.org/10.1115/1.2833765. Jang, Y.H., Barber, J.R., 2011. Frictional energy dissipation in materials containing cracks. J. Mech. Phys. Solids 59 (3), 583e594. http://dx.doi.org/10.1016/ j.jmps.2010.12.010. Klarbring, A., Ciavarella, M., Barber, J.R., December 2007. Shakedown in elastic contact problems with coulomb friction. Int. J. Solids Struct. 44 (25e26), 8355e8365. http://dx.doi.org/10.1016/j.ijsolstr.2007.06.013. Law, S.S., Wu, Z.M., Chan, S.L., 2006. Analytical model of a slotted bolted connection element and its behaviour under dynamic load. J. Sound. Vib. 292 (3), 777e787. http://dx.doi.org/10.1016/j.jsv.2005.09.028. Limmer, L., Nowell, D., Hills, D.A., 2001. A combined testing and modelling approach to the prediction of the fretting fatigue performance of splined shafts. P. I. Mech. Eng. G J. Aer. 215 (2), 105e112. http://dx.doi.org/10.1243/0954410011531808. Melan, E., 1936. Theorie statisch unbestimmter systeme aus ideal-plastischem baustoff. In: Sitzungsber. Akad. Wiss. Wien, Abt. 2A, 145, pp. 195e218. Mindlin, R.D., Deresiewicz, H., 1953. Elastic spheres in contact under varying oblique forces. J. Appl. Mech. 21, 327e344. Nowell, D., Hills, D.A., 1987. Mechanics of fretting fatigue tests. Int. J. Mech. Sci. 29 (5), 355e365. http://dx.doi.org/10.1016/0020-7403(87)90117-2. Powell, M.J.D., 1978. A fast algorithm for nonlinearly constrained optimization calculations. In: Numerical Analysis. Springer, pp. 144e157. http://dx.doi.org/ 10.1007/BFb0067703. Rajasekaran, R., Nowell, D., 2006. Fretting fatigue in dovetail blade roots: experiment and analysis. Tribol. Int. 39, 1277e1285. http://dx.doi.org/10.1016/ j.triboint.2006.02.044. Szolwinski, M.P., Farris, T.N., 1996. Mechanics of fretting fatigue crack formation. Wear 198 (1), 93e107. http://dx.doi.org/10.1016/0043-1648(96)06937-2. Thaitirarot, A., Flicek, R.C., Hills, D.A., Barber, J.R., 2014. The use of static reduction in the finite element solution of two-dimensional frictional contact problems. P. I. Mech. Eng. C J. Mech. 228 (9), 1474e1487. http://dx.doi.org/10.1177/ 0954406213509086. Xi, N.S., Zhong, P.D., Huang, H.Q., Yan, H., Tao, C.H., 2000. Failure investigation of blade and disk in first stage compressor. Eng. Fail. Anal. 7 (6), 385e392. http:// dx.doi.org/10.1016/S1350-6307(99)00045-X.
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
Determination of the shakedown limit for large, discrete frictional systems R.C. Flicek a, *, D.A. Hills a, J.R. Barber b, D. Dini c a
Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, United Kingdom Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109-2125, USA c Department of Mechanical Engineering, Imperial College London, South Kensington Campus, Exhibition Road, London SW7 2AZ, United Kingdom b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 25 March 2014 Accepted 3 August 2014 Available online 13 August 2014
It has been shown that complete frictional contacts subjected to cyclic loads can shake down such that the interface becomes fully stuck after some number of cycles (even if partial-slip occurs at the onset of loading). However, if the amplitude of the cyclic load is greater than a particular value (the shakedown limit), it is impossible for shakedown to occur. In this paper, we examine a numerical approach for determining the shakedown limit for elastic frictional systems subjected to quasi-static loads using a discrete formulation. We then use this technique to determine the shakedown limit for a finite element model of a (coupled) complete contact with ~50,000 total degrees of freedom and ~250 along the contact interface. Finally, we compare the calculated value of the shakedown limit to a series of over 1000 transient simulations and investigate the influence of initial conditions on steady state frictional energy dissipation. The results demonstrate that the dissipative properties of complete contacts can be highly dependent on the initial residual slip displacement state. © 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Shakedown limit Cyclic loading Fretting fatigue
1. Introduction It is extremely common for frictional contacts in engineering systems to be subjected to a combination of a static load (to form the contact) and a time varying cyclic load (e.g. due to vibrations). Some examples include riveted and bolted joints (Abad et al., 2012; Law et al., 2006), dove-tail connections in jet engines (Rajasekaran and Nowell, 2006; Xi et al., 2000), and spline couplings (Limmer et al., 2001; Banerjee and Hills, 2006; Banerjee et al., 2008) to name just a few. Furthermore, it is well known that even when the contact interface appears nominally to be stuck, these loading conditions frequently result in small zones of micro-slip (or partialslip) along the contact interface (where the relative tangential displacement is generally of the order 25100 mm (Szolwinski and Farris, 1996)). The effect of micro-slip on system performance is two-fold: on the one hand, it results in the well known damage process of fretting fatigue, which significantly reduces the service life of components (Farris et al., 2000); but on the other hand, it provides a source of damping, which can be beneficial for system performance (Law et al., 2006).
* Corresponding author. E-mail address: robert.fl[email protected] (R.C. Flicek). http://dx.doi.org/10.1016/j.euromechsol.2014.08.001 0997-7538/© 2014 Elsevier Masson SAS. All rights reserved.
Much of the previous work on micro-slip and fretting fatigue has focused on incomplete (or non-conforming) contacts (e.g. (Nowell and Hills, 1987; Hills et al., 2012)), which are characterized by the edge of contact being smooth such that the bodies do not ‘conform’ in the absence of applied loads. This is probably due to their relevance in engineering practice but also because closed form solutions exist for many relevant loading cases, e.g. (Mindlin €ger, 1998; Ciavarella, 1998). For these and Deresiewicz, 1953; Ja contacts, the contact area is a function of the applied normal load, and the contact pressure distribution smoothly tends to zero as the contact edge is approached. As a direct consequence of this, a zone of micro-slip will almost always arise at the edge of an incomplete contact that is subjected to cyclic loads (Barber et al., JanuaryeFebruary 2008). Complete (or conforming) contacts, such as the contact shown in Fig. 1, are characterized by the two contacting bodies ‘conforming’ in the absence of the applied loads. For this class of contacts, the edge of contact is sharp, and the contact area either does not change in the presence of applied loads or becomes smaller (due to the edges ‘lifting off’).1 Because sharp features result in a large
1 For this reason, this type of contact is sometimes referred to as a receding contact (Dundurs and Stippes, 1970).
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 1. The complete contact formed between a square elastic punch and an elastically similar half-plane showing: the static normal load, P, the time-varying load, s(t), the size of the contact, 2a, and the (x,y) coordinate set.
stress concentration in their vicinity, complete contacts are usually avoided when possible. However, these contacts do arise in some critical applications such as in the spline couplings that join the central shafts in jet engines (Limmer et al., 2001; Banerjee and Hills, 2006). Numerical analysis of complete contacts subjected to cyclic loads sometimes predicts that the contact shakes down (Churchman and Hills, June 2006; Banerjee and Hills, 2006): that is, after the first few loading cycles, a favourable residual displacement state is developed that inhibits slip for all subsequent cycles. This behaviour clearly has some similarities to cyclic plasticity problems, which has prompted the question as to whether a theorem analogous to Melan's plastic shakedown theorem (Melan, 1936) could be developed for frictional contacts. In other words, if there exists a residual displacement distribution that would inhibit slip at all times during load cycle (i.e. a safe displacement distribution), is the existence of this displacement distribution sufficient to guarantee that the contact will shake down? This question has been answered unequivocally by Klarbring et al. (December 2007) for discrete contact problems and by Barber et al. (JanuaryeFebruary 2008) for continuum problems. These authors have proven that the existence of a safe displacement distribution is both a necessary and sufficient condition to guarantee shakedown if the contact is uncoupled: that is, if relative tangential contact displacements have no influence on the distribution of normal contact traction. Thus, for uncoupled contacts, there exists a load factor, i.e. a ratio of cyclic load amplitude to static load (e.g. s/P in Fig. 1), above which shakedown is impossible and below which it is guaranteed; hence the contact's steady state response never depends on initial conditions. On the other hand, if the contact is coupled, the existence of a safe displacement distribution is still a necessary condition for shakedown to be possible, but it is no longer sufficient. In fact, for coupled contacts, there almost always exists a range of load factors for which the steady state response depends on the initial residual displacement distribution (Klarbring et al., December 2007; Barber et al., JanuaryeFebruary 2008). Ahn et al. (December 2008) have demonstrated this behaviour for a two-node system and have shown that for coupled contacts there exists: 1. a load factor below which shakedown is guaranteed to occur, denoted l1 2. a load factor above which shakedown is impossible, denoted l2. Note that these two limits on the applied load may be zero for many contacts, i.e. l1 ¼ l2 ¼ 0, (e.g. for incomplete contacts), implying that shakedown is impossible. Moreover, even if l2 > 0, it
243
is still possible that l1 ¼ 0, implying that although shakedown is sometimes possible, whether it actually occurs is always dependent on initial conditions. Ahn et al. (December 2008) also presented a method for calculating both l1 and l2. Unfortunately, their approach involves solving a system of equations that becomes combinatorially more complex as the number of degrees of freedom in the system increases. For instance, Jang and Barber (2011) were able to apply this approach to determine both l1,l2 for a system of 10 cracks, but applying this analysis to a significantly larger frictional system would be computationally prohibitive. As it is common for the finite element models used in practice to incorporate one or two orders of magnitude more nodes along the frictional interface than are considered by Jang and Barber (2011), it is of practical interest to develop a more computationally efficient approach. € rkman and Klarbring (1987) have presented an efficient Bjo approach for calculating l2, which we will henceforth refer to as the shakedown limit. To do this, these authors formulate the problem in discrete form, and show that the shakedown limit calculation can be posed as a Linear Programming problem for which several efficient solution algorithms exist. These authors then apply this calculation to two example problems, and also compare the results to numerical results obtained using the finite element method. The aim of this paper is: (i) to re-visit the Linear Programming approach for calculating the shakedown limit and (ii) to examine the influence of initial conditions on the level of frictional dissipation present in the steady state. This is of practical importance because energy dissipation is associated with damage processes such as fretting fatigue and wear, which are themselves linked to component life. In particular, we wish to investigate whether energy dissipation is significantly influenced by initial conditions or if their effect can safely be neglected. The structure of the paper is as follows. In x 2, we formulate the discrete contact problem. In x 3, we outline how to determine when the onset of slip or separation will occur for a contact subjected to monotonic loading. In x 4, we describe how the shakedown limit calculation can be posed as a Linear Programming optimization problem for which several solution algorithms already exist. In x 5, these calculations are applied to the example (coupled) complete contact shown in Fig. 1. Over 1000 numerical simulations of the time-evolution of the contact are then performed under a wide range of cyclic loads and initial displacement states, and the influence of initial conditions on steady state frictional energy dissipation is investigated. Finally, the implications of these results are discussed in x 6.
2. Formulation Let us consider some two-dimensional complete contact geometry, such as that shown in Fig. 1, which has been discretized using the finite element method such that there are N nodes in potential frictional contact along the interface. We can write the reaction forces, r, and the relative displacements, u, at the contact nodes as
r ¼ ½q; pT ¼ ½q1 ; …; qN ; p1 ; …; pN T
(1a)
u ¼ ½v; wT ¼ ½v1 ; …; vN ; w1 ; …; wN T ;
(1b)
where qi,pi are the shear, normal reactions and vi,wi are the tangential, normal relative displacements, and where i2{1,…,N}. Here we adopt the convention that pi is positive in compression (the negative y-direction), wi positive for a positive gap (the
244
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
positive y-direction), and qi,vi are positive in the positive x-direction with respect to the x,y coordinate system shown in Fig. 1. We can write an expression for the time-evolution of the contact reactions as
rðtÞ ¼ KuðtÞ þ r w ðtÞ;
(2)
where K is a symmetric 2N 2N contact stiffness matrix and rw accounts for the influence of the external loads on the contact reactions. More specifically, rw represents the reactions, r, that would be generated by the external loads if all the contact displacements were fixed to be zero, i.e. if u ¼ 0. Note that both K and rw can be obtained from a standard finite element model using the static reduction procedure described by Thaitirarot et al. (2014). Also note that the contact stiffness matrix, K, is positive semi-definite (and singular) if either contacting body has rigid body degrees of freedom; otherwise K is positive definite (and non-singular).
qðtÞ pðtÞ
¼
A B
BT C
vðtÞ wðtÞ
þ Fs
qws pws
þ Fc ðtÞ
qwc pwc
;
(5)
where q,p,v,w are defined in equation (1), and K is partitioned into sub-matrices corresponding to the following forceedisplacement relationships: tangentialetangential, A; normaletangential, B; and normalenormal, C. Note that as K is symmetric and positive semidefinite (or positive definite), so too are A and C, whereas B is generally not symmetric. We also define the load factor, l, as
l¼
Fc ðtÞjmax Fc ðtÞjmin Fcmax Fcmin ¼ ; Fs Fs
(6)
where Fcmin , Fcmax are the minimum, maximum values of the cyclic load, respectively. Note that the load factor is only defined for cyclic loading regimes. 3. First violation of the stick condition
2.1. The Coulomb friction law Before defining the friction law, we note that the contacting bodies cannot inter-penetrate and that the normal contact reactions cannot be tensile:
wi 0;
pi 0:
(3a)
Furthermore, if there is a positive gap at node i, its normal and shear reactions must be zero. Conversely, if the normal reaction at node i is compressive, its gap must be zero:
wi > 00pi ¼ qi ¼ 0
(3b)
pi > 00wi ¼ 0:
(3c)
Finally, we assume that the Coulomb friction condition applies at each node, and this requires that
jqi j fpi ;
The load ratio, Fc/Fs, at which slip or separation first occurs under monotonic loading is always dependent on the initial residual slip displacement distribution, v(0). Nevertheless, this information may be helpful when planning or analysing experiments, particularly when the case of a null initial slip distribution, i.e. v(0) ¼ 0, is considered. Moreover, this calculation is trivially simple to carry out. Thus, suppose we have a complete contact with v(0) ¼ 0 subject to a loading regime in which the static load, Fs, is first applied and held constant, and the time-varying load, Fc, is subsequently applied monotonically. We can determine the load ratio at which slip or separation first occurs by checking when the stick condition is first violated. In other words, we search for the lowest value of Fc/Fs that violates the condition pi 0 or jqij < fpi at any node, i. Recall that for complete contacts, w ¼ 0 in the absence of external loads by definition. Thus we can use equation (5) to write these conditions as
pws i þ
Fc wc p 0 Fs i
(7a)
Fc wc Fc wc qi < f pws p : þ i Fs Fs i
(7b)
f >0
(3d)
jqi j < fpi 0v_i ¼ 0
(3e)
qws i þ
0 < jqi j ¼ fpi 0sgnðv_i Þ ¼ sgnðqi Þ;
(3f)
Note that these conditions are only valid for predicting first violations, so once slip, for example, is predicted to initiate from equation (7b), separation is no longer accurately predicted from equation (7a).
where f is the coefficient of friction and a superposed dot denotes the derivative with respect to time, t. Thus equation (3d) defines the admissible range of shear reactions; equation (3e) states that if the magnitude of qi is less than a critical value, node i is stationary; and equation (3f) states that if jqij is equal to a critical value, slip occurs in the direction that opposes the shear reaction at that node.
2.2. Loading regime We consider a loading regime comprised of a static component, Fs, superposed with a quasi-statically applied time-varying (cyclic) component, Fc(t), where Fs,Fc are scalars that account for the strength of the applied load. The effect of these applied loads on contact reactions can be expressed as
rw ðtÞ ¼ Fs r ws þ Fc ðtÞr wc ;
(4)
where rws,rwc are the rw vectors corresponding to a unit application of Fs,Fc, respectively. With this notation, equation (2) can be rewritten as
4. The shakedown limit It is perhaps most natural to think about frictional shakedown as it might arise in the context of an experiment: namely, there is some test configuration that is subjected to a particular combination of loads, and the coefficient of friction (which we will assume remains constant) is determined by the material properties, surface treatment, etc. Thus if we wanted to determine the shakedown limit experimentally, we might apply cyclic loads of various amplitudes, vary the initial displacement state, and use digital image correlation to determine the maximum cyclic load amplitude for which the contact shakes down. In contrast, the approach we describe here for calculating the shakedown limit is, in some sense, the reverse of this hypothetical experimental approach. To pose the calculation of the shakedown limit as a Linear €rkman and Klarbring, 1987), we Programming problem as in (Bjo begin by choosing a particular cyclic loading regime in which the cyclic load ranges from some minimum value, Fcmin , to some
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
maximum value, Fcmax . We then assume that the load factor (see equation (6)) corresponding to this load range is the shakedown limit if the coefficient of friction is some currently unknown value, denoted fmin. Our aim therefore becomes to determine fmin: that is, the coefficient of friction below which it is impossible for the contact to shake down for the specified load range. We can then perform this calculation for a series of load ranges to determine the shakedown limit for a range of friction coefficients. Note that if the load range is sufficiently large, it may be impossible for the contact to shake down irrespective of the coefficient of friction, making our initial assumption false. This can occur if there does not exist any residual slip displacement that inhibits separation throughout the load cycle. In this case, the calculation will simply fail to converge, implying shakedown is impossible for the specified load range.
To determine fmin, first recall that if a contact is to shake down, separation cannot occur at any time in the load cycle, i.e. w(t) ¼ 0. Therefore, for a contact that has shaken down, we can use equation (5) to write the ratio of the (element-wise) absolute value of the shear reactions to the normal reactions, R, as
Rðv; Fc Þ ¼
Note that whereas we have formulated the problem as a €rkman and ‘minimax’ constrained optimization problem, Bjo Klarbring (1987) formulate it as a constrained maximization. However, both approaches are equivalent, and the reader is € rkman and Klarbring, 1987) for a directed to the original paper (Bjo discussion of the mathematical structure of this optimization. 4.2. Improving convergence One additional constraint that is often helpful for improving convergence of the optimization in equation (11) is to eliminate the potential for rigid body motion in the slip displacement distribution. This can be done by removing one node's slip displacement from the optimization and optimizing against a subset of the slip displacements, v*, where
v ¼ ½v1 ; …; vN1 T :
4.1. Calculation
jqðFc Þj jAv þ Fs qws þ Fc qwc j ¼ : pðFc Þ Bv þ Fs pws þ Fc pwc
(8)
Notice that we have written the contact reactions as a function of Fc as opposed to a strict function of time, t. We can do this because every load cycle for a shaken down contact is identical, and the reactions simply vary with the applied loads. Also notice that v is not a function of Fc because, by definition, slip cannot occur at any point during the load cycle for a contact that has shaken down. According to the friction law defined in x 2.1, slip will occur if any element of R becomes equal to the coefficient of friction. As we are considering the case of quasi-static loading, R will take on its maximum value at one (or both) of the extreme points (in time) of the load cycle, i.e. at Fcmin or Fcmax. We can construct a vector, S(v), of length 2N that is comprised of the contact reaction ratios, R, at the extreme points in the load cycle as
245
(12)
The slip displacement vN, which is not optimized against, can then be set to
vN ¼
N1 X
vi
(13)
i¼1
to force the sum of all slip displacements to be zero. Note that it is not important which node's slip displacement is removed from the optimization.2 5. Example application
(11)
We now intend to illustrate how these calculations can be used in practice by applying them to the example problem shown in Fig. 1. This figure shows the complete contact formed between two elastic continuum bodies: a square punch and an elastically similar half-plane. As our approach uses the discrete formulation, the first step in our analysis is to discretize the contact, which we have done using the commercial finite element software ABAQUS/CAE. Of course once we discretize these bodies, we are no longer studying the same problem in a rigorous sense. However, in practice it is usually assumed that a continuum is well represented by a discrete system with a sufficiently refined mesh. The details of the finite element model we have used are in Appendix A, but here it suffices to state that both bodies have elastic modulus E ¼ 200 GPa, Poission's ratio n ¼ 0.3, the model assumes plane strain with unit depth, the contact interface comprises 128 uniformly spaced nodes, and the full model incorporating the punch and the half-plane comprises 50,040 degrees of freedom, where every node in the model is associated with two degrees of freedom. Note that the finite element model we used does not take advantage of the symmetry about the vertical centreline of the contact. This is simply because we generated these results from a set of models we created to study several different loading regimes, and some of these loads do not maintain this horizontal symmetry, e.g. when a uniform shear traction is applied to the top of the punch. To determine the sensitivity of the contact solution to the discretization, we performed a mesh convergence study in which we varied both the number of contact nodes and the total number of nodes in the model. The results of this study suggest that a model with 128 contact nodes provides a converged result at the interior
where j2{1,…,2N}. Several algorithms exist for the solution of such problems, e.g. (Brayton et al., 1979; Powell, 1978). Note that the solution to this problem, vopt, generally is not unique.
2 We tested this by removing nodes at various positions along the contact from the optimization and specifying their slip displacement according to equation (13); the solution was insensitive to which node was chosen.
) ( R v; Fcmin SðvÞ ¼
: R v; Fcmax
(9)
Thus further slip will be inhibited throughout the load cycle if
f SðvÞjmax ;
(10)
i.e. if the coefficient of friction is greater than the maximum element value of S. To determine fmin, we simply need to optimize v to determine an optimal slip distribution, vopt. For our purposes, an optimal slip distribution is one that: (i) prevents separation at all times in the load cycle and (ii) requires the lowest coefficient of friction to prevent slip from occurring. From equation (10), we see that (ii) requires that vopt minimizes the maximum value of S while (i) requires that the condition p 0 is simultaneously satisfied. Therefore, finding vopt is a standard ‘mini-max’ optimization problem of the form
min max Sj ðvÞ v
j
such that p 0;
246
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
of the contact but that the solution close to the corner is not fully converged. However, we found that obtaining a converged result at the corners with a uniform mesh along the contact would be computationally impractical, so a more complex biased mesh would be required. As our primary object is to demonstrate the technique, we feel that this simple mesh is sufficient. However, note that the technique presented in x 4.1 can be applied to larger problems without modification. The next step in this analysis is to apply the static reduction technique described by Thaitirarot et al. (2014), reducing the finite element model to a discrete system of the N contact nodes alone. In doing this, we obtain the contact stiffness matrix, K, and the rw vectors corresponding to each applied load, P,s(t). Note that P corresponds to the mean contact pressure, and s(t) corresponds to the mean bulk load; hence both loads have units of [FL2], where F,L denote force, length, respectively. Thus we can write the applied loads for this example problem in the form given by equation (4) as
rw ðtÞ ¼ PrP þ sðtÞr s ;
(14)
where rP,rs have units of [F/(FL2)] or simply [L2] and correspond to the application of a unit load in the directions shown in Fig. 1. Finally, the load factor is defined as
l¼
sðtÞjmax sðtÞjmin smax smin ¼ ; P P
that is implied for monotonically increasing/decreasing s, which corresponds to a bulk tension/compression, respectively. Note that we adopt the convention that leading edge slip refers to the edges of the punch slipping outwards, whereas trailing edge slip refers to the edges of the punch slipping inwards. The first thing to notice from this figure is that if f < 0.48, the first violation of the stick condition occurs before s is applied; that is, leading edge slip initiates from the contact edge on the application of P alone. It is worth mentioning that an asymptotic analysis of the contact edge in the continuum formulation predicts that edge slip will occur if f < 0.543 (see (Hills and Dini, January 2004)). The discrepancy between these results is due to the mesh being insufficiently fine near the contact edges. In any case, when f > 0.48, the entirety of the contact interface remains stuck on the application of P. If s is then monotonically applied in tension, slip initiates in the trailing edge sense when s/ P~2, depending on the coefficient of friction. Conversely, if s is monotonically applied in compression, slip initiates in the leading edge sense within the range 2.3 < s/P < 0. Note that although slip initiates somewhat inwards from the contact edge when f > 0.48 (in both tension and compression), the point of first slip always remains quite close to the contact edge. Also notice that for f < 1, the first violation of the stick condition always results from slip, not separation.
(15)
where smin and smax are the minimum and maximum values of s, respectively. 5.1. First violation of the stick condition We can now determine the load ratio, s/P, at which slip or separation first occurs for this contact when v(0) ¼ 0. To do this, we assume that P is first raised to some value and then held constant, and that s(t) is subsequently applied monotonically (in tension or compression). We then check for violations of the conditions specified in equation (7), and the results of this calculation are shown in Fig. 2 on a plot of s/P vs. f. The way to interpret this figure is to pick a coefficient of friction, which corresponds to some vertical line, and to view the behaviour
5.2. The shakedown limit We now intend to calculate the shakedown limit, l2, and to compare it with the results of a series of transient simulations. The loading history we consider is shown in Fig. 3: namely, P is first applied and held constant, and s(t) is then varied within the range 0 s(t) smax. In other words, we consider the case of ‘zero-topeak’ loading in cyclic tension; thus the load factor is defined as l ¼ smax/P. Also note that the initial slip displacement condition, i.e. v(0), need not be null; indeed, the transient simulations we carried out explored a wide range of v(0) distributions.3 We calculated l2 for the loading regime described above in MATLAB R2013a using the fminimax function, and the results are shown in Fig. 4 on a plot of l vs. f. Also shown in this figure are the results of 347 ‘marching in time’ transient simulations,4 which we carried out using the Gauss-Seidel algorithm described by Ahn and Barber (2008). Simulations that eventually shook down are marked with a circle, and those that failed to shake down are marked with a dot. The main thing to recognize about this figure is that it demonstrates that the calculated values of l2 are consistent with results of all the transient simulations we performed; that is, every transient simulation for l > l2 failed to shake down. 5.3. Frictional energy dissipation All the transient simulations we performed were subjected to the loading regime shown in Fig. 3 until a steady state response was reached, which we determined using a criterion based on the frictional energy dissipation per cycle, W. Specifically, the contact was determined to have reached a steady state response if: (i) W was similar in size to the smallest meaningful value of W that can be represented with a double precision floating point number or
Fig. 2. The slip-stick behaviour of the contact shown in Fig. 1 when P is first applied to a contact with a null initial slip displacement distribution, and s is then applied monotonically in tension (positive) or compression (negative).
3 Note that we varied v(0) by assuming a particular displacement distribution and then simply scaled this distribution. This approach enables any distribution to be used, but we mostly assumed either a linear or sinusoidal distribution or the optimal distribution obtained from the calculation of l2. 4 It is difficult to identify some of these simulations in this and subsequent figures since we ran multiple simulations with the same l but with different v(0).
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 3. The cyclic loading regime considered for the shakedown limit calculation.
(ii) the change in dissipation from the previous load cycle was less than 0.01%. We calculated the frictional energy dissipation per cycle as
W¼
I N X i¼1
display a value of zero; thus, we plot simulations that shook down (for which WSS ¼ 0) as the arbitrarily small value of 108. Both Figs. 4 and 5 demonstrate that the range of l for which shakedown is conditional (i.e. where some simulations shook down while others did not) can be quite large. For example, Fig. 5 shows that simulations failed to shake down when l was as low as 66%, 50%, and 62% the value of l2 when f was 0.3, 0.6, and 0.9, respectively. In addition, Figs. 4 and 5 suggest that the size of the conditional region (as a percentage of l2) remains quite similar over the range of f we have considered. However, it is almost certain that the conditional region is at least somewhat larger than is suggested by these figures because it is unlikely that any of the v(0) distributions we tested were the absolute ‘worst’ v(0) possible (in terms of encouraging shakedown). 5.4. Optimal initial residual slip displacement
qi ðtÞv_i ðtÞdt;
(16)
cycle
where the time-integration was performed numerically with a piecewise-linear approximation in each load increment. Note that W is a non-negative quantity because equation (3) states that the sign of qi always opposes that of v_i for v_i s0. We performed a convergence study and found a load increment of 0.05 (with respect to s) was sufficient to obtain a converged result; thus this load increment was used for all the transient simulations we performed. As this model is two-dimensional, W is the dissipation per load cycle per unit depth; hence its units are [(FL)/L] or simply [F]. Thus we can write the dissipation per cycle in the steady state in dimensionless form, WSS, which is given by
WSS ¼
247
E
W; a2 P 2 1 n2
(17)
where 2a is the contact-width. To illustrate the effect that initial conditions can have on WSS, we plot WSS against l in Fig. 5 for 873 transient simulations carried out at three coefficients of friction: 0.3, 0.6, and 0.9. Each of these plots also includes a vertical line corresponding to l2 at the relevant coefficient of friction. Note that as the vertical scale on these plots is logarithmic, they cannot
Fig. 4. The (calculated) shakedown limit and the results of 347 transient simulations that were run until a steady state was reached. Transient simulations that shook down are shown with a circle, and those that did not shake down are shown with a dot.
Fig. 5 illustrates that the frictional energy dissipated by a complete contact after several load cycles can be very sensitive to the initial residual slip displacement state, v(0). Hence, it is of practical interest to determine which sorts of initial displacement conditions result in a more/less dissipative contact. Conveniently, we automatically obtain an optimal slip distribution for promoting shakedown when calculating l2. Fig. 6 shows three of these distributions in dimensionless form plotted against normalized position along the contact interface, x/a, for f ¼ 0.3,0.6,0.9, where the dimensionless slip displacement is given by
vdm ¼
E
v: aP 1 n2
(18)
Although the information provided by Fig. 6 is useful, it is clear that pre-loading a contact with distributions such as these might be difficult in practice. However, a linear distribution of v(0) can be achieved by the simple expedient of applying an appropriate bulk stress, s, before the constant normal load, P, is applied. Hence we might ask: what is the best/worst initial condition for promoting shakedown that can be achieved by this procedure? To address this question, we performed a series of simulations with different bulk pre-loads, and the results of these simulations are shown in Fig. 4. In addition, we show the effect of a (dimensionless) bulk pre-load, sdm, on steady state dissipation, WSS, in Fig. 7 for the example case of f ¼ 0.6 and l ¼ 3, where the dimensionless bulk load is given by sdm ¼ s/P. Fig. 7 illustrates that shakedown can be encouraged by applying a tensile bulk pre-load prior to the formation of the contact, or it can be inhibited with a compressive bulk pre-load. In fact, even though l ¼ 3 is 97% the value of l2 for f ¼ 0.6, Fig. 7 demonstrates that shakedown can still be achieved by applying a bulk pre-load of sdm 1.5. Although it is perhaps surprising how large a range of WSS can be excited using a bulk pre-load alone, it should not be surprising that tensile/compressive bulk pre-loads promote/inhibit shakedown. This is because all the optimal slip distributions shown in Fig. 6 consist mainly of trailing edge slip, and trailing/leading edge slip is generated by a tensile/compressive bulk load. Also notice that the effect of the bulk pre-load on WSS ‘saturates’ quite quickly and that the largest effect occurs near sdm~1. This saturation effect is not surprising because if excessively large bulk pre-loads are applied, there is simply a larger amount of slip in the first cycle, but the steady state remains largely unchanged. Also, note that the reason the largest effect on WSS is observed for tensile pre-loads (and not about zero pre-load) is probably because we consider a zero-to-peak loading regime, where the s(t) cycles about the mean bulk load of smax/2, as opposed to a fully reversing loading regime, where s(t) cycles about zero.
248
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Fig. 5. The dimensionless frictional energy dissipated per load cycle in the steady state, WSS, plotted against load factor, l, for 873 transient simulations at three coefficients of friction: (a) 0.3, (b) 0.6, and (c) 0.9.
5.5. Comparison of results €rkman and Klarbring (1987) apply the shakeIn their work, Bjo down limit optimization calculation to two example problems. They also perform a series of finite element transient simulations (until a steady state is reached) with a null initial displacement condition and determine the cyclic load above which shakedown does not occur, which can be regarded as an estimate of the shakedown limit.
Fig. 6. The optimal residual slip displacement distributions plotted in dimensionless form against normalized position along the contact, x/a.
These authors then find that this estimate of the shakedown limit gives similar results to the optimization calculation in most cases. This result could imply that l1~l2 for their example problems, so initial conditions are unimportant. Conversely, it may imply that a null initial displacement condition is a near optimal displacement condition for promoting shakedown.
Fig. 7. The effect of a (dimensionless) bulk pre-load, sdm, on the dimensionless energy dissipated per cycle in the steady state, WSS, for the example case of f ¼ 0.6 and l ¼ 3.
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
However, it seems unlikely for a null initial displacement distribution to be an optimal distribution for both examples problems €rkman and Klarbring (1987). Thus, this may be an considered by Bjo indication that l1~l2 for the two examples considered by these authors, and this would be in sharp contrast to the results shown in Figs. 4 and 5, which demonstrate that contacts may fail to shake down significantly below the shakedown limit. However, more transient simulations would need to be performed with different initial conditions to investigate the size of the conditional region for € rkman and Klarbring (1987). the contacts examined by Bjo 6. Discussion The aim of the analysis we have presented is to bring attention to a computationally efficient way to calculate the shakedown limit. This work has practical relevance because frictional contacts in engineering structures are very frequently subjected to a combination of static and cyclic loads resulting, at least transiently, in partial slip conditions. Although it is well known that partial-slip plays a significant role in determining component performance (Farris et al., 2000; Barber, 2011), it is not always clear if a given configuration will result in these conditions after several load cy€ rkman and cles. The Linear Programming approach presented by Bjo Klarbring (1987) provides an efficient way to determine the load range above which partial-slip is guaranteed. The desirability of partial-slip is application specific and generally involves a trade off between the positive effects of damping and the detrimental effects of fretting damage. Still, in many cases, the aim is to eliminate or minimize the extent of partial-slip and the fretting damage with which it is associated. For engineers with this aim, the results presented in Figs. 4 and 5 may be a bit disheartening; these figures illustrate that the dissipative properties of practical complete contact geometries, which invariably exhibit significant levels of coupling even if both bodies are elastically similar, are highly sensitive to initial conditions. Moreover, partial-slip can persist in the steady state even for load factors that are well below the shakedown limit. Although the calculation we have presented does not provide sufficient information to determine if a given load level is ‘safe’ (i.e. that shakedown is guaranteed), it does provide some very valuable information on how best to install components to encourage or discourage shakedown. This information is illustrated in Fig. 6, which shows the optimal residual slip displacement distributions for three coefficients of friction. Notice that this information is automatically obtained from the calculation of l2. Furthermore, the nature of the optimal slip distribution for this particular configuration suggests that shakedown can be encouraged by pre-loading the contact with bulk tension, and this is confirmed in Fig. 7. This suggests that for some configurations there may be a simple procedure for encouraging/inhibiting shakedown in practice. It is not clear how much can be inferred about other configurations from the small set of results we have considered here. This is meant both with regard to changes in contact geometry and changes in the types of loads that are applied, e.g. if we considered a shear load instead of a bulk load. However, one effect we expect will apply generally is that the size of the conditional region (see x 5.3) will increase as the level of coupling increases. This is simply because shakedown is never conditional for the uncoupled case, so one might expect the size of the conditional region to be related to how much coupling is present. This prompts the question as to whether the level of coupling can be quantified and used to predict the size of the conditional region. Some preliminary work has been done on this subject (Brake et al., September 2013), but more work is required to determine the relationship between metrics of frictional coupling and the size of the conditional region.
249
7. Conclusion We have revisited a technique for efficiently calculating the € rkman and Klarbring (1987) that can be shakedown limit due to Bjo applied to coupled, discrete frictional systems with a large number of degrees of freedom. This calculation is then applied to an example complete contact problem between a square punch and an elastically similar half-plane, and the influence of initial conditions on steady state energy dissipation is examined. To do this, we first discretize the problem using the finite element method and then apply a static reduction technique to reduce the problem to a discrete frictional system with 128 contact nodes. We then compare the calculated shakedown limit to a series of transient simulations that were run until a steady state was reached. These results verify that the shakedown limit calculation is correct. They also demonstrate that the dissipative properties of coupled complete contacts can be highly dependent on the initial residual slip displacement state. For instance, some transient simulations failed to shake down for load factors as low as 66%, 50%, and 62% the value of the shakedown limit when coefficient of friction was 0.3, 0.6, and 0.9, respectively.
Acknowledgements R. C. Flicek would like to thank Rolls-Royce plc and the Technology Strategy Board for financial support under the programme SILOET-II. The authors wish to thank Rolls-Royce plc for granting permission to publish this work. Appendix A. The finite element model We modelled the punch and the half-plane using ABAQUS/CAE version 6.11-1 as two-dimensional, plane-strain geometries, which incorporated a combination of CPE4R and CPE3 elements. The square punch was modelled with no external boundary conditions, and each of its edges were of length 2a. The model of the punch was comprised of 5316 nodes, and it had 128 uniformly spaced nodes along the contact interface. Because a true half-plane cannot be modelled using the finite element method, we approximated the half-plane as a large rectangular block of width 20a and height 10a. This block had two boundary conditions applied to it: (i) the displacement in the y-direction was constrained to be zero along the entirety of its bottom edge and (ii) the displacement in the xdirection was constrained to be zero at one node at the centre of the bottom edge of the block as shown in Fig. 1. The model of the block was comprised of 19,704 nodes. Like the punch, the block also had 128 uniformly spaced nodes along the contact interface, which was centred about its top edge. Thus the full model incorporated a total of 50,040 degrees of freedom with 256 degrees of freedom along the contact interface.
References Abad, J., Franco, J.M., Celorrio, R., Lez aun, L., 2012. Design of experiments and energy dissipation analysis for a contact mechanics 3d model of frictional bolted lap joints. Adv. Eng. Softw. 45 (1), 42e53. http://dx.doi.org/10.1016/ j.advengsoft.2011.09.021. Ahn, Y.J., Barber, J.R., October-November 2008. Response of frictional receding contact problems to cyclic loading. Int. J. Mech. Sci. 50 (10e11), 1519e1525. http://dx.doi.org/10.1016/j.ijmecsci.2008.08.003. Ahn, Y.J., Bertocchi, E., Barber, J.R., December 2008. Shakedown of coupled twodimensional discrete frictional systems. J. Mech. Phys. Solids 56 (12), 3433e3440. http://dx.doi.org/10.1016/j.jmps.2008.09.003. Banerjee, N., Hills, D.A., 2006. Analysis of stickeslip and contact-edge behaviour in a simplified fretting fatigue test. J. Strain Anal. Eng. Des. 41 (3), 183e192. http:// dx.doi.org/10.1243/03093247JSA83.
250
R.C. Flicek et al. / European Journal of Mechanics A/Solids 49 (2015) 242e250
Banerjee, N., Dini, D., Hills, D.A., 2008. Frictional complete contacts subject to shear and bulk tension. P. I. Mech. Eng. C J. Mech. 222 (12), 2301e2309. http:// dx.doi.org/10.1243/09544062JMES1076. Barber, J.R., 2011. Frictional systems subjected to oscillating loads. Ann. Solid Struct. Mech. 2 (2e4), 45e55. http://dx.doi.org/10.1007/s12356-011-0017-5. Barber, J.R., Klarbring, A., Ciavarella, M., January-February 2008. Shakedown in frictional contact problems for the continuum. C. R. Mec. 336 (1e2), 34e41. http://dx.doi.org/10.1016/j.crme.2007.10.013. € rkman, G., Klarbring, A., 1987. Shakedown and residual stresses in frictional Bjo systems. In: Gladwell, G.M.L., Ghonem, H., Kalousek, J. (Eds.), Contact Mechanics and Wear of Rail/Wheel Systems II: Proceedings of the 2nd International Symposium, pp. 27e39. Brake, M.R., Flicek, R.C., Hills, D.A., September 2013. Development of a coupling metric to assess the shakedown limits for a contact interface. In: 5th World Tribology Congress, Torino, Italy. Brayton, R.K., Director, S.W., Hachtel, G.D., Vidigal, L., 1979. A new algorithm for statistical circuit design based on quasi-newton methods and function splitting. IEEE Trans. Circuits Syst. 26 (9), 784e794. http://dx.doi.org/10.1109/ TCS.1979.1084701. Churchman, C.M., Hills, D.A., June 2006. General results for complete contacts subject to oscillatory shear. J. Mech. Phys. Solids 54 (6), 1186e1205. http:// dx.doi.org/10.1016/j.jmps.2005.12.005. Ciavarella, M., 1998. The generalized Cattaneo partial slip plane contact problem. I e Theory. Int. J. Solids Struct. 35 (18), 2349e2362. http://dx.doi.org/10.1016/ S0020-7683(97)00154-6. Dundurs, J., Stippes, M., 1970. Role of elastic constants in certain contact problems. J. Appl. Mech. 37 (4), 965e970. http://dx.doi.org/10.1115/1.3408725. Farris, T.N., Szolwinski, M.P., Harish, G., 2000. Fretting in aerospace structures and materials. In: Hoeppner, D.W., Chandrasekaran, V., Elliott, C.B. (Eds.), Fretting Fatigue: Current Technology and Practice, ASTM STP, vol. 1367. ASTM, West Conshohocken, PA, pp. 523e537. Hills, D.A., Dini, D., January 2004. What level of friction guarantees adhesion in a complete contact? J. Strain Anal. Eng. Des. 39 (5), 549e551. http://dx.doi.org/ 10.1243/0309324041896425. Hills, D.A., Thaitirarot, A., Barber, J.R., Dini, D., 2012. Correlation of fretting fatigue experimental results using an asymptotic approach. Int. J. Fatigue 43, 62e75. http://dx.doi.org/10.1016/j.ijfatigue.2012.02.006.
J€ ager, J., 1998. A new principle in contact mechanics. J. Tribol. 120 (4), 677e684. http://dx.doi.org/10.1115/1.2833765. Jang, Y.H., Barber, J.R., 2011. Frictional energy dissipation in materials containing cracks. J. Mech. Phys. Solids 59 (3), 583e594. http://dx.doi.org/10.1016/ j.jmps.2010.12.010. Klarbring, A., Ciavarella, M., Barber, J.R., December 2007. Shakedown in elastic contact problems with coulomb friction. Int. J. Solids Struct. 44 (25e26), 8355e8365. http://dx.doi.org/10.1016/j.ijsolstr.2007.06.013. Law, S.S., Wu, Z.M., Chan, S.L., 2006. Analytical model of a slotted bolted connection element and its behaviour under dynamic load. J. Sound. Vib. 292 (3), 777e787. http://dx.doi.org/10.1016/j.jsv.2005.09.028. Limmer, L., Nowell, D., Hills, D.A., 2001. A combined testing and modelling approach to the prediction of the fretting fatigue performance of splined shafts. P. I. Mech. Eng. G J. Aer. 215 (2), 105e112. http://dx.doi.org/10.1243/0954410011531808. Melan, E., 1936. Theorie statisch unbestimmter systeme aus ideal-plastischem baustoff. In: Sitzungsber. Akad. Wiss. Wien, Abt. 2A, 145, pp. 195e218. Mindlin, R.D., Deresiewicz, H., 1953. Elastic spheres in contact under varying oblique forces. J. Appl. Mech. 21, 327e344. Nowell, D., Hills, D.A., 1987. Mechanics of fretting fatigue tests. Int. J. Mech. Sci. 29 (5), 355e365. http://dx.doi.org/10.1016/0020-7403(87)90117-2. Powell, M.J.D., 1978. A fast algorithm for nonlinearly constrained optimization calculations. In: Numerical Analysis. Springer, pp. 144e157. http://dx.doi.org/ 10.1007/BFb0067703. Rajasekaran, R., Nowell, D., 2006. Fretting fatigue in dovetail blade roots: experiment and analysis. Tribol. Int. 39, 1277e1285. http://dx.doi.org/10.1016/ j.triboint.2006.02.044. Szolwinski, M.P., Farris, T.N., 1996. Mechanics of fretting fatigue crack formation. Wear 198 (1), 93e107. http://dx.doi.org/10.1016/0043-1648(96)06937-2. Thaitirarot, A., Flicek, R.C., Hills, D.A., Barber, J.R., 2014. The use of static reduction in the finite element solution of two-dimensional frictional contact problems. P. I. Mech. Eng. C J. Mech. 228 (9), 1474e1487. http://dx.doi.org/10.1177/ 0954406213509086. Xi, N.S., Zhong, P.D., Huang, H.Q., Yan, H., Tao, C.H., 2000. Failure investigation of blade and disk in first stage compressor. Eng. Fail. Anal. 7 (6), 385e392. http:// dx.doi.org/10.1016/S1350-6307(99)00045-X.