# Nonsmooth contact dynamics for the numerical simulation of collisions in musical string instruments

Clara Issanchou, Vincent Acary, Franck Pérignon, Cyril Touzé, Jean-Loic Le Carrou

► To cite this version:

Clara Issanchou, Vincent Acary, Franck Pérignon, Cyril Touzé, Jean-Loic Le Carrou. Nonsmooth contact dynamics for the numerical simulation of collisions in musical string instruments. Journal of the Acoustical Society of America, 2018, 143 (5), pp.1-13. 10.1121/1.5039740. hal-01806133v2

HAL Id: hal-01806133

https://ensta-paris.hal.science/hal-01806133v2

Submitted on 15 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonsmooth contact dynamics for the numerical simulation of collisions in musical string instruments

Clara Issanchou,<fnref n="1" />, a) Vincent Acary,<fnref n="2" /> Franck Pérignon,<fnref n="2" /> Cyril Touzé,<fnref n="3" />, b and Jean-Loïc Le Carrou<fnref n="4" />

Collisions in musical string instruments play a fundamental role in explaining the sound production in various instruments such as sitars, tanpuras and electric basses. Contacts occurring during the vibration provide a nonlinear effect which shapes a specific tone due to energy transfers and enriches the hearing experience. As such, they must be carefully simulated for the purpose of physically-based sound synthesis. Most of the numerical methods presented in the literature rely on a compliant modeling of the contact force between the string and the obstacle. In this contribution, numerical methods from nonsmooth contact dynamics are used to integrate the problem in time. A Moreau-Jean time-stepping scheme is combined with an exact scheme for phases with no contact, thus controlling the numerical dispersion. Results for a two-point bridge mimicking a tanpura and an electric bass are presented, showing the ability of the method to deal efficiently with such problems while invoking, as compared to a compliant approach, less modelling parameters and a reduced computational burden.

©2018 Acoustical Society of America. [http://dx.doi.org(DOI number)]

[XYZ]

Pages: 1–13

## I. INTRODUCTION

Collisions are of prime importance in musical acoustics for explaining the particular timbre of a number of instruments ranging from strings (a typical example being that of indian instruments such as sitar, tanpura and veena) to drums (e.g. snare drum) (Bilbao, 2012; Fletcher and Rossing, 1998; Raman, 1921). In these examples, the role of contacts is to alter the frequency content due to a nonlinear, nonsmooth interaction that generates high frequencies and contributes to enrich the hearing experience. This effect is particularly prominent in the case of the sitar, where a curved bridge contributes to significantly modify the frequency content of the string vibration, see e.g. (Bilbao et al., 2015; Mandal and Wahi, 2015; Siddiq, 2012; Vyarasayani et al., 2009), and in the case of the tanpura and its particular bridge (Chatziioannou and van Walstijn, 2015; Issanchou et al., 2017; Valette et al., 1991). Other examples that attracted interest in the recent years concern the case of string/frets interactions (for e.g. guitar or electric bass) (Bilbao

and Torin, 2015; Issanchou et al., 2018; Trautmann and Rabenstein, 2004) and snare drum where metal wires are in contact with a vibrating membrane (Bilbao, 2012; Bilbao et al., 2015).

Earlier studies on vibrating strings with contacts derive a number of analytical results, mostly in the 1980s, see e.g. (Amerio, 1978; Cabannes, 1987; Citrini, 1991; Schatzman, 1980). To overcome the limitations of these approaches, where the string needs to be perfect without stiffness, thus discarding the dispersive effect which has been shown to be of prime importance in the case of the sitar and tanpura (Chatziioannou and van Walstijn, 2015; Issanchou et al., 2017; Siddiq, 2012), recent research efforts concentrate toward the development of efficient, robust and accurate numerical methods in order to simulate musical strings encountering an obstacle during their vibration. Most of the methods presented in the last years use a regularisation in order to treat numerically the contact force, see e.g. the energy-conserving schemes proposed by Bilbao et al. (Bilbao et al., 2015; Desvages and Bilbao, 2015; Ducceschi et al., 2016) and by van Walstijn et al. (Chatziioannou and van Walstijn, 2015; van Walstijn and Bridges, 2016; van Walstijn et al., 2016), the modal approach proposed in (Issanchou et al., 2017) or the approach followed in (Inácio et al., 2006) to model the interaction between a *puja* (exciting stick) and a Tibetan bowl. In all these studies the contact

<footnote n="1">¹Sorbonne Université, CNRS, quipe LAM, Institut Jean Le Rond d'Alembert, 4 place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="2">²Université Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 655, Avenue de l'Europe, 38334 Saint Ismier Cedex, France</footnote>

<footnote n="3">³IMSIA, ENSTA ParisTech-CNRS-EDF-CEA, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France</footnote>

<footnote n="4">⁴Sorbonne Université, CNRS, Equipe LAM, Institut Jean Le Rond d'Alembert, 4 place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote>aAlso at: IMSIA, ENSTA ParisTech-CNRS-EDF-CEA, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France</footnote>

<footnote>b) cyril.touze@ensta-paristech.fr; Corresponding author.</footnote>
force is modeled using a power-law method with two pa-
rameters defining the stiffness of the repelling force, thus
allowing to cover a wide range of contact laws, from soft
collisions, see e.g. contact between felt and hammer in
piano (Boutillon, 1988), between mallet and membrane
in a kettledrum (Rhaouti et al., 1999) or between finger
and fretboard (Bilbao and Torin, 2015), to hard contacts.
In this article, we will name as compliant approach the
methods using a regularisation to express the repelling
force. In this framework, power-law are often used but
other functions may also be selected. Also, the parame-
ters defining the contact force may have a physical basis,
see e.g. (Goldsmith, 2001) for some examples. On the
other hand, ad hoc values may be used as numerical free
parameters, so that a penalty approach is at hand.

On the other hand, a large body of research has
been dedicated to the development of nonsmooth nu-
merical methods in order to deal efficiently with nu-
merical challenges posed by contact and friction forces.
These methods rely on specific assumptions (e.g. no in-
terpenetration is allowed between the contacting bod-
ies) and use mathematical tools from the measure the-
ory, differential inclusions and complementarity systems.
The first developments have been pioneered by Jean and
Moreau, see e.g. (Jean, 1999; Jean and Moreau, 1987),
continued by numerous investigations (Doyen et al., 2011;
Janin and Lamarque, 2001; Paoli and Schatzman, 2002),
and are now summarized in reference books (Acary and
Brogliato, 2008; Studer, 2009). Nonsmooth methods
have been successfully applied in a variety of contexts
ranging from granular media (Renouf et al., 2004), geo-
materials (Jean, 1995), multibody dynamics (Chen et al.,
2013) to realistic simulations of hair motions and living
systems (Acary et al., 2014; Bertails-Descoubes et al.,
2011). In the field of vibration, the method has been ap-
plied to rotor/casing contacts in (Meingast et al., 2014)
as well as to string vibrations in (Ahn, 2007), where the
study was however limited to the case of a perfect string
without stiffness and a frictionless contact. In the area of
musical acoustics, the action of a grand piano has been
recently simulated efficiently by using a nonsmooth ap-
proach (Thorin et al., 2017).

As remarked by a number of investigators, numer-
ical integration for contact dynamics is generally time-
consuming due to the high-frequency content gener-
ated (Doyen et al., 2011; Issanchou et al., 2017), leading
to consider very small time steps in order to achieve con-
vergence. In this context, nonsmooth numerical meth-
ods, as avoiding the costly step of finding the zeros of
a nonlinear function with a Newton-Raphson approach
and using efficient numerical methods to solve a linear
complementarity problem, should decrease the compu-
tational burden (Acary and Brogliato, 2008). On the
modeling point of view, nonsmooth methods are partic-
ularly appealing when one only needs an efficient nu-
merical method to repel a vibrating structure from an
obstacle with a simple representation of the dissipation
phenomena at contact. As such, nonsmooth methods de-
scribe the contact law with a single parameter: the coef-

