ABSTRACT
A threedimensional coupled vibroacoustic finite element model for physicsbased simulations of the sound generation by mallet percussion instruments in the time domain is discussed in the present paper. The mechanical model takes the orthotropic material properties of the wooden sound bars and the nonlinear nature of the interaction force between the mallet head and the sound bar into account while the acoustical model considers radiation into an unbounded domain. A direct coupling of the sound bars, acoustical cavity resonators, and the excitation by a mallet is considered with exploiting the modal basis to reduce the number of degrees of freedom of the system. Both the mechanical and acoustical models are validated by comparing them to measurements performed on an Orff xylophone. A case study shows the capabilities of the coupled model, including the analysis of the energy balance, the effect of tuning the resonator, and the excitation of the torsional modes of the sound bar.
Mallet percussion instruments belong to the family of tuned idiophones and produce sound by the radiation of their struck sound bars.^{1} The strike by the mallet excites the vertical bending modes (see the y axis in Fig. 1) of vibration of the sound bar, and the eigenfrequencies of these modes are the dominant partials of the radiated sound. Because the sound bar in itself is an inefficient radiator, tubular or cavity resonators are used in different instruments for reinforcing the radiated sound.
As a bar with a uniform cross section produces nonharmonic eigenfrequencies, the vertical bending modes of the sound bars are tuned by a characteristic undercut to attain a harmonic ratio of $1:4:10$ or $1:3:9$ of the frequencies of the first three partials. The tuning of marimba and xylophone bars was investigated in detail by Bork^{2,3} and many practical aspects were explored. The properties of mallets were also investigated by Bork^{4} by examining the shock spectrum using force transducers. OrduñaBustamente^{5} and Petrolito and Legge^{6} developed numerical methods for optimizing the undercuts to attain a given harmonic ratio of the partials. Recently, Beaton and Scavone^{7} examined iterative and genetic algorithms for the optimization.
Numerical simulations of xylophone bars based on a onedimensional finite difference method were introduced by Chaigne and Doutaut^{8} where the nonlinear interaction of the mallet and the sound bar was also examined in detail by both experiments and simulations. Then, the model was extended by Doutaut et al.^{9} by a tubular resonator that was also described as a unidimensional system. The latter contribution assumed that the backcoupling effect of the resonator on the bar is negligible, which is a specific property of xylophone bars in the higher frequency range, and cannot be generalized to all mallet percussion instruments as also remarked by the authors. Henrique and Antunes^{10} presented a modal simulation of the vibration of sound bars, taking nonlinear interaction phenomena into account; however, acoustical radiation was not considered.
Creating a physicsbased numerical model of the sound production by a mallet percussion instrument involves various challenges. The wooden sound bar has to be represented by an orthotropic material model as shown by Bork et al.^{11} As the instrument radiates into an open or semiopen space, the sound radiation model has to be able to treat unbounded domains. Finally, the nonlinear nature of the interaction between the mallet and sound bar renders solutions in the frequency domain infeasible. Whereas this paper is restricted to the collision model of Chaigne and Doutaut,^{8} the numerical modeling of collisions in a musical instrument is examined in a much wider scope by Bilbao et al.^{12}
In this contribution, a finite element (FE) approach that is capable of combating the challenges mentioned above is presented. The objective of this paper is to introduce the model, compare it with measurements, and show its capabilities by means of a practical case study. The proposed approach relies on the threedimensional (3D) finite element method (FEM), which is applied for modeling both the sound bar and its surrounding acoustical environment. The mechanical model takes the orthotropic material properties and viscoelastic losses into account and is capable of computing all mode shapes of the sound bars. Acoustical radiation into the free field is modeled by means of infinite elements. A direct coupling between the two subsystems is computed by making use of the conservative interpolation on noncompatible meshes. The modal basis is exploited for efficiently reducing the size of the system matrices. Finally, the system is transformed into the time domain and solved by a time stepping scheme. The validation measurements and the case study are performed on an Orff xylophone.
The remainder of the paper is structured as follows. Section II introduces the FE model, including the mechanical and acoustical systems, their coupling, and the applied time stepping scheme. The mechanical and acoustical FE models are validated by comparing the simulation results to measurements in Sec. III. A case study is presented in Sec. IV, highlighting the capabilities of the proposed coupled model. Section V concludes the paper with a short discussion and an outlook on possible further examinations.
A. Mechanical model
A 3D mechanical model, based on the Navier–Cauchy equations for linear elasticity, is considered for describing the vibrations of the sound bars. Whereas the vibration of the vertical bending modes can also be represented using a simple, onedimensional beam model,^{8,9} the 3D model is also capable of capturing horizontal and torsional modes. The equation of motion in the solid is expressed as
where $\mathit{\sigma}$ is the stress tensor, ρ is the density of the material, u is the displacement, x is the spatial coordinate, and ${\Omega}_{\mathrm{m}}$ is the mechanical domain, i.e., the sound bar in our case. The angular frequency is ω and timeharmonic functions with a time dependence of ${\mathrm{e}}^{\mathrm{j}\omega t}$ are assumed with $\mathrm{j}$ being the imaginary unit. The body forces are neglected hereafter.
$$\nabla \cdot \mathit{\sigma}(\mathit{x})={\omega}^{2}\rho \mathit{u}(\mathit{x}),\hspace{1em}\mathit{x}\in {\Omega}_{\mathrm{m}},$$  (1) 
The stress–strain relationship is given by the generalized Hooke's law, using a Newtonian linear viscoelastic model^{13} of intrinsic losses as
where $\mathsf{C}$ is the stiffness tensor of the material, η is the coefficient of viscous losses, and $\mathit{\epsilon}$ is the strain tensor, which is defined as
$$\mathit{\sigma}=\mathsf{C}:\left(\mathit{\epsilon}+\mathrm{j}\omega \eta \mathit{\epsilon}\right),$$  (2) 
$$\mathit{\epsilon}=\frac{1}{2}\left[\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^{\mathrm{T}}\right].$$  (3) 
The discretization of the weak form of Eqs. (1)–(3) by a standard Galerkin FEM leads to an algebraic system of equations, which are written as (see, e.g., Ref. 14)
The matrices ${\mathbf{K}}_{\mathrm{m}},\hspace{0.17em}{\mathbf{C}}_{\mathrm{m}}$, and ${\mathbf{M}}_{\mathrm{m}}$ are the mechanical stiffness, damping, and mass system matrices, respectively. The matrix ${\mathbf{G}}_{\mathrm{m}}$ results from a boundary integral appearing in the integral equation of the weak form. The subscript “m” refers to the mechanical subsystem hereafter. The vector u contains the coefficients of the displacement for each degree of freedom of the discretized system. Similarly, the vector f holds the coefficients of the discretized force distribution acting on the surface of the body.
$${\mathbf{K}}_{\mathrm{m}}\mathbf{u}+\mathrm{j}\omega {\mathbf{C}}_{\mathrm{m}}\mathbf{u}{\omega}^{2}{\mathbf{M}}_{\mathrm{m}}\mathbf{u}={\mathbf{G}}_{\mathrm{m}}\mathbf{f}.$$  (4) 
The applied viscoelastic model leads to a proportional damping provided that the properties of the material are homogeneous such that
From Eq. (2), ${a}_{\mathrm{m}}=\eta $ and ${b}_{\mathrm{m}}=0$ are found. It is noted that a nonzero constant ${b}_{\mathrm{m}}$ may be applied for incorporating the losses due to the air load into the mechanical model such as in Ref. 8. As a direct coupling with the acoustical field is computed here, the latter option is not used. Special properties of proportional damping are exploited in the modal analysis of the system.
$${\mathbf{C}}_{\mathrm{m}}={a}_{\mathrm{m}}{\mathbf{K}}_{\mathrm{m}}+{b}_{\mathrm{m}}{\mathbf{M}}_{\mathrm{m}}.$$  (5) 
B. Acoustical model
Propagation of acoustical waves is described by the Helmholtz equation. As the sound field is excited by the vibrating sound bar, the Euler equation is also used for prescribing the boundary conditions for the acoustical system. The two equations read as
where p is the sound pressure, $k=\omega /c$ is the acoustical wave number with c denoting the speed of sound, ρ_{0} is the equilibrium density of air, $\mathit{v}(\mathit{x})$ is the particle velocity, and ${\Omega}_{\mathrm{a}}$ is the acoustical domain, which contains both the interior of the resonator and the semiopen exterior space that surrounds the instrument.
$${\nabla}^{2}p(\mathit{x})+{k}^{2}p(\mathit{x})=0,\hspace{0.5em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\mathit{x}\in {\Omega}_{\mathrm{a}}$$  (6) 
$$\nabla p(\mathit{x})+\mathrm{j}\omega {\rho}_{0}\mathit{v}(\mathit{x})=0,\hspace{1em}\mathit{x}\in {\Omega}_{\mathrm{a}},$$  (7) 
The Helmholtz equation (6) can be discretized in a finite domain using a standard Galerkin type FE procedure. However, as the resonator and the sound bar radiate into the free acoustical field, special techniques are required to treat the unbounded domain in the formulation. The general approach is to truncate the computational domain and impose special boundary conditions at the truncated boundary ${\Gamma}_{\infty}$. There are various solutions available for constructing such boundary conditions, and among these, the infinite element method is used here, based on the formulation given by Astley et al.^{15,16} The acoustical FE system can be written in matrix form as
where ${\mathbf{K}}_{\mathrm{a}},\hspace{0.17em}{\mathbf{C}}_{\mathrm{a}}$, and ${\mathbf{M}}_{\mathrm{a}}$ are the acoustical stiffness, damping, and mass system matrices, respectively. Similar to Eq. (4), the matrix ${\mathbf{G}}_{\mathrm{a}}$ appears on the righthand side from a surface integral term. The vectors p and v contain the coefficients of the discretized sound pressure and the coefficients of the surface normal component of the particle velocity, respectively.
$${\mathbf{K}}_{\mathrm{a}}\mathbf{p}+\mathrm{j}\omega {\mathbf{C}}_{\mathrm{a}}\mathbf{p}{\omega}^{2}{\mathbf{M}}_{\mathrm{a}}\mathbf{p}=\mathrm{j}\omega {\mathbf{G}}_{\mathrm{a}}\mathbf{v},$$  (8) 
C. Mallet–sound bar interaction
In this section, the mallet–sound bar interaction is discussed, following Chaigne and Doutaut.^{8} The interaction model is based on Hertz's law of contact. Figure 1 illustrates the deformation of the mallet head during the contact and the resulting interaction force. Vibrations are also induced in the handle of the mallet; however, they are not significant during the contact due to the short interaction time and are neglected.^{8}
The compression of the mallet head ${\delta}_{\mathrm{M}}$ during the interaction with the sound bar is expressed from the interaction force F using Hertz's law of contact written for two colliding spheres as
where the constant D is calculated from the elastic material properties of the mallet head and the bar, and ${R}_{\mathrm{B}}$ and ${R}_{\mathrm{M}}$ denote the radii of the two spheres. The subscrpits “B” and “M” refer to the bar and the mallet, respectively. The equivalent radius of the sound bar considered as a sphere can be taken as infinite and, hence, the interaction force is expressed from the compression of the mallet head as
where ${K}_{\mathrm{M}}$ is the equivalent stiffness of the mallet head.
$${\delta}_{\mathrm{M}}={\left[{F}^{2}{D}^{2}\left(\frac{1}{{R}_{\mathrm{B}}}+\frac{1}{{R}_{\mathrm{M}}}\right)\right]}^{1/3},$$  (9) 
$$F=\frac{{\sqrt{R}}_{\mathrm{M}}}{D}{\delta}_{\mathrm{M}}^{3/2}={K}_{\mathrm{M}}{\delta}_{\mathrm{M}}^{3/2},$$  (10) 
It is assumed that the interaction force acts at a single degree of freedom (along the y axis at a single node of the mesh) B of the sound bar. This simplification is reasonable as the duration of the contact is dominant in determining the frequency content of the excitation and the size of the contact area has a negligible effect. Due to the finite resolution of the mesh, the effective area of the excitation ${S}_{\text{eff}}$ is nonzero and found as the sum of the Bth row of the matrix ${\mathbf{G}}_{\mathrm{m}}$. When the mallet–sound bar interaction is considered, the vector f in Eq. (4) is an allzero vector except for the single degree of freedom B where the force over the unit area of ${f}_{\mathrm{B}}=F/{S}_{\text{eff}}$ is present.
The equation of motion for the mallet head reads as
where the dot denotes the time derivative and ${u}_{\mathrm{M}}$ is the displacement of the mallet head. The interaction force acting on the mallet head ${F}_{\mathrm{M}}$ is expressed using Eq. (10) as
where ${u}_{\mathrm{B}}$ is the vertical displacement of the bar at the point of contact.
$${m}_{\mathrm{M}}{\ddot{u}}_{\mathrm{M}}={F}_{\mathrm{M}},$$  (11) 
$${F}_{\mathrm{M}}=\{\begin{array}{cc}{K}_{\mathrm{M}}{\left{u}_{\mathrm{B}}{u}_{\mathrm{M}}\right}^{3/2},\hfill & \text{if}\hspace{1em}{u}_{\mathrm{M}}<{u}_{\mathrm{B}},\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$  (12) 
D. The coupled system
To arrive at a fully coupled FE system, twoway interactions between the sound bar and sound field are taken into account in the following manner. First, the acoustical pressure exerts force on the surface of the sound bar in its normal direction. Second, the motion of the bar appears as a velocity boundary condition for the acoustical system. Only the motion in the normal direction generates acoustical waves, and the normal particle velocity is taken to be equal to the normal vibration velocity of the bar. Other interactions, such as the friction between the bar and surrounding air, are neglected.
The FE system of equations with the forcing terms appearing on the righthand sides of Eqs. (4) and (8) read as
$${\mathbf{K}}_{\mathrm{m}}\mathbf{u}+\mathrm{j}\omega {\mathbf{C}}_{\mathrm{m}}\mathbf{u}{\omega}^{2}{\mathbf{M}}_{\mathrm{m}}\mathbf{u}={\mathbf{G}}_{\mathrm{m}}\left(\mathbf{f}+{\mathbf{A}}_{\mathrm{m}}\mathbf{p}\right),$$  (13a) 
$${\mathbf{K}}_{\mathrm{a}}\mathbf{p}+\mathrm{j}\omega {\mathbf{C}}_{\mathrm{a}}\mathbf{p}{\omega}^{2}{\mathbf{M}}_{\mathrm{a}}\mathbf{p}=\mathrm{j}\omega {\mathbf{G}}_{\mathrm{a}}\left(\mathrm{j}\omega {\mathbf{A}}_{\mathrm{a}}\mathbf{u}\right).$$  (13b) 
(13)
The coupling matrices ${\mathbf{A}}_{\mathrm{m}}$ and ${\mathbf{A}}_{\mathrm{a}}$ are used for computing the pressure force and displacement in each nodal location of the mechanical and acoustical meshes, respectively. As the two meshes are not necessarily compatible, ${\mathbf{A}}_{\mathrm{m}}$ and ${\mathbf{A}}_{\mathrm{a}}$ are evaluated by means of conservative interpolation as illustrated in Fig. 2. The procedure consists of two steps. First, in the inverse mapping step [see Fig. 2(a)], the coordinates of each node of one mesh are computed in the local coordinate system of the elements of the other mesh. Then, in the second step, the shape functions are evaluated using the local coordinates, resulting in interpolation weights stored in the matrices ${\mathbf{A}}_{\mathrm{m}}$ and ${\mathbf{A}}_{\mathrm{a}}$. Finally, the function value at the interpolation point is found by summing the nodal values multiplied by the interpolation weights, i.e., a multiplication by ${\mathbf{A}}_{\mathrm{m}}$ or ${\mathbf{A}}_{\mathrm{a}}$ [see Fig. 2(b)]. The procedure preserves the surface integral of the interpolated quantities.^{17}
Rearranging the unknowns to the lefthand side and writing Eq. (13) as a block system gives
$$\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{m}}& {\mathbf{G}}_{\mathrm{m}}{\mathbf{A}}_{\mathrm{m}}\\ \mathbf{0}& {\mathbf{K}}_{\mathrm{a}}\end{array}\right]\left\{\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right\}+\mathrm{j}\omega \left[\begin{array}{cc}{\mathbf{C}}_{\mathrm{m}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_{\mathrm{a}}\end{array}\right]\left\{\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right\}{\omega}^{2}\left[\begin{array}{cc}{\mathbf{M}}_{\mathrm{m}}& \mathbf{0}\\ {\mathbf{G}}_{\mathrm{a}}{\mathbf{A}}_{\mathrm{a}}& {\mathbf{M}}_{\mathrm{a}}\end{array}\right]\left\{\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{G}}_{\mathrm{m}}\mathbf{f}\\ \mathbf{0}\end{array}\right\}.$$  (14) 
E. Modal approach
The modal basis is utilized for reducing the degrees of freedom of the coupled system. The mechanical modes are found as the eigenvectors of the undamped system by solving the generalized eigenvalue problem
for the eigenvalues ${\omega}_{\mathrm{m}}^{2}$ and eigenvectors ${\psi}_{\mathrm{m}}$. The mechanical eigen angular frequencies are ${\omega}_{\mathrm{m}}$. The damping factors ${\xi}_{\mathrm{m}}$ and the damped eigenfrequencies ${\omega}_{\mathrm{m}}^{\star}$ are found by using the proportional damping model [Eq. (5)] as
$$\left({\mathbf{K}}_{\mathrm{m}}{\omega}_{\mathrm{m}}^{2}{\mathbf{M}}_{\mathrm{m}}\right){\psi}_{\mathrm{m}}=\mathbf{0}$$  (15) 
$${\xi}_{\mathrm{m}}=\frac{1}{2}\left({a}_{\mathrm{m}}{\omega}_{\mathrm{m}}+\frac{{b}_{\mathrm{m}}}{{\omega}_{\mathrm{m}}}\right),\hspace{0.17em}\hspace{0.17em}{\omega}_{\mathrm{m}}^{\star}=\sqrt{1{\xi}_{\mathrm{m}}^{2}}.$$  (16) 
Whereas the mechanical damping matrix ${\mathbf{C}}_{\mathrm{m}}$ is given by the proportional model [Eq. (5)], the acoustical damping matrix ${\mathbf{C}}_{\mathrm{a}}$ depends on the geometrical arrangement of the infinite elements and is not proportional to the mass and stiffness matrices, which renders the computation of the acoustical modes more involved. To compute the acoustical eigenmodes ${\psi}_{\mathrm{a}}$ and the corresponding eigen angular frequencies ${\omega}_{\mathrm{a}}$, the following generalized eigenvalue problem^{18} needs to be solved:
The resulting eigenvalues are complex conjugate pairs^{18} ${\lambda}_{\mathrm{a}}={\xi}_{\mathrm{a}}{\omega}_{\mathrm{a}}\pm \mathrm{j}{\omega}_{\mathrm{a}}$, where ${\xi}_{\mathrm{a}}$ is the damping factor. The eigenvectors ${\psi}_{\mathrm{a}}$ represent propagating acoustical waves.
$$\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{a}}& {\mathbf{C}}_{\mathrm{a}}\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left\{\begin{array}{c}{\psi}_{\mathrm{a}}\\ {\lambda}_{\mathrm{a}}{\psi}_{\mathrm{a}}\end{array}\right\}={\lambda}_{\mathrm{a}}\left[\begin{array}{cc}\mathbf{0}& {\mathbf{M}}_{\mathrm{a}}\\ \mathbf{I}& \mathbf{0}\end{array}\right]\left\{\begin{array}{c}{\psi}_{\mathrm{a}}\\ {\lambda}_{\mathrm{a}}{\psi}_{\mathrm{a}}\end{array}\right\}.$$  (17) 
Using the eigenmodes, the coefficient vectors of both the displacement u and sound pressure p are written, respectively, in truncated modal bases
where the number of mechanical and acoustical modes are ${n}_{\mathrm{m}}$ and ${n}_{\mathrm{a}}$, respectively. The matrices ${\mathbf{\Psi}}_{\mathrm{m}}$ and ${\mathbf{\Psi}}_{\mathrm{a}}$ contain the eigenmodes in their columns, whereas ${\mathbf{q}}_{\mathrm{m}}$ and ${\mathbf{q}}_{\mathrm{a}}$ hold the corresponding modal weights. The truncation limit is chosen by prescribing an upper limit for the maximal eigenfrequency kept in the modal representation.^{19}
$$\mathbf{u}\approx \sum _{i=1}^{{n}_{\mathrm{m}}}{\mathbf{\psi}}_{\mathrm{m},i}{q}_{\mathrm{m},i}={\mathbf{\Psi}}_{\mathrm{m}}{\mathbf{q}}_{\mathrm{m}},$$  (18a) 
$$\mathbf{p}\approx \sum _{i=1}^{{n}_{\mathrm{a}}}{\mathbf{\psi}}_{\mathrm{a},i}{q}_{\mathrm{a},i}={\mathbf{\Psi}}_{\mathrm{a}}{\mathbf{q}}_{\mathrm{a}},$$  (18b) 
(18)
The coupled system of equations (14) is converted into the modal basis by applying the modal representation [Eq. (18)]. By introducing the block matrices,
and the vectors
the coupled modal system is written as
$$\mathbf{K}\hspace{0.17em}=\left[\begin{array}{cc}{\mathbf{\Psi}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{K}}_{\mathrm{m}}{\mathbf{\Psi}}_{\mathrm{m}}& {\mathbf{\Psi}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{G}}_{\mathrm{m}}{\mathbf{A}}_{\mathrm{m}}{\mathbf{\Psi}}_{\mathrm{a}}\\ \mathbf{0}& {\mathbf{\Psi}}_{\mathrm{a}}^{\mathrm{H}}{\mathbf{K}}_{\mathrm{a}}{\mathbf{\Psi}}_{\mathrm{a}}\end{array}\right],$$  (19a) 
$$\mathbf{C}\hspace{0.17em}=\left[\begin{array}{cc}{\mathbf{\Psi}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{C}}_{\mathrm{m}}{\mathbf{\Psi}}_{\mathrm{m}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Psi}}_{\mathrm{a}}^{\mathrm{H}}{\mathbf{C}}_{\mathrm{a}}{\mathbf{\Psi}}_{\mathrm{a}}\end{array}\right],$$  (19b) 
$$\mathbf{M}\hspace{0.17em}=\left[\begin{array}{cc}{\mathbf{\Psi}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{M}}_{\mathrm{m}}{\mathbf{\Psi}}_{\mathrm{m}}& \mathbf{0}\\ {\mathbf{\Psi}}_{\mathrm{a}}^{\mathrm{H}}{\mathbf{G}}_{\mathrm{a}}{\mathbf{A}}_{\mathrm{a}}{\mathbf{\Psi}}_{\mathrm{m}}& {\mathbf{\Psi}}_{\mathrm{a}}^{\mathrm{H}}{\mathbf{M}}_{\mathrm{a}}{\mathbf{\Psi}}_{\mathrm{a}}\end{array}\right],$$  (19c) 
$$\mathbf{q}=\left\{\begin{array}{c}{\mathbf{q}}_{\mathrm{m}}\\ {\mathbf{q}}_{\mathrm{a}}\end{array}\right\}\hspace{1em}\hspace{0.17em}\text{and}\hspace{1em}\hspace{0.17em}\mathbf{g}=\left\{\begin{array}{c}{\mathbf{\Psi}}_{\mathrm{m}}^{\mathrm{T}}{\mathbf{G}}_{\mathrm{m}}\mathbf{f}\\ \mathbf{0}\end{array}\right\},$$  (19d) 
$$\mathbf{Kq}+\mathrm{j}\omega \mathbf{Cq}{\omega}^{2}\mathbf{M}\mathbf{q}=\mathbf{g}.$$  (20) 
Despite that the modal system matrices have diagonal and full blocks as well, it is still less demanding to solve the modal system, thanks to the efficient reduction of the size of the system by the modal basis [Eq. (18)].
F. Time stepping scheme
As the displacement to the force relation [Eq. (10)] is nonlinear, the system cannot be directly solved in the frequency domain. Therefore, a suitable time stepping algorithm is introduced by writing Eq. (20) in the time domain using the inverse Fourier transform and coupling with Eq. (11) as
The lower indices of the vectors denote that Eq. (21) is written for the n + 1th time step. The coupled Eq. (21) is solved by applying the secondorder accurate Newmark time stepping scheme.^{20}
$$\mathbf{M}{\ddot{\mathbf{q}}}_{n+1}+\mathbf{C}{\stackrel{\u0307}{\mathbf{q}}}_{n+1}+\mathbf{K}{\mathbf{q}}_{n+1}={\mathbf{g}}_{n+1},$$  (21a) 
$${\mathrm{M}}_{\mathrm{M}}{\ddot{u}}_{\mathrm{M},n+1}={F}_{\mathrm{M},n+1}.$$  (21b) 
(21)
The following computations are performed in each time step. At the beginning of the n + 1th time step, it is assumed that the state of motion of both the mallet and sound bar is known. First, an approximation of the contact force ${\stackrel{\u0303}{F}}_{\mathrm{M}}$ is evaluated by substituting the result of the previous time step ${u}_{\mathrm{B},n}$ and ${u}_{\mathrm{M},n}$ into Eq. (12). The resulting force is inserted into Eq. (19d) to attain ${\mathbf{g}}_{n+1}$. Then, the time stepping scheme is applied in a predictor–corrector approach to the bar in the following steps:
(1)  Prediction step
 
