© The Author(s), Published by EDP Sciences, 2025

1 Introduction

The relations between the sound and playing characteristics of finished violins and the properties of their components (geometry, materials) have remained elusive for as long as they have been scientifically studied, since the pioneering work of Carleen Hutchins and Frederik Saunders in the 1950s which led to the proposal of guidelines for makers [1]. To this end, numerous efforts have been put in the last few decades, on developing finite element (FE) models of violins that simulate the main features of the instrument [2–9]. Parametric models have been also developed to study the influence of the construction parameters [10] or material choice [11] on the vibrational behaviour. Few of these computational studies, however, present a comparison with experiments.

The lack of experimental validation comes from different factors. Violins are usually made by hand, by makers who need to consider not only objective features of the wood, but also aesthetic and economic factors. The nature of the handwork makes geometric reproduction of one instrument a rather difficult task. Moreover, the arching and the variations in thickness all along the plates, which have an impact on the static and dynamic behavior, are not trivial to model using computer aided design (CAD) based FEM. In addition, violins are made in a very variable material (spruce for the top, maple for the back) and assembled with glue which is an organic material prepared by the luthiers in different ways and which adds uncertainty to the manufacturing process, especially when making the instrument body. These factors explain as well why it has been extremely difficult to find relationships between, on the one hand, the construction parameters and the vibratory behavior of the sub-structural components and, on the other hand, the vibro-acoustical behavior and the sound qualities of the finished instrument.

In such context, the Bilbao project [12] aimed at exploring the influence of a single geometrical parameter by controlling in an unprecedented manner the other parameters and the construction process. The chosen parameter was the thickness of the plates, as the heritage of old instruments has all kinds of graduation schemes for which we lack knowledge of how the playing and sound characteristics are affected and how the decisions were made [13]. Six violins were thus built with different thickness patterns: three instruments with normal backs, each paired with a thin, normal, or thick top; similarly, three violins with normal tops, each paired with a thin, normal, or thick back. The two examples of normal top paired with normal back served as a control.

A high level of accuracy was achieved by using, when it was possible, automated routing process with computer numeric control (CNC); by carefully selecting the wood used (from the same trees); and by controlling and monitoring the construction of all the sub-structures and the assembly process in order to minimise the various sources of variations. Despite this careful control, some variability due to the wood properties of all the different parts and the geometrical tolerances was inevitable. Therefore, creating a FE model is proposed to initially aid in investigating the relative influence of this unavoidable variability versus the controlled variations in thickness within the framework of this project. Subsequently, it can be utilized more broadly as a predictive tool for instrument makers. To this end, the FE model should be quite refined and would need to be validated at different stages of the construction. As the second goal of the Bilbao project was to follow the vibratory behavior from the substructures to the assembled instruments, various measurements (including modal analyses) were conducted at different stages of the construction. The large amount of measurements available represent therefore an unprecedented data set for validating the FE model all along the assembly of the sub-structural components, and this for many instruments, while such validation has so far been only done on a single instrument. Indeed, after a pioneering study [14] on a disassembled factory instrument (which geometry could only be known by hand measurements), [15] modelled a violin body that had been specifically made and assembled for the study, and whose components had been 3D scanned and measured in terms of vibratory behavior.

The long-term goal of this project is thus to propose accurate models of the Bilbao violins that are validated through test-model correlation at all steps in order to be able to differentiate uncertainties from deliberate modifications (in thickness like for the real violins but eventually for other geometric parameters too). An outcome would entail establishing thresholds for the amount of geometric modifications necessary to reliably evaluate their impact on real instruments.

Similar to the construction of a violin, developing the full FE model of the six Bilbao violins requires a large number of steps, in order to model progressively all sub-structural components and their assembly. This paper focuses on the modelling of the top and back plates just after they were CNC routed. The batch of plates and the experimental data that will be used for validating the model are presented in Section 2.1 while Section 2.2 presents the methodology used to build the model and to optimise the wood’s material properties. Section 3 provides the resulting optimised properties and illustrate the very good agreement between the simulated and the experimental modes. Finally, the conclusions are presented in Section 4.

2 Materials and methods

2.1 Materials

2.1.1 The batch of 12 plates

