ORIGINAL RESEARCH article

Front. Signal Process., 21 May 2025

Sec. Audio and Acoustic Signal Processing

Volume 5 - 2025 | https://doi.org/10.3389/frsip.2025.1525044

This article is part of the Research TopicSound Synthesis through Physical ModelingView all 5 articles

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

  • 1Department of Music Acoustics (IWK), University of Music and Performing Arts Vienna, Vienna, Austria
  • 2Centre for Interdisciplinary Research in Music and Sound, School of Electronics, Electrical Engineering, and Computer Science, Queen’s University Belfast, Belfast, United Kingdom

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 m undergoing a tangential friction force F (Figure 1). The mass is excited by a bow moving with velocity vb which is modeled as a harmonic oscillator with bow hair mass mh, stiffness Kh [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̇F,(1)
mhη̈=KhηΓhη̇F,(2)
v=u̇vbη̇,(3)

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.

Figure 1
www.frontiersin.org

Figure 1. 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=σ0z+σ1ż,(4)
ż=v1αz,vzzssv,(5)

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,vz0αmz,v,vz>0,(6)

where

αmz,v=0,|z|zbaᾱmv,z,zba<|z|<|zssv|1,|z||zssv|,(7)

and where zbaμCfN/σ0 is the breakaway displacement

ᾱmz,v=12sinπθz,v+1,θz,v=|z|12|zssv|+zba|zssv|zba,(8)

and zss is the steady-state displacement for constant velocities

zssv=fNσ0μC+μSμCe|v/vS|p,v0fNσ0μC+μSμCe|v/vS|p,v<0,(9)

with vS>0 and p1. When p=2, vS 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 fN is the normal force applied by the bow. For small bristle displacements, when |z|zba, α(v,z)=0, and consequently ż=v, a purely elastic and reversible regime is entered, referred to as “pre-sliding” (sticking). For larger displacements—that is, when zba<|z|<|zss(v)|—some bristles start to break, and a mixed elasto-plastic sliding occurs. Finally, for |z||zss(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 15, the energy balance of the system is considered.

Multiplying Equation 1 with u̇, Equation 2 with η̇, and Equation 5 with σ0z, summing up and using the relation in Equation 3 yields the energy balance

müu̇+κuu̇Ḣr+mhη̈η̇+Khηη̇Ḣh+σ0zżḢb=vbFPγu̇2QrΓhη̇2Qhσ1v2+αz,vvzzssvσ0zσ1vQb,

where

Hr=m2u̇2+κ2u20istheoscillatorkinetic&potentialenergy,Hh=mh2η̇2+Kh2η20isthebowhairkinetic&potentialenergyandHb=σ02z20istheenergystoredinthebristles.

P is the power supplied by the bow via the friction force, and Qr, Qh, and Qb are dissipation terms, corresponding to the oscillator, the bow hair, and the bristles, respectively. This leads to the following conservation law:

Hr+Hh+Hb+P+Qr+Qh+Qbdt=H0,(10)

where H(0) is the initial system energy (if any). Qr and Qh are trivially non-negative, and Qb is non-negative whenever α1. For α=1, and for high—relative to σ0—values of σ1, Qb 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:

σ14σ0zssv|v|.(11)

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— Qb. 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

żt=vt1αzt,vtztzssvt,(12)

and

signżt=1,zt>zssvt0,zt=zssvtorvt=0.1,zt<zssvt(13)

This means that for v(t)>0, if z(t) reaches the maximal value of zss(v(t)) which is zss(0)=μSfNσ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 zss(v(t)), which equals limv0zss(v(t))=μSfNσ0, it cannot decrease any further but has to either rise or stay constant. Therefore, |z(t)|μSfNσ0 (see Figure 2 for visualization). The boundedness of z by itself does not, however, imply stability.

Figure 2
www.frontiersin.org

Figure 2. Behavior of the bristle displacement z(t). Left: pink color indicates regions where z(t) is decreasing (sign(ż(t))=1) and yellow is where z(t) is increasing (sign(ż(t))=1). Blue line is the region where ż(t)=0, and the black dashed line denotes the maximum and minimum value of zss. Right: Example trajectory in the vz 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 vb as described in Section 2, by Equations 15, where bristle damping σ1 is not constant but varies with v and fN as follows:

σ1v=μCfNv2+ϵ2,ϵ=μCfNσ̄1,(14)

such that 0σ1(v)σ̄1 and σ1(0)=σ̄1.

A similar velocity-dependent damping term σ1(v)=σ̄1e(v/vd2) was introduced by Olsson (1996) for the LuGre model, motivated by the need to reproduce certain friction phenomena. Provided that the additional free parameter vd 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 15 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 15 is passive.

Proof: All terms in the energy balance (Section 2.1) are trivially non-negative, apart from Qb. For α(v,z)<1, Qb is always non-negative. Therefore, only the case when α(z,v)=1, which happens when sign(z)=sign(v) and |z||zss(v)|, is considered here. In that case, we have

Qb=σ1vv2+|z||v||zssv|σ0|z|σ1v|v|σ1vv2+|z||v||zssv|σ0|zssv|σ1v|v|σ1vv2+|z||v||zssv|μCfNσ1v|v|=σ1vv2+|z||v||zssv|μCfN1|v|v2+ϵ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 15 together with Equation 14 can be written as an autonomous system of equations.

u̇=ϕ,(15)
η̇=v+vbϕ,(16)
ϕ̇=1mκuγϕσ0zσ1vv1αz,vzzssv,(17)
v̇=κuγϕm+KhηΓhv+vbϕmhm+mhmmhσ0z+σ1vv1αz,vzzssv,(18)
ż=v1αz,vzzssv.(19)

Let x=(u,η,ϕ,v,z), then the above can be written as

ẋ=Sx,(20)

where S=(su,sη,sϕ,sv,sz) with the functions s corresponding to Equations 1519. Given initial conditions x(t0)=(u0,η0,ϕ0,v0,z0), existence and uniqueness of the solution to Equation 20 is guaranteed whenever S(x) is Lipschitz continuous (Barreira and Valls, 2012).

Definition 1:. A vector-valued function S(x) is Lipschitz continuous if there exists an L0, called the “Lipschitz constant,” such that for all x,y in the domain of S

|SxSy|L|xy|.(21)

If S(x)=(s1(x),,sK(x)), then the norm of S(x) is defined as

|Sx|2=i=1K|six|2.

The norm of x is similarly defined.

Proposition 2: Consider a bowed mass m described by Equations 15 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 su,,sz being Lipschitz continuous with respect to each variable u, η, ϕ, v, and z.

By passivity, all the variables are bounded. Lipschitz continuity of su and sη is trivial. Similarly, Lipschitz continuity with respect to u, η, and ϕ is trivially satisfied for all the functions su,sη,sϕ,sv,sz. 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,vzvzssv(22)

is Lipschitz continuous.

(i) Let v be fixed. Then, for z1,z2,

g̃z1,vg̃z2,v|v||zssv||αz1,v||z1z2|+|z2||αz1,vαz2,v||v||zssv||z1z2|+π2μSfNσ0θz1,vθz2,v,

where the last inequality follows from the Lipschitz continuity of the sin function with Lipschitz constant equal to 1, with |sin(π2θ(z1,v))sin(π2θ(z2,v))|π2|θ(z1,v)θ(z2,v)|. The function θ is also Lipschitz continuous with respect to z,

θz1,vθz2,v=||z1||z2|||zssv|zba|z1z2||zssv|zbaσ0μCfNσ0zba|z1z2|.

Let v be bounded by Bv, then

g̃z1,vg̃z2,vBvσ0μCfN+πBvσ0μSμCμCfNσ0zbaLz|z1z2|.

Hence, g̃(z,v) is Lipschitz continuous with respect to z with Lipschitz constant Lz.

(ii) Now, let z be fixed and z>0. Then

g̃z,v=0,v0αz,vzvzssv,v0

is continuous and everywhere differentiable except for v=0, with

g̃v=0,v<0αvzvzss+αzzsszvżsszss2,v>0,

and

αv=π2cosπθżsszba|z|zsszba2.

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 v1,v2>0

g̃z,v1g̃z,v2Lv1|v1v2|,Lv1=maxg̃v,(23)

and g̃ is Lipschitz continuous for v>0, with Lipschitz constant Lv(1), and trivially for v<0, with any Lv(2)0 as a Lipschitz constant. Now let v1>0 and v20, then

g̃z,v1g̃z,v2|v1v2|=g̃z,v1|v1v2|g̃z,v1|v1|=|z||zssv1|μSμC.

Let Lv(3)=μS/μC. Hence, g̃ is Lipschitz continuous with respect to v with Lipschitz constant Lv=max{Lv(1),Lv(2),Lv(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 15 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 15. The proposed discretization is energy conserving.

4.1 Numerical preliminaries

The approximations to u(t) at points nΔt are denoted as un, where Δt is the time step. The variable z is approximated at an interleaved grid—that is, zn+12 denotes the approximation to z(t) at time t=(n+12)Δt. The following centered, second-order accurate discretization operators are defined:

δtun=un+12un12Δt=u̇nΔt+OΔt2,δtun=un+1un12Δt=u̇nΔt+OΔt2,μtun=un+12+un122=unΔt+OΔt2,μtun=un+1+un12=unΔt+OΔt2.

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

δtun=unun1Δt=u̇nΔt+OΔt,δt+un=un+1unΔt=u̇nΔt+OΔt,μtun=un+un12=unΔt+OΔt,μt+un=un+1+un2=unΔt+OΔt.

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

δttun=δt+δtun=un+12un+un1Δt2=ünΔt+OΔt2,μttun=μt+μtun=un+1+2un+un14=unΔt+OΔt2,

and several useful identities can be constructed, including

δtunδttun=δt+12δtun2,δtunμttun=δt+12μtun2,δtunun=δt+12unun1,δtznμtzn=δt+12(zn122.

4.2 Discretization

For simplicity of notation, let

gz,v=v1αz,vzzssv.(24)

The system of Equations 15 is discretized as follows:

mδttun=κunγδtunFn,(25)
mhδttηn=KhμttηnΓhδtηnFn,(26)
vn=δtunvbn+δtηn,(27)
Fn=σ0μtzn+σ1vnδtzn,(28)
δtzn=gμtzn,vn.(29)

In the case of the original Dupont model, σ1(vn)=σ̄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<2m/κ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 rn=μtzn. Using the discretization operators, one can establish identities

δttun=2Δtδtunδtun,δtzn=2Δtrnzn12,μttun=Δt2δtunδtun+un.(30)

Substituting into Equations 25, 26, and 28, one obtains

m2Δt+γδtun=κun+m2Δtδtunσ0rnσ1vn2Δtrnzn12,mh2Δt+KhΔt2+Γhδtηn=Khηn+mh2Δt+KhΔt2δtηnσ0rnσ1vn2Δtrnzn12.

Considering un, un1, ηn, ηn1, and zn12 to be known, let pn=δtun and sn=δtηn. Using Equation 27, rn can be expressed as a function of vn

rn=Avnvn+Bvn,(31)

where

Avn=m2Δt+γmh2Δt+KhΔt2+Γh2Δtm+mh+KhΔt2+γ+Γhσ0+2Δtσ1vn1,Bvn=Avnκun+m2Δtpnm2Δt+γ+Khηn+mh2Δt+KhΔt2snmh2Δt+KhΔt2+Γhvbn+2Δtσ1vnzn12σ0+2Δtσ1vn.

For the original Dupont model, A(vn) and B(vn) are independent of vn, and rn is linearly related to vn. Substituting Equation 31 into Equation 29 and using the middle identity in Equation 30 yields a nonlinear equation in the unknown vn,

2ΔtBvnAvnvnzn12=gBvnAvnvn,vn.

An iterative solver, such as the Newton–Raphson method, can be applied to

Gvn=gBvnAvnvn,vn2ΔtBvnAvnvnzn12,(32)

in order to solve for vn such that G(vn)=0. Once vn is known, the friction force Fn can be calculated and the variables un+1, ηn+1, and zn+12 can be updated as follows:

Fn=BvnAvnvnσ0+σ1vn2ΔtBvnAvnvnzn12,zn+12=2BvnAvnvnzn12,ηn+1=2Δtmh2Δt+KhΔt2+ΓhKhηn+mh2Δt+KhΔt2sn+mh2Δt+KhΔt2+Γh2Δtηn1Fn,un+1=un1+2Δtvn+vbnηn+1ηn12Δt.

The solution of the nonlinear equation G(vn)=0, for G(vn) 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 δtun, Equation 26 with δtηn, and Equation 29 with σ0μtzn yields

mδtunδttun+κδtununδt+Hrn=δtunFnγδtun2Qrn,mhδtηnδttηn+Khδtηnμttηnδt+Hhn=δtηnFnΓhδtηn2Qhn,σ0μtznδtznδt+Hbn=σ0μtzngμtzn,vn.

Summing up the above and using relation Equation 27 yields

δt+Hrn+Hhn+Hbn=PnQrnQhnQbn,(33)

where

Hrn=m2δtun2+κ2unun10,Hhn=mh2δtηn2+Kh2μtηn20,Hbn=σ02zn1220,Pn=vbnFn,Qbn=σ0vn2μtzn+gμtzn,vnσ1vnσ0μtzn0,(34)

P is externally supplied power and Qr, Qh, and Qb are dissipation terms with Qr and Qh 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 μtzn and v by vn.

The energy balance in Equation 33 induces the following discrete conservation law (Chatziioannou and van Walstijn, 2015):

En=Hn+1+Δti=0nPi+Qri+Qhi+Qbi=H0,(35)

where Hn=Hrn+Hhn+Hbn. 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 en=EnH0.

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 vb (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 6466 in the Appendix. The energy conservation error en (where, in this case H0=0) is also shown in Figure 3.

Figure 3
www.frontiersin.org

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 1
www.frontiersin.org

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

Table 2
www.frontiersin.org

Table 2. 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.

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 fN 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
www.frontiersin.org

Figure 4. 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 (fs=44100 Hz).

Figure 5
www.frontiersin.org

Figure 5. Nonlinear function G(v) for the refined model (left) and the original Dupont model (right). OS is an oversampling factor, with sampling rate fs=OSfs̄ and fs̄=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 û that is obtained with 1024 times oversampling, as

eg=uû2û2(36)

Figure 6
www.frontiersin.org

Figure 6. 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 ũ(x,t) be real-valued functions defined over an interval [0,L] and for time t0, with an inner product and a norm defined as

u,ũ=0Lux,tũx,tdx,u2=u,u.

Using the subscripts x and t to denote differentiation with respect to space (x) and time (t), respectively, the following identities hold.

tu,u=ddt12u2,(37)
u,xũ=xu,ũ+uL,tũL,tu0,tũ0,t,(38)

where ddt 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):

ρAt2u=Tx2uEIx4u2γ0ρAtu+2γ1ρAtx2uf,(39)
PTt2w=KTx2w2γ2PTtw+rf,(40)
Dht2η=KhηΓhtηf,(41)

where ρ is the material density, A=πr2 is the cross-sectional area of the string with radius r, T is the tension of the string, E is Young’s modulus, I=πr4/4 are the area moment of inertia, and γ0 and γ1 represent frequency independent and frequency dependent damping. In addition, PT denotes the polar moment of inertia, KT torsional stiffness, and γ2 is a torsional damping coefficient. The bow stick is regarded as a rigid frame moving at a given velocity vb and supporting a ribbon of compliant bow-hair of density Dh ([kg/m]) with distributed spring and damping constants Kh and Γh, respectively. The relative bow–string velocity is then expressed as

v=turδtwvbtη.(42)

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

ux,t|x=0,L=0,x2u|x=0,L=0,(43)

and for the torsional movement we assume fixed boundary conditions

wx,t|x=0,L=0.(44)

A fourth equation, needed to close the system, describing the friction force f is

fz,v=σ0z+σ1vtz,(45)

where σ0 and σ1 are now distributed stiffness and distributed damping, respectively, and tz is the time derivative of z. It is related to v through

tz=v1αz,vzzssv,(46)

where the adhesion map α is defined in Equation 6 and zss 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 tu, tw, tη, and σ0z, respectively, and summing the results. Utilizing Equation 37 and the identity in Equation 38 repeatedly, the energy balance follows:

dHdt=QP+B0L,(47)

where

H=Hr+Hw+Hh+Hb,Q=Qr+Qw+Qh+Qb,B=Br+Bw,

and P(t)=vb,f is the power supplied by the bow. The string energy coming from the transverse and torsional motions, Hr and Hw, respectively, bow hair energy Hh, and bristle energy Hb are given by

Hrt=ρA2tu22+T2xu22+EI2x2u220,Hwt=PT2tw22+KT2xw220,Hht=Dh2tη22+Kh2η220,Hbt=σ02z220.

The dissipated energies in the string (Qr and Qw), bow hair (Qh), and bristles (Qb) are given by

Qrt=2γ0ρAtu22+2γ1ρAxtu220,Qwt=2γ2PTtw220,Qht=Γhtη220,Qbt=σ1vv,v+σ0zσ1vv,αz,vzvzssv,

and the boundary terms are

Brt=TtuxuEItux3u+EItxux2u+2γ1ρAtutxu,Bwt=KTtwxw.

Under simply supported boundary conditions, Br vanishes. Similarly, Bw vanishes due to fixed boundary conditions for the torsional movement of the string.

Given that H0, the passivity of the system may be assessed by observing the dissipated energy expressions. More precisely, passivity is guaranteed if Qb0. Let Wb be the width of the bow, then

Qbt=Wbσ1vv2+αz,vvzzssvσ0zσ1vvdt0

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 t0. Let dN={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 uln. Let un=[u0n,,uNn]T. For two vectors un and ũn, the discrete inner product and norm on dN are defined as

un,ũndN=Δxl=0Nulnũln,undN2=un,undN.

Other domains that differ from dN by removing endpoints and will later be used are

dN̄=lΔx:l=1,,N,dN̲=lΔx:l=0,,N1,dN̲̄=lΔx:l=1,,N1.(48)

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:

δtun,undN=δt+12un,un1dN,δtun,δttundN=δt+12δtundN2.

Spatial forward, backward, and central discretization operators are defined as

δx+uln=ul+1nulnΔx,δxuln=ulnul1nΔx,δxuln=ul+1nul1n2Δx,(49)

and can be used to obtain approximations to higher-order partial differential operators:

δxx=δxδx+,δxxxx=δxxδxx.

Like in the continuous case, the following relation can be derived:

un,δxxũndN=δx+un,δx+ũndN̲+uNnδx+ũNnu0nδx+ũ0n.

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

The model is discretized in time and space with functions uln 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/fs (in s) with fs the sampling rate (in Hz) and nN, 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Δt22,(50)

where τ=c2Δt2+4γ1Δt with c=T/ρA are the wave speed and κ=EI/ρA is a stiffness coefficient. The grid points are dN={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ΔxT, where the grid spacing ΔxT (in m) for the torsional movement of the string must satisfy (Bilbao, 2009):

ΔxTcTΔt,(51)

where cT=KT/PT is the torsional wave speed. The grid points are dNT={lΔxT:l=0,,NT}, where NT=L/ΔxT; hence the total number of grid points is NT+1. Torsional waves travel much faster than transverse waves; therefore, the number of grid points NT is much smaller than N. The spatial discretization operators associated with grid dNT will be denoted with an upper T superscript—for example, δx+Twln=wl+1nwlnΔxT.

For a point xm under the bow, the interpolation vectors ixm and iTxm interpolate the string displacement at position xm for the transverse and torsional motion, respectively. ixm is a row vector of size N+1 that multiplies the column vector un=[u0n,,uNn]T, and iTxm is a row vector of size NT+1 that multiplies the column vector wn=[w0n,,wNTn]T. The simplest interpolation is the one of 0th-order where ilxm=1 for l=xm/Δx, and iT,lxm=1 for l=xm/ΔxT, respectively, and zeros elsewhere. For the definition of higher orders of interpolation vectors, see Bilbao (2009). On the other hand, a spreading vector jxm is a column vector that distributes the friction force around the bowing point xm on the grid lΔx. Similarly, jTxm is a spreading vector that distributes the friction force around the bowing point xm on the grid lΔxT. The spreading and interpolation vectors are related through

jxm=1ΔxixmT,jTxm=1ΔxTiTxmT.

To simplify the notation, we first divide Equation 39 by ρA, then discretize it to obtain

δttun=c2δxxunκ2δxxxxun2γ0δtun+2γ1δtδxxunMρA1Jfn,(52)

where J=[jx1||jxM] is an (N+1)×M matrix with mth column being jxm, m=1,,M, and fn=[f1n,,fMn]T is a column vector with M rows where each row describes the friction for a point xm of the string that is in contact with the bow. The friction force is discretized using an interleaved grid, with

fmn=fμtzmn,vmn,wherefμtzmn,vmn=σ0μtzmn+σ1vmngμtzmn,vmn,

with

gμtzmn,vmn=vmn1αμtzmn,vmnμtzmnzssvmn,

for m=1,,M. For simplicity of notation, let gn denote the M×1 column vector with entries gmn

gmn=gμtzmn,vmn,

then

δtzn=gn.(53)

Assuming simply supported ends, the boundary conditions imply

u0n=uNn=0,u1n=u1n,uN+1n=uN1n,(54)

for all nN.

Similarly, by dividing Equation 40 by PT, it is discretized as follows:

δttTwn=cT2δxxTwn2γ2δtwn+rMPTJTfn,(55)

where JT=[jTx1||jTxM] is an (NT+1)×M matrix with the mth column being jTxm, m=1,,M. Assuming fixed ends, the boundary conditions imply

w0n=wNTn=0,(56)

for all nN.

Discretization of the equation governing bow hair displacement is performed as in the lumped case:

Dhδttηn=KhμttηnΓhδtηn1Mfn,(57)

where ηn=[η1n,,ηMn]T. Here, for each point xm under the bow, the bow hair compliance is computed. Then, if the bow velocity at time nΔt is vbn, the relative velocities at points of the string in contact with the bow are discretized as

vn=IδtunrITδtwnvbnδtηn,(58)

where I=[ix1;;ixM], IT=ΔxΔxT[iTx1;;iTxM]T are M×(N+1) and M×(NT+1) matrices with the mth row being ixm and iTxm, respectively, for m=1,,M. Let EM be an M by M identity matrix and

L=ξSMρAIJ+r2ξTMPTITJT+ξHMEM,

where ξS=(2Δt+2γ0)1, ξT=(2Δt+2γ2)1, and ξH=(2ΔtDh+Δt2Kh+Γh)1. Utilizing the expression in Equation 30 for the discrete operators δtt and μt, Equation 58 becomes

vn+vbn=Lfn+sn,(59)

where sn is an M×1 column vector with entries smn=ξSamnξTbmn+ξHcmn, where

amn=ixmc2δxxunκ2δxxxxun+2γ1δtδxxun+2Δtδtun,bmn=ΔxΔxTiTxmrcT2δxxTwn+2rΔtδtwn,cmn=Δt2Kh+2ΔtDhηmnηmn1ΔtKhηmn.

In order to update the system variables, the vectors vn and μtzn must first be computed. A system of 2M equations is formed using Equation 59 and relation in Equation 53 together with δtzn=2Δtμtznzn12.

vn+Lfnsn+vbn=0,(60)
gn2Δtμtznzn12=0.(61)

The system can then be solved for vn and μtzn using an iterative solver. Once vn and μtzn are known, un+1, wn+1, and ηn+1 can be updated. First, the variables related to the friction force and the friction force itself are computed,

zn+12=2μtznzn12,gn=zn+12zn12Δt,fn=σ0μtzn+Ξ1ngn,

where Ξ1n is an M×M diagonal matrix with diagonal terms being {σ1(vmn)}m=1M. Then, the string and bow hair variables are updated as

un+1=2ΔtξSAun+Bun12ΔtξSMρAJfn,wn+1=2ΔtξTCwn1γ2Δtwn1+r2ΔtξTMPTJTfn,ηn+1=D2DhΔt2Kh2ηnDDhΔt2+Kh4Γh2Δtηn1DMfn,

where

A=2E+2γ1Δt+c2Δt2δxxκ2Δt2δxxxx,B=γ0Δt1E2γ1Δtδxx,C=2ET+cT2Δt2δxxT,D=DhΔt2+Kh4+Γh2Δt1,

with E an (N+1)×(N+1) identity matrix and ET an (NT+1)×(NT+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 μtzn as

fn=σ0μtzn+2ΔtΞ1nμtznzn12,

the bristle displacement can be expressed as

μtzn=ΞnL1vn+vbn+L1sn+2ΔtΞ1nzn12,(62)

where Ξn is an M×M diagonal matrix with entries (σ0+2Δtσ1(vmn))1, m=1,,M on the diagonal. Plugging (62) into (61) and gn, only Equation 61 must be solved for vn. zn+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 δtun, an inner product of Equation 55 with δtwn, an inner product of Equation 57 with δtηn, and an inner product of σ0Mμtzn with δtzn. Using summation by parts identities as well as boundary conditions leads to

δt+Hn=QnPn+Bn,(63)

where H, the total numerical energy, is defined as H=Hr+Hw+Hh+Hb and, assuming that the stability conditions in Equations 50 and 51 are satisfied,

Hrn=ρA2δtundN2+T2δx+un,δx+un1dN̲+EI2δxxun,δxxun1dN̲̄0,Hwn=PT2ΔxΔxTδtwndNT2+KT2ΔxΔxTδx+Twn,δx+Twn1dNT̲0,Hhn=Dh2δtηnTδtηn+Kh2μtηnTμtηn0,Hbn=1Mσ02zn12Tzn120.

The total energy lost due to damping is defined as Q=Qr+Qw+Qh+Qb with

Qrn=2γ0ρAδtundN22γ1ρAδtun,δtδxxundN0,Qwn=2γ2PTΔxΔxTδtwndNT20,Qhn=ΓhδtηnTδtηn0,Qbn=1MΞ1nvnTvn+σ0μtznTonΞ1nvnTon0,

where on is a vector with entries omn=α(μtzmn,vmn)μtzmnzss(vmn) for m=1,,M. The non-negativity of the dissipation term Qbn can be shown similarly to the lumped case. The energy supplied by the bow P and boundary terms B=Br+Bw is

Pn=1Mvbnm=1Mfmn,Brn=TδtuNnδx+uNnδtu0nδx+u0n+EIδtuNnδx+δxxuNnδtu0nδx+δxxu0nEIδtδx+uNnδxxuNnδtδx+u0nδxxu0n,Bwn=KTδtwNnδx+TwNnδtw0nδx+Tw0n.

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:

Qbn=1Mm=1Mσ1vmnvmn2+αμtzmn,vmnvmnμtzmnzssvmnσ0μtzmnσ1vmnvmn0.

Non-negativity of Qbn 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 vb is reached. The physical model parameters used to generate these signals are given in Table 2.

Figure 7
www.frontiersin.org

Figure 7. 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.

Figure 8
www.frontiersin.org

Figure 8. 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% (fN=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., and Valls, C. (2012). Ordinary differential equations: qualitative theory. America, American Mathematical Society.

Google Scholar

Bilbao, S. (2009). Numerical sound synthesis. J. Wiley and Sons. doi:10.1002/9780470749012

CrossRef Full Text | Google Scholar

Bilbao, S., Torin, A., and Chatziioannou, V. (2015). Numerical modeling of collisions in musical instruments. Acta Acust. United Acust. 101, 155–173. doi:10.3813/AAA.918813

CrossRef Full Text | Google Scholar

Chabassier, J., and 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. doi:10.1016/j.cma.2010.04.013

CrossRef Full Text | Google Scholar

Chatziioannou, V., and 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.

Google Scholar

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.

Google Scholar

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

Google Scholar

Deuflhard, P. (2011). Newton methods for nonlinear problems: affine invariance and adaptive algorithms. Springer Sci. and Bus. Media 35.

Google Scholar

de Wit, C., Olsson, H., Åström, K., and Lischinsky, P. (2024). A new model for control of systems with friction. IEEE Trans. Autom. Control 40, 419–425. doi:10.1109/9.376053

CrossRef Full Text | Google Scholar

Ducceschi, M., and Bilbao, S. (2022). Simulation of the geometrically exact nonlinear string via energy quadratisation. J. Sound Vib. 534, 117021. doi:10.1016/j.jsv.2022.117021

CrossRef Full Text | Google Scholar

Dupont, P., Armstrong, B., and Hayward, V. (2000). Elasto-plastic friction model: contact compliance and stiction. Proc. 2000 Acc. IEEE Chic., 1072–1077 vol.2. doi:10.1109/ACC.2000.876665

CrossRef Full Text | Google Scholar

Dupont, P., Hayward, V., Armstrong, B., and Altpeter, F. (2002). Single state elasto-plastic friction models. IEEE Trans. Autom. Control 47, 787–792. doi:10.1109/TAC.2002.1000274

CrossRef Full Text | Google Scholar

Falaize, A., and 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. doi:10.1007/s11071-024-10438-9

CrossRef Full Text | Google Scholar

Galluzzo, P. (2004). On the playability of stringed instruments. Ph.D. thesis. doi:10.17863/CAM.14046

CrossRef Full Text | Google Scholar

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

Google Scholar

Hueso, J., Martínez, E., and Torregrosa, J. (2009). Modified Newton’s method for systems of nonlinear equations with singular Jacobian. J. Comput. Appl. Math. 224, 77–83. doi:10.1016/j.cam.2008.04.013

CrossRef Full Text | Google Scholar

Lampis, A., Mayer, A., and Chatziioannou, V. (2024). Assessing playability limits of bowed-string transients using experimental measurements. Acta Acust. 8, 44. doi:10.1051/aacus/2024034

CrossRef Full Text | Google Scholar

Maestre, E., Spa, C., and Smith, J. (2014). A bowed string physical model including finite-width thermal friction and hair dynamics. ICMC.

Google Scholar

Matusiak, E., and Chatziioannou, V. (2024). Elasto-plastic friction modeling toward reconstructing measured bowed-string transients. J. Acoust. Soc. Am. 156, 1135–1147. doi:10.1121/10.0028228

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

Pitteroff, R., and 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.

Google Scholar

Pitteroff, R., and 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.

Google Scholar

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

Google Scholar

Serafin, S., Avanzini, F., and 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.

Google Scholar

Smith, J., and Woodhouse, J. (1999). The tribology of rosin. J. Mech. Phys. Solids 48, 1633–1681. doi:10.1016/S0022-5096(99)00067-8

CrossRef Full Text | Google Scholar

van Walstijn, M., Chatziioannou, V., and 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. doi:10.1016/j.jsv.2023.117968

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Willemsen, S., Bilbao, S., and 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.

Google Scholar

Woodhouse, J. (2003). Bowed string simulation using a thermal friction model. Acta Acustica united Acustica 89, 355–368. doi:10.25144/18216

CrossRef Full Text | Google Scholar

Appendix

The modal expansion for the displacement of the string, assuming simply supported boundary conditions, is

ux,t=i=1NMqitψix,whereψix=siniπLx,

where NM 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,

ρAt2u=Tx2uEIx4u2γ0ρAtu+2γ1ρAtx2u,

and an inner product with ψn is taken. Since ψi,ψn=L2δi,n, we obtain

mq̈n=κnqnγnq̇n,

where

m=ρAL2,(64)
κn=nπ22LT+EInπL2,(65)
γn=ρALγ0+γ1nπL2.(66)

Therefore, by setting NM=1, the dynamics of the system reduce to a damped harmonic oscillator.

Keywords: bowed mass, bow–string interaction, friction, passivity, finite differences, energy methods, numerical stability, elasto-plastic

Citation: Matusiak E, Chatziioannou V and Van Walstijn M (2025) Numerical modelling of elasto-plastic friction in bow–string interaction with guaranteed passivity. Front. Signal Process. 5:1525044. doi: 10.3389/frsip.2025.1525044

Received: 08 November 2024; Accepted: 11 March 2025;
Published: 21 May 2025.

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

Copyright © 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.

*Correspondence: Ewa Matusiak, bWF0dXNpYWtAbWR3LmFjLmF0

Disclaimer: 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.