ficient of restitution, instead of the two parameters used
in the power-law approach as done in e.g. (Bilbao et al.,
2015; Chatziioannou and van Walstijn, 2015; Issanchou
et al., 2017) in the conservative case, and additional pa-
rameters to model the dissipation following a Hunt and
Crossley approach. Let us note that it is also difficult
with the power-law approach to obtain large dissipation
rates at contact. For this purpose, more complicated
compliant models must be introduced including plastic-
ity effects (Nguyen and Brogliato, 2014). Less modelling
parameters could be seen as an advantage, but also as a
drawback if one prefers to have more degrees of freedom
in order to represent different contact laws with a large
variety of stiffness and/or damping to account for a phys-
ical reality having perceptual effects. A typical example
is that of the guitar where very different parameters are
used for the string/fret and finger/fingerboard collisions
(Bilbao and Torin, 2015). Note however that nonsmooth
methods can also be conjugated with a material descrip-
tion of a soft structure including nonlinear stiffness and
damping to account for these effects and tune them at
ease, see e.g. the recent modeling of felt in piano action
proposed in (Thorin et al., 2017).

In this contribution, a nonsmooth numerical ap-
proach based on a Moreau-Jean time-stepping scheme is
adapted to the case of a string vibrating against a stiff
obstacle with hard contacts, for the specific purpose of
musical acoustics. The examples have been purposely
chosen to restrict our discussion to hard contacts where
the nonsmooth method should be an interesting alterna-
tive to compliant approach. Standard nonsmooth meth-
ods are known to produce numerical dispersion (Yoong
et al., 2018), which is a specific issue in the field of mu-
sical acoustics. During non-contacting phases, an exact
scheme, introduced in (Bilbao, 2009) and used for a com-
pliant contact approach in (Issanchou et al., 2017; van
Walstijn and Bridges, 2016), is considered in order to
ensure a numerical integration method with controlled
dispersion. The efficiency of the scheme is demonstrated
on two examples which are compared either to a com-
pliant approach or to experimental results. First, the
case of a two-point bridge mimicking a tanpura is shown.
Results are then extended to a fretted electric bass hav-
ing at most 20 contact points. The accuracy of the two
different modeling options (nonsmooth vs compliant ap-
proach) are discussed and the computational burdens of
these methods are compared.

II. MODEL

A. Vibrating string against a unilateral obstacle

We consider a stiff string of length *L* (m), mass
per unit length *μ* (kg·m⁻¹) and tension *T* (N). Its
Young’s modulus *E* (Pa) and moment of inertia *I* de-
fine its stiffness. The string vibrates against a unilateral
obstacle, the profile of which is described by *g*(*x*) (see
Fig. 1). Eq. (1) describes the displacement of the string
<figure><img src="image_2.png" /><figcaption>FIG. 1. Scheme of a string vibrating against a unilateral obstacle.</figcaption></figure>

along (Oz):

$$ \mu u_{tt} - T u_{xx} + EI u_{xxxx} = f(x, t), \quad (1) $$

where the subscript t (respectively x) refers to a partial derivative with respect to time (respectively space). The right-hand side f(x, t) refers to the contact force per unit length.

The string is simply supported at its endpoints, so that $\forall t \in \mathbb{R}^+$: $u(0, t) = u(L, t) = u_{xx}(0, t) = u_{xx}(L, t) = 0$. In order to achieve a fine representation of eigenfrequencies and damping parameters, we employ the following modal description of the string:

$$ u(x, t) = \sum_{j=1}^{N_m} q_j(t) \phi_j(x), \quad (2) $$

where $\phi_j(x) = \sqrt{\frac{2}{L}} \sin\left(\frac{j\pi x}{L}\right)$ and $q_j$ is the jth modal amplitude.

Inserting the modal expansion (2) into the equation of motion (1), multiplying by another mode shape and integrating over the string length, and then adding losses to the final modal equations, one obtains:

$$ \mu(\ddot{\mathbf{q}} + \Omega^2 \mathbf{q} + 2\mathbf{Y}\dot{\mathbf{q}}) = \mathbf{F}, \quad (3) $$

where $\mathbf{q} = [q_1, q_2, \dots, q_{N_m}]^T$ is the vector of unknown modal amplitudes. $\Omega$ and $\mathbf{Y}$ are diagonal matrices such that $\Omega$ contains radian frequencies $\Omega_{jj} = \omega_j = 2\pi\nu_j$, $\nu_j$ being the jth eigenfrequency, and $\mathbf{Y}$ contains damping coefficients $\mathbf{Y}_{jj} = \sigma_j$. These quantities follow a model presented in (Issanchou et al., 2017; 2018; Paté et al., 2014; Valette and Cuesta, 1993). We briefly recall the expression of eigenfrequencies and damping parameters in Appendix. The vector $\mathbf{F}$ contains modal forces and represents the projection of the contact force onto mode shapes, its entries read $F_p = \int_0^l f(x, t)\phi_p(x)dx$, for $p \in \{1, \dots, N_m\}$.

## B. Contact force

We choose to model the contact condition by a unilateral constraint in order to avoid the interpenetration of the solids. Defining $\eta(x, t) = u(x, t) - g(x)$ as the distance between the string and the obstacle, this constraint is given by:

$$ \eta(x, t) \ge 0, \forall t \in \mathbb{R}^+. \quad (4) $$

In order to satisfy Eq. (4), the reaction force introduced in Eq. (1) has to be positive and vanishes only if the constraint is not active (open contact). The complete model is the Signorini law which can be simply written with the following notation (Acary and Brogliato, 2008; Signorini, 1933):

$$ 0 \le \eta(x, t) \perp f(x, t) \ge 0, \quad (5) $$

meaning that either $\eta(x, t) = 0$ and $f(x, t) \ge 0$, or $f(x, t) = 0$ and $\eta(x, t) \ge 0$.

Second order dynamics with unilateral constraints involves velocity jumps. The usual setting is to consider that the velocity map $t \mapsto u_t(x, t)$ is a right continuous function, i.e. $u_t^+(x, t) = u_t(x, t)$, and of bounded variation. The notation $u_t^+(x, t)$ (respectively $u_t^-(x, t)$) stands for the right limit (respectively the left limit) of $t \mapsto u_t(x, t)$ at time $t$, i.e. $u_t^+(x, t) = \lim_{s \to t} u_t(x, s)$ (respectively $u_t^-(x, t) = \lim_{s \to t} u_t(x, s)$). Since the function is of bounded variation, these limits exist.

In order to define the solution, especially in a discrete (finite-dimensional) system, a condition on the velocity after an impact must be specified. To this purpose, we choose the Newton impact law (Acary, 2016):

$$ \dot{\eta}^+(x, t) = -\rho \dot{\eta}^-(x, t) \quad \text{if} \quad \eta(x, t) = 0, \quad (6) $$

where the coefficient of restitution $\rho \in [0, 1]$ defines the string behaviour at impact instants.

In the sequel, the Signorini law is formulated at the velocity level, which allows one to explicitly control contact losses in taking into account the Signorini condition together with the impact law. Therefore, the contact law writes:

$$ \begin{cases} 0 \le \dot{\eta}^+(x, t) + \rho \dot{\eta}^-(x, t) \perp f(x, t) \ge 0 & \text{if } \eta(x, t) = 0 \\ f(x, t) = 0 & \text{otherwise.} \end{cases} \quad (7) $$

