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

Authors: Osama M. Abuzeid; Anas N. Al-Rabadi; Hashem S. Alkhaldi
Journal: Mathematical Problems in Engineering
 (2011)

---

## 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.


---

## Body

## 1. Introduction

In this paper, Section1 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, σm2, and  (3) the root mean square of curvatures, σk2. 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σm2 and σk2 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: G(ω)~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 H (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 D, and the Sierpinski gasket (triangle) [6–10]. One may note that the representations of the fractal dimension D and the Hurst index H in [6–8] differ from those in [9, 10]; in [9, 10], the representations for D and H are based on the fractional Brownian motion (fBm), while in [6–8] D and H 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.Common fractals: (a) the middle-third Cantor set, (b) the von Koch curve, (c) the Weierstrass-Mandelbrot functionC in the range 0≤x≤3 (p=1.5 and γ=0.5), where p and γ are two numerical parameters (cf. (1.17)), where the trend of the function is ~x2, (d) trajectories of a fractional Brownian process for different H and D, 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.
(a)(b)(c)(d)(e)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 X be a compact metric space and O be the totality of open balls in X. The Hausdorff s-measure of a subset S⊂X which is defined for s≥0 as the following limit:
(1.1)mH(S,s)=lim⁡σ→0inf⁡G∈O{∑V∈G(dim⁡V)s:S⊆⋃V∈GV,diam⁡V≤δ}.HereG is finite or denumerable subset of O. It was proven that there exists a value s0 such that:
(1.2)mH(S,s)={∞,fors<s0,0,fors>s0.The Hausdorff dimension of the setS, denoted by dimHS, is the number s0 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.LetE be the Euclidean dimension of the space in which a set S is embedded. For δ>0, let N(δ) be the smallest number of E-dimensional balls or cubes of diameter d needed to cover the set S. The box counting dimension or box dimension, denoted by dimBS, can be defined if the following limit exists:
(1.3)dimBS=lim⁡δ→0+log⁡N(δ)-log⁡δ.It can be proven thatdimBS does not change if one takes N(δ) to be either the smallest number of δ-cubes that cover S; or the number of δ-mesh cubes that intersect S; or the smallest number of sets of diameter at most d that cover S; or the largest number of disjoint δ-balls with centers in S. Unfortunately, the box dimension is not always equal to the Hausdorff dimension. For example, the set S={0,1,1/2,…,1/n,…} has unequal values for the Hausdorff and box dimensions for dimHS≠dimBS=1/2. However, it can be proven that dimHS≤dimBS.As a simple alternative to the Hausdorff measure, we can introduce thes-measure ms of a set as the following limit:
(1.4)ms(S)=lim⁡δ→0+N(δ)δs
and define the box dimension as the value s=D such that ms(S) has a jump from 0 to ∞ similar to the behavior of mH(S,s) 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 D when the model considers scales R such that δ*<R<Δ*, where δ* and Δ* are the upper and lower cutoffs for the fractal representation. To get the value D of the dimension, the considered region is discretized into cubes with side length δ*. Then the smallest number of E-dimensional cubes needed to cover the cluster (N(δ*)) is counted. One says that the cluster is fractal if the numbers N(δ*) satisfy the so-called number-radius relation for different sizes of the considered region of the cluster R as follows:
(1.5)N(δ*)≈(Rδ*)D,forδ*<R<Δ*.The value ofD is estimated as the slope of linear growth of ln(N(δ*)) plotted against ln(R). The power D 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 setS it is proven that dimBS 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 mappingM of a plane π onto a plane π′ is called a similarity mapping with coefficient λ>0, or simply a similarity, when the following property holds: if {A,B} are any two points of π, and {A′,B′} are their images under M, then |A′B′|=λ|AB| [35]. It is known that any similarity transformation of a plane is a homogeneous (isotropic) dilation of coordinates {x′=λx,z′=λz} up to a rotation and translation. A set S 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 S′=λS.In practice, it is impossible to verify that all statistical moments of the two distributions are identical. Frequently, a setS is said to be self-similar if only a few moments do not change under scaling [36]. A one-to-one mapping M 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 x′=λx and z′=λHz. Here z is a graph of a surface profile and H 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 S does not change under the action of the Lipschitz homeomorphism L as follows:
(1.6)dimHS=dimHL(S).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 of1-D (fBm) VH(x) of index H, where x is taken as the time and z is the random variable of the single valued function VH(x) with the following property:
(1.7)〈[VH(x+δ)-VH(x)]2〉~δ2H,for0<H<1,
where 〈〉 denotes averaging over the ensemble, and H is the Hurst index. The scaling behavior of the different traces, VH(x), is characterized by a particular H which relate the typical change in Δz(x), where z(x)+VH(x), is the trace of the fBm, and the change in the spatial coordinate Δx by the simple scaling law [36, 39, 40]:
(1.8)Δz(x)~ΔxH.It is known that, with probability equal to “1”, the following holds [34]:
(1.9)dimHVH(x)=dimBVH(x)=2-H.
The autocorrelation function is one of the main tools for studying statistical models of rough surfaces. The autocorrelation function R(δ) of the profile is expressed as:(1.10)R(δ)=lim⁡T→∞12T∫-TT[z(x+δ)-z̅][z(x)-z̅]dx=〈[z(x+δ)-z̅][z(x)-z̅]〉,
which also can be expressed as:
(1.11)R(δ)=lim⁡T→∞12T∫-TTz(x+δ)z(x)dx-(z̅)2,
where z̅ is the average value of the profile function z(x).Another tool for the characterization of surfaces is the spectral density functionG(ω) which is the Fourier transform of R(δ):
(1.12)G(ω)=2π∫0∞R(δ)cos⁡ωδdδ,R(δ)=2π∫0∞G(ω)cos⁡ωδdω.In general, the following is accepted in fractional Brownian motion [2] (i)
If the auto-correlation functionR(δ) of the profile z(x) satisfies R(0)-R(δ)~δ2(2-s), then it is reasonable to expect that the box dimension of the graph z(x) is equal to s, note that one can find R(0)-R(δ)~δ2H for the fBm defined by (1.7).(ii)
If the profilez(x) has spectral density:
(1.13)G(ω)~1ωυ,
then it is reasonable to expect that the box dimension of the graph z(x) is equal to (5-υ)/2 [2]. The above conclusions are valid for mathematical models of the profile, for which the relation 2(5-s)=υ-1 or υ=5-2s 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 (mn) of the spectral density G(ω) provide a useful description of the surface roughness:
(1.14)mn=∫ω0∞ωnG(ω)dω,
where (ω0=2π/λ0) is the wave number corresponding to the profile length λ0. It is possible to show that m0 is the variance of heights (r.m.s height) of the surface, m1 is the variance of slopes (r.m.s slope), and m3 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: (1.15)W(x;p)=∑n=0∞p-γnh(pnx),forp>1,0<γ<1,
where h is a bounded Hölder function of order greater than β. The following complex generalization of the W(x;p) was considered:(1.16)W̃(x;p)=∑n=-∞∞p-γn[(1-eipnt)eiΦn],forp>1,0<γ<1,
where Φn 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) functionC and the Takagi-Hopson function T:(1.17)C(x;p)=∑n=-∞∞p-γn(1-cos⁡⁡pnx),forp>1,0<γ<1,(1.18)T(x;p)=∑n=-∞∞p-γn|pnx-[pnx+12]|,forp>1,0<γ<1,
have been studied in numerous papers [3, 18, 35, 37]. By direct calculations, one may obtain:
(1.19)|W̃(x+δ;p)-W̃(x;p)|~δγ,
which is similar to the behavior of (1.7) of fractional Brownian motion. The box dimension of the Weierstrass function graphs is D=2-γ and it is believed that their Hausdorff dimension is the same [34, 43]. Currently, the only known bounds for the Hausdorff dimensions are D-(c/logp)≤dimHgraphC≤D, provided that p is large and constant c is large enough [25]. It is possible to calculate the spectral density of the WM function W̃(x;p) as follows:(1.20)G(ω)=∑n=-∞∞δ(ω-pn)p2(2-D)n,
where δ is the Dirac delta function. Some arguments for approximating this discrete spectral density by a continuous spectral density G̃(ω)~1/ω5-2D, whose exponent (5-2D) is in agreement with (1.13) with respect to the box dimension were suggested [3]. The following truncated WM function:
(1.21)W̃1(x;p)=A(D-1)∑n=n1∞p(D-2)ncos⁡2πpnx
is often used for fractal characterization of the surface topography [4, 37, 38]. Here n1 is an integer number, which corresponds to the low cut-off frequency of the profile, and A is the so-called characteristic length scale of the profile. The number n1 depends on the length L of the sample and is given by (pn1=1/L) and the parameter A determines the position of the spectral density along the logG axis. It was stated that both parameters A and D of the function W1 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 D are not unique and depend on instruments or resolution of a given instrument [2]. Evidently, the function C(x;p) is not homogeneous. Nevertheless, it exhibits the property C(pkx;p)=pkγC(x;p), with k∈Z where Z is the set of all integers which looks similar to the definition of a homogeneous function hd of degree d, that is, hd(λx)=λdhd(x) for λ>0.Thus, the graph of the functionC(x;p) near any point x0 is repeated in scaling form near all points pkx0, k∈Z. 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 bd(pkx;p)=pkdbd(x;p), with k∈Z where d is degree of homogeneity. As examples of 1-dimensional fractal PH-curves we can consider the graphs of functions b1 and b2 with degrees d=1 and d=2, respectively:
(1.22)b0(x;p)=x-γC(x;p),b1(x;p)=xb0(x;p),b2(x;p)=x2b0(x;p).Because of (1.6), these functions have the same Hausdorff dimension as the WM function C(x;p) whose box-dimension is D. Another consequence is that the WM function C(x;p), with C(x;p)~x2-D 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].


