Fractal Geometry-Based Hypergeometric Time Series Solution to the Hereditary Thermal Creep Model for the Contact of Rough Surfaces Using the Kelvin-Voigt Medium

Abstract

This paper aims at constructing a continuous hereditary creep model for the thermoviscoelastic contact of a rough punch and a smooth surface of a rigid half-space. The used model considers the rough surface as a function of the applied load and temperatures. The material of the rough punch surface is assumed to behave as Kelvin-Voigt viscoelastic material. Such a model uses elastic springs and viscous dashpots in parallel. The fractal-based punch surface is modelled using a deterministic Cantor structure. An asymptotic power law, deduced using approximate iterative relations, is used to express the punch surface creep which is a time-dependent inelastic deformation. The suggested law utilized the hypergeometric time series to relate the variables of creep as a function of remote forces, body temperatures, and time. The model is valid when the approach of punch surface and half space is in the order of the size of the surface roughness. The closed-form results are obtained for selected values of the system parameters; the fractal surface roughness and various material properties. The obtained results show good agreement with published experimental results, and the methodology can be further extended to other structures such as the Kelvin-Voigt medium within electronic circuits and systems.

1. Introduction

Surface topography plays a significant role in tribology, that is, in problems of friction, wear, lubrication, and contact [1]. Therefore, the problem of analysis of rough surfaces attracts the attention of engineers and applied mathematicians. Historically, the following engineering parameters, statistical in nature, were used for the characterization of surface roughness: (1) the root mean square of the heights, 𝜎, (2) the root mean square of slopes, 𝜎²_𝑚, (3) and the root mean square of curvatures, 𝜎²_𝑘. However, it was realized that the topography of engineered surfaces is too complex to be described completely by a few statistical parameters. Thus, it was found that roughness has a multiscale nature and requires sophisticated mathematical techniques for its description.

First attempts to model the distribution of heights of surface asperities utilize the classical random field theory which assumed that the functions of surface model are differentiable. In particular, this implies that limiting values for 𝜎²_𝑚 and 𝜎²_𝑘 should exist as the sample interval tends to “0” [1]. However, it turned out that such limiting behavior is in contradiction with the results of advanced investigations of surfaces. For example, the exponential behavior of the autocorrelation function implies that the engineering parameters should tend to infinity rather than to constant values when the sampling interval is infinitely reduced [2]. Furthermore, it was shown that the profiles of a large number of both natural and artificial surfaces have the following form of the spectral density function 𝐺(𝜔)∼1/𝜔^𝜐 where 𝜐≈2, and 𝜔 is the spatial frequency (see (1.13)). It follows from this that all wavelengths are equally represented in the profile and that there exists no characteristic scale; in other words, after arbitrary magnification roughness looks like before. Moreover, it was found that the values of engineering parameters depend on the measurement scale, that is, these parameters are scale dependent [2, 3].