The viability Lemma of Moreau (Moreau, 1999) ensures that the condition at the velocity level (7) implies the condition at the position level (5) if the constraint (4) is satisfied at the initial time.

Since the system has some discontinuities in the velocity $\mathbf{w} = \dot{\mathbf{q}}$, the acceleration $\dot{\mathbf{w}}$ is not everywhere defined in the classical sense. A differential measure $d\mathbf{w}$ is associated with the velocity $\mathbf{w}$ and plays the role of the acceleration. If we assume furthermore that $\mathbf{w}$ is of special bounded variation (i.e. $\mathbf{w}$ may be decomposed into a sum of an absolutely continuous function and a jump function), the differential measure can be decomposed with respect to the Lebesgue measure $dt$ as:

$$ d\mathbf{w} = \ddot{\mathbf{q}}dt + (\mathbf{w}^+ - \mathbf{w}^-)d\tilde{\nu}, \quad (8) $$

where $d\tilde{\nu}$ is a discrete measure of the form $\sum_i a_i \delta_{s_i}$ with given sequences $\{a_i\}$ and $\{s_i\}$ of real numbers. The notation $\delta_s$ refers to the Dirac measure supported at time $s$ (see e.g. (Moreau, 1988) for details on differential measures). In other words, $\dot{\mathbf{w}}(t) = \ddot{\mathbf{q}}(t)$ almost everywhere
and at impact instants $t^*$ we have $d\mathbf{w} = (\mathbf{w}^+ - \mathbf{w}^-)\delta_{t^*}$. Similarly, the reaction force in the modal space is defined by a vector measure:

$$d\mathbf{I} = \mathbf{F}dt + \mathbf{P}d\tilde{\nu}, \qquad (9)$$

where the vector $\mathbf{F}$ corresponds to a modal continuous contact force and the vector $\mathbf{P}$ is a modal contact impulse corresponding to velocity jumps. In terms of differential measure, the modal equations read as:

$$\mu(\mathbf{d}\mathbf{w} + \mathbf{\Omega}^2\mathbf{q}dt + 2\mathbf{\Upsilon}\dot{\mathbf{q}}dt) = d\mathbf{I}. \qquad (10)$$

By substituting (8) and (9) in (10), we remark that Eq. (3) is satisfied $dt$-almost everywhere. At the instants of discontinuities, we obtain the impact equation:

$$\mu(\mathbf{w}^+ - \mathbf{w}^-) = \mathbf{P}. \qquad (11)$$

The contact condition at the velocity level is also reformulated in terms of measures as follows:

$$\left\{ \begin{array}{ll} 0 \le \dot{\eta}^{+}(x,t) + \varrho\dot{\eta}^{-}(x,t) \perp di \ge 0 & \text{if } \eta(x,t) = 0 \\ di = 0 & \text{otherwise,} \end{array} \right. \qquad (12)$$

where $di$ is the reaction force in the physical space. The relation between $\eta$, $\dot{\eta}$ and $\mathbf{q}$, $\dot{\mathbf{q}}$ is given by the relations between the quantities in the physical space and the modal space. Let us introduce the column vector $\phi(x)$ containing the $N_m$ first string modes, defined by $\phi_j = \phi_j(x)$, $\forall j \in \{1, ..., N_m\}$. From (2), we get that $u(x,t) = \phi^T(x)\mathbf{q}(t)$ where $\mathbf{q} = [q_1(t), ..., q_{N_m}(t)]^T$ and therefore:

$$\eta(x,t) = \phi^T(x)\mathbf{q}(t) - g(x) \text{ and } \dot{\eta}(x,t) = \phi^T(x)\dot{\mathbf{q}}(t). \qquad (13)$$

By duality, we also have:

$$dI_j = \int_0^L \phi_j di \, dx. \qquad (14)$$

Altogether, the dynamics is given by the following measure differential complementarity problem:

$$\left\{ \begin{array}{ll} \dot{\mathbf{q}} = \mathbf{w} & (15a) \\ \mu(\mathbf{d}\mathbf{w} + \mathbf{\Omega}^2\mathbf{q}dt + 2\mathbf{\Upsilon}\dot{\mathbf{q}}dt) = d\mathbf{I} & (15b) \\ \mathbf{q}(0) = \mathbf{q}_0, \mathbf{w}(0) = \mathbf{w}_0 & (15c) \\ \eta = \phi^T\mathbf{q} - g & (15d) \\ \dot{\eta} = \phi^T\dot{\mathbf{q}} & (15e) \\ dI_j = \int_0^L \phi_j di \, dx & (15f) \\ 0 \le \dot{\eta}^+ + \varrho\dot{\eta}^- \perp di \ge 0 \text{ if } \eta \le 0. & (15g) \end{array} \right.$$

### III. NUMERICAL METHOD

In this section, we present a discretisation of the problem presented in Section II which is able to handle

nonsmooth contacts. The distinctive feature of the proposed scheme is to combine an exact method for the linear (non-contacting) part of the equations of motion with a Moreau-Jean time-stepping approach to handle impulses and velocity jumps. The resulting scheme thus prevents dispersion when the string freely vibrates, which represents an improvement of already existing nonsmooth time-stepping methods.

In the rest of the paper, the variable $t^n$ denotes the discrete time $t^n = n\Delta t$, where $\Delta t$ is the time step. The spatial grid is defined by $x_i = i\Delta x$, $\forall i \in \{0, ..., N\}$, where $\Delta x = \frac{L}{N}$ is the spatial step. Due to selected boundary conditions, $u(x_0, t) = 0$ and $u(x_N, t) = 0 \forall t \in \mathbb{R}^+$, so that only the interior points of the grid are considered.

As prescribed in a Moreau-Jean scheme, Eq. (15b) is integrated over $(t^n, t^{n+1}]$. One thus obtains:

$$\mu(\mathbf{w}^+(t^{n+1}) - \mathbf{w}^+(t^n)) + \mu\mathbf{\Omega}^2 \int_{t^n}^{t^{n+1}} \mathbf{q}dt \\ + 2\mu\mathbf{\Upsilon} \int_{t^n}^{t^{n+1}} \dot{\mathbf{q}}dt = \int_{(t^n, t^{n+1}]} d\mathbf{I}. \qquad (16)$$

The collision term on the right-hand side of Eq. (16) is treated with the Moreau-Jean scheme, leading to the following discrete approximation:

$$P_i^{n+1} \approx \int_{(t^n, t^{n+1}]} dI_i, \qquad (17)$$

where $\mathbf{P}^{n+1} = [P_1^{n+1}, ..., P_{N-1}^{n+1}]^T$ is a contact impulse over the time interval, which will be defined later through a complementarity condition. The Moreau-Jean scheme is generally used with a $\theta$-method to approximate the stiffness and damping terms in Eq. (16), as prescribed for example in (Acary and Brogliato, 2008; Jean, 1999; Jean and Moreau, 1987; Moreau, 1999). However, for mechanical vibratory systems expressed in the modal basis and with application to musical acoustics where the problem of numerical dispersion is particularly stringent, the exact scheme used in (Bilbao, 2009; Issanchou et al., 2017) shall be considered in order to improve the discretisation of the stiffness and damping terms. Using the following notations for discrete approximations, $w_i^n \approx w_i^+(t^n)$, and $q_i^n \approx q_i(t^n)$, it reads:

$$\int_{t^n}^{t^{n+1}} q_i dt \approx \frac{\Delta t q_i^n}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t}, \qquad (18)$$

$$\int_{t^n}^{t^{n+1}} \sigma_i \dot{q}_i dt \approx \frac{\sigma_i^* \Delta t w_i^n}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t}, \qquad (19)$$
where the introduced terms $\gamma_i$ and $\sigma_i^*$ are such that:

