Dynamic contact interactions of fractal surfaces

Applied Surface Science 392 (2016) 872–882

Contents lists available at ScienceDirect

Applied Surface Science journal homepage: www.elsevier.com/locate/apsusc

Dynamic contact interactions of fractal surfaces Tamonash Jana, Anirban Mitra, Prasanta Sahoo ∗ Department of Mechanical Engineering, Jadavpur University, Kolkata-700032, India

a r t i c l e

i n f o

Article history: Received 9 June 2016 Received in revised form 29 August 2016 Accepted 8 September 2016 Keywords: Fractal surface Contact Nonlinear vibration Tangent modulus Yield strength

a b s t r a c t Roughness parameters and material properties have significant influence on the static and dynamic properties of a rough surface. In the present paper, fractal surface is generated using the modified two-variable Weierstrass-Mandelbrot function in MATLAB and the same is imported to ANSYS to construct the finite element model of the rough surface. The force-deflection relationship between the deformable rough fractal surface and a contacting rigid flat is studied by finite element analysis. For the dynamic analysis, the contacting system is represented by a single degree of freedom spring mass-damper-system. The static force-normal displacement relationship obtained from FE analysis is used to determine the dynamic characteristics of the rough surface for free, as well as for forced damped vibration using numerical methods. The influence of fractal surface parameters and the material properties on the dynamics of the rough surface is also analyzed. The system exhibits softening property for linear elastic surface and the softening nature increases with rougher topography. The softening nature of the system increases with increase in tangent modulus value. Above a certain value of yield strength the nature of the frequency response curve is observed to change its nature from softening to hardening. © 2016 Elsevier B.V. All rights reserved.

1. Introduction A variety of mechanical elements, such as gear, cam-follower mechanism, rolling element bearings have dynamic contact at their contact surfaces. It is also well known that all surfaces are inherently associated with roughness. So, it can be inferred that vibration (or dynamics) at the contacting rough surface has major influence on the fatigue and wear performance of such components. Hence the study of dynamic contact interactions at rough surfaces has always been an area of research interest. The modeling of the rough surface during the dynamic analysis can be done using fractal geometry which was first introduced by Mandelbrot [1–4]. Since a rough surface is proved to be a random non-stationary process [5] due to its scale dependent nature, it can be characterized by fractal geometry. The applicability of fractal geometry as an efficient model of rough surface was presented by several researchers [6–8]. The equation for generating fractal surface, i.e. the Weierstrass–Mandelbrot equation was modified by Berry and Lewis [9] and Ausloos and Berman [10]. Few general distribution functions involving fractal parameters and the corresponding contact model were presented [11,12]. It was also found that the fractal roughness parameters are dependent upon the mean separation of

∗ Corresponding author. E-mail addresses: [email protected], [email protected] (P. Sahoo). http://dx.doi.org/10.1016/j.apsusc.2016.09.025 0169-4332/© 2016 Elsevier B.V. All rights reserved.

two contacting surfaces [13,14]. Analytical approach on the real contact area and deformation of elastic and elastic-plastic contact of fractal surfaces was carried out by Yan and Komvopoulos [15]. Kogut and Jackson [16] conducted a comparison between contact mechanics results obtained with statistical and fractal approaches tocharacterize surface topography. Several researchers have undertaken the study of normal contact stiffness and contact area of fractal surfaces in elastic and elastic-plastic regime [17–22]. In the last two decades, Finite Element Analysis has evolved as an efficient tool for contact analysis. Hyun et al. [23] and Pei et al. [24] adopted a 3D finite element analysis for elastic and elasto-plastic contact between rough surfaces with a range of self-affine fractal scaling behaviour. Influences of material and surface parameters on elastic and elastic-plastic contact interactions was analyzed through finite element analysis by several researchers [25–28]. In most of research works on fractal contact analysis available in the literature, the contact is assumed to be equivalently modeled rough to flat [8,15,19,26] or rough to spherical contact [25]. But analysis with rough to rough contact is found to be comparatively less in number. In recent years some researches have taken up contact analysis for rough to rough contact interfaces as the domain of their work [29–31]. It has been seen that when it comes to dynamic contact interactions between rough surfaces, Hertzian theory has been employed most extensively for modeling. First significant theoretical work on contact dynamics was by Nayak [32], who developed some of

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

Nomenclature A0 c D E E F Fs Fˆ f fmax G g k L Ls M m n P P∗ t u u(max) z zs  ı ı∗    ω ωs

