Modelling elastic structures with strong nonlinearities with application to stick–slip friction

Robert Szalai


An exact transformation method is introduced that reduces the governing equations of a continuum structure coupled to strong nonlinearities to a low-dimensional equation with memory. The method is general and well suited to problems with isolated discontinuities such as friction and impact at point contact. It is assumed that the structure is composed of two parts: a continuum but linear structure and finitely many discrete but strong nonlinearities acting at various contact points of the elastic structure. The localized nonlinearities include discontinuities, e.g. the Coulomb friction law. Despite the discontinuities in the model, we demonstrate that contact forces are Lipschitz continuous in time at the onset of sticking for certain classes of structures. The general formalism is illustrated for a continuum elastic body coupled to a Coulomb-like friction model.

1. Background

One of the greatest concerns of engineers is modelling friction and impact. These two strong nonlinearities occur in many mechanical structures, e.g. underplatform dampers of turbine blades [1,2], tyre models [3], or in general any jointed structure [4,5]. The most common way of modelling such systems is to take a finite-dimensional approximation assuming that the omitted dynamics has only a small effect on the overall result. Such models can then be analysed using the well-established theory of non-smooth dynamical systems [6,7]. The application of this theory to engineering structures provides a great insight, even though not all phenomena could be experimentally confirmed [8]. Recent results however indicate significant deficiencies that question the predictive power of finite-dimensional non-smooth models. It was shown that when a rigid constraint becomes slightly compliant in a friction-type system small-scale instabilities develop [9]. This means that refining the model by including more degrees of freedom can lead to qualitatively disagreeing results. The solution can also become non-deterministic [10,11] or non-unique in forward time for an otherwise well-specified initial condition. Therefore, a better modelling framework is necessary that either eliminates inconsistency and non-determinism or at least provides a hint about the physical mechanism that causes such behaviour.

The most apparent problem with finite-dimensional non-smooth models of mechanical systems is that they use rigid body dynamics to describe the motion. This includes finite mode approximation of elastic structures, where each mode has a non-zero modal mass [12]. When two contacting elastic bodies slip, and then suddenly stick their contact points will experience a jump in acceleration. In the case of a finite mass at the contact point, the contact force has a jump, too. In reality, however, the mass of the contact point or contact surface is zero, which implies that the contact force must be continuous. For this reason, a standard finite mode description of elastic bodies is qualitatively inaccurate. The continuity of the contact force should be preserved by mechanical models.

In this paper, we propose a formalism that helps better understand and perhaps solve the above problems. We investigate mechanical systems that consist of linear elastic structures coupled at isolated points of contact with strong nonlinearities. This class of problems include mechanisms with Coulomb-like friction models. The dynamics of impact is considered in a companion paper [13]. Our formalism accounts for the zero mass of the contact point without artificially introducing coupling springs as in Melcher et al. [14] to regularize the problem. To achieve such model reduction and still provide an exact description, we introduce memory terms. Within this new framework, the dynamics is described by a low-dimensional delay equation. We show that in our formulation contact forces are Lipschitz or continuous for certain classes of structures when the solution transitions onto the discontinuity surface. The new formalism also leads to well-defined dynamics, because small perturbations of the reduced model do not affect the qualitative features of the dynamics in general.

Time-delayed models have already been in use when modelling friction [15,16]. In these cases, however the delay parameters are fitted to experimental observations. We hope that through our theory these empirical parameters can gain physical meaning.

The paper is organized as follows. In §2, we present our general mechanical model. In §3, we describe our model reduction technique, discuss the convergence of the method and its implications to non-smooth systems. The derivation of the memory term is illustrated through the examples of a pre-tensed string and a cantilever beam. In §5, we present the example of a bowed string. We demonstrate the properties of the transformed equation of motion, in particular its convergence as the number of vibration modes goes to infinity.

2. Mechanical model

The mechanical model of our structure is divided into two parts, a linear elastic body and several discrete non-smooth nonlinearities that are coupled to the continuum structure. First, we describe our assumptions on the linear but infinite-dimensional part of the model, and then we explain how non-smooth nonlinearities are coupled to the system. The description is sufficiently general to describe friction oscillators and impact phenomena. For simplicity, we only focus on a single elastic structure, but our framework is trivially extensible to mechanisms involving multiple linear structures coupled at (strongly) nonlinear joints.

We assume that the displacement of a material point χ of the structure at time t is represented by u(χ,t). We also assume that the motion u(χ,t) can be expressed as a series Embedded Image2.1where ψk(χ) are three-dimensional vector-valued functions depending on the spatial coordinates of the structure only. The generalized coordinates xk can be arranged into a vector Embedded Image to simplify the notation. As in infinite dimensions not all norms are equivalent, we assume that Embedded Image, which turns Embedded Image into a Banach space [17]. Owing to linearity the governing equations can be written as Embedded Image2.2where C and K are the damping and stiffness matrices, respectively, both assumed being multiplied by the inverse mass matrix from the left. The forcing term fe(t) acts as a placeholder for the non-smooth part of the system and will be replaced with specific terms. Equation (2.2) allows for internal resonances. We assume that these resonances are restricted to arbitrarily large but finite-dimensional subspaces of the state space.

