Numerical modelling of elasto-plastic friction in bow–string interaction with guaranteed passivity

Abstract

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.

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 undergoing a tangential friction force (Figure 1). The mass is excited by a bow moving with velocity which is modeled as a harmonic oscillator with bow hair mass , stiffness [kg/s²], and damping [kg/s]. The bow hair displacement relative to the rigid bow is denoted by . The motion of the two masses can be then described bywhere denotes the displacement of the mass attached to a spring of stiffness [kg/s²] and damping [kg/s]. The relative velocity between bow and mass is given by Equation 3.

FIGURE 1

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 , the average bristle deflection, and by , the relative velocity between the string and the bow, the model is described as follows.where in [kg/s²] and 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 aswhereand where is the breakaway displacementand is the steady-state displacement for constant velocitieswith and . When , is referred to as “Stribeck velocity.” In the numerical experiments in Sections 4.4 and 6.4, is used. Here, and denote static and dynamic friction coefficients, respectively, and is the normal force applied by the bow. For small bristle displacements, when , , and consequently , a purely elastic and reversible regime is entered, referred to as “pre-sliding” (sticking). For larger displacements—that is, when —some bristles start to break, and a mixed elasto-plastic sliding occurs. Finally, for , all bristles break, and a purely plastic regime is achieved—the string slips under the bow. In that situation, . 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 , Equation 2 with , and Equation 5 with , summing up and using the relation in Equation 3 yields the energy balancewhere is the power supplied by the bow via the friction force, and , , and are dissipation terms, corresponding to the oscillator, the bow hair, and the bristles, respectively. This leads to the following conservation law:where is the initial system energy (if any). and are trivially non-negative, and is non-negative whenever . For , and for high—relative to —values of , 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 to guarantee passivity for the elasto-plastic friction model. The condition reads:

However, without knowing the maximal relative velocity, , it is not possible to set the value for . 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 .

Falaize and Roze (2024) employed the elasto-plastic friction model in the framework of port-Hamiltonian systems and arrived at the same dissipation term— . 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 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 can be thought of as a parametric curved defined by and whose rate of change isandThis means that for , if reaches the maximal value of which is , it cannot rise any further but has to either decrease or stay constant. Similarly, for , if reaches the minimal value of , which equals , it cannot decrease any further but has to either rise or stay constant. Therefore, (see Figure 2 for visualization). The boundedness of by itself does not, however, imply stability.

FIGURE 2

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 bowed with velocity as described in Section 2, by Equations 1–5, where bristle damping is not constant but varies with and as follows:such that and .

A similar velocity-dependent damping term was introduced by Olsson (1996) for the LuGre model, motivated by the need to reproduce certain friction phenomena. Provided that the additional free parameter 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 values, stays closer to the original constant than with the exponential formula of Olsson.

3.1 Energy balance

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

Proposition 1: Let be a friction force as defined in Equation 4 with bristle damping 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 . For , is always non-negative. Therefore, only the case when , which happens when and , is considered here. In that case, we have

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.

Let , then the above can be written aswhere with the functions corresponding to Equations 15–19. Given initial conditions , existence and uniqueness of the solution to Equation 20 is guaranteed whenever is Lipschitz continuous (Barreira and Valls, 2012).

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 at points are denoted as , where is the time step. The variable is approximated at an interleaved grid—that is, denotes the approximation to at time . The following centered, second-order accurate discretization operators are defined:

In addition, one may define non-centered, first-order accurate discretization operators

The composition of first-order accurate operators results in second-order accurate operators,

and several useful identities can be constructed, including

4.2 Discretization

For simplicity of notation, letThe system of Equations 1–5 is discretized as follows:

In the case of the original Dupont model, . 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):

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 . Using the discretization operators, one can establish identities

Substituting into Equations 25, 26, and 28, one obtainsConsidering , , , , and to be known, let and . Using Equation 27, can be expressed as a function of whereFor the original Dupont model, and are independent of , and is linearly related to . Substituting Equation 31 into Equation 29 and using the middle identity in Equation 30 yields a nonlinear equation in the unknown ,

