Front. Signal Process. Frontiers in Signal Processing Front. Signal Process. 2673-8198 Frontiers Media S.A. 1525044 10.3389/frsip.2025.1525044 Signal Processing Original Research Numerical modelling of elasto-plastic friction in bow–string interaction with guaranteed passivity Matusiak et al. 10.3389/frsip.2025.1525044 Matusiak Ewa 1 * Chatziioannou Vasileios 1 Van Walstijn Maarten 2 1 Department of Music Acoustics (IWK), University of Music and Performing Arts Vienna, Vienna, Austria 2 Centre for Interdisciplinary Research in Music and Sound, School of Electronics, Electrical Engineering, and Computer Science, Queen’s University Belfast, Belfast, United Kingdom

Edited by: Augusto Sarti, Polytechnic University of Milan, Italy

Reviewed by: Sorin Vlase, Transilvania University of Brașov, Romania

Silviu Nastac, Dunarea de Jos University, Romania

 *Correspondence: Ewa Matusiak, matusiak@mdw.ac.at 21 05 2025 2025 5 1525044 08 11 2024 11 03 2025 Copyright © 2025 Matusiak, Chatziioannou and Van Walstijn. 2025 Matusiak, Chatziioannou and Van Walstijn

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

In order to simulate the function of bowed string instruments, it is necessary to model the frictional interaction between the bow hair and the vibrating string. This is possible using an elasto-plastic friction model, which has previously succeeded in reproducing experimental data captured on a monochord setup. In this study, this elasto-plastic model is refined to guarantee passivity, and a stable numerical scheme is derived that inherits the energy balance of the underlying continuous model. The approach presented considers a finite-width bow, thus spreading the bow–string interaction over an area. The compliance of the bow hair and the torsional motion of the string are also taken into account. A sound example and animations of the string motion are provided to demonstrate the behavior of the model.

 bowed mass bow–string interaction friction passivity finite differences energy methods numerical stability elasto-plastic P34852 Austrian Science Fund10.13039/501100002428 section-at-acceptance Audio and Acoustic Signal Processing 1 Introduction

In order to simulate bow–string interaction, it is crucial to accurately model the friction between string and bow. Several friction models—both static and dynamic—have been developed in recent decades. The static models were obtained after measuring the coefficient of friction either in a steady-state (see friction curve model suggested by Smith and Woodhouse (1999)) or in a transient part of the waveform (see friction curve model suggested by Galluzzo, 2004). Among existing dynamic friction models used to simulate string vibrations are an elasto-plastic friction model developed by Dupont et al. (2002) and first applied in bow–string simulations by Serafin et al. (2003) and a thermal model introduced by Woodhouse (2003) where the temperature of rosin is considered; it is implemented using digital waveguides in Maestre et al. (2014). Regardless of the modeling choice, it is essential to use a guaranteed-passive model along with a discretization method that preserves the energy-conservation properties of the continuous system (e.g., Chabassier and Joly, 2010; Bilbao et al., 2015; Desvages, 2018; Ducceschi and Bilbao, 2022; van Walstijn et al., 2024).

This study examines the application of the elasto-plastic friction model by Dupont et al. (2000) and Dupont et al. (2002) to numerical modeling of bow–string interaction. Their formulation is a refinement of the LuGre friction model (de Wit et al., 2024) which encountered drift at low sliding velocities. The idea behind the model is that two sliding surfaces are irregular at the microscopic level; their interaction is modelled as a bundle of elastic bristles, with each bristle contributing to the overall frictional force. This model has already been applied to bowed strings: it was implemented using a finite difference method in Willemsen et al. (2019), where it was applied to point-bowing a stiff string, and in Matusiak and Chatziioannou (2024), where it was applied to a finite-width bow model. A comparison between this elasto-plastic model and the thermal friction model introduced in Smith and Woodhouse (1999) was conducted in Serafin (2004) using both a digital waveguide implementation and a finite difference method locally under the bow. That study focused on highlighting the differences and similarities between these two dynamic friction models. It was demonstrated in Matusiak and Chatziioannou (2024) that, given the right set of parameters, the elasto-plastic friction model is able to reconstruct the steady-state as well as the transient of a measured waveform.

Although the elasto-plastic friction model has been successfully implemented in the above studies, the passivity of the model has not been thoroughly investigated. This is important in order to guarantee the stability of numerical simulations. In a recent study by Falaize and Roze (2024) on nonlinear interaction modelling, the Dupont model was incorporated in a Port-Hamiltonian formulation, which—besides the part related to the elastic bristles—can be shown to be passive. Passivity, however, is not guaranteed for the coupled bow–string interaction model, since the dissipation term associated with the bristle displacement can become negative for certain parameter values.

In this study, we re-examine this issue and propose a refined model which is shown to be passive for any model parameters. Such a refinement can also demonstrate the existence and uniqueness of solutions to the underlying differential equations. Furthermore, in order to numerically implement the model in a stable manner, an energy preserving discretization scheme is derived. This is an improvement on the implementation proposed in Willemsen (2021), where the numerical bristle energy was not guaranteed to be non-negative.

The paper is structured as follows. In Section 2, the elasto-plastic model is introduced in the context of a lumped bowed mass, and energy analysis is performed in the continuous domain to reveal the lack of passivity of the original elasto-plastic model. Section 3 presents a refined version of the model that preserves passivity and guarantees the uniqueness of the solution. Section 4 is concerned with the numerical formulation of the bowed mass model; discretization using the finite difference method and energy analysis in the discrete setting are performed, followed by numerical experiments and comparison of the two models. In Section 5, the friction model is applied to the problem of bowing a string with a bow of finite width. Transverse and torsional waves on the string are accounted for. Section 6 presents the numerical model for the bow–string interaction, including analysis of the discrete energy, and some concluding remarks are given in Section 7.

 2 Elasto-plastic friction model

It is helpful to first present the problem in its simpler form, which involves modelling the bowing of a lumped mass connected to a stiff spring. This represents a specific case of the bowed string model, restricted solely to the transverse movement of a string that is bowed at a single contact point and considering only its fundamental mode of vibration (see Appendix).

Consider a bowed mass m undergoing a tangential friction force F (Figure 1). The mass is excited by a bow moving with velocity v b which is modeled as a harmonic oscillator with bow hair mass m h , stiffness K h [kg/s2], and damping Γ h [kg/s]. The bow hair displacement relative to the rigid bow is denoted by η = η ( t ) . The motion of the two masses can be then described by m u ̈ = − κ u − γ u ̇ − F , m h η ̈ = − K h η − Γ h η ̇ − F , v = u ̇ − v b − η ̇ , where u = u ( t ) denotes the displacement of the mass attached to a spring of stiffness κ [kg/s2] and damping γ [kg/s]. The relative velocity v between bow and mass is given by Equation 3.

Sketch of the bowed lumped mass system.

Following Dupont et al. (2002), an elasto-plastic model is used to simulate the tangential friction force. This model assumes that the two surfaces—in this case the bow hair and the mass—are irregular at the microscopic level, and their contact is modelled through an ensemble of elastic bristles, each contributing to the total friction load. The bristles are modelled as damped stiff springs, and when the strain exceeds a certain breakaway threshold, the bristles break, and the two surfaces begin to slide. Denoting by z , the average bristle deflection, and by v , the relative velocity between the string and the bow, the model is described as follows. F = σ 0 z + σ 1 z ̇ , z ̇ = v 1 − α z , v z z ss v , where σ 0 in [kg/s2] and σ 1 in [kg/s] are bristle stiffness and damping, respectively. The change of rate with which the bristles stretch or contract is related to the relative velocity through the adhesion map α , defined as α z , v = 0 , v z ≤ 0 α m z , v , v z > 0 , where α m z , v = 0 , | z | ≤ z ba α ̄ m v , z , z ba < | z | < | z ss v | 1 , | z | ≥ | z ss v | , and where z ba ≤ μ C f N / σ 0 is the breakaway displacement α ̄ m z , v = 1 2 sin π θ z , v + 1 , θ z , v = | z | − 1 2 | z ss v | + z ba | z ss v | − z ba , and z ss is the steady-state displacement for constant velocities z ss v = f N σ 0 μ C + μ S − μ C e − | v / v S | p , v ≥ 0 − f N σ 0 μ C + μ S − μ C e − | v / v S | p , v < 0 , with v S > 0 and p ≥ 1 . When p = 2 , v S is referred to as “Stribeck velocity.” In the numerical experiments in Sections 4.4 and 6.4, p = 2 is used. Here, μ S and μ C denote static and dynamic friction coefficients, respectively, and f N is the normal force applied by the bow. For small bristle displacements, when | z | ≤ z ba , α ( v , z ) = 0 , and consequently z ̇ = v , a purely elastic and reversible regime is entered, referred to as “pre-sliding” (sticking). For larger displacements—that is, when z ba < | z | < | z ss ( v ) | —some bristles start to break, and a mixed elasto-plastic sliding occurs. Finally, for | z | ≥ | z ss ( v ) | , all bristles break, and a purely plastic regime is achieved—the string slips under the bow. In that situation, α ( z , v ) = 1 . This model for a friction force was first developed in Dupont et al. (2000) and used for simulating friction in various industrial applications.

