# Viscoelastic Contact Simulation under Harmonic Cyclic Load.

1. Introduction

Important engineering applications involving products like automotive belts and tires, seals, or biomedical devices require accurate prediction of tribological processes between viscoelastic materials such as elastomers or rubbers. Considering that a closed-form description of the viscoelastic contact is difficult to achieve due to complexity of the emerging equations, numerical simulation presents itself as a worthy substitute, capable of assisting the design of tribologically competent products using viscoelastic materials.

Recent developments in numerical resolution of elastic contact problems encouraged a new approach to the problem of viscoelastic contact. Kozhevnikov et al. [9] advanced a new algorithm for the indentation of a viscoelastic half-space based on the Matrix Inversion Method (MIM). Chen et al. [10] developed a new semi-analytical method (SAM) for contact modeling of polymer-based materials with complicated properties and surface topography. These authors [11] studied the multi-indentation of a viscoelastic half-space by rigid bodies using a two-scale iterative method (TIM).

From the point of view of algorithmic complexity, the semi-analytical method (SAM) employed herein, derived from the boundary element method (BEM), has significant computational advantages over the finite element method (FEM), as it requires a 2D spatial discretization only (i.e., the meshing of the potential contact surface), whereas a FEM simulation entails the 3D meshing of the entire contacting bodies. According to this review [16], the computational efficiency of SAM greatly exceeds that of FEM; for example, a 3D SAM contact simulation can be conducted with a computational effort comparable to that of a 2D finite element analysis. In this paper, the algorithmic computational efficiency is optimized by employing state-of-the-art numerical techniques: (1) the conjugate gradient method (CGM), with its superlinear rate of convergence, is used for the resolution of the emerging linear system of equations, whereas (2) the Discrete Convolution Fast Fourier Transform (DCFFT) technique [17] is engaged in the rapid computation of discrete convolution products.

2. Viscoelastic Constitutive Law and Associated Contact Solutions

In the framework of linear theory of viscoelasticity [18], the material exhibits a linear stress-strain relationship; that is, an increase in stress by a constant factor leads to an increase in the strain response by the same factor. The response functions to excitations conveyed by Heaviside step functions are referred to as the material functions of the viscoelastic body, namely, the creep compliance and the relaxation modulus, which are both functions of time t. The creep compliance function [phi](f) describes the viscoelastic strain response to a unit step change in stress, and the relaxation modulus [psi](t), conversely, describes the stress response to a unit step change in strain. With these functions, the linear viscoelastic response to various sequences of stress or strain is assessed, according to the Boltzmann hereditary integral, by either of the two Volterra integral equations:

[mathematical expression not reproducible] (1)

where s and e are the tensors of deviatoric stress and deviatoric strain, respectively. Consequently, (1) can be regarded as the superposition of a series of loading histories consisting in infinitesimal changes in strain or stress, respectively, applied separately in a window of observation [0, t], chosen so that initially (i.e., at t = 0) the viscoelastic body was undisturbed.

Analog mechanical models, constructed from linear springs and dashpots, arranged in series or in parallel, are convenient tools to model the linear viscoelastic behavior under uniaxial loading. The combination rules for these basic units state that creep compliances add in series and relaxation moduli in parallel.

The ideal spring, also referred to as the Hooke model or the ideal solid, is the elastic element in which the force is proportional to the extension. By identifying force with stress and elongation with strain and according to Hooke's law, [sigma](t) = E[epsilon](t), with E being the pertinent (i.e., longitudinal or shear) elasticity modulus. The dashpot, also referred to as the Newton model or the perfect liquid, is the viscous element in which the force is proportional to the rate of extension. According to the Newton equation, [mathematical expression not reproducible] is the rate of strain and [eta] is the viscosity coefficient. Both Hooke and Newton models represent limiting cases of viscoelastic bodies.

A branch constituted by a spring in parallel with a dashpot is known as the Kelvin-Voigt model, whereas a branch constituted by a spring in series with a dashpot is known as the Maxwell model. The differential equation for the Maxwell model can be expressed as [18] [mathematical expression not reproducible], under the assumption that [epsilon](t) = [sigma](t) = 0 for t < [t.sub.0]. The creep compliance function for the viscoelastic half-space described by a Maxwell model consisting in a spring of shear modulus G in series with a dashpot of viscosity [eta] results as [19] [PHI] (t) = 1/(2G)(1 + t/[tau]), where r denotes the relaxation time: [tau] = [eta]/G.

Ting's formalism [6, 7] can be employed to describe the indentation of the Maxwell viscoelastic half-space by a rigid spherical punch of radius R in a step loading W(t) = [W.sub.0]H(t), where H(t) denotes the Heaviside step function and [W.sub.0] is fixed but otherwise arbitrarily chosen. The equation for the pressure distribution achieved at time t after the first contact, at the radial coordinate r, results as

[mathematical expression not reproducible] (2)

where a(t) denotes the contact radius at time t:

a(t) = ([3R[W.sub.0][PHI](t)/8).sup.1/3], (3)

and Re() denotes the real part of its complex argument. This partially analytical solution is a more computationally friendly form of the corresponding equation derived in [20].

The Maxwell model accounts well for relaxation but handles badly both creep (model creeps without bound at constant rate; therefore it is also referred to as the Maxwell fluid) and recovery (only elastic deformation is recovered, and this is done instantaneously). The Kelvin-Voigt model, on the other hand, handles creep and recovery fairly well but does not account for relaxation. Moreover, the latter model exhibits no instantaneous elastic response; consequently, the elasticity modulus is formally infinite, and contact pressure results to be infinite at the contact boundary in the beginning of the loading process, as demonstrated by Ting's model [6, 7] implementation reported in [19]. Therefore, it can be asserted that the assumptions of the Kelvin model make it inappropriate for contact analyses.

The Maxwell and Kelvin models are adequate for qualitative or conceptual analyses, but quantitative representation of the behavior of real materials requires an increase in the number of parameters. The generalized Maxwell model is composed of a number of Maxwell models and an isolated spring in parallel, whereas the generalized Kelvin model consists in a number of Kelvin units plus an isolated spring in series. The Standard Linear Solid Model, also known as the Zener model, can be represented as a spring of shear modulus G in series with a Kelvin model of parameters [G.sub.K] and [[eta].sub.K] or as a spring in parallel with a Maxwell model. By adopting the former representation, the differential equation for the Zener model results as [18]