As it is very hard to fully CNC cut violin plates with specific thickness patterns, only the outside surface of the plates was actually CNC cut for the Bilbao project. The graduation was then done manually, but performed entirely on the inside surface leaving the outside surface untouched. The batch used in this study corresponds to the 12 plates in their initial stage, just after they were CNC routed, using a computer aided design file that was derived from a medical CT-scan of the Huberman violin, an instrument crafted by Antonio Stradivari in 1713. For practical reasons linked to how violins are built, the thickness of the tops was constant (3.5 mm) while the thickness of the backs varied between 3 and 6 mm.

Wood pieces were matched as closely as possible – for both the tops and the backs independently –, in terms of sound speed (therefore specific rigidity) and density (see Sect. 2.1.3). The six spruce wedges were selected as the closest among the nine horizontally adjacent wedges available from the same tree, while the maple wedges were chosen as vertically adjacent wedges coming from three logs. This led to six very similar top plates (labelled from A to F) and six very similar back plates (labelled similarly).

2.1.2 Experimental modal analysis

Experimental modal analysis was conducted in situ at the workshop of one of the authors (RJR), just after the 12 plates were cut, and before they were graduated. The impulse responses were obtained using an accelerometer and a manual impact hammer. The plates were laid on a thin silicon net so the boundary conditions can be considered as free. The accelerometer (Dytran, 0.2g, series 3224) was fixed, located on the longitudinal middle line, a few centimeters off the vertical central line, towards the tail block (see Fig. 1, left). This location was chosen with the foreknowledge of typical violin plate modes (e.g. [16]) so the accelerometer was not on any nodal lines for at least the first 10 modes. The plates were hit with a miniature impact hammer (PCB, model 086E80), five times on each of the 70 and 43 points spread out evenly on the surface of the top and back plates respectively (see Fig. 1, right). The locations (and their number) were selected as a trade-off between providing a good surface coverage while keeping a reasonable duration for the measurements.

Figure 1 Left: Setup for performing the modal analysis of the free plates. Right: Grid of the points at which the transfer functions were measured, on the top (left, 70 points) and the back (right, 43 points) plates.

The acquisition as well as the modal identification were conducted with an in-house software, entirely coded by one of the authors (GS). The latter was done using the global singular value decomposition technique [17, 18]. The eigenfrequencies and interpolated modal shapes of the first six modes obtained for the top A and back A are provided as examples in Figure 2, while the others are provided in Appendix B (Figs. B1–B10).

Figure 2 Interpolated modal shapes of the first six modes and corresponding eigenfrequencies measured on free top A (a) and back A (b). The displacement along the Z axis (see Fig. 3) is displayed.

Figure 3 Actual mesh of the top (left) and back (right) plates of the Bilbao violins, generated with tetrahedral elements in COMSOL.

Figure 4 Sensitivity analysis matrix for a 1% variation on each of relevant material parameters (density and elastic constants). Results for top C (left) and back C (right). The values are normalized for each mode. A high value is associated with a large influence and the sign indicates the direction in which the eigenfrequencies are affected.

Figure 5 (a) Distribution of simulated (red dots) and experimental points (blue circles), (b) Distribution of the experimental points (blue circles) and their closest simulated neighbours (red circles), highlighting the alignment between them.

Figure 6 MAC matrix for back B before optimisation (a) and after the optimisation (b) using only εfreq for the cost function.

Figure 7 MAC matrix for back B after the optimisation using the cost function C.

Figure 8 Modal shapes of the first six simulated modes and corresponding eigenfrequencies for the free top A (a) and back A (b). The displacement along the Z axis is displayed.

Figure 9 MAC matrices for top A before optimisation (a) (with the initial values of the wood properties described in Sect. 2.2.1); and after optimisation (b).

2.1.3 Measured wood’s mechanical properties

Determining the properties of wood without resorting to destructive testing presents significant challenges [19]. This is particularly true for violin wood, given its orthotropic mechanical properties, which entails variable behaviors across different axes and the evolution of properties in response to hygrometry. Additionally, the material may exhibit local variations in density and rigidities, influenced by factors such as annual rings (for temperate species), fiber orientation, and the arrangement of different cells, predominantly tracheids for spruce and maple (for the latter, vessels and multiseriate rays are also present and have an impact on density and anisotropy) [20–24]. Density was estimated on the wedges, i.e. the quarter-sawn trapezoidal prisms that serve as a basis for cutting the plates. The Bilbao project team (led by two of the authors, UI and GS) employed two measurement methods. The first measurement was geometric and relied on measuring the dimensions of the wedge (and thus estimating the volume) and the weight after the wedge was kept in a climatic chamber (a closed cabinet with UV-A and incandescent lamps, with a controlled moisture and temperature, respectively set at 44% air humidity and 33 °C) to reach an equilibrium moisture content of 8%, according to [25]. This value is considered standard to make string instruments, as the wood has to be worked with a lower value of moisture content than the one that is reached in usual room conditions. The second measurement was done by the so-called water displacement method. Both approaches yielded similar densities and an average value is provided in Table 1 for the 12 wedges, the six spruce (Picea abies) that were used for the top plates and the six maple (Acer pseudoplatanus) that were used for the back plates.

