# Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar

Chandrika P Vyasarayani, Stephen Birkett, John Mcphee

► To cite this version:

Chandrika P Vyasarayani, Stephen Birkett, John Mcphee. Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar. Journal of the Acoustical Society of America, 2009, 125 (6), pp.3673-3682. 10.1121/1.3123403. hal-01525495

HAL Id: hal-01525495

https://hal.science/hal-01525495

Submitted on 21 May 2017

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
# Modeling the dynamics of a vibrating string with a finite distributed unilateral constraint: Application to the sitar

Chandrika P. Vyasarayani, Stephen Birkett, and John McPhee

Department of Systems Design Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

The free vibration response of an ideal string impacting a distributed parabolic obstacle located at its boundary has been analyzed, the goal being to understand and simulate a sitar string. The portion of the string in contact with the obstacle is governed by a different partial differential equation (PDE) from the free portion represented by the classical string equation. These two PDEs and corresponding boundary conditions, along with the transversality condition that governs the dynamics of the moving boundary, are obtained using Hamilton's principle. A Galerkin approximation is used to convert them into a system of nonlinear ordinary differential equations, with lower mode-shapes parametrized with respect to the location of the moving boundary as basis functions. This system is solved numerically and the behavior of the string studied from simulations. The advantages and disadvantages of the proposed method are discussed in comparison to the penalty approach for simulating wrapping contacts. Simulations with bridge-string parameters consistent with the configuration of a real sitar show that any degree of obstacle wrapping may occur during normal playing. Finally, the model is used to investigate the mechanism behind the generation of the buzzing tone in a sitar.

## I. INTRODUCTION

The sitar (Fig. 1) is an Indian musical instrument¹ with plucked strings that can interact with a shallow curved ledge situated underneath the vibrating length of the string at one of its boundaries. The bridge, which includes the ledge, as well as grooves to constrain the strings in their correctly spaced lateral positions for playing, is typically carved from a piece of bone and rests on small wooden feet in contact with the soundboard of the instrument (Fig. 2). The tone of the sitar, and other instruments of Indian origin with a similar bridge design, such as the veena and tambura,¹ is markedly different from that of Western plucked stringed instruments such as a guitar. The interaction of the string with the obstacle generates high frequency components and creates a characteristic buzzing sound.

This mechanism is also not unknown in Western instruments. The Medieval and Renaissance bray harp, for example, has small bray-pins which provide a metal surface for the vibrating string to impact close to one end, increasing the upper partial content in the tone and providing a means for the harp to be audible in larger spaces and in ensemble with other instruments.² The arpichordium stop common on some virginals (a plucked string keyboard instrument) imitates the bray harp by soft metal (lead or brass) pins that can be bent so as to lie close to the vibrating strings about 15 mm from the termination point.³

The phenomenon common to all of these instruments is caused by the presence of a physical obstacle which alters the behavior of a vibrating string by interacting with it close

to one of the termination points. Raman⁴ gave a detailed description of bridge geometry for the sitar, veena, and tambura. He noted that even though a veena string may be plucked at a node the corresponding vibration mode appears in the response. Raman concluded that this phenomenon is a consequence of the interaction of the string with the bridge. Various approaches for modeling the interaction between a vibrating string and an obstacle have been presented in the literature. Amerio and Prouse,⁵ Schatzman,⁶ Burridge et al.,⁷ and Cabannes⁸⁻¹⁰ used the method of characteristics and energy conservation for simulating the impact between a string and a rigid obstacle. Cabannes¹⁰ noted that modeling the case of a string not initially at rest is an open problem. Ahn¹¹ used a finite element approach and Newton’s kinetic coefficient of restitution to simulate the string and its impact with the bridge. Han and Grosenbaugh¹² and Taguti¹³ simulated the impact using a penalty approach combined with finite difference discretization for string motion, while Vyasarayani et al.¹⁴ used the penalty approach with a modal representation of string motion. A different direction was taken by Krishnaswamy and Smith¹⁵ who considered the curved sitar-bridge as a point obstacle and applied digital wave-guides and finite difference methods to obtain the solution for a rigid impact. Velette¹⁶ analyzed the tambura string interaction with a distributed obstacle modeled as a unilateral point constraint considering completely plastic impact.¹⁷

When the obstacle is located near the termination of the vibrating continuum, it is possible to model the dynamics using a moving boundary problem. Fung and Chen¹⁸ proposed this method and studied perfect wrapping of a flexible cantilever beam on a circular rigid foundation as a moving boundary problem. The sitar-bridge-string problem is closely related to this work, but differs from it in several respects: (i)
<figure><img src="image_2.png" /><figcaption>FIG. 1. Sitar, a stringed instrument of Indian origin.</figcaption></figure>

the vibrating continuum is a string instead of a beam; (ii) the string can have a non-point contact length at static equilibrium; (iii) the contacting boundary is of finite length, thus limiting the maximum amount of wrapping around the obstacle; and (iv) the obstacle geometry is closer to a parabola than a circle.

The present paper describes a new modal formulation for solving the string contact problem for a finitely terminated parabolic obstacle, obtaining the equation of motion of the string during wrapping motion from Hamilton's principle. The developed model is used to study the general behavior of the string motion. Advantages and disadvantages of the proposed method for simulating wrapping contacts are discussed in comparison to the penalty approach. The parameter space for which the simulation model approximates the sitar-bridge-string interaction is analyzed and the mechanism behind the generation of the buzzing tone in a sitar is investigated.

## II. MATHEMATICAL MODELING

A schematic representation of the bridge-string geometry is shown in Fig. 3. The bridge is a finite obstacle defined by a parabolic surface for $X$ values between $\Gamma_1$ and $\Gamma_2$. The string has fixed termination points on the $X$-axis at $X=L$, and on the parabola at $X=\Gamma_1$.

String motion can be divided into three distinct phases: phase-I motion occurs when there is no contact with the obstacle; partial wrapping on the obstacle is called phase-II motion; a completely wrapped string is considered to be in phase-III motion. In this section we assume, for simplicity and without loss of generality, that $\Gamma_1=0$ and $\Gamma_2=B$. The equation of motion governing the dynamics of the string during each of the three phases is derived, as well as the corresponding switching conditions as the string passes between the phases. This approach is more general than that of Burridge et al.⁷ who only considered phase-II motion because the parabola extends indefinitely below the string.

