Physical modelling techniques for the dynamical characterization and sound synthesis of historical bells

Capable of maintaining characteristics practically intact over the centuries, bells are musical instruments able to provide important and unique data for the study of musicology and archaeology essential to understand past manufacturing and tuning techniques. In this research we present a multidisciplinary approach based on both direct and reverse engineering processes for the dynamical characterization and sound synthesis of historical bells which proven particularly useful to extract and preserve important information for Cultural Heritage. It allows the assessment of the bell’s 3D morphology, sound properties and casting and tuning techniques over time. The accuracy and usefulness of the developed techniques are illustrated for three historical bells, including the oldest recognized bell in Portugal, dated 1287, and two eighteenth century bells from the Mafra National Palace carillons (Portugal). The proposed approach combines non-invasive up-to-date imaging technology with modelling and computational techniques from vibration analysis, and can be summarized in the following steps: (1) For the diagnosis of existing bells, a precise assessment of the bell geometry is achieved through 3D scanning technologies, used for the field measurement and reconstruction of a 3D geometry model of each bell; (2) To access the modal properties of the bells, for any given (at the design stage) or measured geometry, a finite element model is built to compute the significant frequencies of the bell partials, and the corresponding modal masses and modeshapes. In the case of existing bells, comparison of the computed modes with those obtained from vibrational data, through experimental modal identification, enables the validation (or otherwise correction) of the finite element model; (3) Using the computed or experimentally identified modes, time-domain dynamical responses can be synthesized for any conceivable bell, providing realistic sounds for any given clapper and impact location. Although this study primarily aimed to better understand the morphology and sounds of historical bells to inform their conservation/preservation, this technique can be also applied to modern instruments, either existing or at design stages. To a larger extent, it presents strong potential for applications in the bell industry, namely for restoration and re-tuning, as well as in virtual museology.


Introduction
The existence of bells dates back to the Bronze age, c.a. 2nd millenium BC [1], and their morphology has taken many forms over time, depending on aspects such as culture, function, knowledge and technology. Intrinsically related to the history, arts and science of every culture, bells carry important and unique value of cultural heritage. Through the creative efforts of gifted crafters, primitive bells gradually became capable musical objects fulfilling cultural criteria for music performance. If not physically altered, bells maintain their characteristics over time. Their morphology and tuning can be taken as representative of where the instrument was built or located, providing a testimony of the casting and tuning tecniques of the time they were cast.
Apart from shape or size, all bells follow the same operation and physical principles. Their sound results from the vibrational response to a clapper excitation and depends essentially on three main aspects: (1) the bell 3D geometry, (2) the bell material and (3) the excitation conditions (e.g. clapper mass and geometry, impact velocity and contact location). These principles have captured the attention of bell-founders and scientists for centuries. However, only from the 1980s, with the significant growth of technology and computational power, significant steps have been made on physical modelling and simulations of musical instruments' dynamical behaviour. Using computational and experimental techniques, essential information has been obtained regarding vibrational modes of bells [2], the influence of the bell geometry on the bell sound [3], the analysis of bell material [4][5][6], dynamical responses of bells to an excitation [7,8], as well as the development of physics-based sound synthesis techniques [9,10]. Several studies have included 3D analysis of bell vibrations [11][12][13][14]. However, as far as we know, previous research on bell sound and their vibrational features has not yet included precise 3D geometrical data, that can only be obtained through cuttingedge 3D imaging technologies, such as structured-light scanning, a technique based on distortion of known projected light patterns.
In this paper we propose a general physical modelling approach, which takes into account all the underlying physical principles of the bell functioning, in order to characterize its dynamical behaviour and generate realistic synthesized sounds for virtually any bell excited by a clapper. A significant aspect of this research is the use of 3D imaging technologies, namely structured-light scanning, for obtaining the bell geometry, which combined with finite elements computations [15] and modal methods, allows a highly accurate approximation of the dynamical behaviour of these musical instruments. Ultimately, this approach allows a morphological and dynamical characterization of bells, thus contributing to the preservation of important information of the bell sound properties.
The techniques used throughout this work are presented in detail, and their usefulness and effectiveness are illustrated by results from several case studies on historical bells. After providing an overview of the research workflow ("Research workflow" section), the three main stages of this approach are presented in a thorough manner. Section "Geometry assessment" presents the techniques for the measurement of the bell 3D geometry. After describing the used 3D imaging technology, we compare two metrology methods used to infer the full 3D geometry of a thirteenth century bell. The "Modal analysis" section presents the developed modal computation techniques for assessing the modal parameters of bells. Two different modal analysis approaches were used and their effectiveness is illustrated for a eighteenth century bell from the Witlockx carillon at the Mafra National Palace (MNP). Finally, the "Sound synthesis" section addresses the dynamical system formulation and timedomain numerical simulations, in order to perform the sound synthesis of a bell excited by a clapper. To validate our approach, we compare the simulated and measured bell responses for two differently sized bells from the Witlockx carillon at the MNP (see "Illustrative results" section). Additionaly, several simulations on these bells are used to illustrate the power and potentialities of the proposed approach.