Nominal contact area Damping coefficient Fractal dimension Elastic modulus of deformable surface Composite elastic modulus of equivalent rough surface Amplitude of harmonic excitation Restoring force of spring Normalized amplitude of harmonic excitation Frequency index Maximum frequency index Fractal roughness Gravitational acceleration Nonlinear stiffness of spring/deformable surface Sample length of fractal surface Cut off length of fractal surface Number of superimposed ridges to construct the fractal surface Mass of the block Exponent value in the power law force n displacement relationship equation P ∗ = k ı∗ Applied normal load on rough surface Normalized applied normal load, P ∗ = A PE 0 Time Normalized vertical displacement Maximum normalized vertical displacement (above which contact loss occurs) Vertical displacement of block Static vertical displacement of block due to its own weight Frequency density Displacement of contacting rough surface Normalized displacement of contacting rough surface, ı∗ = ı/L Poisson’s ratio Damping ratio Normalized time,  = tωs Harmonic forcing frequency Natural frequency of system at the static equilibrium position

the theoretical groundwork necessary for detailed physical explanations of experimentally observed phenomena in vibratory point contact by modeling a single-degree-of-freedom dynamic system and obtained analytical solutions by using a single term harmonic balancing method (HBM) [33]. Multi term HBM was also used by Ma et al. [34] while studying sphere-plane (Hertzian) contact model analytically as well as experimentally. Besides HBM, Method of Multiple Scales (MMS) [33] is used extensively by researchers for solving nonlinear differential equation of motion. Recently, Xiao et al. [35] used exact method, MMS and HBM to determine the natural frequency of undamped free vibration of a mass interacting with Hertzian contact stiffness. Hess and Soom [36,37] studied nonlinear vibrations at Hertzian contact as well as at the contact region formed between rough surfaces excited by the dynamic component of an externally applied normal load by MMS. Similar Hertzian contact problem for sub-harmonic and super-harmonic resonance of order two was investigated by Perret-Liaudet [38,39] using MMS. The condition for contact loss was taken into consideration for the first time in this work. Dynamic analysis of Hertzian contact using analytical and numerical methods was presented by number of researchers including Sabot et al. [40] and Perret-Liaudet

873

and Sabot [41]. Experimental analysis of sphere-plane contact with sinusoidal, random excitation and with subharmonic and superharmonic excitation was also carried out by some researchers [42–45]. However, very few research works is found on dynamic analysis of multi-asperity contact. Xiao et al. [46] recently modeled an elastic fractal surface contacting with a rigid flat surface, analyzed its forcedisplacement relation and studied its effect on the free and forced vibration responses of the system. Tian and Xie [47] investigated the dynamic contact stiffness at the interface between a vibrating rigid sphere and a semi-infinite transversely isotropic viscoelastic solid with an oscillating force superimposed onto a static compressive force. In the present paper the effect of variation of material properties of the contacting fractal surfaces on free and forced vibration responses of the system is analyzed. A three-dimensional rough fractal surface is constructed using a modified two-variable Weierstrass-Mandelbrot function [15] and the force-deflection relationship of the rough surface contacting with a rigid flat surface is determined using the finite element analysis. The current study considers nominally smooth contact which are assumed to be frictionless, as the effect of friction arising out of asperity interactions and atomic friction are neglected. The equation of motion of the system is represented by the Helmholtz–Duffing equation using a third-order Taylor series expansion which is dependent on the power value of the force-normal displacement relationship. For different post-elastic material properties such as tangent modulus (assuming bilinear model) and yield strength values, the force displacement relationship is determined. The natural frequency of the undamped free vibration of the system is determined by numerical quadrature method and for the forced damped system, harmonic response amplitude is calculated numerically using Runge-Kutta method. 2. Fractal surface modeling A realistic rough surface can be modeled by a 3D fractal surface topography [3]. In the present paper the fractal surface is generated with the help of modified two-variable Weierstrass-Mandelbrot function expressed in the following form

 D−2  z(x, y) = L

G L

 ×

ln  M

fmax M 1/2  

=0



cos  − cos

2 

 f

x2 + y L

 (D−3)f

f =0

1/2 2 × cos



tan

y −1 x

  − + M



(1)

where, L is the sample length, and G is a fractal height scaling parameter independent of frequency within the scale range.D,M and  are fractal dimension, number of superimposed ridges to construct the surface and frequency density, respectively. For physical rough surfaces the value of fractal dimension is found to be between 2 and 3 [15]. Random phase angle is represented by . Frequency index is represented byf . Frequency index has a lower limit at zero for a truncated series of the height function and the upper limit is given by,



fmax = int



log L/Ls log 



(2)

Where, int [.] represents the maximum integer value of the number within the bracket and Ls is the cut off length. The values of sample length and cut-off length are set to be 1 × 10−6 m and 1.5 × 10−7 m respectively. The parametersM,  and fmax have values of 10, 1.5 and 5 respectively. Low values of Gand high values of Dsignify smoother and highly dense surface profile and vice-versa. Using the above

874

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

Table 1 CLA and RMS roughness values for different fractal dimension (D) and fractal roughness (G) values. G (m)

D

CLA Roughness (m)

RMS Roughness (m)

1.36 × 10−11 1.36 × 10−11 1.36 × 10−11 1.36 × 10−11 1.36 × 10−10 1.36 × 10−12 1.36 × 10−13

2.3 2.4 2.5 2.6 2.4 2.4 2.4