(2)  Solution step
with
 
(3)  Correction step

$\Delta t$ denotes the time step size, which is not necessarily uniform throughout the simulation. The constants γ and β define a linear scaling between a fully explicit ($\beta =0,\gamma =0$) and a fully implicit ($\beta =1/2,\gamma =1$) scheme. A common choice that is also applied here is $\gamma =1/2$ and $\beta =1/4$. With these settings, the scheme is unconditionally stable.^{20}
The predictor step for the mallet reads
With ${\mathbf{q}}_{n+1}$ already attained, ${u}_{\mathrm{B},n+1}$ and ${\stackrel{\u0303}{u}}_{\mathrm{M}}$ are substituted into Eq. (12) to attain the interaction force ${F}_{\mathrm{M},n+1}$ and Eq. (21b) is applied to get ${\ddot{u}}_{\mathrm{M},n+1}$. Finally, the corrector step for the mallet head is performed,
$${\stackrel{\u0303}{\stackrel{\u0307}{u}}}_{\mathrm{M}}={\stackrel{\u0307}{u}}_{\mathrm{M},n}+\Delta t(1\gamma ){\ddot{u}}_{\mathrm{M},n},$$  (26a) 
$${\stackrel{\u0303}{u}}_{\mathrm{M}}={u}_{\mathrm{M},n}+\Delta t{\stackrel{\u0307}{u}}_{\mathrm{M},n}+\frac{\Delta {t}^{2}}{2}(12\beta ){\ddot{u}}_{\mathrm{M},n}.$$  (26b) 
$${\stackrel{\u0307}{u}}_{\mathrm{M},n+1}={\stackrel{\u0303}{\stackrel{\u0307}{u}}}_{\mathrm{M}}+\Delta t\gamma {\ddot{u}}_{\mathrm{M},n+1},$$  (27a) 
$${u}_{\mathrm{M},n+1}={\stackrel{\u0303}{u}}_{\mathrm{M}}+\Delta {t}^{2}\beta {\ddot{u}}_{\mathrm{M},n+1}.$$  (27b) 
The computationally demanding step of the scheme is solving Eq. (23). As the linear system of equations needs to be solved in a large number of time steps with the same coefficients ${\mathbf{M}}^{\star}$, it is useful to precompute and store the inverse of ${\mathbf{M}}^{\star}$. This precomputation is possible for the modal system [Eq. (20)], whereas it is prohibitively expensive in the case of the original system [Eq. (14)].
The accurate simulation of the mallet–sound bar interaction necessitates very small time steps due to the interaction force [Eq. (12)] depending on the displacement of both the sound bar and mallet head. After the contact, during free vibration, much larger time steps are affordable. Therefore, the time domain simulations are performed in two stages with two different time step sizes.
All numerical procedures introduced in this section were programmed in an inhouse MATLABbased software tool for FEs. Both the mechanical and acoustical models employ a nodal formulation using isoparametric elements with linear shape functions. Sixthorder Lagrange polynomials are used in the radial direction in the infinite element formulation.^{15} The computation of the eigenvalues and eigenvectors in Eqs. (15) and (17) is performed using a restarted Arnoldi iteration process^{21} implemented in the builtin MATLAB routine eigs.^{22}
The mechanical and acoustical FE models are validated by comparisons with measurements carried out on a Sonor Orff primary series soprano xylophone. The instrument has 13 rosewood sound bars, tuned in a diatonic scale covering the range $\mathrm{C}4$–$\mathrm{A}5$. To allow for playing different scales, the sound bars of notes F4, B4, and $\mathrm{F}5$ can be replaced by $\mathrm{F}\u266f4,\hspace{0.17em}\mathrm{B}\u266d4$, and $\mathrm{F}\u266f5$.
The sound produced by striking the sound bars with the mallet is reinforced by cavity resonators. There are a total of six boxshaped resonators with an isosceles trapezoid base, covered by 2–3–2–2–2–2 sound bars from the lowest to the highest note. Each resonator has a rectangular opening with rounded corners at the center of its top plate. The resonator bodies are made of plastic. The sound bars lay on a felt support attached to the top plate of the resonators.
A. Material properties
To investigate the material properties of the rosewood sound bars, eight untuned samples of the same material were measured. Each sample had a length of 270 mm, a width of 31 mm, and a thickness of 16 mm. The latter two dimensions are the same as those of the actual sound bars. The average density of the samples was found as $\rho =1116\hspace{0.17em}{\text{kg}/\mathrm{m}}^{3}$ with a standard deviation of $20\hspace{0.17em}{\text{kg}/\mathrm{m}}^{3}$.
The vibrations of the samples were measured by striking them using a felt covered plastic mallet with the bars resting on wedge shaped felt supports. A GRAS AF46 1/2 in. condenser microphone (Holte, Denmark) measured the radiated sound pressure. After amplification by a B and K 2690–0S2 Nexus amplifier (Nærum, Denmark) and bandpass filtering to 500 Hz–10 kHz by an Ithaco 4213 filter (Ithaca NY), the spectrum of the microphone signal was measured by an HP 35670A FFT Dynamic Signal Analyzer (Santa Rosa, CA). For one measurement, the spectra of 20 strikes were averaged. The responses to strikes at different positions (top, side, edge) were recorded, allowing the identification of different modes in the spectra.
The rosewood bars are modeled by a homogeneous orthotropic material with the following properties, based on Ref. 23. Elastic moduli, ${E}_{L}=24.0\hspace{0.17em}\text{GPa},\hspace{0.17em}{E}_{T}/{E}_{L}=0.086,{E}_{R}/{E}_{L}=0.097$; shear moduli, ${G}_{TR}/{E}_{L}={G}_{RL}/{E}_{L}=0.147$ and ${G}_{LT}/{E}_{L}=0.097$; Poisson ratios, ${\nu}_{TR}=0.303,{\nu}_{RL}=0.077$, and ${\nu}_{LT}=0.428$. Here, the subscripts L, T, and R refer to the longitudinal, tangential, and radial directions of the material, which correspond to the coordinates along the lengths (x), widths (z), and thicknesses (y) of the sound bars, respectively. For the sake of comparison, a substitute isotropic material with a Young modulus E = 24 GPa and Poisson ratio $\nu =0.4$ was also used in the FE models. A regular brick mesh consisting of hexahedra elements with a maximum edge length of 2 mm was created. It was verified that further increasing the resolution of the mesh does not affect the eigenfrequencies in the frequency range of interest.
Table I shows the measured eigenfrequencies of the first few modes of the untuned sample bars compared to the FE results. In case of the FEM, the eigenfrequencies are evaluated by solving Eq. (15). The measured eigenfrequencies of the samples also varied to some extent, therefore, data of one sample whose eigenfrequencies were close to the average values are shown. As seen, the orthotropic FEM gives a good agreement with the measured values for all three types of bending modes with the error increasing slightly with the frequency. The isotropic material model provides an acceptable approximation for the vertical bending modes, but horizontal and torsional modes cannot be computed accurately. This justifies the use of the more complex orthotropic material model.
Mode  #  Measurement  FEM, isotropic  FEM, orthotropic  