[mathematical expression not reproducible] (4)

yielding the following creep compliance function:

[mathematical expression not reproducible] (5)

The Zener model exhibits instantaneous elastic strain when stress is instantly applied; if the stress is held constant, the strain creeps towards a limit, whereas, under constant strain, the stress relaxes towards a limit. Moreover, when stress is removed, instantaneous elastic recovery occurs, followed by gradual recovery towards vanishing strain. Using Ting's formalism [6, 7], the pressure distribution in the step loading of a Zener half-space by a rigid spherical punch of radius R results as

[mathematical expression not reproducible] (6)

with a(t) being the contact radius given by (3).

These basic models, having only one relaxation time, are capable of providing only qualitative description of viscoelastic behavior, whereas precise quantitative assessments require more parameters, related to the naturally occurring spectrum of relaxation times of the real viscoelastic material. Such a goal can be accomplished by using a complex model such as the generalized Wiechert model, which consists in several Maxwell units and a free spring, connected in parallel. The shear relaxation modulus function of the Wiechert model can be expressed as [18]

[PSI](t) = [g.sub.[infinity]] + [n.summation over (i=1)] [g.sub.i] exp (-t/[[tau].sub.i]), (7)

where [g.sub.[infinity]] is the long-term modulus (longitudinal or shear) once the material is totally relaxed and [[tau].sub.i] and [g.sub.i], with [[tau].sub.i] = [[eta].sub.i]/[g.sub.i], are the relaxation time and the spring stiffness of each Maxwell subunit. The naturally occurring spectrum of relaxation times of a viscoelastic material can be described by including as many exponential terms as needed. Relation (7) is also referred to as a Prony series. The Prony series of a viscoelastic material is usually obtained by a one-dimensional relaxation test, in which the viscoelastic material is subjected to a sudden strain that is kept constant, while measuring the stress response over time. The initial stress is related to the purely elastic response of the material. Later on, the stress relaxes due to the viscous effects in the viscoelastic material. Mathematical description employing the Prony series can be achieved by fitting the experimental data to (7) by adjusting the model parameters [g.sub.[infinity]], [g.sub.i], and [[tau].sub.i].

The constitutive law of the real viscoelastic material considered in this paper is that of the polymethyl methacrylate (PMMA), a thermoplastic polymer whose mechanical properties were studied extensively by Ramesh Kumar and Narasimhan [21]. These authors obtained experimentally the PMMA relaxation modulus data under uniaxial compression in a window of observation of 1000 s. Based on their results, the two-term Prony series of PMMA results as [10]

[mathematical expression not reproducible] (8)

which is the modulus relaxation function of the material (expressed in terms of longitudinal modulus), depicted in Figure 1.

In the time domain, the creep compliance and relaxation modulus are not reciprocal (like in the purely elastic case); that is, [psi](t)[PHI](f) [not equal to] 1. However, in the Laplace transform domain, the following relation applies [18] to their transforms:

[bar.[PSI]] (s) [bar.[PHI]] (s) = 1/[s.sup.2], (9)

where s is the variable in the Laplace transform domain. The latter equation can be used to derive the creep compliance function of PMMA by computing first its Laplace transform [bar.[PHI]] (s) from (9) and by subsequently applying inverse Laplace transform to obtain [PHI](T) in the time domain, leading to

[mathematical expression not reproducible] (10)

The latter equations express the PMMA creep compliance in terms of shear modulus by using a Poisson's ratio v = 0.38 [21]. The semi-analytical solution of the contact problem involving linear viscoelastic materials described by complex rheological models like the Prony series is discussed in the following sections.

3. Contact Model

The contact model employed in this paper is based on the general contact model developed by Johnson [20] for the elastic domain. The contact equations, as well as the imposed assumptions and limitations, are repeated here for clarity, and the newly established dependencies, related to the viscoelastic constitutive law described in the previous section, are then discussed in detail.

The set of equations and inequalities governing the contact process are written in a Cartesian coordinate system with [x.sub.1] and [x.sub.2] axes laying in the common plane of contact (i.e., the plane passing through the initial point of contact, which separates best the two contacting surfaces). The two contacting solids are subjected to a normal force aligned with x3 axis, compressing the two bounding surfaces together. As opposed to a time-independent purely elastic contact process, in which the final state depends only on loading level, the viscoelastic contact state depends on time, as well as on the loading history, due to the memory effect of the viscoelastic materials, thus adding a third time parameter to the elastic contact model.

The static force equilibrium relates the normal force W to the pressure distribution p at any time in the observation window [0, [t.sub.0]]. To keep the number of independent parameters to a minimum, the contact is assumed to be frictionless, meaning that shear tractions cannot be sustained at the contact interface. Moreover, the problem is considered as quasistatic, meaning that the inertia forces due to deformation are negligible:

[mathematical expression not reproducible] (11)

The equations of the surface of separation between the two contacting bodies yield the geometrical conditions of deformation in the normal direction:

[mathematical expression not reproducible] (12)

where hi is the gap between the undeformed (i.e., initial, at time t = 0) surfaces, h is the gap between the deformed surfaces, u is the relative normal displacement, and to is the rigid-body approach.

The contact model is completed with the boundary conditions and constraints, also referred to as the complementarity conditions in the literature of the elastic contact. The latter equations are required as the contact area is not known a priori and consequently must be found in an iterative manner by a trial-and-error approach. The gap h between the deformed contacting surfaces vanishes on the contact area, as no interpenetration of the contacting solids is allowed in the frame of elasticity. On the other hand, the gap must be positive outside the contact area, where there is clearance between the contacting bodies. In the same manner, pressure is positive on the contact area and vanishes outside the contact area. These boundary conditions and constraints that must be satisfied simultaneously can be expressed as

[mathematical expression not reproducible] (13)

[mathematical expression not reproducible] (14)

The difficulty in solving the contact model (11)-(14) stems from the fact that neither the contact area nor the pressure distribution is known in advance. An iterative approach is therefore needed, involving a trial-and-error approach, in which a contact region is assumed, and the pressure distribution is then computed based on this assumption. If all constraints in the contact model are verified by the obtained solution, the contact problem solution is achieved. This solution is unique based on the theorem of uniqueness of solution of the elastostatic problem. Otherwise, the contact area is adjusted and a new pressure distribution is computed with the new guess. This iterative approach requires that the response of the contacting material, that is, the displacement induced by the surface tractions, is computed for arbitrary contact area and pressure distribution. The latter computation can only be achieved semi-analytically, and therefore a problem discretization is imposed to perform the numerical analysis of the contact process.

