Recent Advancements in Fractal Geometric-Based Nonlinear Time Series Solutions to the Micro-Quasistatic Thermoviscoelastic Creep for Rough Surfaces in Contact

Abstract

To understand the tripological contact phenomena, both mathematical and experimental models are needed. In this work, fractal mathematical models are used to model the experimental results obtained from literature. Fractal geometry, using a deterministic Cantor structure, is used to model the surface topography, where recent advancements in thermoviscoelastic creep contact of rough surfaces are introduced. Various viscoelastic idealizations are used to model the surface materials, for example, Maxwell, Kelvin-Voigt, Standard Linear Solid and Jeffrey media. Such media are modelled as arrangements of elastic springs and viscous dashpots in parallel and/or in series. Asymptotic power laws, through hypergeometric series, were used to express the surface creep as a function of remote forces, body temperatures and time. The introduced models are valid only when the creep approach of the contact surfaces is in the order of the size of the surface roughness. The obtained results using such models, which admit closed-form solutions, are displayed graphically for selected values of the systems' parameters; the fractal surface roughness and various material properties. Results obtained showed good agreement with published experimental results, where the utilized methodology can be further extended to the utilization for the contact of surfaces within micro- and nano-electronic devices, circuits and systems.

1. Introduction

In this paper, Section 1 presents the practical implications and motivations for our work and the basic background materials on fractal geometry from mathematical and physical point of views. Section 2 presents the fractal contact model of the middle-third Cantor structure which is utilized to obtain the new results introduced in this paper. The continuous elastic model used in this paper is presented in Section 3. Section 4 presents the effect of temperature on viscoelastic materials and the Arrhenius relation. Section 5 presents the elastic-viscoelastic correspondence principle. Results and discussions for the various viscoelastic models are presented in Section 6. Conclusions and future works are presented in Section 7.

1.1. Motivations and Practical Implications

Contact problem is central to solid mechanics, as contact is the principal method of applying loads to a deformable body and the resulting stress concentration is often the most critical point in the body [1, 2]. Therefore, the determination of specific pressure and deformation in the contact area of two deformable bodies is of a major concern in many practical applications [3, 4]. One of the most important features of the contacting surface is its contact stiffness, that is, its ability to bear a load. Therefore, determination of surface deformations plays a vital role in the assembling of machines; in many instances the deformations are responsible for “sinking” effect in contact surfaces, occurring in bolted joints, also for slacking of nuts [1, 3]. Resilient deformations in contact joints plays extremely important role in machinery dynamics. It often affects considerably stiffness also vibration damping of complex mechanical systems. As far as dynamics of machines is concerned, much significance is given to stabilized characteristics of contacting joints that reflect their force displacement relationships [2, 4].

Real contacting surfaces are rough, leading to the concentration of contact in a cluster of microscopic actual contact areas. The level or roughness depends on the surface preparation but even the most thoroughly polished surfaces show irregularities. As a consequence, at incipient approach of two surfaces, contact will be imperfect and electrical, thermal and mechanical properties will be affected, and in general in a different way as if surfaces are assumed smooth. Also in many applications involving contact mechanics, the tribological performance of surfaces can be enhanced by covering them with thin layers of soft materials. The dissimilar elastic properties of those layers to those of the substrate are of particular interest owing to their increasing applications in industry [1, 3, 4]. Many of these applications include overlay soft lubricants to hard surfaces, sealing rubbers and viscoelastic load bearing members as, that is, in the motor vehicle industry.

Many mechanical joints such as bolt joints, press fits, springs, and dies work under creep and/or stress relaxation condition; therefore, in many cases the characteristics of creep and/or relaxation stability are the determining characteristics of the material, and the reliability of many structures incorporating such elements are closely linked to the creep and/or relaxation stability of the material used. Joint elements carry in service static stresses caused by initial tightening, and cyclic stresses generated by operational loads. Such stresses are responsible for the strains which cause microscopic changes in the dimensions of joint elements. These changes are responsible for the reduction in the initial tightening, and a possible fatigue failure. The danger of failure during creep and/or stress relaxation becomes real if the part has stress concentrations which sharply reduce the deformability of materials.

In electrical engineering, an electrical connector must maintain good contact force throughout the functional life of the product. In creep and/or stress relaxation, the failure mode is a loss of contact force, and an inability for a deflected contact to return to its original position. However, the normal force necessary to maintain electrical contact will no longer be present. Therefore, creep and/or stress relaxation can be interpreted as a constant increase in contact resistance across the contact interface, eventually leading to an open circuit [1].

The mechanical implications are manifold and involve issues of friction, wear and fatigue on one hand and elements like bearings and gears on the other. Also the automotive industry has started to specify creep and/or stress relaxation tests for critical sealing products in their vehicles. Accordingly the evaluation of the time-dependent stress and strain is important in order to make clear the fundamental mechanism of those phenomena which spoil the mechanical functions and shorten their lifecycle. For these reasons, contact of machined surfaces has been and still are subjected to analytical and experimental studies.