shape  f (Hz)  f (Hz)  e (%)  f (Hz)  e (%)  
Vertical  1  1040  1040  0.0  1027  −1.3 
2  2728  2805  2.8  2712  −0.6  
3  4976  5339  7.3  5019  0.9  
4  7320  8511  16.3  7755  9.7  
Horizontal  1  1796  1940  8.0  1842  2.6 
2  4160  4968  19.4  4281  2.9  
3  6968  8918  28.0  7025  0.8  
4  9380  13 411  43.0  9795  4.4  
Torsional  1  2064  3886  88.3  2123  2.9 
2  4232  7789  84.1  4314  1.9  
3  6504  11 726  80.3  6624  1.8  
4  8544  15 710  83.9  9085  6.3 
B. Modal behavior of the sound bars
The sound bars are tuned by tuning cuts, which have the typical shape observed on marimba bars. The length of the sound bars changes from 325 mm to 228 mm along the scale. The tuning cuts that are made by handwork using a grindstone with a 125 mm diameter reduce the thickness of the bar at its middle to 6.8–8.2 mm. In case of the lower notes, the tuning cut has a straight section in the middle of the bar whose relative length decreases along the scale and completely disappears in the case of the $\mathrm{A}5$ bar. The meshes of the sound bars are created by subjecting a brick mesh of hexahedra elements (edge length 2 mm) to an appropriate geometrical transformation, leading to structured meshes consisting of $17\hspace{0.17em}500$–$25\hspace{0.17em}000$ nodes.
As the first three vertical bending modes are of key importance from a musical point of view, Table II shows the comparison of the measured and simulated eigenfrequencies of these modes. It was found that using the reference material properties gives large deviations from the measured frequencies for a number of sound bars, indicating that these properties vary from bar to bar. To mitigate this uncertainty, the following approach was pursued. The elastic moduli of the material were scaled by a factor of $E\prime /E$ such that the frequency of the first vertical mode matches the measured value within $\pm 0.5\hspace{0.17em}\text{Hz}$. In case of the orthotropic model, the scaling affects all elastic and shear moduli to the same extent. The accuracy of the FE model can then be assessed by comparing the frequencies of the second and third modes.
Measurement  FEM, isotropic  FEM, orthotropic  