48.221 × 10−9 16.794 × 10−9 5.9186 × 10−9 2.059 × 10−9 42.088 × 10−9 6.6883 × 10−9 2.6470 × 10−9

55.799 × 10−9 19.394 × 10−9 6.8175 × 10−9 2.3747 × 10−9 48.692 × 10−9 7.7110 × 10−9 3.0573 × 10−9

mentioned parameter values centre line average (Ra ) and root mean square roughness (Rq ) values for a generated 50 × 50 grid fractal surface for different combinations of G and D values are presented in Table 1. The surface heights z (x, y)at different (x, y)points are evaluated in MATLAB with the value of  generated by the random number generator of MATLAB. The surface function given by Eq. (1) possesses a scale-invariant (fractal) behaviour [9] only within a finite range of length scales, outside of which, the surface topography can be represented by a deterministic function. It may be noted here that the smallest length corresponds to the instrument resolution and the upper length to the length of the profile. Since frequencies outside the range determined by the lower and upper wavelengths do not contribute to the observed profile, self-similarity is satisfied at all scales only approximately. The fractal roughness G is a height scaling parameter independent of frequency within the scale range where fractal power-law behaviour is observable. Physically, higher G values correspond to rougher (less dense) surface topographies. The magnitude of the fractal dimension D determines the contribution of high and low frequency components in the surface function z(x,y). Thus, high values of D indicate that high-frequency components are more dominant than low-frequency components in the surface topography profile. The physical significance of D is the extent of space occupied by the rough surface, i.e., larger D values correspond to denser profiles (smoother topography). The surface height function given by Eq. (1) is continuous, nondifferentiable, scale-invariant within the range determined by the upper and lower wavelengths used in the truncated series, and selfaffine asymptotically according to the analysis of Blackmore and Zhou [12]. The self-affinity implies that as the surface is repeatedly magnified, more and more surface features appear and the magnified image shows a close resemblance to that of the original surface obtained at a different scale. These properties make the function given by Eq. (1) suitable for constructing surfaces possessing topographies closely resembling the actual surfaces with the same fractal parameters D and G. 3. Finite element modeling The surface heights z (x, y)as per the supplied x and y values generated in MATLAB are imported to ANSYS 14.5 as keypoints. Joining the keypoints, lines are generated and subsequently from the lines, areas are constructed. Finally a cubic solid block of dimension 1 ␮m × 1 ␮m × 1 ␮m is generated, where the top surface is the generated rough surface profile. The rough surface on the top is set as the CONTACT surface. A rigid flat surface is constructed to touch the contact surface from top which is identified as the TARGET surface. The rough deformable solid body is meshed using 3D solid element SOLID187 which is a higher order 3D, 10-node element, having a quadratic displacement behaviour. The TARGET surface is meshed with TARGE170 element which comes in contact with the CONTA174 elements overlaying the SOLID187 elements. Fig. 1(a) shows the constructed fractal surface simulated for D = 2.3, G = 1.36 × 10−11 m in MATLAB with 50 × 50 grid and Fig. 1(b) shows

Fig. 1. Weierstrass-Mandelbrot fractal surface for D = 2.3, G = 1.36 × 10−11 m (a) simulated in MATLAB with 50 × 50 grid and (b) meshed model of the surface in ANSYS.

the meshed model of fractal surface having the same G and D values.The mesh convergence study of the FE model is given in the Section 5.2. A pilot node is chosen at the centre of the rigid target surface where external load/displacement is to be applied. All the nodes at z = 0 are set fixed and the rigid plane is allowed to move in the z direction only. The contact with the rigid plane is realized using surface-to-surface contact elements that use the augmented Lagrangian method [48]. In the current method a downward displacement is imparted on the rigid surface incrementally with 100 sub-steps to come in contact and deform the deformable rough surface. The normal reaction force  at the contact  elements is recorded. The non-dimensional load P ∗ = P/A0 E  ) vs. non-dimensional nor-



mal displacement ı∗ = ı/L relationship in the form

 n

P ∗ = k ı∗



curves are fitted in a power law

(3)

k is a  constant parameter with positive value. E/ where  2

E = E/ 1 −  is the equivalent composite Young’s modulus, P is the applied load and A0 is the nominal contact area.

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

875

The equation of motion can be expressed in normalized form using the following non-dimensional variables. c u = zzs , ˝ = ωωs ,  = ωs t and  = 2mω s Natural frequency at the static equilibrium position,ωs =

nkzsn−1 m is

utilized to normalize the external forcing frequency. u,˝, and  are normalized displacement, frequency, time and damping ratio respectively. After normalizing, the dimensionless equation motion becomes u + 2u +

 1

1 (n + 1)n − 1 = F cos(t) n n

(6)

where, the non-dimensional excitation force is represented as Fˆ = F/mg and ’ denotes differentiation w.r.t. non-dimensional time. The expression (u + 1)n is expanded into a third order Taylor series expansion and hence the following equation is obtained. u + 2u + u + ˛2 u2 + ˛3 u3 =