2.1 Energy balance

In order to investigate the passivity of the model given by Equations 1–5, the energy balance of the system is considered.

Multiplying Equation 1 with u ̇ , Equation 2 with η ̇ , and Equation 5 with σ 0 z , summing up and using the relation in Equation 3 yields the energy balance m u ̈ u ̇ + κ u u ̇ ⏟ H ̇ r + m h η ̈ η ̇ + K h η η ̇ ⏟ H ̇ h + σ 0 z z ̇ ⏟ H ̇ b = − v b F ⏟ P − γ u ̇ 2 ⏟ Q r − Γ h η ̇ 2 ⏟ Q h − σ 1 v 2 + α z , v v z z ss v σ 0 z − σ 1 v ⏟ Q b , where H r = m 2 u ̇ 2 + κ 2 u 2 ≥ 0   i s t h e o s c i l l a t o r k i n e t i c & p o t e n t i a l e n e r g y , H h = m h 2 η ̇ 2 + K h 2 η 2 ≥ 0   i s t h e b o w h a i r k i n e t i c & p o t e n t i a l e n e r g y a n d H b = σ 0 2 z 2 ≥ 0 i s t h e e n e r g y s t o r e d i n t h e b r i s t l e s . P is the power supplied by the bow via the friction force, and Q r , Q h , and Q b are dissipation terms, corresponding to the oscillator, the bow hair, and the bristles, respectively. This leads to the following conservation law: H r + H h + H b + ∫ P + Q r + Q h + Q b d t = H 0 , where H ( 0 ) is the initial system energy (if any). Q r and Q h are trivially non-negative, and Q b is non-negative whenever α ≠ 1 . For α = 1 , and for high—relative to σ 0 —values of σ 1 , Q b can become negative, which violates passivity and is unphysical (Dupont et al., 2000). Therefore, without a certain condition on the relationship between the stiffness and damping coefficients, passivity cannot be ensured. Following Olsson (1996), who analyzed a simpler version of the Dupont model—the so-called LuGre friction model (de Wit et al., 2024)—Willemsen (2021) derived a condition on σ 1 to guarantee passivity for the elasto-plastic friction model. The condition reads: σ 1 ≤ 4 σ 0 z ss v | v | .

However, without knowing the maximal relative velocity, v , it is not possible to set the value for σ 1 . In Section 3, a different condition on the bristle damping term is proposed that is less restrictive and does not require any knowledge of the limits of v .

Falaize and Roze (2024) employed the elasto-plastic friction model in the framework of port-Hamiltonian systems and arrived at the same dissipation term— Q b . For certain parameter choices, the dissipation matrix defined in Falaize and Roze (2024) is semi positive-definite, but for that property to generally hold, a refinement of σ 1 is needed.

 2.2 Boundedness of the bristle displacement

Boundedness of bristle displacement was first shown in Olsson (1996) for a LuGre friction model, and subsequently in Dupont et al. (2000) for an elasto-plastic model, by defining a positive definite Lyapunov function and an invariant set of solutions. Here a slightly different approach is presented.

The bristle displacement z ( t ) can be thought of as a parametric curved defined by z ( t ) and v ( t ) whose rate of change is z ̇ t = v t 1 − α z t , v t z t z ss v t , and sign z ̇ t = − 1 , z t > z ss v t 0 , z t = z ss v t or v t = 0 . 1 , z t < z ss v t This means that for v ( t ) > 0 , if z ( t ) reaches the maximal value of z ss ( v ( t ) ) which is z ss ( 0 ) = μ S f N σ 0 , it cannot rise any further but has to either decrease or stay constant. Similarly, for v ( t ) < 0 , if z ( t ) reaches the minimal value of z ss ( v ( t ) ) , which equals lim v → 0 − z ss ( v ( t ) ) = − μ S f N σ 0 , it cannot decrease any further but has to either rise or stay constant. Therefore, | z ( t ) | ≤ μ S f N σ 0 (see Figure 2 for visualization). The boundedness of z by itself does not, however, imply stability.

Behavior of the bristle displacement z ( t ) . Left: pink color indicates regions where z ( t ) is decreasing ( sign ( z ̇ ( t ) ) = − 1 ) and yellow is where z ( t ) is increasing ( sign ( z ̇ ( t ) ) = 1 ) . Blue line is the region where z ̇ ( t ) = 0 , and the black dashed line denotes the maximum and minimum value of z ss . Right: Example trajectory in the v − z plane for model parameters as in Table 1.

 3 Refined elasto-plastic model

This section presents a refined version of the elasto-plastic friction model that addresses the passivity of the system. Consider a mass m bowed with velocity v b as described in Section 2, by Equations 1–5, where bristle damping σ 1 is not constant but varies with v and f N as follows: σ 1 v = μ C f N v 2 + ϵ 2 , ϵ = μ C f N σ ̄ 1 , such that 0 ≤ σ 1 ( v ) ≤ σ ̄ 1 and σ 1 ( 0 ) = σ ̄ 1 .

A similar velocity-dependent damping term σ 1 ( v ) = σ ̄ 1 e − ( v / v d 2 ) was introduced by Olsson (1996) for the LuGre model, motivated by the need to reproduce certain friction phenomena. Provided that the additional free parameter v d chosen is sufficiently small, the condition in Equation 11 is satisfied, in turn guaranteeing passivity of the LuGre model. This dependency on an external parameter is avoided in the refinement proposed here (Equation 14). Furthermore, for large v values, σ 1 stays closer to the original constant σ ̄ 1 than with the exponential formula of Olsson.

3.1 Energy balance

With σ 1 now being a function of the relative velocity v , the passivity of the system in Equations 1–5 is guaranteed.

Proposition 1: Let F be a friction force as defined in Equation 4 with bristle damping σ 1 defined in Equation 14. Then, the system in Equations 1–5 is passive.

Proof: All terms in the energy balance (Section 2.1) are trivially non-negative, apart from Q b . For α ( v , z ) < 1 , Q b is always non-negative. Therefore, only the case when α ( z , v ) = 1 , which happens when sign ( z ) = sign ( v ) and | z | ≥ | z ss ( v ) | , is considered here. In that case, we have Q b = σ 1 v v 2 + | z | | v | | z ss v | σ 0 | z | − σ 1 v | v | ≥ σ 1 v v 2 + | z | | v | | z ss v | σ 0 | z ss v | − σ 1 v | v | ≥ σ 1 v v 2 + | z | | v | | z ss v | μ C f N − σ 1 v | v | = σ 1 v v 2 + | z | | v | | z ss v | μ C f N 1 − | v | v 2 + ϵ 2 > 0 .  

3.2 Existence and uniqueness of the solution

The passivity of the system can be shown to guarantee the existence and uniqueness of the solution. By introducing a variable ϕ , the bowed mass system of Equations 1–5 together with Equation 14 can be written as an autonomous system of equations. u ̇ = ϕ , η ̇ = v + v b − ϕ , ϕ ̇ = 1 m − κ u − γ ϕ − σ 0 z − σ 1 v v 1 − α z , v z z ss v , v ̇ = − κ u − γ ϕ m + − K h η − Γ h v + v b − ϕ m h − m + m h m ⋅ m h σ 0 z + σ 1 v v 1 − α z , v z z ss v , z ̇ = v 1 − α z , v z z ss v .

Let x = ( u , η , ϕ , v , z ) , then the above can be written as x ̇ = S x , where S = ( s u , s η , s ϕ , s v , s z ) with the functions s corresponding to Equations 15–19. Given initial conditions x ( t 0 ) = ( u 0 , η 0 , ϕ 0 , v 0 , z 0 ) , existence and uniqueness of the solution to Equation 20 is guaranteed whenever S ( x ) is Lipschitz continuous (Barreira and Valls, 2012).