When C and K matrices are simultaneously diagonalizable the equation of motion can be written in the form of Embedded Image2.3where Ω=diag(ω1,ω2,…) and D=diag(D1,D2,…). In the unforced case (fe=0), the vector components of x on the left-hand side of equation (2.3) are decoupled, which means that the homogeneous equation can be solved for each xk independently. Therefore, xk(t) and ψk(χ) are called the modes and mode shapes of the system, respectively, with ωk natural frequencies and Dk damping ratios [12]. System (2.3) are called modal equations describing the motion through the modal coordinates x. Our results do not require that the equations of motion assume form (2.3), however, the examples later illustrate that certain calculations are greatly simplified in this case.

In order to take into account the coupling of the contact points to non-smooth nonlinearities, we need to characterize their motion. For simplicity, we assume point contact. Let us denote the motion of the jth contact point at χj along the direction of vector ϕj by Embedded Image2.4as illustrated in figure 1. We call the positions yj(t) and the velocities Embedded Image of the contact points resolved variables. We assume M contact points, thus we define Embedded Image. Substituting (2.1) into (2.4), we obtain the motion of the contact points as a function of the solution of equation (2.2) Embedded Image2.5where Embedded Image2.6Vectors nj are assumed to be linearly independent, spanning an M-dimensional subspace of Embedded Image.

Figure 1.

A linear elastic structure with contact points χj. Each contact point has a three-dimensional motion, however, we project that motion to vectors ϕj to obtain a scalar valued resolved variable yj. If we need to resolve more than one direction of the motion of the contact point χj, we attach multiple labels χj=χj+1=χj+2 to the same point. Then the motion is projected by the linearly independent vectors ϕj,ϕj+1,ϕj+2 to yield the resolved variables yj,yj+1,yj+2. By this definition, we ensure a one-to-one mapping between indices of labels χj and vectors ϕj. (Online version in colour.)

The nonlinearities are incorporated into the model through the forcing term fe(t). We assume that the nonlinearities only depend on the resolved variables. They are also piecewise continuous with isolated discontinuities. Thus, the contact forces acting at contact points χj in the direction ϕj are written as fj(y(t),t). Summing up all the nonlinearities completes our model description by providing the right-hand side of equation (2.3) in the form of Embedded Image2.7

3. Reduction of the mechanical model

We aim to reduce the number of dimensions of our mechanical model (2.7) to an equation that only involves the 2M number of resolved variables all contained in the vector y. To achieve this, we use the Mori–Zwanzig formalism [18] to arrive at a 2M-dimensional first-order delay equation. The solution of the reduced system agrees with the solution of the full system for the resolved variables, whereas the rest of the variables are discarded. Our method can be viewed as a way of producing a Green's function for only parts of the system. In this formalism, a convolution with the resolved variables represents the effect of the eliminated variables on the dynamics of the resolved variables.

To simplify our calculation, we transform equation (2.2) into a first-order form of Embedded Image3.1where Embedded Image3.2We already have a way of obtaining the resolved variables from the generalized coordinates x through a dot product with vectors nj as shown in equation (2.5). In a matrix–vector notation, we write the conversion as Embedded Image3.3To obtain our reduced model, we construct a projection matrix S that acts on the generalized coordinates and has a 2M-dimensional range. In order to do that, we also define a lifting operator in the form of Embedded Image3.4which is a linear mapping from Embedded Image to Embedded Image. Note that W does not depend on the physical system, it can be chosen to suit the reduction procedure. By combining matrices V and W, we obtain our projections S=WV and Q=IS on the condition that wk satisfy wkvl=0 if kl and wkvk=1. This constraint can also be expressed as I=VW (note the order of the two matrices).

For technical reasons, the choice of W needs to be further restricted. Therefore, we assume that the range of W is invariant under R, that is, Embedded Image3.5Equivalently, wj can be constructed as linear combinations of 2M number of eigenvectors of R, because eigenvectors are invariant by definition. If Embedded Image are some eigenvectors of R and Embedded Image is the matrix with Embedded Image in its columns, the R invariant lifting operator becomes Embedded Image. The matrix Embedded Image is invertible if Embedded Image are not in the kernel of V. Also, both V and Embedded Image have maximum rank, hence their product which is a 2M×2M matrix has maximal rank, too.

In the case of modal equations (2.3), the columns of W can be explicitly constructed in the form of Embedded Image3.6Because each of the four blocks of R are diagonal, condition (3.5) holds when the mj vectors are chosen such that they only have M number of non-zero components Embedded Image3.7

Using the results in appendix A and expression (2.7) of the forcing term, we find that equation (3.1) coupled to the contact forces can be transformed into Embedded Image3.8where Embedded Image3.9 Embedded Image3.10 and Embedded Image3.11The integral in (3.8) is understood in the Riemann–Stieltjes sense (see the electronic supplementary material), z(s) is the initial condition of equation (3.1) at time s and eRt is the fundamental matrix of equation (3.1). The matrix A is bounded if condition (3.5) holds, while H(t)z(s) is bounded if the initial condition z(s) is bounded, too. To obtain the memory kernel Lj(τ), one needs to solve first-order system (3.1) for M different initial conditions vM+j. This solution may not be bounded for most mechanical systems, however, the right-hand side of equation (3.8) stays bounded under general conditions. In subsequent sections, we derive an alternative equation that overcomes this issue.