The Bilbao team measured as well the sound velocity with different methods, for the original wedges and then for the CNC cut plates. For the wedges, they only measured the velocity along the wood longitudinal L direction (along the grain), using two methods:

• by measuring the frequency f of the first longitudinal (compression) mode of the wedge and approximating it as a beam of length L so c = 2fL. This was done by hitting on one side with a hammer and measuring on the other side (along the length of the wedge) with a microphone.

• by using a Lucchimeter [26, 27].

The longitudinal Young modulus E_L was then obtained with the following formula:

$$E_L=\rho c^2.$$(1)

For the cut plates, they measured the sound velocities in both the longitudinal L (along the X axis, as defined in Fig. 3) and radial R (across the grain, along the Y axis) directions of the cut plates using the Lucchimeter. However, as such measurements rely on longitudinal (compression) waves, the relationship between the sound velocities and the Young’s moduli along the grain (E_L = E_X) and across the grain (E_R = E_Y) is not straightforward for an arched thin plate and equation (1) is no longer valid. So only the values measured on the wedges are reported in Table 1.

Values derived from Lucchimeter measurements are consistently higher than the ones derived from the frequency of the first longitudinal/compression mode, with a discrepancy going up to almost 3 GPa in some cases, confirming previous observations [28].

2.2 Method

2.2.1 Model description

In this section, the FE element model used for simulating is described. The simulations were conducted using a three-dimensional (3D) model developed with COMSOL Multiphysics® software. The geometry files used for CNC carving were meshed with tetrahedral elements with quadratic interpolation. A number of at least six elements with linear shape function per wavelength is usually considered. In our case, we defined mesh size to satisfy a condition of at least six elements with quadratic shape function (11 nodes) per wavelength. The first six modes that will be used for the optimisation (see Sect. 2.2.3) are below 1000 Hz, leading to a wavelength lower than about 180 mm in the longitudinal direction (L) and about 90 mm in the radial (R) direction. This wavelength calculation is based on the wave speed of flexural waves in the L and R directions of a 3 mm thickness wood plate, with E_L = 13000 MPa, E_R = 1000 MPa and a density ρ = 400 kg · m⁻³. This implies that the elements should be smaller than 15 mm, which is the case. The mesh is comprised of nearly 40000 and 30000 elements for the top and the back respectively and displayed on Figure 3. The exact same mesh was used for the six same plates (tops and backs respectively).

The quality of the elements was evaluated with the skewness metric, which measures the distortion of tetrahedral elements compared to an ideal shape, and is commonly accepted in the literature for assessing the quality of tetrahedral meshes [29]. For the top plates, the minimum element quality (skewness) is 0.1967, with an average element quality of 0.6626. For the back plates, the minimum element quality is 0.1921, with an average element quality of 0.6631. The total volume corresponding to all the 3D elements of the mesh was checked to be equivalent to the volume of the plates (190 cm³ for the top and 213 cm³ for the back) estimated from the CAD files, showing that the plates were properly meshed without holes or incomplete filling.

2.2.2 Wood’s mechanical properties

The material has been modeled assuming linear elastic orthotropic behaviour. While the geometry of the plates was accurately known, the wood properties were largely unknown. They have therefore been estimated by optimising the simulated eigenfrequencies so they match as closely as possible the measured ones. As initial values are always needed to start any non linear optimisation process (which is what is used here, see Sect. 2.2.5), they were set to the values obtained by using the laws proposed by Guitard [30] and modified by Viala [5] for spruce and maple tonewood, based on a sample of 60 specimens. The Guitard laws were indeed more general for soft and hard woods and did not take into account the specificity of spruce and maple species neither the fact they are selected for instrument making (lack of defects and resin pockets, relatively low angle of the crystalline cellulose micro-fibrils in the S2 layer of the secondary wall). These relations are provided in Table 2. The Poisson’s ratios and E_T were given the following nominal values: ν_LR = ν_RT = ν_TR = ν_LT = 0.3 and ν_RL = ν_TL = 0.03 (spruce and maple), $E_{\mathrm{T}}^{\mathrm{s}}=520$ MPa (spruce) and $E_{\mathrm{T}}^{\mathrm{m}}=766$ MPa (maple) [20].