A vector-valued function S ( x ) is Lipschitz continuous if there exists an L ≥ 0 , called the “Lipschitz constant,” such that for all x , y in the domain of S | S x − S y | ≤ L | x − y | .

If S ( x ) = ( s 1 ( x ) , … , s K ( x ) ) , then the norm of S ( x ) is defined as | S x | 2 = ∑ i = 1 K | s i x | 2 . The norm of x is similarly defined.

Proposition 2: Consider a bowed mass m described by Equations 1–5 with bristle damping given by Equation 14. Then, for given initial conditions, there exists a solution to the bowed mass system that is unique.

Proof: By the Picard–Lindelöf theorem (Barreira and Valls, 2012), a system of ordinary differential equations S has a global unique solution if S is Lipschitz continuous with respect to x with a Lipschitz constant not depending on x . This is equivalent to all the functions s u , … , s z being Lipschitz continuous with respect to each variable u , η , ϕ , v , and z .

By passivity, all the variables are bounded. Lipschitz continuity of s u and s η is trivial. Similarly, Lipschitz continuity with respect to u , η , and ϕ is trivially satisfied for all the functions s u , s η , s ϕ , s v , s z . Moreover, it is straightforward to verify that σ 1 ( v ) defined in Equation 14 is Lipschitz continuous. It remains to show that g ̃ z , v = α z , v z v z ss v is Lipschitz continuous.

(i) Let v be fixed. Then, for z 1 , z 2 , g ̃ z 1 , v − g ̃ z 2 , v ≤ | v | | z ss v | | α z 1 , v | | z 1 − z 2 | + | z 2 | | α z 1 , v −   α z 2 , v | ≤ | v | | z ss v | | z 1 − z 2 | + π 2 μ S f N σ 0 θ z 1 , v − θ z 2 , v ,

where the last inequality follows from the Lipschitz continuity of the sin function with Lipschitz constant equal to 1, with | sin ( π 2 θ ( z 1 , v ) ) − sin ( π 2 θ ( z 2 , v ) ) | ≤ π 2 | θ ( z 1 , v ) − θ ( z 2 , v ) | . The function θ is also Lipschitz continuous with respect to z , θ z 1 , v − θ z 2 , v = | | z 1 | − | z 2 | | | z ss v | − z ba ≤ | z 1 − z 2 | | z ss v | − z ba ≤ σ 0 μ C f N − σ 0 z ba | z 1 − z 2 | . Let v be bounded by B v , then g ̃ z 1 , v − g ̃ z 2 , v ≤ B v σ 0 μ C f N + π B v σ 0 μ S μ C μ C f N − σ 0 z ba ⏟ L z | z 1 − z 2 | . Hence, g ̃ ( z , v ) is Lipschitz continuous with respect to z with Lipschitz constant L z .

(ii) Now, let z be fixed and z > 0 . Then g ̃ z , v = 0 , v ≤ 0 α z , v z v z ss v , v ≥ 0

is continuous and everywhere differentiable except for v = 0 , with ∂ g ̃ ∂ v = 0 , v < 0 ∂ α ∂ v z v z ss + α z z ss − z v z ̇ ss z ss 2 , v > 0 , and ∂ α ∂ v = π 2 cos π θ z ̇ ss z ba − | z | z ss − z ba 2 . The partial derivative of g ̃ with respect to v is bounded and positive, since all the elements are bounded and positive. Therefore, by the mean value theorem for v 1 , v 2 > 0 g ̃ z , v 1 − g ̃ z , v 2 ≤ L v 1 | v 1 − v 2 | , L v 1 = max ∂ g ̃ ∂ v , and g ̃ is Lipschitz continuous for v > 0 , with Lipschitz constant L v ( 1 ) , and trivially for v < 0 , with any L v ( 2 ) ≥ 0 as a Lipschitz constant. Now let v 1 > 0 and v 2 ≤ 0 , then g ̃ z , v 1 − g ̃ z , v 2 | v 1 − v 2 | = g ̃ z , v 1 | v 1 − v 2 | ≤ g ̃ z , v 1 | v 1 | = | z | | z ss v 1 | ≤ μ S μ C . Let L v ( 3 ) = μ S / μ C . Hence, g ̃ is Lipschitz continuous with respect to v with Lipschitz constant L v = max { L v ( 1 ) , L v ( 2 ) , L v ( 3 ) } . A similar proof holds for z < 0 .

Therefore, S is Lipschitz continuous with respect to x , and a unique solution to Equation 20 is guaranteed.

The original Dupont model was not shown to be passive, as there was no guarantee that v is bounded. In that case, the system in Equations 1–5 has a unique solution only locally.

 4 Numerical formulation

In the present section, a finite difference numerical scheme is utilized to discretize the system in Equations 1–5. The proposed discretization is energy conserving.

4.1 Numerical preliminaries

The approximations to u ( t ) at points n Δ t are denoted as u n , where Δ t is the time step. The variable z is approximated at an interleaved grid—that is, z n + 1 2 denotes the approximation to z ( t ) at time t = ( n + 1 2 ) Δ t . The following centered, second-order accurate discretization operators are defined: δ t ○ u n = u n + 1 2 − u n − 1 2 Δ t = u ̇ n Δ t + O Δ t 2 , δ t ⋅ u n = u n + 1 − u n − 1 2 Δ t = u ̇ n Δ t + O Δ t 2 , μ t ○ u n = u n + 1 2 + u n − 1 2 2 = u n Δ t + O Δ t 2 , μ t ⋅ u n = u n + 1 + u n − 1 2 = u n Δ t + O Δ t 2 .

In addition, one may define non-centered, first-order accurate discretization operators δ t − u n = u n − u n − 1 Δ t = u ̇ n Δ t + O Δ t , δ t + u n = u n + 1 − u n Δ t = u ̇ n Δ t + O Δ t , μ t − u n = u n + u n − 1 2 = u n Δ t + O Δ t , μ t + u n = u n + 1 + u n 2 = u n Δ t + O Δ t .

The composition of first-order accurate operators results in second-order accurate operators, δ t t u n = δ t + δ t − u n = u n + 1 − 2 u n + u n − 1 Δ t 2 = u ̈ n Δ t + O Δ t 2 , μ t t u n = μ t + μ t − u n = u n + 1 + 2 u n + u n − 1 4 = u n Δ t + O Δ t 2 ,