The contact model reviewed herein was also used extensively in the simulation of contact scenarios involving history-dependent processes like plasticity [24], wear [25], or friction [26] by adding an external loop in which the load was applied incrementally. In the latter contact scenarios, the time parameter does not need to be considered explicitly as long as the history of the contact process is properly simulated (i.e., the load is applied with sufficiently small increments). The present work attempts to link the contact model to the theory of viscoelastic behavior, in which the material properties depend explicitly on time. Manipulation of the existing semianalytical solution in the elastic domain, aiming to achieve a viscoelastic contact algorithm, is detailed in the following sections.

4. Viscoelastic Displacement

A surface distribution of normal tractions, such as the pressure resulting from a mechanical contact process, induces a displacement field whose knowledge is essential in solving the contact problem and in performing stress analysis in the contacting bodies. Although the limiting boundary of a real solid is intrinsically rough, computational contact mechanics employ the half-space assumption, allowing for the use of fundamental solutions (i.e., the Green functions) derived in the theory of linear elasticity for a semi-infinite body bounded by a plane surface. For this approximation to remain valid, the slope of the contact geometry must remain small throughout the contact region. The normal displacement field [u.sup.(e)] generated in a linear elastic and isotropic solid by a distribution of normal tractions p([x.sub.1], [x.sub.2]) is computed by applying the superposition principle to the Green function [G.sup.(e)] ([x.sub.1], [x.sub.2]) for the elastic half-space derived by Boussinesq [27]:

[mathematical expression not reproducible] (15)

where [mathematical expression not reproducible] is the normal displacement induced at a point of coordinates ([x.sub.1], [x.sub.2]) by a unit concentrated force acting in origin along direction of [[??].sub.3]], and v and G are Poisson's ratio and the shear modulus of the elastic half-space.

Lee and Radok [3] obtained the contact radius in the viscoelastic spherical contact problem by applying the hereditary integral operator of type (2) to the Hertz (i.e., purely elastic) solution in which the elastic compliance 1/(2G) was replaced by the viscoelastic creep compliance [PHI](f). This course of action is justified by the classic method for solving the linear viscoelastic problems of stress analysis, which is based on the concept of associated elastic problem [1, 2]. Capitalizing on the fact that basic integral equations for stress analysis in viscoelastic materials reduce in the Laplace transform domain to the type of integral equations describing stresses in elastic materials, it has been shown [1, 2] that a viscoelastic problem has an associated elastic problem, to which the former reduces after removal of time dependency by application of the Laplace transform. Consequently, if the boundary conditions are time-independent, a solution in the frequency domain is identical in form to the corresponding elastic solution. This technique of deriving viscoelastic solutions from their elastic counterpart is also referred to as the correspondence principle. The indentation of a viscoelastic half-space by a rigid indenter cannot generally be solved in this manner, as the contact problem features time-dependent boundary conditions, which impede transfer to Laplace domain. When applying this technique to the contact radius formula in the associated Hertz elastic problem, the viscoelastic solution holds true as long as the contact radius increases monotonically with time, but additional manipulations are required when the time-dependent contact area is an arbitrary function of time, as shown in [6, 7].

In this paper, the same technique of replacing the elastic contact compliance with the viscoelastic creep compliance function is applied to (15):

[mathematical expression not reproducible], (16)