2.2.3 Sensitivity analysis

The optimisation is based on matching the simulated eigenfrequencies with the experimental frequencies for a selection of modes. Including too many modes in the optimisation can make it too cumbersome, while increasing the risk of mode inversions and more local minima. It was therefore decided to run it on only a limited number of modes, which would be chosen after checking the sensitivity of these modes to the material parameters. Computing this sensitivity is essential for the reciprocal reason: only the wood properties that are influential on the selected eigenfrequencies can be obtained by optimisation in a reliable way. Therefore, a sensitivity analysis was conducted to determine the influence on the six first frequencies of each of the 10 parameters of the wood as modeled in Section 2.2.2. This was done by varying each parameter by 1% and determining, for each mode, the induced relative variation of the eigenfrequency. Then, for each mode, these 10 relative variations were normalised so their sum equalled 1. The resulting values are shown in Figure 4 for one top plate (left) and one back plate (right). E_L, E_R, G_LR and G_RT were found to be moderately to strongly influential and so these four parameters were included in the optimisation procedure, detailed in the following section. In contrast, the influence of E_T, G_TL and the three Poisson’s ratios appeared very weak, so weak that, except for G_TL they are not even represented in Figure 4. They will therefore not be included in the optimisation and will be given the initial values provided in Section 2.2.1.

Finally, Figure 4 shows that while four modes could be enough for optimising E_L, E_R and G_LR, six modes are needed for G_RT. The optimisation was thus based on the first six modes and is detailed hereafter.

2.2.4 Density evaluation

While, among all wood parameters, the density has the largest influence on the vibratory response of the plate, the system is underconstrained: a change in density can easily be compensated by a change in elastic parameters (as they have opposite effects, see Fig. 4) and so a solution can be found by optimisation for any density with a small error. Therefore, the density needs to be fixed a priori. As the density measurements were made on the wedges and not the plates, and knowing that the density can be heterogeneous, the density was estimated for each plate. The weight recorded by the Bilbao team was divided by the volume of the object defined in the geometry files (190 cm³ for the top and 213 cm³ for the back, the back being thicker in the central area). The estimations are provided in Table 3.

While the estimated values for the backs are in good agreement with the measured densities (with a mean relative error of 1.2%), the estimated values for the tops significantly deviate from the measured ones, with a mean relative variation of 10%. One hypothesis for this discrepancy lays in the density difference between sapwood and heartwood for spruce, sapwood being denser [31].

To make a violin plate, the wedge is divided in two and the two pieces are glued together so the sapwood is in the center and the heart wood on the outside. When a violin plate is carved in this new piece, more heartwood is thus removed than sapwood which may increase the averaged density of the plate compared to the density measured on the quarter.

It has to be noted that the coefficient of variation of tops and backs is close to 1.8% and 2.7% respectively. The usual variability of tonewood is higher, even when specially selected for violin making (from 7% to 13% [20, 32, 33]). This suggests a relatively high similarity between the selected woods, as intended by the Bilbao team.

These density estimations were then used to determine the initial values of the wood properties (given in Table B1) for the optimisation process, using the laws provided in Section 2.2.2. Note that the initial values for E_L correspond quite closely to the measured values based on the first compression eigenfrequency. This suggests that measurements with the Lucchimeter can overestimate by 1–3 GPa the longitudinal Young Modulus, confirming previous observations [28].

2.2.5 Optimisation, cost function and modal assurance criterion (MAC)

The optimisation was done by using the Finite Element Model Updating Approach (FEMU) [34] via Matlab’s fminsearchbnd function which allows to specify a range for each of the material parameters to be optimised. We have chosen reference ranges from the literature that are provided in Tables 4 and 5. We believe that, as there is an error on the measured frequencies anyway (i.e. hard to quantify), it is better to have a solution within the typical range for wood properties even if a better match could possibly be obtained with values outside this range (but would be physically very unlikely).

The optimisation process relies on minimising a cost function. Traditionally, this cost function is based on the error that represents the average relative discrepancy between the simulated frequencies (f_sim,i) and the experimental frequencies (f_exp,i) over a certain number of modes [35], here the first n modes as justified in the previous section. This ε_freq quantifies how closely the simulated frequencies match the experimental ones:

$${\epsilon }_{\mathrm{freq}}=\frac{100}{n}\sum _{i=1}^n\mathrm{}\frac{|f_{\exp ,i}-f_{\mathrm{sim},i}|}{f_{\exp ,i}}.$$(2)

However, a low error on the eigenfrequencies is not enough to validate the optimisation. Indeed, the simulated mode shapes should be similar to the measured ones too. The Mode Assurance Criterion (MAC) is a widely used tool to assess the correspondence between experimentally obtained vibration mode vectors and those resulting from numerical simulations [36]. The MAC provides a numerical index between 0 and 1, where a value of 1 indicates a perfect match between the vibration modes that are compared, and is defined as follows:

$$\mathrm{MAC}({\varphi }_i,{\psi }_j)=\frac{|{\varphi }_i^T{\psi }_j{|}^2}{\left({\varphi }_i^T{\varphi }_i\right)\left({\psi }_j^T{\psi }_j\right)}$$(3)

where ϕ_i and ψ_j represent the mode vectors of the i-th experimental modal shape and the j-th simulated modal shape, respectively. As the number of simulation points is much higher than the number of experimental points, the nearest simulation point is searched for each experimental point to compute the MAC matrices, which leads to an average error of 0.8 mm (see Fig. 5), far below the wavelengths in the considered frequency range (8–50 cm) (see Appendix A).

Neglecting the MAC during optimisation could lead to a scenario where the eigenfrequencies appear to be well matched, but the mode shapes differ significantly. Such a difference could compromise the modal integrity of the system, leading to an erroneous interpretation of the vibrational behavior. Figure 6 illustrates, with the example of back B, what can happen when the optimisation is solely based on minimising ε_freq: though ε_freq has slightly decreased from 2% to 1.6%, the MAC matrix has actually worsened.

We have therefore integrated the proximity of the simulated modes compared to the experimental modes in three ways in the cost function.

First, for the part of the cost function C_freq that is related to the match in terms of eigenfrequencies, only the modes for which the MAC is above 0.7 were included:

$$\begin{array}{cc}{C}_{\mathrm{freq}}=\frac{100}{N_I}\sum _{i\in I}\mathrm{}\frac{|f_{\exp ,i}-f_{\mathrm{sim},i}|}{f_{\exp ,i}}& \mathrm{with} I=\left\{i\in \left[\mathrm{1,6}\right]\right| \mathrm{MAC}({\varphi }_i,{\psi }_i)>0.7\}\end{array}$$(4)

and N_I the length of I.

Second, a penalty C_mode based directly on the MAC matrix was added. It was defined to include both the diagonal elements which should be as close to 1 as possible, to ensure a good match between the experimental and the simulated mode shapes, and the non diagonal elements, which should be as small as possible, to minimize cross-modal influences and ensure distinct modal identification [37]:

$$C_{\mathrm{mode}}=\alpha \cdot D+\beta \cdot N$$(5)

with

$$D=\sum _{i=1}^6\mathrm{}(1-\mathrm{MAC}({\varphi }_i,{\psi }_i){)}^2$$(6)

and

$$N=\sum _{i=1}^6\mathrm{}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^6\mathrm{}\mathrm{MAC}({\varphi }_i,{\psi }_j{)}^2.$$(7)

The coefficients α and β adjust the relative importance of the penalties D and N, allowing for fine-tuning of the optimization process to balance mode matching accuracy and minimization of cross-modal contamination.

Third, a penalty C_number was added to reflect a potential decrease in the modeling accuracy when N_I, the number of matched modes (MAC > 0.7), gets too far from 6. It was therefore defined as

$$C_{\mathrm{number}}=\gamma \left(\frac{1}{N_I}-\frac{1}{6}\right).$$(8)

The total cost function used in the optimisation was therefore:

$$C=C_{\mathrm{freq}}+C_{\mathrm{mode}}+C_{\mathrm{number}}$$(9)

α, β and γ were calibrated to 1, 5 and 60 respectively. This was done empirically to ensure comparable weights of the three penalties as the interplay of these three components is needed for a good match between the experimental results and the simulation. It is indeed essential to achieve a good frequency match without distorting the MAC, maintaining the integrity of the mode shapes.

Figure 7 shows the MAC matrix of back B after optimisation with the total cost function C. While the final value of ε_freq is very similar to the value obtained before (1.8%), the MAC matrix has not been degraded.

3 Results and discussion

3.1 Identification of the elastic constants E_L, E_R, G_LR and G_RT