(a) An alternative form of the reduced equations

It was mentioned before that Lj(τ) may not be a bounded function. Therefore, we integrate the convolution in equation (3.8) by parts, so that the integral of Lj(τ) and the derivative of the contact force fj appear in the rewritten equation. However, to be able to integrate Lj(τ) and ensure that the integral does not grow out of bound, we need to decompose Lj(τ) into a constant and an oscillatory term Embedded Image3.12where Embedded Image3.13As we have a formal expression for Lj(τ) in the form of equation (3.11), we can calculate the integral Embedded Image3.14From the simple rule of differentiating an exponential, we find that Embedded Image. Using this formula twice in equation (3.14), we are left with Embedded Image3.15Assuming that Embedded Image is bounded (see §3d) limit (3.13) becomes Embedded Image3.16Note that the first M element of VR−1vM+j is the static displacement of the contact points under a static unit load at contact point χj. Indeed, by expanding equation R(ξ,η)T=(0,nj)T, we get Kξ=nj and η=0, where ξ is the static displacement of the generalized coordinates, hence, nlξ=[VR−1vM+j]l is the static displacement of the lth contact point.

Our definition (3.12) of Embedded Image implies that either Embedded Image exponentially if the system is dissipative or Embedded Image oscillates with zero mean. Therefore, we define the integral Embedded Image3.17The boundedness and smoothness of (3.17) is discussed in §3d. By virtue of (3.17) the forcing terms in (3.8) can be integrated by parts as follows: Embedded Image3.18which has only bounded terms if the derivative (d/dt)[fj(t)] exists and is bounded almost everywhere and condition (3.33) in §3d holds. Using this transformation, the governing equation becomes Embedded Image3.19where Embedded Image.

It is important to note that according to the theory of delay equations [19] if the terms in equation (3.19) are slightly perturbed, the solution of (3.19) changes only slightly under general conditions. This is a clear advantage over the finite-dimensional description where perturbation generally leads to qualitatively disagreeing solutions [9].

(b) Vibrations of a pre-tensed string

In this section, we illustrate the calculation of A and Embedded Image for a pre-tensed string without bending stiffness. We keep the calculation as general as possible so that it applies to systems with a single contact point that can be written in the form of (2.3). The schematic of the string is shown in figure 2a, whose motion is described by the following equation: Embedded Image3.20where c is the wave speed, u(t,ξ) is the deflection of the string, t stands for time and 0≤ξ≤1 is the coordinate along the string. The boundary conditions u(0,t)=0 and u(1,t)=0 express that there is no movement at the two ends of the string. Equation (3.20) indicates damping, which explicitly appears in the mode decomposed system (2.3) with non-zero damping ratios. The string is forced at ξ=ξ by a contact force Fc(t), which is represented by the Dirac delta function in equation (3.20).

Figure 2.

(a) Schematic of a string and (b) a beam. The displacement of material points is described by u(ξ,t), which represents motion in the direction of the arrows. (Online version in colour.)

The first step is to provide a mode decomposition in the form of equation (2.3), so that the vibration of the string is expressed by (2.1), where the mode shapes are the scalar valued Embedded Image. The natural frequencies of the system are ωk=ckπ and we assume uniform damping Embedded Image for all modes, which gives us Embedded Image3.21The vibration at the contact point can be expressed as a linear combination of all the modes u(ξ,t)=nx(t), where Embedded Image. The resolved variables therefore are y1(t)=nx(t) and Embedded Image. We also choose Embedded Image in formula (3.6), thus Embedded Image3.22

Next, we calculate the function eRtv2 that appears in expression (3.11) of L(τ). The equation whose solution we are seeking is Embedded Image, which is expanded into Embedded Image3.23The initial conditions that correspond to z(0)=v2 are z1,k(0)=0 and z2,k(0)=[n]k. The solution for the modes are Embedded Image3.24Without assuming the form of z1,k(t) and evaluating formula (3.11), we find that Embedded Image3.25Eventually substituting (3.24) into (3.25) in the conservative case (Dk=0) [L(t)]2 of our system (3.20) becomes Embedded Image3.26which is a divergent Fourier series, therefore equation (3.8) cannot be used to describe the dynamics. The constant term in equation (3.25) regardless of the damping ratios assumes the form Embedded Image3.27Note that this is Embedded Image times the static displacement of the string under unit load at ξ=ξ. Using the expressions for z1,k(t) and integrating Embedded Image as per definition (3.17), we get Embedded Image3.28Assuming that Dk=D are constant, the right limit of [L1(τ)]2 becomes Embedded Image3.29The detailed calculation of (3.29) can be found in the electronic supplementary material, which also indicates the boundedness of [L1(τ)]2. The graph of [L1(τ)]2 for c=1 is illustrated in figure 3a for both the conservative and the damped case.