and several useful identities can be constructed, including δ t ⋅ u n δ t t u n = δ t + 1 2 δ t − u n 2 , δ t ⋅ u n μ t t u n = δ t + 1 2 μ t − u n 2 , δ t ⋅ u n u n = δ t + 1 2 u n u n − 1 , δ t ○ z n μ t ○ z n = δ t + 1 2 ( z n − 1 2 2 .

4.2 Discretization

For simplicity of notation, let g z , v = v 1 − α z , v z z ss v . The system of Equations 1–5 is discretized as follows: m δ t t u n = − κ u n − γ δ t ⋅ u n − F n , m h δ t t η n = − K h μ t t η n − Γ h δ t ⋅ η n − F n , v n = δ t ⋅ u n − v b n + δ t ⋅ η n , F n = σ 0 μ t ○ z n + σ 1 v n δ t ○ z n , δ t ○ z n = g μ t ○ z n , v n .

In the case of the original Dupont model, σ 1 ( v n ) = σ ̄ 1 . In order to ensure numerical stability for the lumped mass model (Equation 25), the following stability condition is obtained using frequency domain analysis (Bilbao, 2009): Δ t < 2 m / κ 1 / 2 .

This discretization choice for the equation of motion of the lumped mass is motivated by the discretization that is carried out for the string model in the distributed case in Section 6. On the other hand, the stiffness term for the bow hair in Equation 26 is discretized using an averaging operator in order to avoid introducing a further stability condition.

For simpler notation, let r n = μ t ○ z n . Using the discretization operators, one can establish identities δ t t u n = 2 Δ t δ t ⋅ u n − δ t − u n , δ t ○ z n = 2 Δ t r n − z n − 1 2 , μ t t u n = Δ t 2 δ t ⋅ u n − δ t − u n + u n .

Substituting into Equations 25, 26, and 28, one obtains m 2 Δ t + γ δ t ⋅ u n = − κ u n + m 2 Δ t δ t − u n − σ 0 r n − σ 1 v n 2 Δ t r n − z n − 1 2 , m h 2 Δ t + K h Δ t 2 + Γ h δ t ⋅ η n = − K h η n + m h 2 Δ t + K h Δ t 2 δ t − η n − σ 0 r n − σ 1 v n 2 Δ t r n − z n − 1 2 . Considering u n , u n − 1 , η n , η n − 1 , and z n − 1 2 to be known, let p n = δ t − u n and s n = δ t − η n . Using Equation 27, r n can be expressed as a function of v n r n = − A v n v n + B v n , where A v n = m 2 Δ t + γ m h 2 Δ t + K h Δ t 2 + Γ h 2 Δ t m + m h + K h Δ t 2 + γ + Γ h σ 0 + 2 Δ t σ 1 v n − 1 , B v n = A v n − κ u n + m 2 Δ t p n m 2 Δ t + γ + − K h η n + m h 2 Δ t + K h Δ t 2 s n m h 2 Δ t + K h Δ t 2 + Γ h − v b n + 2 Δ t σ 1 v n z n − 1 2 σ 0 + 2 Δ t σ 1 v n . For the original Dupont model, A ( v n ) and B ( v n ) are independent of v n , and r n is linearly related to v n . Substituting Equation 31 into Equation 29 and using the middle identity in Equation 30 yields a nonlinear equation in the unknown v n , 2 Δ t B v n − A v n v n − z n − 1 2 = g B v n − A v n v n , v n .

An iterative solver, such as the Newton–Raphson method, can be applied to G v n = g B v n − A v n v n , v n − 2 Δ t B v n − A v n v n − z n − 1 2 , in order to solve for v n such that G ( v n ) = 0 . Once v n is known, the friction force F n can be calculated and the variables u n + 1 , η n + 1 , and z n + 1 2 can be updated as follows: F n = B v n − A v n v n σ 0 + σ 1 v n 2 Δ t B v n − A v n v n − z n − 1 2 , z n + 1 2 = 2 B v n − A v n v n − z n − 1 2 , η n + 1 = 2 Δ t m h 2 Δ t + K h Δ t 2 + Γ h − K h η n + m h 2 Δ t + K h Δ t 2 s n + m h 2 Δ t + K h Δ t 2 + Γ h 2 Δ t η n − 1 − F n , u n + 1 = u n − 1 + 2 Δ t v n + v b n − η n + 1 − η n − 1 2 Δ t .

The solution of the nonlinear equation G ( v n ) = 0 , for G ( v n ) as defined in Equation 32 plays a key role in the algorithm described above. It was shown in Section 3 that in the continuous case, the refined elasto-plastic model has a unique solution. However, this property is not immediately transferable to the numerical case. It remains an open problem to find a threshold on the time step Δ t that would guarantee a unique solution.

 4.3 Numerical energy balance

The stability of the numerical scheme is analyzed by investigating whether the discrete energy balance preserves the passivity of the underlying continuous system.

Multiplying Equation 25 with δ t ⋅ u n , Equation 26 with δ t ⋅ η n , and Equation 29 with σ 0 μ t ○ z n yields m δ t ⋅ u n δ t t u n + κ δ t ⋅ u n u n ⏟ δ t + H r n = − δ t ⋅ u n F n − γ δ t ⋅ u n 2 ⏟ Q r n , m h δ t ⋅ η n δ t t η n + K h δ t ⋅ η n μ t t η n ⏟ δ t + H h n = − δ t ⋅ η n F n − Γ h δ t ⋅ η n 2 ⏟ Q h n , σ 0 μ t ○ z n δ t ○ z n ⏟ δ t + H b n = σ 0 μ t ○ z n g μ t ○ z n , v n .

Summing up the above and using relation Equation 27 yields δ t + H r n + H h n + H b n = − P n − Q r n − Q h n − Q b n , where H r n = m 2 δ t − u n 2 + κ 2 u n u n − 1 ≥ 0 , H h n = m h 2 δ t − η n 2 + K h 2 μ t − η n 2 ≥ 0 , H b n = σ 0 2 z n − 1 2 2 ≥ 0 , P n = v b n F n , Q b n = σ 0 v n 2 μ t ○ z n + g μ t ○ z n , v n σ 1 v n − σ 0 μ t ○ z n ≥ 0 , P is externally supplied power and Q r , Q h , and Q b are dissipation terms with Q r and Q h trivially non-negative. The proof of passivity in the numerical formulation goes line by line as in the continuous case, with z being substituted by μ t ○ z n and v by v n .

The energy balance in Equation 33 induces the following discrete conservation law (Chatziioannou and van Walstijn, 2015): E n = H n + 1 + Δ t ∑ i = 0 n P i + Q r i + Q h i + Q b i = H 0 , where H n = H r n + H h n + H b n . This is the discrete equivalent of Equation 10. The conservation of this quantity, subject to machine precision, can be assessed by monitoring the energy conservation error e n = E n − H 0 .

 4.4 Numerical experiments

To demonstrate the behavior of the model, simulation results are shown in Figure 3. For these simulations, the bow accelerates from 0 at 3.439 m/s2 until it reaches the steady-state value v b (given in Table 1) and then remains constant. The model parameters (Table 1) are set to values found in Matusiak and Chatziioannou (2024) and Pitteroff and Woodhouse (1998a) considering this lumped model hypothesis. The values were obtained for the fundamental mode of a vibrating string (Table 2) according to Equations 64–66 in the Appendix. The energy conservation error e n (where, in this case H 0 = 0 ) is also shown in Figure 3.

Comparison of the original (orange) and refined (blue) elasto-plastic friction model applied to simulating a bowed lumped mass. Model parameters are given in Table 1. Plotted signals are the mass displacement, bow hair displacement, relative velocity, and energy error, as defined in Equation 35. Bottom right plot shows the deviation of the bristle damping coefficient, in the case of the refined model, from the constant value σ ̄ 1 .

Table with parameter values used to generate signals in Figure 3.

Parameter Value Parameter Value
v b Bow velocity [m/s] 0.3439 m Mass [kg] 0.0028
f N Bow force [N] 1.6403 κ Spring stiffness [N/m] 1,055.7
σ 0 Bristle stiffness [N/m] 1 ⋅ 1 0 5 γ Spring damping [kg/s] 0.0095
σ ̄ 1 Bristle damping [kg/s] 0.5
v S Stribeck velocity [m/s] 0.228 m h Bow hair mass [kg] 0.0042
μ C Dynamic friction [-] 0.5071 K h Bow hair stiffness [N/m] 48,297
μ S Static friction [-] 1.0207 Γ h Bow hair damping [kg/s] 57.674

Table with parameter values used to generate the signals in Figure 7. Bow-hair mass, stiffness, and damping are as in the lumped case (see Table 1), divided by the width of the bow hair.

String parameters Value Bow parameters Value
L String length [m] 0.7 a b Bow acceleration [m/s2] 0.8722
r String radius [m] 5 ⋅ 1 0 − 4 v b Bow velocity [m/s] 0.3439
T String tension [N] 149.74 f N Bow force [N/m] 2.3433
ρ Material density [kg/m3] 10,128 W b Bow-hair width [m] 0.01
E Young’s modulus [Pa] 1.37 ⋅ 1 0 10
c Wave speed [m/s] 137.2 σ 0 Bristle stiffness [N/m2] 3.186 ⋅ 1 0 5
γ 0 Freq. independent damping ( s − 1 ) 1.537 σ 1 Bristle damping [kg/(m ⋅ s)] 0.0027
γ 1 Freq. dependent damping [ m 2 /s] 0.0087 v S Stribeck velocity [m/s] 0.228
K T Torsional stiffness [N ⋅ m 2 ] 3.03 ⋅ 1 0 − 4 μ C Dynamic friction [-] 0.5071
P T Polar moment of inertia [kg ⋅ m] 4.2 ⋅ 1 0 − 10 μ S Static friction [-] 1.0207
γ 2 Torsional damping [1/s] 0.0172 β Bow position [-] 0.0786

For this parameter set, the refined model generates signals that are nearly identical to those generated by the original Dupont model. The latter has been used to simulate a string bowed by a finite-width bow and was validated against experimental measurements for the case of a monochord played by a bow (Matusiak and Chatziioannou, 2024). Therefore, it is possible to deduce that the refined model can also reliably resynthesize measured signals.

The difference between the two models comes into play when the Dupont model violates the passivity condition (Equation 34); the two models then behave quite differently (Figure 4). To generate this figure, the bristle stiffness σ 0 was reduced to 500 N/m, the bristle damping σ 1 was increased to 3 kg/s, and the bow force f N was increased to 0.25 N. It can be observed that, in the case of the Dupont model, the total energy loss may become negative, which violates passivity. This results in the friction force not closely following the underlying steady state friction curve. By allowing the bristle damping to vary (bottom left of Figure 4), passivity is guaranteed and the friction force trajectory remains close to the steady-state curve.

Comparison of the behavior of the original (orange) and refined (blue) elasto-plastic friction models when the passivity condition of the original model is violated. The total energy loss of the original model becomes negative, violating passivity, and the friction force deviates significantly from the underlying theoretical steady-state friction curve.

As discussed in Section 4.2, a further issue with this modeling approach is whether the system possesses a unique solution. This can only be shown for the refined model in the continuous case. For the discrete case, the uniqueness of the solution could only be demonstrated empirically. The nonlinear function G ( v ) in Equation 32 is plotted in Figure 5 for both the refined and the Dupont model for increasing sampling rates. Model parameters are as in Table 1, except for σ ̄ 1 = 100 kg/s. The nonlinear function is plotted for time instance t = 0.0299 s. It can be observed that while, for the refined model, G ( v ) has a single root, this is not the case for the original model. Furthermore, the existence of multiple roots in the latter case cannot be avoided by increasing the sampling rate (i.e., oversampling towards the continuous case does not alleviate this issue). While a strict upper bound for Δ t guaranteeing uniqueness is not yet available for the refined model, it has been empirically observed, for a large set of parameter values, that a unique solution exists even for sampling rates lower than the audio sampling rate ( f s = 44100 Hz).

Nonlinear function G ( v ) for the refined model (left) and the original Dupont model (right). OS is an oversampling factor, with sampling rate f s = O S ⋅ f s ̄ and f s ̄ = 44100 Hz.

Furthermore, it is possible to observe that for both models, the derivative of G ( v ) may become equal (or approximately equal) to 0 for certain values of v . While this may hinder the convergence of the Newton–Raphson method, there are alternative approaches that may be used to approximate the root of G ( v ) (e.g., Deuflhard, 2011; Hueso et al., 2009).

Finally, the convergence of the refined model is illustrated in Figure 6. A global error is defined, assuming a reference signal u ̂ that is obtained with 1024 times oversampling, as e g = ∑ u − u ̂ 2 ∑ u ̂ 2

Error of the refined model for different time steps, as defined in Equation 36. Dashed line indicates second-order convergence, while dash-dotted line indicates first order convergence. Model parameters are taken from Table 1.

 5 Distributed system

The insights obtained while studying the bowed-mass system are now applied to a distributed system—a string bowed with a finite width bow. In the following, let u ( x , t ) and u ̃ ( x , t ) be real-valued functions defined over an interval [ 0 , L ] and for time t ≥ 0 , with an inner product and a norm defined as ⟨ u , u ̃ ⟩ = ∫ 0 L u x , t u ̃ x , t   d x , ‖ u ‖ 2 = ⟨ u , u ⟩ . Using the subscripts x and t to denote differentiation with respect to space ( ∂ x ) and time ( ∂ t ) , respectively, the following identities hold. ⟨ ∂ t u , u ⟩ = d d t 1 2 ‖ u ‖ 2 , ⟨ u , ∂ x u ̃ ⟩ = − ⟨ ∂ x u , u ̃ ⟩ + u L , t u ̃ L , t − u 0 , t u ̃ 0 , t , where d d t is the total derivative with respect to time.

5.1 String model

The governing equations for the motion of a string excited by a bow are the equations describing transverse (Equation 39) and torsional (Equation 40) waves. They are coupled through the distributed friction force f ([N/m]), and this force in turn is linked to the bow hair displacement η (Equation 41). The bow hair is modeled as a harmonic oscillator (Pitteroff and Woodhouse, 1998a). The friction force is modelled according to the elasto-plastic friction model. The partial differential equations describing the motion of the bowed string are (Bilbao, 2009; Pitteroff and Woodhouse, 1998b): ρ A ∂ t 2 u = T ∂ x 2 u − E I ∂ x 4 u − 2 γ 0 ρ A ∂ t u + 2 γ 1 ρ A ∂ t ∂ x 2 u − f , P T ∂ t 2 w = K T ∂ x 2 w − 2 γ 2 P T ∂ t w + r f , D h ∂ t 2 η = − K h η − Γ h ∂ t η − f , where ρ is the material density, A = π r 2 is the cross-sectional area of the string with radius r , T is the tension of the string, E is Young’s modulus, I = π r 4 / 4 are the area moment of inertia, and γ 0 and γ 1 represent frequency independent and frequency dependent damping. In addition, P T denotes the polar moment of inertia, K T torsional stiffness, and γ 2 is a torsional damping coefficient. The bow stick is regarded as a rigid frame moving at a given velocity v b and supporting a ribbon of compliant bow-hair of density D h ([kg/m]) with distributed spring and damping constants K h and Γ h , respectively. The relative bow–string velocity is then expressed as v = ∂ t u − r δ t w − v b − ∂ t η .

This model simplifies string damping and omits body coupling. This simplified approach is favored in this case, as including these additional factors would not enhance the presentation of the friction model.

Assuming simply supported ends, the boundary conditions for the transverse movement of the string are u x , t | x = 0 , L = 0 , ∂ x 2 u | x = 0 , L = 0 ,

and for the torsional movement we assume fixed boundary conditions w x , t | x = 0 , L = 0 .

A fourth equation, needed to close the system, describing the friction force f is f z , v = σ 0 z + σ 1 v ∂ t z , where σ 0 and σ 1 are now distributed stiffness and distributed damping, respectively, and ∂ t z is the time derivative of z . It is related to v through ∂ t z = v 1 − α z , v z z ss v , where the adhesion map α is defined in Equation 6 and z ss is a steady-state displacement function given by the Stribeck curve (Equation 9). The damping term σ 1 ( v ) is defined in Equation 14.

 5.2 Energy analysis

The time derivative of the total energy of the combined transverse and torsional movement of the string, bow hair, and bristle energy may be derived by taking an inner product of Equations 39, 40, 41, and 46 with ∂ t u , ∂ t w , ∂ t η , and σ 0 z , respectively, and summing the results. Utilizing Equation 37 and the identity in Equation 38 repeatedly, the energy balance follows: d H d t = − Q − P + B 0 L , where H = H r + H w + H h + H b , Q = Q r + Q w + Q h + Q b , B = B r + B w , and P ( t ) = ⟨ v b , f ⟩ is the power supplied by the bow. The string energy coming from the transverse and torsional motions, H r and H w , respectively, bow hair energy H h , and bristle energy H b are given by H r t = ρ A 2 ‖ ∂ t u ‖ 2 2 + T 2 ‖ ∂ x u ‖ 2 2 + E I 2 ‖ ∂ x 2 u ‖ 2 2 ≥ 0 , H w t = P T 2 ‖ ∂ t w ‖ 2 2 + K T 2 ‖ ∂ x w ‖ 2 2 ≥ 0 , H h t = D h 2 ‖ ∂ t η ‖ 2 2 + K h 2 ‖ η ‖ 2 2 ≥ 0 , H b t = σ 0 2 ‖ z ‖ 2 2 ≥ 0 . The dissipated energies in the string ( Q r and Q w ), bow hair ( Q h ) , and bristles ( Q b ) are given by Q r t = 2 γ 0 ρ A ‖ ∂ t u ‖ 2 2 + 2 γ 1 ρ A ‖ ∂ x ∂ t u ‖ 2 2 ≥ 0 , Q w t = 2 γ 2 P T ‖ ∂ t w ‖ 2 2 ≥ 0 , Q h t = Γ h ‖ ∂ t η ‖ 2 2 ≥ 0 , Q b t = ⟨ σ 1 v v , v ⟩ + σ 0 z − σ 1 v v , α z , v z v z ss v ,

and the boundary terms are B r t = T ∂ t u ∂ x u − E I ∂ t u ∂ x 3 u + E I ∂ t ∂ x u ∂ x 2 u + 2 γ 1 ρ A ∂ t u ∂ t ∂ x u , B w t = K T ∂ t w ∂ x w . Under simply supported boundary conditions, B r vanishes. Similarly, B w vanishes due to fixed boundary conditions for the torsional movement of the string.

Given that H ≥ 0 , the passivity of the system may be assessed by observing the dissipated energy expressions. More precisely, passivity is guaranteed if Q b ≥ 0 . Let W b be the width of the bow, then Q b t = ∫ W b σ 1 v v 2 + α z , v v z z ss v σ 0 z − σ 1 v v d t ≥ 0 since, by Proposition 1, the expression under the integral is positive if σ 1 ( v ) is defined as in Equation 14. Therefore, the refined elasto-plastic friction model results in a passive system, also for this distributed system.

 6 Distributed system—Numerical formulation 6.1 Operators and identities

Let u ( x , t ) be a function defined over an interval [ 0 , L ] and for t ≥ 0 . Let d N = { l Δ x : l = 0 , … , N } be a N + 1 discrete spatial domain corresponding to [ 0 , L ] , with Δ x = L / N . The approximations to u ( x , t ) at points ( l Δ x , n Δ t ) are denoted as u l n . Let u n = [ u 0 n , … , u N n ] T . For two vectors u n and u ̃ n , the discrete inner product and norm on d N are defined as ⟨ u n , u ̃ n ⟩ d N = Δ x ∑ l = 0 N u l n u ̃ l n ,   ‖ u n ‖ d N 2 = ⟨ u n , u n ⟩ d N . Other domains that differ from d N by removing endpoints and will later be used are d N ̄ = l Δ x : l = 1 , … , N , d N ̲ = l Δ x : l = 0 , … , N − 1 , d N ̄ ̲ = l Δ x : l = 1 , … , N − 1 .

The time difference and averaging operators introduced in Sections 4.1, 4.2 are valid in their implementation to grid functions. Similarly, as in the continuous case, the following identities hold: ⟨ δ t ⋅ u n , u n ⟩ d N = δ t + 1 2 ⟨ u n , u n − 1 ⟩ d N , ⟨ δ t ⋅ u n , δ t t u n ⟩ d N = δ t + 1 2 ‖ δ t − u n ‖ d N 2 . Spatial forward, backward, and central discretization operators are defined as δ x + u l n = u l + 1 n − u l n Δ x , δ x − u l n = u l n − u l − 1 n Δ x , δ x ⋅ u l n = u l + 1 n − u l − 1 n 2 Δ x , and can be used to obtain approximations to higher-order partial differential operators: δ x x = δ x − δ x + ,   δ xxxx = δ x x δ x x . Like in the continuous case, the following relation can be derived: ⟨ u n , δ x x u ̃ n ⟩ d N = − ⟨ δ x + u n , δ x + u ̃ n ⟩ d N ̲ + u N n δ x + u ̃ N n − u 0 n δ x + u ̃ 0 n .

6.2 Discretization

Finite-difference schemes for the string in isolation and the bowed string have been described in studies such as Bilbao (2009) and Pitteroff and Woodhouse (1998b). In order to fix the notation, let x B L and x B R be the positions on the string of the inner and outer bow edges, respectively, with the center of the bow lying at x B . Let M be a desired number of grid points under the bow, denoted by x m .

The model is discretized in time and space with functions u l n that are approximations of u ( x , t ) at points ( l Δ x , n Δ t ) . Discretization in time is performed with t = n Δ t , where Δ t = 1 / f s (in s ) with f s the sampling rate (in Hz) and n ∈ N , and in space with x = l Δ x , where the grid spacing Δ x (in m ) for the transverse movement of the string must satisfy the following stability condition (Bilbao, 2009): Δ x ≥ τ + τ 2 + 16 κ 2 Δ t 2 2 , where τ = c 2 Δ t 2 + 4 γ 1 Δ t with c = T / ρ A are the wave speed and κ = E I / ρ A is a stiffness coefficient. The grid points are d N = { l Δ x : l = 0 , … , N } , where N = ⌊ L / Δ x ⌋ ; hence, the total number of grid points is N + 1 . For the torsional movement of the string, the discretization in time is performed as for the transverse motion while in space with x = l Δ x T , where the grid spacing Δ x T (in m ) for the torsional movement of the string must satisfy (Bilbao, 2009): Δ x T ≥ c T Δ t , where c T = K T / P T is the torsional wave speed. The grid points are d N T = { l Δ x T : l = 0 , … , N T } , where N T = ⌊ L / Δ x T ⌋ ; hence the total number of grid points is N T + 1 . Torsional waves travel much faster than transverse waves; therefore, the number of grid points N T is much smaller than N . The spatial discretization operators associated with grid d N T will be denoted with an upper T superscript—for example, δ x + T w l n = w l + 1 n − w l n Δ x T .

For a point x m under the bow, the interpolation vectors i x m and i T x m interpolate the string displacement at position x m for the transverse and torsional motion, respectively. i x m is a row vector of size N + 1 that multiplies the column vector u n = [ u 0 n , … , u N n ] T , and i T x m is a row vector of size N T + 1 that multiplies the column vector w n = [ w 0 n , … , w N T n ] T . The simplest interpolation is the one of 0 th -order where i l x m = 1 for l = ⌊ x m / Δ x ⌋ , and i T , l x m = 1 for l = ⌊ x m / Δ x T ⌋ , respectively, and zeros elsewhere. For the definition of higher orders of interpolation vectors, see Bilbao (2009). On the other hand, a spreading vector j x m is a column vector that distributes the friction force around the bowing point x m on the grid l Δ x . Similarly, j T x m is a spreading vector that distributes the friction force around the bowing point x m on the grid l Δ x T . The spreading and interpolation vectors are related through j x m = 1 Δ x i x m T , j T x m = 1 Δ x T i T x m T .

To simplify the notation, we first divide Equation 39 by ρ A , then discretize it to obtain δ t t u n = c 2 δ x x u n − κ 2 δ xxxx u n − 2 γ 0 δ t ⋅ u n + 2 γ 1 δ t − δ x x u n − M ρ A − 1 J f n , where J = [ j x 1 | … | j x M ] is an ( N + 1 ) × M matrix with m th column being j x m , m = 1 , … , M , and f n = [ f 1 n , … , f M n ] T is a column vector with M rows where each row describes the friction for a point x m of the string that is in contact with the bow. The friction force is discretized using an interleaved grid, with f m n = f μ t ○ z m n , v m n , where f μ t ○ z m n , v m n = σ 0 μ t ○ z m n + σ 1 v m n g μ t ○ z m n , v m n , with g μ t ○ z m n , v m n = v m n 1 − α μ t ○ z m n , v m n μ t ○ z m n z ss v m n , for m = 1 , … , M . For simplicity of notation, let g n denote the M × 1 column vector with entries g m n g m n = g μ t ○ z m n , v m n , then δ t ○ z n = g n .

Assuming simply supported ends, the boundary conditions imply u 0 n = u N n = 0 , u − 1 n = − u 1 n , u N + 1 n = − u N − 1 n , for all n ∈ N .

Similarly, by dividing Equation 40 by P T , it is discretized as follows: δ t t T w n = c T 2 δ x x T w n − 2 γ 2 δ t ⋅ w n + r M P T J T f n , where J T = [ j T x 1 | … | j T x M ] is an ( N T + 1 ) × M matrix with the m th column being j T x m , m = 1 , … , M . Assuming fixed ends, the boundary conditions imply w 0 n = w N T n = 0 , for all n ∈ N .

Discretization of the equation governing bow hair displacement is performed as in the lumped case: D h δ t t η n = − K h μ t t η n − Γ h δ t ⋅ η n − 1 M f n , where η n = [ η 1 n , … , η M n ] T . Here, for each point x m under the bow, the bow hair compliance is computed. Then, if the bow velocity at time n Δ t is v b n , the relative velocities at points of the string in contact with the bow are discretized as v n = I δ t ⋅ u n − r I T δ t ⋅ w n − v b n − δ t ⋅ η n , where I = [ i x 1 ; … ; i x M ] , I T = Δ x Δ x T [ i T x 1 ; … ; i T x M ] T are M × ( N + 1 ) and M × ( N T + 1 ) matrices with the m th row being i x m and i T x m , respectively, for m = 1 , … , M . Let E M be an M by M identity matrix and L = ξ S M ρ A I J + r 2 ξ T M P T I T J T + ξ H M E M , where ξ S = ( 2 Δ t + 2 γ 0 ) − 1 , ξ T = ( 2 Δ t + 2 γ 2 ) − 1 , and ξ H = ( 2 Δ t D h + Δ t 2 K h + Γ h ) − 1 . Utilizing the expression in Equation 30 for the discrete operators δ t t and μ t ○ , Equation 58 becomes v n + v b n = − L ⋅ f n + s n , where s n is an M × 1 column vector with entries s m n = ξ S a m n − ξ T b m n + ξ H c m n , where a m n = i x m c 2 δ x x u n − κ 2 δ xxxx u n + 2 γ 1 δ t − δ x x u n + 2 Δ t δ t − u n , b m n = Δ x Δ x T i T x m r c T 2 δ x x T w n + 2 r Δ t δ t − w n , c m n = Δ t 2 K h + 2 Δ t D h η m n − η m n − 1 Δ t − K h η m n . In order to update the system variables, the vectors v n and μ t ○ z n must first be computed. A system of 2 M equations is formed using Equation 59 and relation in Equation 53 together with δ t ○ z n = 2 Δ t μ t ○ z n − z n − 1 2 . v n + L f n − s n + v b n = 0 , g n − 2 Δ t μ t ○ z n − z n − 1 2 = 0 . The system can then be solved for v n and μ t ○ z n using an iterative solver. Once v n and μ t ○ z n are known, u n + 1 , w n + 1 , and η n + 1 can be updated. First, the variables related to the friction force and the friction force itself are computed, z n + 1 2 = 2 μ t ○ z n − z n − 1 2 , g n = z n + 1 2 − z n − 1 2 Δ t , f n = σ 0 μ t ○ z n + Ξ 1 n g n , where Ξ 1 n is an M × M diagonal matrix with diagonal terms being { σ 1 ( v m n ) } m = 1 M . Then, the string and bow hair variables are updated as u n + 1 = 2 Δ t ξ S A u n + B u n − 1 − 2 Δ t ξ S M ρ A J f n , w n + 1 = 2 Δ t ξ T C w n − 1 − γ 2 Δ t w n − 1 + r 2 Δ t ξ T M P T J T f n , η n + 1 = D 2 D h Δ t 2 − K h 2 η n − D D h Δ t 2 + K h 4 − Γ h 2 Δ t η n − 1 − D M f n , where A = 2 E + 2 γ 1 Δ t + c 2 Δ t 2 δ x x − κ 2 Δ t 2 δ xxxx , B = γ 0 Δ t − 1 E − 2 γ 1 Δ t δ x x , C = 2 E T + c T 2 Δ t 2 δ x x T , D = D h Δ t 2 + K h 4 + Γ h 2 Δ t − 1 , with E an ( N + 1 ) × ( N + 1 ) identity matrix and E T an ( N T + 1 ) × ( N T + 1 ) identity matrix.

Note that Equations 60 and 61 can be reduced to solving just one equation, as was performed in the case of the bowed lumped mass. Using Equation 60 and writing the friction force in terms of μ t ○ z n as f n = σ 0 μ t ○ z n + 2 Δ t Ξ 1 n μ t ○ z n − z n − 1 2 ,

the bristle displacement can be expressed as μ t ○ z n = − Ξ n L − 1 v n + v b n + L − 1 s n + 2 Δ t Ξ 1 n z n − 1 2 , where Ξ n is an M × M diagonal matrix with entries ( σ 0 + 2 Δ t σ 1 ( v m n ) ) − 1 , m = 1 , … , M on the diagonal. Plugging (62) into (61) and g n , only Equation 61 must be solved for v n . z n + 1 / 2 can then be computed from Equation 62. This approach increases computational efficiency for point bowing when M = 1 , but with more points under the bow it involves matrix inversion, which is computationally expensive.

 6.3 Numerical energy and stability condition

An energy balance for the discretized scheme follows from a discrete inner product of Equation 52 with δ t ⋅ u n , an inner product of Equation 55 with δ t ⋅ w n , an inner product of Equation 57 with δ t ⋅ η n , and an inner product of σ 0 M μ t ○ z n with δ t ○ z n . Using summation by parts identities as well as boundary conditions leads to δ t + H n = − Q n − P n + B n , where H , the total numerical energy, is defined as H = H r + H w + H h + H b and, assuming that the stability conditions in Equations 50 and 51 are satisfied, H r n = ρ A 2 δ t − u n d N 2 + T 2 δ x + u n , δ x + u n − 1 d N ̲ + E I 2 δ x x u n , δ x x u n − 1 d N ̲ ̄ ≥ 0 , H w n = P T 2 Δ x Δ x T δ t − w n d N T 2 + K T 2 Δ x Δ x T δ x + T w n , δ x + T w n − 1 d N T ̲ ≥ 0 , H h n = D h 2 δ t − η n T δ t − η n + K h 2 μ t − η n T μ t − η n ≥ 0 , H b n = 1 M σ 0 2 z n − 1 2 T z n − 1 2 ≥ 0 . The total energy lost due to damping is defined as Q = Q r + Q w + Q h + Q b with Q r n = 2 γ 0 ρ A δ t ⋅ u n d N 2 − 2 γ 1 ρ A δ t ⋅ u n , δ t − δ x x u n d N ≥ 0 , Q w n = 2 γ 2 P T Δ x Δ x T δ t ⋅ w n d N T 2 ≥ 0 , Q h n = Γ h δ t ⋅ η n T δ t ⋅ η n ≥ 0 , Q b n = 1 M Ξ 1 n v n T v n + σ 0 μ t ○ z n T o n − Ξ 1 n v n T o n ≥ 0 , where o n is a vector with entries o m n = α ( μ t ○ z m n , v m n ) μ t ○ z m n z ss ( v m n ) for m = 1 , … , M . The non-negativity of the dissipation term Q b n can be shown similarly to the lumped case. The energy supplied by the bow P and boundary terms B = B r + B w is P n = 1 M v b n ∑ m = 1 M f m n , B r n = T δ t ⋅ u N n δ x + u N n − δ t ⋅ u 0 n δ x + u 0 n + E I δ t ⋅ u N n δ x + δ x x u N n − δ t ⋅ u 0 n δ x + δ x x u 0 n − E I δ t ⋅ δ x + u N n δ x x u N n − δ t ⋅ δ x + u 0 n δ x x u 0 n , B w n = K T δ t ⋅ w N n δ x + T w N n − δ t ⋅ w 0 n δ x + T w 0 n .

For simply supported and fixed boundary conditions, the boundary term B vanishes.

The bowed string model is passive, and the argument follows that in the continuous case, where the integral becomes the sum: Q b n = 1 M ∑ m = 1 M σ 1 v m n v m n 2 + α μ t ○ z m n , v m n v m n μ t ○ z m n z ss v m n σ 0 μ t ○ z m n − σ 1 v m n v m n ≥ 0 . Non-negativity of Q b n is guaranteed, since all terms inside the sum are positive, as was the case for the lumped model.

 6.4 Numerical experiments

Simulated signals using the refined bow–string interaction model are shown in Figure 7. In this case, the bow starts in contact with the string, and the force is kept constant while the bow is accelerated from rest with a chosen acceleration value until the steady-state bow velocity v b is reached. The physical model parameters used to generate these signals are given in Table 2.

Simulated bridge force using the refined elasto-plastic model, along with various energy components. Total energy is the sum of the transverse and torsional energies of the string as well as the bristle and bow hair energies. The energy error is shown on the bottom right. Model parameter values are given in Table 2.

The original Dupont friction model was recently applied to the distributed case of bowing a string in Matusiak and Chatziioannou (2024). The performance of the original elasto-plastic model was evaluated by simulating a Guettler diagram and comparing it to a Guettler diagram obtained from measurements (Lampis et al., 2024). A Guettler diagram is generated by choosing a fixed location on the string and bowing from rest with an accelerating bow. The bow force while bowing is kept constant. Bowing is then repeated for different accelerations and bow forces that are incremented in small steps. The number of periods required to reach Helmholtz motion is visualized via the color of each pixel (Figure 8). White corresponds to 0 periods (perfect transient) and black to 20 periods, with intermediate transient lengths resulting in different grayscale values. Guettler (2002) proposed such a diagram as a measure of assessing the playability of a bowed string by measuring the length of the transient necessary to arrive at Helmholz motion.

Guettler diagrams using Dupont (left) and refined models (right) with string parameters as in Table 2. The color bar refers to the number of periods required for Helmholtz motion to be achieved (transient length).

The refined model was utilized to simulate a Guettler diagram with the same parameters as in Matusiak and Chatziioannou (2024). The two diagrams are shown in Figure 8 for comparison. The chaotic nature of the frictional interaction manifests itself via the patchiness of the playability regions. Neighboring pixels may correspond to largely different transient durations, indicating sensitivity to small changes in bow force and acceleration. Therefore, the diagrams slightly differ due to the small underlying numerical differences of the two models. Qualitative observations regarding the playability of the system are nevertheless the same for both models.

 6.5 Supplementary material

A sound example of the synthesis of a fast sautillé passage is included in the supplementary material. Simulation of this bow stroke style exposes the transient behavior of the proposed bow–string model and offers the reader a preliminary aural impression of it. For synthesis of the G2 notes, the bow velocity and force were varied periodically over time according to patterns similar to those observed in Demoucron (2008). Convolution with a measured cello impulse response was applied to the bridge force signal to render a more realistic audio signal.

Furthermore, animations of the string motion are provided for the case shown in Figure 7, where Helmholtz motion is achieved, as well as for a case where the bow force is reduced by 50% ( f N = 1.17 N), resulting in a double-slip pattern. The matlab code for the bowed-string simulations is available at doi:10.5281/zenodo.15341818.

 7 Conclusion

An elasto-plastic friction model has been investigated from the point of view of energy conservation in the continuous domain and discretized using a finite difference scheme. The model was first analyzed in the setting of a bowed lumped mass. A refinement of the model has been suggested that guarantees passivity for elasto-plastic friction with the Stribeck effect and simultaneously leads to the existence and uniqueness of the solution. A numerical scheme has been derived that respects the energy balance of the underlying continuous model, thus leading to a guaranteed passive model and hence to stable simulations.

Based on this refined version of the elasto-plastic model, simulations of a bowed string were revisited, including bow compliance, string torsion, and a finite bow width. While results are similar to those previously obtained using the original Dupont model, the refined model presented is proven to be guaranteed passive, as is the case for the lumped system.

The main limitation of the proposed implicit scheme is that the Jacobian of the nonlinear function to be solved iteratively at each time step can become singular, which—even when applying deliberately modified versions of the iterative solver (e.g., Hueso et al. 2009)—can lead to the necessity of a huge number of iterations. In practice, this often necessitates heavy oversampling, especially when driving the model with articulation parameters (e.g., bowing force) that vary across over time. A logical future research direction is therefore to develop numerical schemes that sidestep the need for an iterative solver, as has been achieved recently for numerical simulation of various other nonlinear phenomena in musical instruments, including collisions (e.g., van Walstijn et al., 2024).

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

 Author contributions

EM: conceptualization, investigation, methodology, software, and writing – original draft. VC: funding acquisition, investigation, project administration, validation, and writing – review and editing. MV: investigation, validation, and writing – review and editing.

 Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded in whole or in part by the Austrian Science Fund (FWF) [10.55776/P34852]. For open access purposes, the authors have applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsip.2025.1525044/full#supplementary-material

References Barreira L. Valls C. (2012). Ordinary differential equations: qualitative theory. America, American Mathematical Society. Bilbao S. (2009). Numerical sound synthesis. J. Wiley and Sons. 10.1002/9780470749012 Bilbao S. Torin A. Chatziioannou V. (2015). Numerical modeling of collisions in musical instruments. Acta Acust. United Acust. 101, 155–173. 10.3813/AAA.918813 Chabassier J. Joly P. (2010). Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: application to the vibrating piano string. Comput. Methods Appl. Mech. Eng. 199, 2779–2795. 10.1016/j.cma.2010.04.013 Chatziioannou V. van Walstijn M. (2015). “Discrete-time conserved quantities for damped oscillators,” in Proceedings of the third Vienna talk on music acoustics (Vienna, AT: Department of Music Acoustics - Wiener Klangstil), 135–139. Demoucron M. (2008). On the control of virtual violins physical modelling and control of bowed string instruments (Stockholm: Royal Institute of Technology). Ph.D. thesis. Desvages C. (2018). Physical modelling of the bowed string and applications to sound synthesis (The University of Edinburgh). Ph.D. thesis. Deuflhard P. (2011). Newton methods for nonlinear problems: affine invariance and adaptive algorithms. Springer Sci. and Bus. Media 35. de Wit C. Olsson H. Åström K. Lischinsky P. (2024). A new model for control of systems with friction. IEEE Trans. Autom. Control 40, 419–425. 10.1109/9.376053 Ducceschi M. Bilbao S. (2022). Simulation of the geometrically exact nonlinear string via energy quadratisation. J. Sound Vib. 534, 117021. 10.1016/j.jsv.2022.117021 Dupont P. Armstrong B. Hayward V. (2000). Elasto-plastic friction model: contact compliance and stiction. Proc. 2000 Acc. IEEE Chic., 1072–1077 vol.2. 10.1109/ACC.2000.876665 Dupont P. Hayward V. Armstrong B. Altpeter F. (2002). Single state elasto-plastic friction models. IEEE Trans. Autom. Control 47, 787–792. 10.1109/TAC.2002.1000274 Falaize A. Roze D. (2024). Generic passive-guaranteed nonlinear interaction model and structure-preserving spatial discretization procedure with applications in musical acoustics. Nonlinear Dyn. 112, 3249–3275. 10.1007/s11071-024-10438-9 Galluzzo P. (2004). On the playability of stringed instruments. Ph.D. thesis. 10.17863/CAM.14046 Guettler K. (2002). On the creation of the helmholtz motion in bowed strings. Acta Acust. united Acust. 88, 970–985. Hueso J. Martínez E. Torregrosa J. (2009). Modified Newton’s method for systems of nonlinear equations with singular Jacobian. J. Comput. Appl. Math. 224, 77–83. 10.1016/j.cam.2008.04.013 Lampis A. Mayer A. Chatziioannou V. (2024). Assessing playability limits of bowed-string transients using experimental measurements. Acta Acust. 8, 44. 10.1051/aacus/2024034 Maestre E. Spa C. Smith J. (2014). A bowed string physical model including finite-width thermal friction and hair dynamics. ICMC. Matusiak E. Chatziioannou V. (2024). Elasto-plastic friction modeling toward reconstructing measured bowed-string transients. J. Acoust. Soc. Am. 156, 1135–1147. 10.1121/10.0028228 Olsson H. (1996). Control systems with friction (Lund Institute of Technology LTH). Ph.D. thesis. Pitteroff R. Woodhouse J. (1998a). Mechanics of the contact area between a violin bow and a string. Part I: reflection and transmission behaviour. Acta Acust. United Acust. 84, 543–562. Pitteroff R. Woodhouse J. (1998b). Mechanics of the contact area between a violin bow and a string. Part II: simulating the bowed string. Acta Acust. United Acust. 84, 744–757. Serafin S. (2004). The sound of friction: real-time models, playability and musical applications (Stanford: Stanford University). Ph.D. thesis. Serafin S. Avanzini F. Rocchesso D. (2003). “Bowed string simulations using an elasto-plastic friction model,” in Proceedings of the Stockholm Music Accoustics Conference, USA, 14–15 June 2023. Smith J. Woodhouse J. (1999). The tribology of rosin. J. Mech. Phys. Solids 48, 1633–1681. 10.1016/S0022-5096(99)00067-8 van Walstijn M. Chatziioannou V. Bhanuprakash A. (2024). Implicit and explicit schemes for energy-stable simulation of string vibrations with collisions: refinement, analysis, and comparison. J. Sound Vib. 569, 117968. 10.1016/j.jsv.2023.117968 Willemsen S. (2021). “Real-time simulation of musical instruments using finite-difference time-domain methods,”.Aalb. Univ. Ph.D. thesis. 10.54337/aau451025772 Willemsen S. Bilbao S. Serafin S. (2019). Real-time implementation of an elasto-plastic friction model applied to stiff strings using finite-difference schemes. Proceedings of the International Conference on Digital Audio Effects, USA, 8-10 Sept. 2021, 40–46. Woodhouse J. (2003). Bowed string simulation using a thermal friction model. Acta Acustica united Acustica 89, 355–368. 10.25144/18216 Appendix

The modal expansion for the displacement of the string, assuming simply supported boundary conditions, is u x , t = ∑ i = 1 N M q i t ψ i x , where ψ i x = sin i π L x , where N M is the number of modes, and it is normally set according to the relevant frequency range. To isolate a single mode of vibration, u ( x , t ) is substituted into the equation governing the transverse motion of the string, ρ A ∂ t 2 u = T ∂ x 2 u − E I ∂ x 4 u − 2 γ 0 ρ A ∂ t u + 2 γ 1 ρ A ∂ t ∂ x 2 u , and an inner product with ψ n is taken. Since ⟨ ψ i , ψ n ⟩ = L 2 δ i , n , we obtain m q ̈ n = − κ n q n − γ ~ n q ̇ n , where m = ρ A L 2 , κ n = n π 2 2 L T + E I n π L 2 , γ ~ n = ρ A L γ 0 + γ 1 n π L 2 . Therefore, by setting N M = 1 , the dynamics of the system reduce to a damped harmonic oscillator.