Transient modelling of flexible belt drive dynamics using the equations of a deformable string with discontinuities

ABSTRACT Dynamics of a belt drive is analysed using a non-linear model of an extensible string at contour motion, in which the trajectories of particles of the belt are predetermined. The equations of string dynamics at the tight and slack spans are considered in a fixed domain by transforming into a spatial frame. Assuming the absence of slip of the belt on the surface of the pulleys, we arrive at a new model with a discontinuous velocity field and concentrated contact forces. Finite difference discretization allows numerical analysis of the resulting system of partial differential equations with delays. Example solution for the acceleration of a belt drive and an investigation of its frequency response depending on the velocity are presented and discussed.


Introduction
In the present paper we numerically study transient dynamics of a belt drive based on the model, which has earlier been presented by the authors in [1][2][3]. Velocity and strain of the extensible belt vary between the tight and the slack sides of the drive, which leads to the sliding friction between the belt and the pulleys. This model of creep of the belt near the points, where it leaves the pulleys, was first studied by Reynolds [4] and is adopted in the engineering and technical literature, see e.g. [5][6][7]. Assumptions, which are common for the technical literature, ignore dynamic effects due to the acceleration of the particles of the belt; other hypotheses concerning redistribution of tension forces between the tight and the slack sides of the drive are often involved.
A whole body of literature deals with the dynamics of axially moving free spans of the belt, which are considered decoupled from the segments in contact with pulleys by assuming particular boundary conditions, see [8][9][10]. Modelling the belt by string mass and momentum. The range of applicability of this idealized contact model to real structures is yet to be justified in the future work. Further treatment of more sophisticated and practically relevant models, featuring transverse vibrations of the belt with varying contact zones, bending and shear stiffness of the belt as well as dry friction contact law is foreseen. The analysis will rely on the present strategy of spatial description of the entire belt, and the results will be presented in subsequent publications.
In section 3 of the present work we shortly recall the results of the analysis of steady operation of belt drives, presented by the authors in [2]. The benefits of the spatial description and of the idealized model of belt-pulley contact are demonstrated in section 4, in which we outline a strategy for simple and consistent transient analysis, earlier sketched in [3] n the form of a mathematical model without any discussion of possible numerical analysis. Subsequent sections of the paper present a novel numerical scheme for the solution of the system of partial differential equations with delays in the boundary conditions. Based on the finite difference discretization with respect to the spatial variable, the scheme has been implemented in Wolfram Mathematica 1 environment, which was proven to be efficient for numerical modeling of problems of structural mechanics [27].
An earlier version of this work has been partially presented in the proceedings of 8th Vienna International Conference on Mathematical Modelling -MATHMOD 2015 [28]. The present contribution includes additional numerical investigation concerning the harmonic excitation and frequency response of a belt drive in the vicinity of a stationary regime.

