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

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 introducedmodels are valid onlywhen 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 microand nano-electronic devices, circuits and systems.


Introduction
In this paper, Section 1 presents the practical implications and motivations for our work and the basic background materials on fractal geometry from mathematical and physical point of views.Section 2 presents the fractal contact model of the middle-third Cantor structure which is utilized to obtain the new results introduced in this paper.The continuous elastic model 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.

Fractals Background
Surface topography plays a significant role in tribology, that is, in problems of friction, wear, lubrication and contact 1 .Therefore, the problem of analysis of rough surfaces attracts the attention of engineers and applied mathematicians.Historically, the following engineering parameters, statistical in nature, were used for the characterization of surface roughness: 1 the root mean square of the heights, σ, 2 the root mean square of slopes, σ 2 m , and 3 the root mean square of curvatures, σ 2  k .However, it was realized that the topography of engineered surfaces is too complex to be described completely by a few statistical parameters.Thus, it was found that roughness has a multiscale nature and requires sophisticated mathematical techniques for its description.
First attempts to model the distribution of heights of surface asperities utilize the classical random field theory assumed that the functions of surface model are differentiable.In particular, this implies that limiting values for σ 2  m and σ 2 k 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.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 selfaffine 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.

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: Here G is finite or denumerable subset of O.It was proven that there exists a value s 0 such that: The Hausdorff dimension of the set S, denoted by dim H S, is the number s 0 such that 1.2 holds.Unfortunately, the calculation of the Hausdorff dimension of mathematical objects often demands a lot of effort.Even to find some estimations of the dimension, it is necessary to overcome a number of rather complex mathematical difficulties 25 .This issue called for the use of other definitions of dimension which are useful in applied mathematics for the characterization of fractal objects.One such alternative is the box dimension 2 .The analytical calculation of the box dimension is usually easier since the corresponding definition of this dimension involves coverings by spheres of equal radii.
Let E 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 Mathematical Problems in Engineering cover the set S. The box counting dimension or box dimension, denoted by dim B S, can be defined if the following limit exists: It can be proven that dim B S 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 dim H S / dim B S 1/2.However, it can be proven that dim H S ≤ dim B S.
As a simple alternative to the Hausdorff measure, we can introduce the s-measure m s of a set as the following limit: and define the box dimension as the value s D such that m s S has a jump from 0 to ∞ similar to the behavior of m H 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 .

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: The value of D 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 set S it is proven that dim B S 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.

Self-Similarity and Self-Affinity of Surfaces
Let us recall that a one-to-one mapping M 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 set S 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 λ H z.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 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.

Brownian Surfaces and Random Fractals
Fractional Brownian processes are widely used in creating computer-generated surfaces, in particular landscapes.For example, a profile can be constructed as a graph of 1 − D fBm V H x of index H, where x is taken as the time and z is the random variable of the single valued function V H x with the following property: where denotes averaging over the ensemble, and H is the Hurst index.The scaling behavior of the different traces, V H x , is characterized by a particular H which relate the typical change in Δz x , where z x V H x , is the trace of the fBm, and the change in the spatial coordinate Δx by the simple scaling law 36, 39, 40 : It is known that, with probability equal to "1", the following holds 34 : 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:

R δ lim
which also can be expressed as:

R δ lim
where z is the average value of the profile function z x .
Another tool for the characterization of surfaces is the spectral density function G ω which is the Fourier transform of R δ : G ω cos ωδ dω.

1.12
In general, the following is accepted in fractional Brownian motion 2 .
i If the auto-correlation function R δ 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 profile z x has spectral density: 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 m n of the spectral density G ω provide a useful description of the surface roughness: where ω 0 2π/λ 0 is the wave number corresponding to the profile length λ 0 .It is possible to show that m 0 is the variance of heights r.m.s height of the surface, m 1 is the variance of slopes r.m.s slope , and m 3 is the variance of curvatures r.m.s curvature 42 .

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: where h is a bounded H ölder function of order greater than β.The following complex generalization of the W x; p was considered: 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 function C and the Takagi-Hopson function T : 1.17 have been studied in numerous papers 3, 18, 35, 37 .By direct calculations, one may obtain: 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: is often used for fractal characterization of the surface topography 4, 37, 38 .Here n 1 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 n 1 depends on the length L of the sample and is given by p n 1 1/L and the parameter A determines the position of the spectral density along the log G axis.It was stated that both parameters A and D of the function W 1 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 p k x; p p kγ C x; p , with k ∈ Z where Z is the set of all integers which looks similar to the definition of a homogeneous function Thus, the graph of the function C x; p near any point x 0 is repeated in scaling form near all points p k x 0 , 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 socalled parametric-homogeneous PH functions introduced 2, 15