1

F cos(t) n

(7)

(n−1)(n−2)

Fig. 2. Model of the (a) deformable body coming in contact with rigid flat surface and (b) the SDOF dynamic model representation.

4. Dynamic behaviour of the rough surface

4.2. Undamped free vibration

4.1. Dynamic model The present dynamic system consists of a deformable body of mass m with a rough surface at the top and a rigid flat surface capable of vertical motion. While in motion, the rigid flat comes in contact with the rough surface causing deflection at its multiple asperities (Fig. 2(a)). However, the rough surface is considered to be deflected with a bulk stiffness k for the entire surface. This dynamic system is equivalently modeled as a single degree of freedom spring-mass-damper system (Fig. 2(b)) and it is considered to be valid for low frequency behaviour only [6]. The mass of the block is equal to the mass of the deformable body whereas the spring stiffness is taken as the bulk stiffness of the surface. It should be pointed out that the spring is non-linear in nature. On the other hand, damping model is considered to be linear and viscous, with linear damping co-efficient c. The spring undergoes a deflection of zs due to the weight of the block and the static equilibrium position of the system is denoted by ’o’ as shown in Fig. 2(b). The system, under the influence of a harmonic external excitation, executes oscillation about this equilibrium position and z indicates the normal displacement of the system corresponding to the equilibrium. The force displacement relation for the nonlinear spring is expressed as Fs = k(z + zs )n where Fs is the restoring force and k is the nonlinear stiffness of the spring. At any time instant  t  the equation of motion of the system can be written as m¨z + c z˙ + Fs − mg = F cos (ωt)

The undamped natural vibration characteristics of the rough surface-rigid flat contact can also be analyzed by setting the damping factor and the forcing amplitude to zero. Substituting  = 0 and Fˆ = 0 in Eq. (7) the following expression is obtained. u +

 1 (n + 1)n − 1 = 0 n

n 2 1 u + (u + 1)n+1 − u = C0 2 n+1

n

m¨z + c z˙ + k(z + zs ) − mg = F cos (ωt)

(5)

(9)

where, C0 is the constant of integration which can be calculated from the initial conditions. For the present analysis following initial conditions are used (␶ = 0) = 0 u (␶ = 0) = 0 The value of C0 obtained from the above condition is C0 =

1 (u0 + 1)n+1 − u0 n+1

(10)

Putting the value of C0 in Eq. (9),



u =

2 n (n + 1)



(u0 + 1)n+1 − (u + 1)n+1 + (u − u0 ) (n + 1)



(11)

The natural time period ␶0 can be obtained by carrying out the integration



u(max)

0 = 2

du |u |

(12)

u(min)

1/n

aszs = mg/k . This equation of motion is only valid when the rough surface is in contact with the rigid flat surface, i.e. z ≥ −zs . Beyond this limit contact loss occurs. Substituting the expression for normal contact force incorporating the parameters k and n, the equation of motion becomes

(8)

Integrating Eq. (8) w.r.t non-dimensional time ()

(4)

where, F is the amplitude of the harmonic excitation, while ω is the forcing frequency. The static deflection of the spring is expressed



and ˛3 = where,˛2 = (n−1) 2 6 The governing differential equation describes linearly damped forced vibration of the contact interface under external harmonic excitation. As the equation contains both quadratic and cubic nonlinearity terms, it can be classified as a Duffing-Helmholtz type equation. This equation is solved numerically using Runge-Kutta method to obtain the frequency response curve for the vibrating rough surface.

 =2



u(max)

n (n + 1) 2

 u(min)

du (u0 + 1)n+1 − (u + 1)n+1 + (u − u0 ) (n + 1) (13)

876

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

Fig. 3. Comparison of the present method with the results of Buzio et al. [49] and Chatterjee and Sahoo [28].

The normalized natural frequency is obtained as. ˝0 = 2 /0 As mentioned earlier, contact loss between the rough surface and the rigid flat occurs whenu < −1. Hence, the minimum value of the initial displacement, u(min) is −1. Setting the value of u/ in Eq. (11) equal to zero, the maximum value of the initial displacement before contact loss is obtained as, u(max) = (n + 1)1/n − 1

(14)

The value of the integration of Eq. (13) is computed numerically using adaptive Simpson quadrature method and the normalized natural frequency is calculated. 5. Results and discussion 5.1. Validation problem Validation of the finite element modeling between a deformable rough surface and a rigid flat is obtained by comparison with experimental results provided by Buzio et al. [49] who has studied the load deflection behaviour between atomic force microscope (AFM) probe and nanostructured carbon films. The material property was taken as, equivalent Young’s modulus, (E/ ) = 0.88 GPa, Poisson’s ratio () = 0.3 and yield strength ( ys ) = 16.07 MPa. The surface properties obtained from the AFM were fractal dimension (D) = 2.3, fractal roughness (G) = 1.55 × 10−10 m and nominal contact area (A0 ) = 4.8 × 10−10 m. With the same material and surface properties, SOLID185, TARGE170 and CONTA174 elements are used to construct the solid, target and contact elements respectively during the FE analysis using ANSYS. The result (Fig. 3) is in favourable agreement with the experimental result [49] as well as the result furnished by Chatterjee and Sahoo [28].