Equations of belt dynamics at contour motion
A two-pulley belt drive is schematically presented in Fig. 1. At contour motion, the material points of the belt follow the trajectory, which is plotted by a thick line. The material particles travel across the boundaries of the contact regions 1. . . 4. This makes spatial description of the dynamics of the belt as an extensible string advantageous in the present case. In this section we recall and revise the equations of string dynamics, which were derived in [2].
The considered trajectory can be parametrically specified by the dependence of the position vector R on the arc coordinate σ with the following geometrical rela-tions: in which k(σ) is the curvature, τ is the unit tangent vector and n is the unit vector of normal to the trajectory. In Fig. 1 the contour coordinate is looped, and the values σ = 0 and σ = σ 4 correspond to the same point 4. The domains 0 ≤ σ ≤ σ 1 and σ 2 ≤ σ ≤ σ 3 answer respectively to the driving pulley A with k = k A and the driven one B with k = k B . The slack span σ 1 ≤ σ ≤ σ 2 and the tight one σ 3 ≤ σ ≤ σ 4 are straight with k = 0. At spatial description, our aim is to find the variation of the material coordinate s of the string for a given point in space over the time t: In the reference configuration s = S(σ, 0) = σ. Stretch of the string equals to ∂σ/∂s ≡ ∂ s σ, and the strain measure is expressed using the rule for an inverse function: we denote the derivative with respect to the spatial coordinate σ with a prime. The tension force Q is related to the strain by a linear constitutive law in which Q 0 is the tension in the reference configuration (pre-tension) and b is the stiffness of the belt. The velocity of a particle of the belt is the material time derivative Transforming to the spatial description, we find in which the derivative with respect to time t at a given point in space σ = const is denoted with a dot. The equation of balance of forces [15] ∂ s (Qτ ) + q = ρw (7) features the vector of the external force q and mass ρ per unit length of the belt in the reference configuration ds. At spatial description, components of the vector of acceleration w ≡ ∂ t v take the form Introducing the force factorq and belt density m per unit length of the contour in space dσ according toq ≡ qS , m ≡ ρS (9) and projecting onto the local basis vectors τ and n, we rewrite the equations of balance (7) in a spatial form: The equation of balance of mass in a spatial description reads in which we integrate with respect to the material and spatial coordinates over the same domain, s −,+ = S(σ −,+ , t), and ds = S dσ. At the driven pulley with the radius k −1 B the loading moment M B (t) is computed as follows:

Steady operation
In a steady regime the velocity, force and strain at a given point in space do not change in time. We have v = v(σ), Q = Q(σ), ε = ε(σ), and from (3) it follows that Now using (6) we find S: This relation between the material and spatial coordinates at steady motions of strings was established earlier by Nordenholz and O'Reilly [16], and was applied in the analysis of steady operation of belt drives by Rubin [13]. From (11) we find the "mass flow rate" mv = ρc = const. The length of the belt in the reference configuration L equals to the length of the contour, and the variation of σ in the range from 0 to L ≡ σ 4 answers to the growth of s from 0 to L. The condition S(L) − S(0) = L determines the constant c (this is the "condition of compatibility" in the terminology of Rubin [13]).
The strain and the tension force can now be expressed through the field of velocity v: The tangential projection of the balance equations in (10) leads tõ At a straight free span q = 0, and from (16) follows as Q depends on v according to (15). At the tight span we have and at the slack one see Fig. 1. At steady motion the moment at the driven pulley can be related to the velocities at the free spans. Substituting (16) in (12) and integrating, we obtain At the pulleys we assume perfect contact: The seeming contradiction between (12), (20) and (21) is resolved by admitting discontinuities and corresponding concentrated forces at the points σ = σ i : only the "outermost" zones on the pulleys are loaded. For the loading moment we get According to (16), the concentrated forces result from the jumps of the tension force and of the velocity in the corresponding points; see Irschik [26] for a discussion. Traditional treatment of the creep model [4,13] assumes that the belt is contracting near the point σ 1 , where it leaves the driving pulley A, and extending near the point σ 3 , where it leaves the driven one B. For a friction belt drive a continuity condition needs to be imposed in the points σ 2 and σ 4 , where the belt hits the pulleys. Indeed, according to the idealized model of contact the particles of the belt adhere to the surface of the pulleys at these points, and any jump here would have been transported further. Thus we conclude that in a steady regime the field v(σ) is piecewise constant with two points of discontinuity σ 1 and σ 3 . The velocities at the free spans coincide with the velocities of the surfaces of the pulleys: To simplify the subsequent analysis, from now on we assume the drive to be symmetric with Then from the compatibility condition S(L) − S(0) = L (the material length of the belt remains constant) and from the relation between the material and spatial coordinates (14) we find as the velocity is piecewise constant in the domain of integration. The working regime is determined by the ratio of the velocities at the tight and at the slack sides and for the moments at the pulleys we find The coefficient of efficiency of the drive is more power is lost at higher transmitted moments. This simple expression for the coefficient of efficiency is a consequence of the relation (25), which is based upon the idealized contact model. The analysis with the finite slip zones results in more complicated expressions, see, e.g., [13].
In the technical literature it is commonly assumed that the mean value of the tension forces at the free spans equals to the pre-tension (see e.g. [5,7]). However, it is easy to see that and the difference grows at higher transmitted moments. Moreover, further additional analysis shows that the proposition does not generally hold for unsymmetric belt drives, when k A = k B .