<figure><img src="image_3.png" /><figcaption>FIG. 3. (Color online) Bridge-string geometry and different phases of string motion (exaggerated for clarity). The right string termination lies on the $X$-axis; the left termination lies on the parabolic bridge surface as shown.</figcaption></figure>

## A. Phase-I motion

During phase-I motion the string is governed by the classical string equation given as follows:

$$ \rho\Lambda\frac{\partial^2 Y_1}{\partial t^2} - T\frac{\partial^2 Y_1}{\partial X^2} = 0, \qquad (1) $$

with boundary conditions

$$ Y_1(0,t) = Y_1(L,t) = 0, \qquad (2) $$

where $Y_1$ is the transverse displacement of the string, $L$ is the length of the string, $T$ is the string tension, $\rho$ is the density, $A$ is the cross-sectional area, $X$ is the co-ordinate along the length, and $t$ is the time. The shape of the obstacle is assumed to be a parabola as in the work of Burridge et al.⁷ and the geometry can be analytically represented as

$$ Y_B(X) = A_p X(B - X). \qquad (3) $$

We substitute the following non-dimensional parameters into the equation of motion to facilitate analysis:

$$ y_1 = \frac{Y_1}{h}, \quad x = \frac{X}{L}, \quad \tau = t \sqrt{\frac{T}{\rho A L^2}}, \qquad (4) $$

where $h=A_p B^2/4$ is the height of the obstacle. The equation of motion after substituting the non-dimensional variables is as follows:

$$ \frac{\partial^2 y_1}{\partial \tau^2} - \frac{\partial^2 y_1}{\partial x^2} = 0, \qquad (5) $$

with boundary conditions

<figure><img src="image_4.png" /><figcaption>FIG. 2. Sitar-bridge. The profile can be approximated by a parabolic curve. The base of the parabola defines the $x$-axis passing through the far string termination point and parallel to the reference line shown, which runs along the neck of the instrument. See text for detailed dimensions.</figcaption></figure>
$$y_1(0, \tau) = y_1(1, \tau) = 0. \qquad (6)$$

A solution to Eq. (5) is assumed to be of the form

$$y_1(x, \tau) = \sum_{j=1}^{\infty} \phi_j(x) \eta_j(\tau). \qquad (7)$$

In Eq. (7), $\phi_j(x) = \sqrt{2} \sin(j\pi x)$ are mass-normalized modes of the string and $\eta_j(\tau)$ are modal co-ordinates. Substituting Eq. (7) into Eq. (5), multiplying by $\phi_k(x)$, integrating over the domain, and simplifying the resulting equation by using orthogonality relations result in a set of uncoupled ordinary differential equations of the following form:

$$\ddot{\eta}_j(\tau) + \omega_j^2 \eta_j(\tau) = 0, \qquad (8)$$

where $\omega_j = j\pi$ are the natural frequencies of the string. The modal initial conditions corresponding to physical initial conditions of $y_1(x, 0)$ and $\dot{y}_1(x, 0)$ are as follows:

$$\eta_j(0) = \int_0^1 y_1(x, 0) \phi_j(x) dx \quad \text{and} \\ \dot{\eta}_j(0) = \int_0^1 \dot{y}_1(x, 0) \phi_j(x) dx. \qquad (9)$$

After specifying the system parameters and initial conditions, the modal equations of motion given by Eq. (8) can be numerically integrated using any standard numerical scheme. When the string starts to contact the obstacle, the equation of motion given by Eq. (8) is no longer valid. In Sec. II B we derive the equation of motion during wrapping motion of the string around the obstacle by using Hamilton's principle with a moving boundary.

## B. Phase-II motion

The Hamiltonian with the distributed spatial constraint¹⁸ can be written as follows:

$$H = \delta \int_{t_1}^{t_2} \left[ \int_{0}^{\Gamma_{-}} (\Pi + \lambda(X)G(X))dX + \int_{\Gamma_{+}}^{L} \Pi dX \right] dt = 0, \qquad (10)$$

where $\Pi$ is the Lagrangian density function defined as

$$\Pi = \frac{1}{2} \rho A \dot{Y}_2(X, t)^2 - \frac{1}{2} T \left( \frac{\partial Y_2(X, t)}{\partial X} \right)^2, \qquad (11)$$

and $G(X)$ is the gap function defined as $G(X) = Y_2(X, t) - Y_B(X)$ and $\Gamma$ is the wrapped string length. $Y_2$ is the displacement of the string during wrapping motion and $\lambda(X)$ is the distributed constraint force. Substituting Eq. (11) into the Hamiltonian [Eq. (10)] and simplifying

$$H = \delta \int_{t_1}^{t_2} \int_{0}^{\Gamma_{-}} \left( \frac{1}{2} \rho A \dot{Y}_2(X, t)^2 \right) dX dt \\ - \delta \int_{t_1}^{t_2} \int_{0}^{\Gamma_{-}} \left( \frac{1}{2} T \left( \frac{\partial Y_2(X, t)}{\partial X} \right)^2 \right) dX dt \\ + \delta \int_{t_1}^{t_2} \int_{0}^{\Gamma_{-}} (\lambda(X)G(X)) dX dt$$

$$+ \delta \int_{t_1}^{t_2} \int_{\Gamma_{+}}^{L} \left( \frac{1}{2} \rho A \dot{Y}_2(X, t)^2 \right) dX dt \\ - \delta \int_{t_1}^{t_2} \int_{\Gamma_{+}}^{L} \left( \frac{1}{2} T \left( \frac{\partial Y_2(X, t)}{\partial X} \right)^2 \right) dX dt = 0. \qquad (12)$$

After evaluating the variations in Eq. (12),