Fig. 4. Mesh convergence plot at non-dimensional displacement values of (a) 0.02 and (b) 0.03.

elements and the material is taken as linearly elastic with values of elastic modulus and the Poisson’s ratio of 200 GPa and 0.3 respectively. The mesh convergence plots for two different nondimensional displacement values are given in Fig. 4(a) and (b). 5.3. Force-deflection characteristics of the rough surface

5.2. Mesh convergence study The mesh configuration for the finite element analysis is determined by increasing the number of divisions of the unmeshed lines i.e. increasing the mesh density in steps. The mesh density is increased until the difference of the value of contact force of two successive steps is less than 0.1%. As mentioned earlier, the FE model was constructed by SOLID187, CONTA174 and TARGE170

Force-deflection relation of the deformable rough surface in contact with the rigid flat is analyzed varying the surface properties. The fractal dimension D is varied from 2.3 to 2.6 and fractal roughness parameter G from 1.36 × 10−10 m to 1.36 × 10−13 m. The material is taken as infinitely linear elastic i.e. no yielding occurs up to infinite stress. The elastic modulus and the Poisson’s ratio of the material is taken as 200 GPa and 0.3 respectively. The variation of

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

877

Fig. 6. Plots of non-dimensional force versus non-dimensional displacement of rough surface for varying (a) tangent modulus (b) yield strength. Fig. 5. Plots of non-dimensional force versus non-dimensional displacement of rough surface for varying (a) fractal dimension (D) and (b) fractal roughness (G).

force-deflection behaviour in non-dimensional plane corresponding to varying D and G values are shown in Fig. 5(a) and (b) respectively. The k and n values from Eq. (2) for different G and D values are shown in Table 2. The n value is a measure of nonlinearity. It can be seen that the nonlinearity increases for rougher surface topography i.e. higher G value and lower D value and vice versa. Keeping the surface properties fixed at D = 2.4 and G = 1.36 × 10−11 m, the material is considered to yield and behave according to a bilinear model after yielding. Yield strength is considered to be fixed at 250 MPa and the value of the tangent

Table 2 k and n values for different fractal dimension (D) and fractal roughness (G) values. G (m)

D

k

n

1.36 × 10−11 1.36 × 10−11 1.36 × 10−11 1.36 × 10−11 1.36 × 10−10 1.36 × 10−12 1.36 × 10−13

2.3 2.4 2.5 2.6 2.4 2.4 2.4

2.5751 2.7927 1.9254 1.5275 5.4703 1.9939 1.5909

1.74723 1.5372 1.2719 1.1421 2.1092 1.2934 1.1638

878

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

Fig. 7. Comparison with undamped free vibration results of Xiao et al. [46].

modulus is varied as 100 GPa, 60 GPa and 10 GPa. The elastic modulus and Poisson’s ratio is set fixed at the previously mentioned value. Similarly finite element analysis is applied to obtain the force-displacement relationship to different tangent modulus values. The results are presented in Fig. 6(a). It can be seen that for higher tangent modulus value the surface exhibits higher nonlinearity. The elastic modulus, Poisson’s ratio and the tangent modulus is set fixed at 200 GPa, 0.3 and 10 GPa respectively and yield strength is varied as 250 MPa, 560.8 MPa, 911.5 MPa, 1265.3 MPa and 1619 MPa. These yield strength values cover a range of steels used in engineering applications [27,50]. The result of the FE analysis is showed in Fig. 6(b).For different tangent modulus and yield strength values, the coefficient and power values (k and n) in the Eq. (3) are given in Table 3. The power value (n) in Eq. (3) is found to decrease with the increase of yield stress value. An important phenomenon is observed for the yield strength values 560.8 MPa, 911.5 MPa, 1265.3 MPa, and 1619 MPa where the n value is found to be less than 1. 5.4. Undamped free vibration results For the validation of the dynamic analysis carried out in the present work, the non-dimensional natural frequency versus nondimensional initial displacement plots for fractal dimension (D) values 2.3 and 2.7 at undamped-free condition are compared with results of Xiao et al. [46] in Fig. 7. From the figure it is evident that the two sets of results are similar in trend. From Eq. (14) it can be seen that with increase of the roughness of the surface i.e. with surTable 3 kand n values for different tangent modulus (TM) and yield strength (YS) values. YS(MPa)

TM (GPa)

k

n

N/A 250 250 250 560.8 911.5 1265.3 1619

Linear Elastic 100 60 10 10 10 10 10

2.7927 1.2666 0.71652 0.06957 0.05301 0.05112 0.05416 0.05944

1.5372 1.4989 1.453 1.1057 0.8816 0.7709 0.7174 0.6914