1.2. Fractals Background

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, 𝜎^2𝑚, and (3) the root mean square of curvatures, 𝜎^2𝑘. 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 assumed that the functions of surface model are differentiable. In particular, this implies that limiting values for 𝜎^2𝑚 and 𝜎^2𝑘 should exist as the sample interval tends to “0” [2]. 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 [3]. 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 (cf. (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 [3, 4].

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 [5]. Figure 1 shows a picture of several popular fractals, that is, the middle-third Cantor set, the von Koch curve, graph of the Weierstrass-Mandelbrot function, trajectories of a fractional Brownian process for different Hurst indices 𝐻 (i.e., the Hurst parameter which is the key parameter of the fractal surface that describes the smoothness of the surface) and for different fractal dimensions 𝐷, and the Sierpinski gasket (triangle) [6–10]. One may note that the representations of the fractal dimension 𝐷 and the Hurst index 𝐻 in [6–8] differ from those in [9, 10]; in [9, 10], the representations for 𝐷 and 𝐻 are based on the fractional Brownian motion (fBm), while in [6–8] 𝐷 and 𝐻 are separately and independently represented based on a fractal model which is called the generalized Cauchy (GC) process. The work in [6] provides a new bound of fractal time series, the work in [7] uses the GC process for the internet traffic modeling, and the work in [8] addresses the simulation of the GC process.

(a)

(b)

(c)

(d)

(e)

Figure 1 

Common fractals: (a) the middle-third Cantor set, (b) the von Koch curve, (c) the Weierstrass-Mandelbrot function 𝐶 in the range 0≤𝑥≤3 (𝑝=1.5 and 𝛾=0.5), where 𝑝 and 𝛾 are two numerical parameters (cf. (1.17)), where the trend of the function is ~𝑥2, (𝑑) trajectories of a fractional Brownian process for different 𝐻 and 𝐷, and (e) the Sierpinski gasket (triangle), where the four small diagrams show the point of departure of the construction, then its first three stages, while the large diagram shows an advanced stage.

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 [11].

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 [12], where it is important to mention that (GW) model is a nonscale-invariant [13]. Currently, the development of models of contact between nominally flat fractal rough surfaces presented for the Cantor profile is an active area of research [14]. Various contact problems utilizing Cantor profile were considered [15–24]. 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 [15].

It is accepted that fractal dimension is not a compressive geometric parameter that could characterize alone the behavior of contacting rough bodies [2]. 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.2.1. Mathematical Formulation 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 [9]. 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 𝐺 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 [25]. 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 [2]. 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 𝑆; or 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 no fractal curves that are originated in other sciences such as in biology, engineering, physics, quantum physics and computing [26–33].

1.2.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 𝛿_∗ [2]. 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: 𝑁𝛿_∗≈𝑅𝛿_∗^𝐷,for𝛿_∗<𝑅<Δ_∗.(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 [34], 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.2.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 |𝐴^𝐵^|=𝜆|𝐴𝐵| [35]. 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 [36]. 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 [35]. 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 𝑧 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 [2, 15]. 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 [4, 37, 38]. 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.2.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−𝐷 (fBm) 𝑉_𝐻(𝑥) of index 𝐻, where 𝑥 is taken as the time and 𝑧 is the random variable of the single valued function 𝑉_𝐻(𝑥) with the following property: 𝑉^𝐻(𝑥+𝛿)−𝑉_𝐻(𝑥)²∼𝛿^2𝐻,for0<𝐻<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 [36, 39, 40]: Δ𝑧(𝑥)∼Δ𝑥^𝐻.(1.8)

It is known that, with probability equal to “1”, the following holds [34]: 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 expressed as:𝑅(𝛿)=lim_𝑇→∞12𝑇^𝑇−𝑇𝑧(𝑥+𝛿)−𝑧𝑧(𝑥)−𝑧𝑑𝑥=𝑧(𝑥+𝛿)−𝑧𝑧(𝑥)−𝑧,(1.10) which also can be expressed as: 𝑅(𝛿)=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𝜋^∞0𝑅(𝛿)=2𝑅(𝛿)cos𝜔𝛿𝑑𝛿,𝜋^∞0𝐺(𝜔)cos𝜔𝛿𝑑𝜔.(1.12)

In general, the following is accepted in fractional Brownian motion [2] (i)If the auto-correlation function 𝑅(𝛿) of the profile 𝑧(𝑥) satisfies 𝑅(0)−𝑅(𝛿)∼𝛿^2(2−𝑠), then it is reasonable to expect that the box dimension of the graph 𝑧(𝑥) is equal to 𝑠, 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 𝑧(𝑥) is equal to (5−𝜐)/2 [2]. The above conclusions are valid for mathematical models of the profile, for which the relation 2(5−𝑠)=𝜐−1 or 𝜐=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 [41]. The moments (𝑚_𝑛) of the spectral density 𝐺(𝜔) provide a useful description of the surface roughness: 𝑚_𝑛=^∞𝜔₀𝜔^𝑛𝐺(𝜔)𝑑𝜔,(1.14) where (𝜔₀=2𝜋/𝜆₀) is the wave number corresponding to the profile length 𝜆₀. It is possible to show that 𝑚₀ is the variance of heights (r.m.s height) of the surface, 𝑚₁ is the variance of slopes (r.m.s slope), and 𝑚₃ is the variance of curvatures (r.m.s curvature) [42].

1.2.5. Weierstrass Functions and the Modeling of Rough Surfaces

A number of researchers have used the Weierstrass-type functions for fractal modeling of surface roughness [4, 37, 38] and Fractal modeling applications such as in quantum computing [26, 27]. The real Weierstrass-type function can be defined as: 𝑊(𝑥;𝑝)=^∞𝑛=0𝑝^−𝛾𝑛ℎ(𝑝^𝑛𝑥),for𝑝>1,0<𝛾<1,(1.15) where ℎ is a bounded Hölder function of order greater than 𝛽. The following complex generalization of the 𝑊(𝑥;𝑝) was considered:𝑊(𝑥;𝑝)=^{∞𝑛=−∞}𝑝^−𝛾𝑛1−𝑒^{𝑖𝑝𝑛𝑡}𝑒^𝑖Φ𝑛,for𝑝>1,0<𝛾<1,(1.16) where Φ_𝑛 are arbitrary phases [35].

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 𝐶 and the Takagi-Hopson function 𝑇:𝑇(𝐶𝑥;(𝑝𝑥);=𝑝)=^{∞𝑛=−∞}_{∞𝑛=−∞}𝑝^−𝛾𝑛(1−cos𝑝^𝑛𝑥),for𝑝>1,0<𝛾<1,(1.17)𝑝^−𝛾𝑛|||𝑝_𝑛𝑥−𝑝^𝑛𝑥+12|||,for𝑝>1,0<𝛾<1,(1.18) have been studied in numerous papers [3, 18, 35, 37]. 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 [34, 43]. Currently, the only known bounds for the Hausdorff dimensions are 𝐷−(𝑐/log𝑝)≤dim_𝐻graph𝐶≤𝐷, provided that 𝑝 is large and constant 𝑐 is large enough [25]. 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−2𝐷) is in agreement with (1.13) with respect to the box dimension were suggested [3]. The following truncated WM function: 𝑊₁(𝑥;𝑝)=𝐴^(𝐷−1)^∞𝑛=𝑛₁𝑝^{(𝐷−2)𝑛}cos2𝜋𝑝^𝑛𝑥(1.21) is often used for fractal characterization of the surface topography [4, 37, 38]. Here 𝑛₁ is an integer number, which corresponds to the low cut-off frequency of the profile, and 𝐴 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 𝐴 determines the position of the spectral density along the log𝐺 axis. It was stated that both parameters 𝐴 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 𝐴 and 𝐷 are not unique and depend on instruments or resolution of a given instrument [2]. Evidently, the function 𝐶(𝑥;𝑝) is not homogeneous. Nevertheless, it exhibits the property 𝐶(𝑝^𝑘𝑥;𝑝)=𝑝^𝑘𝛾𝐶(𝑥;𝑝), with 𝑘∈𝑍 where 𝑍 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 [2, 15] 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 [11]. Subsequently, several approaches were used to analyze the contact interaction between the soft layer and the indenting object surface [44–48]. These methods are based upon the 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 [49].

Various models for the approach of the fractal punches were considered and utilized [14, 16–21, 23, 24]. In the previously cited published works, different constitutive relations were considered: (1) linear elastic material [14], (2) rigid-perfectly plastic material [17], (3) elastic-perfectly plastic material [16], (4) linear viscoelastic creep model via Maxwell medium [18], (5) linear viscoelastic creep-contact model via Kelvin-Voigt medium [19], (6) linear viscoelastic creep model via Jeffreys’ type material [20], (7) linear viscoelastic creep model via standard linear solid (SLS) material [21], (8) linear thermo-viscoelastic relaxation model via Maxwell medium [22], (9) linear thermo viscoelastic creep-contact model via Kelvin-Voigt medium [23], and (10) linear thermo viscoelastic creep model via standard linear solid (SLS) material [24].

2. The Fractal Contact Model

In this work, the problem of contact between a nominal flat surface with a half-space is studied. The considered surface is constructed on the basis of a deterministic fractal; the Cantor structure. The contact problem for that surface is analyzed, assuming that the results hold for all problems with surfaces of the same fractal dimension. The Cantor structure is constructed by joining the segments obtained at successive stages of the construction of a Cantor set to one another, Figures 2 and 3, where 𝐿₀ corresponds to the nominal surface length, and ℎ₀ is equal to twice the r.m.s height of the roughness.

Figure 2 

The fractal middle-third Cantor structure 𝑠=2, where 𝐸0 is the initiator step, {𝐸1,𝐸2,𝐸3} are the other generated step of cantor structure, 𝐿's are the lengths of the 𝐸's steps, ℎ's are the heights of 𝐸's steps, and 𝐹 is the applied load.

Figure 3 

The fractal middle-third Cantor structure 𝑠=3, where 𝐸0 is the initiator step, {𝐸1,𝐸2} are the other generated step of cantor structure, 𝐿's are the lengths of the 𝐸's steps, ℎ's are the heights of 𝐸's steps, and 𝐹 is the applied load.

It is established that not all engineering surfaces are fractals in their nature from mathematical point of view, but when surfaces do possess fractal behavior, it is oftentimes valid only within a certain length range, and this is fractals from the physical point of view [50]. Therefore there will be some restrictions over the scaling factors that constitute Cantor structure, that is, 𝑐_𝑥 and 𝑐_𝑧, where this will be imposed in the sequel.

At each stage of 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 be easily shown that the horizontal length and recess depth of the (𝑖+1)th step are:𝐿ℎ_𝑖+1=𝑐^𝑥−1𝐿_𝑖=𝑐^{𝑥−(𝑖+1)}𝐿₀,(2.1)_𝑖+1=𝑐^𝑧−1ℎ_𝑖=𝑐^{𝑧−(𝑖+1)}ℎ₀,(2.2) respectively, in which the surface is assumed to be smooth in a direction perpendicular to the plane which is shown. This restriction is not very important because it is possible to construct a fractal Cantor surface perpendicular to the plane too [14]. The length of the constructed contour, after 𝑖 iterations, is found to be:𝐿_(𝑖)=𝐿₀+2ℎ⁰2𝑐^{𝑧−(𝑖+1)}−12𝑐^𝑧−1−1.(2.3)

It is obvious that only for 𝑐_𝑧≤2 the contour of the surface of the structure will be a fractal, because in this case the length 𝐿_(𝑖) tends to an infinite limit. It is known that most rough surfaces have a self-affine scaling structure which implies that length scales change by different amounts in different directions, [9, 10]. This is also evident with the case of the Cantor structure in Figure 2. At the 𝑖th generation, the Cantor structure contains 𝑁=𝑠^𝑖 segments, each of length [14]: 𝛿_𝑖=1𝑠𝑐_𝑥^𝑖𝐿₀,(2.4) where the parameter 𝑠 corresponds to the number of asperities on a repeating segment [16, 17].

In Figure 2 we have 𝑠=2, and in Figure 3 we use 𝑠=3. Therefore, by changing the value of 𝑠, an infinite number of different structures could be constructed.

The profile of the surface in Figures 2 and 3 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 produced as: Δ𝑥_𝑖+1=𝑠𝑐_𝑥⁻¹Δ𝑥_𝑖,(2.5) while in the vertical direction, the corresponding fluctuations Δ𝑧_𝑖 in the 𝑖th generation can be defined by considering the probability of obtaining the value of: 𝑧_𝑖=𝑐^𝑧−𝑖ℎ₀.(2.6)

The fluctuation Δ𝑧_𝑖 in the 𝑖th generation can be obtained by assuming that Δ𝑧_𝑖 scales as the expected value 𝑧_𝑖𝑃(𝑧_𝑖) [14] where: Δ𝑧_𝑖∝𝑧_𝑖𝑃𝑧_𝑖,(2.7) and 𝑃(𝑧_𝑖) is the probability of obtaining the value 𝑧_𝑖 which is given by: 𝑃𝑧_𝑖=𝐿^𝑖−𝐿_𝑖+1𝐿₀,(2.8) where it is also produced as: 𝑃𝑧_𝑖=𝑐_𝑥^−𝑖1−𝑐^𝑥−1.(2.9)

Thus, the expected value of the fluctuation in the (𝑖+1)th generation is related to the expected value of the fluctuation in the 𝑖th generation through: 𝑧_𝑖+1𝑃𝑧_𝑖+1=𝑐_𝑥𝑐_𝑧^−𝑖𝑧_𝑖𝑃𝑧_𝑖.(2.10) Hence, we have: Δ𝑧_𝑖+1=𝑐_𝑥𝑐_𝑧^−𝑖Δ𝑧_𝑖.(2.11)

Substituting (2.5) and (2.11) into (1.8), we get:Δ𝑧_𝑖+1Δ𝑧_𝑖=Δ𝑥^𝑖+1Δ𝑥_𝑖^2−𝐷,(2.12) from which the self-affine fractal dimension for the contour of the Cantor structure is:𝐷=1+ln𝑠ln𝑠𝑐_𝑥−lnln𝑠𝑐𝑐^𝑧_𝑥=1+𝐷^𝑐−ln𝑐^𝑧ln𝑠𝑐_𝑥,for1<𝐷<2,(2.13) where (0<𝐷_𝑐<1) is the fractal dimension of the middle-third Cantor set.

3. The Continuous Elastic Model

The main assumptions utilized for the theoretical assessment of the model are as follows.(1)The interactions between asperities are not taken into account. This assumption is confirmed by engineering practice; motivation for this assumption is that the height of the tallest asperities in contact is roughly 50–80 times smaller than the distance of its closest neighbor [51]. (2)The two rough surfaces in contact are modelled as an equivalent single rough surface in contact with a rigid smooth surface. The contact condition where the rigid surface is infinitely resistant to compression is referred to as the Signorini nonpenetration condition [52, 53]. (3)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 Figures 2 and 3, and their deformation is assumed to be linear elastic [54]. (4)The analysis in this study is based upon the Radok’s technique [55] 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. This approach can be applied to the contact problem provided that the loading program is such that the contact area is increasing throughout [44]. Thus the simplest approach to this problem is to find the compliance in cases where a corresponding solution for purely elastic material is known.(5)It is also assumed that, with reference to Figures 2 and 3, there exist a series of one-dimensional elastic bars distributed in a way 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, [18–24].

Let 𝐹₃ 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 so forth, or in general: Δ𝐹_𝑖+1=𝐸𝐿⁰𝑐_𝑧𝑐_𝑧−1𝑐^𝑧−𝑖𝑐^𝑥−𝑖.(3.1)

Also, let 𝐹_𝑖+1 be the limit force for protrusion of the (𝑖+1)th generation. We assume that when the limit load is reached, the surface approaches a distance Δ𝑢_𝑖+1 equal to the difference between the heights protrusion of 𝑖th and (𝑖+1)th generations. We will assume that 𝐸₂ is deflected by a distance 𝑢₂=ℎ₁, 𝐸₁ is deflected by a distance 𝑢₁=ℎ₀, 𝐸₀ is deflected a distance 𝑢₀=𝑐_𝑧ℎ₀, and so Δ𝑢₁=𝑢₀−𝑢₁=ℎ₀(𝑐_𝑧−1) and Δ𝑢₂=ℎ₀𝑐^𝑧−1(𝑐_𝑧−1). Thus: Δ𝑢_𝑖+1=ℎ₀𝑐_𝑧−1𝑐^𝑧−𝑖.(3.2)

The above-mentioned assumptions are sufficient to determine the dependence of the limit load 𝐹 on the approach 𝑢. We will use the fact that the surfaces were approached by an amount Δ𝑢_𝑖+1 when the limit remote increases from 𝐹_𝑖+1 to 𝐹_𝑖, then, for load effect only: Δ𝐹_𝑖+1Δ𝑢_𝑖+1=𝐸𝐿⁰𝑐_𝑧ℎ₀𝑐^𝑥−𝑖.(3.3)

In the limit when 𝑖→∞, the following asymptotic behavior is obtained as: 𝜀=𝑢ℎ₀=𝑐^𝑧𝜒𝐸𝐿𝐹₀^1/(𝜒+1),(3.4) where 𝜀=𝑢/ℎ₀ is the strain due to the remote applied force and 𝜒=ln𝑐_𝑥/ln𝑐_𝑧.

4. The Effect of Temperature on Viscoelastic Materials and the Arrhenius Relation

Temperature has a dramatic influence on the rates of viscoelastic response, and it is often necessary to adjust a viscoelastic analysis for varying temperature. This strong dependence of temperature can also be useful in experimental characterization. If for instance 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 “thermo-rheologically simple” [56]. For such simple materials, the effect of lowering the temperature is simply to shift the viscoelastic response (which is plotted against log time) to the right without a 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 𝑇 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 [57]: log𝑎_𝑇=𝑄2.303𝑅𝑇1−𝑇1_ref,(4.1) where 𝑄 is the activation energy (J/mol), 𝑅 is the gas constant (J/mol K), and 𝑇 is the temperature (K).

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

One way to simplify this function is to separate it into three functions of stress, time and temperature as follow [53, 58]: 𝜀=𝑓₁(𝜎)𝑓₂(𝑡)𝑓₃(𝑇).(4.3)

Temperature has a significant effect on the creep of materials. In some steels, it is found that the temperature has a more pronounced effect than the strain rate [59]. 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. The Arrhenius relation is a simple, but remarkably accurate, formula for the temperature dependence. According to the Arrhenius law [58], the temperature dependence is given as: 𝑓₃(𝑇)=𝐵exp𝑄𝑅𝑇,(4.4) where 𝐵 is a constant given by:𝐵=exp−𝑄𝑅𝑇₀.(4.5)

When combining (4.4) and (4.5), it could be shown that at the reference temperature 𝑇₀ the function 𝑓₃(𝑇) is equal to unity, and the creep is a function of the stress and time only.

5. The Elastic-Viscoelastic Correspondance Principle

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 [55]. This approach can be applied to the contact problem provided that the loading program is such that the contact area is increasing throughout [44].

Various linear viscoelastic models were employed to describe the viscoelastic behavior of the compliant layer (which is the assumption (1.3)). Such models are an arrangement of spring and dashpot in parallel and/or series as shown in Figure 4 [18–24] in which 𝜂's and 𝐸's are the Newtonian viscosity and the elastic modulus, respectively.

(a)

(b)

(c)

(d)

Figure 4 

Various linear viscoelastic models, where 𝜂is the Newtonian viscosity, 𝐸 is the elastic modulus, 𝑢 is the approach and 𝐹 is the applied external remote load.

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[𝑈(𝑡)], that is, 𝐹=𝑃₀[𝑈(𝑡)]. In this study the viscoelastic asperities forms a contact with a rigid half-space under a normal step creep load, the friction between the viscoelastic asperities and the rigid plane is assumed to be negligible.

Tables 1–4 show, for various viscoelastic Maxwell, Kelvin-Voigt, SLS and Jeffreys' models, the time-dependent constitutive equation, viscoelastic operator corresponds to the modulus of elasticity 𝐸, and the modified creep model, where (𝑢(𝑡)/ℎ₀) is the relative approach, (𝑃₀/𝐿₀) is the applied stress per unit depth, 𝑇 is the bulk temperature, 𝑇₀ is the reference temperature, 𝜂 is the Newtonian viscosity, 𝐸_𝑣 is the elastic modulus, 𝑄 is the activation energy, 𝑅 is the ideal gas constant, and (𝜏≡𝜂/𝐸_𝑣) is a characteristics parameter with units of time called the retardation time [18–24]. The effect of the fractal dimension 𝐷 appears through the constant 𝜒 which combines the two scaling parameters: 𝑐_𝑥 and 𝑐_𝑧, that is, (𝜒=ln𝑐_𝑥/ln𝑐_𝑧) (cf. (2.13)). It is also clear in Tables that simple constant of proportionality between stress and strain does no longer exist.

Table 1 

Maxwell viscoelastic medium.

Table 2 

Kelvin-Voigt viscoelastic medium.

Table 3 

Standard Linear Solid viscoelastic medium.

Table 4 

Jeffrey viscoelastic medium.

6. Results and Discussion

Various continuous models for the creep approach of a viscoelastic surface, modeled by a fractal Cantor structure, in contact with a perfectly rigid foundation are listed in Tables 1–4. The models present an approximate closed form solution for the approach, (𝑢(𝑡)/ℎ₀), of the fractal surface as a function of the applied load, (𝑃₀/𝐿₀). Three sets of numerical parameters are required to evaluate such models as follows.(1)Those related to the material properties, that is, 𝐸_𝑣's, 𝜂's, and 𝜏: these parameters are usually adjusted to fit with the experimental results [54], and it could take different values related to the solids, that is, the magnitudes of such parameters are not unique among different publications.(2)Those related to the surface texture, that is, 𝑐_𝑥, 𝑐_𝑧, and 𝑠 that characterize the Cantor structure of the rough surface and relate through the fractal dimensions 𝐷 and 𝐷_𝑐, where 𝐷_𝑐 is the fractal dimension of the cantor set and (𝐷_𝑐=0.631) [9, 10].(3)Those extracted from the experimental results [60], that is, ℎ₀ which corresponding to twice the rms height. Another experimental result showed that the fractal dimension 𝐷 of a ground stainless steel surfaces is 𝐷=1.5 [16]. Values for 𝑃₀ and 𝐿₀ are not assumed, instead values for (𝜎=𝑃₀/𝐿₀) (for unit depth) are assumed.

It is important to mention here that the surface parameters like the fractal dimension 𝐷 along with the parameters 𝐿₀ and ℎ₀ can be determined experimentally from the surface of material specimen. The fractal dimension 𝐷 can be directly obtained from the slope of the structure function of surface profile or surface that can be presented by fractional Brownian process. The 𝐿₀ corresponds to the profile length (usually taken as autocorrelation length), and so the length 𝐿₀ could be taken to be in the range from 10 μm to 30 μm depending on the way how the surface is produced [16, 17]. In fact this range gives an idea about how small the size of voids when constructing the middle-third Cantor structure is.

The comparison-based experimental results cited in this paper are extracted from the work in [60]. These experiments were conducted on a carbon steel (0.45% carbon) specimen with surface roughness resulting from different finishing processes: face turning, grinding, and bead-blasting. The experiments were carried out 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. The validity of the analytical models presented in Tables 1–4 is motivated by comparing its results with the above-extracted experimental results [16, 60].

Analytical creep data for the steel considered in the form of strain, (𝑢/ℎ₀), versus the nondimensional time record, (𝑡/𝜏), are shown in Tables 1–4 and Figures 5–7.

(a) Maxwell medium

(b) Kelvin-Voigt medium

(c) Standard linear solid medium

(d) Jeffrey medium

Figure 5 

Non-dimensional time-strain curves for different stresses and fractal dimensions.

(a) Maxwell medium, 𝐷=1.24

(b) Maxwell medium, 𝐷=1.5

(c) Standard Linear Solid medium, 𝐷=1.4

(d) Standard Linear Solid medium, 𝐷=1.5

(e) Jeffrey medium, 𝐷=1.4

(f) Jeffrey medium, 𝐷=1.5

(g) Kelvin-Voigt medium, 𝐷=1.4

(h) Kelvin-Voigt medium, 𝐷=1.5

(i) Maxwell medium

(j) Kelvin-Voigt medium

Figure 6 

Experimental results and analytical isochronous stress-strain creep curves for constant room temperature and various (𝑡/𝜏).

(a) Kelvin-Voigt Medium, 𝐷=1.4, 𝑡/𝜏≈44, Δ𝑇=0, Δ𝑇=100, and Δ𝑇=200

(b) Kelvin-Voigt Medium, 𝐷=1.5, 𝑡/𝜏≈42, Δ𝑇=0, Δ𝑇=100, and Δ𝑇=200

(c) Standard Linear Solid, 𝐷=1.4, 𝑡/𝜏≈1680, Δ𝑇=0, Δ𝑇=100, and Δ𝑇=200

(d) Standard Linear Solid, 𝐷=1.5, 𝑡/𝜏≈1100, Δ𝑇=0, Δ𝑇=100, and Δ𝑇=200

Figure 7 

Experimental results and Analytical Isochronous Stress-Strain creep curves for constant (𝑡/𝜏) and various bulk temperatures.

6.1. Maxwell Model

For the Maxwell model to be representative of a solid rather than a fluid, the modulus 𝜂 (MPa⋅sec) must be large, comparable in magnitude with the modulus 𝐸_𝑣 (MPa), that is, (𝜏≡𝜂/𝐸_𝑣) is assumed to be of the order 1 (units of time; sec) [18]. In the confluent hypergeometric equation in Table 1, the value of 𝑛 should be large enough for the result to converge to a stable value; numerically it is assumed that (𝑛=550). The parameter (𝑐_𝑥=1.5) is held fixed, and two different values of the parameter {𝑐_𝑧=1.155,𝑐_𝑧=1.54} (𝑠=2) are used to yield two different fractal dimensions of the Cantor structure, 𝐷=1.5 and 𝐷=1.24, respectively. The depth ℎ₀ was taken as 9.7𝜇m which corresponds to twice r.m.s height obtained in [60] for the bead-blasted surface.

Stress-strain curves for various constant values of the time ratio (𝑡/𝜏), called isochronous stress-strain curves, in comparison with the experimental results are depicted in Figures 5–7 [18–24]. For (𝐷=1.24) the nondimensional time duration, (𝑡/𝜏), required to get an agreement between the experimental results and the isochronous curves varies between 100 and 200, while it varies between 10 and 40 for the fractal dimension (𝐷=1.5).

It is clear that a longer duration of agreement occurs for lower fractal dimension, which could be explained from the definition of the fractal dimension. Lower fractal dimensions means less roughness and consequently the bulk material dominates, in the contrary, higher fractal dimensions means excessive roughness and consequently the asperities deformation dominate.

Because of the periodicity of the Cantor set model, it undergoes the same construction procedure at each hierarchical level producing contact areas that are all of the same size. Therefore, it is doubtful that this model will provide an exact simulation of the deformation of random rough surface. However, the model does admit an analytic solution and, as proposed in [14], the specific character of a 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.

6.2. Kelvin-Voigt Model

In this study the modulus 𝐸_𝑣 is selected to be 130 GPa and the retardation time, (𝜏≡𝜂/𝐸_𝑣), is selected to be of the order “1” [61]. The value of 𝑛 is taken large enough for the result to converge to an accepted accuracy; it is assumed that (𝑛=550). The Cantor structure, shown in Figures 2 and 3, is built from the middle-third Cantor set where the parameter (𝑐_𝑥=1.5) is held fixed, giving a Cantor set dimension 𝐷_𝑐=0.63093 [9, 10]. Two different values of the parameter {𝑐_𝑧=1.155,𝑐_𝑧=1.29} with 𝑠=2 are used to yield two different fractal dimensions of the Cantor structure, 𝐷=1.5 and 𝐷=1.4, respectively. It was pointed out [14] that, only for (𝑐_𝑧≤2) the profile of the surface of the contacting body is fractal. The fractal dimension, 𝐷, of a ground stainless steel surfaces is 𝐷=1.5 [16, 43]. The value of ℎ₀, which corresponds to twice r.m.s height for the ground surface, is taken as 6.6 μm [60].

Figures 5–7 show a set of numerical creep data obtained applying the equations presented in Table 2. In these figures, the strain (𝑢/ℎ₀) is plotted versus the nondimensional time record of (𝑡/𝜏) for different ambient temperatures. As might be expected, higher strain rates occur for higher temperatures for a constant stress. Figure 5 shows the creep data for two different sets of fractal dimension; (1.4, 1.5) and (1.24, 1.5). It is clear that the higher strain rates results for higher fractal dimension. This 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.

Isochronous curves, for various constant values of the dimensionless time ratio (𝑡/𝜏) in comparison with the experimental results [60], conducted at room temperature, are depicted in Figures 6 and 7. These figures show a good agreement between the present 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, where it is clear that the relatively longer duration of agreement occurs for lower fractal dimension which could be attributed to the same reason mentioned above. The mathematical model shows instability when (𝑡/𝜏) exceeds 45 [19, 23].

6.3. Standard Linear Solid (SLS) Model

For SLS medium, as an illustrative example, the elastic moduli 𝐸_𝑣1 and 𝐸_𝑣2 are selected to be equal, they are given the value of 130 GPa, and the Newtonian viscosity, 𝜂₂, is selected to be of the same order, and consequently 𝜏=1 [61]. Three sets of numerical parameters are required to evaluate the equations in Table 3.(1)Those related to the material properties, that is, 𝐸_𝑣1, 𝐸_𝑣2, 𝜂₂, and 𝜏, these parameters are usually adjusted to fit with the experimental results [60], and it could take different values related to the solids, that is, the magnitudes of such parameters are not unique among different publications.(2)Those related to the surface texture, that is, 𝑐_𝑥, 𝑐_𝑧 and 𝑠 that characterize the Cantor structure of the rough surface and relate through the fractal dimensions 𝐷 and 𝐷_𝑐, where 𝐷_𝑐 is the fractal dimension of the cantor set and (𝐷_𝑐=0.631) [9, 10]. The parameter (𝑐_𝑥=1.5) is held fixed, and two different values of the parameter 𝑐_𝑧=1.15 and 𝑐_𝑧=1.29 with 𝑠=2 are used to yield two different fractal dimensions of the Cantor structure, (𝐷₁=1.5) and (𝐷₂=1.4), respectively (cf. (2.13)).(3)Those extracted from the experimental results [60], that is, ℎ₀=6.6𝜇m of the ground surface which corresponding to twice the rms height. Another experimental result showed that the fractal dimension 𝐷 of a ground stainless steel surfaces is (𝐷=1.5) [16]; values for 𝑃₀ and 𝐿₀ are not assumed, instead values for (𝜎=𝑃₀/𝐿₀), assuming unit depth, are assumed.

Figures 6 and 7 present the isochronous creep curves for different nondimensional time durations (𝑡/𝜏) accompanied by experimental results available in the literature [60]. The close agreement shown proves the validity of the analytical model presented in this work. Since a standard linear solid constitutive relation is adopted, the ratio (𝑡/𝜏) as seen in Figures 6 and 7 is large enough, which is in accordance with the creep of real metallic materials. The model presented is used in order to conduct an analytical investigation of the solution of contact of a viscoelastic problem, but nevertheless it might turn out not to be too far from reality since tests on the actual contact area of steel surfaces show that they contain waviness and tortuosity at different scales [22, 24].

6.4. Jeffrey Model

To synthesize the required fractal surface, parameter (𝑐_𝑥=1.5) is held fixed, and two different values of the parameter 𝑐_𝑧, that is, {𝑐_𝑧=1.155,𝑐_𝑧=1.29} yields two different dimensions of the Cantor structure, that is, 𝐷=1.5 and 𝐷=1.4, respectively. The double r.m.s. value of the heights of roughness ℎ₀ is taken as 9.7 μm [60].

The material parameters {𝜂₁,𝜂₂,𝐸_𝑣2} are usually identified such that the mathematical model derived can fit the experimental results. For the Jeffreys' model to be representative of a solid rather than a fluid, the moduli 𝜂₁ and 𝜂₂ (MPa⋅sec) must be large, comparable in magnitude with the modulus 𝐸_𝑣2 [11]. In this work, and as an empirical example, the retardation time 𝜏 is assumed to be of the order “1”, and 𝐸_𝑣=210 GPa [61]. The activation energy 𝑄 is taken as 7.3 KJ/mol [62]. Figures 5‐ 7 show a set of numerical creep data that are obtained by applying the equations presented in Table 4.

Figures 6 and 7 show the isochronous curves for creep data for different nondimensional time duration, (𝑡/𝜏) in comparison with the experimental results available in the open literature [20]. This model is used in order to conduct an exact analytic investigation of the solution of contact problems, but nevertheless it might turn out not to be too far from reality, since tests on the actual contact area of steel surfaces show that they contain waviness and tortuosity at different scales. It is also clear that quantitative reproduction of real material behavior requires a nonlinear viscoelastic stress-strain relation.

7. Conclusions and Future Work

This paper presents recent advancements in the creep contact of various linear viscoelastic models; Maxwell model, thermal and isothermal Kelvin-Voigt model, Jeffrey model, and thermal and isothermal standard linear solid (SLS) model. Such models were used to study the effect of the normal creep load, the surface quality, and/or the bulk temperature on the creep process of a nominally flat surface in contact with a rigid half space. Since these models are constituted from combinations of dashpots and springs, that is, not real materials, it is not expected to predict precisely the real material behavior for all situations. In fact they fail to predict the tertiary region, which means that quantitative reproduction of real material behavior requires a nonlinear viscoelastic stress-strain constitutive relation.

The effects of temperatures on the behavior of such mathematical models are shown to follow trend similar to the isothermal curve. Such behavior is a consequence of assuming thermorheologically simple materials. Also, because of the periodicity of the Cantor set model, it undergoes the same construction procedure at each hierarchical level producing contact areas that are all of the same size. Therefore, it is doubtful that these models will provide an exact simulation of the deformation of random rough surface. However, these models do admit an analytical solution, where the specific character of a fractal models 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 solutions obtained here provide further insight into the effect that surface texture has on the deformation process, and it also provides indications of the effect that different surface forming processes (since they produce different surface textures) may have on subsequent surface deformation.

For future work, further advancements in the study of rough surfaces might be achieved if the random field based on the generalized Cauchy (GC) correlation structure will be used. In addition, further investigation for the prediction of surface roughness, similar to the work in [40] which discusses the predictability of fractal time series, will be performed. Also, it is intended to extend the methodology which is used in this paper for the applications of contact-surfaces within the micro- and nanoelectronic (and electrical) devices and circuits, such as resistors, capacitors and inductors that usually occur in electronic manufacturing. Further investigation of fractals within quantum computing applications will also be conducted.