An iterative solver, such as the Newton–Raphson method, can be applied toin order to solve for such that . Once is known, the friction force can be calculated and the variables , , and can be updated as follows:

The solution of the nonlinear equation , for 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 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 , Equation 26 with , and Equation 29 with yields

Summing up the above and using relation Equation 27 yieldswhere is externally supplied power and , , and are dissipation terms with and trivially non-negative. The proof of passivity in the numerical formulation goes line by line as in the continuous case, with being substituted by and by .

The energy balance in Equation 33 induces the following discrete conservation law (Chatziioannou and van Walstijn, 2015):where . 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 .

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/s² until it reaches the steady-state value (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 (where, in this case ) is also shown in Figure 3.

FIGURE 3

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 was reduced to 500 N/m, the bristle damping was increased to 3 kg/s, and the bow force 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.

FIGURE 4

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 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 kg/s. The nonlinear function is plotted for time instance s. It can be observed that while, for the refined model, 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 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 ( Hz).

FIGURE 5

Furthermore, it is possible to observe that for both models, the derivative of may become equal (or approximately equal) to 0 for certain values of . While this may hinder the convergence of the Newton–Raphson method, there are alternative approaches that may be used to approximate the root of (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 that is obtained with 1024 times oversampling, as

FIGURE 6

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 and be real-valued functions defined over an interval and for time , with an inner product and a norm defined asUsing the subscripts and to denote differentiation with respect to space and time , respectively, the following identities hold.where 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 ([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):where is the material density, is the cross-sectional area of the string with radius , is the tension of the string, is Young’s modulus, are the area moment of inertia, and and represent frequency independent and frequency dependent damping. In addition, denotes the polar moment of inertia, torsional stiffness, and is a torsional damping coefficient. The bow stick is regarded as a rigid frame moving at a given velocity and supporting a ribbon of compliant bow-hair of density ([kg/m]) with distributed spring and damping constants and , respectively. The relative bow–string velocity is then expressed as

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

and for the torsional movement we assume fixed boundary conditions

A fourth equation, needed to close the system, describing the friction force iswhere and are now distributed stiffness and distributed damping, respectively, and is the time derivative of . It is related to throughwhere the adhesion map is defined in Equation 6 and is a steady-state displacement function given by the Stribeck curve (Equation 9). The damping term 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 , , , and , respectively, and summing the results. Utilizing Equation 37 and the identity in Equation 38 repeatedly, the energy balance follows:whereand is the power supplied by the bow. The string energy coming from the transverse and torsional motions, and , respectively, bow hair energy , and bristle energy are given byThe dissipated energies in the string ( and ), bow hair , and bristles are given by

and the boundary terms areUnder simply supported boundary conditions, vanishes. Similarly, vanishes due to fixed boundary conditions for the torsional movement of the string.

Given that , the passivity of the system may be assessed by observing the dissipated energy expressions. More precisely, passivity is guaranteed if . Let be the width of the bow, thensince, by Proposition 1, the expression under the integral is positive if 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 be a function defined over an interval and for . Let be a discrete spatial domain corresponding to , with . The approximations to at points are denoted as . Let . For two vectors and , the discrete inner product and norm on are defined asOther domains that differ from by removing endpoints and will later be used are

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:Spatial forward, backward, and central discretization operators are defined asand can be used to obtain approximations to higher-order partial differential operators:Like in the continuous case, the following relation can be derived:

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 and be the positions on the string of the inner and outer bow edges, respectively, with the center of the bow lying at . Let be a desired number of grid points under the bow, denoted by .

The model is discretized in time and space with functions that are approximations of at points . Discretization in time is performed with , where (in ) with the sampling rate (in Hz) and , and in space with , where the grid spacing (in ) for the transverse movement of the string must satisfy the following stability condition (Bilbao, 2009):where with are the wave speed and is a stiffness coefficient. The grid points are , where ; hence, the total number of grid points is . For the torsional movement of the string, the discretization in time is performed as for the transverse motion while in space with , where the grid spacing (in ) for the torsional movement of the string must satisfy (Bilbao, 2009):where is the torsional wave speed. The grid points are , where ; hence the total number of grid points is . Torsional waves travel much faster than transverse waves; therefore, the number of grid points is much smaller than . The spatial discretization operators associated with grid will be denoted with an upper superscript—for example, .

For a point under the bow, the interpolation vectors and interpolate the string displacement at position for the transverse and torsional motion, respectively. is a row vector of size that multiplies the column vector , and is a row vector of size that multiplies the column vector . The simplest interpolation is the one of -order where for , and for , respectively, and zeros elsewhere. For the definition of higher orders of interpolation vectors, see Bilbao (2009). On the other hand, a spreading vector is a column vector that distributes the friction force around the bowing point on the grid . Similarly, is a spreading vector that distributes the friction force around the bowing point on the grid . The spreading and interpolation vectors are related through

To simplify the notation, we first divide Equation 39 by , then discretize it to obtainwhere is an matrix with ^th column being , , and is a column vector with rows where each row describes the friction for a point of the string that is in contact with the bow. The friction force is discretized using an interleaved grid, withwithfor . For simplicity of notation, let denote the column vector with entries then

Assuming simply supported ends, the boundary conditions implyfor all .

Similarly, by dividing Equation 40 by , it is discretized as follows:where is an matrix with the ^th column being , . Assuming fixed ends, the boundary conditions implyfor all .

Discretization of the equation governing bow hair displacement is performed as in the lumped case:where . Here, for each point under the bow, the bow hair compliance is computed. Then, if the bow velocity at time is , the relative velocities at points of the string in contact with the bow are discretized aswhere , are and matrices with the ^th row being and , respectively, for . Let be an by identity matrix andwhere , , and . Utilizing the expression in Equation 30 for the discrete operators and , Equation 58 becomeswhere is an column vector with entries , whereIn order to update the system variables, the vectors and must first be computed. A system of equations is formed using Equation 59 and relation in Equation 53 together with .The system can then be solved for and using an iterative solver. Once and are known, , , and can be updated. First, the variables related to the friction force and the friction force itself are computed,where is an diagonal matrix with diagonal terms being . Then, the string and bow hair variables are updated aswherewith an identity matrix and an 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 as

the bristle displacement can be expressed aswhere is an diagonal matrix with entries , on the diagonal. Plugging (62) into (61) and , only Equation 61 must be solved for . can then be computed from Equation 62. This approach increases computational efficiency for point bowing when , 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 , an inner product of Equation 55 with , an inner product of Equation 57 with , and an inner product of with . Using summation by parts identities as well as boundary conditions leads towhere , the total numerical energy, is defined as and, assuming that the stability conditions in Equations 50 and 51 are satisfied,The total energy lost due to damping is defined as withwhere is a vector with entries for . The non-negativity of the dissipation term can be shown similarly to the lumped case. The energy supplied by the bow and boundary terms is

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

The bowed string model is passive, and the argument follows that in the continuous case, where the integral becomes the sum:Non-negativity of 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 is reached. The physical model parameters used to generate these signals are given in Table 2.

FIGURE 7

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.

FIGURE 8

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% ( 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).

References

• 1

  BarreiraL.VallsC. (2012). Ordinary differential equations: qualitative theory. America, American Mathematical Society.

  • Google Scholar
  • View reference in article

• 2

  BilbaoS. (2009). Numerical sound synthesis. J. Wiley and Sons. 10.1002/9780470749012

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BilbaoS.TorinA.ChatziioannouV. (2015). Numerical modeling of collisions in musical instruments. Acta Acust. United Acust.101, 155–173. 10.3813/AAA.918813

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  ChabassierJ.JolyP. (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

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  ChatziioannouV.van WalstijnM. (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.

  • Google Scholar
  • View reference in article

• 6

  DemoucronM. (2008). On the control of virtual violins physical modelling and control of bowed string instruments (Stockholm: Royal Institute of Technology). Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 7

  DesvagesC. (2018). Physical modelling of the bowed string and applications to sound synthesis (The University of Edinburgh). Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 8

  DeuflhardP. (2011). Newton methods for nonlinear problems: affine invariance and adaptive algorithms. Springer Sci. and Bus. Media35.

  • Google Scholar
  • View reference in article

• 9

  de WitC.OlssonH.ÅströmK.LischinskyP. (2024). A new model for control of systems with friction. IEEE Trans. Autom. Control40, 419–425. 10.1109/9.376053

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  DucceschiM.BilbaoS. (2022). Simulation of the geometrically exact nonlinear string via energy quadratisation. J. Sound Vib.534, 117021. 10.1016/j.jsv.2022.117021

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  DupontP.ArmstrongB.HaywardV. (2000). Elasto-plastic friction model: contact compliance and stiction. Proc. 2000 Acc. IEEE Chic., 1072–1077 vol.2. 10.1109/ACC.2000.876665

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  DupontP.HaywardV.ArmstrongB.AltpeterF. (2002). Single state elasto-plastic friction models. IEEE Trans. Autom. Control47, 787–792. 10.1109/TAC.2002.1000274

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FalaizeA.RozeD. (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

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  GalluzzoP. (2004). On the playability of stringed instruments. Ph.D. thesis. 10.17863/CAM.14046

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  GuettlerK. (2002). On the creation of the helmholtz motion in bowed strings. Acta Acust. united Acust.88, 970–985.

  • Google Scholar
  • View reference in article

• 16

  HuesoJ.MartínezE.TorregrosaJ. (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

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  LampisA.MayerA.ChatziioannouV. (2024). Assessing playability limits of bowed-string transients using experimental measurements. Acta Acust.8, 44. 10.1051/aacus/2024034

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  MaestreE.SpaC.SmithJ. (2014). A bowed string physical model including finite-width thermal friction and hair dynamics. ICMC.

  • Google Scholar
  • View reference in article

• 19

  MatusiakE.ChatziioannouV. (2024). Elasto-plastic friction modeling toward reconstructing measured bowed-string transients. J. Acoust. Soc. Am.156, 1135–1147. 10.1121/10.0028228

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  OlssonH. (1996). Control systems with friction (Lund Institute of Technology LTH). Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 21

  PitteroffR.WoodhouseJ. (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.

  • Google Scholar
  • View reference in article

• 22

  PitteroffR.WoodhouseJ. (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.

  • Google Scholar
  • View reference in article

• 23

  SerafinS. (2004). The sound of friction: real-time models, playability and musical applications (Stanford: Stanford University). Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 24

  SerafinS.AvanziniF.RocchessoD. (2003). “Bowed string simulations using an elasto-plastic friction model,” in Proceedings of the Stockholm Music Accoustics Conference, USA, 14–15 June 2023.

  • Google Scholar
  • View reference in article

• 25

  SmithJ.WoodhouseJ. (1999). The tribology of rosin. J. Mech. Phys. Solids48, 1633–1681. 10.1016/S0022-5096(99)00067-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  van WalstijnM.ChatziioannouV.BhanuprakashA. (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

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  WillemsenS. (2021). “Real-time simulation of musical instruments using finite-difference time-domain methods,”.Aalb. Univ.Ph.D. thesis. 10.54337/aau451025772

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  WillemsenS.BilbaoS.SerafinS. (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.

  • Google Scholar
  • View reference in article

• 29

  WoodhouseJ. (2003). Bowed string simulation using a thermal friction model. Acta Acustica united Acustica89, 355–368. 10.25144/18216

  • CrossRef
  • Google Scholar
  • View reference in article