Transient dynamics: mathematical model
For both free spans with q = 0 a single partial differential equation follows from the equalities (3)-(10):S Indeed, expressing the tangential acceleration component w τ via S using (6) and (8), and then finding the force Q by substituting (3) in (4), we find that the tangential component of the balance equation (10) becomes equivalent to (30). This wave equation can be classified as quasilinear, as in the material description for σ( The spatial form (30) is more advantageous as the boundary conditions are then placed at fixed points σ i . As discussed before (23), the particles of the belt immediately adhere to the pulley surface at the points σ 2 and σ 4 . In these points we demand the velocity of the belt to be continuous, which leads to the boundary conditions for the equations of dynamics at the free spans (30) Discussing the remaining conditions at the points σ 1 and σ 3 , we focus our attention at the conditions for the driven pulley B. The equations for the driving pulley A are to be formulated by replacing the indices accordingly.
Moving together with the surface of the pulleys, the particles of the belt preserve constant distances to their neighbours. Thus, neither the strain nor the tension force vary in a given material point during the contact stage. In the point σ 3 we observe the material points from the point σ 2 with a certain delay: the particle with the material coordinate S 3 (t) came into contact with the pulley at the time instance t − τ B (t). Between the time points t − τ B and t the particle has moved from σ 2 to σ 3 together with the pulley surface. This means that during this time span the pulley has rotated to the known angle: Differentiating with respect to time, we find the relation of the delay time to the angular velocity of the driven pulley: .
The problem is closed by the equations for the pulleys The moments of inertia of the pulleys are denoted as I. External moments M ext A , M ext B are either given functions of time or may be governed by additional equations for a control loop. Computing the moments M A and M B , which act on the pulleys from the side of the belt, we again focus the attention at the driven one and substitute (10) in (12): The fields m, v andv have a discontinuity of the first kind at σ 3 . "Concentrated" terms under the integral appear only due to v : In the domain σ 2 < σ < σ 3 we havev = k −1ω B , v = 0, and the integral (37) results in the following three terms: Here S 3 is again determined by the solution in the point 2 with a delay. In the expression for the jump of (S ) −1 the value at σ 3 + 0 is taken from the solution on the span 34, and the stretch at σ 3 − 0 is "transported" by the pulley: A simple interpretation can be given to all three terms at the right hand side of (39). The first term is just the moment due to the tension forces and is taken into account by all technical approaches. The second term represents additional inertia of the pulley due to the attached part of the belt. And the third non-trivial term is the reactive force at the point, where the velocity undergoes a jump.
It should be noted, that the expression (39) could be derived using relations of jump of mass and momentum known from the literature. Following [26], we write the concentrated force, which acts on the belt in the point σ 3 : The force F 3 contributes to the moment M B as in (22), and the jump of mv 2 evidently leads to the third term (reactive force) at the right hand side in (39).

Numerical scheme
A central point of the simulation is the numerical solution of the partial differential equation (30). A system of ordinary differential equations, which can then be integrated over time using conventional methods, follows after a finite difference discretization with respect to the spatial coordinate σ [27,29]. We divide both free spans 12 and 34 into n segments and consider values of the unknown function at the equidistantly distributed nodes σ (i) , i = 0 . . . n. The outmost nodes correspond to the end points of the domain, i.e., for the slack span and for the tight span. We project the equation (30) onto the internal nodes of the grid i = 1 . . . n − 1 by using central finite difference approximation for the spatial derivatives and obtain n − 1 ordinary differential equations for each free span; here h = σ (i+1) − σ (i) is the step of the grid. There are n unknowns S (1) (t), . . . , S (n) (t), as the values S (0) at the starting nodes of the grid σ 1 and σ 3 are directly available from the boundary conditions (33), which couples the solutions at both free spans. At the end points of the grid the particles immediately adhere to the surface of the pulleys, which means boundary conditions of the form The first equality in (38) provides a relation between S and v, and for the stretch at the end point of the grid we use a left-sided finite difference approximation: The conditions (46) complete the set of differential equations for each free span. Delay times τ A and τ B are additional unknowns in the problem with their own differential equations (35). Finally, for the angular velocities of the pulleys ω A and ω B we integrate (36) over time. The resulting system of equations features 2(n − 1) + 2 + 2 unknown functions of time.
According to (39) and (40), the moments M A and M B depend on the solution at the time instants t − τ A and t − τ B owing to both (33) and (35). In the numerical simulations we need to solve a system of delay differential equations with variable delay times. An implementation in Wolfram Mathematica has been developed on the basis of the regular routine NDSolve for solving systems of ordinary differential equations. After each step of the time integration, the values S 2 ,Ṡ 2 , S 4 ,Ṡ 4 , ω A and ω B are saved into a data array for later use. The values of the unknowns at the time instants t − τ A and t − τ B at the right-hand sides of the differential equations are then interpolated over this array. Stretches (S ) −1 in (40) are expressed viaṠ using the first equality in (38) as the velocities at the pulleys are known from (46). The accuracy of the time integration is justified by the convergence of solutions with respect to the spatial discretization (at growing n).
It remains to notice, that the delays in the differential equations mean that the solution depends on the "pre-history" at the time period t < 0; conventional initial conditions at the time instant t = 0 are insufficient in the present case. It can easily be explained by the fact that the motion of the system is determined by the initial distribution of S(σ, 0) not only at the free spans, but also on the pulleys. Because the domains along the pulley surfaces are not explicitly included in the numerical scheme, they influence the solution with delays. As an alternative, one could avoid the delays by including the whole range 0 ≤ σ ≤ L in the finite difference scheme with an advection equation in the contact domains. Conducted numerical experiments have shown, however, that this approach not only increases the number of unknowns, but also it does significantly affect the stability and the accuracy of the time integration. Indeed, it is known that the upwind difference schemes for the advection model introduce severe numerical diffusion/dissipation in the solution where large gradients exist [29], and its stability regarding the time integration is guaranteed only when the time step size is small enough. Moreover, matching various difference schemes at the boundaries between two domains is an additional non-trivial task.
Assuming the drive to be freely moving at a constant speed, such that v = v 0 at t < 0, we find all necessary functions for this "pre-historic" period of time: Starting from a moving state with v 0 > 0, we avoid difficulties with the infinite delay times in the beginning of the simulation. This initial velocity may be chosen small enough compared to the working regime.

Results of simulations: start up
The present model features dynamic effects, which are ignored by traditional approaches to modelling dynamics of belt drives; namely, the acceleration and deceleration of the particles of the belt when leaving the pulleys as well as the axial vibrations of the belt in the free spans, coupled with the oscillations of the pulleys. To make these effects prominent, we considered high rotation speeds and light pulleys. Choosing Young's modulus of the belt E = 10 8 , its volumetric material density ρ 3 = 1200 and its cross-section as a square with the size a = 5 × 10 −3 (all values are in SI units here and below), we compute its properties as follows: The radius of the pulleys and their moment of inertia are such that I is 20 times larger than the moment of inertia of the part of the belt around the half of a pulley. The distance between the centers of the pulleys is At t < 0, the structure is moving with the velocity v 0 = 10. Starting from the steady unloaded motion (48) at t = 0, we apply an accelerating moment at the driving pulley M ext A = M = 100. After a certain period of timet we apply the same, but oppositely acting loading moment at the other pulley, such that Depending on the size of the time intervalt, the system accelerates to a certain speed level, and then keeps moving in a steady mode after the decay of the oscillations. The structure is very light with short characteristic times, and we choose a value of the acceleration timet = 0.06 to make dynamic effects visible. The computed time histories of particularly interesting characteristics of the solution are presented in Fig. 2-6. We used n = 60 segments of the finite difference discretization, as further refinement does not lead to noticeable changes in the presented results. The angular velocities of both pulleys grow during the acceleration stage, and after the drive is balanced by a loading moment at the driven pulley, the system quickly arrives at a steady regime of motion. The high ratio between the stationary angular velocities η = ω A /ω B ≈ 1.66 means high levels of deformations of the belt, and indeed at the tight span we have strain ε ≈ 0.33, while at the slack span ε ≈ −0.2. Fast rotation of the pulleys leads to very short delay times τ A,B at the steady regime, which is clear from (34). Rotational moments, which act on the pulleys from the    (27), while the constant c ≈ 142.8 follows from (26). During the acceleration stage, the velocity varies almost linearly along the free spans, and at t =t the difference between the two ends is approximately 2%. It remains to notice, that the reached velocities are not too far below the critical value b/ρ ≈ 288.7, at which the physical transport speed equals to the wave speed and the solution of the transport problem undergoes a divergence instability [30]. Indeed, increasinḡ t further we observed "exploding" solutions as soon as the theoretically predicted critical speed level was exceeded.
In Fig. 7-8 we have shown the same results for the short initial period of time, such that the starting stage of the motion can be analyzed. The moment at the driven pulley appears after it is reached by a wave travelling through the free spans. The middle points of the free spans react to the applied driving moment after a two times shorter period of time, and several waves travelling across the belt can clearly be seen.

Results of simulations: response to a harmonic excitation
In the next step of numerical simulation, we again bring the system into the stationary motion regime by applying two moments one after another with a particular time intervalt. But now we wait until the decay of oscillations, induced by the acceleration phase, and then add a small harmonic excitation at the driven pulley such, that after a certain period of time the external moments are  Fig. 9. After a transient stage the system arrives in a regime of motion with small harmonic oscillations in the angular velocities of the pulleys as depicted in Fig. 10. These results correspond to the acceleration timet = 0.0075, such that the velocity of the stationary regime is lower than in the previous example and it takes longer time for the vibrations after the acceleration stage to decay. The harmonic excitation is applied at the time instant t = 50t = 0.375 and has a frequency λ = 315. Varying the excitation frequency along a certain range and running the time integration, we determined the steady-state amplitude A of the harmonic part of the angular velocity of the driving pulley ω A . This provides us with a frequency response of the system on the considered type of excitation. The computations were repeated for different velocities of the background stationary regime, which was achieved by changing the acceleration timet in (52). The resulting curves for three values oft are shown in Fig. 11. A simple analytical prediction for the resonance frequency of a system at rest can be computed using a simple model with two degrees of freedom, in which the belt acts as two springs, connecting the upper and lower points of the pulleys: in the denominator we accounted for the added inertia from the part of the belt on a pulley and also added half of the inertia of the free span of the belt; see (49)-(51) for the numerical values of the parameters. While the resonance peaks of slowly moving systems with smaller values oft are close to λ * , the reasons for the non-monotonous behavior of the maximal vibration amplitudes on the velocity of the system are yet t=0.00375 = * =314.9 Figure 11. Frequency response: steady-state amplitude of the angular velocity of the driving pulley depending on the frequency of excitation for different velocity levels of the stationary regime; the resonance frequency for a non-moving system is depicted by a vertical dashed line to be investigated in the future work.

Conclusion
We have presented a mathematical model for transient analysis of dynamics of a belt drive under the assumption of contour motion of the belt. A numerical scheme for solving the resulting partial differential equations with solution dependent delay times has been developed and discussed. The presented simulation results concerning steady acceleration with subsequent decay of oscillations are justified by the theoretical analysis of stationary regimes of motion. The numerical study of harmonic response in the vicinity of a steady motion shall be accomplished by an analytical investigation with the help of perturbation methods. This shall provide a solid basis for testing the strategies of automatic control for belt drives at high and rapidly changing speeds. The demonstrated mathematical model shall also serve as a basis and as a reference for the further ongoing research in this field, which will make use of the Eulerian description and account for more sophisticated effects, such as transverse vibrations of the belt, its bending and shear stiffness and finite creep zone.