f_{1}  f_{2}  f_{3}  $E\prime /E$  e_{2}  e_{3}  $E\prime /E$  e_{2}  e_{3}  
Note  (Hz)  (Hz)  (Hz)  —  (%)  (%)  —  (%)  (%) 
$\mathrm{C}4$  264  1048  2656  0.685  −5.5  −7.3  0.715  −7.3  −10.1 
$\mathrm{D}4$  296  1176  2788  0.695  −5.6  −3.2  0.723  −7.6  −6.2 
$\mathrm{E}4$  324  1320  3196  1.000  −2.3  3.6  1.050  −4.5  0.5 
$\mathrm{F}4$  349  1400  3376  0.950  0.2  3.8  0.997  −2.2  0.9 
$\mathrm{F}\u266f4$  369  1480  3376  1.000  0.0  7.6  1.048  −2.4  4.7 
$\mathrm{G}4$  392  1570  3530  1.000  −0.3  8.5  1.052  −2.6  5.6 
$\mathrm{A}4$  440  1762  3876  1.065  1.0  9.3  1.120  −1.6  6.5 
$\mathrm{B}\u266d4$  464  1864  4104  1.145  3.7  12.0  1.208  0.8  9.4 
$\mathrm{B}4$  496  1976  4216  1.155  3.5  11.8  1.220  0.7  9.4 
$\mathrm{C}5$  520  2096  4304  0.990  2.1  8.7  1.047  −0.7  6.6 
$\mathrm{D}5$  584  2300  4368  0.945  2.9  11.9  1.002  0.2  10.0 
$\mathrm{E}5$  659  2523  4733  1.008  6.7  12.2  1.076  4.0  10.7 
$\mathrm{F}5$  705  2755  5176  0.955  6.5  6.6  1.032  4.1  5.9 
$\mathrm{F}\u266f5$  742  2816  4984  0.920  6.0  8.5  0.997  4.1  5.0 
$\mathrm{G}5$  784  2984  5360  0.868  3.7  4.7  0.940  1.7  3.7 
$\mathrm{A}5$  880  3168  5308  0.820  3.2  11.0  0.884  1.0  9.1 
Mean/abs. errors  0.950  3.3  8.2  1.007  2.8  6.7 
As seen from Table II, the average scaling factors for the elastic moduli are close to one, i.e., the reference material represents the average material properties quite well. However, in the case of the $\mathrm{C}4$ and $\mathrm{D}4$ sound bars, the elastic moduli needed to be decreased by $\approx 30\%$, whereas an increase of $\approx 20\%$ was necessary in the case of the $\mathrm{B}\u266d4$ and $\mathrm{B}4$ bars. Overall, a good agreement of the measurement and FEM is found. Except for the $\mathrm{C}4$ and $\mathrm{D}4$ sound bars, the results of the orthotropic model are slightly closer to the measurements with the average absolute errors for the second and third vertical bending eigenfrequencies being 2.8% ($47\hspace{0.17em}\text{cent}$) and 6.7% ($112\hspace{0.17em}\text{cent}$), respectively.
No systematic error is visible in the data presented in Tables I and II. Thus, the observed discrepancies are attributed to the variations of the material properties and not to errors of the applied numerical techniques.
Figure 3 shows the relative eigenfrequencies of the first few modes, computed by the orthotropic FEM, normalized by the frequency of the first vertical mode of each sound bar. The attempt to keep the $1:4:10$ ratio of the first three vertical modes is clearly visible. The $1:4$ ratio is achieved quite well in the whole range, whereas the $1:10$ ratio decreases gradually to $1:7$ in the upper range $\mathrm{B}4$–$\mathrm{A}5$. This result is explained by the fact that to keep the widths and thicknesses of the bars constant along the scale, the length of the undercut must be decreased. As also visible in Fig. 3, the line of the second torsional mode crosses that of the third vertical mode at the $\mathrm{D}5$ note. Such overlapping of vertical and torsional modes may lead to undesired coupling effects and a poor quality sound for the given sound bar.
C. Resonators
The resonance properties of the six cavity resonators were measured by removing all sound bars and exciting the resonators using an external loudspeaker driven by a chirp signal. Only one resonator was excited one time, the others were closed by covering them using a sheet of plastic. A reference signal was recorded outside the resonator to compensate for the frequency response of the speaker. The response of the resonator was measured inside the cavity using a Knowles FG23629P16 miniature electret microphone (Itasca, IL). In the acoustical FE model, the source is an ideal point source, and the response was evaluated at the center of the bottom of the resonator.
The dimensions and the measured eigenfrequencies are listed in Table III. A very short neck of 3 mm length is formed at the opening of each resonator, whose length is corrected by the end correction $\Delta l$ to attain the effective length of the neck ${l}_{\text{eff}}$. The low frequency length correction value is found as $\Delta {l}_{0}=2\times 0.82\sqrt{S/\pi}$ with S denoting the area of the opening. The values presented in Table III also take the frequency dependence of the correction into account, following Ref. 24.
Dimensions  Meas.  Theory  FEM  