yielding the viscoelastic displacement generated by a known history of pressure p([x'.sub.1], [x'.sub.2],t') in a window of observation [0,t]:

[mathematical expression not reproducible] (17)

Unlike its counterpart expressing the contact radius in a Hertz-type viscoelastic contact, the displacement equation (17) does not require additional manipulations and can be used in conjunction with any history of boundary conditions. Interchanging differentiation and integration in (17) yields

[mathematical expression not reproducible] (18)

In a viscoelastic contact problem, the contact area and the pressure distribution are not known in advance and, moreover, keep changing during the contact process, as the response of the viscoelastic material also changes with time. Consequently, integral (18) must be evaluated for various loading histories, implying integration of arbitrary functions over arbitrary domains. The semianalytical treatment of these equations to attain a computationally friendly form is detailed in the following sections.

5. Problem Discretization

Neither integral (15) for the elastic case nor (18) for the viscoelastic framework can be computed analytically for general contact geometry, loading history, or material properties. Important research efforts were dedicated to obtaining the solution of these integrals in a semi-analytical manner [17, 28]. The principle of the semi-analytical method consists in considering all continuous distributions as piecewise constant functions, uniform within each surface element in a uniformly spaced rectangular mesh established in the common plane of contact. Control points must be chosen for all the elementary cells of the grid (the centers of the cells make good candidates), and all continuous parameters are evaluated in these representing points, resulting in a digitized counterpart for each continuous distribution. This discretization encourages a simplified notation taking as arguments the indexes of the cells rather than the continuous coordinates. For example, p(i,j), with i and j integers, denotes the pressure value computed in the center of the cell (i,j), and p(i, j) = p([x.sup.(i).sub.1], [x.sup.(i).sub.2]), where [x.sup.(i).sub.1] and [x.sup.(j).sub.2] are coordinates of the center of the cell (i, j). Consequently, pressure is assumed to be uniform in any rectangular patch and therefore can be factored outside the integral operator in (15). The integral of the Green function [G.sup.(e)]([x.sub.1], [x.sub.2]), taken over the elementary patch of side lengths [[DELTA].sub.1] and [[DELTA].sub.1] 2 along directions of [mathematical expression not reproducible], respectively, yields the influence coefficient [28] for the elastic displacement [K.sup.(e)]:

[mathematical expression not reproducible] (19)

which expresses the normal displacement induced in the observation cell (i,j) by a uniform pressure of magnitude 1/([[DELTA].sub.1] [[DELTA].sub.2]) Pa acting in the cell (k,l). The closed-form solution of the double integral (19) was derived by Love [29]:

[mathematical expression not reproducible] (20)

Within this framework, the semi-analytical counterpart of (15) results as

[mathematical expression not reproducible] (21)

where [N.sub.1] and [N.sub.2] denote the numbers of grids along directions of [mathematical expression not reproducible], respectively. The discrete double convolution in (21) can be performed for any imposed pressure distribution. Optimum algorithmic efficiency is achieved using the DCFFT algorithm advanced by Liu et al. [17]. The reduction of the order of computation stems from the convolution theorem, which states that the convolution operation reduces to an element-wise product in the Fourier transform domain. The semi-analytical displacement computation using (21) together with the DCFFT technique is now widely used in computational contact mechanics. In this paper, a generalization of this technique to the case of viscoelastic behavior is proposed.

As equations describing the purely elastic model are intrinsically time-independent, spatial discretization is adequate to circumvent the continuous integration in (15). The additional integration over the time span in which the body was loaded in (17) requires an additional temporal mesh capable of simulating the memory effect specific to viscoelastic materials (i.e., the property that the current state depends upon all previous states succeeded from the first loading). This temporal discretization should be chosen so that at t _0 the body was undisturbed, and the time increment [[DELTA].sub.t] should be small enough so that, during each step, the problem parameters can be assumed to be constant. A piecewise constant law is thus imposed along the temporal axis, adding a third parameter to the notation implemented in the purely elastic model. For example, p(i, j,k) is the discrete counterpart of p([x.sub.1], [x.sub.2],t), denoting the pressure in the elementary cell (i, j) in the spatial mesh, achieved after k time increments, where t = k[[DELTA].sub.t], with k = 1... [N.sub.t]. These assumptions regarding the temporal variation of model parameters authorize the substitution of the partial derivative [partial derivative]p([x'.sub.1], [x'.sub.2],t')dt'/[partial derivative]t' in (18) with the finite difference p(i, j, k) - p(i, j, k- l). Discretization of (18) in the time domain yields

[mathematical expression not reproducible] (22)

It should be noted that, in the latter equation, the simplified notation related to problem digitization was used only for the temporal parameter, whereas continuous coordinates are employed for spatial localization. As pressure is uniform within each elementary patch, it can be factored outside the spatial integral operator as in the purely elastic model, allowing for a viscoelastic influence coefficient [K.sub.(v)] defined similarly to its elastic counterpart [K.sup.(e)] in (19):

[mathematical expression not reproducible] (23)

The latter equation employs notation related to discretized parameters in both spatial and temporal dimensions. The influence coefficient (23) expresses the displacement observed after k time steps in the elementary patch (i, j) of the spatial mesh, due to a uniform pressure of 1/([[DELTA].sub.1] [[DELTA].sub.2]) Pa which acted in the patch (l, m) in the nth time step after the reference moment in which the body was undisturbed, with n [less than or equal to] [N.sub.t], k [less than or equal to] [N.sub.t], and n [less than or equal to] k. The semi- analytical counterpart of (18), discrete in both time and space dimensions, can thus be expressed as

[mathematical expression not reproducible] (24)

where 1 = l... [N.sub.1], j = 1... [N.sub.2], and k = 1... [N.sub.t]. This equation clearly shows that the memory effect is considered explicitly in the displacement computation, as pressure distributions in all previous states (i.e., in all previous time increments), together with the current pressure, are needed to evaluate the current displacement. It is noteworthy that the contribution of the historical pressures can be separated from the contribution of the current pressure, which leads to the algorithmic strategy described in the following section.

6. Algorithm Overview

The viscoelastic contact simulation is achieved by constructing a series of elastic contact problems with boundary conditions that match the ones of the viscoelastic contact problem at each new time increment in the temporal discretization. This approach is based on the fact that, provided the compatibility and internal equilibrium equations are satisfied instantaneously, any elastic solution to the contact problem befits an instantaneous viscoelastic solution. The set of equations and inequalities (11)--(15) describe in fact a purely elastic contact process, whereas substitution of (15) with relation (18) accomplishes the algorithm generalization to the viscoelastic constitutive law.

Due to the robust nature of the semi-analytical formulation for the purely elastic contact process, history-dependent processes can also be simulated using this algorithm by applying the load in small increments, assuring that the loading path is properly reproduced. In the case of viscoelastic materials, however, the time parameter appears explicitly, and the contact parameters are expected to vary even when the load level is kept constant.

Essentially, the nodal pressures result as the solution of the system of equations arising from digitization of (12). The latter is essentially non-linear due to the presence of the rigid-body approach [omega] but can be linearized by adopting an estimate for this parameter based on the current pressure. Once initial guess values for pressure and the contact area are adopted, (12) can be considered for every cell in the contact area, and the resulting set of equations, all having vanishing gap h, provide an estimation of [omega]. In other words, the rigid-body approach, the pressure distribution, and the contact area are iterated simultaneously, and the true value of [omega] is found when all contact parameters converge. The advantage of this approach consists in linearization of (12), which can thus be solved using appropriate numerical methods for linear systems of equations. According to this review [30] on the subject of numerical methods applied to contact mechanics, the best candidate is the CGM, providing the fastest rate of residual decrease.

Equation (12) applied to all cells in the computational domain generates discrete equations, but only the ones related to cells in the contact area form the system to be solved. Additional difficulty arises as the contact area is not known in advance, and therefore both pressure distribution and contact area have to be iterated simultaneously. Consequently, the size of the linear system, which is directly linked to the number of cells in the contact area, may also vary during the iterative process. At any iteration, cells with negative pressure are excluded from the contact area, while cells with negative gap (i.e., cells where the contacting bodies are predicted to overlap) are (re)included. These adjustments are requested by the contact complementarity conditions in (13) and (14) and force the computed nodal pressures to comply to the boundary conditions of the contact problem. As depicted in Figure 2, these restrictions are imposed outside the core of the CGM-type iteration. The resizing of the linear system leads to reconsideration of the descent directions and the descent steps previously computed in the CGM residual minimization process.

In order to use the CGM, which is proven to converge only for systems having a symmetric and positive definite matrix, it is not convenient to include the static equilibrium equation (11) in the system. In this manner, the system matrix is in fact the influence coefficients matrix, which is symmetric and positive definite (as a diagonally dominant matrix). It should be noted that the diagonal entries in the influence coefficients matrix (for any time t) are reserved for the coefficients expressing the contribution of the pressure located in each cell to the displacement in the same cell. Moreover, (20) suggests that the influence coefficients decay rapidly with the distance between the excitation (i.e., pressure) and the effect (i.e., displacement); therefore the influence coefficients matrix is diagonally dominant for fine meshes. The static equilibrium equation (11) is imposed during each iteration of the CGM algorithm in an additional correction of the system solution outside the CGM core, as depicted in Figure 2. This correction consists in a modification of the nodal pressures proportional to the ratio of the current load to the imposed load.

The iterative process stops when such a system is found that its solution in pressure features only positive entries and the gap distribution takes only positive values on the entire computational domain. The finding of a set of equations of type (12) whose solution additionally satisfies both boundary constraints and static equilibrium marks the achievement of simultaneously converging solutions for pressure, contact area, and rigid-body approach. For each new time interval, a new set of such contact parameters has to be found, as the viscoelastic displacement is updated at each new time increment to account for the new pressure history.

At any iteration and for each newly considered time interval, the computational domain P should fully enclose the contact area; otherwise, the contact process simulation should be restarted with a larger estimate. As discussed in [28], if the predicted contact area reaches the boundaries of P, spurious pressure concentrators due to end effects compromise the solution precision. This is especially true for rough contact scenarios, in which the specimen roughness sample should be properly handled by extending the geometry information with smooth convex shapes. For a detailed description of the algorithm to solve the frictionless normal contact problem, the reader is referred to [15, 28].

The contact solver is adapted in this paper to allow more general boundary conditions, thus addressing an important part of indentation scenarios, in which the rigid-body approach (i.e., the displacement history), rather than the load (i.e., the surface stress history), is imposed. The contact solver advanced by Polonsky and Keer [28] is load-driven (LD) (i.e., the applied normal force is expected as input) but can be modified for displacement-driven (DD) scenarios (i.e., in which the load is not specified, but the rigid-body approach is imposed). A distinctive feature of the original solver [28] is that the normal approach is not explicitly computed. This can be achieved as the algorithm [28] minimizes the residual r(i, j, k) computed from a modified interference equation of type (12):

r (i, j, k) = hi (i, j, k) + u (i, j, k), (i, j) [member of] P, k = 1... [N.sub.t], (25)

in which the normal approach is not directly considered. As rigid-body approach is independent of spatial localization (but dependent of time), its contribution is intrinsically accounted for when the residual is corrected by its mean value [28]. In the DD formulation, on the other hand, the normal approach is known as input data, so the interference computationshouldbeconductedasin thetheoreticalmodel. The system residual matches the gap h between the deformed surfaces and vanishes on the contact area:

r (i, j, k) = hi (i, j, k) + u (i, j, k) - [omega] (k), (i, j) [member of] P, k = 1... [N.sub.t]. (26)

As a result, the residual correction by its mean value in the original algorithm [28] becomes redundant and should be removed. The lack of normal load as input, on the other hand, compromises the adoption of the guess value for the pressure distribution, as the mean pressure can no longer be computed. In our numerical simulations, no convergence problems occurred with nonvanishing uniform pressure distribution as initial guess. The algorithm sequence performing pressure correction with respect to the imposed load also becomes redundant and should be eliminated. The fact that the correct value of the rigid-body approach is known from the start coerces pressure to converge without any additional algorithm modifications. The speed of convergence for the modified DD version of the contact solver was found to be of the same order of magnitude as its LD counterpart. The flowchart for the generalized algorithm, comprising both common and specific steps pertinent to each type of boundary conditions, is depicted in Figure 2.

It was shown in the previous section that the instantaneous viscoelastic displacement can be derived based on the knowledge of all previous contact states. Capitalizing on this result, the simulation of the viscoelastic contact process can be achieved by solving a series of sequential contact states related to a temporal discretization fine enough so that the memory effect specific to viscoelastic materials is accurately captured. The main algorithm steps are detailed as follows:

(1) Acquire the input of the viscoelastic contact problem: the loading history W(t), t [member of] [0, [t.sub.0]], the initial contact geometry hi([x.sub.1], [x.sub.2]), and the contact compliance, that is, the creep compliance function [PHI](f) of the viscoelastic material.

(2) Adopt the computational domain P and the spatial and temporal discretization parameters: [N.sub.1], [[DELTA].sub.1], [N.sub.2], [[DELTA].sub.2], [N.sub.t], and [[DELTA].sub.t].

(3) Digitize all problems inputs. Obtain hi(i, j), (i, j) [member of] P,W(k), [PHI](k), k = 1... [N.sub.t].

(4) Compute the array of influence coefficients K(i, j, k), (i, j) [member of] P, k = 1 ... [N.sub.t], which can be stored as a three-dimensional array, having two dimensions related to the spatial dimensions and the third dimension for the temporal discretization. In the kth temporal step, with k [less than or equal to] [N.sub.t], only the first k layers of this array are used.

(5) Solve an initial contact state (k = 1) with the contact compliance [PHI] (1) and obtain the related pressure distribution p(i, j, 1), (i, j) [member of] P. In this particular case, the surface displacement is simply the convolution between pressure and the influence coefficients, with no memory effect to account for:

[mathematical expression not reproducible] (27)

(6) Apply a new time increment (i.e., increase k) and solve the kth instantaneous contact state to get p(i, j, k), (i, j) [member of] P. This can be achieved as all historical pressures p(i, j, 1) ... p(i, j, k - 1), (i, j) [member of] P are known at this point. The product between the influence coefficients and the historical pressure distributions can be included in (12) as an initial state (i.e., superimposed to hi). In this manner, the structure of the model (11)-(14) remains unchanged, and the same type of algorithm can be applied to solve each new time step in the viscoelastic contact problem. This technique is similar to the one used in the semi-analytical resolution [24] of elastic-plastic contact problems, in which the residual displacement due to the plastic region is added to the contact geometry, forming a modified initial state.

(7) If the end of the simulation window is reached (i.e., k = [N.sub.t]), stop program execution and export the computed data. Otherwise, go to step 6.

Based on (24), it follows that the computational complexity of the proposed algorithm is of order O([N.sup.2.sub.1][N.sup.2.sub.2][N.sup.2.sub.t]) when direct multi-summation is applied. The algorithmic efficiency can be increased dramatically by implementing the DCFFT technique [17] in computation of viscoelastic displacement, leading to an improved order of operations of 0([N.sup.2.sub.t][N.sub.1][N.sub.2] log ([N.sub.1] [N.sub.2])). Practically, a viscoelastic contact scenario with a 256 x 256 spatial discretization and 200 time steps can be solved on a 3 GHz quad-core processor in a matter of minutes.

In this paper, an additional conditional branch is added to allow for computation of displacement which persists (for a limited or an indefinite time period, according to the imposed rheological model) after contact opening. The semianalytical equation (24) for viscoelastic displacement computation holds true even after contact opening by assuming vanishing pressure distributions for the time steps subsequent to load removal and until the load is resumed. However, an additional algorithm conditional branch should be added, as the contact solver cannot handle vanishing loads directly. The flowchart of the extended algorithm, proposed for the simulation of contact problems involving linear viscoelastic materials, with general boundary conditions and arbitrary loading programs, is presented in Figure 3. In that figure, spatial localization is omitted for brevity and integer superscripts are used to index the referred time increment.

7. Code Validation and Results

The newly proposed algorithm was benchmarked against the partially analytical results derived by Ting [6, 7] for the Maxwell and Zener viscoelastic half-spaces indented by a rigid spherical punch of radius R in a step loading W(t) = [W.sub.0]H(t), where H(t) denotes the Heaviside step function and [W.sub.0] is fixed but otherwise arbitrarily chosen. It should be noted that whereas Ting's framework requires the punch to be axisymmetric and leads to tedious algebraic manipulations due to the five possible algorithm branches when the contact area does not vary monotonically with time, the algorithm proposed in this paper can readily handle irregular contact geometry, generalized boundary conditions, and complex material properties.

The pressure distributions predicted by the newly advanced computer program at several times after the first contact, depicted with discrete symbols in Figures 4 and 5, match well the data resulting from numerical computation of relations (2) and (6), respectively, displayed using continuous lines. Dimensionless pressure distribution and radial coordinate are defined as ratio to Hertz central pressure [p.sub.H] and contact radius [a.sub.H], both corresponding to the maximum loading level [W.sub.0]. The good agreement between the two data sets legitimizes the novel algorithm for viscoelastic contact simulation.

It should be noted that the contact model for the Maxwell material predicts a contact radius that grows indefinitely in indentation under constant load. The solution should, however, be rejected for large strain.

In order to test the novel DD formulation for the linear viscoelastic domain, the simulation of a Zener viscoelastic half-space with [G.sub.K] = G, indented by a rigid sphere of radius R, was conducted following an imposed history of rigid-body approach, as resulting from Ting's model applied to a step loading, that is, described by the equation [omega](t) = [a.sup.2](t)/R, with the contact radius a(t) given by (3). The predicted history of pressure distribution over a [0, 5[tau]] time range, similar to the one depicted in Figure 5, was compared to that generated by the numerical computation of (6). A good agreement between the two data sets was found, legitimizing the new displacement-driven contact model.

8. Dynamic Mechanical Analysis (DMA)

The inducement of a steady-state undulated response in a viscoelastic specimen (e.g., a uniaxial bar) during a steady-state oscillation test is often used [31] as a method for assessing the mechanical properties of the viscoelastic material. The DMA methodology employs the so-called complex modulus of the viscoelastic material, consisting in a real part, that is, the storage modulus related to the elastic energy storage, and an imaginary part, that is, the loss modulus related to the internal energy dissipation. This complex modulus results by applying the Fourier transform to the Boltzmann hereditary integral of type (1). The phase angle of the complex modulus is defined as the hysteresis angle of the strain response falling behind the stress excitation. In a DMA procedure, the response in amplitude and the phase angle of the viscoelastic specimen to various excitation frequencies are employed to compute the complex modulus in the frequency domain. The latter modulus can be further used [32, 33] in the characterization of the contact behavior of the viscoelastic material under harmonic load indentation.

An approximate analytical solution in terms of the viscoelastic material complex modulus for the indentation depth amplitude in a frictionless spherical indentation under sinusoidal load was derived by Huang et al. [34]. The latter solution is, however, based on the tentative viscoelastic contact solution [3] that assumes a monotonically increasing contact area during the loading history. The latter limitation can be considered too severe unless, in a loading program consisting in a sinusoidal load superposed on a carrier step load, the indentation depth amplitude due to the oscillation is negligible compared to that due to the step load. The aforementioned limiting assumptions can be released in the numerical simulation, and the analysis of the hysteresis energy loss under the cyclic load can be conducted with a precision challenged only by the spatial and temporal discretization errors.

The Zener viscoelastic contact performance under harmonic load is first analyzed in this section using the newly proposed semi-analytical model. The indentation loading history is sinusoidal, and the phase is chosen so that at t = 0 the load vanishes together with its first derivative:

W (t) = [W.sub.0]/2 x [1 + sin (2[pi]ft/[tau] + 3[pi]/2)] H (t), (28)

where f(Hz x s) denotes the number of oscillations (cycles) that occur each time interval equal to [tau]. The latter parameter and the load amplitude are fixed but otherwise can be arbitrarily chosen. By varying f (i.e., by varying the excitation frequency with respect to [tau]), the effect of loading frequency on contact parameters can be assessed. Dimensionless time parameter is defined as ratio to [tau].

Further observations can be made based on analysis of contact response in the steady-state regime in an imposed frequency spectrum of the excitation load (with f ranging between 0.1 and 1.5 Hz x s). The numerical simulations suggest that the amplitude of contact radius decays with increasing excitation frequency, while the amplitude of peak pressure increases, as shown in Figure 10. Another interesting feature is that, at the same loading level, the contact is more conformal on the descending branch of the sinusoid, as proven by the pressure profiles depicted in Figure 11. The gap is more apparent at smaller loading frequencies.

The force-displacement curves depicted in Figure 12 prove that the contact process is hysteretic. The out-of-phase response between load and indentation depth is the source of energy damping in the viscoelastic material. The energy loss occurring each cycle is quantifiable by the area enclosed by the hysteresis loop. As this loop is obtained in discrete form, fine temporal discretization is required to obtain its detailed description. The loops presented in Figure 12 were obtained by imposing 100 temporal steps per loading cycle, and the simulation time was extended until a steady-state regime was reached (a maximum of 10 cycles were necessary for the case f = 1.5 Hz x s). The simulation results suggest that, for the investigated spectrum of excitation frequency and for the considered rheological model, the energy loss decays with the increasing of the excitation frequency.

A study of the variation of energy dissipation under various frequencies of the excitation load is conducted for the PMMA material, involving the spherical indentation of a PMMA half-space in a loading program consisting in a sinusoidal load with an amplitude of [W.sub.A], fixed but otherwise arbitrarily chosen, superposed to a step carrier load of [W.sub.0] = 2[W.sub.A]. The oscillation was applied after the contact stabilization under the step load. Considering the creep behavior of the contact condition obtained for the rough contact of a PMMA specimen in [15], 500 s were allowed for the contact to reach a steady state:

W(t) = [[W.sub.0] + [W.sub.A] sin ([[omega].sub.f]t) H (t - 500)] H (t). (29)

By conducting the numerical simulation under an extended frequency spectrum of the excitation load, the stabilized viscoelastic contact behavior and the energy loss per cycle were obtained without any simplifying assumptions. The predicted hysteresis loops are similar to those presented in Figure 12, but the energy loss per cycle, related to the internal viscous friction in the viscoelastic material, was found to possess a maximum at a specific harmonic load frequency, as depicted in Figure 13. The numerically computed areas [delta]H of the hysteresis loops enclosed by the force-displacement curves, normalized by the maximum value, suggest that the energy loss rises to a peak value when [[omega].sub.f] increases from 0 to 0.0031[s.sup.-1] and then decays with further increase of the loading frequency to an asymptotic value. This behavior may be attributed to the use of two relaxation times in the PMMA relaxation curve, as opposed to the single relaxation time in the Zener rheological model. This result proves the ability of the proposed numerical model to simulate and predict the contact performance of viscoelastic materials with complicated mechanical properties, as resulting from actual experimental measurements.

9. Conclusions

A robust algorithm for the resolution of linear viscoelastic contact scenarios is advanced in this paper by generalizing an existing well-known numerical solution for the contact of linear elastic bodies. The method relies on applying the correspondence principle together with the Boltzmann hereditary integral to achieve numerical computation of viscoelastic displacement, based on a temporal discretization in which every new state is assessed considering all previous time increments, thus simulating the memory effect specific to viscoelastic materials. Compared to other existing methods, the new approach can readily handle arbitrary contact geometry, arbitrary creep compliance (which can be inputted in discrete form as resulting from experimental data), and arbitrary loading programs, leading to contact radii described by increasing or decreasing functions of time. Whereas the classical viscoelastic contact solution still requires numerical treatments that can often raise convergence issues, the newly proposed algorithm is based on the conjugate gradient method, whose convergence is guaranteed. The computational efficiency is increased using a well-established method for rapid calculation of convolution products.

The contact solver assessing the contact parameters from the geometrical condition of deformation is enhanced to allow for more general boundary conditions. The solution for the displacement-driven contact is achieved by modifying the algorithm for the load-driven contact scenario. Another model extension deals with computation of post-unloading displacement, allowing for contact openings during the loading program.

The predictions of the newly advanced numerical program are benchmarked against existing results for a linear viscoelastic half-space described by the Maxwell and Zener rheological models, indented by a spherical punch in step loading.

The contact simulations under harmonic cyclic load suggest that a steady-state, out-of-phase harmonic response is attained after a few cycles. For the Zener theoretical model, the amplitude of contact radius, as well as the energy loss per cycle, decays with increasing excitation frequency, while the amplitude of peak pressure increases. The contact is found to be more conformal on the descending branch of the loading profile. The numerical computation of the energy dissipation in polymethyl methacrylate under an extended frequency spectrum of the excitation load suggests that a maximum is attained at a specific frequency. Further increase of the harmonic load frequency leads to an asymptotic value of the area enclosed by the resulting hysteresis loop.

The performed numerical simulations prove the ability of the newly proposed algorithm to assist the dynamic mechanical analysis in derivation of the complex modulus of the viscoelastic material, providing an important theoretical tool for the study of the viscoelastic contact with fewer simplifying assumptions than other existing partially analytical solutions.

Data Availability

The data supporting the PMMA constitutive law are from previously reported studies, which have been cited. The results of the numerical simulations that support the findings of this study are available from the corresponding author upon request.

https://doi.org/10.1155/2018/9432894

Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was partially supported from the project "Integrated Center for Research, Development and Innovation in Advanced Materials, Nanotechnologies, and Distributed Systems for Fabrication and Control" (Contract no. 671/09.04.2015) and Sectoral Operational Program for Increase of the Economic Competitiveness co-funded from the European Regional Development Fund.

References

[1] E. H. Lee, "Stress analysis in visco-elastic bodies," Quarterly of Applied Mathematics, vol. 13, pp. 183-190,1955.

[2] J. R. Radok, "Visco-elastic stress analysis," Quarterly of Applied Mathematics, vol. 15, pp. 198-202,1957

[3] E. H. Lee and J. R. M. Radok, "The contact problem for viscoelastic bodies," Journal of Applied Mechanics, vol. 27, no. 3, pp. 438-444, 1960.

[4] S. C. Hunter, "The Hertz problem for a rigid spherical indenter and a viscoelastic half-space," Journal of the Mechanics and Physics of Solids, vol. 8, pp. 219-234,1960.

[5] W. H. Yang, "The contact problem for viscoelastic bodies," Journal of Applied Mechanics, vol. 33, no. 2, pp. 395-401,1966.

[6] T. C. Ting, "The contact stresses between a rigid indenter and a viscoelastic half-space," Journal of Applied Mechanics, vol. 33, no. 4, p. 845,1966.

[7] T. C. Ting, "Contact problems in the linear theory of viscoelasticity," Journal of Applied Mechanics, vol. 35, no. 2, pp. 248-254, 1968.

[8] J. A. Greenwood, "Contact between an axisymmetric indenter and a viscoelastic half-space," International Journal of Mechanical Sciences, vol. 52, no. 6, pp. 829-835, 2010.

[9] I. F. Kozhevnikov, J. Cesbron, D. Duhamel, H. P. Yin, and F. Anfosso-Ledee, "A new algorithm for computing the indentation of a rigid body of arbitrary shape on a viscoelastic half space," International Journal of Mechanical Sciences, vol. 50, no. 7, pp. 1194-1202, 2008.

[10] W. W. Chen, Q. J. Wang, Z. Huan, and X. Luo, "Semi-analytical viscoelastic contact modeling of polymer-based materials," Journal of Tribology, vol. 133, no. 4, Article ID 041404,10 pages, 2011.

[11] I. F. Kozhevnikov, D. Duhamel, H. P. Yin, and Z.-Q. Feng, "A new algorithm for solving the multi-indentation problem of rigid bodies of arbitrary shapes on a viscoelastic half-space," International Journal of Mechanical Sciences, vol. 52, no. 3, pp. 399-409, 2010.

[12] S. Spinu and D. Gradinaru, "Semi-analytical computation of displacement in linear viscoelastic materials," IOP Conference Series: Materials Science and Engineering, vol. 95, p. 012111, 2015.

[13] S. Spinu and D. Cerlinca, "A numerical solution to the cattaneo-mindlin problem for viscoelastic materials," IOP Conference Series: Materials Science and Engineering, vol. 145, p. 042033, 2016.

[14] S. Spinu, "Numerical Simulation of Viscoelastic Contacts. Part I. Algorithm Overview," Journal of the Balkan Tribological Association, vol. 21, no. 2, pp. 269-280, 2015.

[15] S. Spinu and D. Cerlinca, "Modelling of rough contact between linear viscoelastic materials," Modelling and Simulation in Engineering, vol. 2017, Article ID 2521903, 2017

[16] M. Renouf, F. Massi, N. Fillot, and A. Saulot, "Numerical tribology of a dry contact," Tribology International, vol. 44, no. 7-8, pp. 834-844, 2011.

[17] S. Liu, Q. Wang, and G. Liu, "A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses," Wear, vol. 243, no. 1-2, pp. 101-111, 2000.

[18] S. P. Marques and G. J. Creus, Computational viscoelasticity, SpringerBriefs in Applied Sciences and Technology, Springer, Heidelberg, 2012.

[19] F. C. Ciornei, On the hertzian circular dynamic contact in the linear viscoelastic domain [Ph.D. thesis], University of Suceava, 2004 (Romanian).

[20] K. L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, UK, 1985.

[21] M. V. Ramesh Kumar and R. Narasimhan, "Analysis of spherical indentation of linear viscoelastic materials," Current Science, vol. 87, no. 8, pp. 1088-1095, 2004.

[22] J. J. Kalker and Y. van Randen, "A minimum principle for frictionless elastic contact with application to non-Hertzian half-space contact problems," Journal of Engineering Mathematics, vol. 6, pp. 193-206, 1972.

[23] H. Yu, Z. Li, and Q. Jane Wang, "Viscoelastic-adhesive contact modeling: application to the characterization of the viscoelastic behavior of materials," Mechanics of Materials, vol. 60, pp. 55-65, 2013.

[24] V. Boucly, D. Nelias, and I. Green, "Modeling of the rolling and sliding contact between two asperities," Journal ofTribology, vol. 129, no. 2, pp. 235-245, 2007.

[25] L. Gallego, D. Nelias, and S. Deyber, "A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III," Wear, vol. 268, no. 1-2, pp. 208-222, 2010.

[26] S. Spinu and M. Glovnea, "Numerical Analysis of Fretting Contact between Dissimilar Elastic Materials," Journal of the Balkan Tribological Association, vol. 18, no. 2, pp. 195-206,2012.

[27] J. Boussinesq, Application despotentielsal'etudedel'equilibreet dumouvementdessolideselastiques, Reed. A. Blanchard, Paris, 1969.

[28] I. A. Polonsky and L. M. Keer, "A numerical method for solving rough contact problems based on the multi-level multisummation and conjugate gradient techniques," Wear, vol. 231, no. 2, pp. 206-219, 1999.

[29] A. E. Love, "The stress produced in a semi-infinite solid by pressure on part of the boundary," Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences, vol. 228, no. 659-669, pp. 377-420,1929.

[30] J. Allwood, "Survey and performance assessment of solution methods for elastic rough contact problems," Journal of Tribology, vol. 127, no. 1, pp. 10-23, 2005.

[31] H. F. Brinson and L. C. Brinson, Polymer Engineering Science and Viscoelasticity: An Introduction, Springer, New York, NY, USA, 2008.

[32] A. Le Gal, X. Yang, and M. Kluppel, "Evaluation of sliding friction and contact mechanics of elastomers based on dynamicmechanical analysis," The Journal of Chemical Physics, vol. 123, no. 1, Article ID 014704, 2005.

[33] A. Le Gal and M. Kluppel, "Investigation and modelling of rubber stationary friction on rough surfaces," Journal of Physics: Condensed Matter, vol. 20, no. 1, p. 015007, 2008.

[34] G. Huang, B. Wang, and H. Lu, "Measurements of viscoelastic functions of polymers in the frequency-domain using nanoindentation," Mechanics of Time-Dependent Materials, vol. 8, no. 4, pp. 345-364, 2004.

Sergiu Spinu (iD) (1, 2)

(1) Department of Mechanics and Technologies, Stefan cel Mare University of Suceava, 13th University Street, 720229 Suceava, Romania

(2) Integrated Center for Research, Development and Innovation in Advanced Materials, Nanotechnologies, and Distributed Systems for Fabrication and Control (MANSiD), Stefan cel Mare University of Suceava, Suceava, Romania

Correspondence should be addressed to Sergiu Spinu; sergiu.spinu@fim.usv.ro

Received 29 January 2018; Accepted 16 April 2018; Published 20 May 2018

Caption: Figure 1: Relaxation modulus function of PMMA.

Caption: Figure 2: Generalized algorithm flowchart for the instantaneous (i.e., with k fixed but arbitrary) contact state.

Caption: Figure 3: Algorithm flowchart for the transient viscoelastic contact.

Caption: Figure 4: Comparison of pressure history in the step loading spherical indentation of a Maxwell half-space.

Caption: Figure 5: Comparison of pressure history in the step loading spherical indentation of a Zener half-space.

Caption: Figure 6: Evolution of contact radius in the indentation of a Maxwell half-space under specific loading programs.

Caption: Figure 8: Contact stabilization under harmonic cyclic load: (a) f = 0.1 Hz-s, simulation time t = 40[tau] (i.e., 4 loading cycles); (b) f = 1 Hz-s, simulation time t = 8[tau] (i.e., 8 loading cycles).

Caption: Figure 9: Contact stabilization, 7th and 8th loading cycles, f = 0.1 Hz x s, simulation time t = 8[tau] (i.e., 8 loading cycles).

Caption: Figure 10: Amplitude of peak contact pressure and contact radius versus excitation frequency.

Caption: Figure 11: Pressure profiles corresponding to the same loading level, f = 0.6 Hz x s.