1.22
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 ∼ x 2−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 stressstrain 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

The Fractal Contact Model
In this work, the problem of contact between a nominal flat surface with a half-space is studied.The considered surface is constructed on the basis of a deterministic fractal; the Cantor structure.The contact problem for that surface is analyzed, assuming that the results hold for all problems with surfaces of the same fractal dimension.The Cantor structure is constructed by joining the segments obtained at successive stages of the construction of a Cantor set to one another, Figures 2 and 3, where L 0 corresponds to the nominal surface length, and h 0 is equal to twice the r.m.s height of the roughness.
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, c x and c z , 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/c x times the length of the initial segment where c x > 1.The depth of the recesses measured from the last step at the i 1 th construction step of the fractal surface is 1/c z times less than the depth of the ith step where c z > 1.From this it can be easily shown that the horizontal length and recess depth of the i 1 th step are: 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
It is obvious that only for c z ≤ 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 s i segments, each of length 14 : where the parameter s corresponds to the number of asperities on a repeating segment 16, The fractal middle-third Cantor structure s 2, where E 0 is the initiator step, {E 1 , E 2 , E 3 } are the other generated step of cantor structure, L's are the lengths of the E's steps, h's are the heights of E's steps, and F is the applied load.In Figure 2 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 Figures 2 and 3 can be considered as a certain graph of a step function.It can be seen that, during an iterative step in constructing the surface, scaling in the horizontal direction is produced as: while in the vertical direction, the corresponding fluctuations Δz i in the ith generation can be defined by considering the probability of obtaining the value of:

Mathematical Problems in Engineering
The fluctuation Δz i in the ith generation can be obtained by assuming that Δz i scales as the expected value z i P z i 14 where: and P z i is the probability of obtaining the value z i which is given by: where it is also produced as: 2.9 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: Hence, we have:

2.11
Substituting 2.5 and 2.11 into 1.8 , we get: 12 from which the self-affine fractal dimension for the contour of the Cantor structure is: where 0 < D c < 1 is the fractal dimension of the middle-third Cantor set.

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; c z h 0 in Figures 2 and 3, and their deformation is assumed to be linear elastic 54 .
4 The analysis in this study is based upon the Radok's technique 55 of replacing the elastic constants in the elastic solution by the corresponding integral or differential operators, which appear in the stress-strain relations for linear viscoelastic materials.This approach can be applied to the contact problem provided that the loading program is such that the contact area is increasing throughout 44 .Thus the simplest approach to this problem is to find the compliance in cases where a corresponding solution for purely elastic material is known.
5 It is also assumed that, with reference to Figures 2 and 3, there exist a series of onedimensional elastic bars distributed in a way that the distance from the initiator step E 0 to the generated step E 3 is indicated by h 0 , from E 1 to E 3 is indicated by h 1 , from E 2 to E 3 is indicated by h 2 , and so forth, 18-24 .
Let F 3 be the force required to compress E 3 until E 2 , F 2 be the force require to compress E 3 and E 2 until E 1 , and F 1 be the force required to compress E 3 , E 2 and E 1 until E 0 , and so forth, or in general: 3.1 Also, let F i 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 Δu i 1 equal to the difference between the heights protrusion of ith and i 1 th generations.We will assume that E 2 is deflected by a distance u 2 h 1 , E 1 is deflected by a distance u 1 h 0 , E 0 is deflected a distance u 0 c z h 0 , and so Δu 1 u 0 − u 1 h 0 c z − 1 and Δu 2 h 0 c −1 z c z − 1 .Thus:

3.2
The above-mentioned assumptions are sufficient to determine the dependence of the limit load F on the approach u.We will use the fact that the surfaces were approached by an amount Δu i 1 when the limit remote increases from F i 1 to F i , then, for load effect only:

3.3
In the limit when i → ∞, the following asymptotic behavior is obtained as: where ε u/h 0 is the strain due to the remote applied force and χ ln c x / ln c z .

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 "thermorheologically simple" 56 .For such simple materials, the effect of lowering the temperature is simply to shift the viscoelastic response which is plotted against log time to the right without a change in shape.This is equivalent to increasing the relaxation time τ, without changing the relaxation modulus.
A time-temperature shift factor a T 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 T ref .If the creep time obeys an Arrhenius relation, the shift factor can be shown to be 57 : 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: ε f σ, t, T .