$$H = \int_{t_1}^{t_2} \int_{0}^{\Gamma_{-}} \left( T \frac{\partial^2 Y_2(X, t)}{\partial X^2} - \rho A \frac{\partial^2 Y_2(X, t)}{\partial t^2} + \lambda(X) \right) \times \delta Y_2(X, t) dX dt + \int_{t_1}^{t_2} \int_{0}^{\Gamma_{-}} (\delta \lambda(X) G(X)) dX dt \\ + \int_{t_1}^{t_2} \int_{\Gamma_{+}}^{L} \left( T \frac{\partial^2 Y_2(X, t)}{\partial X^2} - \rho A \frac{\partial^2 Y_2(X, t)}{\partial t^2} \right) \delta Y_2(X, t) dX dt \\ + \int_{t_1}^{t_2} \left( T \frac{\partial Y_2(\Gamma_{+}, t)}{\partial X} \delta Y_2(\Gamma_{+}, t) \right) dt \\ - \int_{t_1}^{t_2} \left( T \frac{\partial Y_2(\Gamma_{-}, t)}{\partial X} \delta Y_2(\Gamma_{-}, t) \right) dt \\ + \int_{t_1}^{t_2} \left( T \frac{\partial Y_2(0, t)}{\partial X} \delta Y_2(0, t) \right) dt \\ - \int_{t_1}^{t_2} \left( T \frac{\partial Y_2(L, t)}{\partial X} \delta Y_2(L, t) \right) dt = 0. \qquad (13)$$

It should be noted that the virtual variables $\delta Y_2(\Gamma, t)$ and $\delta \Gamma$ are unspecified, but they are related due to the presence of the geometrical constraint. After some tedious algebra, the relation between the virtual variables can be shown to be

$$\delta Y_2(\Gamma, t) = 2A_p(B - 2\Gamma)\delta\Gamma. \qquad (14)$$

The reader can find a more detailed derivation on relating the virtual variables in the work of Fung and Chen,¹⁸ for a similar problem. Substituting Eq. (14) into Eq. (13), and utilizing the fact that the virtual variables are arbitrary and must vanish at the boundaries results in the following equations of motion that must be satisfied at all time:

$$T \frac{\partial^2 Y_2(X, t)}{\partial X^2} - \rho A \frac{\partial^2 Y_2(X, t)}{\partial t^2} + \lambda(X) = 0, \quad 0 < X < \Gamma_{-}, \qquad (15)$$

$$T \frac{\partial^2 Y_2(X, t)}{\partial X^2} - \rho A \frac{\partial^2 Y_2(X, t)}{\partial t^2} = 0, \quad \gamma_+ < X < L, \qquad (16)$$

with the boundary conditions

$$Y_2(0, t) = 0, \quad Y_2(\Gamma_{-}, t) = A_p \Gamma_{-} (B - \Gamma_{-}), \qquad (17)$$

$$Y_2(L, t) = 0, \quad Y_2(\Gamma_{+}, t) = A_p \Gamma_{+} (B - \Gamma_{+}), \qquad (18)$$

since $\delta\Gamma = \delta\Gamma_{-} = \delta\Gamma_{+}$. The transversality condition can be written as
$$ \frac{\partial Y_2(\Gamma(t),t)}{\partial X} = A_P(B - 2\Gamma(t)). \quad (19) $$

The transversality condition is the necessary condition that must be satisfied for the variations to vanish at the free boundary $\gamma$. The physical implication of the transversality condition is the enforcement of the string slope to be equal to that of the slope of the parabola at the point of separation ($\Gamma$). We substitute the same non-dimensional parameters given in Eq. (4) along with $\gamma=\Gamma/L$, $b=B/L$, and $\alpha = 4L^2/B^2$ into Eqs. (15)–(19). It should be noted that the solution of Eq. (15) is the geometry of the parabola, as the string in the domain $0<x<\gamma_-$ perfectly wraps around the obstacle. So we have to solve for the motion of the string in the domain $\gamma_+<x<1$; thus the equation of motion reduces to the following moving boundary problem:

$$ \frac{\partial^2 y_2(x, \tau)}{\partial x^2} - \frac{\partial^2 y_2(x, \tau)}{\partial \tau^2} = 0, \quad \gamma_+ \le x \le 1, \quad (20) $$

$$ y_2(\gamma_+, \tau) = \alpha \gamma_+(b - \gamma_+), \quad y_2(1, \tau) = 0. \quad (21) $$

We need one further equation that should be solved for obtaining the separation point (moving boundary), which comes from the transversality condition

$$ \frac{\partial y_2(\gamma, \tau)}{\partial x} = \alpha(b - 2\gamma). \quad (22) $$

To transform the non-homogenous boundary conditions given by Eq. (21) into homogenous boundary conditions, the following transformation defining $y_3$ is substituted into Eqs. (20)–(22):

$$ y_2(x, \tau) = y_3(x, \tau) + s(x, \tau), \quad (23) $$

where

$$ s(x, \tau) = \frac{\alpha \gamma_+(b - \gamma_+)}{(\gamma_+ - 1)(x - 1)}. \quad (24) $$

The transformed equations are shown below:

$$ \frac{\partial^2 y_3(x, \tau)}{\partial x^2} - \frac{\partial^2 y_3(x, \tau)}{\partial \tau^2} = - \frac{\partial^2 s(x, \tau)}{\partial \tau^2}, \quad \gamma_+ < x < 1, \quad (25) $$

with boundary conditions

$$ y_3(\gamma_+, \tau) = 0, \quad y_3(1, \tau) = 0. \quad (26) $$

The transversality condition now becomes

$$ \frac{\partial y_3(\gamma_+, \tau)}{\partial x} = \alpha(b - 2\gamma_+) - \frac{\alpha\gamma_+(b - \gamma_+)}{(\gamma_+ - 1)}. \quad (27) $$

Substituting a solution of the following form:

$$ y_3(x, \tau) = \sum_{i=1}^{N} \psi_i(x, \tau) \beta_i(\tau) \quad (28) $$

into Eq. (25), where $\psi_i(x, \tau)$ are the mass-normalized mode-shapes parametrized with respect to the moving boundary $\gamma_+(\tau)$. They are obtained by solving for the mode-shapes of Eq. (25):

$$ \psi_i(x, \tau) = \sqrt{\frac{2}{(1 - \gamma_+)}} \sin\left(j\pi \frac{(x - \gamma_+)}{(1 - \gamma_+)}\right). \quad (29) $$

