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

Robert Szalai

## Abstract

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 2.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 to simplify the notation. As in infinite dimensions not all norms are equivalent, we assume that , which turns into a Banach space [17]. Owing to linearity the governing equations can be written as 2.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 2.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 2.4as illustrated in figure 1. We call the positions yj(t) and the velocities of the contact points resolved variables. We assume M contact points, thus we define . Substituting (2.1) into (2.4), we obtain the motion of the contact points as a function of the solution of equation (2.2) 2.5where 2.6Vectors nj are assumed to be linearly independent, spanning an M-dimensional subspace of .

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 2.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 3.1where 3.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 3.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 3.4which is a linear mapping from to . 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, 3.5Equivalently, wj can be constructed as linear combinations of 2M number of eigenvectors of R, because eigenvectors are invariant by definition. If are some eigenvectors of R and is the matrix with in its columns, the R invariant lifting operator becomes . The matrix is invertible if are not in the kernel of V. Also, both V and 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 3.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 3.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 3.8where 3.9 3.10 and 3.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 3.12where 3.13As we have a formal expression for Lj(τ) in the form of equation (3.11), we can calculate the integral 3.14From the simple rule of differentiating an exponential, we find that . Using this formula twice in equation (3.14), we are left with 3.15Assuming that is bounded (see §3d) limit (3.13) becomes 3.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 implies that either exponentially if the system is dissipative or oscillates with zero mean. Therefore, we define the integral 3.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: 3.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 3.19where .

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 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: 3.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 . The natural frequencies of the system are ωk=ckπ and we assume uniform damping for all modes, which gives us 3.21The vibration at the contact point can be expressed as a linear combination of all the modes u(ξ,t)=nx(t), where . The resolved variables therefore are y1(t)=nx(t) and . We also choose in formula (3.6), thus 3.22

Next, we calculate the function eRtv2 that appears in expression (3.11) of L(τ). The equation whose solution we are seeking is , which is expanded into 3.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 3.24Without assuming the form of z1,k(t) and evaluating formula (3.11), we find that 3.25Eventually substituting (3.24) into (3.25) in the conservative case (Dk=0) [L(t)]2 of our system (3.20) becomes 3.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 3.27Note that this is times the static displacement of the string under unit load at ξ=ξ. Using the expressions for z1,k(t) and integrating as per definition (3.17), we get 3.28Assuming that Dk=D are constant, the right limit of [L1(τ)]2 becomes 3.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 . (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 3.30The natural frequencies of (3.30) are determined by the equation , which can be approximated by 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 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 . 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

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 3.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 . Even if we know that , the bound of cannot be directly estimated in the straightforward way, because . However, the two examples in the previous sections show that can be bounded, but not necessarily smooth. In case of the string and without damping, is a piecewise-smooth function, while appears to be continuous but non-differentiable for the undamped Euler–Bernoulli beam. If damping is introduced, 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, 3.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, 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 are bounded if there are such that 3.33

Smoothness of 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 3.34The stronger condition that guarantees the smoothness of for τ>0 becomes 3.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 3.36

Condition (3.33) has a mechanical meaning. When the structure is forced at contact point χj with , γ>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 , we require that the decaying forcing 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 , which implies that 3.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 . In particular, for the two examples of the string and the beam we have 3.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 4.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 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, 4.2Without restricting generality, we assume that . 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 4.3and the vector field must be tangential to Σ, that is, 4.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 4.5Equation (4.5) involves the history of the contact force, which is either 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 This 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 . 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 , 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 4.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 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 on the left-hand side, that is, 4.7At the onset of stick at t the initial condition is . As all the terms in (4.7) are bounded 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 4.8which is by definition αΓ(α) times the α-fractional integral of  [25]. We assume that the stick phase starts at time t. Using the rules of fractional integration, we find that 4.9By separating the singular component of equation (4.5), we get 4.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 4.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 5.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: 5.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 5.3that is, the discretized version of . 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 5.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 , the damping ratios are , 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.

## Acknowledgements

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: A1where , and . Assume matrices and 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 A2Assume that the solution of 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 A3Substituting this result into the second term on the right-hand side of (A2), we get A4Note that Sz(t)=WVz(t)=Wy(t) and project (A4) using V, to get A5which 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 A6where A7 A8 A9 and A10If 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 A11which means that A12With this result, the forcing term and L(τ) become A13

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  and B1for λ>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 B2After inversion the component of (RλI)n that corresponds to Jk becomes B3The norm of this Jordan block can be estimated by B4Note that |λkλ|≥|ℜλkλ|. As λ>0≥ℜλk, one can find an Mk such that B5Considering this estimate for all Jordan blocks, we find that B6where . This proves (B1).

To conclude the proof, we show that . 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 such that it has P number of non-zero components. This guarantees that for any zP, , hence . It is also clear that owing to its construction, thus we have shown that .

## Appendix C. Boundedness and smoothness of

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

We use the inverse Laplace transform [22] to obtain C1where . Using integration rules for the Laplace transform and multiplying (C1) by vectors vM+j and v from the left and right, respectively, we get C2The equality makes sense if the integral converges. Evaluating the inverse operator in (C2), we find that C3Note 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 C4The integral of inverse Laplace transform (C2) converges for all t≥0 if C5and for j,l=1,…,M. Indeed, by estimating the bound we get C6This implies that 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, , 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 C7can be estimated by C8This means that the derivative of 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 C9for which proves that is smooth for τ>0.

## Footnotes

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