$$ \gamma_i = \frac{2}{\omega_i^2 \Delta t^2} - \frac{A_i}{1 + e_i - A_i}, \quad (20) $$

$$ \sigma_i^* = \left( \frac{1}{\Delta t} + \frac{\omega_i^2 \Delta t}{2} - \gamma_i \frac{\omega_i^2 \Delta t}{2} \right) \frac{1 - e_i}{1 + e_i}, \quad (21) $$

with:

$$ A_i = e^{-\sigma_i \Delta t} \left( e^{\sqrt{\sigma_i^2 - \omega_i^2} \Delta t} + e^{-\sqrt{\sigma_i^2 - \omega_i^2} \Delta t} \right), \quad (22) $$

$$ e_i = e^{-2\sigma_i \Delta t}. \quad (23) $$

Finally, one obtains the following time-stepping scheme from (16) for the update of the modal velocity vector using (18) and (19):

$$ w_i^{n+1} - w_i^n + \frac{\Delta t \omega_i^2}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t} q_i^n + \frac{2 \Delta t \sigma_i^*}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t} w_i^n = \frac{1}{\mu} P_i^{n+1}, \quad (24) $$

In matrix form, it writes:

$$ \mathbf{w}^{n+1} - \mathbf{w}^n + \Delta t \mathbf{C} \mathbf{q}^n + \Delta t \mathbf{\check{C}} \mathbf{w}^n = \frac{1}{\mu} \mathbf{P}^{n+1}, \quad (25) $$

where $\mathbf{C}$ and $\mathbf{\check{C}}$ are diagonal matrices with entries:

$$ C_{ii} = \frac{\omega_i^2}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t}, $$

$$ \check{C}_{ii} = \frac{2\sigma_i^*}{1 + (1 - \gamma_i) \frac{\omega_i^2 \Delta t^2}{2} + \sigma_i^* \Delta t}. $$

The update of the modal amplitudes is given by the discretisation of (15)(a) as:

$$ q_i^{n+1} = q_i^n + \Delta t w_i^{n+1} \quad (26) $$

We introduce a matrix $\mathbf{S}$ containing the $N_m$ first string modes, defined by $S_{ij} = \phi_j(x_i), \forall (i, j) \in \{1, ..., N-1\} \times \{1, ..., N_m\}$. Defining the physical discrete displacement vector $\mathbf{u} = [u(x_1, t), ..., u(x_{N-1}, t)]^T$, one easily obtains the following relationship: $\mathbf{u} = \mathbf{S}\mathbf{q}$ (Issanchou et al., 2017). Denoting the displacement of the string at contact points by $\mathbf{u}_c = \mathbf{S}_c\mathbf{q}$, its velocity by $\mathbf{v}_c = \mathbf{S}_c\mathbf{w}$ and introducing the velocity without contact $\mathbf{w}_{\text{free}}^{n+1}$ as:

$$ \mathbf{w}_{\text{free}}^{n+1} = \mathbf{w}^n - \Delta t \mathbf{C} \mathbf{q}^n - \Delta t \mathbf{\check{C}} \mathbf{w}^n, \quad (27) $$

one obtains the following linear equation between $\mathbf{v}_c^{n+1}$ and $\mathbf{p}_c^{n+1}$, with $\mathbf{P}^n = \Delta x \mathbf{S}_c^T \mathbf{p}_c^n$:

$$ \mathbf{v}_c^{n+1} = \mathbf{S}_c \mathbf{w}_{\text{free}}^{n+1} + \mathbf{W}_{cc} \mathbf{p}_c^{n+1}, \quad (28) $$

where the Delassus' matrix $\mathbf{W}_{cc}$ is computed as:

$$ \mathbf{W}_{cc} = \frac{\Delta x}{\mu} \mathbf{S}_c \mathbf{S}_c^T. \quad (29) $$

The complementarity condition (15g) is discretised in a fully implicit way as prescribed for the Moreau-Jean scheme. For a contact point indexed by $c_i$, one has to solve:

$$ \left\{ \begin{array}{ll} \mathbf{p}_{c_i}^{n+1} = 0 & \text{if } \eta_{c_i}^n > 0 \\ 0 \le \mathbf{v}_{c_i}^{n+1} + \varrho \mathbf{v}_{c_i}^n \perp \mathbf{p}_{c_i}^{n+1} \ge 0 & \text{if } \eta_{c_i}^n \le 0. \end{array} \right. \quad (30) $$

Let us denote the index set of active contact points as $\bar{c} = \{c_i | \eta_{c_i}^n \le 0\}$. The impulse $\mathbf{p}_{\bar{c}}^{n+1}$ is the solution of:

$$ \left\{ \begin{array}{l} \mathbf{v}_{\bar{c}}^{n+1} = \mathbf{S}_{\bar{c}} \mathbf{w}_{\text{free}}^{n+1} + \mathbf{W}_{\bar{c}\bar{c}} \mathbf{p}_{\bar{c}}^{n+1} \\ 0 \le \mathbf{v}_{\bar{c}}^{n+1} + \varrho \mathbf{v}_{\bar{c}}^n \perp \mathbf{p}_{\bar{c}}^{n+1} \ge 0 \end{array} \right. \quad (31) $$

which is a Linear Complementarity Problem (LCP). A LCP consists in finding the vectors $\mathbf{z}, \mathbf{\lambda}$ such that:

$$ \left\{ \begin{array}{l} \mathbf{z} = \mathbf{W}\mathbf{\lambda} + \mathbf{a} \\ 0 \le \mathbf{z} \perp \mathbf{\lambda} \ge 0 \end{array} \right. \quad (32) $$

for some given matrix $\mathbf{W}$ and given vector $\mathbf{a}$. Let us recall that the LCP has a unique solution for all $\mathbf{a}$ if the matrix $\mathbf{W}$ is a positive definite matrix (Cottle et al., 1992). It is straightforward to identify in (32) the vector $\mathbf{a}$ as $\mathbf{a} = \mathbf{S}_{\bar{c}}\mathbf{w}_{\text{free}}^{n+1} + \varrho \mathbf{v}_{\bar{c}}^n$ and the matrix $\mathbf{W} = \mathbf{W}_{\bar{c}\bar{c}}$. In our application, the Delassus matrix $\mathbf{W}$ is obviously symmetric positive definite and a unique solution is ensured. If the string vibrates against a single point obstacle, the Delassus matrix reduces to a positive scalar and the LCP can be solved directly as:

$$ p_{\bar{c}}^{n+1} = \max \left( 0, -\frac{1}{W_{\bar{c}\bar{c}}} (\mathbf{v}_{\bar{c},\text{free}}^{n+1} + \varrho \mathbf{v}_{\bar{c}}^n) \right). \quad (33) $$

In the case of a distributed obstacle, a LCP solver has to be employed. In this article, we employ the numerical solvers provided by the Siconos software (Acary et al., 2016) which is dedicated to the modelling and simulation of nonsmooth dynamical systems. In Siconos, several numerical algorithms are implemented to solve LCPs. Depending on the number of constraints and the required accuracy, pivoting or projected successive overrelaxation (PSOR) techniques may be used (Cottle et al., 1992). Simulations in the present paper are based on the use of a standard pivoting technique for solving LCP, known as Lemke's method (Cottle et al., 1992) to get a high precision solution at the machine accuracy.

# IV. NUMERICAL RESULTS

## A. One contact point: the case of tanpura

In this section, the case of a two point bridge is considered, which consists in placing a point obstacle near a boundary. The string parameters are selected as in (Issanchou et al., 2017) for the purpose of comparison. Geometric and material relevant properties are recalled in Table I. The point obstacle is located at 6 mm from the boundary $x = 0$, in agreement with the range of positions
mentioned in (Valette and Cuesta, 1993) to mimick the tanpura bridge with a two point approximation.

<table><caption>TABLE I. Electric guitar string properties.</caption><thead><tr><th>L (m)</th><th>d (mm)</th><th>T (N)</th><th>μ (kg.m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1.002</td><td>0.43</td><td>180.5</td><td>1.17 × 10<sup>-3</sup></td></tr></tbody></table>

The initial position of the string is a smoothed centered triangle with a maximal amplitude $u_{0,max} = 1.8$ mm and no velocity. The nonsmooth numerical procedure is compared to the power-law approach described in (Issanchou et al., 2017), where the contact force is modeled as $f = K[\eta]_+^\alpha$, where $[\eta]_+ = \frac{1}{2}[\eta + |\eta|]$, with $K = 10^{13}$ and $\alpha = 1.5$ for simulations. In order to be in line with the assumptions retained in (Issanchou et al., 2017), we select an equal number of modes $N_m$ and grid points $N$, specifically $N_m = N - 1$, and set $N = 1002$ for the simulations. For completeness, the outcomes of the nonsmooth method are compared for both $\varrho = 0$ and $\varrho = 1$.

A convergence study has been led using the stringent criterion introduced in (Issanchou et al., 2017) to control the long-term behaviour. The criterion inspects relative errors of the numerical results on a 3 seconds simulation and imposes the error to be less than $10^{-1}$. Following that, the sampling frequency is selected as $F_s = 2$ MHz for both compliant and nonsmooth cases, since it has been found that the convergence is obtained at that sampling rate for the two methods.

The linear modal characteristics of the string are also selected as in (Issanchou et al., 2017), eigenfrequencies and damping ratios resulting from an experimental identification on a real guitar string. The identification has been realised up to 7200 Hz, then models are employed to input linear characteristics into the simulation up to the desired maximal frequency. Models used for linear characteristics are recalled in Appendix, and selected values resulting from model fitting shown in (Issanchou et al., 2017) are given in Table II.

<table><caption>TABLE II. Model parameters for the electric guitar string.</caption><thead><tr><th>B</th><th>δ<sub>ve</sub></th><th>Q<sub>te</sub><sup>-1</sup></th></tr></thead><tbody><tr><td>1.78 × 10<sup>-5</sup></td><td>4.5 × 10<sup>-3</sup></td><td>2.03 × 10<sup>-4</sup></td></tr></tbody></table>

Fig. 2 shows the string displacement at a point located at 1 cm from $x = L$. Associated sound files corresponding to the displacement resampled at 44.1 kHz are available as supplementary materials¹. Its appears that selecting $\varrho = 0$ or $\varrho = 1$ has no significant macroscopic effect on the simulation result. Both values give a temporal displacement extremely similar to that obtained with the penalty approach. In all cases, a crenel shape appears which evolves in time due to effects of dispersion and damping. As analysed for the power-law approach in (Issanchou et al., 2017) where simulation outcomes are

compared to experimental signals, we can thus conclude that the nonsmooth method also retrieves the essential features of the string vibration against a point obstacle, with the same global accuracy as that obtained with the penalty approach.

Spectrograms of numerical signals with the penalty approach and the nonsmooth method with $\varrho = 1$ are shown in Fig. 3. No significant difference is observed when $\varrho = 0$, therefore we do not show the spectrogram in this case. In both cases, there is no missing mode despite the centered initial condition, due to energy transfers induced by collisions. Moreover, both methods accurately recover the descending formant observed experimentally. Such formants constitute a distinctive feature of musical instruments such as tanpuras.

We now explore the local behaviour of the simulated string displacement at the obstacle position $u_c$ in Fig. 4, with the nonsmooth method, $\varrho = 0$ and $\varrho = 1$, as well as with the penalty approach. Zooms focus on the first two contact periods. Despite similar global waveshapes, the string has local specific and distinguishable behaviours, around the obstacle, when colliding with it. When $\varrho = 0$, as expected since no back velocity is enforced, the string sticks to the obstacle for a short period of time when colliding, and then leaves it when the wave propagates back (i.e. around $3.8 \times 10^{-3}$ s for the first contact time, then around $8.8 \times 10^{-3}$ s for the second contact). When $\varrho = 1$, the string bounces off the obstacle with a magnitude of about $10^{-7}$ m and the signal takes the form of one time step oscillations until the string leaves the obstacle for a longer time around $3.8 \times 10^{-3}$ s. With the penalty approach, smoother oscillations appear with similar amplitudes. Interestingly, the envelope of oscillations of the nonsmooth signal closely fits oscillations of the compliant signal. The detailed behaviour of the string, computed with the nonsmooth method and different values of $F_s$, is presented in Fig. 5. Increasing the sampling frequency does not significantly change the global string behaviour, however the local behaviour at the obstacle position is affected. Note also that a very small penetration occurs with the nonsmooth method, a numerical artifact linked to the time discretisation. This numerical penetration is reduced and small bounces can be observed when $\varrho = 1$, after the string reaches the obstacle. When $\varrho = 0$, the string sticks to the obstacle at the height given by the first discrete point where contact condition is activated, defining the depth of this numerical penetration.

Finally, nonsmooth and power-law approaches give very similar results in the case of a point obstacle, and no significant difference is observed between simulations with $\varrho = 0$ and $\varrho = 1$ in the nonsmooth case, except when focusing on the obstacle position. This can be related to the analytic solutions for a string with a single contact point where a sticking time is obtained as long as the wave has not come back, see e.g. (Cabannes, 1984; 1987) and (Doyen et al., 2011; Yoong et al., 2018) for a similar case with bars.

We now discuss time costs of the nonsmooth ($\varrho = 1$) and power-law ($\alpha = 1.5$, $K = 10^{13}$) approaches. Com-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Displacement of the string vibrating against a point obstacle near a boundary, taken at 1 cm from the extremity $x = L$. Comparison between the penalty approach, $F_s = 2$ MHz (green line) and the nonsmooth method, $F_s = 2$ MHz, $\rho = 0$ (black dashed line) and $\rho = 1$ (red line).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Spectrograms of the string displacement in the case of a point obstacle near a boundary (dB, with a 70 dB dynamic) (a) power-law method, $F_s = 2$ MHz (b) nonsmooth method, $F_s = 2$ MHz, $\rho = 1$.</figcaption></figure>

putation costs, excluding initialisation, saving data and energy computation, are presented in Table III. For the compliant method, they are presented either with a spatial grid which is fully included in the matrix $S$, meaning that its dimensions are $(N - 1)^2$, or with a matrix $S$ which only includes the obstacle position, meaning that its dimensions are $1 \times N - 1$. The former case strictly corresponds to the scheme presented in (Issan-chou et al., 2017), while the latter case (modified compliant method) contains an improvement where the number of grid points $N$ is not anymore related to the number of modes $N_m$. This further step is described in (van Wal-stijn and Bridges, 2016), and is also implemented for the nonsmooth method. All computations presented in Table III were led with MATLAB on a single CPU with a clock at 2.4 GHz.

One can observe that taking into account only the contact point in the modified compliant method allows one to obtain a significant gain in computational times since a factor two is present as compared to the power-law method with a full $S$ matrix. A more significant gain can be obtained due to the nonsmooth method: computation costs are 12 to 37 times smaller. This is due in particular to the reduced number of matrix vector products as well as a reduced time for solving the LCP (straightforward in the point obstacle case) compared to Newton-Raphson iterations. Note also that computation times has been found to be once again divided by a factor two by running the same code with the SICONOS software in a optimized C/C++ implementation (Acary et al., 2016).

In this section, we discussed the behaviour of a string vibrating against a point obstacle, mimicking the case of
<figure><img src="image_5.png" /><figcaption>FIG. 4. Displacement of the string vibrating at the obstacle position. Comparison between the numerical simulation with the penalty approach, $F_s = 2$ MHz (green line) and the nonsmooth method with $\rho = 0$ (black dashed line, black line in zooms) and $\rho = 1$ (red line), $F_s = 2$ MHz.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Displacement of the string at the obstacle position. Comparison between the numerical nonsmooth simulation with $F_s = 2$ MHz ($\rho = 0$: black dashed line, $\rho = 1$: red line) and $F_s = 8$ MHz ($\rho = 0$: blue dashed line, $\rho = 1$: magenta line): (a) first overall string/obstacle contact (b) second overall string/obstacle contact.</figcaption></figure>

a tanpura. We observed that the nonsmooth simulation in the case of a two point bridge was very close to the simulation led with the power-law method, both being consistent with experimental data presented in (Issanchou et al., 2017). Moreover, it seems that the global behaviour of the string does not rely on the coefficient of restitution, even though the local string behaviour at the contact point does. In the following, we investigate a case implying multiple contact points.

## B. Multiple contact points: application to the electric bass

In this section, the complexity is increased by considering a string vibrating against the fretboard of an electric bass with up to 20 contact points. Electric basses are known for their modern playing techniques such as pop and slap where numerous contacts are intentionally provoked in the transient attack, giving a peculiar bright and percussive sound (Bacon, 2013; Issanchou et al., 2018).
<table><caption>TABLE III. Computation times with the compliant and nonsmooth methods (implemented in Matlab), in seconds, for the simulation of a one second signal in the case of a point obstacle.</caption><tbody><tr><th colspan="5">Compliant approach, full spatial grid included in S</th></tr><tr><th><i>F</i><sub>s</sub> (kHz)</th><td>44.1</td><td>88.2</td><td>176</td><td>1000</td></tr><tr><th>Newton-Raphson</th><td>30</td><td>55</td><td>105</td><td>362</td></tr><tr><th>Other operations</th><td>175</td><td>352</td><td>710</td><td>3964</td></tr><tr><th>Total time</th><td>205</td><td>407</td><td>815</td><td>4326</td></tr><tr><th colspan="5">Compliant approach, S at the obstacle position only</th></tr><tr><th><i>F</i><sub>s</sub> (kHz)</th><td>44.1</td><td>88.2</td><td>176</td><td>1000</td></tr><tr><th>Newton-Raphson</th><td>7.7</td><td>12</td><td>29</td><td>140</td></tr><tr><th>Other operations</th><td>100</td><td>196</td><td>343</td><td>1592</td></tr><tr><th>Total time</th><td>108</td><td>208</td><td>372</td><td>1732</td></tr><tr><th colspan="5">Nonsmooth method</th></tr><tr><th><i>F</i><sub>s</sub> (kHz)</th><td>44.1</td><td>88.2</td><td>176</td><td>1000</td></tr><tr><th>LCP</th><td>1.4</td><td>5.3</td><td>3.9</td><td>21</td></tr><tr><th>Computation in (27)</th><td>1.4</td><td>2.6</td><td>5.2</td><td>38</td></tr><tr><th>Other operations</th><td>4.3</td><td>8.0</td><td>13</td><td>83</td></tr><tr><th>Total time</th><td>7.1</td><td>16</td><td>22</td><td>142</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIG. 6. A string vibrating against a bass guitar fretboard represented by the function *g*.</figcaption></figure>

Fig. 6 shows the guitar neck profile with its 20 frets together with the initial condition imposed to the string. A G-string of an electric bass is considered for simulations, its properties are given in Table IV. The string is plucked at 64 cm from the nut with a maximal initial amplitude $u_{0,max} = 3.6$ mm. As in the previous section and according to the study led in (Issanchou et al., 2018), measured values of the linear characteristics are selected to fit with the model up to about 3400 Hz, then models are employed with parameters reported in Table V. In order to perform comparisons with the power-law approach and with experimental results obtained in (Issanchou et al., 2018), 863 modes are retained for simulations. The convergence of data is determined according to the criterion presented in (Issanchou et al., 2017; 2018) applied over a period of time recovering all string/neck collisions. This implies large values of $F_s$, nevertheless smaller values may be selected for other applications. For instance, a convergence study on the spectral content, associated to a perception evaluation of signals with different $F_s$,

may lead to smaller values of $F_s$ with a satisfactory sound rendition. This is however not the purpose of the present paper.

<table><caption>TABLE IV. Electric bass string properties.</caption><thead><tr><th><i>L</i> (m)</th><th><i>d</i> (mm)</th><th><i>d</i><sub>core</sub> (mm)</th><th><i>T</i> (N)</th><th>μ (kg.m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.863</td><td>1.14</td><td>0.43</td><td>191.6</td><td>6.69 × 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>TABLE V. Models parameters for the electric bass string installed on an electric bass.</caption><thead><tr><th><i>B</i></th><th>δ<sub>ve</sub></th><th><i>Q</i><sub>te</sub><sup>-1</sup></th></tr></thead><tbody><tr><td>3.5 × 10<sup>-5</sup></td><td>0.01</td><td>6 × 10<sup>-6</sup></td></tr></tbody></table>

Temporal results are shown in Fig. 7, with a comparison of the nonsmooth approach with the experimental result already shown in (Issanchou et al., 2018), and then with the compliant method. A spectral comparison is provided due to the spectrograms of the displacements shown in Fig. 8. Related sound files corresponding to the displacement resampled at 44.1 kHz are available as supplementary materials². The presence of multiple contacts makes the dynamics more complex but leads once again to similar results for the two limit values of the restitution coefficient, namely $\varrho = 0$ and $\varrho = 1$. Simulation results show that the temporal signals are almost coincident. Consequently only the case $\varrho = 1$ is presented. It thus appears that the dynamics of collisions, from a macroscopic point of view, is still dominated by waves going back and forth in inter-fret intervals.

Comparing outcomes of the displacement obtained with the nonsmooth method with experimental results shows an excellent agreement, both in the global shape of the time series and in the details shown in zooms in Fig. 7. The occurrence of a substantial high frequency content from the very first collisions is clearly ascertained and retrieved, and the long term behaviour remains similar. A few discrepancies appear, however as explained in (Issanchou et al., 2018) they are probably mainly due to (i) uncertainties related to the experimental measurement of the neck profile and, to a lesser extent, to the measured properties of the string; (ii) the assumption that the fretboard is a rigid obstacle. Comparing the two numerical methods shows that they behave very similarly, with only slight discrepancies.

Focusing on spectrograms, a similar structure between experimental and numerical data can be observed. In particular, energy transfers appear in the transient attack between 0 and 0.09 s, this time window corresponding to the occurrence of contacts between the string and the fretboard. Once again, this collision period of time, together with specific reinforced spectral zones (e.g. around 4500 Hz), is particularly well retrieved by numerical simulations. From a perceptive point of view, associ-
<figure><img src="image_8.png" /><figcaption>FIG. 7. Displacement of the string vibrating against the neck of a bass guitar, taken at 9 mm from the extremity $x = L$. Comparison between the experimental signal (blue line) and the nonsmooth method, $\rho = 1$, $F_s = 4$ MHz (red line), and between the nonsmooth method, $\rho = 1$, $F_s = 4$ MHz (red line) and the compliant method, $F_s = 8$ MHz (green line).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Spectrograms of the string displacement on an electric bass (dB, with a 70 dB dynamic) (a) experimental signal (b) compliant method, $F_s = 8$ MHz (c) nonsmooth method, $F_s = 4$ MHz, $\rho = 1$.</figcaption></figure>

ated sounds, given as supplementary materials, are very close to the ear.

Finally, the global dynamics of the string colliding with frets is well described by both compliant and nonsmooth approaches. In the studied configuration, both numerical methods produce very similar results at the measurement point considered, selected from experiments. Furthermore, the value of the restitution

coefficient $\rho$ does not substantially affect simulations, meaning that, in the present case, the nonsmooth approach would be free of parameters to be adjusted.

The simulation of 0.1 s of signal (time interval recovering all contacts) at $F_s = 1$ MHz lasts 5 minutes with Matlab while it lasts only 30 s with Siconos. Numerical costs are thus found to be divided by a factor 10: once
again, an important gain in computation can be achieved with the nonsmooth method on Siconos.

## V. CONCLUSION

A nonsmooth approach to the numerical simulation of musical strings colliding with an obstacle has been presented. The main feature consists in combining an exact scheme for the vibratory (linear) part and a Moreau-Jean scheme for the contact force, embedded in a modal description of the dynamics. This results in a computationally efficient scheme preventing numerical dispersion during free flight phases. The method has been tested on two cases involving either one contact point (a two point bridge mimicking a tanpura) or multiple contacts (the case of an electric bass). In both cases, it has been found that the nonsmooth approach is very accurate and compares well with experiments and simulations led with a compliant method. A significant gain in computational times is obtained with the proposed nonsmooth method. Interestingly, for the two test cases, simulation results are relatively insensitive to the choice of the value of the restitution coefficient, so that the model could finally be considered as free of parameters to describe collisions. This particular behaviour should however be limited to a few number of pointwise contacts where, from a macroscopic point of view, the behaviour during collisions is governed by waves going back and forth in the string. In contrast, the coefficient of restitution should have an effect when considering a large number of contact points or continuous contact regions (e.g. a fretless bass), where sticking may occur on a whole interval. Note also that the nonsmooth method has been tested in specific examples where a hard collision is at hand and no specific physical detail on the contact is needed to have efficient simulations accounting for the rich contact dynamics of the string. For other cases encountered in musical acoustics where the contact is softer, e.g. mallet/membrane or finger/string interactions, the nonsmooth approach may also be used, either with a refined continuous description of the deformable bodies, involving an increase of the computational burden to take into account the whole dynamics; or with a localised behaviour law to account for the felt as done for example in (Thorin et al., 2017). These points are thus left for further studies where a more complete comparison of nonsmooth and compliant methods could be given.

## APPENDIX: LINEAR CHARACTERISTICS

A complete model for the computation of eigenfrequencies is given in (Issanchou et al., 2018) which includes the tension and the stiffness of the string, together with the potential mobility at one extremity, as observed e.g. for string instruments such as bass and guitar. The eigenfrequencies are then described by:

$$ \nu_j = j \frac{c}{2L} \left( 1 + \frac{B j^2}{2} + \frac{\mu c}{j \pi} \mathrm{Im}(Y_{nut}(\omega_{0,j})) \right), \quad (A.1) $$

where $c = \sqrt{\frac{T}{\mu}}$ is the wave velocity of the ideal string, $B = \frac{\pi^2 EI}{TL^2}$ is the inharmonicity coefficient and $Y_{nut}$ is the mobility at the nut, evaluated at $\omega_{0,j} = j \frac{\pi c}{L}$. In this paper, the mobility is taken as zero in the case of a two point bridge while it takes measured values (see (Issanchou et al., 2018)) in the case of an electric bass.

Damping parameters are modeled through the quality factor $Q_j = \pi \nu_j / \sigma_j$, such that:

$$ Q_j^{-1} = Q_{j,air}^{-1} + Q_{j,ve}^{-1} + Q_{te}^{-1} + \frac{\mu c^2}{\pi L \nu_j} \mathrm{Re}(Y_{nut}(\omega_j)), \quad (A.2) $$

where subscripts air, ve and te respectively refer to losses due to air friction, viscoelastic and thermoelastic effects. Their detailed expressions are given by:

$$ Q_{j,air}^{-1} = \frac{jc}{2L\nu_j} \frac{R}{2\pi\mu\nu_j}, \quad (A.3) $$

$$ Q_{j,ve}^{-1} = \frac{4\pi^2\mu E_{core} I_{core} \delta_{ve}}{T^2} \frac{\nu_{0,j}^3}{\nu_j}, \quad (A.4) $$

where $R = 2\pi\eta_{air} + 2\pi d_{eq}\sqrt{\pi\eta_{air}\rho_{air}\nu_j}$, with $\eta_{air}$ and $\rho_{air}$ the dynamic viscosity coefficient and the air density respectively. Usual values $\eta_{air} = 1.8 \times 10^{-5} \text{ kg.m}^{-1}\text{s}^{-1}$ and $\rho_{air} = 1.2 \text{ kg.m}^{-3}$ are selected. The employed equivalent diameter is expressed as:

$$ d_{eq} = \frac{\pi}{2} \left( \left( 1 + \frac{2}{\pi} \right) 2r_{winding} + 2r_{core} \right), \quad (A.5) $$

where $r_{winding}$ and $r_{core}$ are the radii of the string winding and core respectively.

Viscoelastic effects are assumed to be located in the string core so that the Young's modulus $E_{core}$ and the moment of inertia $I = \pi r_{core}^4 / 4$ are related to the core only. The notation $\nu_{0,j}$ refers to the value $\frac{jc}{2L}$ and $Q_{te}^{-1}$ is taken as a constant.

## References

<footnote n="1">¹See Supplementary materials at [URL will be inserted by AIP].</footnote>

<footnote n="2">²See Supplementary materials at [URL will be inserted by AIP].</footnote>

<footnote>Acary, V. (2016). “Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact,” ZAMM 96(5), 585–603.</footnote>

<footnote>Acary, V., Brémond, M., Huber, O., and Pérignon, F. (2016). “Siconos. a software for modeling and simulation of nonsmooth dynamical systems” http://siconos.gforge.inria.fr, http://github.com/siconos/siconos.</footnote>

<footnote>Acary, V., and Brogliato, B. (2008). Numerical methods for nonsmooth dynamical systems. Applications in mechanics and electronics, 35 of Lecture Notes in Applied and Computational Mechanics (Springer Verlag).</footnote>

<footnote>Acary, V., de Jong, H., and Brogliato, B. (2014). “Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems,” Physica D: Nonlinear Phenomena 269(Supplement C), 103–119.</footnote>
Ahn, J. (2007). "A vibrating string with dynamic frictionless impact," Applied Numerical Mathematics 57, 861-884.

Amerio, L. (1978). "Continuous solutions of the problem of a string vibrating against an obstacle," Rendiconti del Seminario Matematico della Università di Padova 59, 67-96.

Bacon, T. (2013). "Electric bass guitar," Grove Music Online. Oxford Music Online. Oxford University Press., https://doi.org/10.1093/gmo/9781561592630.article.45729.

Bertails-Descoubes, F., Cadoux, F., Davier, G., and Acary, V. (2011). "A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies," ACM Transactions on Graphics 30(1), 6:1-6:14.

Bilbao, S. (2009). *Numerical sound synthesis: finite difference schemes and simulation in musical acoustics* (Wiley, Chichester).

Bilbao, S. (2012). "Time domain simulation of the snare drum," Journal of the Acoustical Society of America 131, 914-925.

Bilbao, S., and Torin, A. (2015). "Numerical modeling and sound synthesis for articulated string/fretboard interaction," J. Audio Eng. Soc 63, 336-347.

Bilbao, S., Torin, A., and Chatzioioannou, V. (2015). "Numerical modeling of collisions in musical instruments," Acta Acustica united with Acustica 101, 155-173.

Boutillon, X. (1988). "Model for piano hammers: Experimental determination and digital simulation," Journal of the Acoustical Society of America 83, 746-754.

Cabannes, H. (1984). "Periodic motions of a string vibrating against a fixed point-mass obstacle: II," Mathematical Methods in the Applied Sciences 6, 55-67.

Cabannes, H. (1987). "Motion of a string vibrating against a rigid fixed obstacle," Mathematical Modelling 8, 364-365.

Chatzioioannou, V., and van Walstijn, M. (2015). "Energy conserving schemes for the simulation of musical instrument contact dynamics," Journal of Sound and Vibration 339, 262-279.

Chen, Q.-Z., Acary, V., Virlez, G., and Brüls, O. (2013). "A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints," International Journal for Numerical Methods in Engineering 96(8), 487-511.

Citrini, C. (1991). "Vibrating strings with obstacles: the analytic study," Advances in kinetic theory and continuum mechanics 249-259.

Cottle, R. W., Pang, J., and Stone, R. E. (1992). *The Linear Complementarity Problem* (Academic Press, Inc., Boston, MA).

Desvages, C., and Bilbao, S. (2015). "Two-polarisation finite difference model of bowed strings with nonlinear contact and friction forces," in *International Conference on Digital Audio Effects, DAFx-15*, Trondheim.

Doyen, D., Ern, A., and Piperno, S. (2011). "Time-integration schemes for the finite element dynamic Signorini problem," SIAM Journal on Scientific Computing 33(1), 223-249.

Ducceschi, M., Bilbao, S., and Desvages, C. (2016). "Modelling collisions of nonlinear strings against rigid barriers: conservative finite difference schemes with application to sound synthesis," in *International Congress on Acoustics ICA, Buenos Aires*.

Fletcher, N. H., and Rossing, T. D. (1998). *The physics of musical instruments, second edition* (Springer New York).

Goldsmith, W. (2001). *Impact* (Dover publications, Mineola).

Inácio, O., Henrique, L., and Antunes, J. (2006). "The dynamics of tibetan singing bowls," Acta Acustica united with Acustica 92(4), 637-653.

Issanchou, C., Bilbao, S., Le Carrou, J.-L., Touzé, C., and Doaré, O. (2017). "A modal-based approach to the nonlinear vibration of strings against a unilateral obstacle: simulations and experiments in the pointwise case," Journal of Sound and Vibration 393, 229-251.

Issanchou, C., Le Carrou, J.-L., Touzé, C., Fabre, B., and Doaré, O. (2018). "String/frets contacts in the electric bass sound: simulations and experiments," Applied Acoustics 129, 217-228.

Janin, O., and Lamarque, C.-H. (2001). "Comparison of several numerical methods for mechanical systems with impacts," International Journal for numerical methods in engineering 51, 1101-1132.

Jean, M. (1995). "Frictional contact in collections of rigid or deformable bodies: numerical simulation of geomaterial motions," in *Mechanics of Geomaterial Interfaces*, edited by A. Selvadurai and M. Boulon, **42** of *Studies in Applied Mechanics* (Elsevier), pp. 463–486.

Jean, M. (1999). "The non-smooth contact dynamics method," Computer Methods in Applied Mechanics and Engineering 177(3), 235-257.

Jean, M., and Moreau, J.-J. (1987). *Dynamics in the presence of unilateral contacts and dry friction: a numerical approach*, 151-196 (Springer Vienna, Vienna).

Mandal, A. K., and Wahi, P. (2015). "Natural frequencies, modes and modal interactions for strings vibrating against an obstacle: relevance to Sitar and Veena," Journal of Sound and Vibration 338, 42-59.

Meingast, M. B., Legrand, M., and Pierre, C. (2014). "A linear complementarity problem formulation for periodic solutions to unilateral contact problems," International Journal of Non-Linear Mechanics 66, 18-27.

Moreau, J. (1988). "Bounded variation in time," in *Topics in Non-smooth Mechanics*, edited by J. Moreau, P. Panagiotopoulos, and G. Strang, Birkhäuser, Basel, pp. 1-74.

Moreau, J. (1999). "Numerical aspects of the sweeping process," Computer Methods in Applied Mechanics and Engineering 177, 329-349 special issue on computational modeling of contact and friction, J.A.C. Martins and A. Klarbring, editors.

Nguyen, N. S., and Brogliato, B. (2014). *Multiple Impacts in Dissipative Granular Chains*, **72** of *Lecture Notes in Applied and Computational Mechanics* (Springer Verlag).

Paoli, L., and Schatzman, M. (2002). "A numerical scheme for impact problems i: the one-dimensional case," SIAM Journal on Numerical Analysis 40(2), 702-733.

Paté, A., Le Carrou, J.-L., and Fabre, B. (2014). "Predicting the decay time of solid body electric guitar tones," Journal of the Acoustical Society of America 135(5), 3045-3055.

Raman, C. (1921). "On some indian stringed instruments," Indian Association for the Cultivation of Science (7), 29-33.

Renouf, M., Dubois, F., and Alart, P. (2004). "A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media," Journal of Computational and Applied Mathematics 168(1), 375-382.

Rhaouti, L., Chaigne, A., and Joly, P. (1999). "Time-domain modeling and numerical simulation of a kettledrum," Journal of the Acoustical Society of America 135(6), 3545-3562.

Schatzman, M. (1980). "A hyperbolic problem of second order with unilateral constraints: the vibrating string with a concave obstacle," Journal of Mathematical Analysis and Applications 73, 138-191.

Siddiq, S. (2012). "A physical model of the nonlinear sitar string," Archives of acoustics 37(1), 73-79.

Signorini, A. (1933). "Sopra alcune questioni di statica dei sistemi continui," *Annali della Scuola Normale Superiore di Pisa* 2, 231-251.

Studer, C. (2009). *Numerics of unilateral contacts and friction*, **47** of *Lecture Notes in Applied and Computational Mechanics* (Springer Verlag, Berlin Heidelberg).

Thorin, A., Boutillon, X., Lozada, J. and Merlhiot, X. (2017). "Non-smooth dynamics for an efficient simulation of the grand piano action," *Meccanica* 52(11-12), 2837-2854.

Trautmann, L., and Rabenstein, R. (2004). "Multirate simulations of string vibrations including nonlinear fret-string interactions using the functional transformation method," EURASIP Journal on Applied Signal Processing 7, 949-963.

Valette, C., and Cuesta, C. (1993). *Mécanique de la corde vibrante* (Hermès, Paris).

Valette, C., Cuesta, C., Besnainou, C., and Castellengo, M. (1991). "The tempura bridge as a recursive wave generator," Acustica 74, 201-208.

van Walstijn, M., and Bridges, J. (2016). "Simulation of distributed contact in string instruments: a modal expansion approach," in *European Signal Processing Conference (EUSIPCO)*, pp. 1023-1027.

van Walstijn, M., Bridges, J., and Mehes, S. (2016). "A real-time synthesis oriented tanpura model," in *International Conference on Digital Audio Effects, DAFx-16*, Brno.

Vyasarayani, C. P., Birkett, S., and McPhee, J. (2009). "Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar," Journal of the