The partial differential equation (25) after substituting the solution (28) and using the orthogonality property of $\psi_j$ reduces to the following coupled ordinary differential equations with time dependent coefficients:

$$ \ddot{\beta}_i(\tau) + \left[ \sum_{j=1}^{N} 2C_{ij}(\gamma_+, \dot{\gamma}_+) \right] \dot{\beta}_i(\tau) + \left[ \left( \frac{j\pi}{(1 - \gamma_+)} \right)^2 + \sum_{j=1}^{N} D_{ij}(\gamma_+, \dot{\gamma}_+, \ddot{\gamma}_+) \right] \beta_i(\tau) = -E_{ij}(\gamma_+, \dot{\gamma}_+, \ddot{\gamma}_+). \quad (30) $$

In the above equation, $C_{ij}$, $D_{ij}$, and $E_{ij}$ are defined as follows:

$$ C_{ij}(\gamma_+, \dot{\gamma}_+) = \int_{\gamma_+}^{1} \psi_i(x, \tau) \dot{\psi}_j(x, \tau) dx = J_{1ij}(\gamma) \dot{\gamma}, \quad (31) $$

$$ D_{ij}(\gamma_+, \dot{\gamma}_+, \ddot{\gamma}_+) = \int_{\gamma_+}^{1} \psi_i(x, \tau) \ddot{\psi}_j(x, \tau) dx = J_{1ij}(\gamma) \ddot{\gamma} + J_{2ij}(\gamma) \dot{\gamma}^2, \quad (32) $$

$$ E_{ij}(\gamma_+, \dot{\gamma}_+, \ddot{\gamma}_+) = \int_{\gamma_+}^{1} \ddot{s}_i(x, \tau) \psi_j(x, \tau) dx = J_{3ij}(\gamma) \dot{\gamma}^2 + J_{4ij}(\gamma) \ddot{\gamma}, \quad (33) $$

where $J_{1ij}$, $J_{2ij}$, $J_{3ij}$, and $J_{4ij}$ are functions of $\gamma$ only. The transversality condition [Eq. (27)], after substituting the solution given by Eq. (28), becomes

$$ \sum_{i=1}^{N} i\pi\sqrt{2}\beta_i(\tau) = (1 - \gamma_+)^{3/2} \left( \alpha(b - 2\gamma_+) - \frac{\alpha\gamma_+(b - \gamma_+)}{(\gamma_+ - 1)} \right). \quad (34) $$

Differentiating the above equation twice, we get

$$ \sum_{i=1}^{N} i\pi\sqrt{2}\ddot{\beta}_i(\tau) = H_1(\gamma_+)\ddot{\gamma}_+(t) + H_2(\gamma_+)\dot{\gamma}_+^2, \quad (35) $$

where $H_1$ and $H_2$ are functions of $\gamma$. Equations (30) and (35) can be solved simultaneously for $\beta_i(\tau)$ and $\gamma(\tau)$ to predict the motion of the string during the contact phase. The above method of satisfying the second derivative of the displacement constraint instead of displacement constraint directly, thus converting the constraint equation into a differential equation, is a very well known procedure in the field of multibody dynamics.¹⁹

### C. Switching conditions between phase-I and phase-II motions

Let $\tau_{c1}$ be the time at which the string in phase-I comes in contact with the obstacle. The string in phase-I motion contacts the obstacle when the slope of the string at $x=\gamma_+ = 0$ matches with the slope of the obstacle. At the event of contact, we have
$$y_1(x, \tau_{c1}) = y_2(x, \tau_{c1}), \quad (36)$$

which can be further written as follows:

$$\sum_{i=1}^{N} \eta_{i}(\tau_{c1}) \phi_{i}(x) = \sum_{i=1}^{N} \psi_{i}(x, \tau_{c1}) \beta_{i}(\tau_{c1}) + \frac{\alpha \gamma_{+} (b - \gamma_{+})}{(\gamma_{+} - 1)} (x - 1). \quad (37)$$

Since $\gamma_+=0$ at $\tau=\tau_{c1}$ and $\phi_i(0)=\psi_i(0, \tau_{c1})=0$, we have

$$\eta_i(\tau_{c1}) = \beta_i(\tau_{c1}). \quad (38)$$

Since the velocity distributions in phase-I motion and phase-II motion should also be equal at transferring time $\tau_{c1}$, we have

$$\sum_{i=1}^{N} \dot{\eta}_i(\tau_{c1}) \phi_i(x) = \sum_{i=1}^{N} (\dot{\psi}_i(x, \tau_{c1}) \beta_i(\tau_{c1}) + \psi_i(x, \tau_{c1}) \times \dot{\beta}_i(\tau_{c1})) - b\alpha(x-1)\dot{\gamma}_+. \quad (39)$$

Multiplying both sides of Eq. (39) with $\psi_j(x, \tau_{c1})$ and integrating over the domain result in the following:

$$\begin{aligned} \dot{\eta}_i(\tau_{c1}) \sum_{i=1}^{N} \int_{0}^{1} \phi_i(x) \psi_j(x, \tau_{c1}) dx \\ = \beta(\tau_{c1}) \sum_{i=1}^{N} \int_{0}^{1} \dot{\psi}_i(x, \tau_{c1}) \psi_j(x, \tau_{c1}) dx \\ + \dot{\beta}_i(\tau_{c1}) \sum_{i=1}^{N} \int_{0}^{1} \psi_i(x, \tau_{c1}) \psi_j(x, \tau_{c1}) dx \\ - b\alpha \dot{\gamma}_+ \int_{0}^{1} (x-1) \psi_j(x, \tau_{c1}) dx. \end{aligned} \quad (40)$$

After some simplifications using orthogonality relations, the above equation reduces to

$$\dot{\eta}_i(\tau_{c1}) = \dot{\beta}_i(\tau_{c1}) + \left[ \beta(\tau_{c1}) \sum_{j=1}^{N} J_{1ij}(0) + \sqrt{2} \left( \frac{b\alpha}{i\pi} \right) \right] \dot{\gamma}_+. \quad (41)$$

The above set of $N$ equations contain $N+1$ unknowns, so we need one further equation to solve for $\dot{\gamma}_+$ which can be obtained by differentiating the transversality condition given by Eq. (34) with respect to time