Fig. 8. Plots of non-dimensional natural frequency vs. non-dimensional initial displacement of rough surfaces for varying (a) fractal dimension (D) (b) fractal roughness (G).

faces having higher value of n, the maximum displacement value for which the contact loss occurs, decreases. For example, for the surface with G = 1.36 × 10−11 m and D = 2.3, u(max) = 0.7832 and for G = 1.36 × 10−11 m and D = 2.6, u(max) = 0.9484 . The values of the non-dimensional natural frequency for non-dimensional initial displacement ranging from the minimum to the maximum values for the variation of fractal dimension (D) and fractal roughness parameter (G) values are shown in Fig. 8(a) and (b) respectively. In this case the material is considered to be linearly elastic upto infinite stress. The non-dimensional natural frequency reaches its maximum value of ˝0 = 1at the static equilibrium position ofu = 0. As the initial displacement increases or decreases from the static equilibrium position, the natural frequency decreases.

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

879

Fig. 10. Comparison between frequency response curves for surfaces having D values 2.3 and 2.6.

Fig. 9. Plots of non-dimensional natural frequency vs. non-dimensional initial displacement of rough surfaces for varying (a) tangent modulus (b) yield strength. Fig. 11. Comparison between frequency response curves for surfaces having G = 1.36 × 10−10 m and G = 1.36 × 10−13 m.

Results for variation of tangent modulus of the material of deformable surface are put forward in Fig. 9(a) and the results for the variation of yield strength are shown in Fig. 9(b). The corresponding n values used for the analysis are given in Table 3. A significant phenomenon is observed for the yield strength values 560.8 MPa, 911.5 MPa, and 1619 MPa, where the n values are found to be less than 1. For these values the nature of the initial displacement-natural frequency curves are vertically reversed w.r.t. the curve for yield strength 250 MPa. For these yield stress values, the non-dimensional natural frequency has its minimum value of ˝0 = 1at the static equilibrium position. For any positive or negative value within the limits the non-dimensional natural frequency is greater than 1, i.e. natural frequency at that position

is greater than the natural frequency at the static equilibrium position. 5.5. Forced damped vibration results The force damped dynamic characteristics of the system is represented by the normalized excitation frequency versus harmonic response amplitudeplots. The comparison between the forced damped characteristics for the surfaces having D = 2.3 and 2.6 is presented in Fig. 10. The value of G in both cases is taken as 1.36 × 10−11 m. Similarly the comparison between the frequency response curves for the surfaces with G = 1.36 × 10−10 m and 1.36 × 10−13 m is shown in Fig. 11. For both the cases the value

880

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

Fig. 12. Frequency response curves for rough surface material with yield stress 250 MPa and varying tangent modulus 10 GPa and 100 GPa.

of D is taken as 2.4. It can be observed that for linearly elastic surface, the system exhibits softening nature while vibrating. Rougher is the surface topography, higher is the softening nature of the system.The frequency response curves for the deformable surfaces having yield strength 250 MPa and varying tangent modulus values as 10 GPa, and 100 GPa are presented in Fig. 12. It can be observed that for higher tangent modulus, nonlinearity increases and the vibrating system tends to be more softening in nature. The jump up and jump down frequencies are indicated for each case. Now, keeping the tangent modulus fixed at 10 GPa, the yield strength is varied. For the yield strength values of 560.8 MPa, 911.5 MPa, 1265.3 MPa and 1619 MPa the frequency response curves are showed in Fig. 13 along with the comparison between the curves for the minimum and maximum yield strengths. Significantly hardening nature of the vibrating rough surface is observed for the yield stress values mentioned above while softening behaviour is observed for yield strength of 250 MPa. For higher yield strength value, the hardening nature of the system increases with higher nonlinearity. The present results indicate that, in order to predict the dynamic contact interactions of surfaces, adequate surface characterization in terms of the fractal nature must be made. The scale-independent fractal parameters D and G to be used in the analysis can be experimentally determined by some standard techniques such as structure function method, power spectrum method etc. Though the range of D and G values used in the present analysis refers to the same typically observed on real rough surfaces, the present study does not properly represent the real contact situations. This is due to the fact that the present analysis does not include the effects of adhesion, interfacial friction and tangential loading of the solids. The present results and trends are thus applicable when the magnitudes of tangential traction and adhesion at asperity microcontacts are insignificant. However, future studies will attempt to evaluate these effects. The present approach predicts the dynamic contact behaviour of fractal surfaces that agrees qualitatively well with observations. But, fractal surface structures span over a large range of length scales. Thus meshing surface structures in such a way that allows computational efficiency and realistic asperity morphologies is problematic. This poses the major disadvantage of using finite element analysis in fractal contact mechanics.

Fig. 13. Frequency response curves for rough surface material with tangent modulus 10 GPa and yield strength (a) 560.8 MPa, (b) 911.5 MPa, (c) 1265.3 MPa and (d) 1619 MPa and (e) comparison between 250 MPa and 1619 MPa.

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