Figure 3.

(a) Graph of [L1(τ)]2 for the string equation (3.20) with at ξ=0.4, c=1. The non-smooth dark (blue) line represents the conservative case and the light smooth (red) corresponds to the damped case with Graphic. (b) Graph of the function [L1(τ)]2 when truncating the series (3.28) at 20,40,80,160 terms, denoted by cross symbols, triangles, squares and circles, respectively. (Online version in colour.)

Note that the convolution kernel L1(τ) for Dk=0 is periodic, and therefore the delay that occurs as the effect of nonlinearities is infinite. If damping is introduced L1(τ) decays in time so that the past of the system will have a smaller effect. This is illustrated by the light smooth (red) line in figure 3a. For non-zero damping, as an approximation one can truncate the delay to a finite time-interval. Truncation is a reasonable choice for most practical purposes, but it is not quite clear what are the theoretical implications [20].

It is worth noting that equation (3.20) in the conservative case can be solved using D’Alembert's formula that also leads to a delay-differential equation [21], which is similar to (3.19).

(c) Euler–Bernoulli cantilever beam

We choose the Euler–Bernoulli beam as our second example to illustrate the theory. This model can support waves of infinite speed, therefore its physical validity is questionable. Nevertheless, it is worth investigating how this property of the Euler–Bernoulli beam translates into the properties of the memory kernel. The non-dimensional governing equation and boundary conditions are Embedded Image3.30The natural frequencies of (3.30) are determined by the equation Embedded Image, which can be approximated by Embedded Image for ωk sufficiently large. Therefore, ωk≈(π/2)2. We use this estimate as a starting point to numerically find more accurate ωk values. The mode shapes at the free end of the beam assume the values given by vector n=(2,−2,2,−2,…)T. We choose Embedded Image in formula (3.6), so that W satisfies our assumption (3.5).

The general formulae that were derived in §3b still apply to equation (3.30) with the appropriate ωk, Dk and [n]k values. In particular, we use (3.28) to plot the memory kernel in figure 4. The graph of [L1(τ)]2 in figure 4a shows that the quadratically growing natural frequencies make the function non-smooth in the conservative case. When damping is introduced, the function becomes smooth for τ>0. Figure 4b shows that for 0≤τ≪1 [L1(τ)]2 grows like a power curve τp, 0<p<1. Therefore, [L1(τ)]2 is not differentiable at τ=0.

Figure 4.

(a) The graph of [L1(τ)]2 for the Euler–Bernoulli cantilever beam. The dark and serrated (blue) curves correspond to the conservative case and the light and smooth (red) curves represent the damped system with Graphic. The conservative case illustrates the lack of smoothness of [L1(τ)]2 for τ>0. (b) Illustrates that [L1(τ)]2 is not differentiable at τ=0. (Online version in colour.)

(d) The convergence of Embedded Image

So far we derived the reduced equation of motion (3.19) without any consideration whether intermediate terms are well defined. To make the analysis rigorous, we introduce infinite-dimensional vector spaces Embedded Image3.31that contain the solutions of the second (2.2) and first-order system (3.1), respectively. One can check that in the previous two examples vM+jZ, because their norm is infinite. This also implies that Embedded Image. Even if we know that Embedded Image, the bound of Embedded Image cannot be directly estimated in the straightforward way, because Embedded Image. However, the two examples in the previous sections show that Embedded Image can be bounded, but not necessarily smooth. In case of the string and without damping, Embedded Image is a piecewise-smooth function, while Embedded Image appears to be continuous but non-differentiable for the undamped Euler–Bernoulli beam. If damping is introduced, Embedded Image becomes smooth for τ>0 and discontinuity or non-differentiability occurs only at τ=0 in both examples.

Boundedness and smoothness depends on the eigenvalues of R, which are directly related to the natural frequencies and damping ratios of system (2.2). We assume that all the eigenvalues of R are in the left half of the complex plane including the imaginary axis, that is, Embedded Image3.32This guarantees that eRt is a strongly continuous semigroup with ∥eRt∥≤C0 as shown in appendix B. Therefore, the solutions of (2.2) with bounded initial conditions are also bounded. However, Embedded Image are calculated as integrals of solutions with possibly unbounded initial conditions. If these unbounded solutions decay sufficiently fast, their integral becomes bounded. Appendix C shows that Embedded Image are bounded if there are Embedded Image such that Embedded Image3.33

Smoothness of Embedded Image for τ>0 can be guaranteed if we replace (3.32) and (3.33) with stronger assumptions. The eigenvalues of R now need to be in a sector of the complex left-half plane to guarantee that high-frequency oscillations decay with a rate proportional to their frequency, so that infinite slopes cannot develop for τ>0. Therefore, in formal terms, we assume that there exists a D0>0 such that the eigenvalues of R satisfy Embedded Image3.34The stronger condition that guarantees the smoothness of Embedded Image for τ>0 becomes Embedded Image3.35The proof of this fact can be found in appendix C, which follows standard techniques from the theory of operator semigroups [22].