4.2
One way to simplify this function is to separate it into three functions of stress, time and temperature as follow 53, 58 : 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: where B is a constant given by: When combining 4.4 and 4.5 , it could be shown that at the reference temperature T 0 the function f 3 T is equal to unity, and the creep is a function of the stress and time only.

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 stressstrain 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 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.
Creep test is a widely used standard test, wherein a force P 0 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 P 0 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.
Tables 1-4 show, for various viscoelastic Maxwell, Kelvin-Voigt, SLS and Jeffreys' models, the time-dependent constitutive equation, viscoelastic operator corresponds to the modulus of elasticity E, and the modified creep model, where u t /h 0 is the relative approach, P 0 /L 0 is the applied stress per unit depth, T is the bulk temperature, T 0 is the reference temperature, η is the Newtonian viscosity, E v is the elastic modulus, Q is the activation energy, R is the ideal gas constant, and τ ≡ η/E v 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: c x and c z , that is, χ ln c x / ln c z cf.2.13 .It is also clear in Tables that simple constant of proportionality between stress and strain does no longer exist.

Results and Discussion
Various continuous models for the creep approach of a viscoelastic surface, modeled by a fractal Cantor structure, in contact with a perfectly rigid foundation are listed in Tables 1-4.The models present an approximate closed form solution for the approach, u t /h 0 , of the fractal surface as a function of the applied load, P 0 /L 0 .Three sets of numerical parameters are required to evaluate such models as follows.
1 Those related to the material properties, that is, E v '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.
Jeffrey viscoelastic medium Figure 4: 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.
Model Constitutive equation hypergeometric function called Kummer's equation written as: where Γ is the gamma function, M μ , ν x is the Whittaker's hypergeometric function and equals to: x and 1 F 1 c; d; x is the Kummer's confluent hypergeometric function which could be expressed as 1 x n n! , or: 2 Those related to the surface texture, that is, c x , c z , and s that characterize the Cantor structure of the rough surface and relate through the fractal dimensions D and D c , where D c is the fractal dimension of the cantor set and D c 0.631 9, 10 .
3 Those extracted from the experimental results 60 , that is, h 0 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 P 0 and L 0 are not assumed, instead values for σ P 0 /L 0 for unit depth are assumed.
It is important to mention here that the surface parameters like the fractal dimension D along with the parameters L 0 and h 0 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 L 0 corresponds to the profile length usually taken as autocorrelation length , and so the length L 0 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 .

Model Constitutive equation
The modified Thermal Creep Model where Φ 2 is the confluent hypergeometric function of two variables Appell function which is defined by: x m y n m!n! , and Analytical creep data for the steel considered in the form of strain, u/h 0 , versus the nondimensional time record, t/τ , are shown in Tables 1-4 and Figures 5-7.

