A model for contact and friction between beams under large deformation and sheaves

This work focuses on the modeling of contact between sheaves and flexible axially moving beams. A two-dimensional beam finite element is employed, based on the absolute nodal coordinate formulation (ANCF) with an improved selective reduced integration for the virtual work of elastic and viscous damping forces. For the efficient modeling of contact between flexible axially moving beams and sheaves in systems such as belt-drives or reeving systems, a number of newly developed algorithms is presented. The computation of normal contact is based on a penalty formulation using a spring-damper model, while for the efficient contact detection a bounding box which fits the exact dimensions of the finite elements is employed. For the detection and computation of contact, the beam elements are divided into linear segments. The modeling of tangential contact is based on a bristle model for friction extended for being compatible with an implicit time integration. A numerical example of a belt drive showed good convergence and agreement with analytical solutions.


Introduction
Mechanical systems including contact between highly slender beams and circular objects are met in ropeway systems [1-3] cranes and elevators [4,5], belt drives and conveyor belts [6].Computing the normal and tangential contact between circular objects and axially moving beams under large deformations is not straightforward.As stiffness may be high, the formulation should be compatible with an implicit time integrator.Especially in the case that we aim to solve industrial-scale problems or to perform an optimization through parameter variations the need to save simulation time becomes stronger which leads us to a coarser mesh refinement.To meet these goals, we focus on a segment method, which allows coarse space-wise discretization and an implicit integration which can cope with large time steps.
Various contact detection algorithms for rods undergoing large deformations exist in the literature often focusing on beam to beam contact.A pioneer work on beam to beam contact modeling was the one of Wriggers and Zavarise [7] which has introduced the so-called point-to-point formulation modeling the contact force as a point force at the closest point between two curves.Other approaches, e.g.[8], are using the line-to-line modeling which assumes distributed contact forces and is most commonly based on the mortar element method, [9].A few works on beam to beam contact are using Gauss-points-to-segment modeling, e.g.[10].Likewise, works on modeling the contact between flexible ropes and circular objects often use a point-wise contact modeling using points along the beam with respect to which they define the relative position of circle and rope, the indentation and they use them to apply the contact forces, e.g., [5,[11][12][13].The present work introduces an algorithm for tracking the circle's and rope's relative position which is based on linear segments to subdivide the beam's geometry which is modeled by third order polynomials.A similar algorithm was used in [14] where a flexible belt is discretized by 4-node shell linear elements.The proposed algorithm overcomes the limitation of methods based on the use of integration points to detect the contact with a circle if it lies on the interval between two nodes.
Existing literature for modeling contact forces is based on continuous approaches either using a continuous force model [5,[15][16][17][18], or the unilateral constraint methodology.The latter is based on the linear complementary problem [19][20][21] or the differential variational inequality [22].The continuous force model or penalty method has widely been used in the literature either using linear spring elements [11,12], or nonlinear spring elements [5].In the current work, the penalty method is preferred as it is compatible to an existing multibody code and because for soft materials such as belts, the contact stiffness is no limiting factor.
For modeling the tangential contact force, many works, e.g.[11,12,17,18], are using a regularization of the Coulomb friction law that gives a continuous tangential force, denoted as tri-linear Coulomb friction law or creep-rate-dependent friction, as introduced by Rooney and Deravi [23].A variation of this model is used in [13].A work on contact modeling in reeving systems, [5], employs a bristle friction model which tracks the slip and stick regions on the contact region and updates the tangential contact force accordingly.This model is originally proposed for control-systems applications [24], see also review paper [25].In the present work, a bristle friction model has been significantly extended for being compatible with the use of linear segments and an implicit time integration.
For the numerical modeling of beams undergoing large deformations in contact with circular objects, many researchers have used the Absolute Nodal Coordinate Formulation (ANCF), e.g., [13].Kerkkänen et al. [12] showed how an ANCF beam model is suitable for modeling belt-drives, while Dufva et al. [11] explored the effect of bending stiffness in the beam model when modeling structures that include contact of pulleys and beams.The ANCF beam model used in the current work is based on an extension for the original ANCF element [26], inferred using corrected expressions for the strain energy and bending moment relations.
The novelty of the present work can be summarized by an improved implementation of ANCF elements with a selective reduced integration, a contact modeling which uses linear segments for dividing higher order interpolation functions, the use of a bounding box which allows effective contact detection and a bristle model with history variables which tracks accurately the stick and slip zones and reproduces the solution from classical belt theory.

ANCF beam framework
For the numerical modeling of the beam, we are using the ANCF beam element in the extended version of Ref. [26] with slight improvements as shown in the following.An element with nodes N j 1 and N j 2 is described by the generalized element coordinates, q = q ( j 1 ) T q ( j 2 ) T T , in which each node N j is related with containing the position vector of the two nodes and their spatial derivatives, r .Throughout this paper, the abbreviation () = ∂() ∂ x is used.The geometry of an ANCF beam element with reference length L is defined with the use of the local coordinate x ∈ [0, L] as r( x, t) = S( x)q(t) and r ( x, t) = S ( x)q(t) , (1) in which S is the matrix of the shape functions defined as The shape functions used are the third-order polynomials above which guarantee C2 continuity for the elements.For the derivation of the equations of motion, we are using Lagrange kinematics and the virtual work of applied, elastic and visco-elastic forces.The virtual work of elastic forces is derived following the definition of the axial strain and curvature proposed in [27].The virtual work of elastic forces is derived using the axial strain ε = r − 1 and the material measure of curvature K = r ×r r 2 , in which r × r is the cross product of planar vectors as defined in [26].The virtual work of elastic forces reads in which δW ea and δW eb are the virtual works of axial and bending forces, while δε and δ K are the variations of ε and K , E is the Young's modulus, E A the axial rigidity for a cross section A and E I the flexural rigidity for second moment of area I .The elastic forces are integrated numerically using selective reduced integration based on Gauss quadrature numerical integration rules.Thus, the integrals for δW ea and δW eb are computed separately allowing an effective combination of two Gauss quadrature integration rules.The implemented integration schemes, that is the combined integration rules for the virtual works of axial and bending forces, make use of the Gauss Legendre (GL) and Gauss Lobatto (LO) integration rules as shown in Table 1.These three integration schemes are available in the multibody dynamics simulation code Exudyn 1 which is used for implementing the numerical example performed in Sect. 4.
Integration points and weights for the aforementioned integration rules can be found in [28].If we denote by (A), the Gauss quadrature rules used for δW ea and by (B) the ones used for δW eb , and m, n are the numbers of integration points used for δW ea and 1 https://github.com/jgerstmayr/EXUDYN. δW eb , respectively, we can rewrite (4) as in which x i and w i are the integration points and weights, respectively.For example, ε(x ) is the axial strain at the integration point i of (A m ) integration rule.
In previous works employing the above-described ANCF beam element, [26] the integration schemes a and b, see Table 1, had been considered.For the investigations performed in this paper, the newly introduced scheme c is selected.This scheme resolves the problem of oscillations occurring in the axial strain and requires less computational efforts.The improvement of this integration scheme follows from the fact that axial and bending strain terms are evaluated at completely disjointed locations.In specific, the axial strain terms are evaluated at 0, L 2 and L, while the bending terms at 3 .This allows axial strains to freely follow the bending terms at L 2 ± L 2 1 3 , while the axial strains are almost independent from bending terms at 0, L 2 and L. However, note that the numerical integration is highly reduced in this case which might cause spurious modes in some applications, [29].
For the virtual work of viscous damping forces, we are using the approach of [27], with the damping parameters d ε and d K .The time derivatives ε and and The virtual work of the viscous damping forces in ( 6) is integrated according to (5).

Normal force model
For modeling the normal contact force, we are using the model of Lankarani and Nikravesh [30]: in which g represents the gap between the beam element and v n is the normal relative velocity.The exponent n in the first term originally comes from the Hertz contact law and is related with the geometry and the material of the contacting surfaces, [23].For the case of our belt drive where contact occurs between a flat belt and cylinders, we use n = 1 throughout in our implementation.In the second term, d c is a function of g (d c ∝ g n ) which makes the second term nonlinear, [30].Although this model is proven to be suitable for contact force models, for the contact models we examine where the normal relative velocity is very small we use a constant d c , as a simplification.
The algorithm for defining g and v n is described in the following section.

Contact geometry and kinematics
Bounding box for ANCF A boxed search algorithm is used in order to improve the performance.This boxed search specifies which beam elements are potential candidates for having contact with the circular objects based on their relative location.
The proposed boxed search algorithm is based on the use of a bounding box the dimensions of which fit exactly the beam finite elements.
Since the shape functions defined in (3) are thirdorder polynomials (1) can be rewritten as with a 0 , a 1 ,... and b 0 , b 1 ,... being constants for the bounding box calculation.The coordinates of the i th bounding box are r x,min,i , r y,min,i , r x,max,i , r y,max,i see Fig. 1, in which and r y,min,i , r y,max,i accordingly defined.For finding the extrema of r x , we compute its derivative, and by setting r x equal to 0 and we get the two roots x1 and x2 .Hence, we exclude the according r x ( x1,2 ) from ( 13).The y-coordinates of the bounding box are computed likewise, using ( 11)-( 13) for r y .Detailed search For computing the relative position of a beam element with a circular object, cubic interpolation functions for the beam are subdivided using n s linear segments, see Fig. 2b.Then, we examine the relative position of the center C i of the circle, located in p c i , and the linear segment s i , lying between p i and p i+1 of a beam element, see Fig. 2a.
Using the vectors v s = p i+1 −p i and v p = p c i −p i , the projection of vector v p on v s is defined as, Fig. 1 The dimensions of the bounding box are fitting exactly the beam elements Fig. 2 in which we used the substitution ρ = We distinguish three cases with respect to ρ: p i+1 and the distance, d, can be computed as Therefore, the gap between the circle and the segment used for defining contact and for calculating f n , see (9), is We approximate the velocity of the point P by interpolating linearly the velocities of p i and p i+1 as follows: The normal and tangential gap velocity which are used for the computation of normal and tangential contact forces are and respectively.In ( 18) and (19), n denotes the normalized vector of d, d = p p − p c i in which p p the position vector of P, and t denotes the perpendicular normalized vector, see Fig. 2.

Tangential force model
Existing tangential contact models, e.g.[13,17], are based on a regularized Coulomb friction that makes use only of a velocity penalty.Such models suffer from certain limitations as they are not suitable for static or quasi-static simulations where no relative velocity exists.
In this section, we describe an extended tangential model which makes use of a bristle stiffness μ k .The incorporation of this bristle stiffness is inspired by the LuGre model introduced in [24] which reproduces spring-like behavior for small relative displacements during sticking.Similar to the LuGre model, we consider a sticking phase, in which the relative displacement is constrained by a virtual spring.The tangential force applied in sticking phase reads: (20) in which v t was defined in (19), μ v is an optional penalty factor (can be also understood as friction damping) and Δx stick is the displacement of the contact point while sticking, see Fig. 3.When f (lin) t > μ• f n , with μ being the friction coefficient, the friction spring breaks and the tangential force is modeled according to Coulomb friction.The resulting tangential model can be summarized as follows Although a bristle friction model was already proposed for modeling the contact in reeving systems [5], the model introduced in the present work is significantly different.First, for computing Δx stick we make Fig. 3 Contact stiffness visualized with a spring use of a set of history variables.In this way, we avoid the introduction of an additional dynamic parameter for the deflection of the bristles, [24].Furthermore, the proposed model has been extended from a point-tocircle contact to a segment-to-circle contact.Finally, it has been adjusted for being compatible with an implicit time integration.Calculation of relative displacement during sticking Assuming, a material point P on a segment s i of a beam element being in contact with the circle, see Fig. 4, we want to track the relative position of P and circle's local frame.The reference position of P on the segment is x s, p = ρ L seg , in which ρ is the one defined in 2.3, with L seg being the segment's length.For defining the relative position of P to the circle, we use the arc length Therefore, we have no relative motion if: Fig. 4 The current sticking position results from the sum of the position of P relative to the beginning of the segment and the circle's frame Defining the current sticking position x cur Stick as the displacement of contact point while sticking Δx stick reads: in which the function floor() is a standardized version of rounding, available in C and Python programming languages.In (25), x stick Re f is equal to x cur Stick when sticking first starts, see Fig. 3, and can be thought as the reference length of the virtual spring which constrains the tangential motion.The algorithm for updating x cur Stick and x stick Re f is described in the next paragraph.Note that the vertical offset from the beams' centreline (resulting by the nonzero width of the beam), as well as the axial stretch are not considered in the calculation of (un-)winding of the segment.While a reduction in segment length reduces this effect because less (un-)winding occurs, a more consistent computation of this effect would require an integration of relative motion, as stretch slightly influences the local changes of the relative sticking position.
Post Newton step As mentioned, the implemented tangential model considers two possible states, sticking and slipping, that lead to different computation formulas for the tangential force.Although the computation of tangential forces occurs during the Newton method of the implicit integrator the switching from the one state to the other is not allowed inside the Newton method.Thus, we introduce a Post-Newton step (PNS()) which computes a set of data variables: In this data set x gap = g as defined in (16), and for x sli p , there are three cases: a. x sli p = −2: undefined, used for initialization b. x sli p = 0: sticking c. x sli p = ±1: slipping, sign defines slipping direction For initializing this data set, we are using [0, −2, 0].An error (in N) is calculated during every PNS(): e PNS = e n PNS + e t PNS (27) in which e n PNS corresponds to the calculation of the normal contact force and e t PNS to the calculation of the tangential contact force.For computing e n PNS , if the gap of the previous PNS(), x gap,last PNS , had a different sign as the current gap, set while otherwise e n PNS = 0.For computing e t PNS , if stick-slip-state of previous PNS(), x sli p,last P N S , is different from current x sli p , set while otherwise e t PNS = 0.The PNS(), with all operations performed for every segment, s i , is summarized in Fig. 5.

Contact forces computation in the Newton step
The calculation of contact occurs during Newton iterations, and it is performed only if the PNS() determined contact in any segment, for efficiency.For a segment s i , if x gap,s i ≤ 0, the steps are followed: 1. We compute the contact force, f n , according to (9). 2. If x sli p = 1, -and if the friction stiffness μ k = 0 or if x sli p = −2, we set Δx stick = 0, -else, we compute the current sticking position, x cur Stick according to (24) and following the Δx stick from (25) using as x stick Re f the updated data variable of the last PNS().-using the tangential velocity from (19), the tangential force is computed by (20).
3. If x sli p = 1, the tangential force is set as Comparison of proposed tangential contact model and LuGre model To compare the agreement of the proposed tangential contact model with LuGre model, we reproduce a numerical example described in [24].To briefly describe the experimental set up of this example, a mass, M = 1 kg, is connected to a spring, of stiffness K = 2 N m −1 , the tip of which has a constant velocity, v = 0.1 m s −1 .The LuGre model, described in [24], is characterized by six parameters, the bending stiffness of the bristle σ 0 = 10 5 N m −1 , the bending damping of the bristle σ 1 = 10 5/2 N s m −1 , the viscous stiffness σ 2 = 0.4 N s m −1 , the Coulomb friction level F c = 1 N, the level of the stiction force F s = 1.5 N and the Stribeck velocity v s = 0.001 m s −1 .We reproduce the example with the proposed tangential contact model using μ k = 10 5 and μ v = 10 5/2 .We set viscous stiffness to zero in both models.We measure the displacement and the velocity of the mass as well as the friction force and we compare the results obtained by the LuGre and the proposed model, see Figs. 6 and 7.While the LuGre model converges to the proposed model for increasing bending stiffness of the bristle, the proposed model manifests its superior performance as it reproduces the expected friction force without the need for an internal parameter.Note that in the LuGre model the shifting from sticking to sliding depends on the internal dynamic parameter σ 0 and converges to the level of the stiction force F s = 1.5 N for increasing σ 0 , see Fig. 7.

Mapping contact forces on the beam elements and sheaves
Contact forces f i with i ∈ [0, n cs ] are applied at the points p i , and they are computed for every contact segment.For every contact computation, first all contact forces at segment points are set to zero.If there is contact in a segment s i , e.g., gap state x gap ≤ 0, see Fig. 8, contact forces f s i are computed per segment, and added to every force f i , initialized with zero, at segment points according to while in case x gap > 0 nothing is added.Note that for the computation of the normal and tangent vector a second option which defines the normal as a vector pointing from the center of the circle to the segment point has been investigated.In case of segments that are short as compared to the circle radius, the second option is more consistent and produces tangents only in circumferential direction, which improves accuracy for coarse meshes.However, the first option, which uses n s i , t s i , vectors vertical and tangential to the segment, see Fig. 8, lead to always good approxi-mations for normal directions, irrespective of short or extremely long segments as compared to the circle.The forces f i are then applied to the ANCF cable element at according segment points.The forces on the circle are computed as the total sum of all segment contact forces, and additional torques on the circle's rotation simply follow from  connected to the ground.We apply a lateral force f = 200 N at the center of the circle.In Fig. 10, we plot the absolute error for the y−coordinate of the position of the center of the circle and the axial force measured at the middle of beam.The error is calculated as the absolute difference from a reference solution calculated for 128 elements and 12 linear segments.As the numbers of elements and linear segments increase, we observe convergence.

Numerical experiment
Description of set-up The presented methods are applied in the numerical example of a belt drive consisting of an elastic belt and two identical pulleys, P 1 and P 2 , see Fig. 11.Although the overall set up of the system is similar to [13], some key-parameters have Fig. 11 Belt drive with two pulleys been changed for making the system suitable for convergence analysis and capable to reach a steady state motion.Therefore, we find it necessary to describe the system here.The belt is modeled by cubic ANCF elements described in the previous section.The pulleys are modeled as rigid bodies mounted to the ground through revolute joints.Dimensions and properties of belt and pulleys are given in Table 2.The initial as well as the reference length of the belt centreline is the length that results from the geometry of the system: The stress-free length of the belt centreline is chosen as lb = 0.95•l b , resulting in a pre-stretch of ε re f = −0.05,which is added to all finite elements before the static  computation for initial conditions and thus, results in a pre-tension of While the reference configuration of all beams is straight, the initial values are computed from a static computation, in order to avoid vibrations at the beginning of the dynamic simulation.Dynamic simulation During the dynamic simulation, the angular velocity of P 1 is prescribed, see Fig. 12.In order to reach a steady state, a torque is added to P 2 in the time range t = [1, 1.5] seconds: Fig. 12 Angular velocity of P 2 for varying number of elements and comparison with angular velocity of P 1 The example was simulated in the multibody dynamics code Exudyn [31].For the time integration, an implicit second-order time integration method (trape- zoidal rule) is used, similar to Newmark's method [32], but without numerical damping and with an index 2 constraint reduction, which allows stable integration of constraints.The discontinuous iteration for the PNS() uses a tolerance of 10 −3 and a maximum of 5 iterations, which is not reached in case of smaller time steps.The belt drive has been computed with different parameter sets (e.g., different number of elements), starting from a nominal parameter set as well as a reference parameter set, which uses finer space and time discretisation as compared to the nominal test case, see Table 3.

Results
The results' section is split into three parts: -results of angular velocities and torques over time, which reflect the overall behavior of the belt drive, -results of axial forces and contact stresses over the belt arc length, s l , see Fig. 11.Note that the quantities shown w.r.t. the belt length are given in terms of the curve parameter sl, see Fig. 11, and -comparisons with analytical solutions for the belt drive.

Results over simulation time
The overall behavior of the belt can be seen in Fig. 12 which shows the angular velocities of both pulleys, noting that the velocity of P 1 is prescribed and thus, it is the same for all variations.The angular velocity of pulley P 1 speeds up at 0.05s, while due to belt elasticity, slip and inertia effects, the angular velocity of pulley P 2 shows vibrations and a smaller velocity in steady state.We observe almost no effects for varying number of elements for over 60 elements.
Fig. 13 Torque at P 1 for varying number of elements and comparison with torque load at P 2 (not including torque due to damping) The torque in pulley P 1 , which is highly influenced by the torque in pulley P 2 , is shown in Fig. 13 for different numbers of elements.Note that there is an angular velocity-proportional (damping) d P2 acting on pulley P 2 as well.It is also clearly seen that vibrations in the P 2 angular velocity are merely damped out until t = 1s, while there is an additional drop hereafter due to torque τ P2 .Torques in Fig. 13 show some oscillatory behavior for the coarse discretization with 60 elements.Results along the belt arc length We compare and evaluate the axial velocity, the tension, tangential contact stresses, normal contact stresses along the belt arc length varying the number of finite elements as well as varying the quantities n cs , μ, k c and μ k .In Figs.14,15,16,17,18,19,20,21,22,23 one can distinguish four different zones for the belt arc length:  Figures 14,15,16,17 show the convergence of the proposed contact algorithms w.r.t.number of elements.This number significantly influences the axial velocity, see Fig. 16, and the tangential contact stresses, see Fig. 17.Compared to that, the normal contact stresses only show oscillations for 60 elements, see Fig. 14, and only minor influences of the number of elements are seen in Fig. 15.
We evaluate the influence of the normal contact stiffness k c (which is coupled to normal contact damping, d c ), the tangential contact stiffness μ k , the number of segments n s and the friction coefficient μ for otherwise nominal parameters as shown in Figs. 18, 19, 20, 21.We observe that increasing the number of linear segments from 4 to 8 has almost no effect.However, using fewer segments would worsen the geometric approximation of the circle for small number of finite elements and was therefore avoided.In Fig. 19, the friction coefficient shows a high influence on the exponential force drop in the contact slipping zones of the pulleys as known from the belt drive theory.
In Fig. 22, we plot the axial forces at t = 1 s for the different integration schemes explained in 2.1.We observe oscillatory behavior for integration schemes a and b, while the solution obtained using integration scheme c is much smoother in all regions of belt drive.Comparison with analytical solution For comparing the numerical results with available analytical values, we use the solution at time t = 2.45 s, at which the system has reached the steady state as shown in Fig. 23.As known, from classic belt drive literature [33], if F 1 is the axial force on the tight side and F 2 is the axial force on the slack side, see Table 4, the drop of span forces while the creep computed from the relative elongation of the tight and slack spans is Both values are relative values, but they agree well with 3.3% error.We plot the axial forces in the contact zone of P 1 as resulted from the simulation, and we compare it with the analytical solution calculated by Euler-Eytelwein also known as capstan formula, [33]: (42) for φ ∈ [0, β] in which with corresponding region on the arc length of the belt see Fig. 24.The maximum absolute error is 13.60 N, and the maximum relative error is equal to 2.5%.

Conclusions
In this paper, a model for contact and friction for flexible beams and sheaves and a numerical example have been provided.We have shown that the proposed contact model consisting of an improved ANCF implementation with a selective reduced integration, a contact detection based on the subdivision of beam elements into linear segments, a bounding box for effective contact detection and a bristle friction model with history variables can be successfully applied for a numerical example of a belt drive, as the one demonstrated here.The example of the belt drive showed convergence of the implemented methods for increasing number of elements and for other varying parameters such as the number of segments and the dry friction.Moreover, we saw agreement of the obtained axial forces in the contact zone of the first pulley of the belt drive with the analytical solution from classical belt theory.Finally, the resulted beam axial forces along the belt arc length showed improved behavior for the proposed integration scheme for the virtual work of elastic and viscous damping forces.As a consequence, we conclude that we can use the developed methods for an efficient contact modeling for systems consisting of flexible beams and sheaves such as belt drives and reeving systems.

Fig. 5 Fig. 6 Fig. 7 Fig. 8
Fig.5 Flowchart for PNS().Returns the set of data variables [x gap , x sli p , x stick Re f ] and e P N S

) 3
Convergence in the static friction-less caseWe examine the example of a circle of radius r = 0.35 m moving toward a flexible beam of Young's modulus E = 2 10 8 N m 2 , length L = 1 m, density

Fig. 9 9 .
Fig. 9 Quasi-static experiment of circle in contact with beam

Fig. 10
Fig. 10 Convergence of error for the y-coordinate of the position of the center of the circle, on the left, and for the axial force at the middle of the length of the beam, on the right

Fig. 14 Fig. 15 Fig. 16
Fig.14 Normal contact stresses along the length of the belt for time t = 2.45 s for varying number of elements

Fig. 18 Fig. 19 Fig. 20 Fig. 21 Fig. 22 Fig. 23
Fig.18 Normal contact stresses along the length of the belt for time t = 2.45 s for varying parameters

Table 1
Implemented integration rules in Exudyn

Table 4
24proximated numerical values at t = 2.45 s for overall validation of contact algorithm Note that some of the measured values are not constant along the free spans and thus are approximated and rounded.Furthermore, we only consider absolute values here to simplify computations Fig.24Comparison of axial forces in the contact zone of P 1 and analytical solution is related to torques τ P1 and τ P2 by (F 1 − F 2 )r = 47.857Nm≈τP1=47.68Nm(38)whichresults to a relative error equal to 0.37%.Note, also, that τ P2 at time t = 2.45 s results from the velocity-proportional damping d P2 and the added load of 25 N m given by (37):τ P2 = 25 + d P2 • ω P2 = 47.68 N m , (39) which agrees with τ P1 , as measured at t = 2.45 s.The creep computed by the velocities of the spans is