## 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, σm2, and  (3) the root mean square of curvatures, σk2. 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σm2 and σk2 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: G(ω)~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 H (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 D, and the Sierpinski gasket (triangle) [6–10]. One may note that the representations of the fractal dimension D and the Hurst index H in [6–8] differ from those in [9, 10]; in [9, 10], the representations for D and H are based on the fractional Brownian motion (fBm), while in [6–8] D and H 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.Common fractals: (a) the middle-third Cantor set, (b) the von Koch curve, (c) the Weierstrass-Mandelbrot functionC in the range 0≤x≤3 (p=1.5 and γ=0.5), where p and γ are two numerical parameters (cf. (1.17)), where the trend of the function is ~x2, (d) trajectories of a fractional Brownian process for different H and D, 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.
(a)(b)(c)(d)(e)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 X be a compact metric space and O be the totality of open balls in X. The Hausdorff s-measure of a subset S⊂X which is defined for s≥0 as the following limit:
(1.1)mH(S,s)=lim⁡σ→0inf⁡G∈O{∑V∈G(dim⁡V)s:S⊆⋃V∈GV,diam⁡V≤δ}.HereG is finite or denumerable subset of O. It was proven that there exists a value s0 such that:
(1.2)mH(S,s)={∞,fors<s0,0,fors>s0.The Hausdorff dimension of the setS, denoted by dimHS, is the number s0 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.LetE be the Euclidean dimension of the space in which a set S is embedded. For δ>0, let N(δ) be the smallest number of E-dimensional balls or cubes of diameter d needed to cover the set S. The box counting dimension or box dimension, denoted by dimBS, can be defined if the following limit exists:
(1.3)dimBS=lim⁡δ→0+log⁡N(δ)-log⁡δ.It can be proven thatdimBS does not change if one takes N(δ) to be either the smallest number of δ-cubes that cover S; or the number of δ-mesh cubes that intersect S; or the smallest number of sets of diameter at most d that cover S; or the largest number of disjoint δ-balls with centers in S. Unfortunately, the box dimension is not always equal to the Hausdorff dimension. For example, the set S={0,1,1/2,…,1/n,…} has unequal values for the Hausdorff and box dimensions for dimHS≠dimBS=1/2. However, it can be proven that dimHS≤dimBS.As a simple alternative to the Hausdorff measure, we can introduce thes-measure ms of a set as the following limit:
(1.4)ms(S)=lim⁡δ→0+N(δ)δs
and define the box dimension as the value s=D such that ms(S) has a jump from 0 to ∞ similar to the behavior of mH(S,s) 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 D when the model considers scales R such that δ*<R<Δ*, where δ* and Δ* are the upper and lower cutoffs for the fractal representation. To get the value D of the dimension, the considered region is discretized into cubes with side length δ*. Then the smallest number of E-dimensional cubes needed to cover the cluster (N(δ*)) is counted. One says that the cluster is fractal if the numbers N(δ*) satisfy the so-called number-radius relation for different sizes of the considered region of the cluster R as follows:
(1.5)N(δ*)≈(Rδ*)D,forδ*<R<Δ*.The value ofD is estimated as the slope of linear growth of ln(N(δ*)) plotted against ln(R). The power D 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 setS it is proven that dimBS 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 mappingM of a plane π onto a plane π′ is called a similarity mapping with coefficient λ>0, or simply a similarity, when the following property holds: if {A,B} are any two points of π, and {A′,B′} are their images under M, then |A′B′|=λ|AB| [35]. It is known that any similarity transformation of a plane is a homogeneous (isotropic) dilation of coordinates {x′=λx,z′=λz} up to a rotation and translation. A set S 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 S′=λS.In practice, it is impossible to verify that all statistical moments of the two distributions are identical. Frequently, a setS is said to be self-similar if only a few moments do not change under scaling [36]. A one-to-one mapping M 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 x′=λx and z′=λHz. Here z is a graph of a surface profile and H 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 S does not change under the action of the Lipschitz homeomorphism L as follows:
(1.6)dimHS=dimHL(S).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 of1-D (fBm) VH(x) of index H, where x is taken as the time and z is the random variable of the single valued function VH(x) with the following property:
(1.7)〈[VH(x+δ)-VH(x)]2〉~δ2H,for0<H<1,
where 〈〉 denotes averaging over the ensemble, and H is the Hurst index. The scaling behavior of the different traces, VH(x), is characterized by a particular H which relate the typical change in Δz(x), where z(x)+VH(x), is the trace of the fBm, and the change in the spatial coordinate Δx by the simple scaling law [36, 39, 40]:
(1.8)Δz(x)~ΔxH.It is known that, with probability equal to “1”, the following holds [34]:
(1.9)dimHVH(x)=dimBVH(x)=2-H.
The autocorrelation function is one of the main tools for studying statistical models of rough surfaces. The autocorrelation function R(δ) of the profile is expressed as:(1.10)R(δ)=lim⁡T→∞12T∫-TT[z(x+δ)-z̅][z(x)-z̅]dx=〈[z(x+δ)-z̅][z(x)-z̅]〉,
which also can be expressed as:
(1.11)R(δ)=lim⁡T→∞12T∫-TTz(x+δ)z(x)dx-(z̅)2,
where z̅ is the average value of the profile function z(x).Another tool for the characterization of surfaces is the spectral density functionG(ω) which is the Fourier transform of R(δ):
(1.12)G(ω)=2π∫0∞R(δ)cos⁡ωδdδ,R(δ)=2π∫0∞G(ω)cos⁡ωδdω.In general, the following is accepted in fractional Brownian motion [2] (i)
If the auto-correlation functionR(δ) of the profile z(x) satisfies R(0)-R(δ)~δ2(2-s), then it is reasonable to expect that the box dimension of the graph z(x) is equal to s, note that one can find R(0)-R(δ)~δ2H for the fBm defined by (1.7).(ii)
If the profilez(x) has spectral density:
(1.13)G(ω)~1ωυ,
then it is reasonable to expect that the box dimension of the graph z(x) is equal to (5-υ)/2 [2]. The above conclusions are valid for mathematical models of the profile, for which the relation 2(5-s)=υ-1 or υ=5-2s 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 (mn) of the spectral density G(ω) provide a useful description of the surface roughness:
(1.14)mn=∫ω0∞ωnG(ω)dω,
where (ω0=2π/λ0) is the wave number corresponding to the profile length λ0. It is possible to show that m0 is the variance of heights (r.m.s height) of the surface, m1 is the variance of slopes (r.m.s slope), and m3 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: (1.15)W(x;p)=∑n=0∞p-γnh(pnx),forp>1,0<γ<1,
where h is a bounded Hölder function of order greater than β. The following complex generalization of the W(x;p) was considered:(1.16)W̃(x;p)=∑n=-∞∞p-γn[(1-eipnt)eiΦn],forp>1,0<γ<1,
where Φn 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) functionC and the Takagi-Hopson function T:(1.17)C(x;p)=∑n=-∞∞p-γn(1-cos⁡⁡pnx),forp>1,0<γ<1,(1.18)T(x;p)=∑n=-∞∞p-γn|pnx-[pnx+12]|,forp>1,0<γ<1,
have been studied in numerous papers [3, 18, 35, 37]. By direct calculations, one may obtain:
(1.19)|W̃(x+δ;p)-W̃(x;p)|~δγ,
which is similar to the behavior of (1.7) of fractional Brownian motion. The box dimension of the Weierstrass function graphs is D=2-γ and it is believed that their Hausdorff dimension is the same [34, 43]. Currently, the only known bounds for the Hausdorff dimensions are D-(c/logp)≤dimHgraphC≤D, provided that p is large and constant c is large enough [25]. It is possible to calculate the spectral density of the WM function W̃(x;p) as follows:(1.20)G(ω)=∑n=-∞∞δ(ω-pn)p2(2-D)n,
where δ is the Dirac delta function. Some arguments for approximating this discrete spectral density by a continuous spectral density G̃(ω)~1/ω5-2D, whose exponent (5-2D) is in agreement with (1.13) with respect to the box dimension were suggested [3]. The following truncated WM function:
(1.21)W̃1(x;p)=A(D-1)∑n=n1∞p(D-2)ncos⁡2πpnx
is often used for fractal characterization of the surface topography [4, 37, 38]. Here n1 is an integer number, which corresponds to the low cut-off frequency of the profile, and A is the so-called characteristic length scale of the profile. The number n1 depends on the length L of the sample and is given by (pn1=1/L) and the parameter A determines the position of the spectral density along the logG axis. It was stated that both parameters A and D of the function W1 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 D are not unique and depend on instruments or resolution of a given instrument [2]. Evidently, the function C(x;p) is not homogeneous. Nevertheless, it exhibits the property C(pkx;p)=pkγC(x;p), with k∈Z where Z is the set of all integers which looks similar to the definition of a homogeneous function hd of degree d, that is, hd(λx)=λdhd(x) for λ>0.Thus, the graph of the functionC(x;p) near any point x0 is repeated in scaling form near all points pkx0, k∈Z. 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 bd(pkx;p)=pkdbd(x;p), with k∈Z where d is degree of homogeneity. As examples of 1-dimensional fractal PH-curves we can consider the graphs of functions b1 and b2 with degrees d=1 and d=2, respectively:
(1.22)b0(x;p)=x-γC(x;p),b1(x;p)=xb0(x;p),b2(x;p)=x2b0(x;p).Because of (1.6), these functions have the same Hausdorff dimension as the WM function C(x;p) whose box-dimension is D. Another consequence is that the WM function C(x;p), with C(x;p)~x2-D 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].