Maxwell Model
For the Maxwell model to be representative of a solid rather than a fluid, the modulus η MPa • sec must be large, comparable in magnitude with the modulus E v MPa , that is, τ ≡ η/E v 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 c x 1.5 is held fixed, and two different values of the parameter {c z 1.155, c z 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 h 0 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 where Φ 2 is the confluent hypergeometric function of two variables Appell function which is defined by: x m y n m!n! where a k is Pochhammer's notation, that is, a 0 1, a 1 a, a k a k−1 a k − 1 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.

Kelvin-Voigt Model
In this study the modulus E v is selected to be 130 GPa and the retardation time, τ ≡ η/E v , 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 c x 1.5 is held fixed, giving a Cantor set dimension D c 0.63093 9, 10 .Two different values of the parameter {c z 1.155 , c z 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 c z ≤ 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 h 0 , which corresponds to twice r.m.s height for the ground surface, is taken as 6.6 μm 60 .
Figures 5-7 show a set of numerical creep data obtained applying the equations presented in Table 2.In these figures, the strain u/h 0 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 .

Standard Linear Solid (SLS) Model
For SLS medium, as an illustrative example, the elastic moduli E v1 and E v2 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, E v1 , E v2 , η 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.result showed that the fractal dimension D of a ground stainless steel surfaces is D 1.5 16 ; values for P 0 and L 0 are not assumed, instead values for σ P 0 /L 0 , assuming unit depth, are assumed.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 .

Jeffrey Model
To synthesize the required fractal surface, parameter c x 1.5 is held fixed, and two different values of the parameter c z , that is, {c z 1.155, c z 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 h 0 is taken as 9.7 μm 60 .
The material parameters {η 1 , η 2 , E v2 } 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 E v2 11 .In this work, and as an empirical example, the retardation time τ is assumed to be of the order "1", and E v 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.
Figures 6 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.

Conclusions and Future Work
This paper presents recent advancements in the creep contact of various linear viscoelastic models; Maxwell model, thermal and isothermal Kelvin-Voigt model, Jeffrey model, and thermal and isothermal standard linear solid SLS model.Such models were used to study the effect of the normal creep load, the surface quality, and/or the bulk temperature on the creep process of a nominally flat surface in contact with a rigid half space.Since these models are constituted from combinations of dashpots and springs, that is, not real materials, it is not expected to predict precisely the real material behavior for all situations.In fact they fail to predict the tertiary region, which means that quantitative reproduction of real material behavior requires a nonlinear viscoelastic stress-strain constitutive relation.
The effects of temperatures on the behavior of such mathematical models are shown to follow trend similar to the isothermal curve.Such behavior is a consequence of assuming thermorheologically simple materials.Also, because of the periodicity of the Cantor set model, it undergoes the same construction procedure at each hierarchical level producing contact areas that are all of the same size.Therefore, it is doubtful that these models will provide an exact simulation of the deformation of random rough surface.However, these models do admit an analytical solution, where the specific character of a fractal models has little effect on the asymptotic behavior of the process, and the fractal dimension D 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.

H 2 H 5 HFigure 1 :
Figure 1: Common fractals: a the middle-third Cantor set, b the von Koch curve, c the Weierstrass-Mandelbrot function C in the range 0 ≤ 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 ∼ x 2 , 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.

19 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/ log p ≤ dim H graph C ≤ 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: which strictly satisfy the equation b d p k x; p p kd b d 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 b 1 and b 2 with degrees d 1 and d 2, respectively: b 0 x; p x −γ C x; p , b 1 x; p xb 0 x; p , b 2 x; p x 2 b 0 x; p .

Figure 3 :
Figure 3: The fractal middle-third Cantor structure s 3, where E 0 is the initiator step, {E 1 , E 2 } 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 6 :
Figure 6: Experimental results and analytical isochronous stress-strain creep curves for constant room temperature and various t/τ .

2 20 F 16 F
Those related to the surface texture, that is, c x , c z and s that characterize the Cantor structure of the rough surface and relate through the fractal dimensions D and D c , where D c is the fractal dimension of the cantor set and D c 0.631 9, 10 .The parameter c x 1.5 is held fixed, and two different values of the parameter c z 1.15 and c z 1.29 with s 2 are used to yield two different fractal dimensions of the Cantor structure, D 1 1.5 and D 2 1.4 , respectively cf. 2.13 .3 Those extracted from the experimental results 60 , that is, h 0 6.6 μm of the ground surface which corresponding to twice the rms height.Another experimental 0 /L 0 (MPa) Displacement u (μm) a Kelvin-Voigt Medium, D 1.4, t/τ ≈ 44, ΔT 0, ΔT 100, and ΔT 200 /L 0 (MPa) Displacement u (μm) b Kelvin-Voigt Medium, D 1.5, t/τ ≈ 42, ΔT 0, ΔT 100, and ΔT 200 0 (MPa) Displacement u (μm) c Standard Linear Solid, D 1.4, t/τ ≈ 1680, ΔT 0, ΔT 100, and ΔT 200 0 (MPa) Displacement u (μm) Experimental d Standard Linear Solid, D 1.5, t/τ ≈ 1100, ΔT 0, ΔT 100, and ΔT 200

Figure 7 :
Figure 7: Experimental results and Analytical Isochronous Stress-Strain creep curves for constant t/τ and various bulk temperatures.

Figures 6
Figures 6 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 Figures6 and 7is 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 scales22, 24 .
, 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 .

Table 2 :
Kelvin-Voigt viscoelastic medium.Viscoelastic Operator that Corresponds to the Modulus of Elasticity E

Table 3 :
Standard Linear Solid viscoelastic medium.

Table 4 :
Jeffrey viscoelastic medium.Viscoelastic Operator that Corresponds to the Modulus of Elasticity E