*Proceeding Paper* **Computational Intelligence Model of Orally Disintegrating Tablets: An Attempt to Explain Disintegration Process †**

**Jakub Szl ˛ek \*, Adam Pacławski, Natalia Czub and Aleksander Mendyk**

Department of Pharmaceutical Technology and Biopharmaceutics, Jagiellonian University Medical College, Medyczna 9, 30-688 Kraków, Poland; adam.paclawski@uj.edu.pl (A.P.); natalia.czub@doctoral.uj.edu.pl (N.C.); aleksander.mendyk@uj.edu.pl (A.M.)

**\*** Correspondence: j.szlek@uj.edu.pl; Tel.: +48-12-620-56-04

† Presented at the 2nd International Electronic Conference on Applied Sciences, 15–31 October 2021; Available online: https://asec2021.sciforum.net/.

**Abstract:** We obtained a curated database based on the database published elsewhere. Chemical descriptors were introduced as characteristics of active pharmaceutical ingredients (APIs). We used H2O AutoML platform in order to develop a Deep Learning model and SHAP method to explain its predictions. Obtained results were satisfactory with *NRMSE* of 8.1% and R2 of 0.84. Finally, we identified critical parameters affecting the process of disintegration of directly compressed ODTs.

**Keywords:** machine learning model; computational intelligence; AutoML; orally disintegrating tablets; ODTs; disintegration time

#### **1. Introduction**

Traditional tablets are not an ideal drug dosage form. Many groups of patients, e.g., pediatric or geriatric patients, have problems with swallowing or simply are not willing to take tablets. As a consequence, all these factors may reduce a patient's compliance. In order to overcome the inconvenience of conventional tablet use, orally disintegrating tablets (ODTs) were introduced into the drug market. One of the methods of preparing ODTs is direct compression, which is cost efficient and simple. It involves comparatively fewer stages than compression preceded by wet or dry granulation. In brief, powders are grinded, if necessary, and blended. Then, the mixture is compressed into the tablets. Although the process is quite simple, there are many factors that influence the characteristics of the ODTs. Among the crucial factors is the disintegration time.

One of the methods used to solve problems with many factors, where the hypothesis governing the phenomenon is unknown or the whole process is complex, is machine learning (ML). Automated machine learning (AutoML) is currently in focus as a branch of ML automating the time-consuming, iterative tasks of model development. AutoML enables machine-driven building of large-scale, high-performance, and superb predictability models with minimal human intervention.

The Motivation of this study was a limited knowledge of relationships between excipients, APIs, and process parameters of direct compression and their influence on disintegration of ODTs. Knowing such behavior would enhance the design and development of novel drug dosage forms. In this work we applied a concept of AutoML-based heuristic model development for prediction of the disintegration time based on the quantitative and qualitative composition of powder mixtures.

#### **2. Materials and Methods**

Our database was built based on the database presented by Han et al. [1]. First, we curated the existing database [1], neglecting any unclear or uncertain data records. We put emphasis on the occurrence of the ODTs' characteristic and process parameters,

**Citation:** Szl ˛ek, J.; Pacławski, A.; Czub, N.; Mendyk, A. Computational Intelligence Model of Orally Disintegrating Tablets: An Attempt to Explain Disintegration Process. *Eng. Proc.* **2021**, *11*, 24. https:// doi.org/10.3390/ASEC2021-11163

Academic Editors: Nunzio Cennamo and Takayoshi Kobayashi

Published: 15 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

such as tablet hardness, thickness, and dimension of tablet press die. Moreover, we performed a literature survey in order to enhance the database. Scopus® database was searched for publications fulfilling the following criteria: the direct compression method of ODTs should be used in processing, the amount of all excipients should be present, tablet characteristics (hardness, thickness, and die dimension) should be present, and the compendial disintegration test should be applied (Ph. Eur. or USP).

After data scrapping, we calculated APIs' two-dimensional (2D) molecular descriptors using mordred-descriptor v.1.2.1a1 Python package [2] and included them in the curated database. Excipients' types and amounts were encoded in a topological manner. The only output was the time needed for the disintegration of tablets.