Research aims
The main aim of this research is to develop a multidisciplinary approach, which combines up-to-date imaging technology with physical modelling techniques for the dynamical characterization and sound synthesis of historical bells. This methodology allows access to the bell morphology and the significant modal parameters of bells, thus providing and preserving information on sound qualities such as the bell tuning, the effects of asymmetry and ornamentation on the warble [16] and the partials decay times.
Finally, the developed sound synthesis approach allows a set of realistic and comprehensive simulations of the vibrational responses of bells of any size and shape to be performed, either existing or at the stage of design. The developed methodology has potential to be applied in many fields, including bell design, bell restoration, re-tuning, virtual museology and Cultural Heritage preservation.

Research workflow
A flowchart of the global research strategy applied on this study is presented in Fig. 1 and is based on three stages: (1) geometry assessment, (2) modal analysis and (3) sound synthesis.

Stage 1-geometry assessment
A precise assessment of a bell profile is a challenging issue even for modern bell-founders. The usage of calipers, plaster moulds, or laser beams are some of the currently used methodologies for that purpose. In the specific case of carillon bells, given their variable sizes and their frequently hard-to-access locations, these techniques can be painstaking, unreliable, and in some cases, unfeasible. Moreover, because of unpredictable aspects inherent to the bell casting process, the bell geometry varies not only along its height but also slightly along its radius. These radius variations translate in symmetry breaks with important consequences for the bell sound. If perfectly axi-symmetric, bells normal modes would occur in degenerate pairs with identical frequencies and possibly distorted modeshapes. However, because of the mentioned symmetry imperfections, modal pairs become non-degenerate with slightly different frequencies. This causes beats in the radiated sound, a phenomenon known as warble among campanologists [16]. For these reasons, in order to achieve a realistic physical model of the structure, an accurate and detailed assessment of the full 3D geometry of the bell in the order of tenths of millimetre is essential.
Today's 3D imaging technologies combine refined computational methods of in-depth image acquisition with processing algorithms, providing fast, precise and versatile tools with numerous advantages for the assessment of complex geometries. In this work we use a high resolution, portable 3D scanner to assess the full 3D geometry of bells. Because many musically important aspects of the bell sound are related to fine details of the bell geometry, this tool can provide comprehensive ways to characterize and describe historical bells, as well as giving information on casting and tuning techniques.

Stage 2-modal analysis
In a second part of the work, the objective is to assess the modal parameters relevant for the dynamical characterization and sound synthesis of bells (i.e. modal frequencies, modal damping values and modeshapes). Despite the existence of a wide literature on bells, there is still a lack of research addressing the full assessment of these parameters, notwithstanding their significance. The identification of the ratios between the first five partials, at least, is of crucial importance for asserting the bell tuning quality, as bell-founders commonly agree that they should fall into the specific ratios of 0.5 (hum-note), 1 (prime), 1.2 (tierce), 1.5 (fifth), 2 (nominal) [17,18]. On the other hand, as mentioned above, an accurate characterization of close modal pairs can contribute to a more efficient control of the warble. Moreover, the knowledge of modal damping provides important information on the energy dissipation processes, allowing an objective quantification of the partials' decay rate [19]. Finally, through the modeshapes, we can examine the spatial vibration patterns associated with each modal frequency, thus assessing relevant spatial information on the nature of the respective modes, which can be particularly useful to correct the tuning of bells through local structural modifications. In practice, these parameters can be asserted either experimentally or numerically. In this work both approaches are considered and can be used individually or combined, depending on the available information.