Resonator  V  S  ${l}_{\text{eff}}$  ${f}_{\text{res}}$  ${f}_{\text{res}}$  Error  ${f}_{\text{res}}$  Error 
(cm^{3})  (cm^{2})  (mm)  (Hz)  (Hz)  (%)  (Hz)  (%)  
1  1830.4  33.56  55.8  305.5  313.2  2.5  326.0  6.7 
2  1560.0  49.91  66.1  386.5  380.2  −1.6  389.5  0.8 
3  627.8  34.72  54.8  559.0  548.8  −1.8  554.0  −0.9 
4  531.0  37.62  55.9  641.5  615.0  −4.1  638.5  −0.5 
5  415.5  40.52  56.6  740.5  717.5  −3.1  751.5  1.5 
6  259.9  46.32  55.9  1 007.5  975.5  −3.2  1 031.5  2.4 
Table III also compares the measured, theoretical, and simulated values of the eigenfrequencies of the six resonators. As seen, both the theory and numerical models agree quite well with the measurements. In most cases, the FE result is closer to the measured value except in the case of the first resonator where the largest error is observed.
As also reported by Bork,^{3} the presence of the sound bars detunes the resonators to a significant extent. To investigate this effect, the resonators were also measured with all of the sound bars in place and damped using plasticine. In the FE simulation, the bars are cut out of the acoustical domain and all of their boundaries are regarded as perfectly rigid.
Table IV shows the simulated and measured detuning effects of the sound bars. In the simulations, the sound bars were added gradually, keeping a symmetrical arrangement in the case of the middle resonators. It is interesting to observe that the greatest changes of frequencies occur when the second sound bar is added. This is explained as the second bar reduces the relative amount of the open surface at the opening of the resonator to the greatest extent. Adding the neighboring sound bars still has a significant effect, whereas adding further bars results only in minor shifts of the frequency. As seen, the detuning predicted by the FEM is very close to that measured with resonator 2 being the only exception. The detuning can be as large as 15% ($240\hspace{0.17em}\text{cent}$), and this effect cannot be calculated theoretically. Therefore, a numerical approach is necessarily applied.
FEM  Meas.  

Res.  0 bars  1 bar  2 bars  3 bars  4 bars  5 bars  Detuning^{a}  
${f}_{\text{res}}$ (Hz)  $\Delta {f}_{\text{res}}$  (%)  
1  326.0  320.3  310.6  306.8  305.7  —  −6.2  −6.1  
2  389.5  383.7  —  365.4  —  357.5  −8.2  −14.9  
3  554.0  543.5  523.8  —  506.5  —  −8.6  −9.2  
4  638.5  623.7  594.9  —  570.2  —  −10.7  −12.8  
5  751.5  730.4  690.1  —  654.5  —  −12.9  −13.8  
6  1 031.5  998.1  920.3  885.1  872.1  —  −15.5  −15.5 
^{a}Change of the eigenfrequency due to the presence of sound bars. For the measurements, the setups without and with the sound bars are compared. For the FEM, the case with the highest number of bars (four or five) is compared to the zero bars case.
Besides the resonance frequencies, the damping of the natural resonances is also important from the viewpoint of the instrument design. Measuring the damping on the instrument was found to be challenging due to the inevitable excitation of the vibrations of the plastic walls of the resonator. Nevertheless, as can be seen in Table III, the resonators have relatively large opening surfaces and, hence, the damping due to radiation losses becomes too high, especially in the higher frequency range. Furthermore, the amplification of two notes having a major second or minor third distance by the same resonator necessitates a flat resonance curve, which results in too low a gain. Therefore, in the case studies presented in Sec. IV, a revised resonator design is examined. Additional illustrations of the xylophone model examined in this section, together with a sensitivity study of material parameters, are presented in the supplementary material.^{26}
In the sequel, time domain simulations on a revised xylophone model are discussed. The model shown in Fig. 4 contains only the $\mathrm{F}4$ ($\approx 349\hspace{0.17em}\text{Hz}$) note of the 13 sound bars and resonators. The geometry consists of (1) a Helmholtz cavity resonator with a circular opening (25 mm diameter), a very short neck (3 mm length), and an irregularly shaped base; (2) a hemisphere volume representing the exterior sound field; and (3) five sound bars. As was shown in Sec. III, the presence of neighboring sound bars can have a significant effect on the tuning of the resonator. Thus, besides the middle sound bar, which is struck by the mallet, the neighboring (passive) bars are also cut out from the acoustical domain ${\Omega}_{\mathrm{a}}$.
The acoustical mesh is created using the parametric mesh generation tool Gmsh^{25} and consists of linear tetrahedra. Infinite elements are attached to the truncated boundary ${\Gamma}_{\infty}$ shown in Fig. 4. The acoustical mesh has $\approx 100\hspace{0.17em}000$ nodes and $\approx 400\hspace{0.17em}000$ elements with a nonuniform spatial resolution. The finest element sizes of 1.0–2.5 mm are applied at the neck of the resonator and at the common boundaries with the sound bars, whereas the largest element sizes of 10 mm are used at the interface with infinite elements.
To attain the modal system [Eq. (20)], the mechanical modes are computed using Eq. (15). As the mechanical modes are relatively sparsely spaced in frequency, a few tens of modes are enough to represent the displacement of the sound bar in the frequency range of 0–20 kHz. The numerical evaluation of the acoustical eigenmodes using Eq. (17) is challenging as many nonphysical modes arise inevitably due to the artificial truncation of the computational domain at ${\Gamma}_{\infty}$. Spurious modes are filtered out using two criteria. First, modes that represent waves traveling inward from ${\Gamma}_{\infty}$ are dropped. Second, only modes with realistic damping factors ${10}^{4}<{\xi}_{\mathrm{a}}<{10}^{1}$ are kept. This method significantly reduces the number of modes with 14–22 eigenmodes (7–11 pairs) being kept in the 0–5 kHz frequency range.
Representative values of the most important parameters of the simulation setup are listed in Table V. The coeffient of viscous losses η and the properties of the mallet were set based on the parameters found in Refs. 3, 4, and 8. The sampling frequency of the time stepping algorithm is ${f}_{s}=48\hspace{0.17em}\text{kHz}$ ($\Delta t=20.83\hspace{0.17em}\mu \mathrm{s}$). During the short time of the mallet–sound bar contact (0.4–1.5 ms) a 100–200 times smaller time step size was found to be necessary to keep the energy balance of the discretized system during the interaction. The validation of the proposed time stepping scheme and visualization of the eigenmodes are presented in the supplementary material.^{26}
Component  Parameter  Symbol  Value 