In the case of modal equations (2.3), assumption (3.34) holds if there is a D0 such that for the damping ratios Embedded Image3.36

Condition (3.33) has a mechanical meaning. When the structure is forced at contact point χj with Embedded Image, γ>0, the velocity response yM+l(t) when scaled back with the exponential growth must be bounded independent of the forcing frequency ω, that is |yM+l(t) eγt|≤Mj,l. For smoothness of Embedded Image, we require that the decaying forcing Embedded Image produces a similarly decaying velocity with |yM+l(t) eδωt|≤Mj,l for 0<δ<D0 independent of ω.

In general, it is not straightforward to check whether (3.33) holds. Let us consider modal equations (2.3) without damping and assume that the natural frequencies scale as ωk=ω0kℓ/2, where ℓ=2,3,4,…. We also assume that Embedded Image, which implies that Embedded Image3.37where ψ(⋅) is the logarithmic derivative of the Euler Gamma function [23]. Function ψ has isolated singularities on the negative real axis, otherwise it is bounded. Therefore, a γ>0 of (3.33) can be chosen such that none of the arguments of ψ goes through these singularities. This means that there is an Mj,l such that Embedded Image. In particular, for the two examples of the string and the beam we have Embedded Image3.38Note that for ℓ<2, sum (3.37) is not uniformly bounded, owing to the λ1−2/ℓ term in the denominator.

4. Non-smooth dynamics

We are now in the position to include the strongly nonlinear contact forces (2.7) into reduced model (3.19) and investigate their effect. For the sake of simplicity, in this section we assume a single contact force Fc(y), so that the governing equation becomes Embedded Image4.1where g(t)=H(t)z(s)+L0(ts)Fc(y(s)). The properties of solutions of (4.1) strongly depend on both L1(τ) and Fc(y). We also assume that L1(τ) is smooth for τ>0 as it is outlined in §3d. In the case when Fc(y) is a smooth function, but has a steep gradient the solution of (4.1) can become discontinuous [24]. As we are focusing on the discontinuities of Fc(y), we assume that the smooth branches of Fc(y) have sufficiently small gradients so that slipping motion is continuous.

Our definition of solution at the discontinuities of Fc(y) is based on a mechanical analogy. If the elastic bodies stick together, there is an algebraic constraint that restricts the trajectories to sticking motion and one must be able to calculate the contact force implicitly from equation (4.1). If these contact forces are admissible by physical law, the bodies will stick, otherwise they will continue slipping.

To formalize this definition, we assume that Fc(y) is discontinuous along a smooth surface defined by h(y)=0, which stands for the algebraic constraint of sticking. We call Embedded Image the switching surface. The physical bound of the contact force can be defined as the two limits of Fc(y) on the two sides of Σ, that is, Embedded Image4.2Without restricting generality, we assume that Embedded Image. Alternatively, Fc and F+c can be defined on the switching surface Σ independently of Fc, when one wants to distinguish between static and dynamic friction.

According to our physical interpretation of the solution, when a trajectory reaches the switching surface Σ the trajectory either crosses Σ or becomes part of Σ, which means sticking in the physical sense. The algebraic constraint of sticking is h(y(t))=0. While sticking the contact force Fc(t) must stay within physical bounds Embedded Image4.3and the vector field must be tangential to Σ, that is, Embedded Image4.4so that the solution continues on the switching surface. If such a contact force cannot be found, the solution crosses the switching surface and a discontinuity develops in the contact force. To calculate the contact force Fc(t) that makes the solution restricted to Σ, we substitute (4.1) into (4.4), which yields Embedded Image4.5Equation (4.5) involves the history of the contact force, which is either Embedded Image if h(y(t))≠0 or it is calculated from (4.5).

The question is whether the contact force Fc(t) is well defined during the stick phase by equation (4.5), which is an integral equation for Fc(t). To answer this, we need to consider possible singularities of L1(τ) at τ=0. As L1(τ) is bounded, one can find a maximal 0≤α≤1 and a positive constant C such that ∥L1(τ)∥<α. This is called the Hölder condition and α is the Hölder exponent. If α<1 we can also find a constant L1+ and positive C0 such that Embedded ImageThis means that L1(τ) is a sum of the singular L1+τα and a differentiable function. There are three cases to consider

  • (i) α=1, so that L1(τ) is differentiable. We assume that ∇h(y(t))⋅L(0)≠0.

  • (ii) α=0, so that L1(τ) is discontinuous and Embedded Image. We assume that ∇h(y(t))⋅L1+≠0.

  • (iii) 0<α<1, when L1(τ) is not differentiable, but continuous. Similarly, we assume that ∇h(y(t))⋅L1+≠0.