The values obtained by optimisation using finite element model updating for the four elastic constants are provided in Tables 4 and 5 for the tops and backs respectively, along with the range that has been reported for tonewood spruce and maple [20] (provided in the last line).

The values are very similar across the plates showing again the careful selection of the wood for the Bilbao project. They are in general quite close to the ones expected from Viala’s laws though consistently a bit lower for E_L. This is not understood so far, all the less that the expected values were in really good agreement with the measured values (based on the frequency of the first compression mode).

To test the robustness of the optimisation, different initial values were used, deviating by 15% from Viala’s laws. The resulting optimised elastic constants were exactly the same.

The microfibril angle (MFA) of the S2 layer of the secondary wall of tracheids (defined in relation to the longitudinal direction) is partially responsible for the stiffness as microfibrils of cellulitis are the stiffer parts of a wood cell, with a young modulus close to 13 GPa: the lower the angle, the higher the stiffness in the longitudinal direction, and consequently, the lower the stiffness in the radial and tangential directions [38–40]. Therefore having E_L in the lower bound will generally lead E_R (and E_T) to be in the upper bound, for a given density.

3.2 Simulated modes

The results of the optimisation process in terms of eigenfrequencies and mode shapes are shown in Figure 8 for the top A (Fig. 8a) and back A (Fig. 8b). The comparison with the measured data (Fig. 2) reveals a strong similarity in frequencies and modal shapes for both plates. This is the case for the other 10 plates (see Figs. B1–B10 in Appendix B) and will be quantified below.

3.2.1 Eigenfrequencies error

The eigenfrequency error as calculated with equation (2) is provided before and after optimisation for each plate in Table 7. The average error is around 0.9 and 1.2% for tops and backs respectively, which is lower than the measurement precision (repeatability error and humidity variations) [41].

It is interesting to note that ${ϵ}_{\mathrm{freq}}^0$ is already quite low for the 12 plates. This shows the validity of the laws (Table 2) established by one of the authors (RV) for a large corpus of tonewood samples [5] and that the initial values derived from these laws (Table B1) can be used for estimating the unknown parameters and lead to a good prediction of the eigenfrequencies.

It is also possible to predict the behavior of the plates outside the identification domain [42–47], as illustrated by the results presented in Appendix C.

3.2.2 Modal assurance criterion

Diagonal MAC values, i.e. the similarity between the simulated and experimental i-th modes are provided for the six first modes, for each of the 12 plates in Table 8. The values are overall quite high, especially for the top plates. The minimum MAC value for the top plates is above 0.74 (mode 6 of top D), which represents a relatively good match between experiments and simulations. For the back plates, however, the minimum value is much lower, 0.41 (mode 6 for back E) which is rather low. But all other values are above 0.58. This difference between the tops and the backs may be explained by the fact that the experimental grid for the tops is much finer (70 points) than for the backs (43 points).

MAC matrices obtained before (with the initial values defined in Sect. 2.2.1) and after optimisation are presented in Figure 9, only for top A as Figures 5 and 6 already show an example for a back plate, and the MAC matrices are very similar for the other plates. Thanks to the MAC penalty, the MAC matrices are not degraded by the optimisation, but are barely improved. Therefore, when optimisation is not possible, using Viala’s laws offer an excellent alternative, leading to a good prediction for both the eigenfrequencies (as seen in Table 7) and mode shapes.

4 Conclusions

A FE model can be a potential tool for violin makers to predict the effect of small variations which would be difficult (if not impossible) to study in practice. Indeed, this would be a natural follow-up of the current trend that can be observed in the violin making community that is acquainted with acoustics: with the growing availability of free and easy to use FEM online platforms, some violin makers have already started using FE models to predict the vibratory behavior of small pieces like, for instance, the bridge. Many little details of the assembly process may be difficult to model but the model only needs to be accurate within the variability due to hand work and natural variations in wood. The Bilbao project offers an unprecedented database to calibrate the FE model at various stages during the construction of a violin. This paper presented the first stage of building this model, on the 12 plates, after they were CNC carved, before graduation. The geometry was obtained from the CAD file used for the CNC routing while an optimisation procedure was developed to optimise the wood’s material parameters in order to match the simulated eigenfrequencies of the first six modes with the experimental ones. The optimisation used a cost function based on the eigenfrequencies as well as the MAC values in order to ensure the proximity of the simulated mode shapes with the experimental ones. The obtained material parameters are close to the expected values (though consistently lower for E_L), the error on the eigenfrequencies is very small and the match in terms of modal shapes is good for the backs and excellent for the tops as shown by the MAC matrices. The difference between the two types of plates may be attributed to a smaller number of experimental measuring points for the backs. It has been shown as well that using Viala’s laws lead to a good (and even very good in some cases) prediction and can offer thus a useful alternative to optimisation when the latter is not possible. However, one should be aware that some wood can deviate strongly from these laws, especially for low densities. This is significant because violin makers often prefer lighter woods, such as torrefied maple, where the flames can alter the mechanical properties, making predictions challenging [48].