A Computational experiment was performed according to the scheme presented in Figure 1. In brief, a preprocessed database was passed to the Python script [3] performing at the first stage a feature selection. Then, the final model was built according to a 10-fold cross-validation scheme. All available algorithms in H2O implementation of AutoML were used [4]: Distributed Random Forest (DRF), Extremely Randomized Trees (XRT), Generalized Linear Model (GLM), Extreme Gradient Boosting Machine, (XGBoost), Gradient Boosting Machine (GBM), Deep Learning (fully connected multi-layer artificial neural network), and Stacked Ensemble models.

**Figure 1.** Scheme of computational experiment design.

Model performance was assessed according to the 10-fold cross-validation (10-CV) and expressed by three goodness of fit metrics: root-mean-square error (*RMSE*), normalized root-mean-square error (*NRMSE*), and coefficient of determination (*R*2). For reference, please see Equations (1)–(3).

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{n} (pred\_i - obs\_i)^2}{n}} \tag{1}$$

where *obsi* and *predi* are observed and predicted values, *i* is the data record number, and *n* is the total number of records.

$$NRMSE = \frac{RMSE}{obs\_{max} - obs\_{min}} \times 100\% \tag{2}$$

where *RMSE* is the root-mean-square error and *obsmax* and *obsmin* are observed minimal and maximal values.

$$R^2 = 1 - \frac{SS\_{res}}{SS\_{tot}} = 1 - \frac{\sum\_{i=1}^{n} \left(pred\_i - obs\right)^2}{\sum\_{i=1}^{n} \left(obs\_i - obs\right)^2} \tag{3}$$

where *R*<sup>2</sup> is the coefficient of determination, *SSres* is the sum of squares of the residual errors, *SStot* is the total sum of the errors, *obsi* and *predi* are observed and predicted values, and obs is the arithmetical mean of observed values.

Predictions of the best model were explained with the use of another Python wrapper [5] implementing among others, including the SHapley Additive exPlanations (SHAP) method by Lundberg et al. [6].

#### **3. Results**

Each record of the curated database represented one formulation of ODTs. It consisted of 633 chemical descriptors encoding API, 28 inputs encoding amounts of excipients, and 9 inputs characterizing drug dosage form. A single, independent variable was the disintegration time. The database consisted of 243 records (formulations), of which only 52 records (~21%) overlapped the Han et al. database [1].

In the feature selection stage, inputs' number was reduced to 39, among which there were 28 inputs (amount of 2-hydroxypropyl-beta-cyclodextrin [%], Aerosil [%], Amberlite IRP 64-69 [%], API [%], beta-cyclodextrin [%], calcium silicate [%], camphor [%], colloidal silicon dioxide [%], croscarmellose sodium [%], crospovidone [%], cyclodextrin methacrylate [%], Eudragit EPO [%], hydroxy propyl methyl cellulose [%], lactose [%], low-substituted hydroxy propyl cellulose [%], magnesium stearate [%], mannitol [%], microcrystalline cellulose [%], Poloxamer 188 [%], polyvinyl alcohol [%], polyvinylpyrrolidone [%], pregelatinized starch [%], sodium bicarbonate [%], sodium carboxymethyl starch [%], sodium lauryl sulphate [%], sodium starch glycolate [%], sodium stearyl fumarate [%], and talc [%]) responsible for encoding the quantity of excipients and API, 8 molecular descriptors characterizing API (API Geary autocorrelation of lag 7 weighted by ionization potential, API topological charge index of order 7, API Geary autocorrelation—lag 7/weighted by polarizabilities, API modified information content index, API Moran autocorrelation of lag 4 weighted by polarizability, API negative logarithm of the partition (oil/water) coefficient, API number of 12-membered rings (includes counts from fused rings), API number of 8-membered fused rings containing heteroatoms (N, O, P, S, or halogens)), and 3 inputs characterizing the drug dosage form (diameter of die or tablet [mm], hardness of ODT [N], thickness of ODT [mm]). A list of selected features along with their type and relative importance is presented in Table 1.


**Table 1.** The First 15 selected features and their relative importance.

The best results were obtained by a Deep Learning (DL) model, which had *RMSE* = 10.9, *NRMSE* = 8.1%, and *R*<sup>2</sup> = 0.84. The model had two hidden layers with 100 neurons in each layer and a rectifier with dropout as an activation function. A plot of predicted versus observed disintegration values is presented in Figure 2.