In case (i), when L1(τ) is continuously differentiable on Embedded Image, the integral term can be expressed using L(τ) as in equation (3.8). This is the case when the governing equations are finite dimensional or nj have finite norms. Therefore, the same dynamical phenomena should occur as in finite-dimensional systems, which cannot be resolved by our method. Applying (4.4) to equation (3.8), we find that the contact force obeys the integral equation Embedded Image4.6Owing to the differentiability of L1(τ) its derivative L0(τ), and therefore L(τ) must be bounded. When calculating the contact force by equation (4.6), the resulting Embedded Image can be discontinuous at the onset of the stick phase. Another cause of singularity is when ∇h(y(t))⋅L(0)=0, which can occur in case of the two-fold singularity [10].

The pre-tensed string model falls into case (ii). Owing to the discontinuity of L1(τ), equation (4.5) can be rearranged as a delay-differential equation with Embedded Image on the left-hand side, that is, Embedded Image4.7At the onset of stick at t the initial condition is Embedded Image. As all the terms in (4.7) are bounded Embedded Image must also be bounded. Therefore, Fc(t) is a Lipschitz continuous function of time when the solution gets restricted to Σ and all throughout the stick phase. It remains to be investigated what are the dynamical consequences when ∇h(y(t))⋅L1+=0 and whether the uniqueness of solution is preserved through such a singularity.

The Euler–Bernoulli beam falls into case (iii). First, we note that Embedded Image4.8which is by definition αΓ(α) times the α-fractional integral of Embedded Image [25]. We assume that the stick phase starts at time t. Using the rules of fractional integration, we find that Embedded Image4.9By separating the singular component of equation (4.5), we get Embedded Image4.10As we assumed that ∇h(y(t))⋅L1+≠0, it follows that all the terms on the right-hand side of (4.10) are bounded by C1. Therefore, we fractional integrate (4.10) with 1−α exponent exactly as in (4.9) and get Embedded Image4.11This means that if ∇h(y(t))⋅L1+≠0, Fc(t) is Hölder continuous with exponent 1−α. We note that Hölder continuity implies continuity in the traditional sense, therefore the friction force is continuous during the transition from slip to stick.

5. Stick–slip motion of a bowed string

To see our theory applied to a mechanical system, consider the example of a bowed string in figure 5a. We consider the same equation of motion as in §3b, where we derived all the necessary ingredients of the reduced model apart from the contact force. Because L1(τ) has a discontinuity at τ=0, this example falls into case (ii) of §4.

Figure 5.

(a) Schematic of a bowed string. The bow is pulled with a constant velocity v0, while the string exhibits a stick–slip vibration generated by the friction between the bow and the string. (b) Graph of the Coulomb-like friction force. (Online version in colour.)

To complete the model, we define the contact force Fc of equation (4.1) as the friction force between the bow and the string. We assume that the string is being bowed at ξ=ξ with velocity v0 that generates the friction force Embedded Image5.1where vrel is the relative velocity of the string and the bow. The graph of the friction force function can be seen in figure 5b. In this example, the static friction force is within the interval [−μ,μ]. The relative velocity between the string and the bow is expressed as a function of a resolved variable vrel=y2(t)−v0. We also use the relative velocity to define the switching surface Σ by h(y(t))=vrel. Therefore, the contact force of equation (4.1) becomes Fc(y(t))=fc(h(y(t))).

(a) Numerical method

We use a simple explicit Euler method to approximate the solutions of (4.1) and (4.5). We assume that time is quantized in ε chunks, so that yq=y(), fc,q=fc(), where q=0,1,2,…. In the case of slipping, the only unknown is the state variable yq that is calculated using the following formula: Embedded Image5.2where the friction force fc,q=Fc([yq]2v0) is used. The integration is approximated by the rectangle rule. For just illustrating the theory, such a crude approximation is sufficient while for better accuracy and efficiency higher order methods, such as the Runge–Kutta [26] method could be used. In our calculations, we keep the step size reasonably short at ε=5×10−4.

If the relative velocity h(yq)=[yq]2v0 of the string and the bow changes sign, there are two possibilities. Either the trajectory crosses the switching surface Σ or it will stay on Σ satisfying the equation Embedded Image5.3that is, the discretized version of Embedded Image. To test which case applies, we substitute (5.2) into (5.3) and solve for the friction force fc,q that would hold the string and the bow together, which becomes Embedded Image5.4While calculating fc,q, the relative velocity is forced to be zero. If the calculated friction force satisfies −μfc,qμ, the bow and the string stick together. For the stick phase of motion, we use equation (5.2) to advance the solution together with this dynamic friction force of equation (5.4). During the stick phase, using the friction force from equation (5.4) theoretically keeps the relative velocity zero. However, numerical errors can cause a small drift. To prevent the drift, equation (5.2) is solved only for the position variable, while the velocity is fixed at [yq]2=v0.

(b) Numerical results

To illustrate the properties of the dimension-reduced equation (4.1), we calculated a typical stick–slip trajectory starting at the initial condition z(0)=y1(0)w1+y2(0)w2 with y1(0)=−2.9224 and y2(0)=−2.7668. The parameters of the friction force in equation (5.1) are μ=4, κ=0.32, σ=1, the speed of the bow is Embedded Image, the damping ratios are Embedded Image, k=1,…,N and the wave speed on the string is c=1. The string is bowed at ξ=0.4.1 We solved equation (2.3) using Matlab's ode113 solver and applying the solution strategy described by Piiroinen & Kuznetsov [27]. Equation (2.3) was truncated at N=160 modes. Reduced equation (4.1) was solved using the algorithm described in §5a. The results of the simulations for the reduced and the full model, shown in figure 6a,b, are nearly indistinguishable, because the only approximations are within the numerical methods. The dark (blue) curves denote the solution of (4.1) and the (almost invisible) light (red) curves underneath represent the solution of (2.3).