Numerical approach
Numerically, the bells modal parameters are assessed, in this work, through the finite elements method (FEM). The basic idea of this technique is to obtain a numerical model that approximates the dynamical behaviour of the structure by using the fundamental equations of dynamics. The continuous structure is discretized into a mesh constituted by a finite number of elements with mass and stiffness properties, which are interconnected at points called nodes. Ultimately, by taking into consideration the continuity relations between neighbouring elements and the conditions of the structure fixation (boundary conditions), the modal parameters of the structure can be computed by solving an eigen-problem. This method has been widely used in the field of campanology for analysing bell modes [15], for exploring new bell designs, including the major third bells [20], as well as for several other purposes on the bellfounding industry. However, apart from recent studies [21,22], reverse engineering techniques, combining 3D imaging technologies with FEM for the study of bells are still scarcely existent in the literature. In this work we explore this powerful combination to achieve a virtual model as close as possible to the real bell. This constitutes a more accurate and efficient approach, which incorporates several steps ranging from the bell geometry (using 3D scans) to the model parameters computations (using FEM).
A clear advantage of using FEM for the modal computations comes from the fact that one can easily assess the modal frequencies and the 3D modeshapes of the structure. However, despite quite straightforward and capable of providing very satisfying results, this method still has some limitations. From the bell 3D geometry, its frequency ratios can be understood by assuming a general bell bronze. However, for an accurate assessment of the modal frequencies values, a precise knowledge of the material properties (mass density, Young's modulus and Poisson coefficient) is needed. Also, the computation of the modal damping values is still a challenging issue. These limitations can be overcome, in the case of existent bells, through the use of experimental techniques.

Experimental approach
The musically important modal parameters (modal frequencies, modal damping values and modeshapes) can also be assessed experimentally, from vibrational measurements. To that end, several experimental approaches were proposed, which can be divided into frequencydomain or time-domain identification techniques [23]. In this work we use a recently developed experimental technique [24], based on a robust time-domain modal identification algorithm, the Eigensystem Realization Algorithm (ERA) [25], which has proven suitable for assessing the modal parameters of complex structures with close frequencies. The basic idea is to adjust some parameters of a mathematical model of the bell response (in this case a time-domain impulsive response) to the vibrational measurements obtained experimentally in real bells. From the ones that best suit the model we can identify the bell modal parameters, including the modal damping values.
Despite the referred advantages, the experimental approach cannot apply for the case of non-existent or non-functional bells. Besides, performing a full vibratory assessment of a bell is a painstaking process, which implies considerable work in practice. A detailed presentation of both numerical and experimental methods used in this work is given in the "Modal analysis" section.

Stage 3-sound synthesis
During the last decades, research in the field of music acoustics has provided a set of techniques for simulating the real behavior of musical instruments. Generally speaking, sound synthesis methods can be divided into two groups: abstract and physical [26]. The fundamental difference between these two categories is the absence (abstract synthesis) or existence (physical synthesis) of an underlying modelling of the instrument physical behaviour. Despite higher computational complexity, physics-based sound synthesis methods allow to perform physically meaningful computations by simple changes in the model parameters.
Sound synthesis of musical instruments has been investigated extensively and significant efforts have been made to simulate their dynamic responses to an excitation, namely for wind instruments [27], string instruments [28] and percussion instruments [29]. A thorough treatment for the specific case of bells can be found in Roozen-Kroon [7].
In this work, based on the computed and experimentally identified modal parameters, a physics-based synthesis technique capable of simulating the time-domain dynamical responses of a bell to an excitation is developed. The idea is to obtain a physical representation of the musical instrument through a set of partial differential equations. In the case of a bell excited by a clapper, this includes a formulation of the bell's structure dynamics, clapper dynamics and the bell/clapper interaction force. Finally, from the physical modelling, a time-step integration algorithm is used to calculate the temporal responses of the bell to a clapper excitation. A particular feature of the bell sound is that the perceived notes, although resulting from actual modal frequencies, often arise as virtual pitches. This phenomenon, known as strike note [30], relates to physical and subjective aspects which are still a matter of debate among campanologists. Although these aspects are beyond the scope of this study, it is worth noting that the strike note is a psychoacoustic effect that results from the behaviour of the auditory system when reacting to the physical modes, and as such, it is intrinsically taken into account in this work. As a result, realistic sounds for any conceivable bell and clapper excitation conditions can be generated, including from archeological artefacts. A detailed description of this approach is presented in the "Sound synthesis" section.

Methodology
As previously mentioned, in order to obtain a precise assessment of the full 3D nearly axi-symmetrical geometry of a bell, we use 3D scanning techniques. In this work we opted for structured-light scanning, a non-invasive contactless method that infers an unknown 3D geometry of a surface based on the distortion of known projected light patterns [31]. This technique is rapid and allows to assess bells 3D geometry in situ, independently of their sizes even in hard-to-access locations as it is often the case in bell towers. We use a portable and high resolution scanner (model Artec Eva), equipped with a structured light pattern generator and two cameras, capable of capturing and simultaneously processing up to two million geometrical points per second, with an accuracy of up to 0.5 mm. The bell is scanned in several parts capturing 3D frames in a point-and-shoot manner (see Fig. 2a). The device is swept along the bell's inner and outer surfaces and the data is stored in the form of sets of point clouds containing the geometrical information. Figure 2a shows the scanning of a thirteenth century bell with the Artec Eva scanner and Fig. 2b shows the real time monitoring of the scanning process. The scanning is followed by a post-processing of the gathered data, using a dedicated software (Artec Studio 12). During this process the separate scans are cleaned and aligned to obtain the full 3D geometry of the bell. Finally, the bell 3D geometrical information can be exported as a point cloud or a surface polygons mesh to a CAD file. Applicable in situ to objects from a couple of centimetres to several meters, this is a versatile, rapid and non-invasive technique which as proven to be effective for scanning bells and has been tested successfully by the authors in different scenarios as presented throughout this paper.

Illustrative results
To illustrate the effectiveness of the scanning technique, we present results from a case study on a thirteenth century bell (see Fig. 2). Dated 1287, this bell was found during an archaeological excavation next to the church of S. Pedro de Coruche, Portugal, and is considered the oldest known bell in Portugal. Figure 3a and b show the photo from the bell from S. Pedro de Coruche, first as it was found in pieces and second after cleaning and restoring the bell. The bell parts were assembled by means of an hypoxic glue. The bell has a diameter of 0.22 m at its mouth and 0.18 m height to its shoulder, with a total weight of 5.6 kg.
Given its delicate condition, the bell was scanned at the museum in Coruche and the post-processing was performed at the laboratory in a total time of one day. The scanning results are shown in Fig. 3c.
On the left, we can see the colored rendered scan of the outer and inner surfaces of the bell. Once the objects are inevitably subject to transformations over time due to corrosion or even sometimes due to accidents, the storage of the textured 3D geometry of the bell can be of  great relevance for Cultural Heritage preservation purposes. At the center we can see the same geometry without color. This representation, together with the colored one, can provide a complementary view to analyse some details of the bell geometry, highlighting for instance geometry imperfections and blisters from the corrosion of the bell bronze [32]. Finally, on the right, we can see a simplified polygon's mesh of the bell geometry, which can be readily imported to a FEM software.
For comparative reasons we present below the results from a previous work by Debut et al. [9], in which a different approach was used for assessing the 3D geometry of the same bell. Without the possibility of using a scanner at the time, a series of geometrical measurements were made, including diameters and thicknesses of the bell in a total of 936 points. A mesh with 36 points regularly spaced azimuthally at 26 height locations was defined and the bell was placed upside-down in a rotating support. At every ten degrees, absolute position measurements were performed through the usage of a singlepoint optical displacement transducer and a slide gauge was used to assess the thickness of the bell wall at each point (see Fig. 4). The 3D coordinates of each point of the mesh were subsequently computed from the geometrical measurements.
Because of its complexity, the geometry of the bell crown could not be asserted through this technique. The complete process took approximately 4 weeks (a detailed explanation can be found in Debut et al. [9]). Figure 5 shows the bell geometry assessed through this approach. When comparing with the results from Fig. 3c, the scanner advantages are clear. Its capacity of capturing millions of geometrical reference points results in a much more detailed geometrical representation of the object when comparing with the laser/gauge technique. Moreover, the analysis is much faster.
In Fig. 6a and b, measurements stemming from the scanner (in green) are superimposed with the ones from the laser/gauge technique (in red). The superposition of the profile measurements (Fig. 6a) highlights the differences between the two approaches. We can see that apart from the head, the bell radius from the laser/gauge technique are, in this case, smaller than the ones obtained with the scanner. This verifies both for the outer and the inner surface. In Fig. 6b the same results can be observed from a different perspective. The bell radius of the outer and inner surface of the bell are alternately represented at four different heights. The larger two circles correspond to the bell soundbow (at 10 mm heigth) and the smaller two circles (184 mm height) correspond to the point where the laser could no longer be used, because of limitations due to the angle of the bell profile at the shoulder. As we can see, the green circles (measurements from the scanner) are, consistently, slightly larger than the red ones. The same can be observed from the bell profile.
In Fig. 6c and d, the differences between the two 3D geometries are represented. The geometries obtained through both approaches were superimposed, taking the scanner as the reference and the distances between them were computed. In the figure, we can see the geometry of the outer (Fig. 6c) and inner (Fig. 6d) profile of the bell assessed using the laser/gauge technique. The distances in millimeters to the reference geometry are represented through the colour map. This representation allows to visualize what is inside (in blue) and what is outside (in red) of the volume circumscribed by the reference Fig. 4 Setup for the bell geometry description. a Outer surface measurements b horizontal thickness measurements [9] geometry surface. In white are represented the parts where both meshes nearly coincide. Figure 6c shows that apart from the bell head, in general, the measurements of the outside surface with the laser/gauge technique are inside the volume created by the surface measured using the scanner. Figure 6d shows the inner surface, where apart from the lip and top plate, the bell walls with the laser/gauge technique "come out" from the scanner measurements. Ultimately, this means that the bell represented through the laser/gauge technique was in this case slightly smaller than the one stemming from the scanner. From the histograms we can also understand that, despite a general good agreement between both approaches, significant differences up to 1.5 mm can be found.
In conclusion, a comparison between the scanner and the laser/gauge technique demonstrates an improvement of more than one order of magnitude in dimensional precision and an enhancement of spatial coverage by several orders of magnitude. Additionally, only the scanner allows the assessment of the bell crown, a crucial element for the computation of the full volume of the bell. Apart from morphological information, this is also extremely useful for asserting material properties such as the mass density from the bell. Another advantage is that the geometrical data can be obtained in a fraction of the time using the scanner and the bell geometry can be readily exported in CAD files to input FEM software.

Finite elements modelling
As previously mentioned, through the finite elements method, the modal parameters are computed based on a physical model of the bell. In a first step, the geometry of the bell needs to be inserted into a finite elements package. In our work, we import it as a CAD file with the bell surface geometry stemming from the 3D scanning methodology, described in the previous section. This constitutes a significant aspect of this work, as the accuracy of the physical model is strongly dependent on the accuracy of the geometrical data.
From the surface geometry, a 3D volumetric mesh is generated to create a FEM model. Mass and stiffness properties can then be defined for each element of the mesh, based on bell bronze material properties such as the Young's modulus, the density and the Poisson's ratio. Taking into consideration the continuity relations between neighbouring elements and the conditions of the structure fixation (boundary conditions), a physical model of the bell can then be obtained.
Mathematically, most structural dynamics problems of a linear mechanical system may be described by the finite element equation of motion where M and K are the global inertia and stiffness matrices of the system (built by assembling the elementary mass and stiffness matrices that describe the physical properties of each element of the mesh), and y(t) is the vector of physical displacements at location r [33]. Solution of Eq. (1) can be assumed in the following harmonic form  where ϕ is the vector of modal amplitudes of each node of the mesh in every direction of motion. From this assumption, Eq. (1) becomes the generalized eigenvalue problem with the classic formulation from which the system modal frequencies f n = ω n /2π and respective modeshapes ϕ n can thus be obtained.
It is worth noting that the convergence of the FEM solutions to represent the real bell behaviour is dependent on the mesh parameters, namely the type and number of elements, that need to be sufficiently detailed for representing the structure. The choice of these settings is problem-dependent. Based on preliminary convergence tests, we use in this work solid elements and take profit of the spatial resolution of the scan data for obtaining a highly detailed mesh, thus enhancing the approximation solutions. For results see "Illustrative results" section.

Experimental modal identification Experimental procedure
The bell's responses to an excitation are collected at several points along the bell rim. A mesh of 32 points equally spaced at the bell's mouth is defined and an instrumented Fig. 6 a Superposition of the thirteenth century bell profile measurements assessed through the 3D scanning technique (in green) and the profile obtained through the measurements with the optical displacement transducer and the slide gauge (in red). b Inner and outer surfaces of the bell at four different heights (scanner in green and optical displacement transducer and laser/gauge technique in red). c and d Colour representation of the distances in millimeters between the 3D geometries, assessed through the scanning technique and the profile obtained through the measurements with the optical displacement transducer and the slide gauge impact hammer is used for the excitation at each location. Attention needs to be paid when selecting the impact hammer configurations to guarantee an adequate amount of energy for the frequency range of interest. In our case, an impact hammer (Brüel&Kjaer type 8202), instrumented with a force transducer (Brüel&Kjaer type 8200) and equipped with several tips of different stiffnesses is used for the small and mid-range bells. For heavier bells a larger hammer was designed by the authors, capable of exciting a lower frequency range.
The corresponding vibrational radial responses are measured with three piezo-electric accelerometers (Brüel&Kjaer type 4375), positioned at the asymmetric angular coordinates 0°, 23° and 146° along the bell rim. This minimizes the probability of all three being located at nodal points of the same mode of vibration, thus ensuring that a maximum number of modes is detected and identified. Through this strategy the modeshapes at the bell's mouth can also be assessed, constituting a useful criterion for the modal identification process.
Signals are acquired and converted into digital form using an acquisition board (SigLab/Spectral Dynamics, model 20-42). To comply with the Shannon theorem, a sampling frequency of F s = 51200Hz is used, as well as anti-aliasing filters. Signals are 12 s long, resulting in a frequency resolution f = 0.08Hz , and allowing accurate estimates of the dissipation bell properties, since the decay times are of the same order than the duration of the acquisition.

Identification of the modal parameters
As previously mentioned, a software was developed at the laboratory, based on a polyreference multi-modal identification algorithm, that uses the time domain responses for the modal identification. From the acquired time signals, the transfer functions H ij (ω) =Ẍ i (ω) F j (ω) are computed through the Fourier transform of the excitation force F j (ω) and accelerations Ẍ i (ω) measured at the mesh locations i and j . By computing the inverse Fourier transform of H ij (ω) , the corresponding impulse responses h ij (t) are obtained, from which the modal parameters are finally extracted using the ERA algorithm. Based on a mathematical parametric description of the system dynamical behavior, the algorithm adjusts the modal parameters to best fit the measured signals, assuming a sum of damped modal responses as where ω n = 2π f n and ζ n are the modal frequency and damping values of mode n , and A ij n is the respective modal participation factor related to the mode-shape (4) h ij (t) = N n=1 A ij n e −ω n ζ n t sin(ω n 1 − ζ 2 n t) values of ϕ n at the excitation and response locations. The modal parameters ω n and ζ n are extracted from a set of measured impulse responses organized in the form of a generalized Hankel matrix, while the modeshapes are obtained by adjusting, in the least-squares sense, a set of model responses to the corresponding measured signals. A sensitive issue of the identification process is to determine the order N of the mathematical model due to the polluting noise in the measured signals, which perturbs all identification algorithms. To overcome this problem, the number of modes of the model must be oversized, leading to the identification of some modes which are unphysical. The difficulty is then to select the physically relevant modes from the modes identified from the oversized model. For that purpose, we rely on a stability diagram that allows an evaluation of the modal convergence as the order of the model grows, as well as on the identified modeshapes. A detailed presentation of the developed ERA method applied to bell modal identification is given at Debut et al. [24].

Illustrative results
The results from the numerical and experimental approaches are now presented and compared through the modal analysis of a eighteenth century bell from the Witlockx carillon of the Mafra National Palace (see Fig. 7). The bell has a diameter of 0.2 m at its mouth and 0.16 m height to its shoulder, in a total weight of 8.4 kg.
Here, for the numerical approach, the bell geometry was scanned as described in "Geometry assessment" section and then used as input for the finite elements software. In Fig. 7b, on the left, is presented the polygons surface mesh from the 3D scan, where the bell inner and outer profile can be recognized. At the center, the achievable level of detail is illustrated through the full-rendering model.
The finite elements software Cast3M [34] was used for constructing the bell physical model and for performing the computations. A volumetric mesh consisting of 29,070 solid elements (4-node tetrahedral) was generated (Fig. 7b) on the right. For the material properties, the bronze density was deduced from the measured mass and the computed volume of the bell, and the Young's modulus was adjusted until the first modal frequencies from the FEM computations reproduced the experimentally identified results. Finally, values of ρ = 8660 kg/cm 3 and E = 90 GPa were used for the density and Young's modulus respectively, and a typical bronze value of 0.34 was assumed for the Poisson's ratio. The computations were performed assuming free boundary conditions as if the bell was suspended.
The FEM computation results for the first five lower partials are presented in Table 1. As previously  mentioned, normal modes appear in pairs because of the symmetry properties of the bell structure. One of the FEM advantages is the possibility of computing the full 3D modeshapes. Table 1 shows a general view of the modeshapes of each modal pair on the left (X, Y view) and on the right, the relative orthogonal position of the modal pairs is represented from the azimuth plane relative to the bell mouth (X, Z view). Dark blue represents the parts where the bell movement is null and the gradually warmer colours correspond to areas where the displacement is larger, with red representing the maximum displacement. The computed modal frequencies are also presented in the table. The differences between modal pairs highlight the effects of symmetry imperfections on the bell geometry captured by the scanner. As previously mentioned, this translates into beats in the bell sound at a rate equal to the frequency difference between modal pairs. This way, a quantification of the warble can be asserted, which is a relevant indicator of the musical quality of bells [16]. The internal tuning of an individual bell can also be analyzed from the modal frequencies.
One may notice that the frequency ratios are far from the typical targets of 0.5, 1, 1.2, 1.5 and 2, in particular for the three highest modes. This can be explained by the fact that for smaller bells, these high frequency modes die away very quickly once the bell is struck by the clapper, which makes them very difficult to tune. On the other hand, the higher partials are situated in a frequency range where the tuning differences are not easily perceptible by the ear and, for this reason, their tuning becomes of less importance for the bell sound. The bell modal parameters were also assessed experimentally through the approach presented in "Modal analysis" section. The experimentally identified modal parameters are shown on Table 1. The good agreement of both approaches with relative differences between modal frequencies of less than 0.5% validates the FEM methodology and gives us confidence on the effectiveness of the developed techniques. Through the experimental approach it was also possible to identify the modal damping values, an important parameter, namely to account for the energy dissipation phenomena for the subsequent modeling in Stage 3. The identified modal damping values are also presented in the table.
Finally, through the combination of numerical and experimental approaches, the complete set of modal parameters required for the bell characterization and for the time domain dynamical simulations in Stage 3 were successfully assessed.

Methodology
The proposed sound synthesis approach is divided in two main steps. First, the physical modelling of a bell impacted by a clapper is presented. This includes the formulation of the bell dynamics, clapper dynamics and the non-linear interaction between bell and clapper. In a second step, the time-domain simulations of the bell responses to a given clapper excitation are computed. A detailed description of this approach is presented below and several results stemming from parametrical computations are illustrated.

Bell dynamics
We now turn to the bell dynamics formulation. After a modal discretization of the bell dynamics in terms of its modal properties, its physical motion y(t) at a generic position r = (r, θ, z) can be described as where ϕ n (r) are the modeshapes and q n (t) the respective modal amplitudes in every direction of motion.
The partial differential equation of the bell motion is then replaced by a set of N ordinary second-order modal equations which governs the time-dependent modal amplitudes, written in matrix form as: where M , C and K are diagonal matrices of the modal masses, modal damping and modal stiffnesses. The modal forces F (t) are obtained by projecting the external forces f (t) on the modeshapes ϕ n of the modal basis, as where This modal approach is particularly suitable for the objective of the work, as it provides a physical model with a reduced number of degree-of-freedom, and consequently requires less computational efforts for the sound synthesis.

Clapper dynamics
The clapper is modelled as a rigid body of mass m c , moving perpendicularly to the bell wall, neglecting any elastic deformation after the bell/clapper contact. Its normal (radial) motion z(t) is thus represented by the following dynamic equation where f R (r 0 , t) is the radial component of the applied force.

Bell/clapper interaction
The bell/clapper interaction for the excitation must be then formulated. The Hertz theory of elasticity [35] has shown to be particularly suitable for this purpose, describing the local interaction force f c (t) during the contact time as where z(t) is the normal clapper motion, y R (r 0 , t) is the normal motion of the bell at the localized radial impact excitation at location r = r 0 , b is a contact parameter equal to 3/2 for a sphere in contact with an elastic halfspace, and K c is a contact stiffness coefficient which depends on the geometric and elastic properties of the two solids. In the modal framework, we can write the interaction force as: where ϕ R n (r 0 ) is the radial component of the mode shape ϕ(r 0 ) . Despite some limitations, the nonlinearity of the model (through the power parameter b) allows an adequate reproduction of the significant spectral changes observed in bell sounds, when a bell is struck with different impact forces [36]. Finally, from Eqs. (5)-(7), (9) and (10), the modal and physical responses of a bell to a given clapper excitation can be computed. In a last step, the time-domain responses of the bell to a given clapper excitation are computed and sound files are generated.

Time domain dynamical simulations
The system motion is computed by time-step integrating the modal equations presented in "Sound synthesis" section. For that purpose we use the velocity-Verlet scheme [37], which has proven to provide accurate enough results for vibro-impacting systems [29,38]. Assuming initial conditions for the positions, velocities and accelerations of the clapper and bell at a given time t n , this explicit algorithm yields these quantities at the following time step t n+1 .

Illustrative results
We now present the results from several parametric computations to explore the potentialities of the proposed approach. Simulations on two historical Witlockx bells from the Mafra National Palace are used to illustrate and discuss different musical effects.

Witlockx small bell
Concerning the eighteenth century bell from the Witlockx carillon of the MNP described in "Illustrative results" section, for comparative purposes we present below the bell sound recorded with a microphone and the corresponding synthesised sound stemming from the proposed sound synthesis approach. The modal parameters used for the sound synthesis were assessed by combining both numerical and experimental approaches as described in "Modal analysis" section. For the recorded sound, the bell was struck at the soundbow (0.03 m height from the rim) and the microphone was placed at 1 m distance from the rim. Regarding the simulations, both the geometrical and material properties of bell and clapper were taken into account, as well as the excitation conditions assuming a clapper with mass m c = 0.3 kg and ball radius r c = 0.015 m, and an impact velocity of v 0 = 0.1 m/s at a height h c = 0.03 m from the rim. Numerically, one major difficulty arises from the short duration of the contact time, usually in the range 0.2-0.5 ms, which also causes the excitation of a wide frequency range. After convergence tests, a time step in the order of 4.10 -6 s and a very extensive modal basis in the order of 110 modes were used in this work to obtain realistic simulations. Figures 8a and b show the spectrograms of the measured and synthesised sounds, respectively. On the spectrograms it is clear that a large number of partials are excited when the sound starts, dying away at different rates. The highest frequencies decay very quickly and the lower ones slowly emerge from the bright initial sound until only the hum-note remains. Apart from the noise, which is absent on the simulations for obvious reasons, a good agreement can be verified between the two spectrograms. A close look at the higher frequencies reveals some minor differences, namely on the decay rates and the partials amplitudes. A possible reason for this is the fact that the modal damping values were only identified for a limited number of modes and using extrapolated values for the non identified higher modes. Besides, in practice, it is also very difficult to assert the precise impact location which leads to some differences. Additional file 1 and Additional file 2 correspond to the measured and synthetised sound respectively. The results illustrate the achievable level of realism by the sound synthesis, thus validating our proposed approach. Figure 8c shows the simulated bell/clapper interaction force and the corresponding bell and clapper motions. On the top of the figure we can see the simulated impact force assuming an initial velocity v 0 = 0.1 m/s at height h c = 0.03 m. A short contact time of around 3 ms corroborates with typical real impacts, as observed experimentally. Also, note that the impact force is symmetrical as it should. In the middle, in blue, we can see the displacement of the clapper, penetrating the bell wall in a first moment, until it starts to move in the opposite direction at half the contact time, finally distancing itself from the bell at the end of the contact time. At the bottom, in black, the bell displacement is represented. After the initial contact, the bell wall is pushed by the clapper, both moving with the same orientation until the bell wall starts to move back towards its equilibrium position initiating a vibratory movement. Overall, these results show that the quasi-static contact Hertz model is well-suited for modeling the dynamic bell/clapper interaction.
In Fig. 8d, two parametric computations simulating the response of a bell as function of the impact velocity are presented. Due to the non-linear behaviour of the bell/ clapper interaction, a large impact velocity translates into a short contact time, thus being more effective in exciting higher frequencies, ultimately resulting in a brighter sound. For comparative purposes, acceleration signals were normalized relative to their energy contents.

Witlockx large bell
To illustrate the versatility of the full methodology, this subsection presents the use of the above-mentioned procedures for the case of a large Witlockx bell from the same carillon of the Mafra National Palace. The bell has a diameter of 0.89 m at its mouth and height h = 0.69 m to its shoulder and a total mass of 416 kg.
In this case, as for the small bell, the experimental and numerical techniques presented in "Modal analysis" section could be combined and the complete set of modal parameters could be assessed. It is worth noting that there was a time lapse between the experimental modal identification and the bell metrology. Given the carillon supporting structure deteriorating conditions, at the time the scan was made, the bell was already supported and surrounded by scaffolds for safety reasons. Together with the hammers and cables of the automatic playing system, this left us with very little room to work. Thus, an extra work for scanning and post-processing was required in this case, as all these features and nearby structures needed to be cleaned from the scans.
In Fig. 9a, we can see the bell scan result. Aspects such as text and ornamentation highlight the level of detail achievable through this metrology technique, even in such difficult conditions. In Fig. 9b the effect of the impact location on the bell vibratory responses is illustrated. It is important to notice that when a clapper impacts the bell at the vicinity of a normal mode node, this partial is significantly supressed in the sound of the bell. From left to right, five similar impacts were simulated at the respective heights h c = 0.2 h, 0.4 h, 0.6 h, 0.8 h and 0.95 h, with h = 0.69 m being the bell height. Obvious changes on the spectral responses are shown as the impact location varies from the soundbow to the bell head. The spectrograms show a clear variation in the modal amplitude of the several modes depending on the striking point. This illustrates that the level of excitation of the various partials can be selected by striking the bell at different locations. Depending on the clapper excitation point, its energy is distributed differently through the bell modes, thus leading to sounds with different energy levels along the frequency range. As the excitation moves towards the bell head, in general, higher-order modes become more excited whilst the amplitudes of the lower modes decrease, leading to a gradual change of the timbral characteristics. This is illustrated in Additional file 3, where five sounds can be heard. We start by hearing a well-balanced sound (strike at the soundbow) which progressively becomes brighter as the impact location rises vertically.
Another important aspect for the bell sound is the choice of the clapper material. Figure 9c shows the simulated spectrograms of the bell responses varying only the contact stiffness parameters (K c = 1e10 N/m and K c = 1e11 N/m, respectively). The spectrogram on the right shows that using a harder clapper, the higher frequencies are excited in a more efficient way, resulting in a brighter sound. A softer clapper (spectrogram on the left) will act as a low pass filter, resulting in a "sweeter" sound.
Finally, in Fig. 9d, the computed spectra of the bell accelerations for different clapper masses are presented. On the left, the simulated bell responses to a lighter (in red) and a heavier (in blue) clapper, with 1.38 kg and 13.8 kg respectively, are superimposed. Spectra of the bell velocity at the contact location (normalized to their energy) are presented to emphasize the difference in their frequency contents. As we can see, lighter clappers are able to excite a larger number of modes, whilst heavier clappers are more effective for exciting the lower ones. This can be explained by the increase of contact-time due to a larger mass, which restricts the energy transfer to the higher modes. The sounds simulations considering a large (13.8 kg) and a small (1.38 kg) clapper corresponding to Fig. 9d are illustrated in the Additional file 4 and Additional file 5, respectively.

Conclusions
In this paper we present a general methodology for the dynamical characterization and sound synthesis of bells, comprised of three main stages. In a first stage 3D imaging techniques are used for the geometry assessment. In a second stage, numerical and experimental modal analysis techniques are combined to assess the significant modal parameters of bells. Finally, based on the computed modal parameters, this approach includes a physics-based synthesis technique capable of simulating the time-domain dynamical responses of a bell to an excitation. For each stage, the effectiveness and usefulness of each technique are illustrated through the obtained results for three historical bells.
One of the strengths of this methodology is the use of leading-edge imaging technologies to overcome the challenging issue of obtaining the full 3D geometry of a bell. Given that many musically important aspects are strongly related to fine details of the bell geometry, this information is basic for a finite elements model and leads to a notable improvement on the modal computation results. Besides, the combination of finite elements analysis with experimental modal identification enables the full assessment of the significant modal parameters of bells. It is worth noting that, given the high accuracy of the developed techniques, significant information on the bell musical qualities can be obtained, namely on the bell tuning, the effects of asymmetry and ornamentation on the warble and the partials decay times.
Finally, the developed sound synthesis approach allows to perform a set of realistic and comprehensive simulations of the bell vibrational responses. This includes the formulation of the bell dynamics, clapper dynamics and bell/clapper interaction forces. Despite the satisfactory results, some limitations can still be pointed out, namely to the sound synthesis approach. The used local contact model still remains a crude representation of the real interaction dynamics. However, it proved to be a plausible pragmatic approach for representing the bell/clapper interaction. Another aspect is that the used sound synthesis approach still lacks the delicate aspect of the sound