**Figure 2.** Predicted vs. Observed values for disintegration time for Deep Learning model.

Following model development, a procedure of the SHAP method was applied. Then, selected plots were analyzed, and conclusions were drawn (Figure 3).

**Figure 3.** Results of the model's explanation: (**a**) Summary plot of impact on model output and feature value; (**b**) Effect of magnesium stearate amount [%] on the average model's prediction.

#### **4. Discussion**

Based on the obtained prediction metrics (*RMSE*, *NRMSE*, *R*2), it was concluded that the model was satisfactory in terms of generalization. The 10-fold cross-validation technique was used as a golden standard. The mean error of the model was 10.9 (*NRMSE* = 8.1%); therefore it was possible to optimize a formulation with its use. Moreover, in Figure 3a, critical parameters and their impact on the disintegration time were identified. It seems that a high amount of sodium lauryl sulphate, magnesium stearate, Eudragit EPO, and colloidal silicon dioxide could increase the disintegration time of ODTs. On the other hand, a high amount of crospovidone, Aerosil, croscarmellose, sodium starch glycolate, or sodium stearyl fumarate could lead to a decreased disintegration time. Looking more closely at the variable effects, a percolation threshold could be found. For example, in Figure 3b, at a magnesium stearate value of about 1% a reverse in effects could be observed. This observation was consistent with

the findings of previous studies [7]. It is believed that magnesium stearate in higher amounts than 1%, besides the usual action as a lubricant, could form a hydrophobic film around API particles and could prevent water from penetrating into the core of the tablet. Using similar reasoning, the XLogP limit was determined for the API, the value of which will increase the disintegration time of ODTs (Figure 4). The general conclusion is that a more hydrophobic API with XLogP higher than 3.5 would negatively affect the disintegration time by increasing it.

**Figure 4.** Effect of XLogP (calculated partition coefficient of API) on the average model's prediction.

**Author Contributions:** Conceptualization, J.S. and A.M.; methodology, J.S., A.P., N.C. and A.M.; software, J.S.; validation, J.S., A.P., N.C. and A.M.; formal analysis, J.S.; investigation, J.S.; resources, A.M.; data curation, J.S.; writing—original draft preparation, J.S., A.P., N.C. and A.M.; writing—review and editing, J.S., A.P., N.C. and A.M.; visualization, J.S.; supervision, J.S.; project administration, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Uniwersytet Jagiello ´nski Collegium Medicum, grant number N42/DBS/000205.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available upon written request.

**Acknowledgments:** Not applicable.

**Conflicts of Interest:** The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. The authors declare no conflict of interest.

#### **References**


## *Proceeding Paper* **Concise Review of Classical Guitar Modelling Technologies †**

**Alexandre M. Löw 1,\*, Herbert M. Gomes <sup>1</sup> and César M. A. Vasques <sup>2</sup>**


**Abstract:** Analytical modeling and numerical simulation of multiphysics coupled systems is an exciting research area, even when it comes to intrinsically linear or linearized formulations, as is usually the case with coupled vibroacoustic problems. The combined effect of many localized geometrical miss-modeling with significant uncertainty in mechanical characterization of some organic materials yields large discrepancies in the natural frequencies and mode shapes obtained. The main goal of this work is to compare two basic approaches for the modeling of stringed musical instruments in the frequency domain: simplified lumped-parameter analytic modeling, considering only the most influential degrees of freedom, and discretized finite element modal analysis. Thus, the emphasis is on a review of some key references in this field, including previous work by the authors, which may shed light on some of the most relevant issues surrounding this problem.

**Keywords:** classical guitar; lumped-parameter modeling; finite element method; vibroacoustics; frequency domain; musical instruments; luthier

#### **1. Introduction**

Due to the diversity of phenomena involved in the modeling of musical instruments and the inherent difficulties in formulating vibroacoustics problems, including the various underlying possible dissipating effects, and despite the availability of some classical and computational approaches, this task remains a significant challenge after many years. It is amazing that despite the huge global popularity of plucked stringed instruments like the guitar, scientific research into the piano and violin began much earlier, paving the way for the application of similar principles to the modeling and study of classic guitars.