## 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 X be a compact metric space and O be the totality of open balls in X. The Hausdorff s-measure of a subset S⊂X which is defined for s≥0 as the following limit:
(1.1)mH(S,s)=lim⁡σ→0inf⁡G∈O{∑V∈G(dim⁡V)s:S⊆⋃V∈GV,diam⁡V≤δ}.HereG is finite or denumerable subset of O. It was proven that there exists a value s0 such that:
(1.2)mH(S,s)={∞,fors<s0,0,fors>s0.The Hausdorff dimension of the setS, denoted by dimHS, is the number s0 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.LetE be the Euclidean dimension of the space in which a set S is embedded. For δ>0, let N(δ) be the smallest number of E-dimensional balls or cubes of diameter d needed to cover the set S. The box counting dimension or box dimension, denoted by dimBS, can be defined if the following limit exists:
(1.3)dimBS=lim⁡δ→0+log⁡N(δ)-log⁡δ.It can be proven thatdimBS does not change if one takes N(δ) to be either the smallest number of δ-cubes that cover S; or the number of δ-mesh cubes that intersect S; or the smallest number of sets of diameter at most d that cover S; or the largest number of disjoint δ-balls with centers in S. Unfortunately, the box dimension is not always equal to the Hausdorff dimension. For example, the set S={0,1,1/2,…,1/n,…} has unequal values for the Hausdorff and box dimensions for dimHS≠dimBS=1/2. However, it can be proven that dimHS≤dimBS.As a simple alternative to the Hausdorff measure, we can introduce thes-measure ms of a set as the following limit:
(1.4)ms(S)=lim⁡δ→0+N(δ)δs
and define the box dimension as the value s=D such that ms(S) has a jump from 0 to ∞ similar to the behavior of mH(S,s) 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 D when the model considers scales R such that δ*<R<Δ*, where δ* and Δ* are the upper and lower cutoffs for the fractal representation. To get the value D of the dimension, the considered region is discretized into cubes with side length δ*. Then the smallest number of E-dimensional cubes needed to cover the cluster (N(δ*)) is counted. One says that the cluster is fractal if the numbers N(δ*) satisfy the so-called number-radius relation for different sizes of the considered region of the cluster R as follows:
(1.5)N(δ*)≈(Rδ*)D,forδ*<R<Δ*.The value ofD is estimated as the slope of linear growth of ln(N(δ*)) plotted against ln(R). The power D 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 setS it is proven that dimBS 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 mappingM of a plane π onto a plane π′ is called a similarity mapping with coefficient λ>0, or simply a similarity, when the following property holds: if {A,B} are any two points of π, and {A′,B′} are their images under M, then |A′B′|=λ|AB| [35]. It is known that any similarity transformation of a plane is a homogeneous (isotropic) dilation of coordinates {x′=λx,z′=λz} up to a rotation and translation. A set S 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 S′=λS.In practice, it is impossible to verify that all statistical moments of the two distributions are identical. Frequently, a setS is said to be self-similar if only a few moments do not change under scaling [36]. A one-to-one mapping M 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 x′=λx and z′=λHz. Here z is a graph of a surface profile and H 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 S does not change under the action of the Lipschitz homeomorphism L as follows:
(1.6)dimHS=dimHL(S).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 of1-D (fBm) VH(x) of index H, where x is taken as the time and z is the random variable of the single valued function VH(x) with the following property:
(1.7)〈[VH(x+δ)-VH(x)]2〉~δ2H,for0<H<1,
where 〈〉 denotes averaging over the ensemble, and H is the Hurst index. The scaling behavior of the different traces, VH(x), is characterized by a particular H which relate the typical change in Δz(x), where z(x)+VH(x), is the trace of the fBm, and the change in the spatial coordinate Δx by the simple scaling law [36, 39, 40]:
(1.8)Δz(x)~ΔxH.It is known that, with probability equal to “1”, the following holds [34]:
(1.9)dimHVH(x)=dimBVH(x)=2-H.
The autocorrelation function is one of the main tools for studying statistical models of rough surfaces. The autocorrelation function R(δ) of the profile is expressed as:(1.10)R(δ)=lim⁡T→∞12T∫-TT[z(x+δ)-z̅][z(x)-z̅]dx=〈[z(x+δ)-z̅][z(x)-z̅]〉,
which also can be expressed as:
(1.11)R(δ)=lim⁡T→∞12T∫-TTz(x+δ)z(x)dx-(z̅)2,
where z̅ is the average value of the profile function z(x).Another tool for the characterization of surfaces is the spectral density functionG(ω) which is the Fourier transform of R(δ):
(1.12)G(ω)=2π∫0∞R(δ)cos⁡ωδdδ,R(δ)=2π∫0∞G(ω)cos⁡ωδdω.In general, the following is accepted in fractional Brownian motion [2] (i)
If the auto-correlation functionR(δ) of the profile z(x) satisfies R(0)-R(δ)~δ2(2-s), then it is reasonable to expect that the box dimension of the graph z(x) is equal to s, note that one can find R(0)-R(δ)~δ2H for the fBm defined by (1.7).(ii)
If the profilez(x) has spectral density:
(1.13)G(ω)~1ωυ,
then it is reasonable to expect that the box dimension of the graph z(x) is equal to (5-υ)/2 [2]. The above conclusions are valid for mathematical models of the profile, for which the relation 2(5-s)=υ-1 or υ=5-2s 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 (mn) of the spectral density G(ω) provide a useful description of the surface roughness:
(1.14)mn=∫ω0∞ωnG(ω)dω,
where (ω0=2π/λ0) is the wave number corresponding to the profile length λ0. It is possible to show that m0 is the variance of heights (r.m.s height) of the surface, m1 is the variance of slopes (r.m.s slope), and m3 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: (1.15)W(x;p)=∑n=0∞p-γnh(pnx),forp>1,0<γ<1,
where h is a bounded Hölder function of order greater than β. The following complex generalization of the W(x;p) was considered:(1.16)W̃(x;p)=∑n=-∞∞p-γn[(1-eipnt)eiΦn],forp>1,0<γ<1,
where Φn 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) functionC and the Takagi-Hopson function T:(1.17)C(x;p)=∑n=-∞∞p-γn(1-cos⁡⁡pnx),forp>1,0<γ<1,(1.18)T(x;p)=∑n=-∞∞p-γn|pnx-[pnx+12]|,forp>1,0<γ<1,
have been studied in numerous papers [3, 18, 35, 37]. By direct calculations, one may obtain:
(1.19)|W̃(x+δ;p)-W̃(x;p)|~δγ,
which is similar to the behavior of (1.7) of fractional Brownian motion. The box dimension of the Weierstrass function graphs is D=2-γ and it is believed that their Hausdorff dimension is the same [34, 43]. Currently, the only known bounds for the Hausdorff dimensions are D-(c/logp)≤dimHgraphC≤D, provided that p is large and constant c is large enough [25]. It is possible to calculate the spectral density of the WM function W̃(x;p) as follows:(1.20)G(ω)=∑n=-∞∞δ(ω-pn)p2(2-D)n,
where δ is the Dirac delta function. Some arguments for approximating this discrete spectral density by a continuous spectral density G̃(ω)~1/ω5-2D, whose exponent (5-2D) is in agreement with (1.13) with respect to the box dimension were suggested [3]. The following truncated WM function:
(1.21)W̃1(x;p)=A(D-1)∑n=n1∞p(D-2)ncos⁡2πpnx
is often used for fractal characterization of the surface topography [4, 37, 38]. Here n1 is an integer number, which corresponds to the low cut-off frequency of the profile, and A is the so-called characteristic length scale of the profile. The number n1 depends on the length L of the sample and is given by (pn1=1/L) and the parameter A determines the position of the spectral density along the logG axis. It was stated that both parameters A and D of the function W1 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 D are not unique and depend on instruments or resolution of a given instrument [2]. Evidently, the function C(x;p) is not homogeneous. Nevertheless, it exhibits the property C(pkx;p)=pkγC(x;p), with k∈Z where Z is the set of all integers which looks similar to the definition of a homogeneous function hd of degree d, that is, hd(λx)=λdhd(x) for λ>0.Thus, the graph of the functionC(x;p) near any point x0 is repeated in scaling form near all points pkx0, k∈Z. 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 bd(pkx;p)=pkdbd(x;p), with k∈Z where d is degree of homogeneity. As examples of 1-dimensional fractal PH-curves we can consider the graphs of functions b1 and b2 with degrees d=1 and d=2, respectively:
(1.22)b0(x;p)=x-γC(x;p),b1(x;p)=xb0(x;p),b2(x;p)=x2b0(x;p).Because of (1.6), these functions have the same Hausdorff dimension as the WM function C(x;p) whose box-dimension is D. Another consequence is that the WM function C(x;p), with C(x;p)~x2-D 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, Figures2 and 3, where L0 corresponds to the nominal surface length, and h0 is equal to twice the r.m.s height of the roughness.Figure 2
The fractal middle-third Cantor structures=2, where E0 is the initiator step, {E1,E2,E3} 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.Figure 3
The fractal middle-third Cantor structures=3, where E0 is the initiator step, {E1,E2} 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.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, cx and cz, 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/cx) times the length of the initial segment where cx>1. The depth of the recesses (measured from the last step) at the (i+1)th construction step of the fractal surface is (1/cz) times less than the depth of the ith step where cz>1. From this it can be easily shown that the horizontal length and recess depth of the (i+1)th step are:(2.1)Li+1=cx-1Li=cx-(i+1)L0,(2.2)hi+1=cz-1hi=cz-(i+1)h0,
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 i iterations, is found to be:(2.3)L(i)=L0+2h0(2cz-(i+1)-1)(2cz-1-1).It is obvious that only forcz≤2 the contour of the surface of the structure will be a fractal, because in this case the length L(i) 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 ith generation, the Cantor structure contains N=si segments, each of length [14]:
(2.4)δi=(1scx)iL0,
where the parameter s corresponds to the number of asperities on a repeating segment [16, 17].In Figure2 we have s=2, and in Figure 3 we use s=3. Therefore, by changing the value of s, an infinite number of different structures could be constructed.The profile of the surface in Figures2 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:
(2.5)Δxi+1=(scx)-1Δxi,
while in the vertical direction, the corresponding fluctuations Δzi in the ith generation can be defined by considering the probability of obtaining the value of:
(2.6)zi=cz-ih0.The fluctuationΔzi in the ith generation can be obtained by assuming that Δzi scales as the expected value ziP(zi) [14] where:
(2.7)Δzi∝ziP(zi),
and P(zi) is the probability of obtaining the value zi which is given by:
(2.8)P(zi)=(Li-Li+1)L0,
where it is also produced as:
(2.9)P(zi)=(cx)-i(1-cx-1).Thus, the expected value of the fluctuation in the(i+1)th generation is related to the expected value of the fluctuation in the ith generation through:
(2.10)zi+1P(zi+1)=(cxcz)-iziP(zi).
Hence, we have:
(2.11)Δzi+1=(cxcz)-iΔzi.Substituting (2.5) and (2.11) into (1.8), we get:(2.12)Δzi+1Δzi=(Δxi+1Δxi)2-D,
from which the self-affine fractal dimension for the contour of the Cantor structure is:(2.13)D=1+ln⁡sln⁡(scx)-ln⁡czln⁡(scx)=1+Dc-ln⁡czln⁡(scx),for1<D<2,
where (0<Dc<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;czh0 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 Figures2 and 3, there exist a series of one-dimensional elastic bars distributed in a way that the distance from the initiator step E0 to the generated step E3 is indicated by h0, from E1 to E3 is indicated by h1, from E2 to E3 is indicated by h2, and so forth, [18–24].LetF3 be the force required to compress E3 until E2, F2 be the force require to compress E3 and E2 until E1, and F1 be the force required to compress E3, E2 and E1 until E0, and so forth, or in general:
(3.1)ΔFi+1=EL0cz(cz-1)cz-icx-i.Also, letFi+1 be the limit force for protrusion of the (i+1)th generation. We assume that when the limit load is reached, the surface approaches a distance Δui+1 equal to the difference between the heights protrusion of ith and (i+1)th generations. We will assume that E2 is deflected by a distance u2=h1, E1 is deflected by a distance u1=h0, E0 is deflected a distance u0=czh0, and so Δu1=u0-u1=h0(cz-1) and Δu2=h0cz-1(cz-1). Thus:
(3.2)Δui+1=h0(cz-1)cz-i.The above-mentioned assumptions are sufficient to determine the dependence of the limit loadF on the approach u. We will use the fact that the surfaces were approached by an amount Δui+1 when the limit remote increases from Fi+1 to Fi, then, for load effect only:
(3.3)ΔFi+1Δui+1=EL0czh0cx-i.In the limit wheni→∞, the following asymptotic behavior is obtained as:
(3.4)ε=uh0=(czχEFL0)1/(χ+1),
where ε=u/h0 is the strain due to the remote applied force and χ=lncx/lncz.


## 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 factoraT 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 Tref. If the creep time obeys an Arrhenius relation, the shift factor can be shown to be [57]:
(4.1)log⁡aT=Q2.303R(1T-1Tref),
where Q is the activation energy (J/mol), R is the gas constant (J/mol K), and T 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)ε=f(σ,t,T).One way to simplify this function is to separate it into three functions of stress, time and temperature as follow [53, 58]:
(4.3)ε=f1(σ)f2(t)f3(T).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:
(4.4)f3(T)=Bexp⁡(QRT),
where B is a constant given by:(4.5)B=exp⁡(-QRT0).When combining (4.4) and (4.5), it could be shown that at the reference temperature T0 the function f3(T) 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 E's are the Newtonian viscosity and the elastic modulus, respectively.Various linear viscoelastic models, whereηis the Newtonian viscosity, E is the elastic modulus, u is the approach and F is the applied external remote load.
(a)(b)(c)(d)Creep test is a widely used standard test, wherein a forceP0 is suddenly applied at time t=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[U(t)], that is, F=P0[U(t)]. 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.Tables1–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 E, and the modified creep model, where (u(t)/h0) is the relative approach, (P0/L0) is the applied stress per unit depth, T is the bulk temperature, T0 is the reference temperature, η is the Newtonian viscosity, Ev is the elastic modulus, Q is the activation energy, R is the ideal gas constant, and (τ≡η/Ev) is a characteristics parameter with units of time called the retardation time [18–24]. The effect of the fractal dimension D appears through the constant χ which combines the two scaling parameters: cx and cz, that is, (χ=lncx/lncz) (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.
Model Constitutive equation(∂t/E+1/η)F=(∂t)u; the differential operator ∂t≡∂/∂tViscoelastic Operator Corresponds to the Modulus of ElasticityE1E→1E∂t+1/τ∂tThe modified Creep Modelu(t)h0=(bχEP0L0)1/(χ+1)e-t/τ×1F1[2+χ1+χ;1;tτ],1F1[2+χ/1+χ;1;t/τ] is the confluent         hypergeometric function (called Kummer’s equation) written as:  1F1[c;d;x]=1+cdx+c(c+1)x2d(d+1)2!+c(c+1)(c+2)x3d(d+1)(d+2)3!+⋯, or 1F1[c;d;x]=∑n=0∞(c)n(d)nxnn!Table 2
Kelvin-Voigt viscoelastic medium.
Model Constitutive equationε=σEe+η∂t; the differential operator ∂t≡∂/∂tViscoelastic Operator that Corresponds to the Modulus of ElasticityE1E→1Ee11+τ∂tThe modified Creep Modelu(t)h0=1Γ((2+χ)/(1+χ))(bχτEeF0L0)1/(χ+1)τ(2+χ)/2(1+χ)×t-χ/2(1+χ)e-t/2τM-χ/2(1+χ),1/2(1+χ)[t/τ],  where Γ is the gamma function, Mμ,ν(x) is the Whittaker’s hypergeometric function and equals to:  Mμ,ν(x)=x1/2+νe-x/2×1F1[12+ν-μ;2ν+1;x] and 1F1[c;d;x] is the Kummer’s confluent  hypergeometric function which could be expressed as 1F1[c;d;x]=∑n=0∞(c)n(d)nxnn!, or:  1F1[c;d;x]=1+cdx+c(c+1)x2d(d+1)2!+c(c+1)(c+2)x3d(d+1)(d+2)3!+⋯The modified Thermal Creep Modelu(t)h0=(bχτEeP0L0)1/(χ+1)t1/(1+χ)Γ((2+χ)/(1+χ))1F1[11+χ;2+χ1+χ;-tτ]×exp(QR(1T-1T0))Table 3
Standard Linear Solid viscoelastic medium.
Model Constitutive equationε=(1Ev1+1Ev211+τ∂t)σ; the differential operator ∂t≡∂/∂tViscoelastic Operator that Corresponds to the Modulus of ElasticityE1E→1Ev1+1Ev211+τ∂tThe modified Creep Modelu(t)h0=(bχEv1P0L0)1/(χ+1)Φ2(-1χ+1,1χ+1,1,-Ev1+Ev2η2t,-tτ)The modified Thermal Creep Modelu(t)h0=(bχEv1P0L0)1/(χ+1)Φ2(-1χ+1,1χ+1,1,-Ev1+Ev2η2t,-tτ)×exp(QR(1T-1T0)),where Φ2 is the confluent hypergeometric function of two variables (Appell function)                which is defined by:       Φ2(α,β,γ,x,y)=∑m=0∞∑n=0∞(α)m(β)n(γ)m+nxmynm!n!, and (x)n=(x+n-1)!(x-1)!Table 4
Jeffrey viscoelastic medium.
Model Constitutive equationε=(1η1∂t+1Ev211+τ2∂t)σ; the differential operator ∂t≡∂/∂tViscoelastic Operator that Corresponds to the Modulus of ElasticityE1E→1η1∂t+1Ev211+τ2∂tThe modified Creep Modelu(t)h0=(bχη1η2P0L0(η1+η2))1/(χ+1)t1/(χ+1)Γ((χ+2)/(χ+1))×Φ2(-1χ+1,1χ+1;χ+2χ+1,-tτ2,-Ev2η1+η2t),     where Φ2 is the confluent hypergeometric function of two variables (Appell function)                which is defined by:          Φ2(α,β,γ,x,y)=∑m=0∞∑n=0∞(α)m(β)n(γ)m+nxmynm!n!      where (a)k is Pochhammer's notation, that is, (a)0=1, (a)1=a, (a)k=(a)k-1(a+k-1)


## 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 Tables1–4. The models present an approximate closed form solution for the approach, (u(t)/h0), of the fractal surface as a function of the applied load, (P0/L0). Three sets of numerical parameters are required to evaluate such models as follows.(1)
Those related to the material properties, that is,Ev'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,cx, cz, and s that characterize the Cantor structure of the rough surface and relate through the fractal dimensions D and Dc, where Dc is the fractal dimension of the cantor set and (Dc=0.631) [9, 10].(3)
Those extracted from the experimental results [60], that is, h0 which corresponding to twice the rms height. Another experimental result showed that the fractal dimension D of a ground stainless steel surfaces is D=1.5 [16]. Values for P0 and L0 are not assumed, instead values for (σ=P0/L0) (for unit depth) are assumed.It is important to mention here that the surface parameters like the fractal dimensionD along with the parameters L0 and h0 can be determined experimentally from the surface of material specimen. The fractal dimension D 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 L0 corresponds to the profile length (usually taken as autocorrelation length), and so the length L0 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, (u/h0), versus the nondimensional time record, (t/τ), are shown in  Tables 1–4 and Figures 5–7.Non-dimensional time-strain curves for different stresses and fractal dimensions.
(a)
Maxwell medium(b)
Kelvin-Voigt medium(c)
Standard linear solid medium(d)
Jeffrey mediumExperimental results and analytical isochronous stress-strain creep curves for constant room temperature and various (t/τ).
(a)
Maxwell medium,D=1.24(b)
Maxwell medium,D=1.5(c)
Standard Linear Solid medium,D=1.4(d)
Standard Linear Solid medium,D=1.5(e)
Jeffrey medium,D=1.4(f)
Jeffrey medium,D=1.5(g)
Kelvin-Voigt medium,D=1.4(h)
Kelvin-Voigt medium,D=1.5(i)
Maxwell medium(j)
Kelvin-Voigt  mediumExperimental results and Analytical Isochronous Stress-Strain creep curves for constant (t/τ) and various bulk temperatures.
(a)
Kelvin-Voigt Medium,D=1.4, t/τ≈44, ΔT=0, ΔT=100, and ΔT=200(b)
Kelvin-Voigt Medium,D=1.5, t/τ≈42, ΔT=0, ΔT=100, and ΔT=200(c)
Standard Linear Solid,D=1.4, t/τ≈1680, ΔT=0, ΔT=100, and ΔT=200(d)
Standard Linear Solid,D=1.5, t/τ≈1100, ΔT=0, ΔT=100, and ΔT=200


### 6.1. Maxwell Model

For the Maxwell model to be representative of asolidrather than a fluid,the modulus η (MPa·sec) must be large, comparable in magnitude with the modulus Ev (MPa), that is, (τ≡η/Ev) is assumed to be of the order 1 (units of time; sec) [18]. In the confluent hypergeometric equation in  Table 1, the value of n should be large enough for the result to converge to a stable value; numerically it is assumed that (n=550). The parameter (cx=1.5) is held fixed, and two different values of the parameter {cz=1.155,cz=1.54} (s=2) are used to yield two different fractal dimensions of the Cantor structure, D=1.5 and D=1.24, respectively. The depth h0 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 (t/τ), called isochronous stress-strain curves, in comparison with the experimental results are depicted in Figures 5–7 [18–24]. For (D=1.24) the nondimensional time duration, (t/τ), 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 (D=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 D 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 modulusEv is selected to be 130 GPa and the retardation time, (τ≡η/Ev), is selected to be of the order “1” [61]. The value of n is taken large enough for the result to converge to an accepted accuracy; it is assumed that (n=550). The Cantor structure, shown in Figures 2 and 3, is built from the middle-third Cantor set where the parameter (cx=1.5) is held fixed, giving a Cantor set dimension Dc=0.63093 [9, 10]. Two different values of the parameter {cz=1.155,cz=1.29} with s=2 are used to yield two different fractal dimensions of the Cantor structure, D=1.5 and D=1.4, respectively. It was pointed out [14] that, only for (cz≤2) the profile of the surface of the contacting body is fractal. The fractal dimension, D, of a ground stainless steel surfaces is D=1.5 [16, 43]. The value of h0, which corresponds to twice r.m.s height for the ground surface, is taken as 6.6 μm [60].Figures5–7 show a set of numerical creep data obtained applying the equations presented in  Table 2. In these figures, the strain (u/h0) is plotted versus the nondimensional time record of (t/τ) 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 (t/τ) 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 (D=1.4) the nondimensional time duration, (t/τ), 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 D=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 (t/τ) exceeds 45 [19, 23].


### 6.3. Standard Linear Solid (SLS) Model

For SLS medium, as an illustrative example, the elastic moduliEv1 and Ev2 are selected to be equal, they are given the value of 130 GPa, and the Newtonian viscosity, η2, 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,Ev1, Ev2, η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,cx, cz and s that characterize the Cantor structure of the rough surface and relate through the fractal dimensions D and Dc, where Dc is the fractal dimension of the cantor set and (Dc=0.631) [9, 10].The parameter (cx=1.5) is held fixed, and two different values of the parameter cz=1.15 and cz=1.29 with s=2 are used to yield two different fractal dimensions of the Cantor structure, (D1=1.5) and (D2=1.4), respectively (cf. (2.13)).(3)
Those extracted from the experimental results [60], that is, h0=6.6μm of the ground surface which corresponding to twice the rms height. Another experimental result showed that the fractal dimension D of a ground stainless steel surfaces is (D=1.5) [16]; values for P0 and L0 are not assumed, instead values for (σ=P0/L0), assuming unit depth, are assumed.Figures6 and 7 present the isochronous creep curves for different nondimensional time durations (t/τ) 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 (t/τ) 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 (cx=1.5) is held fixed, and two different values of the parameter cz, that is, {cz=1.155,cz=1.29} yields two different dimensions of the Cantor structure, that is, D=1.5 and D=1.4, respectively. The double r.m.s. value of the heights of roughness h0 is taken as 9.7 μm [60].The material parameters{η1,η2,Ev2} 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 η1 and η2 (MPa·sec) must be large, comparable in magnitude with the modulus Ev2 [11]. In this work, and as an empirical example, the retardation time τ is assumed to be of the order “1”, and Ev=210 GPa [61]. The activation energy Q 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.Figures6 and 7 show the isochronous curves for creep data for different nondimensional time duration, (t/τ) 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.


## 6.1. Maxwell Model

For the Maxwell model to be representative of asolidrather than a fluid,the modulus η (MPa·sec) must be large, comparable in magnitude with the modulus Ev (MPa), that is, (τ≡η/Ev) is assumed to be of the order 1 (units of time; sec) [18]. In the confluent hypergeometric equation in  Table 1, the value of n should be large enough for the result to converge to a stable value; numerically it is assumed that (n=550). The parameter (cx=1.5) is held fixed, and two different values of the parameter {cz=1.155,cz=1.54} (s=2) are used to yield two different fractal dimensions of the Cantor structure, D=1.5 and D=1.24, respectively. The depth h0 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 (t/τ), called isochronous stress-strain curves, in comparison with the experimental results are depicted in Figures 5–7 [18–24]. For (D=1.24) the nondimensional time duration, (t/τ), 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 (D=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 D 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 modulusEv is selected to be 130 GPa and the retardation time, (τ≡η/Ev), is selected to be of the order “1” [61]. The value of n is taken large enough for the result to converge to an accepted accuracy; it is assumed that (n=550). The Cantor structure, shown in Figures 2 and 3, is built from the middle-third Cantor set where the parameter (cx=1.5) is held fixed, giving a Cantor set dimension Dc=0.63093 [9, 10]. Two different values of the parameter {cz=1.155,cz=1.29} with s=2 are used to yield two different fractal dimensions of the Cantor structure, D=1.5 and D=1.4, respectively. It was pointed out [14] that, only for (cz≤2) the profile of the surface of the contacting body is fractal. The fractal dimension, D, of a ground stainless steel surfaces is D=1.5 [16, 43]. The value of h0, which corresponds to twice r.m.s height for the ground surface, is taken as 6.6 μm [60].Figures5–7 show a set of numerical creep data obtained applying the equations presented in  Table 2. In these figures, the strain (u/h0) is plotted versus the nondimensional time record of (t/τ) 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 (t/τ) 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 (D=1.4) the nondimensional time duration, (t/τ), 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 D=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 (t/τ) exceeds 45 [19, 23].


## 6.3. Standard Linear Solid (SLS) Model

For SLS medium, as an illustrative example, the elastic moduliEv1 and Ev2 are selected to be equal, they are given the value of 130 GPa, and the Newtonian viscosity, η2, 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,Ev1, Ev2, η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,cx, cz and s that characterize the Cantor structure of the rough surface and relate through the fractal dimensions D and Dc, where Dc is the fractal dimension of the cantor set and (Dc=0.631) [9, 10].The parameter (cx=1.5) is held fixed, and two different values of the parameter cz=1.15 and cz=1.29 with s=2 are used to yield two different fractal dimensions of the Cantor structure, (D1=1.5) and (D2=1.4), respectively (cf. (2.13)).(3)
Those extracted from the experimental results [60], that is, h0=6.6μm of the ground surface which corresponding to twice the rms height. Another experimental result showed that the fractal dimension D of a ground stainless steel surfaces is (D=1.5) [16]; values for P0 and L0 are not assumed, instead values for (σ=P0/L0), assuming unit depth, are assumed.Figures6 and 7 present the isochronous creep curves for different nondimensional time durations (t/τ) 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 (t/τ) 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 (cx=1.5) is held fixed, and two different values of the parameter cz, that is, {cz=1.155,cz=1.29} yields two different dimensions of the Cantor structure, that is, D=1.5 and D=1.4, respectively. The double r.m.s. value of the heights of roughness h0 is taken as 9.7 μm [60].The material parameters{η1,η2,Ev2} 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 η1 and η2 (MPa·sec) must be large, comparable in magnitude with the modulus Ev2 [11]. In this work, and as an empirical example, the retardation time τ is assumed to be of the order “1”, and Ev=210 GPa [61]. The activation energy Q 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.Figures6 and 7 show the isochronous curves for creep data for different nondimensional time duration, (t/τ) 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 dimensionD 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.



---

*Source: 691270-2011-03-22.xml*