Mallet  Mass  ${\mathrm{m}}_{\mathrm{M}}$  $35\times {10}^{3}\hspace{0.17em}\text{kg}$ 
Stiffness  K_{M}  $5\times {10}^{7}\hspace{0.17em}\mathrm{N}/{\mathrm{m}}^{3/2}$  
Initial velocity  ${v}_{\mathrm{M}}(0)$  1 m/s  
Interaction area  ${S}_{\text{eff}}$  3.88 mm^{2}  
Interaction time  ${T}_{\text{int}}$  0.7 ms  
Sound bar  Mass  m_{B}  $134.9\times {10}^{3}\hspace{0.17em}\text{kg}$ 
Coefficient of viscosity  η  $9.4\times {10}^{7}\hspace{0.17em}\mathrm{s}$  
Fundamental  ${f}_{\mathrm{m},1}$  349.0 Hz  
Damping  ${\xi}_{\mathrm{m},1}$  $1.03\times {10}^{3}$  
Number of modes  ${n}_{\mathrm{m}}$  34  
Resonator  Fundamental  ${f}_{\mathrm{a},1}$  349.1 Hz 
Damping  ${\xi}_{\mathrm{a},1}$  $8.8\times {10}^{3}$  
Number of modes  n_{a}  18 
A. The energy balance of the system
To assess the accuracy of the proposed time stepping scheme, it is useful to test the conservation of energy in the system. Furthermore, the analysis of the energy balance of the subsystems may also reveal interesting details of the physics of the process. Each of the three subsystems, i.e., the mallet, sound bar, and acoustical field, can store potential and kinetic energy. Due to viscoelastic losses and sound radiation, the sound bar and acoustical field also dissipate energy.
The potential and kinetic energies of the mallet head, ${E}_{\mathrm{M},p}$ and ${E}_{\mathrm{M},k}$, respectively, are found as
Similarly, the potential and kinetic energies stored in the sound bar, denoted by ${E}_{\mathrm{B},p}$ and ${E}_{\mathrm{B},k}$, are found as volume integrals of the strain and kinetic energy densities. The quantities are evaluated in the FE context using the system matrices,
The energy dissipated by the mechanical damping is the time integral of the dissipated power ${P}_{\mathrm{B}}$
$${E}_{\mathrm{M},p}={\int}_{0}^{{\delta}_{\mathrm{M}}}{F}_{\mathrm{M}}(s)\hspace{0.17em}\mathrm{d}s={\int}_{0}^{{\delta}_{\mathrm{M}}}{K}_{\mathrm{M}}{s}^{3/2}\hspace{0.17em}\mathrm{d}s=\frac{2}{5}{K}_{\mathrm{M}}{\delta}_{\mathrm{M}}^{5/2},$$  (28) 
$${E}_{\mathrm{M},k}=\frac{1}{2}{\mathrm{m}}_{\mathrm{M}}{v}_{\mathrm{M}}^{2}.$$  (29) 
$${E}_{\mathrm{B},p}(t)={\int}_{{\Omega}_{\mathrm{m}}}\sum _{i=1}^{3}\sum _{j=1}^{3}{\sigma}_{ij}{\epsilon}_{ij}\hspace{0.17em}\mathrm{d}V\approx \frac{1}{2}{\mathbf{u}}^{\mathrm{T}}(t){\mathbf{K}}_{\mathrm{m}}\mathbf{u}(t),$$  (30) 
$${E}_{\mathrm{B},k}(t)={\int}_{{\Omega}_{\mathrm{m}}}\frac{1}{2}\rho {\stackrel{\u0307}{u}}^{2}\hspace{0.17em}\mathrm{d}V\approx \frac{1}{2}{\stackrel{\u0307}{\mathbf{u}}}^{\mathrm{T}}(t){\mathbf{M}}_{\mathrm{m}}\stackrel{\u0307}{\mathbf{u}}(t).$$  (31) 
$${W}_{\mathrm{B}}(t)={\int}_{0}^{t}{P}_{\mathrm{B}}(\tau )\hspace{0.17em}\mathrm{d}\tau \approx {\int}_{0}^{t}{\stackrel{\u0307}{\mathbf{u}}}^{\mathrm{T}}(\tau ){\mathbf{C}}_{\mathrm{m}}\stackrel{\u0307}{\mathbf{u}}(\tau )\hspace{0.17em}\mathrm{d}\tau .$$  (32) 
The kinetic and potential energies stored in the acoustical system, ${E}_{\mathrm{A},p}$ and ${E}_{\mathrm{A},k}$, are found as volume integrals of the acoustical kinetic and potential energy densities of the sound field.
where ${\mathbf{M}}_{\mathrm{a}}^{\text{FE}}$ and ${\mathbf{K}}_{\mathrm{a}}^{\text{FE}}$ denote the acoustical mass and stiffness matrices, respectively, without the infinite element contributions. The time domain form of the Euler equation (7) was used and $\mathbf{h}(t)$ was introduced as
$${E}_{\mathrm{A},p}(t)={\int}_{{\Omega}_{\mathrm{a}}}\frac{1}{2}\frac{{p}^{2}(t)}{{\rho}_{0}{c}^{2}}\hspace{0.17em}\mathrm{d}V\approx \frac{1}{2}\frac{{\mathbf{p}}^{\mathrm{T}}{\mathbf{M}}_{\mathrm{a}}^{\text{FE}}\mathbf{p}}{{\rho}_{0}^{2}{c}^{2}},$$  (33) 
$${E}_{\mathrm{A},k}(t)={\int}_{{\Omega}_{\mathrm{a}}}\frac{1}{2}{\rho}_{0}{v}^{2}(t)\hspace{0.17em}\mathrm{d}V\approx \frac{1}{2}\frac{{\mathbf{h}}^{\mathrm{T}}(t){\mathbf{K}}_{\mathrm{a}}^{\text{FE}}\mathbf{h}(t)}{{\rho}_{0}^{2}{c}^{2}},$$  (34) 
$$\mathbf{h}(t)={\int}_{0}^{t}\mathbf{p}(\tau )\hspace{0.17em}\mathrm{d}\tau .$$  (35) 
The radiated acoustical power ${P}_{\mathrm{A}}$ is found as the surface integral of the pressure multiplied by the normal particle velocity. The integration surface is conveniently chosen as ${\Gamma}_{\infty}$. The radiated acoustical energy ${W}_{\mathrm{A}}$ is attained by time integration,
In the numerical framework, the surface integral is attained as a weighted sum of the quantities evaluated at quadrature base points. The pressure and normal particle velocity are obtained at the base points using the conservative interpolation anew as shown in Fig. 2. The time integrals of Eqs. (32), (35), and (36) are evaluated by secondorder accurate schemes.
$${W}_{\mathrm{A}}(t)={\int}_{0}^{t}{P}_{\mathrm{A}}(\tau )\hspace{0.17em}\mathrm{d}\tau ={\int}_{0}^{t}{\int}_{{\Gamma}_{\infty}}p(\tau ){v}_{n}(\tau )\hspace{0.17em}\mathrm{d}S\hspace{0.17em}\mathrm{d}\tau .$$  (36) 
The energy balance is written by summing the energy stored and dissipated in the subsystems as
$$\Sigma E=\underset{\text{Mallet}}{\underbrace{{E}_{\mathrm{M},p}(t)+{E}_{\mathrm{M},k}(t)}}+\underset{\text{Bar}\hspace{1em}\text{with}\hspace{1em}\text{losses}}{\underbrace{{E}_{\mathrm{B},p}(t)+{E}_{\mathrm{B},k}(t)+{W}_{\mathrm{B}}(t)}}+\underset{\text{Sound}\hspace{1em}\text{field}\hspace{1em}\text{and}\hspace{1em}\text{radiation}}{\underbrace{{E}_{\mathrm{A},p}(t)+{E}_{\mathrm{A},k}(t)+{W}_{\mathrm{A}}(t)}}=\text{constant}.$$  (37) 
Figure 5 displays the energy relations obtained in the simulation of the coupled mallet–sound bar–resonator system. Figure 5(a) shows the time histories of the energy stored in the mallet head and sound bar during interaction. At t = 0, only the mallet head has kinetic energy and as it reaches the sound bar, its kinetic energy is transferred to mechanical energy stored in the bar. At the same time, as the head of the mallet is deformed, the mallet head also stores potential energy. When the mallet leaves the sound bar, some of its potential energy is restored as kinetic energy. The mechanical and radiation losses are negligible during the interaction.
The energy relations after the contact, during free vibration and sound radiation, are shown in Figs. 5(b) and 5(c). The only difference between the two cases is the height of the resonator ${h}_{\text{res}}$, which naturally affects its tuning. Figure 5(b) shows the case of an untuned resonator (${h}_{\text{res}}=54.0\hspace{0.17em}\text{mm},{f}_{\mathrm{a},1}=330\hspace{0.17em}\text{Hz}$). The energy transferred from the mallet head to the sound bar during the interaction is slowly dissipated by the viscous losses of the sound bar while a roughly ten times smaller portion of it is radiated into the far field. There is no significant energy stored in the sound field.
If the resonator is well tuned (${h}_{\text{res}}=47.7\hspace{0.17em}\text{mm},{f}_{\mathrm{a},1}=349\hspace{0.17em}\text{Hz}$), the energy relations during free vibration are significantly different compared to the previous case: the initial kinetic energy of the mallet is converted quite efficiently into radiated acoustical energy thanks to the welltuned resonator. It is visible that a nonnegligible amount of energy is also stored in the sound field for $t<200\hspace{0.17em}\text{ms}$. By comparing Figs. 5(b) and 5(c), it is also well observable that the energy stored in the sound bar decays much more quickly if the resonator is well tuned.
The example simulation demonstrates that the implemented numerical framework maintains the energy balance [Eq. (37)] of the discretized quantities with a good precision. The small deviations from the constant in $\Sigma E$, visible in Fig. 5(c) for $t<100\hspace{0.17em}\text{ms}$ is due to the approximation error of the numerical derivatives required for the evaluation of Eq. (36). This example also highlights that the energy relations are significantly affected by the tuning of the resonator. Therefore, the effect of the resonator is investigated further in Sec. IV B.
B. Effects of tuning the resonator
The effect of tuning the resonator is examined by varying the depth of the cavity in the range of 39.4–58.9 mm. For each of the 41 setups, a new acoustical FE model is created and the acoustical system matrices and modes and the vibroacoustic coupling are computed. The tuning range of the models is found as $0.91\le {f}_{\mathrm{a},1}/{f}_{\mathrm{m},1}\le 1.09$, i.e., a range of roughly ±1.5 semitones from the perfect tuning is covered. A strike by the mallet with an initial velocity of ${v}_{\mathrm{M}}(0)=1\hspace{0.17em}\mathrm{m}/\mathrm{s}$ is simulated for a duration of ${\mathrm{T}}_{\text{sim}}=10\hspace{0.17em}\mathrm{s}$. The time history of the radiated sound pressure is exported for a virtual microphone location 0.2 m above the center of the sound bar (see Fig. 4). The running rootmeansquare (rms) value of the sound pressure signal is evaluated using the exponential averaging with a time constant of 15 ms and then the maximal level ${L}_{\text{max}}$ and the −60 dB decay time ${\mathrm{T}}_{60\hspace{0.17em}\text{dB}}$ are determined.
Figure 6(a) shows the resulting amplitudes normalized by the maximal amplitude. As seen, the amplification by a welltuned resonator can be higher than 10 dB. For the sake of comparison, the data points of different measurements by Bork (see Figs. 12 and 13 of Ref. 3) are also shown in Fig. 6(a) with the normalizing of each curve by the maximal amplitude. Whereas these measurements were performed with tubular resonators and different types of sound bars, a fairly good agreement with the results of the FEM are observed. In the experiments, the widest tuning range was measured for the metal bar and, in this case, the mean difference of normalized amplifications between the FEM and experiment is 1.2 dB.
Figure 6(b) displays the decay times. Without acoustical coupling, the decay time of the mechanical vibrations of the sound bar was found to be 3.03 s. The minimum decay time of 0.67 s is found with the perfectly tuned resonator. As the resonator is tuned to the fundamental frequency, its tuning mostly affects the amplitude and decay of the fundamental component. Because the internal losses of the sound bar are greater at higher frequencies [see Eq. (16)], the decay times of higher partials are much smaller than those of the fundamentals. The observed tendencies are also in agreement with those reported by Bork.^{3} Finally, Fig. 6(c) depicts the radiated sound energy ${W}_{\mathrm{A}}$ and the energy dissipated by the internal losses of the sound bar ${W}_{\mathrm{B}}$. As is visible, the dissipated energy increases smoothly at the expense of the radiated sound energy as the resonator is detuned from the perfect tuning.
Mm. 1 contains consecutive samples of the radiated sound with the tuning of the resonator gradually in the range ${f}_{\mathrm{a},1}/{f}_{\mathrm{m},1}=0.96,\dots ,1.00$ in four equal steps.
C. Excitation of torsional modes
Beside the vertical modes being mainly responsible for sound radiation, the 3D model of the sound bar allows for examining the role of torsional modes also. Whereas these modes do not radiate sound efficiently, their excitation may result in a sound of poor quality as they can take away a significant amount of energy. The modal approach offers a straightforward way for analyzing the excitation of the torsional vibrations of the sound bar.
Changing the position or strength of the strike by the mallet both can be expected to have a huge influence on the radiated sound. Due to the nonlinear force–compression relation [Eq. (10)], changing the initial velocity of the mallet head affects both the interaction force and duration of the contact. This effect is visualized in Fig. 7 where the position of the beat is near the end of the bar (point $\mathsf{A}$ in Fig. 9). The results of the FEM are compared to the analytical formulas of Chaigne and Doutaut^{8} and a reasonable agreement is observed with the maximal differences being $\approx 20\%$. Note that the analytical interaction model does not take the vibrations of the sound bar during the contact into account, therefore, some deviations are expected. It was found that the results shown in Fig. 7 are not sensitive to the material model used for modeling the sound bar.
Figure 8 shows the efficiency of the excitation of different modes with the varying of the position of the strike by the mallet. Each colored point on the sound bar represents one of the 684 positions of the strike (i.e., the position of the strike travels through each node of the top surface) and 4 cases with different initial velocities are plotted over the 4 quadrants of the sound bar. Acoustical coupling is not considered in this case and, hence, the arrangement is perfectly symmetric to the two planes indicated by the dashed lines. The energy relations at the beginning of the free vibration phase are shown as ${E}_{\mathrm{B}}={E}_{\mathrm{B},p}+{E}_{\mathrm{B},k}$ and ${E}_{\text{tors}}$ is the sum of the mechanical energy stored in the torsional modes of the sound bar. The top diagram depicts ${E}_{\mathrm{B}}$ normalized by the initial kinetic energy of the mallet head ${E}_{\mathrm{M}}$, whereas the bottom plot displays the ratio ${E}_{\text{tors}}/{E}_{\mathrm{B}}$.
With an increase in the initial velocity of the mallet head, the interaction time is reduced and modes having higher frequencies become also efficiently excited. It is seen in the top diagram that the mallet head can transfer its kinetic energy quite efficiently to the sound bar, except near the nodal lines of the first vertical mode of the bar. Near the latter locations, torsional modes can quite efficiently be excited, especially by stronger hits by the mallet. As the initial velocity increases, the area where the torsional modes are excited becomes larger.
Figure 9 depicts the amount of mechanical energy stored in the first four vertical and torsional modes when striking the sound bar at different positions with an initial velocity of ${v}_{\mathrm{M}}(0)=1\hspace{0.17em}\mathrm{m}/\mathrm{s}$. Positions $\mathsf{A}$ and $\mathsf{B}$ are on the longitudinal axis of the bar and torsional modes are not excited at all by hitting the bar at these locations. Points $\mathsf{B}$ and $\mathsf{C}$ are near the nodal line of the first vertical mode and, hence, the second vertical mode is stronger than the first. In the case of position $\mathsf{C}$, the first torsional mode becomes the strongest. Finally, point $\mathsf{D}$ is at half the length of the bar where the second and fourth vertical modes and the first and third torsional modes have nodal lines, thus, the first and third vertical modes have the most energy.
Sound samples, including acoustical coupling with a resonator of depth ${h}_{\text{res}}=49.2\hspace{0.17em}\text{mm}$ are exported. In Mm. 2, the radiated sound of four samples of striking the sound bar in positions $\mathsf{A},\mathsf{B},\mathsf{C},$ and $\mathsf{D}$ (in this order) shown in Fig. 9 with an initial velocity of 1 m/s are presented. Mm. 3 contains four consecutive sound samples of the radiated sound of striking the sound bar at position $\mathsf{A}$ with the four initial velocities shown in Fig. 8 in increasing order. To highlight the differences in the timbre, the latter four samples are normalized to the same $\text{dBA}$ level.
A coupled vibroacoustic model of the sound generation by mallet percussion instruments was discussed in this paper. The proposed approach relies on the 3D FEM extended by infinite elements for treating the unbounded acoustical domain. The mallet–sound bar interaction model applied is the one proposed by Chaigne and Doutaut.^{8} The modal description was shown to efficiently reduce the size of the system matrices with a few tens of modes providing a suitable representation of both the mechanical vibrations and acoustical field. The orthotropic material model and mechanical and acoustical FEMs were validated by comparisons with measurements performed on an Orff xylophone. The detuning effect of the sound bars covering the resonators was found to be significant and was well reproduced by the FEM.
Simulation of the twoway interactions with the resonator enabled investigating the energy balance of the system. The time stepping solution was found to keep the total energy of the discretized system at a constant level. Effects of tuning the resonator were examined and the numerical model reproduced the experimental results of Bork^{3} both qualitatively and quantitatively quite well with regard to the amplifications and the decay times of the radiated sound. Simulating these effects necessitates the twoway coupling between the sound bar and sound field, which to the best of the authors' knowledge is not included in earlier models. Finally, changing the strength and position of the strike by the mallet allows for inspection of the excitation of the torsional modes of the sound bar, which is facilitated by the applied 3D mechanical model. The interaction times and contact forces produced by the FEM were in reasonable agreement with the analytical model of Chaigne and Doutaut.^{8}
The computational performance of the implemented simulation framework on an average desktop computer is as follows. The most timeconsuming step is the evaluation of acoustical eigenmodes [Eq. (17)], which took $\approx 40$–60 min for the presented simulations, depending on the size of the system. Other preprocessing steps, such as the mesh generation, system matrix assembly, computation of mechanical modes, and the conservative interpolation, took $\approx 15$ min. These steps have to be performed only once for each geometrical arrangement. Then, the time domain simulation of 10 s of sound took $\approx 3$ min. Despite the high sampling rate needed for simulating the mallet–sound bar interaction, the simulation of the short period of the contact took only a few seconds. It is noted that no extra effort has been dedicated to the optimization of the performance of the program codes.
The capabilities of the proposed coupled model can be exploited for investigating other interesting phenomena related to mallet percussion instruments such as the acoustical coupling of neighboring cavity resonators, which was also observed in some of the measurements carried out on the Orff xylophone. The simulation framework may also be applied for creating virtual prototypes of new instrument designs.
ACKNOWLEDGMENTS
P.R. gratefully acknowledges the support of the Bolyai János research grant provided by the Hungarian Academy of Sciences. This work was supported by the ÚNKP–20–5 new national excellence program of the Ministry for Innovation and Technology.
REFERENCES
 1. N. H. Fletcher and T. D. Rossing, The Physics of Musical Instruments, second ed. ( Springer, New York, 1998), pp. 623–648. Google ScholarCrossref
 2. I. Bork, “ Zur Abstimmung Und Kopplung Von Schwingenden Stäben Und Hohlraumresonatoren,” (“On the tuning and coupling of vibrating bars and cavity resonators”) Ph.D. thesis, Technischen Universität CaroloWilhelmina., Braunschweig, 1983. Google Scholar
 3. I. Bork, “ Practical tuning of xylophone bars and resonators,” Appl. Acoust. 46, 103–127 (1995). https://doi.org/10.1016/0003682X(95)93953F, Google ScholarCrossref
 4. I. Bork, “ Measuring the acoustical properties of mallets,” Appl. Acoust. 30, 207–218 (1990). https://doi.org/10.1016/0003682X(90)90044U, Google ScholarCrossref
 5. F. OrduñaBustamante, “ Nonuniform beams with harmonically related overtones for use in percussion instruments,” J. Acoust. Soc. Am. 90(6), 2935–2941 (1991).. https://doi.org/10.1121/1.401768, Google ScholarScitation, ISI
 6. J. Petrolito and K. A. Legge, “ Optimal undercuts for the tuning of percussive beams,” J. Acoust. Soc. Am. 102(4), 2432–2437 (1997). https://doi.org/10.1121/1.419625, Google ScholarScitation, ISI
 7. D. Beaton and G. Scavone, “ Optimization of marimba bar geometry by 3D finite element analysis,” in ISMA2019 International Symposium on Music Acoustics, edited by M. Kob ( Deutsche Gesellschaft für Akustik e.V. (DEGA), Detmold, Germany, 2019), pp. 402–407. Google Scholar
 8. A. Chaigne and V. Doutaut, “ Numerical simulations of xylophones. I. Timedomain modeling of the vibrating bars,” J. Acoust. Soc. Am. 101(1), 539–557 (1997). https://doi.org/10.1121/1.418117, Google ScholarScitation, ISI
 9. V. Doutaut, D. Matignon, and A. Chaigne, “ Numerical simulations of xylophones. II. Timedomain modeling of the resonator and of the radiated sound pressure,” J. Acoust. Soc. Am. 104(3), 1633–1647 (1998). https://doi.org/10.1121/1.424376, Google ScholarScitation, ISI
 10. L. L. Henrique and J. Antunes, “ Optimal design and physical modelling of mallet percussion instruments,” Acta Acust. Acust. 89(6), 948–963 (2003). Google Scholar
 11. I. Bork, A. Chaigne, L.C. Trebuchet, M. Kosfelder, and D. Pillot, “ Comparison between modal analysis and finite element nodelling of a marimba bar,” Acta Acust. Acust. 85(2), 258–266 (1999). Google Scholar
 12. S. Bilbao, A. Torin, and V. Chatziioannou, “ Numerical modeling of collisions in musical instruments,” Acta Acust. Acust. 101(1), 155–173 (2015). https://doi.org/10.3813/AAA.918813, Google ScholarCrossref
 13. F. Irgens, Continuum Mechanics ( Springer, Berlin, 2008), Chap. 7.2 and 9.2. Google Scholar
 14. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: For Solid and Structural Mechanics, sixth ed. ( Elsevier Butterworth–Heinemann, Oxford, 2005). Google Scholar
 15. R. J. Astley, G. J. Macaulay, J.P. Coyette, and L. Cremers, “ Threedimensional waveenvelope elements of variable order for acoustic radiation and scattering: Part I. Formulation in the frequency domain,” J. Acoust. Soc. Am. 103, 49–63 (1998). https://doi.org/10.1121/1.421106, Google ScholarScitation, ISI
 16. R. J. Astley, J.P. Coyette, and L. Cremers, “ Threedimensional waveenvelope elements of variable order for acoustic radiation and scattering: Part II. Formulation in the time domain,” J. Acoust. Soc. Am. 103, 64–72 (1998). https://doi.org/10.1121/1.421107, Google ScholarScitation, ISI
 17. M. Kaltenbacher, M. Escobar, S. Becker, and I. Ali, “ Computational aeroacoustics based on Lighthill's acoustic analogy,” in Computational Acoustics of Noise Propagation in Fluids—Finite and Boundary Element Methods, edited by S. Marburg and B. Nolte ( Springer, Berlin Heidelberg, 2008), Chap. 4, pp. 115–142. Google ScholarCrossref
 18. F. Tisseur and K. Meerbergen, “ The quadratic eigenvalue problem,” SIAM Rev. 43(2), 235–286 (2001). https://doi.org/10.1137/S0036144500381988, Google ScholarCrossref
 19. S. Rubin, “ Improved componentmode representation for structural dynamic analysis,” AIAA J. 13(8), 995–1006 (1975). https://doi.org/10.2514/3.60497, Google ScholarCrossref
 20. N. M. Newmark, “ A method of computation for structural dynamics,” ASCE J. Eng. Mech. Div. 85, 67–94 (1959). https://doi.org/10.1061/JMCEA3.0000098, Google ScholarCrossref
 21. R. Lehoucq and D. Sorenssen, “ Deflation techniques for an implicitly restarted Arnoldi iteration,” SIAM J. Matrix Anal. Appl. 17(4), 789–821 (1996). https://doi.org/10.1137/S0895479895281484, Google ScholarCrossref
 22. MathWorks Inc., “ Online documentation of the builtin function eigs,” available at https://de.mathworks.com/help/matlab/ref/eigs.html (Last viewed 2/9/2021). Google Scholar
 23. D. E. Kretschmann, “ Mechanical properties of wood,” in Wood Handbook: Wood as an Engineering Material ( USDA Forest Products Laboratory, Madison, WI, 2010), Chap. 5, pp. 1–46, available at https://www.fpl.fs.fed.us/documnts/fplgtr/fpl_gtr190.pdf (Last viewed 3/16/2021). Google Scholar
 24. A. Norris and I. C. Sheng, “ Acoustic radiation from a circular pipe with an infinite flange,” J. Sound Vib. 135(1), 85–93 (1989). https://doi.org/10.1016/0022460X(89)907566, Google ScholarCrossref
 25. C. Geuzaine and J.F. Remacle, “ Gmsh: A threedimensional finite element mesh generator with builtin pre and postprocessing facilities,” Int. J. Numer. Methods Eng. 79(11), 1309–1331 (2009). https://doi.org/10.1002/nme.2579, Google ScholarCrossref
 26. See supplementary material at https://www.scitation.org/doi/suppl/10.1121/10.0004216 for additional illustrations, a sensitivity study of material parameters, and the validation of the time stepping scheme. Google Scholar
Please Note: The number of views represents the full text views from December 2016 to date. Article views prior to December 2016 are not included.