This study validates the first step in obtaining a FE model that can be sufficiently close to be used as a tool to quantify real modifications on plates. In particular, it can be used to evaluate the impact of a luthier’s design choices on the vibratory behaviour of the plates as this would usually not be feasible in reality. Indeed, the number of plates that would be required to allow the separation of the effects due to intentional modifications from effects due to unintentional variations would be too large, considering the irreducible uncertainty in geometry and wood parameters.

The next steps in this project include the numerical graduation of the plates’ thickness (taking into account that the thickness of violin plates is not homogeneous [13] and that different graduation patterns were used for the six violins of the Bilbao project) and then their assembly which involves modelling blocks, linings, the effect of gluing, adding ribs and a soundpost. Since we have all the experimental data for the finished violins, this project presents an exceptional opportunity to study how accurate FE modelling can be in simulating complete instruments. From the results on the individual plates, the results look promising.

Acknowledgments

The experimental part of the project was partially funded by Erasmus+. The authors would like to thank Sam Zygmuntowicz, the violin maker in charge of maintaining the Huberman violin, for providing the CT-scan.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

The data are available from the corresponding author on request.

Appendix A

For the optimisation, the first six modes considered are all below 1000 Hz. Below are the detailed calculations used to determine the corresponding wavelengths in the longitudinal (L) and radial (R) directions for spruce, an orthotropic material.

The flexural rigidity in the longitudinal direction is defined as [49, 50]:

$$D_L=\frac{E_L\cdot h^3}{12\cdot (1-{\nu }_{\mathrm{LR}}\cdot {\nu }_{\mathrm{RL}})}$$

where E_L is the Young’s modulus in the longitudinal direction, h is the plate thickness, and ν_LR and ν_RL are the Poisson’s ratios for longitudinal-radial and radial-longitudinal directions, respectively. The flexural wave speed in the longitudinal direction is calculated by:

$$v_{\mathrm{flexural},L}={\left(\frac{D_L\cdot {\omega }^2}{\rho \cdot h}\right)}^{\frac{1}{4}}$$

where ω is the angular frequency and ρ is the material density. The flexural wavelength in the longitudinal direction is then given by:

$${\lambda }_{\mathrm{flexural},L}=\frac{v_{\mathrm{flexural},L}}{f}$$

where f is the frequency. Similarly, the flexural wavelength in the radial direction is given by:

$${\lambda }_{\mathrm{flexural},R}=\frac{v_{\mathrm{flexural},R}}{f}.$$

Two sets of material properties (corresponding to spruce and maple) and eigenfrequencies were used to calculate the flexural wavelengths. The results are summarized in Table A1 for spruce and Table A2 for maple. A homogeneous thickness of the plates of 3.5 mm was considered for each case.

Appendix B

The initial values used for the optimisation are provided in Table B1 while Figures B1–B10 show the measured and simulated mode shapes and eigenfrequencies, for all the plates (except top A and back A that are already provided in the main text).

Figure B1 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate B. (b) Simulated modal shapes and eigenfrequencies of top plate B. Displacement along the Z axis is displayed.

Figure B2 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate B. (b) Simulated modal shapes and eigenfrequencies of back plate B. Displacement along the Z axis is displayed.

Figure B3 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate C. (b) Simulated modal shapes and eigenfrequencies of top plate C. Displacement along the Z axis is displayed.

Figure B4 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate C. (b) Simulated modal shapes and eigenfrequencies of back plate C. Displacement along the Z axis is displayed.

Figure B5 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate D. (b) Simulated modal shapes and eigenfrequencies of top plate D. Displacement along the Z axis is displayed.

Figure B6 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate D. (b) Simulated modal shapes and eigenfrequencies of back plate D. Displacement along the Z axis is displayed.

Figure B7 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate E. (b) Simulated modal shapes and eigenfrequencies of top plate E. Displacement along the Z axis is displayed.