One may wonder, however, to what extent are analytical models truly exact? Though convergence towards the analytic solution is mathematically granted for methods such as the finite or boundary element methods [1], to what extent can we call this result correct? This article is intended to discuss the above questions contrasting two major approaches when it comes to acoustic guitar modelling. On one hand, a quantitative analysis can follow a familiar approach among engineers and more practically oriented researchers: the quest for an equivalent model capable of yielding a sufficiently accurate approximate result, given that model parameters lie inside a broad-enough range of permitted variations this paradigm is eventually called gray box modelling, because it does not rely entirely neither on experimental identification nor on analytical modelling. Usefulness has proven to be the key feature here, with the loss of accuracy perceived as a minor drawback in comparison with some gains expected for a computationally inexpensive model. On the other hand, lies the demand for high-fidelity descriptions of physical reality, where a number of numerical methods could be employed to discretize the system into a set of algebraic linear equations—with finite element (FE) models standing by far as the preferred

**Citation:** Löw, A.M.; Gomes, H.M.; Vasques, C.M.A. Concise Review of Classical Guitar Modelling Technologies. *Eng. Proc.* **2021**, *11*, 25. https://doi.org/10.3390/ ASEC2021-11179

Academic Editor: Filippo Berto

Published: 15 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

option. Thus, bigger descriptions can be obtained, even though they should unavoidably require greater sets of input data, certainly contaminated by some amount of uncertainty. Thus, model updating has become a standard practice both in academic research and industry (commercial software) applications [2]. It is important to note from its definition, however, that the whole model uncertainty remains distributed over the parameters in an unpredictable, possibly overlapped manner, consequently raising serious concerns about the uniqueness of obtained solutions [3]. These issues become evidently more severe for huge complex models with possibly hundreds of free parameters. Thus, it is clear that a delicate trade-off emerges at this point for those situations where a smaller, less detailed model, could still bring useful results while providing valuable gains from narrower bounds on uncertainty propagation.

#### **2. Simplified Guitar Models**

Even as a complex and truly intriguing system, the construction of musical instruments still relies almost solely on traditional knowledge in the great majority of famous lutheries, with an increasing number of academic researchers trying to establish this bridge out of a personal passion [4]. However, things were slightly different for a guitar's close relative: the violin. For a long time, being a learned instrument, it has naturally attracted researchers' attention earlier than its plucked-string cousin. Thus, Schelleng [5] proposed, already in the 1960's, an analogue 2nd order electro-magnetic resistance–inductance–capacitance circuit that could emulate the elasto-acoustic interactions of top plate mechanical properties mass distribution, stiffness and damping—with the air confined inside and surrounding the resonant chamber. Firth [6] adapted the analogue circuit to the dynamics of acoustic guitar, relating electrical quantities to the most clearly important structural and acoustical parameters. While valuable for those familiar with electrical–mechanical analogy (a standard theoretical apparatus in classical control theory) this approach eventually gave place to a rather direct description of the coupled phenomena. Thus, Caldersmith [4] devised, constructed and measured a simple resonator composed by parts with similar functions as those pointed in [6], and Christensen [7], relying upon these concepts and related equations [4–6], presented a simplified 2 degrees of freedom (DOF) model derived from purely mechanical considerations, focusing on the relations between 1st and 2nd coupled modes and their intermediate Helmholtz anti-resonance, where almost all model parameters (with exception of those related to damping) can be obtained from simple measurement of these 3 frequencies. This work greatly simplified the resulting model, thereby bringing timely applicability and physical insight to subsequent developments. In this way, Caldersmith [8] introduced a straightforward yet important modelling refinement by taking into account the contribution of back plate flexibility in the coupled dynamics, resulting in a 3 DOF system for the resonant chamber, and French [9,10] tested a procedure for identification of model parameters and prediction of structural modification based on the sensitivities of eigenvalues for both 2 and 3 DOF models. Finally, the last extension to this simplified modelling approach was presented in [11], where additional mass and motion DOF were considered for the ribs, neglecting, however, any associated stiffness in this region. Following [7], here as well, model parameters can be obtained, with the exception of dissipative coefficients, from simple measurements. At the top and back plate, bending is represented by that of a thin axisymmetric disk corresponding approximately to the region circumscribed by the lower bout, in an average condition between hinged and clamped. In this way, theoretical approximate expressions could be obtained for the enclosed air volume variation as the product of equivalent plate area and displacement.