881

to hardening in nature. Jump up and jump down frequencies are also shown for each case.

References

Fig. 13. (Continued)

6. Conclusion In the current paper, fractal surface is generated from the modified two-variable Weierstrass-Mandelbrot function in MATLAB and the points are imported to ANSYS as keypoints to construct the finite element model of the rough fractal surface. The force-deflection relationship between the deformable rough fractal surface and a rigid flat is determined by FE contact analysis. The contacting system is represented by a single degree of freedom spring-mass-damper system and the static force-displacement relationship is used to determine the dynamic characteristics of the rough surface for free as well as for forced damped vibration using numerical methods. The effect of the variation of fractal surface parameters and the material properties (tangent modulus and yield strength) are also analyzed. For linearly elastic surface, the system is found to be softening in nature and the softening nature increases with the increase of roughness of the surface. Similarly, the softening nature increases with increase in tangent modulus for elastic-plastic surface. A significant phenomenon is found while varying the yield strength. Above a certain value of yield strength the nature of the frequency response curve changes from softening

[1] B.B. Mandelbrot, How long is the coast of Britain, Science 156 (1967) 636–638. [2] B.B. Mandelbrot, Stochastic models for the earth’s relief, the shape and the fractal dimension of coastlines, and the number-area rule for islands, Proc. Natl. Acad. Sci. U. S. A. 72 (1975) 3825–3828. [3] B.B. Mandelbrot, The Fractal Geometry of Nature, W. H Freeman, New York, 1982. [4] B.B. Mandelbrot, Fractal geometry: what is it and what does it do? Proc. R. Soc. Lond A 423 (1989) 3–16. [5] R.S. Sayles, T.R. Thomas, Surface topography as a non-stationary random process, Nature (London) 271 (1978) 431–434. [6] F.F. Ling, Fractals, engineering surfaces and tribology, Wear 136 (1990) 141–156. [7] A. Majumdar, C.L. Tien, Fractal characterization and simulation of rough surfaces, Wear 136 (1990) 313–327. [8] G. Palasantzas, J.T. DeHosson, Self-affine roughness effects on the contact area between elastic bodies, J. Appl. Phys. 93 (2003) 898–902. [9] M.V. Berry, Z.V. Lewis, On the Weierstrass-Mandelbrot fractal function, Proc. R. Soc. Lond A 370 (1980) 459–484. [10] M. Ausloos, D.H. Berman, A multivariate Weierstrass-Mandelbrot function, Proc. R. Soc. Lond A 400 (1985) 331–350. [11] D. Blackmore, J.G. Zhou, A general fractal distribution function for rough surface profiles, SIAM J. Appl. Math. 56 (1996) 1694–1719. [12] D. Blackmore, J.G. Zhou, Fractal analysis of height distributions of anisotropic rough surfaces, Fractals 6 (1998) 43–58. [13] J.C. Chung, J.F. Lin, Variation in fractal properties and non-Gaussian distributions of microcontact between elastic-plastic rough surfaces with mean surface separation, J. Appl. Mech. 73 (2006) 143–152. [14] J.L. Liou, J.F. Lin, A new method developed for fractal dimension and topothesy varying with the mean separation of two contact surfaces, J. Tribol. 128 (2006) 515–524. [15] W. Yan, K. Komvopoulos, Contact analysis of elastic-plastic fractal surfaces, J. Appl. Phys. 84 (1998) 3617–3624. [16] L. Kogut, R.L. Jackson, A comparison of contact modeling utilizing statistical and fractal approaches, J. Tribol. 128 (2006) 213–217. [17] S. Jiang, Y. Zheng, H. Zhu, A contact stiffness model of machined plane joint based on fractal theory, J. Tribol. 132 (2010) 011401. [18] R. Pohrt, V.L. Popov, Investigation of the dry normal contact between fractal rough surfaces using the reduction method, comparison to 3D simulations, Phys. Mesomech 15 (2012) 31–35. [19] R. Pohrt, V.L. Popov, Normal contact stiffness of elastic solids with fractal roug hsurfaces, Phys. Rev. Lett. 108 (2012) 104301. [20] R. Pohrt, V.L. Popov, A.E. Filippov, Normal contact stiffness of elastic solids with fractal rough surfaces for one-and three-dimensional systems, Phys. Rev. E 86 (2012) 026710. [21] X. Miao, X. Huang, A complete contact model of a fractal rough surface, Wear 309 (2014) 146–151. [22] R. Buczkowski, M. Kleiber, G. Starzynski, Normal contact stiffness of fractal rough surfaces, Arch. Mech. 66 (2014) 411–428. [23] S. Hyun, L. Pei, J.F. Molinari, M.O. Robbins, Finite-element analysis of contact between elastic self-affine surfaces, Phys. Rev. E 70 (2004) 026117. [24] L. Pei, S. Hyun, J.F. Molinari, M.O. Robbins, Finite element modeling of elasto-plastic contact between rough surfaces, J. Mech. Phys. Solids 53 (2005) 2385–2409. [25] K. Komvopoulos, N. Ye, Three-dimensional contact analysis of elastic-plastic layered media with fractal surface topographies, J. Tribol. 123 (2001) 632–640. [26] P. Sahoo, N. Ghosh, Finite element contact analysis of fractal surfaces, J. Phys. 40 (2007) 4245. [27] B. Chatterjee, P. Sahoo, Finite element based parametric study of elastic-plastic contact of fractal surfaces, Proc. Eng. 64 (2013) 552–561. [28] B. Chatterjee, P. Sahoo, Element based contact analysis of fractal Surfaces–effect of varying elastic modulus, Proc. Eng. 90 (2014) 116–122. [29] D.A. Hanaor, Y. Gan, I. Einav, Contact mechanics of fractal surfaces by spline assisted discretisation, Int. J. Solids Struct. 59 (2015) 121–131. [30] D.A. Hanaor, Y. Gan, I. Einav, Static friction at fractal interfaces, Tribol. Int. 93 (2016) 229–238. [31] F. Alonso-Marroquin, P. Huang, D.A. Hanaor, E.A. Flores-Johnson, G. Proust, Y. Gan, L. Shen, Static friction between rigid fractal surfaces, Phys. Rev. E 92 (2015) 032405. [32] P.R. Nayak, Contact vibrations, J. Sound Vib. 22 (1972) 297–322. [33] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, Wiley, NewYork, 1979. [34] Q.L. Ma, A. Kahraman, J. Perret-Liaudet, E. Rigaud, An investigation of steady-state dynamic response of a sphere-plane contact interface with contact loss, J. Appl. Mech. 74 (2007) 249–255. [35] H. Xiao, M.J. Brennan, Y. Shao, On the undamped free vibration of a mass interacting with a Hertzian contact stiffness, Mech. Res. Commun. 38 (2011) 560–564. [36] D.P. Hess, A. Soom, Normal vibrations and friction under harmonic loads: part I—Hertzian contacts, J. Tribol. 113 (1991) 80–86.