Figure B8 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate E. (b) Simulated modal shapes and eigenfrequencies of back plate E. Displacement along the Z axis is displayed.

Figure B9 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate F. (b) Simulated modal shapes and eigenfrequencies of top plate F. Displacement along the Z axis is displayed.

Figure B10 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate F. (b) Simulated modal shapes and eigenfrequencies of back plate F. Displacement along the Z axis is displayed.

Appendix C

This appendix provides informations about the predictivity of the model outside the frequency range used for the optimisation. The simulated modes 7–13 of the top plate B are compared to the experimental modes (Fig. C1) and the error on the eigenfrequencies are given in Table C1.

Figure C1 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate B. (b) Simulated modal shapes and eigenfrequencies of Top plate B. Displacement along the Z axis is displayed.

The calculated eigenmodes are outside the training domain; however, the errors remain below 4% in all cases, with some modes showing an almost perfect correspondence, as shown in Table C1. This suggests that optimization over more parameters in the space is possible. Nevertheless, it is important to note that the most significant vibration modes are the low-frequency ones, as they better describe the system’s characteristics due to their lower energy dissipation. Therefore, these modes are prioritized to understand the system’s behavior.

References

All Tables

All Figures

	Figure 1 Left: Setup for performing the modal analysis of the free plates. Right: Grid of the points at which the transfer functions were measured, on the top (left, 70 points) and the back (right, 43 points) plates.
In the text

	Figure 2 Interpolated modal shapes of the first six modes and corresponding eigenfrequencies measured on free top A (a) and back A (b). The displacement along the Z axis (see Fig. 3) is displayed.
In the text

	Figure 3 Actual mesh of the top (left) and back (right) plates of the Bilbao violins, generated with tetrahedral elements in COMSOL.
In the text

	Figure 4 Sensitivity analysis matrix for a 1% variation on each of relevant material parameters (density and elastic constants). Results for top C (left) and back C (right). The values are normalized for each mode. A high value is associated with a large influence and the sign indicates the direction in which the eigenfrequencies are affected.
In the text

	Figure 5 (a) Distribution of simulated (red dots) and experimental points (blue circles), (b) Distribution of the experimental points (blue circles) and their closest simulated neighbours (red circles), highlighting the alignment between them.
In the text

	Figure 6 MAC matrix for back B before optimisation (a) and after the optimisation (b) using only εfreq for the cost function.
In the text

	Figure 7 MAC matrix for back B after the optimisation using the cost function C.
In the text

	Figure 8 Modal shapes of the first six simulated modes and corresponding eigenfrequencies for the free top A (a) and back A (b). The displacement along the Z axis is displayed.
In the text

	Figure 9 MAC matrices for top A before optimisation (a) (with the initial values of the wood properties described in Sect. 2.2.1); and after optimisation (b).
In the text

	Figure B1 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate B. (b) Simulated modal shapes and eigenfrequencies of top plate B. Displacement along the Z axis is displayed.
In the text

	Figure B2 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate B. (b) Simulated modal shapes and eigenfrequencies of back plate B. Displacement along the Z axis is displayed.
In the text

	Figure B3 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate C. (b) Simulated modal shapes and eigenfrequencies of top plate C. Displacement along the Z axis is displayed.
In the text

	Figure B4 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate C. (b) Simulated modal shapes and eigenfrequencies of back plate C. Displacement along the Z axis is displayed.
In the text

	Figure B5 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate D. (b) Simulated modal shapes and eigenfrequencies of top plate D. Displacement along the Z axis is displayed.
In the text

	Figure B6 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate D. (b) Simulated modal shapes and eigenfrequencies of back plate D. Displacement along the Z axis is displayed.
In the text

	Figure B7 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate E. (b) Simulated modal shapes and eigenfrequencies of top plate E. Displacement along the Z axis is displayed.
In the text

	Figure B8 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate E. (b) Simulated modal shapes and eigenfrequencies of back plate E. Displacement along the Z axis is displayed.
In the text

	Figure B9 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate F. (b) Simulated modal shapes and eigenfrequencies of top plate F. Displacement along the Z axis is displayed.
In the text

	Figure B10 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free back plate F. (b) Simulated modal shapes and eigenfrequencies of back plate F. Displacement along the Z axis is displayed.
In the text

	Figure C1 (a) Interpolated modal shapes and corresponding eigenfrequencies measured on the free top plate B. (b) Simulated modal shapes and eigenfrequencies of Top plate B. Displacement along the Z axis is displayed.
In the text