Figure 6.

Solution trajectories of equations (2.3) and (4.1). (a) Displacement of the string and (b) velocity of the string at the contact point. (c) Friction force between the bow and string using equation (4.1). The continuous lines denote slip, the dashed-dotted lines represent sticking motion. (d) The discontinuity of the friction force at the onset of sticking disappears in the continuum limit. N indicates the considered number of modes. (Online version in colour.)

Initially, the solution spends short time intervals on the switching surface, and then settles into a periodic stick–slip motion. The stick phases can be recognized in figure 6b as short horizontal sections at y2=1.5. In figure 6c, the friction force is represented by the dark continuous (blue) lines and the light dash-dotted (green) lines for the slipping and the sticking motion, respectively. The friction force also appears to be discontinuous. To calculate this solution, we did not use the converged L1(τ), instead we used a series of mode truncations shown in figure 3b. On a smaller scale, figure 6d shows that the gap in the friction force (dashed line) between the slipping segment (solid line) and the sticking segment (dash-dotted line) of the friction force vanishes as an increasing number of modes of system (2.3) are considered. As the theory dictates, the gap should vanish in the infinite-dimensional case.

6. Conclusion

In this paper, we considered vibrations of structures that are composed of linear elastic bodies coupled through strongly nonlinear contact forces such as friction. The coupling was assumed to occur at point contacts. We introduced an exact transformation based on the Mori–Zwanzig formalism that reduces the infinite-dimensional system of ordinary differential equations to a description with time delay involving small number of variables. We found that the model reduction technique converges and contact forces become continuous even though the governing equation is discontinuous. We illustrated this novel technique through the example of a bowed string.

Through examples, we found that if natural frequencies scale linearly with the mode number, the contact forces are Lipschitz continuous during the transition from slip to stick. This is the case of the elastic string. If the natural frequencies increase faster than linear, the contact forces are only continuous. The Euler–Bernoulli beam exhibits such a behaviour, but it also allows infinite wave speed, which can be thought of as not physical. In reality however, every structure will have small-scale longitudinal vibration components with linearly scaled frequencies similar to the Timoshenko beam model [13]. We expect that if all the details are considered for a linear structure, the contact forces must always be Lipschitz continuous in time. This finding together with the new form of governing equations could be used in further studies to understand the source of non-deterministic motion [10].

The reduced equations are also structurally stable. Small perturbations to the memory kernel or other terms only deform solutions but do not change their qualitative behaviour as long as the qualitative features of the memory kernel are preserved. This is a clear advantage over finite-dimensional approximation of non-smooth systems, where small perturbations can cause qualitatively different solutions. Therefore, once the qualitative form of the memory kernel is established non-smooth mechanical systems can be approximated more successfully using our description.

How non-smooth phenomena of low-dimensional systems manifest in continuum structures is an open question. For low-dimensional systems, many singularities can occur that lead to chaotic and resonant vibration on invariant polygons [28], the Painleve paradox [11] and other types of discontinuity induced bifurcations. It remains to investigate how these phenomena occur in systems involving elastic structures, and hence equations with memory.

Our theory is developed for linear structures coupled to strong nonlinearities. It is however possible to extend this framework to cases where the underlying structure is nonlinear. For the weakly nonlinear case, the Hartman–Grobman theorem [29,30] guarantees the existence of a transformation that takes any weakly nonlinear system into a linear system about an equilibrium if that system is not undergoing a stability change. This generalization is currently being worked on by the author.

We also assumed point contacts in our derivations. This is a significant simplification because most contact problems occur along a surface. The difficulty arises when one needs to deal with contacting surfaces that slip at one part of the contact surface while stick at others. An interesting question is whether it is possible to develop a similar model reduction technique of such problems to involve only a finite number of variables.


The author thanks Gábor Stépán, who brought his attention to the work of Chorin et al. [18]. He also thanks Jan Sieber, Alan R. Champneys and John Hogan for useful discussion and comments on the manuscript. The detailed comments and suggestions of the anonymous referees significantly improved the quality of the paper and are greatly appreciated.

Appendix A. Model transformation

In this appendix, we show that the infinite-dimensional system (3.1) can be transformed into a finite-dimensional delay equation. The delay equation involves convolution integrals that can be related to Green's functions, but only for part of the system. The procedure is based on the variation-of-parameters formula.