$$\sum_{i=1}^{N} i \pi \sqrt{2} \dot{\beta}_i(\tau) = H_1(\gamma_+) \dot{\gamma}_+ \quad (42)$$

at transferring time $\tau=\tau_{c1}$. The above equation becomes

$$\dot{\gamma} = \frac{1}{H_1(0)} \sum_{i=1}^{N} i \pi \sqrt{2} \dot{\beta}_i(\tau_{c1}). \quad (43)$$

The above algebraic equations (41) and (43) can be solved to obtain $\dot{\beta}_i(\tau_{c1})$ and $\dot{\gamma}_+$. The time of switching can be obtained through the transversality condition as shown below:

$$\sum_{i=1}^{N} i \pi \sqrt{2} \eta_{i}(\tau_{c1}) = \sum_{i=1}^{N} i \pi \sqrt{2} \beta_{i}(\tau_{c1}) = \alpha b. \quad (44)$$

Standard event detection algorithms can be used in the simulation to detect the time at which Eq. (44) holds.

## D. Phase-III motion

The phase-III motion is the same as phase-I motion except that the string is completely wrapped around the obstacle and vibrates between $x=\gamma_+=b$ and $x=1$. The dimensionless equation of motion of the string during phase-III can be written as

$$\frac{\partial^2 y_4(x, \tau)}{\partial x^2} - \frac{\partial^2 y_4(x, \tau)}{\partial \tau^2} = 0, \quad b \le x \le 1, \quad (45)$$

$$y_4(b, \tau) = 0, \quad y_4(1, \tau) = 0, \quad (46)$$

substituting a solution of the following form:

$$y_4(x, \tau) = \sum_{j=1}^{\infty} \varphi_j(x) r_j(\tau) \quad (47)$$

into Eq. (45) and performing standard modal analysis, we get the following uncoupled ordinary differential equations:

$$\ddot{r}_j(\tau) + \omega_{4j}^2 r_j(\tau) = 0, \quad (48)$$

where $\varphi_j(x) = \sqrt{2/(1-b)}\sin(j\pi(x-b)/(1-b))$ are the mass-normalized mode-shapes of the string, $r_j(\tau)$ are the generalized coordinates, and $\omega_{4j} = j\pi/(1-b)$ are the natural frequencies of the string.

## E. Switching conditions between phase-II and phase-III motions

Once the slope of the string at $x=\gamma=b$ matches with the slope of the obstacle the string enters into phase-III motion. Let $\tau_{c2}$ be the instant of switching. The displacement and velocity distributions of the string during the last instant of phase-II will be transferred to the phase-III. This can be mathematically represented as

$$\beta_i(\tau_{c2}) = r_i(\tau_{c2}) \quad (49)$$

and

$$\dot{\beta}_i(\tau_{c2}) = \dot{r}_i(\tau_{c2}). \quad (50)$$

The transfer time $\tau_{c2}$ between phase-II and phase-III can be again obtained from the transversality condition as shown below:

$$\sum_{i=1}^{N} i \pi \sqrt{2} \beta_{i}(\tau) = (1 - \gamma_{+})^{3/2} \left( \alpha(b - 2\gamma_{+}) - \frac{\alpha\gamma_{+}(b - \gamma_{+})}{(\gamma_{+} - 1)} \right). \quad (51)$$

Substituting $\tau=\tau_{c2}$ and $\gamma_+=b$, we get

$$\sum_{i=1}^{N} i \pi \sqrt{2} \beta_{i}(\tau_{c2}) = - \alpha b (1 - b)^{3/2}. \quad (52)$$
## F. Switching conditions between phase-III and phase-II motions

The switching conditions between phase-III and phase-II are similar to that of phase-I and phase-II. Let $\tau_{c3}$ be the transfer time, once the phase-III motion is initiated after time $\tau_{c2}$ the string vibrates downwards between the boundaries b and 1. Once the string starts to move upwards the slope of the string $x=\gamma=b$ matches with the slope of the obstacle and the string enters in to phase-II motion, and then it tries to unwrap itself thus again performing a phase-II motion. Now we try to relate the initial conditions between phase-III and phase-II motions. Following similar procedure as in Sec. II C, we get the following relations:

$$r_i(\tau_{c3}) = \beta_i(\tau_{c3}), \qquad (53)$$

$$\dot{r}_i(\tau_{c3}) = \dot{\beta}_i(\tau_{c1}) + \left[ \beta(\tau_{c1}) \sum_{j=1}^{N} J_{1ij}(b) + \sqrt{2} \frac{b\alpha}{i\pi} \right] \dot{\gamma}_+, \qquad (54)$$

$$\dot{\gamma}_+ = \frac{1}{H_1(b)} \sum_{i=1}^{N} i\pi\sqrt{2}\dot{\beta}_i(\tau_{c3}). \qquad (55)$$

Equation (55) can be substituted into Eq. (54) to eliminate $\dot{\gamma}_+$ and that can be solved for $\dot{\beta}_j(t)$. Once $\dot{\beta}_j(t)$ are known they can be re-substituted into Eq. (55) to get $\dot{\gamma}_+$. The transferring time $\tau=\tau_{c3}$ can again obtained from the transversality condition as shown below:

$$\sum_{i=1}^{N} i\pi\sqrt{2}\beta_i(\tau_{c3}) = -\alpha b(1-b)^{3/2}. \qquad (56)$$

Equations (53)–(55) relate the initial conditions between phase-III and phase-II.

## G. Switching conditions between phase-II and phase-I motions

The switching conditions between phase-II and phase-I are similar to that of phase-II and phase-III. Let $\tau_{c4}$ be the switching time between phase-II and phase-I motions. When the string completely unwraps from the obstacle during phase-II motion the slope of the string at $x=\gamma=0$ matches with the slope of the obstacle and again phase-I motion gets initiated. Following similar procedure as in Sec. II E, we get the following relations:

$$\beta_i(\tau_{c4}) = \eta_i(\tau_{c4}), \qquad (57)$$

$$\dot{\beta}_i(\tau_{c4}) = \dot{\eta}_i(\tau_{c4}). \qquad (58)$$

The transfer time $\tau_{c4}$ between phase-II and phase-III can be again obtained from the transversality condition as shown below:

$$\sum_{i=1}^{N} i\pi\sqrt{2}\beta_i(\tau_{c4}) = \alpha b. \qquad (59)$$

## H. Summary of formulation

Now we have the equations governing the dynamics of the string during the three phases of motion given by Eqs. (8), (30), and (48). The switching conditions between phase-I and phase-II are given by Eqs. (38) and (43). The event of switching can be obtained from Eq. (44). The switching conditions between phase-II and phase-III motions are given by Eqs. (49) and (50). The event of switching can be obtained from Eq. (52). During upward motion of the string, the switching conditions between phase-III and phase-II and the event of switching can be obtained from equations (53)–(56). Finally, the switching condition between phase-II and phase-I and the event of switching are given by Eqs. (57)–(59). In Sec. III, we discuss the results obtained by numerical simulations of the formulated equations.

## III. RESULTS AND DISCUSSION

In this section, we discuss the behavior of the string motion from simulations. For computational simplicity, we consider only a single mode representation of the string. It will be evident shortly that even a one mode approximation of the moving boundary formulation can capture the physics of the problem. This approximation requires the string to be plucked at the center of its unwrapped length.

### A. General behavior of string motion

As we have introduced dimensionless quantities in the equation of motion, the natural frequencies of the completely unwrapped string are integer multiples of $\pi$. For a simulation we need two dimensionless quantities: the relationship between the bridge and string given by $b=B/L$, and the modal amplitude of the initial string configuration given by $\beta_1(0)$. The contact length $\gamma=1-\sqrt{(1-b)}$ for the string at static equilibrium on the parabolic obstacle can be obtained from the transversality condition [Eq. (34)] by setting $\beta_1(0)=0$; the string shape at static equilibrium, a straight line from the contact point to right termination, is given by Eqs. (23) and (24) with $y_3=0$.

Transitions between phases are controlled by two factors: (i) the location of the bridge terminations, as determined by $\gamma_1$ and $\gamma_2$; and (ii) the pluck amplitude, as given by the initial condition $\beta_1(0)$. The value of $\gamma_1$ constrains the left boundary of the vibrating string and shortens the effective speaking length at the phase-I transition, and $\gamma_2$ constrains the potential extent of the string wrapping before it enters into phase-III motion. For each defined pair of bridge terminations, limiting values of the initial condition such that the string only vibrates in phase-II motion can be obtained from the transversality condition [Eq. (34)]:

$$\frac{(1 - \gamma_2)^{3/2}}{\sqrt{2\pi}} \left( \alpha(b - 2\gamma_2) - \frac{\alpha\gamma_2(b - \gamma_2)}{(\gamma_2 - 1)} \right) < \beta_1(0) \\ < \frac{(1 - \gamma_1)^{3/2}}{\sqrt{2\pi}} \left( \alpha(b - 2\gamma_1) - \frac{\alpha\gamma_1(b - \gamma_1)}{(\gamma_1 - 1)} \right). \qquad (60)$$

The significance of these inequalities for sitar-bridge design will be demonstrated subsequently. A parabolic bridge that
<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Variation of initial condition with initial contact length.</figcaption></figure>

extends to the x-axis on both sides, as used in the derivation of equations of motion and switching conditions in Sec. II, corresponds to $\gamma_1=0$ and $\gamma_2=b$, in which case the inequalities constraining the string to phase-II motion simplify to $-4(1-b)^{3/2}/\sqrt{2\pi b} < \beta_1(0) < 1/\sqrt{2\pi b}$. The bridge configuration in Ref. 7 corresponds to $\gamma_1=0$ and $\gamma_2 \to \infty$.

Figure 4 shows the variation of initial condition $\beta_1(0)$ against initial contact length $\gamma(0)$ for a string in phase-II configuration with $b=0.05$, the same as considered by Burridge et al.⁷ The corresponding midpoint deflection $y_2(0.5,0)$ is also shown. Interestingly the initial conditions seem to vary almost linearly with contact length. For this bridge configuration, the phase transitions occur when $\gamma(0)$ is 0 or 0.05, giving $\beta_1(0)$ values of 18.0 and -16.7, and corresponding midpoint deflections of 25.5 and -24.1; these demand very large amplitude initial conditions if the string is to vibrate in phase-I or phase-III motion. Increasing $\gamma_1$ will constrain the $\gamma(0)$ value on the left for the phase-I transition; decreasing $\gamma_2$ will constrain the $\gamma(0)$ value on the right giving the phase-III transition. In this way bridge termination can be used to control the string amplitude required for the phase transitions.

Now we study the free vibrations of the string about the equilibrium state. The equations of motions were solved numerically using MATLAB with ode 23s solver. The in-built event detection algorithm in MATLAB was used for detecting events for switching between the three phases of motion. An absolute and relative tolerance of $10^{-9}$ was used in the numerical simulations. Figure 5 shows the variation of modal coordinate $\beta_1(t)$ with time for four different initial conditions. The corresponding phase space plots are given in Fig. 6. The first initial condition is $\beta_1(0)=24 > 1/\sqrt{2\pi b}$ and thus the string starts its motion in phase-I; the string eventually vibrates in all the three phases of motion, but the asymmetry seen in the phase plot should be noted. All the other initial conditions shown satisfy the inequalities $-4(1-b)^{3/2}/\sqrt{2\pi b} < \beta_1(0) < 1/\sqrt{2\pi b}$ and the initial string condition is in phase-II; in this case, the string remains only in phase-II motion; however, this cannot be concluded in general due to asymmetry. Figure 5 shows that the string starting in phase-I

<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) Variation of modal amplitude $\beta_1$ with time for four different initial conditions. Bridge terminations given by $\gamma_1=0$ and $\gamma_2=b=0.05$.</figcaption></figure>

has a higher oscillation frequency than that starting in phase-II, since some motion in phase-III with the highest frequency occurs during a cycle. The frequency of oscillation for all cases that remain in phase-II motion is essentially the same, possibly due to the near linear relationship of $\beta_1(0)$ and $\gamma(0)$, as shown in Fig. 4. The period of oscillation with no obstacle is 2, and it can be seen that this is increased due to the presence of the bridge.