The fractal approach was introduced as an attempt to give a scale-invariant characterization of surface topography. The idea of fractality of roughness was experimentally verified on real surfaces as well as when applied to mathematically simulated profiles [4]. Figure 1 shows a picture of popular fractals, that is, the middle-third Cantor set, the von Koch curve, graph of the Weierstrass-Mandelbrot function C in the range 0≤𝑥≤3 (𝑝=1.5 and 𝛾=0.5 where p and 𝛾 are two numerical parameters (see (1.17)) where the trend of the function is ∼𝑥², and trajectories of a fractional Brownian process for different Hurst index 𝐻 and fractal dimension 𝐷 [5–7]. The index 𝐻 is the key parameter of the fractal surface which describes the smoothness of the surface.

(a)

(b)

(c)

(d)

Figure 1 

Common fractals: (a) the middle-third Cantor set, (b) the von Koch curve, (c) the Weierstrass-Mandelbrot function $C$ in the range $0\le x\le 3$ ($p=1.5$ and $\gamma =0.5$), and (d) trajectories of a fractional Brownian process for different $H$ and $D$.

Evidently, roughness of the surface of a body has a great influence on stress fields that arise when two deformable bodies are pressed together. Analysis of the effect of roughness on the contact interaction of solids has attracted wide attention [8].

One of the most popular models for studying contact of rough bodies is the Greenwood and Williamson (GW) model based on the use of the Hertz theory [9], where it is important to mention that (GW) model is a nonscale-invariant [10]. Currently, the development of models of contact between nominally flat fractal rough surfaces presented for the Cantor profile is an active area of research [11]. Various contact problems utilizing Cantor profile were considered [12–17]. All these models consider the one-level Cantor profile. It is argued that such profile is simple for analytical analysis. However, it has a minor drawback: all asperities of the profile have one-level character, while all real roughness has a hierarchical structure [17].

It is accepted that fractal dimension is not a compressive geometric parameter that could characterize alone the behavior of contacting rough bodies [1]. Moreover, the employment of the fractal approach in the study of surfaces has several drawbacks. The proposed model can be both fractal and nonfractal depending on values of the structural parameters. Regardless of this, the model profile remains rough and possesses certain self-affine properties. The iterative regular construction of the profile allows us to analyze its prestructures, that is, prefractals, of arbitrary generation.

In this introduction, important and relevant definitions and methods that are attributed to fractal geometry with the application to the modeling of rough surfaces will be fully presented. Furthermore, the important differences between mathematical and utilized physical fractals will be explicitly highlighted.

1.1. Mathematical Definition of Fractals

Mandelbrot stated that a set in a metric space is called a fractal set if the Hausdorff-Besicovitch dimension of the set is greater than its topological dimension [18]. Let 𝑋 be a compact metric space and 𝑂 be the totality of open balls in 𝑋. The Hausdorff 𝑠-measure of a subset 𝑆⊂𝑋 which is defined for 𝑠≥0 as the following limit:𝑚_𝐻(𝑆,𝑠)=lim_𝜎→0inf_𝐺∈𝑂{_𝑉∈𝐺(dim𝑉)^𝑠∶𝑆⊆_𝑉∈𝐺𝑉,diam𝑉≤𝛿}. (1.1) Here G is finite or denumerable subset of 𝑂. It was proven that there exists a value 𝑠₀ such that 𝑚_𝐻(𝑆,𝑠)=∞,for𝑠<𝑠₀,0,for𝑠>𝑠₀. (1.2)

The Hausdorff dimension of the set 𝑆, denoted by dim_𝐻𝑆, is the number 𝑠₀ such that (1.2) holds. Unfortunately, the calculation of the Hausdorff dimension of mathematical objects often demands a lot of effort. Even to find some estimations of the dimension, it is necessary to overcome a number of rather complex mathematical difficulties [19]. This issue called for the use of other definitions of dimension which are useful in applied mathematics for the characterization of fractal objects. One such alternative is the box dimension [1]. The analytical calculation of the box dimension is usually easier since the corresponding definition of this dimension involves coverings by spheres of equal radii.

Let 𝐸 be the Euclidean dimension of the space in which a set 𝑆 is embedded. For 𝛿>0, let 𝑁(𝛿) be the smallest number of 𝐸-dimensional balls or cubes of diameter 𝑑 needed to cover the set 𝑆. The box counting dimension or box dimension, denoted by dim_𝐵𝑆, can be defined if the following limit exists: dim_𝐵𝑆=lim_𝛿→0+log𝑁(𝛿)−log𝛿. (1.3)

It can be proven that dim_𝐵𝑆 does not change if one takes 𝑁(𝛿) to be either the smallest number of 𝛿-cubes that cover 𝑆: the number of 𝛿-mesh cubes that intersect 𝑆; or the smallest number of sets of diameter at most 𝑑 that cover 𝑆; or the largest number of disjoint 𝛿-balls with centers in 𝑆. Unfortunately, the box dimension is not always equal to the Hausdorff dimension. For example, the set 𝑆={0,1,1/2,…,1/𝑛…} has unequal values for the Hausdorff and box dimensions for dim_𝐻𝑆≠dim_𝐵𝑆=1/2. However, it can be proven that dim_𝐻𝑆≤dim_𝐵𝑆.

As a simple alternative to the Hausdorff measure, we can introduce the 𝑠-measure 𝑚_𝑠 of a set as the following limit: 𝑚_𝑠(𝑆)=lim_𝛿→0+𝑁(𝛿)𝛿^𝑠 (1.4) and define the box dimension as the value 𝑠=𝐷 such that 𝑚_𝑠(𝑆) has a jump from 0 to ∞ similar to the behavior of 𝑚_𝐻(𝑆,𝑠) in (1.2).

On the other hand, the difficulties involved with calculating the Hausdorff dimension are the reason for the opinion that the Hausdorff dimension is not generally used in applications in the study of fractal and non fractal curves that are originated in other sciences such as in biology, engineering, physics, quantum physics and computing [20–27].

1.2. Physical Concept of Fractals

Evidently, it is impossible to carry out the scaling procedure for any real physical object down to infinitely small scales. Hence, the mathematical concept of the Hausdorff measure is applicable only to mathematical models of objects rather than to the objects themselves and, of course, the Hausdorff dimension cannot be obtained by experimental procedures. In this sense there are no actual fractal objects in nature. For physical objects, the box dimension cannot be calculated analytically but it is estimated by experimental or numerical calculations. However, various errors can arise during such numerical calculations. There is no canonical definition of physical fractals and there are numerous methods for the practical estimation of the fractal dimension of an object. The cluster fractal dimension is taken as the first example of a physical fractal dimension definition.

Let a whole cluster be imagined as consisting of elementary parts of the size 𝛿_∗ [1]. An object can be modeled as a fractal cluster with dimension 𝐷 when the model considers scales 𝑅 such that 𝛿_∗<𝑅<Δ_∗, where 𝛿_∗ and Δ_∗are the upper and lower cutoffs for the fractal representation. To get the value 𝐷 of the dimension, the considered region is discretized into cubes with side length 𝛿_∗. Then the smallest number of 𝐸-dimensional cubes needed to cover the cluster (𝑁(𝛿_∗)) is counted. One says that the cluster is fractal if the numbers 𝑁(𝛿_∗) satisfy the so-called number-radius relation for different sizes of the considered region of the cluster 𝑅 as follows: 𝑁(𝛿_∗𝑅)≈(𝛿_∗)^𝐷,𝛿_∗<𝑅<Δ_∗. (1.5)

The value of 𝐷 is estimated as the slope of linear growth of ln(𝑁(𝛿_∗)) plotted against ln(𝑅). The power 𝐷 is usually called the cluster dimension or mass dimension.

In literature, various methods were utilized to estimate the fractal dimension of a physical object. However, the notion of fractal dimension is not well-defined in that the relative value does depend on the approach used. Indeed, only for the mathematical box dimension of a fractal set 𝑆 it is proven that dim_𝐵𝑆 is the same when using various specific schemes of covering [28], while for physical fractals the estimations of the fractal dimension inevitably involve various techniques, distinct scale ranges, and various computation rules. Therefore, the obtained values can differ strongly and it is unlikely that they could be fruitfully compared for distinct objects. Thus, even in the case of physical objects of a similar nature, it would be wrong to consider the fractal dimension of these objects as their specific property without referring to the estimation technique involved.

1.3. Self-Similarity and Self-Affinity of Surfaces

Let us recall that a one-to-one mapping 𝑀 of a plane 𝜋 onto a plane 𝜋^′ is called a similarity mapping with coefficient 𝜆>0, or simply a similarity, when the following property holds: if {𝐴,𝐵} are any two points of 𝜋, and {𝐴^′,𝐵^′} are their images under 𝑀, then |𝐴^′𝐵^′|=𝜆|𝐴𝐵| [29]. It is known that any similarity transformation of a plane is a homogeneous (isotropic) dilation of coordinates {𝑥^′=𝜆𝑥,𝑧^′=𝜆𝑧} up to a rotation and translation. A set 𝑆 is called statistically self-similar if under homogeneous scaling with the coefficient 𝜆, where 1>𝜆>0, it is identical from the statistical point of view to the set 𝑆^′=𝜆𝑆.

In practice, it is impossible to verify that all statistical moments of the two distributions are identical. Frequently, a set 𝑆 is said to be self-similar if only a few moments do not change under scaling [30]. A one-to-one mapping 𝑀 of a plane 𝜋 onto a plane 𝜋^′ is called an affine mapping, if the images of any three collinear points are collinear in turn [29]. In general, an affine transformation of a plane may be given in any coordinate system as a nondegenerative linear transformation. In practical studies of rough surfaces, one often considers a particular affine mapping, with anisotropic scaling, that is given coordinate wise by 𝑥^′=𝜆𝑥 and 𝑧^′=𝜆^𝐻𝑧. Here z is a graph of a surface profile and 𝐻 is some scaling exponent.

One says that a fractal is self-affine if it is invariant from the statistical point of view under quasihomogeneous (anisotropic) scaling. It is possible to show that usually a quasihomogeneous transformation is a particular case of Lipschitz homeomorphism [1, 17]. The Hausdorff dimension of a set 𝑆 does not change under the action of the Lipschitz homeomorphism 𝐿 as follows: dim_𝐻𝑆=dim_𝐻𝐿(𝑆). (1.6)

The ideas of self-similarity and self-affinity are very popular in studying surface roughness because experimental investigations show that usually profiles of vertical sections of real surfaces are statistically similar to themselves under repeatedly magnifications; however, the profiles should be scaled differently in the direction of nominal surface plane and in the vertical direction. The self-affine fractals were used in a number of papers as a tool for description of rough surfaces [3, 31, 32]. Two standard examples of self-affine fractals are the trace of the fractional Brownian motion and the Weierstrass function. The former is a statistical fractal while the latter is a deterministic fractal.

1.4. Brownian Surfaces and Random Fractals

Fractional Brownian processes are widely used in creating computer-generated surfaces, in particular landscapes. For example, a profile can be constructed as a graph of 1𝐷(𝑓𝐵𝑚)𝑉_𝐻(𝑥) of index 𝐻, where x is taken as the time and z is the random variable of the single valued function 𝑉_𝐻(𝑥) with the following property: ⟨[𝑉_𝐻(𝑥+𝛿)−𝑉_𝐻(𝑥)]²⟩∼𝛿^2𝐻,0<𝐻<1, (1.7) where ⟨⟩ denotes averaging over the ensemble, and 𝐻 is the Hurst index. The scaling behavior of the different traces, 𝑉_𝐻(𝑥), is characterized by a particular 𝐻 which relate the typical change in Δ𝑧(𝑥), where 𝑧(𝑥)=𝑉_𝐻(𝑥), is the trace of the fBm, and the change in the spatial coordinate Δ𝑥 by the simple scaling law [30, 33, 34]: Δ𝑧(𝑥)∼Δ𝑥^𝐻. (1.8)

It is known that, with probability equal to “1”, the following holds [28]: dim_𝐻𝑉_𝐻(𝑥)=dim_𝐵𝑉_𝐻(𝑥)=2−𝐻. (1.9)

The autocorrelation function is one of the main tools for studying statistical models of rough surfaces. The autocorrelation function 𝑅(𝛿) of the profile is𝑅(𝛿)=lim_𝑇→∞12𝑇^𝑇_−𝑇[𝑧(𝑥+𝛿)−𝑧][𝑧(𝑥)−𝑧]𝑑𝑥=⟨[𝑧(𝑥+𝛿)−𝑧][𝑧(𝑥)−𝑧]⟩ (1.10) or 𝑅(𝛿)=lim_𝑇→∞12𝑇^𝑇_−𝑇𝑧(𝑥+𝛿)𝑧(𝑥)𝑑𝑥−(𝑧)², (1.11) where 𝑧 is the average value of the profile function 𝑧(𝑥).

Another tool for the characterization of surfaces is the spectral density function 𝐺(𝜔) which is the Fourier transform of 𝑅(𝛿): 2𝐺(𝜔)=𝜋^∞₀2𝑅(𝛿)cos𝜔𝛿𝑑𝛿,𝑅(𝛿)=𝜋^∞₀𝐺(𝜔)cos𝜔𝛿𝑑𝜔. (1.12)

In general, it is accepted in fractional Brownian motion that [14]:(i)if the autocorrelation function 𝑅(𝛿) of the profile z(x) satisfies 𝑅(0)−𝑅(𝛿)∼𝛿^2(2−𝑠), then it is reasonable to expect that the box dimension of the graph z(x) is equal to s, note that one can find 𝑅(0)−𝑅(𝛿)∼𝛿^2𝐻 for the fBm defined by (1.7).(ii)if the profile 𝑧(𝑥) has spectral density: 1𝐺(𝜔)∼𝜔^𝜐, (1.13) then it is reasonable to expect that the box dimension of the graph z(x) is equal to (5−𝜐)/2 [1]. The above conclusions are valid for mathematical models of the profile, for which the relation 2(5−𝑠)=𝜐−1𝑜𝑟𝜐=5−2𝑠 holds. The exponent 𝜐 varies typically between 0 and 2. Usually, it is assumed that the same conclusions concerning the box dimension are valid for physical fractals as well. It is shown that real surfaces approximately satisfy the property in (1.13) in wide range of scales [35]. The moments 𝑚_𝑛 of the spectral density 𝐺(𝜔) provide a useful description of the surface roughness: 𝑚_𝑛=^∞_𝜔0𝜔^𝑛𝐺(𝜔)𝑑𝜔, (1.14) where 𝜔₀=2𝜋/𝜆₀ is the wave number corresponding to the profile length 𝜆₀. It is possible to show that 𝑚₀ is the variance of heights (rms height) of the surface, 𝑚₁ is the variance of slopes (rms slope) and 𝑚₃ is the variance of curvatures (rms curvature) [36].

1.5. Weierstrass-Type Functions and Modeling of Rough Surfaces

A number of researchers have used the Weierstrass-type functions for fractal modeling of surface roughness [3, 31, 32] and fractal modeling applications such as in quantum computing [20, 21]. The real Weierstrass-type function can be defined as:𝑊(𝑥;𝑝)=^∞_𝑛=0𝑝^−𝛾𝑛ℎ(𝑝^𝑛𝑥),𝑝>1,0<𝛾<1, (1.15) where h is a bounded Höِlder function of order greater than 𝛽. The following complex generalization of the 𝑊(𝑥;𝑝)was considered:𝑊(𝑥;𝑝)=^∞_𝑛=−∞𝑝^−𝛾𝑛[(1−𝑒^{𝑖𝑝𝑛𝑡})𝑒^𝑖Φ𝑛],𝑝>1,0<𝛾<1, (1.16) where Φ_𝑛 are arbitrary phases [29].

The Weierstrass-type functions are continuous everywhere and differentiable nowhere. In addition, their graphs are curves whose fractal dimension exceeds one. Fractal properties of these functions including the Weierstrass-Mandelbrot (WM) function C and the Takagi-Hopson function T:𝐶(𝑥;𝑝)=^∞_𝑛=−∞𝑝^−𝛾𝑛(1−cos𝑝^𝑛𝑥),𝑝>1,0<𝛾<1,𝑇(𝑥;𝑝)=^∞_𝑛=−∞𝑝^−𝛾𝑛|𝑝^𝑛𝑥−[𝑝^𝑛1𝑥+2]|,𝑝>1,0<𝛾<1, (1.17)(1.18) have been studied in numerous papers [2, 14, 29, 31]. By direct calculations, one may obtain: |𝑊(𝑥+𝛿;𝑝)−𝑊(𝑥;𝑝)|∼𝛿^𝛾, (1.19) which is similar to the behavior of (1.7) of fractional Brownian motion. The box dimension of the Weierstrass function graphs is 𝐷=2−𝛾 and it is believed that their Hausdorff dimension is the same [28, 37]. Currently, the only known bounds for the Hausdorff dimensions are 𝐷−(𝑐/log𝑝)≤dim_𝐻graph𝐶≤𝐷, provided that p is large and constant c is large enough [19]. It is possible to calculate the spectral density of the WM function 𝑊(𝑥;𝑝) as follows:𝐺(𝜔)=^∞_𝑛=−∞𝛿(𝜔−𝑝^𝑛)𝑝^{2(2−𝐷)𝑛}, (1.20) where 𝛿 is the Dirac delta function. Some arguments for approximating this discrete spectral density by a continuous spectral density 𝐺(𝜔)∼1/𝜔^5−2𝐷, whose exponent (5-2D) is in agreement with (1.13) with respect to the box dimension were suggested [2]. The following truncated WM function 𝑊₁(𝑥;𝑝)=𝐴^{∞(𝐷−1)}_𝑛=𝑛1𝑝^{(𝐷−2)𝑛}cos2𝜋𝑝^𝑛𝑥 (1.21) is often used for fractal characterization of the surface topography [3, 31, 32]. Here 𝑛₁ is an integer number, which corresponds to the low cutoff frequency of the profile, and A is the so-called characteristic length scale of the profile. The number 𝑛₁ depends on the length 𝐿 of the sample and is given by 𝑝^𝑛1=1/𝐿 and the parameter A determines the position of the spectral density along the log𝐺 axis. It was stated that both parameters A and 𝐷 of the function 𝑊₁ are scale-invariant characteristics of the roughness. However, the extensive experimental studies of this fractal characterization model showed that the values of parameters A and 𝐷 are not unique and depend on instruments or resolution of a given instrument [1].

Evidently, the function 𝐶(𝑥;𝑝) is not homogeneous. Nevertheless, it exhibits the property 𝐶(𝑝^𝑘𝑥;𝑝)=𝑝^𝑘𝛾𝐶(𝑥;𝑝), with 𝑘∈𝑍 where Z is the set of all integers which looks similar to the definition of a homogeneous function ℎ_𝑑 of degree 𝑑, that is, ℎ_𝑑(𝜆𝑥)=𝜆^𝑑ℎ_𝑑(𝑥) for 𝜆>0.

Thus, the graph of the function 𝐶(𝑥;𝑝) near any point 𝑥₀ is repeated in scaling form near all points 𝑝^𝑘𝑥₀,𝑘∈𝑍. This scaling (self-affine) property was often attributed to fractal features of the graph. However, this discrete scaling property is the main property of the so-called parametric-homogeneous (PH) functions introduced [1, 17] which strictly satisfy the equation 𝑏_𝑑(𝑝^𝑘𝑥;𝑝)=𝑝^𝑘𝑑𝑏_𝑑(𝑥;𝑝), with 𝑘∈𝑍 where 𝑑 is degree of homogeneity. As examples of 1-dimensional fractal PH-curves we can consider the graphs of functions 𝑏₁ and 𝑏₂ with degrees 𝑑=1 and 𝑑=2, respectively: 𝑏₀(𝑥;𝑝)=𝑥^−𝛾𝑏𝐶(𝑥;𝑝),₁(𝑥;𝑝)=𝑥𝑏₀𝑏(𝑥;𝑝),₂(𝑥;𝑝)=𝑥²𝑏₀(𝑥;𝑝). (1.22)

Because of (1.6), these functions have the same Hausdorff dimension as the WM function 𝐶(𝑥;𝑝) whose box-dimension is 𝐷.

Another consequence is that the WM function 𝐶(𝑥;𝑝), with 𝐶(𝑥;𝑝)∼𝑥^2-𝐷 can be used only as an example of fractal profile and it cannot be considered as the general fractal functional model for simulations of the rough surface profiles. The assumption that the WM function represents the general fractal properties of rough profiles can lead to wrong conclusions concerning surface roughness parameters and their distributions.

The solution to the problem of mechanical contact between elastically deforming solids was obtained by Hertz [8]. Subsequently, several approaches were used to analyze the contact interaction between the soft layer and the indenting object surface [38–42]. These methods are based upon Radok's technique of replacing the elastic constants in the elastic solution by the corresponding integral or differential operators, which appear in the stress-strain relations for linear viscoelastic materials. Furthermore, these studies assumed that the surfaces of contacting solids are smooth, excluding from consideration all real solids, which have a certain degree of roughness and waviness regardless of how fine their finish is [43].

Various models for the approach of the fractal punches were considered [11–16]. In the previously cited works, different constitutive relations were considered: (1) linear elastic material [11], (2) rigid-perfectly plastic material [13], (3) elastic-perfectly plastic material [12], (4) linear viscoelastic creep model via Maxwell medium [14], (5) linear viscoelastic creep model via standard linear solid (SLS) material [15], and (6) linear thermoviscoelastic relaxation model via Maxwell medium [16].

The objective of this work is to introduce an alternative approach, using fractal geometry, to study the deformation of a viscoelastic surface as a function of the force applied and the bulk temperature. In this model friction force effect is assumed to be negligible. The development of the fractal model of the rough surface is carried out using fractional Brownian motion in conjunction with Cantor set. The Radok's technique [44] is then used to derive the thermoviscoelastic model from the corresponding elastic model. The main contribution of this work is a mathematical model for the time-dependent-creep of a rough surface (cf. (6.7)). This model relates the creep to time, temperature, external applied load, fractal dimension of the rough surface, and various material properties.

Section 2 presents the fractal model, where the Cantor structure is built and its fractal dimension is presented. Section 3 presents the discrete and continuous elastic model. In Sections 4 and 5, the effect of temperature on the viscoelastic behavior is presented and the Arrhenius's relation is introduced. In Section 6, the elastic viscoelastic correspondence is presented which consists of replacing the elastic constant in the elastic solution by the corresponding integral or differential operators from the viscoelastic stress-strain relations. Also, in Section 6 a new continuous model for the creep contact of a thermovisco-elastic punch is presented. In Section 7, the results obtained from the new model is presentd and compared with an experimental results obtained from literature. In Section 8, conclusions and future work are presented.

2. Fractal Model

The surface profile of the punch, in contact with a rigid half-space, will be constructed on the basis of Cantor set [11]. The contacting surface is constructed by joining the segments obtained at successive stages of the construction of a Cantor set to one another, Figure 2, where 𝐿₀ correspond to the profile nominal length, and ℎ₀ is equal to the twice rms height of the roughness.

Figure 2 

The fractal middle-third Cantor structure, where $E_0$ is the initiator step, $E_1\text{,}$$E_2$, and $E_3$ are the other generated step of cantor structure, $L$'s are the lengths of the $E$'s steps, h's are the heights of $E$'s steps, and F is the applied load.

At each stage of profile construction, the middle section of each initial segment is discarded so that the total length of the remaining segments is 1/𝑎 times the length of the initial segment, where 𝑎>1. The depth of the recesses (measured from the last step) at the (𝑖+1)th construction step of the fractal surface is 1/𝑏 times less than the depth of the 𝑖th step, where 𝑏>1. From this it can easily be shown that the horizontal length and recess depth of the (𝑖+1)th step are, respectively𝐿_𝑖+1=𝑎⁻¹𝐿_𝑖=𝑎^−(𝑖+1)𝐿₀,ℎ_𝑖+1=𝑏⁻¹ℎ_𝑖=𝑏^−(𝑖+1)ℎ₀, (2.1)(2.2) where it is assumed that the surface is smooth in a direction perpendicular to the plane of the page. This restriction is not expected to have a significant effect since it is possible to construct a fractal Cantor surface perpendicular to the plane of the page [11].

At the 𝑖th generation, the Cantor structure contains 𝑁=2^𝑖 segments, each of length 𝛿_𝑖=(2𝑎)^−𝑖𝐿₀ [11]. The profile of the surface in Figure 2 can be considered as a certain graph of a step function.

It can be seen that, during an iterative step in constructing the surface, scaling in the horizontal direction is Δ𝜒_𝑖+1=(2𝑎)⁻¹Δ𝜒_𝑖, while in the vertical direction, the corresponding fluctuations Δ𝑧_𝑖 at the 𝑖th generation can be defined by considering the probability of obtaining the value, 𝑧_𝑖=𝑏^−𝑖ℎ₀.

The fluctuation Δ𝑧_𝑖 at the 𝑖th generation can be obtained by assuming the Δ𝑧_𝑖 scales as the expected value 𝑧_𝑖𝑃(𝑧_𝑖) in which Δ𝑧_𝑖∝𝑧_𝑖𝑃(𝑧_𝑖) [8], where 𝑃(𝑧_𝑖) is the probability of obtaining the value 𝑧_𝑖, that is, 𝑃(𝑧_𝑖)=(𝐿_𝑖−𝐿_𝑖+1)/𝐿₀, and it is found that 𝑃(𝑧_𝑖)=𝑎^−𝑖(1−1/𝑎).

Thus, the expected value of the fluctuation at the (𝑖+1)th generation is related to the expected value of the fluctuation at the 𝑖th generation through 𝑧_𝑖+1𝑃(𝑧_𝑖+1)=(𝑎𝑏)^−𝑖𝑧_𝑖𝑃(𝑧_𝑖). Hence Δ𝑧_𝑖+1=(𝑎𝑏)^−𝑖Δ𝑧_𝑖, and thus (Δ𝑧_𝑖+1/Δ𝑧_𝑖)=(Δ𝜒_𝑖+1/Δ𝜒_𝑖)^2-𝐷, from which the self-affine fractal dimension for the contour of the Cantor structure is derived as: 𝐷=1+ln2−ln2𝑎ln𝑏ln2𝑎=1+𝐷_𝑐−ln𝑏ln2𝑎for(1<𝐷<2), (2.3) where 𝐷_𝑐 is the fractal dimension of the Cantor set (0<𝐷_𝑐<1). Equation (2.3) will be used in the next section in the development of the approach-force model.

3. The Continuous Elastic Model

Qualitatively, two size scales are manifested in the contact problem [8]:(1)the bulk scale, for which the elastic compression would be calculated by the Hertz theory and its limitations, (2)the roughness scale, where the asperities act like a compliant layer on the surface, and so all the deformations are limited in a surface layer which represents all the asperities; 𝑏ℎ₀ in Figure 2, and their deformation is assumed to be linear elastic [45].

In this paper, the approach of the punch of Cantor structure surface and length 𝐿₀ will be considered. It is to be noted that the obtained relation may be applied for all problems with surfaces having the same fractal dimension. The contact between two rough surfaces can be modeled as the contact of an effective surface with a rigid flat surface [10]. Hence, a solution for the deformation of an equivalent surface generated using the Cantor structure can be modified for the problem at hand. The bodies treated in this work will be assumed to be isotropic and homogeneous, and obey linear force-displacement laws. The yield strength 𝜎_𝑦, the modulus of elasticity 𝐸, and coefficient of thermal expansion𝛼, are all assumed to be independent of temperature.

Furthermore, it is assumed, with reference to Figure 2, that there exists a series of one-dimensional elastic bars, distributed in a way such that the distance from the initiator step 𝐸₀ to the generated step 𝐸₃ is indicated by ℎ₀, from 𝐸₁ to 𝐸₃ is indicated by ℎ₁, from 𝐸₂ to 𝐸₃ is indicated by ℎ₂, and so forth, [14–16]. By letting 𝐹₃ be the force required to compress 𝐸₃ until 𝐸₂, 𝐹₂ be the force require to compress 𝐸₃ and 𝐸₂ until 𝐸₁, and 𝐹₁ be the force required to compress 𝐸₃, 𝐸₂, and 𝐸₁ until 𝐸₀, and assuming unit depth, one obtains the following discrete force-displacement relations: 𝐹₃=ℎ₂𝑘₃,𝐹₂=(ℎ₁−ℎ₂)𝑘₂+ℎ₁𝑘₃,𝐹₁=(ℎ₀−ℎ₁)𝑘₁+(ℎ₀−ℎ₂)𝑘₂+ℎ₀𝑘₃, (3.1) where 𝑘_𝑖=𝐸𝐿_𝑖/𝑏ℎ₀ is the stiffness of the 𝑖th step, 𝐸 is the modulus of elasticity of the material used, 𝑏ℎ₀ could be understood from Figure 2, and ℎ_𝑖 and 𝐿_𝑖 can be calculated using (2.1) and (2.2), respectively.

It is to be noted that thermal forces and deflections may arise in heated body either because of a nonuniform temperature distribution, or external constraints, or as a combination of these causes. The problem is assumed to be a steady state one with no internal heat source.

Next, by letting Δ𝐹_𝑖+1=𝐹_𝑖−𝐹_𝑖+1, then from (3.1) one can conclude the general equation for any number of steps as follows: Δ𝐹_𝑖+1=𝐸𝐿₀𝑏(𝑏−1)𝑏^−𝑖𝑎^−𝑖. (3.2)

In order to find a recursive relation as in (3.2) for the approach u, one lets 𝐹_𝑖+1 be the limit force for protrusion of the (𝑖+1)th generation. It is assumed that when the limit load is reached, the punch approaches a distance Δ𝑢_𝑖+1, equals to the difference between the heights protrusion of 𝑖th and (𝑖+1)th generations. Consequently the second generation 𝐸₂ deflects a distance 𝑢₂=ℎ₁, 𝐸₁ deflects a distance 𝑢₁=ℎ₀, and 𝐸₀ deflects a distance 𝑢₀=𝑏ℎ₀, so Δ𝑢₁=𝑢₀−𝑢₁=ℎ₀(𝑏−1), and Δ𝑢₂=ℎ₀𝑏⁻¹(𝑏−1). Accordingly:Δ𝑢_𝑖+1=ℎ₀(𝑏−1)𝑏^−𝑖. (3.3)

The above-mentioned assumptions are sufficient to determine the dependence of the limit load F on the approach u. The effects of the remote load and the bulk temperature will be first studied separately and then superimposed. Using the fact that, when the limit load increases from 𝐹_𝑖+1 to 𝐹_𝑖, the punch is approached by an amount Δ𝑢_𝑖+1, and by utilizing (3.2) and (3.3), the remote load effect is given by the discrete force displacement relation: Δ𝐹_𝑖+1Δ𝑢_𝑖+1=𝐸𝐿₀𝑏ℎ₀𝑎^−𝑖. (3.4)

As 𝑖→∞, (3.4) yields the following asymptotic behavior for the strain 𝜀: 𝜀=(𝑏𝜒𝐸𝐹𝐿₀)^1/(𝜒+1), (3.5) where the strain 𝜀 could be defined as 𝑢/ℎ₀ and 𝜒=ln𝑎/ln𝑏 [14–16]. It is to be noted that, the effect of the applied external load should not exceed a limiting yield load 𝐹_𝑦, where 𝐹_𝑦=𝜎_𝑦𝐿₀.

Since the interest of this work is to consider the viscoelastic behavior, the principle of correspondence [44] will be used in the next section to obtain a viscoelastic model corresponding to the elastic model presented in (3.5).

4. Effect of Temperature

Temperature has a dramatic influence on rates of viscoelastic response, and in practical work it is often necessary to adjust a viscoelastic analysis for varying temperature. This strong dependence of temperature can also be useful in experimental characterization, for example, if a viscoelastic transition occurs too quickly at room temperature, for easy measurement, the experimenter can lower the temperature to slow things down and vice versa.

In some viscoelastic materials, the relation between time and temperature can be described by correspondingly simple models. Such materials are termed “thermorheologically simple” [46]. For such simple materials, the effect of lowering the temperature is simply to shift the viscoelastic response (plotted against log-time) to the right without change in shape. This is equivalent to increasing the relaxation time 𝜏, without changing the relaxation modulus.

A time-temperature shift factor (𝑎_𝑇) can be defined as the horizontal shift that must be applied to a response curve, measured at an arbitrary temperature T in order to move it to the curve measured at some reference temperature 𝑇_ref.

If the creep time obeys an Arrhenius relation, the shift factor can be shown to be [47]:log𝑎_𝑇=𝑄(12.303𝑅𝑇−1𝑇₀), (4.1) where 𝑄 is the activation energy (J/mol), 𝑅 is the gas constant (J/mol⋅K), and T is the temperature (K).

5. Arrhenius Relation

The creep properties of materials are usually described by reference to the dependence of the creep (𝜀) on the applied stress, time and temperature, which may be written as: 𝜀=𝑓(𝜎,𝑡,𝑇). (5.1)

One way to simplify this function is to make it to be separable into three functions of stress, time and temperature as follows [38]: 𝜀=𝑓₁(𝜎)𝑓₂(𝑡)𝑓₃(𝑇). (5.2)

Temperature has a significant effect on the creep of materials. In some steels, it is found that the temperature has a pronounced effect than the strain rate [48].

Thermal forces and deflections may arise in a heated body either because of a nonuniform temperature distribution, or external constraints. The problem is assumed to be a steady state one with no internal heat. Arrhenius relation is a simple, but remarkably accurate, formula for the temperature dependence, where according to Arrhenius law, the temperature dependence is given as [49]: 𝑓₃𝑄(𝑇)=𝐵exp(𝑅𝑇), (5.3) where B is a constant, Q is the activation energy and 𝑅 is the ideal gas constant. The functions 𝑓₁(𝜎) and 𝑓₂(𝑡) will be included in (6.4) in Section 6.

It is clear that at the reference temperature, 𝑇₀, the temperature function 𝑓₃(𝑇) is equal to unity, and the creep will be a function of the stress and time, as it was shown in (5.1). From (5.3) at the reference temperature 𝑇₀, it can be easily shown that the value of the constant B will be: 𝐵=exp(−𝑄𝑅𝑇₀). (5.4)

6. Elastic Viscoelastic Correspondence

The simplest approach to this problem consists of replacing the elastic constant in the elastic solution by the corresponding integral or differential operators from the viscoelastic stress-strain relations [44]. This approach can be applied to the contact problem provided that the loading program is such that the contact area is increasing throughout [38].

A Kelvin-Voigt linear model is employed to describe the viscoelastic behavior of the compliant layer. Such a model is an arrangement of spring and dashpot in parallel, as shown in Figure 3, in which 𝜂and 𝐸_𝑒 are the Newtonian viscosity and the elastic modulus, respectively. The time-dependent force-displacement relation could be written in the operator form using the linear differential time operator 𝜕_𝑡≡𝜕/𝜕𝑡 as shown in (6.1) [50]: 𝜎𝜀=𝐸_𝑒+𝜂𝜕_𝑡. (6.1)

Figure 3 

The linear Kelvin-Voigt model, where $\eta \mathrm{\hspace{0.17em}\hspace{0.17em}}$is the Newtonian viscosity, $E_e$ is the elastic modulus, F is the applied load and u indicates the points to be displaced.

It is clear that simple the constant of proportionality between stress and strain does no longer exist. The viscoelastic operator corresponding to the modulus 𝐸 in (3.5) could be written as [50]: 1𝐸→1𝐸_𝑒11+𝜏𝜕_𝑡, (6.2) where (𝜏≡𝜂/𝐸_𝑒) is a characteristics parameter with units of time called the retardation time.

Creep test is a widely used standard test, wherein a force 𝑃₀ is suddenly applied at time 𝑡=0 on the viscoelastic model and then maintained constant thereafter, while measuring the approach as a function of time. The applied force can be expressed as a function of time with the aid of the unit step function[𝑈(𝑡)]. Thus 𝐹=𝑃₀[𝑈(𝑡)]. (6.3)

By substituting (6.1)–(6.3) in (3.5) one obtains𝑓₁(𝜎)𝑓₂(𝑡)=(𝑢(𝑡)ℎ₀)₁=(𝑏𝜒𝜏𝐸_𝑒𝑃₀𝐿₀)^1/(𝜒+1)𝑡^1/(1+𝜒)Γ((2+𝜒)/(1+𝜒))¹𝐹₁[1;1+𝜒2+𝜒𝑡1+𝜒;−𝜏] (6.4)

Tables of Laplace transforms [51] were utilized to obtain (6.4), where Γ represent the gamma function, and ₁𝐹₁[𝑐;𝑑;𝑥] is the Kummer's confluent hypergeometric function which could be expressed as [52] ₁𝐹₁𝑐[𝑐;𝑑;𝑥]=1+𝑑𝑥+𝑐(𝑐+1)𝑥²+𝑑(𝑑+1)2!𝑐(𝑐+1)(𝑐+2)𝑥³𝑑(𝑑+1)(𝑑+2)3!+⋯ (6.5) or ₁𝐹₁[𝑐;𝑑;𝑥]=^∞_𝑛=0(𝑐)_𝑛(𝑑)_𝑛𝑥^𝑛𝑛!. (6.6)

By substituting the value of the constant B from (5.4) into (5.3), and substituting (3.5) and (5.3) into (5.2), the creep stain as a function of stress, time, and temperature, is obtained as follows: 𝜀=𝑓₁(𝜎)𝑓₂(𝑡)𝑓₃(𝑇)=𝑢(𝑡)ℎ₀=(𝑏𝜒𝜏𝐸_𝑒𝑃₀𝐿₀)^1/(𝜒+1)𝑡^1/(1+𝜒)Γ((2+𝜒)/(1+𝜒))¹𝐹₁[1;1+𝜒2+𝜒𝑡1+𝜒;−𝜏]𝑄×exp(𝑅(1𝑇−1𝑇₀)), (6.7) where 𝑢(𝑡)/ℎ₀ is the approach, 𝑃₀/𝐿₀ is the applied stress per unit depth, T is the bulk temperature, 𝑇₀ is the reference temperature, 𝜏≡𝜂/𝐸_𝑒 is the retardation time, 𝜂is the Newtonian viscosity, 𝐸_𝑒 is the elastic modulus, Q is the activation energy, and 𝑅 is the ideal gas constant. The effect of the fractal dimension 𝐷 appears through the constant 𝜒 which combines the two scaling parameters {𝑎,𝑏}, that is, 𝜒=ln𝑎/ln𝑏 (cf. (2.3)). The results obtained using this analytical model are presented and discussed in the next section.

7. Results and Discussion

A new continuous model for the creep contact of a thermovisco-elastic punch has been presented in (6.7). The model presents an approximate closed form solution for the approach 𝑢(𝑡)/ℎ₀ of the fractal surface as a function of the applied load 𝑃₀/𝐿₀. In order to use (6.7), values for the system parameters 𝐸_𝑒, 𝜂, a and b are needed. The constant parameters a and b, in this equation, characterize the Cantor structure of the rough surface and are related through the fractal dimension 𝐷 of the rough surface. For the Kelvin-Voigt model to be capable of describing the experimental results of various viscoelastic materials, the viscous coefficients, 𝜂(GPa⋅sec), and the elastic modulus 𝐸_𝑒 (GPa), should be selected properly. In this study the modulus 𝐸_𝑒 is selected to be 130 GPa and the retardation time, 𝜏≡𝜂/𝐸_𝑒, is selected to be of the order “1” [53]. The value of n is taken large enough for the result to converge to an accepted accuracy; it is assumed that 𝑛=550.

The Cantor structure, shown in Figure 2, is built from the middle-third Cantor set where the parameter 𝑎=1.5 is held fixed, giving a Cantor set dimension 𝐷_𝑐=0.63093 [18]. Two different values of the parameter 𝑏 (𝑏=1.155 and 𝑏=1.29) yield two different dimensions of the Cantor structure; 𝐷=1.5 and 𝐷=1.4, respectively, which are used to verify the proposed analytical model. It was pointed out [11] that only for 𝑏≤2 the profile of the surface of the contacting body is fractal.

In order to examine the validity of the presented model, results obtained using this model for selected values of the system parameters were compared with those obtained experimentally by the results in [12, 54]. These results displayed the approach-force relation between a flat rough surface and an ideally smooth and rigid counter surface. The specimens used in these experiments were made of carbon steel (0.45 percent carbon). Their surface roughness resulted from different finishing processes; face turning, grinding, and bead-blasting. The experiments were conducted using MATLAB for a wide range of the nominal load, up to 600 MPa. The error in the experimental measurements was determined to be approximately ∓0.5𝜇m for the approach, and ∓5 MPa for the load [54]. Furthermore, the fractal dimension, 𝐷, of a ground stainless steel surfaces is 𝐷=1.5 [12]. The value of ℎ₀, which corresponds to twice rms height for the ground surface, is taken as 6.6𝜇m [54], and the value 𝐷=1.5 [37] is used to calculate the value of the parameter b, using (2.3), for a fixed value of the parameter a.

Figure 4 shows a set of numerical creep data obtained by applying the model which was presented in (6.7). In Figure 4, the strain 𝑢/ℎ₀ is plotted versus the nondimensional time record 𝑡/𝜏 for Δ𝑇=0, with two different fractal dimensions 𝐷=1.5 and 𝐷=1.4, and two applied stresses 𝜎₁=300 MPa and 𝜎₂=500 MPa. As might be expected, higher strain rates occur for higher temperatures for a constant stress.

Figure 4 

Non-dimensional time-strain curves for Kelvin-Voigt model when $\mathrm{\Delta }T=0$, with $D=1.5$ and $D=1.4$, and two applied Stresses ${\sigma }_1=300$ MPa and ${\sigma }_2=500$ MPa.

Figure 5 shows the strain 𝑢/ℎ₀ versus the nondimensional time record 𝑡/𝜏. It shows, also, a set of numerical creep data obtained applying (6.7) for three temperature differences Δ𝑇=0, Δ𝑇=100 and Δ𝑇=200 with two different fractal dimensions 𝐷=1.5 and 𝐷=1.4, and constant applied stress 𝜎₁=300 MPa. The behavior is also similar to the known typical creep curves; it indicates that higher strain rates occur for higher temperatures.

Figure 5 

Non-dimensional time-strain curves for Kelvin-Voigt model when for three temperature differences $\mathrm{\Delta }T=0$, 100 and 200, with $D=1.5$ and $D=1.4$, and constant applied stress $\sigma =300$ MPa.

Figure 6 is similar to Figure 5 in all of its aspects except for the applied stress, where in this case 𝜎=500 MPa.

Figure 6 

Non-dimensional time-strain curves for Kelvin-Voigt model when for three temperature differences $\mathrm{\Delta }T=0$, 100 and 200, with dimensions $D=1.5$ and $D=1.4$, and constant applied stress $\sigma =500$ MPa.

Figures 4–6 show that higher strain rates result for higher fractal dimensions, which could be explained from the definition of the fractal dimension itself. Lower fractal dimensions means less roughness and consequently the bulk material dominates, while higher fractal dimensions means excessive roughness and consequently the asperities deformation dominate.

Figure 7 presents the isochronous stress-strain creep curves accompanied by experimental results available in the literature [54] for nondimensional time durations 𝑡/𝜏≈44 (which is held constant), fractal dimension 𝐷=1.4, and for three temperature differences Δ𝑇=0, Δ𝑇=100 and Δ𝑇=200.

Figure 7 

Experimental results and analytical isochronous Stress-Strain curves for constant $t/\tau \approx 44$ for three temperature differences $\mathrm{\Delta }T=0\text{,}$$\mathrm{\Delta }T=100$ and $\mathrm{\Delta }T=200$, with $D=1.4$.

Figure 8 also presents the isochronous stress-strain creep curves accompanied by experimental results available in the literature [54] which are conducted at the room temperature. The analytical results are shown for the nondimensional time durations 𝑡/𝜏≈42 (which is held constant), fractal dimension 𝐷=1.5, and for temperature differences Δ𝑇=0, Δ𝑇=100 and Δ𝑇=200.

Figure 8 

Experimental results and analytical isochronous Stress-Strain curves for constant $t/\tau \approx 42$ for three temperature differences $\mathrm{\Delta }T=0\text{,}$$\mathrm{\Delta }T=100$ and $\mathrm{\Delta }T=200$, with $D=1.5$.

Figures 7 and 8 show good agreement between the presented proposed model and the experimental data. For 𝐷=1.4 the nondimensional time duration, 𝑡/𝜏, required to get an agreement between the experimental results and the isochronous curves is about 44.2, while it is about 44 for the fractal dimension 𝐷=1.5.

It is clear that the relatively longer duration of agreement occurs for lower fractal dimensions which could be attributed to the same reason mentioned above, that is, lower fractal dimensions means less roughness and consequently the bulk material dominates, while higher fractal dimensions means excessive roughness and consequently the asperities deformation dominate. The mathematical model also shows instability when 𝑡/𝜏 exceeds 45.

8. Conclusions and Future Work

As well known, creep analysis is a nonlinear time-dependent phenomenon. The model which is modified in this work presents a solution to the thermal creep-contact of rough surfaces as a hypergeometric time series. Fractal geometry, via Cantor set, is utilized to model roughness of the creeping contact surfaces. The results obtained by this model turn out not to be too far from reality, since tests, at room temperatures, on the actual contact area of ground metal surfaces show that they contain sets of parallel ragged-edged scratches of different depths.

Since the construction of the Cantor structure is periodic in its nature, it undergoes the same construction procedure at each hierarchical level producing contact areas that are all of the same size. Therefore, the presented analytical model provides an approximate, not exact, simulation of the approach of the viscoelastic rough surfaces, where the presented model shows a fairly good agreement with the available experimental results. As linearity is an inherent assumption, it is not expected from this model to be able to describe exactly the real material behavior; roughness and deformation. A nonlinear thermoviscoelastic stress-strain relation is required for the reproduction of real material behavior. It is also clear that the specific character of the fractal model has little effect on the asymptotic behavior of the process, and the fractal dimension 𝐷 which provides a measure of the rate at which a surface is changing is of most importance. The solution obtained in this work provides further insight into the effect that surface structure has on the deformation process, and it also provides indications of the effect that different surface forming processes may have on the subsequent surface deformation. Furthermore, in the averaged sense, the Cantor structure model appears to provide fairly reasonable results.

For future work, it is intended to extend the methodology which is used in this paper for the application to the contact-surfaces within electronic and electrical devices and circuits such as resistors, capacitors and inductors that arise in electronic manufacturing systems. Furthermore, it is intended to further investigate the applications of fractals in the emerging quantum computing domain.