#### *2.1. Formulation of Simplified Models and Definition of Parameters*

In what follows, the equations of these simplified models discussed above will be presented, trying to keep clarity and conciseness as possible. The meaning of each index adopted in the development of the equations below is shown in Table 1 (for an explanation of guitar components, refer to Figure 1 below). Table 2 presents the parameters appearing in the equations of the 4 DOF model [11]; it should be noted that those appearing in the equations for 2 and 3 DOF models [7,8,10] are also covered in this table, but some table entries do not apply for the simpler models. Damping coefficients are presented separately on the last line because these parameters are to be defined through experimental modal analysis.

**Table 1.** Meaning of the indexes used in the various vibroacoustic models of the classic guitar.




About the damping parameters in Table 2, it is important to emphasize at this point that these coefficients alone would imply, for instance, a purely viscous, or viscoelastic/hysteretic dissipation behaviour, which is clearly a non-physical assumption considering the predicted forms of energy loss in the modelling of vibro-acoustic coupled systems; i.e., the energy balance in this cases should include, at least, some sort of viscous, viscoelastic and acoustic radiation forms of damping [12,13]. Indeed, for those parameters related to energy dissipation, that equivalent modelling approach mentioned earlier is almost always adopted—with damping ratios, loss factors or another dissipation quantifier accounting alone for all forms of energy loss.

**Figure 1.** Acoustic guitar parts. (Source: adapted from original image due to William Crochot, distributed under CC BY-SA 3.0 license; Wikimedia Commons).

#### 2.1.1. Simplified 2 Degrees-of-Freedom Model

Now, suppose the tensions in guitar strings transmit a force of magnitude *f* through the bridge (see Figure 1), that can be represented as acting pointwise, up and down, in the centroid of top plate. Thus, the equations of motion representing the dynamic coupling with the air column DOF displayed in Figure 2 can be written as

$$\begin{cases} m\_t \ddot{\mathbf{x}}\_t + k\_t \mathbf{x}\_t - \Delta p A\_t = f\_\prime \\ m\_a \ddot{\mathbf{x}}\_a - \Delta p A\_a = 0. \end{cases} \tag{1}$$

**Figure 2.** Schematic of the simple 2 DOF model: (**a**) Cut plane view of a guitar resonant chamber 3D CAD model, constructed with shells, planes and lines for subsequent FE analysis. (**b**) Representation as a simple 2 DOF model, considering only the motion of top plate and air column.

Considering the adiabatic, linearized form, of the equations of acoustics, the pressure variation Δ*p* may be expressed in terms of the total volume *V* inside the resonant chamber, the sound velocity *c* and mass density *ρ* of the surrounding air, and of the (idealized) chamber volume increment Δ*V* = *Atxt* + *Aaxa* [14], as

$$
\Delta p = -\rho \mathbf{c}^2 \frac{A\_t \mathbf{x}\_t + A\_a \mathbf{x}\_a}{V}. \tag{2}
$$

Thus, substituting Equation (2) into (1) and writing it in matrix form, i.e., *Mx*¨ + *Kx* = *f* , yields

$$
\begin{bmatrix} m\_t & 0\\ 0 & m\_a \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{x}}\_t\\ \ddot{\mathbf{x}}\_a \end{bmatrix} + \begin{bmatrix} k\_t + \frac{\rho c^2 A\_t^2}{V} & \frac{\rho c^2 A\_t A\_a}{V} \\\\ \frac{\rho c^2 A\_t A\_d}{V} & \frac{\rho c^2 A\_d^2}{V} \end{bmatrix} \begin{bmatrix} \mathbf{x}\_t\\ \mathbf{x}\_a \end{bmatrix} = \begin{bmatrix} f\\ 0 \end{bmatrix}. \tag{3}
$$