In order to understand how the natural frequency of the system changes while it wraps around the obstacle, an instantaneous natural frequency can be defined by taking a square root of the coefficient of $\beta_1(\tau)$ in Eq. (30) and by dropping rate dependent terms. Figure 7 shows the variation of instantaneous natural frequency for the same four initial conditions. It is clear that the natural frequencies are time dependent.

To validate our results with the moving boundary formulation, we compared them to those obtained by the penalty approach¹⁴ and found good agreement when 100 modes were retained in the penalty method. Moreover, the moving

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Phase space for $\beta_1$ and $\dot{\beta}_1$ for four different initial conditions. Bridge terminations given by $\gamma_1=0$ and $\gamma_2=b=0.05$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) Variation of instantaneous natural frequency with time.</figcaption></figure>

boundary formulation simulated the string motion 50 times faster than the penalty approach. During phase-III the slope of the string at $\gamma=b$ is discontinuous, except at the beginning and end configuration. This non-smooth behavior in the string shape at $\gamma=b$ is exactly captured by the formulation and its discontinuous slope is shown in Fig. 8. This particular non-smooth behavior in string slope cannot be captured if the impact is modeled using a penalty approach,¹⁴ for which method the obstacle is assumed as a linear continuum of distributed springs. Usually a series solution in terms of normal modes of the classical string is sought and it is well known that the series solution converges very slowly in the presence of non-smooth displacements and is prone to Gibbs phenomena.²⁰ Capturing such discontinuities in derivatives of the spatial displacement exactly is still a challenging problem with a sparse modeling literature.²⁰⁻²²

In order to investigate the frequency components present in the string shape during phase-II motion, the shape of the

<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Slope of the string during phase-II and III motions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. (Color online) Waterfall plot showing the variation of string shape frequency components with time for motion in phase-II only.</figcaption></figure>

string during its entire motion (obtained by solving the moving boundary problem) is projected on the normal modes of the classical string. Mathematically,

$$y(x,t) = \sum_{j=1}^{N} a_j(t)\sin(j\pi x) \qquad (61)$$

represents the shape of the string including both the wrapped and unwrapped portions. The series in Eq. (61) can also be interpreted as the series solution of the penalty approach, for which¹⁴ the $a_j(t)$ are obtained by solving the differential equations

$$\ddot{a}_j(\tau) + \omega_j^2 a_j(\tau) = \Psi f(a_1, a_2, \dots, a_N), \quad j = 1, 2, \dots, N, \qquad (62)$$

with penalty parameter $\Psi$ and penetration function $f(a_1, a_2, \dots, a_N)$. The moving boundary approach of the present paper is the limiting case of the penalty method for $\Psi \to \infty$ and penetration function tending to zero corresponding to a rigid obstacle. Both the methods should give identical results in the limiting case.

Figure 9 is a waterfall plot obtained from Eq. (61) showing the variation of normalized participation factors $a_j(t)$ with time for the initial condition $\beta_1(0)=1.8$ with string motion constrained to phase-II only. The presence of a particular $a_j(\tau)$ in the waterfall plot means that the corresponding modal oscillator in Eq. (62) must participate in the response if the problem is solved using a penalty approach, and automatically its frequency component will be present in the time response. Figure 9 clearly shows the participation of higher modes during phase-II motion. Considerably more terms in the series [Eq. (61)] are needed around $\tau=0.5$ and $\tau=1.5$ when the shape of the string demands higher mode participation. It can be seen that Eq. (30) is highly coupled and during phase-II motion the modes can exchange energy. As time progresses higher modes will start participating and eventually lead to multiple distributed impacts between the bridge and the string, which will violate the perfect wrapping assumption required in the moving boundary method. The
authors believe that these multiple impacts are responsible for the distinct tone of the sitar. Future work will address this by developing a detailed simulation model that can handle spatial string motions with multiple impacts and friction.²³

## B. Simulation approximating the configuration of a sitar

An initial string deflection of 24 times the height of the bridge is required for phase-I string motion with the bridge-string configuration used in the above simulations, which is the same as that of Burridge et al.⁷ This scenario is clearly impractical for a real sitar, if only because the resulting large amplitude string vibrations would have to pass through the back of the instrument. In reality, the geometry of a sitar-bridge, as shown in Fig. 2, is quite different in several respects: (i) the bridge is terminated on the left at its apex, with downbearing from string back length keeping it fixed there; (ii) the bridge is terminated on the right at a level considerably higher than that of the far string termination; and (iii) the slope of the bridge is very shallow, so the bridge surface remains very close to the string.

Measurements obtained from the sitar shown in Fig. 1 give the following dimensions: $B=300$ mm, $\Gamma_1=150$ mm, $\Gamma_2=173$ mm, and $L=1060$ mm (effective speaking length is 910 mm). We choose the straight neck of the instrument to define the horizontal direction. The x-axis for simulations (Fig. 3) is parallel to this and passes through the far string termination which is 14 mm above the neck reference line. The apex of the bridge at 27 mm above the neck reference line gives $h=13$ mm. Figures 2 and 3 show how the simulation configuration relates to the real sitar-bridge.

The tops of the frets at 11 mm above the neck reference line constrain the maximum vertical displacement of the string if it is to avoid hitting them. The normal plucking point is about 200 mm from the bridge apex, giving a plucking ratio or about 2:9. A typical pluck moves the string about 15 mm horizontally and 2–3 mm vertically.

This sitar string-bridge configuration is approximated for simulations using the following non-dimensional parameters: $L=1$, $b=0.283$, $\gamma_1=b/2=0.142$, and $\gamma_2=0.163$. The initial condition for a one mode solution requires a mid-string plucking point for which the $\beta_1(0)$ value of 0.25 is used. This corresponds to a string raised slightly above the horizontal between the bridge termination and the pluck point, a state that is easily achieved in normal playing.²⁴ Substituting the non-dimensional parameters in Eq. (60) gives phase transitions for $\beta_1$ values of 0.209 (phase-I to phase-II) and -0.172 (phase-II to phase-III), corresponding to midpoint deflections of 0.890 and 0.329, respectively. It can be seen that string motion in all three phases can easily be achieved for the configuration of a real sitar in normal playing, as a result of the geometry of the bridge and its terminations. The results of simulations with the above conditions, shown in the phase plot of Fig. 10 should be contrasted with those shown in Fig. 6.