Consider the following linear forced system: Embedded ImageA1where Embedded Image, Embedded Image and Embedded Image. Assume matrices Embedded Image and Embedded Image such that S=WV is a projection matrix S=S2 with a 2M-dimensional range. S is a projection if and only if VW=I2M, where I2M is the 2M-dimensional identity. Also, define the complementary projection matrix Q=IS and the resolved coordinates y(t)=Vz(t). With this notation, we rewrite equation (A1) into Embedded ImageA2Assume that the solution of Embedded Image can be computed for specific initial conditions. Therefore, the solution of (A2) can formally be expressed using the variation-of-parameters or Dyson's [29] formula as Embedded ImageA3Substituting this result into the second term on the right-hand side of (A2), we get Embedded ImageA4Note that Sz(t)=WVz(t)=Wy(t) and project (A4) using V, to get Embedded ImageA5which is the reduced equation for only the resolved coordinates y(t). Note that RQ eRQτ=(d/dτ) eRQτ, hence the integrals can be rewritten with Riemann–Stieltjes integrals as Embedded ImageA6where Embedded ImageA7 Embedded ImageA8 Embedded ImageA9 and Embedded ImageA10If the range of W is invariant under R, then K(τ)=A const. This occurs because the image of the range of W under R is in the kernel of Q. Consequently, the integral with K(τ) vanishes. The memory kernel L(τ) can be simplified in a similar manner. By the definition of the projection operator Q, we find that QW=0, and consequently QRW=0. This further implies that QRQ=QR(IWV)=QR, and thus (RQ)n=RQRn−1. Differentiating the fundamental matrix eRQt, we get Embedded ImageA11which means that Embedded ImageA12With this result, the forcing term and L(τ) become Embedded ImageA13

Note that this procedure is a simplified version of the Mori–Zwanzig formalism [18,31] for linear systems. Therefore, our procedure can be extended to nonlinear systems. In the nonlinear case, Ay, K(τ)y and L(τ)f become nonlinear functions of y and f, respectively.

Appendix B. Strong continuity of eRt

In order to justify our analysis, we need to show that eRt is a strongly continuous semigroup. We use the Hille–Yosida theorem [22,32], which states that eRt is a strongly continuous semigroup satisfying ∥eRt∥≤M0 if and only if Embedded Image and Embedded ImageB1for λ>0 and n=1,2,3,…, where Z is defined by equation (3.31).

First we show that if R satisfies (3.32) and that each eigenvalue of R has finite multiplicity, then condition (B1) is satisfied. Using a linear transformation matrix, R can be brought into its block diagonal Jordan normal form. Each block in the diagonal of the Jordan normal form corresponds to an eigenvalue λk of R and has size lk. The form of such a block is Embedded ImageB2After inversion the component of (RλI)n that corresponds to Jk becomes Embedded ImageB3The norm of this Jordan block can be estimated by Embedded ImageB4Note that |λkλ|≥|ℜλkλ|. As λ>0≥ℜλk, one can find an Mk such that Embedded ImageB5Considering this estimate for all Jordan blocks, we find that Embedded ImageB6where Embedded Image. This proves (B1).

To conclude the proof, we show that Embedded Image. Again, we use the Jordan normal form. We partition every vector in Z such that z=(z1,z2,…)T, where zk are of the size of a Jordan block. Let zZ and construct Embedded Image such that it has P number of non-zero components. This guarantees that for any zP, Embedded Image, hence Embedded Image. It is also clear that Embedded Image owing to its construction, thus we have shown that Embedded Image.

Appendix C. Boundedness and smoothness of Embedded Image

Definition (3.17) with (3.12) of Embedded Image has two terms both including the expression Embedded Image. Therefore, we only need to consider Υj(τ) in our analysis to show boundedness and smoothness of Embedded Image.

We use the inverse Laplace transform [22] to obtain Embedded ImageC1where Embedded Image. Using integration rules for the Laplace transform and multiplying (C1) by vectors vM+j and v from the left and right, respectively, we get Embedded ImageC2The equality makes sense if the integral converges. Evaluating the inverse operator in (C2), we find that Embedded ImageC3Note that it is sufficient to consider [Υj(τ)]M+l because [Υj(τ)]l is the integral of [Υj(τ)]M+l. Substituting (C3) into (C2), we are left with Embedded ImageC4The integral of inverse Laplace transform (C2) converges for all t≥0 if Embedded ImageC5and for j,l=1,…,M. Indeed, by estimating the bound we get Embedded ImageC6This implies that Embedded Image is bounded.

If stronger conditions (3.34) and (3.35) are satisfied, we can alter the contour of integration so that it is within the left half of the complex plane, Embedded Image, where δ>0 is sufficiently small so that all the eigenvalues of R are on the left of Γδ. We parametrize the contour Γδ by λ=κ e±i(π/2+δ), κ≥0. Using this contour, we find that the derivative Embedded ImageC7can be estimated by Embedded ImageC8This means that the derivative of Embedded Image is also bounded but only for τ>0. Owing to eRτ being a strongly continuous semigroup ([22], Theorem 2.5.2) higher derivatives are also bounded Embedded ImageC9for Embedded Image which proves that Embedded Image is smooth for τ>0.


  • 1 The choice of these parameters was guided by the desire of producing stick–slip motion rather than physical consideration.

  • Received September 4, 2013.
  • Accepted October 25, 2013.


View Abstract