Quantities appearing on the left-hand side of Equation (3) could be viewed separately into two classes: those described in Table 2 depend on geometric and/or material properties of guitar components, and thus shall be treated as adjustable parameters of the equivalent system, while the remaining ones (speed of sound and mass density) are physical properties of the surrounding air, and so it is more reasonable to consider them as system-independent constants. On the other hand, the driving force *f* would be truly impossible to obtain, both experimentally or analytically, because this excitation motion is rather imposed by an effectively distributed conjugate on the interfacing area between the top plate and the bridge; i.e., its much higher stiffness impels it to rotate when forced to and fro by the tension of guitar strings, acting some distance above the plate. Theoretically, once the plate is adequately approximated (to some level of accuracy) as a clamped or hinged plate [11] an equivalent pivot distance could be considered, via static equilibrium of moments, in order to compare the driving force over the top plate with the tension in guitar strings. This approach has not been explored up to this moment, however, to the authors' knowledge. Nevertheless, analytical values for modal parameters can be readily extracted from this matrix equation, and in this way it has been extensively used for model comparison [4,7,9], structural modifications prediction based on eigenvalues sensitivity [9,10] or frequency domain model updating [15].

Unlike the mentioned references, the previous equations were written in undamped form. While this choice was made primarily to simplify the presentation, it has a more profound justification stemming from the authors' previous works [16] that is: given the levels of uncertainty in these simplified models, the use of more complex, non-proportional damping formulations would not cause any perceivable alteration in measurable outputs. Thus, since a proportionally damped analytical system produces the same eigenproperties as its undamped version, and taking into consideration the prominence of modal parameters over the forced response for the envisaged applications, damping factors can be confidently considered here as fine tuning adjustments. Among the myriad of available methods to extract damping coefficients (along with mode shapes and natural frequencies) from measured responses, Subspace-Stochastic Identification (SSI) has proven to be a very convenient, output-only, alternative [16].

#### 2.1.2. Simplified 3 and 4 Degrees-of-Freedom Models

From the definitions presented above for 2 DOF, modelling extensions can be readily devised. Figure 3 shows schematic representations for two of these possible expansions, resulting in 3 and 4 DOF models, respectively. The most obvious approach then is to consider the flexibility of back plate in a symmetrical manner as was made for the top plate. This in turn adds up three more free parameters: the mass, the stiffness and one more modal damping coefficient, and the system can now vibrate with three linearly independent modes. Thus, following [10], the resultant displacements vector and the corresponding mass and stiffness matrices for the 3 DOF model are given as

$$\begin{aligned} \mathbf{x} &= \begin{bmatrix} \mathbf{x}\_t \\ \mathbf{x}\_a \\ \mathbf{x}\_b \end{bmatrix}, \mathbf{M} = \begin{bmatrix} m\_t & 0 & 0 \\ 0 & m\_a & 0 \\ 0 & 0 & m\_b \end{bmatrix}, \\ \mathbf{K} &= \begin{bmatrix} k\_l + \frac{\rho c^2 A\_l^2}{V} & \frac{\rho c^2 A\_l A\_a}{V} & \frac{\rho c^2 A\_l A\_b}{V} \\ \frac{\rho c^2 A\_l A\_a}{V} & \frac{\rho c^2 A\_a^2}{V} & \frac{\rho c^2 A\_a A\_b}{V} \\ \frac{\rho c^2 A\_l A\_b}{V} & \frac{\rho c^2 A\_a A\_b}{V} & k\_b + \frac{\rho c^2 A\_b^2}{V} \end{bmatrix}. \end{aligned} \tag{4}$$
 
$$\begin{bmatrix} \frac{\rho c^2 A\_l A\_b}{V} & \frac{\rho c^2 A\_l A\_b}{V} & k\_l + \frac{\rho c^2 A\_l^2}{V} \\ \frac{\rho c^2 A\_l}{V} & \frac{\rho c^2 A\_l}{V} & \frac{\rho c^2 A\_l}{V} & \frac{\rho c^2 A\_l}{V} \\ \frac{\rho c^2 A\_l}{V} & k\_l & \frac{\rho c^2}{V} & k\_l \\ \frac{\rho c^2 A\_l}{V} & \frac{\rho c^2 A\_l}{V} & \frac{\rho c^2 A\_l}{V} & k\_l \end{bmatrix} \tag{5}$$

**Figure 3.** Schematic of 3 and 4 DOF simplified models; (**a**) 3 DOF: allowing for motions on top and back plates, and air column through the hole. (**b**) 4 DOF: basically allowing for the same motions of the 3 DOF model, plus one for the entire structure (free condition).