<figure><img src="image_11.png" /><figcaption>FIG. 10. (Color online) Phase space for $\beta_1$ and $\dot{\beta}_1$ with two initial conditions and simulation parameters approximating those for a sitar: $b=0.283$, and bridge terminations $\gamma_1=0.142$, $\gamma_2=0.163$. Also shown are results for extra wide 30 mm bridge surface with $\gamma_2=0.170$.</figcaption></figure>

## IV. CONCLUSIONS

A mathematical model of the string wrapping against obstacle at its boundary has been formulated using a moving boundary approach. The formulation includes the distributed behaviour of a rigid bridge obstacle which may be terminated at arbitrary locations on either side. Equations of motion have been provided for the three phases of motion corresponding to the string completely (phase-III), partially (phase-II), or not at all (phase-I) wrapped on the bridge. It is shown that a single mode moving boundary approach can reveal much of the underlying physics, including capturing the non-smooth string shape during phase-II motion. As many as 60 natural frequency components of the string are present in the wrapped string, in particular, during phase-II motion. Thus the model captures the characteristic buzzing behavior of the sitar tone. In the simulations given the string motion has been reasonably well represented using only a single mode, requiring solution of only a single ordinary differential equation (ODE) in phase-I and phase-III, or two coupled ODEs in phase-II. This suggests the applicability of the method to physics-based sound synthesis algorithms. The following conclusions can be drawn on phase-II motion:

* The modal amplitude $\beta_1(\tau)$ decreases as the contact length $\gamma(\tau)$ increases.

* The frequency of oscillation of the string initially in phase-I is higher than that initially in phase-II.

* The frequency of oscillation in phase-II remains constant irrespective of initial amplitude given by modal amplitude $\beta_1(0)$.

<footnote n="1">¹B. C. Deva, *Musical Instruments of India: Their History and Development*, 2nd ed. (Munshiram Manoharlal, New Delhi, India, 1987).</footnote>

<footnote n="2">²C. Fulton, "Playing the late medieval harp," in *A Performer's Guide to Medieval Music*, edited by R. Duffin (Indiana University Press, Bloomington IN, 2002), pp. 346–354.</footnote>

<footnote n="3">³G. O'Brien, *Ruckers: A Harpsichord and Virginal Building Tradition* (Cambridge University Press, Cambridge, England, 1990).</footnote>

<footnote n="4">⁴C. V. Raman, "On some Indian stringed instruments," *Proceedings of the Indian Association for the Cultivation of Science* 7, 29-33 (1922).</footnote>
⁵L. Amerio and G. Prouse, "Study of the motion of a string vibrating against an obstacle," Rend. Mat. **2**, 563-585 (1975).

⁶M. Schatzman, "A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle," J. Math. Anal. Appl. **73**, 138-191 (1980).

⁷R. Burridge, J. Kappraff, and C. Morshedi, "The sitar string, a vibrating string with a one-sided inelastic constraint," SIAM J. Appl. Math. **42**, 1231-1251 (1982).

⁸H. Cabannes, "Motion of a string in the presence of a straight rectilinear obstacle," C.R. Acad. Sc. Paris, Series II **295**, 637-640 (1982).

⁹H. Cabannes, "Motion of a vibrating string in the presence of a convex obstacle: A free boundary problem," C.R. Acad. Se. Paris, Series II **301**, 125-129 (1985).

¹⁰H. Cabannes, "Presentation of software for movies of vibrating strings with obstacles," Appl. Math. Lett. **10**, 79-94 (1997).

¹¹J. Ahn, "A vibrating string with dynamic frictionless impact," Appl. Num. Math. **57**, 861-884 (2007).

¹²S. M. Han and M. A. Grosenbaugh, "Non-linear free vibration of a cable against a straight obstacle," J. Sound Vib. **273**, 337-361 (2004).

¹³T. Taguti, "Dynamics of simple string subject to unilateral constraint: A model analysis of sawari mechanism," Acoust. Sci. & Tech. **29**, 203-214 (2008).

¹⁴C. P. Vyasarayani, S. Birkett, and J. McPhee, "Free vibration response of a plucked string impacting against a spatial obstacle," in Proceedings of the International Conference on Theoretical, Applied, Computational and Experimental Mechanics, Kharagpur, India (2007).

¹⁵A. Krishnaswamy and J. O. Smith, "Methods for simulating string collisions with rigid spatial obstacles," in IEEE Workshop on Applications to Signal Processing to Audio and Acoustics. (2003), pp. 233-236.

¹⁶C. Valette, "The mechanics of vibrating strings," in *Mechanics of Musical Instruments*, edited by J. Kergomard and G. Weinreich (Springer-Verlag, Vienna, 1995), pp. 115-183.

¹⁷In the literature, a perfectly wrapped string is often described as a completely plastic impact.

¹⁸R. Fung and Y. Chen, "Free and forced vibration of a cantilever beam contacting with a rigid cylindrical foundation," J. Sound Vib. **202**, 161-185 (1997).

¹⁹F. Pfeiffer and C. Glocker, *Multibody Dynamics With Unilateral Contacts* (Wiley, New York, 1996).

²⁰A. V. Pesterev and L. A. Bergman, "An improved series expansion of the solution to the moving oscillator problem," ASME J. Vibr. Acoust. **122**, 54-61 (2000).

²¹W. D. Zhu and C. D. Mote, "Dynamics of the pianoforte string and narrow hammers," J. Acoust. Soc. Am. **96**, 1999-2007 (1994).

²²S. Bilbao, "Conservative numerical methods for nonlinear strings," J. Acoust. Soc. Am. **118**, 3316-3327 (2005).

²³C. P. Vyasarayani, S. Birkett, and J. McPhee, "A vibrating sitar string: Modeling the 3d dynamics of a plucked string impacting a spatial obstacle with friction," in Acoustics Paris '08 (2008).

²⁴These comments refer to an open string. Fretting the string shortens its length and lowers the far termination about 2 mm, making only a minor difference to the vertical displacement at the normal plucking point.