882

T. Jana et al. / Applied Surface Science 392 (2016) 872–882

[37] D.P. Hess, A. Soom, Normal vibrations and friction under harmonic loads: part II—rough planar contacts, J. Tribol. 113 (1991) 87–92. [38] J. Perret-Liaudet, Subharmonic resonance of order 2 on a sphere-plane contact, C.R. Acad. Sci. Paris Series II 8 (1997) 443–448. [39] J. Perret-Liaudet, Superharmonic resonance of order two on a sphere-plane contact, C.R. Acad. Sci. Paris Series II 12 (1998) 787–792. [40] J. Sabot, P., Krempf, C. Janolin, Non-linear vibrations of a sphere–plane contact excited by a normal load, J. Sound Vib. 214(1998) 359–375. [41] J. Perret-Liaudet, J. Sabot, Vibro-impacts induced by non-linear resonances in hertzian contacts, in: Dynamics of Vibro-impact Systems, Springer, Berlin Heidelberg, 2016, pp. 1999251–1999260. [42] E. Rigaud, J. Perret-Liaudet, Experiments and numerical results on non-linear vibrations of an impacting Hertzian contact. Part 1: harmonic excitation, J. Sound Vib. 265 (2003) 289–307. [43] J. Perret-Liaudet, E. Rigaud, Experiments and numerical results on non-linear vibrations of an impacting Hertzian contact. Part 2: random excitation, J. Sound Vib. 265 (2003) 309–327. [44] J. Perret-Liaudet, E. Rigaud, Response of an impacting Hertzian contact to an order-2 subharmonic excitation: theory and experiments, J. Sound Vib. 296 (2006) 319–333.

[45] J. Perret-Liaudet, E. Rigaud, Superharmonic resonance of order 2 for an impacting Hertzian contact oscillator: theory and experiments, J. Comput. Nonlinear Dyn. 2 (2007) 190–196. [46] H. Xiao, Y. Shao, M.J. Brennan, On the contact stiffness and nonlinear vibration of an elastic body with a rough surface in contact with a rigid flat surface, Eur. J. Mech. A Solid 49 (2015) 321–328. [47] J. Tian, Z. Xie, Dynamic contact stiffness of vibrating rigid sphere contacting semi-infinite transversely isotropic viscoelastic solid, Acta Mech. Solida Sin. 21 (2008) 580–588. [48] Z.H. Zhong, Finite Element Procedures for Contact-impact Problems, NewYork, Oxford University Press, 1993. [49] R. Buzio, C. Boragno, U. Valbusa, Contac tmechanics and friction of fracta lsurfaces probed by atomic force microscopy, Wear 254 (2003) 917–923. [50] R.L. Jackson, I. Green, A finite element study of elasto-plastic hemispherical contact against a rigid flat, J. Tribol. 127 (2005) 343–354.