The next modeling extension, considering mass and motion for the ribs, is a little bit trickier and thus needs deeper attention. First of all, the motion allowed for the 'ribs' in fact means that the whole guitar is free to move. Thus, this is not exactly about freeing a motion in the same sense as it was done for the back plate, but rather allowing a fully free boundary condition. While the governing equations of the 4 DOF model in [11] are not explicitly presented in matrix form, the matrix form was derived and the model is here written in matrix form as

$$\mathbf{x} = \begin{bmatrix} x\_t \\ x\_a \\ x\_b \\ x\_r \end{bmatrix}, \mathbf{M} = \begin{bmatrix} m\_t & 0 & 0 & 0 \\ 0 & m\_a & 0 & 0 \\ 0 & 0 & m\_b & 0 \\ 0 & 0 & 0 & m\_r \end{bmatrix},$$

$$\mathbf{K} = \begin{bmatrix} k\_l + \frac{\rho c^2 A\_l^2}{V} & \frac{\rho c^2 A\_l A\_a}{V} & \frac{\rho c^2 A\_l A\_b}{V} & -k\_l - \frac{\rho c^2 A\_l^2}{V} + \frac{\rho c^2 A\_l A\_b}{V} - \frac{\rho c^2 A\_l A\_a}{V} \\\\ \frac{\rho c^2 A\_l A\_b}{V} & \frac{\rho c^2 A\_l A\_b}{V} & k\_b + \frac{\rho c^2 A\_b^2}{V} & k\_b + \frac{\rho c^2 A\_b^2}{V} - \frac{\rho c^2 A\_l A\_b}{V} - \frac{\rho c^2 A\_l A\_b}{V} \\\\ -k\_l & 0 & k\_b & k\_l + k\_b \end{bmatrix}. \tag{5}$$

$$\frac{\rho c^2 A\_l A\_b}{V} = \frac{\rho c^2 A\_l A\_b}{V} - \frac{\rho c^2 A\_l A\_b}{V} \qquad \frac{\rho c^2 A\_b^2}{V} - \frac{\rho c^2 A\_b A\_b}{V} - \frac{\rho c^2 A\_l A\_b}{V}$$

#### **3. Finite Element Modelling of Acoustic Guitar Vibroacoustics**

Many factors have contributed to the popularity of FE modelling in the last two or three decades, but most of all, possibly, the greater availability of computational resources. Thus, it is not surprising that nowadays, the bulk of academic research in vibro-acoustics of musical instruments makes use of it to a greater or lesser extent. Brooke [17] was a pioneering work on the application of FE to guitar coupled vibro-acoustics. The simplified models were employed alongside the FE modelling of the top plate and boundary element (BE) model of surrounding air, thus achieving an organic coupling between lumped parameters and domain discretization approaches. Elejabarrieta, Ezcurra and Santamaría [18,19] performed fluid-structure interaction FE analyses, both for fluid and structural domains. The fluid domain was considered only in a small air column, in a similar fashion as was done for the lumped models discussed earlier. Neck and headstock participations were disregarded. 8474 brick elements were used to model the structure, while the mesh on the fluid domain was adjusted to give sufficiently accurate results up to 1 kHz. Comparison of coupled and uncoupled mode shapes confirmed the physical insights brought by lumped models and the importance of Helmholtz frequency to the lower range coupled dynamics. Chaigne and collaborators [20–22] successfully employed FE and finite differences in a full non-linear transient simulation, where considerable effort is put into precise definition of purely mathematical aspects related to the problem's well-posedness and stability of the proposed solution procedure. The fictitious domain method was used to avoid what they called an 'ectoplasm', referring to the rather artificial definition of boundaries in the fluid domain. Structural flexibility of the resonant chamber is considered only within the top plate, disregarding fan bracing. Both these simplifications may leave aside important aspects of top plate dynamic behaviour. It is interesting to note that, despite the clear relevance of such a complete description, these important features are especially focused on by the lumped-parameter models presented above. Other instruments, closely related to the acoustic guitar, were successfully modelled via FE analysis as well, and we pinpoint here references [23], on the Brazilian viola caipira, and [24], on the Colombian bandola. Both present detailed equations and solution procedures for the full coupled FE analysis and solve equations using ANSYS commercial software, but in [24] a convenient symmetrized form for the coupling matrices is introduced and an effective length is considered for the air column [18,19].

#### **4. Results and Discussion**

This section is intended to show some results obtained by the authors with both kinds of modelling techniques discussed up to this point. Updating/identification of the 3 DOF model parameters was performed by Löw [15], following [10], but using some heuristic optimization algorithms instead to minimize the sum of squared differences between analytical and experimental FRFs. Masses were measured with 0.1 g accuracy, and stiffness parameters were previously estimated by static force/displacement testing, as suggested by [6,7,9]. Figure 4 below compares the updated theoretical model with experimental data; a modified particle swarm algorithm was used in error minimization. Experimental data was acquired on a low quality, large-scale manufactured instrument.

For the sake of comparison, some more recent results [16] obtained for a Giannini GWNE15, a Brazilian mid quality, partially hand-made guitar, are presented below, this time parametrically modelled in ANSYS commercial software. The mesh used 17,001 shell elements, considering wood plates orthotropic properties. Fan bracing, struts and neck were modelled with shell elements as well. Figure 5 below shows the FE mesh, and Figure 6 shows structural (in vacuo) modes of top and back plates. Table 3 compares experimental and analytical natural frequencies for both models showed earlier. Bigger discrepancies in FE model are primarily due to not considering the air-structure coupling, so it is difficult to establish mode matching. Even so, it is clear that the greater number of available vibration modes can be used to better represent effects associated with higher frequency range dynamics, such as the timber of instruments, and radiation efficiency.

**Figure 4.** 3 DOF model analytical vs. experimental FRFs: (**a**) top plate; (**b**) forced on top plate and measured on the back plate.


**Table 3.** Natural frequencies of the classic guitar with the updated 3DOF and FE models: analytical/experimental results (relative deviation).

**Figure 5.** (**a**) Giannini catalog photo for model GWNE15; (**b**) 3D CAD model; and (**c**) resultant mesh for numerical analysis of Giannini guitar.

The simplicity and low computational burden of the lumped parameter 3D model permitted an easy application of an optimization routine to curve-fit structural parameters to reproduce low-frequency range dynamic behaviour with good agreement. The phase angles were not taken into account in the error function and some acceptable agreement was verified as well, thus bringing reliability to the obtained parameters values. On the other hand, a numerical approximate solution of partial differential equations makes it easy to visualize mode shapes, obtain dozens of natural frequencies, maybe more than one hundred if one needs and has sufficient computational resources, accounting for global and local contributions as well. It should be mentioned, however, that FE results are largely affected by a significant number of purely numerical issues, such as element distortion (close to singular, or even negative jacobian determinant), and spurious modes due to artificial stiffness in FE formulation.

**Figure 6.** Undamped structural modes: (**a**) in phase top plate mode (75.9 Hz); (**b**) anti phase top plate mode (110.2 Hz); (**c**) first dipole top plate mode (169.3 Hz) and (**d**) quadrupole back plate mode (174.1 Hz).

#### **5. Conclusions**

The use of FE modeling is unquestionably advantageous for systems with a high degree of complexity, coupling diverse physical phenomena, where interesting applications for numerical analysis of musical instruments have been presented. While many of the FE codes are very user-friendly and non-commercial versions may be found, they are not accessible to the majority of the luthiers and the process is far too complex for them. This article describes the classical lumped mass approach and the widely used FE analysis for modelling and simulation of the acoustic guitar, highlighting the advantages and disadvantages of each. There is no undisputed standard for coupled vibroacoustic analysis, especially among users and manufacturers of acoustic guitars.

**Author Contributions:** Conceptualization, all; methodology, all; software, A.M.L. and H.M.G.; validation, all; formal analysis, all; investigation, all; writing—original draft preparation, A.M.L.; writing—review and editing, A.M.L. and C.M.A.V.; funding acquisition, H.M.G. and C.M.A.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** The first two authors gratefully acknowledge the support by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, CAPES, Brazil. The third author gratefully acknowledges the support granted by the Foundation for Science and Technology (FCT) of Portugal, in the scope of the project of the Research Unit on Materials, Energy and Environment for Sustainability (proMetheus), Ref. UID/05975/2020, financed by national funds through the FCT/MCTES.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

