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

1 Introduction

Classical guitars are remarkable musical instruments that have been refined through the expert craftsmanship of luthiers. Historical examples, such as “La Leona", crafted by Antonio de Torres in 1856, continue to shape the craft of guitar making, serving as foundational models and benchmarks [1]. Scientific investigation and development of methodologies to comprehend the functioning of guitars and to define the characteristics of a “good" guitar remain active and evolving research topics. Recent studies [2–4] indicate that evaluation criteria based on acoustic vibration characteristics can be used for this purpose allowing a more scientific approach compared to traditional methods such as tap-tuning.

The sound profile of a guitar is highly distinctive, making it nearly impossible to find two guitars with identical tonal characteristics. The properties of the materials and the geometric design are the primary factors that contribute to this uniqueness. Even with identical geometric designs and closely similar wood species, achieving the same sound profile between two guitars is not guaranteed. However, recent studies [5, 6] on guitar geometry reveal that the guitar geometric design has a more profound impact on sound quality compared to material properties variation.

Global warming has already affected tree habitats [7], causing tone-wood species to become increasingly scarce and expensive. Therefore, replicating the sound profile of an existing guitar model to achieve comparable acoustic characteristics, even within the constraints of available material resources, emerges as an interesting area of inquiry in instrument design. This study offers a robust scientific framework to tackle this challenge, leveraging parametric reduced-order finite element modeling integrated with optimization techniques.

In this study, we focus on the geometry of the guitar’s soundboard, encompassing its shape outline, sound hole, length, and thickness, while excluding any bracing patterns. Utilizing the framework presented in [8], a fully geometrically parameterized virtual soundboard model was employed in numerical optimization to align the vibrational modal parameters – eigenfrequencies and eigenmodes – of two soundboards. This approach allows for precise adjustment to the soundboard geometry, ensuring that the design can be optimized according to the specific vibration characteristics. However, it is acknowledged that vibrational modal parameters alone do not fully capture the acoustic radiation behavior or the perceptual qualities of the instrument. Nevertheless, these modal parameters offer a practical and widely accepted foundation for optimization.

In References [8, 9], geometric shape optimization of a complete guitar soundboard was performed by calculating modal parameters using the commercial finite element software Abaqus [10]. Building upon this work, the present study proposes a new and comprehensive methodology for guitar soundboard shape optimization that goes significantly beyond existing approaches. This methodology uniquely integrates several well-established techniques into a unified framework, including affine representation of finite element modeling, mesh morphing for geometric parameterization, sensitivity-based parameter selection, and Krylov-based moment matching for parametric model order reduction (PMOR).

At the core of this approach is the implementation of a parameterized flat shell element model and the symbolic formulation of an affine parameter-dependent form in Matlab [11], specifically designed to reduce computational power and memory requirements. Due to the limited availability of computationally efficient shell formulations suitable for parametric applications, a new formulation was developed and implemented to meet the demands of this study. In particular, the element formulation was carefully chosen to ensure low memory demands during the affine decomposition process. Additionally, expressing the system matrices in affine form required detailed analytical treatment and the adoption of appropriate strategies to ensure compatibility with parametric model order reduction. Therefore, extensive development efforts were undertaken to construct an accurate and efficient parametric shell model with affine dependency, capable of representing thin and flat-structure behavior. This novel formulation supports the effective application of PMOR, resulting in a substantial reduction in computational cost and time compared to commercial software such as Abaqus.

These advancements are particularly crucial for shape optimization procedures that require multiple and intensive model evaluations. Furthermore, the proposed framework is adaptable to any parameterized, geometrically complex guitar design and can be further extended to include material parameterization. To the best of our knowledge, this work represents one of the first applications of a fully integrated affine finite element and PMOR-based optimization framework for geometric shape optimization in the field of computational luthiery. The unified incorporation of these techniques enables a high-fidelity, resource-efficient approach to soundboard design, marking a significant methodological advancement in the field.

The structure of this paper is as follows. First, the theoretical foundation of the finite element method (FEM) for flat shell elements is presented, followed by the description of the newly implemented geometrically parameterized flat shell element approach. Next, the process of soundboard design parameterization is detailed, incorporating mesh morphing techniques. The PMOR of the implemented flat shell element model is then reviewed. The optimization methodology employed is subsequently discussed. Finally, the verification of the proposed methodology is provided, along with an example application demonstrating the geometric optimization of a soundboard.

2 Reduced order finite element model with geometric parameter dependency

Parametric modeling is extensively utilized in many engineering applications, facilitating parameter studies and optimization analyses for model development. For complex mechanical systems such as guitars, utilizing a parameterized FE model offers a practical and effective framework for advanced designs [12, 13].

In this study, the system is parameterized with vector p ∈ ℝ^d representing the geometric configuration of the guitar’s soundboard. For the spatially discretized and parameterized FE model, the second-order linear differential equation of motion for the undamped structural system with the parametric system mass matrix M(p)∈ℝ^N × N and stiffness matrix K(p)∈ℝ^N × N is given as

$$\begin{array}{c}\hfill \mathit{M}(\mathit{p})\mathit{q}+\mathit{K}(\mathit{p})\mathit{q}=\mathit{f}(\mathit{p})\end{array}$$(1)

with the vector of external forces f ( p )∈ℝ ^N , and the vector q  ∈ ℝ ^N contains the N nodal degrees of freedom of the system.

In the following, the theoretical background regarding the isoparametric finite shell elements modeling and its extension with parametrization and PMOR for dynamical, second-order systems via an interpolatory method are described.

2.1 Shell element formulation

The FEM is a versatile technique for addressing complex mechanical system problems that cannot be solved analytically. With sufficiently fine discretization, FEM can provide a highly accurate approximation of the correct solution. The fundamental derivation and concepts of the FEM can be found in [14, 15].

In this study, the guitar soundboard is discretized using four-node quadrilateral flat shell elements, which offer enhanced computational efficiency. This shell element-based discretization results in fewer degrees of freedom compared to solid brick elements, which typically demand multiple element layers for accurate bending behavior representation. Since the primary soundproduction of a guitar soundboard originates from its bending modes, the flat shell element-based finite element discretization is well suited for this application. The used element formulation accurately captures both bending and membrane actions, making it ideal for thin and flat structures like soundboards. The Reissner-Mindlin shell theory [16], is employed to formulate the shell element. Bilinear shape functions are used for interpolation, following the isoparametric concept.

A bounded domain is considered representing a flat shell mid-surface, which is discretized into n _e finite elements. The finite element solution, indexed with h, of the displacement and rotation vector field d for a four-node quadrilateral flat shell element, is expressed as

$$\begin{array}{c}\hfill \mathit{d}\approx \underset{{\mathit{d}}_{\mathit{h}}}{\underbrace{\left[\begin{array}{c}u\\ v\\ w\\ {\theta }_x\\ {\theta }_y\\ {\theta }_z\end{array}\right]}}=\sum _{i=1}^{n_p}\underset{\widetilde{{\mathit{N}}_{\mathit{i}}}}{\underbrace{\left[\begin{array}{cccccc}{N}_i& 0& 0& 0& 0& 0\\ 0& N_i& 0& 0& 0& 0\\ 0& 0& N_i& 0& 0& 0\\ 0& 0& 0& N_i& 0& 0\\ 0& 0& 0& 0& N_i& 0\\ 0& 0& 0& 0& 0& 0\\ \end{array}\right]}}\underset{{\mathit{q}}_{\mathit{i}}}{\underbrace{\left[\begin{array}{c}{u}_i\\ v_i\\ w_i\\ {\theta }_{\mathit{xi}}\\ {\theta }_{\mathit{yi}}\\ {\theta }_{\mathit{zi}}\end{array}\right]}},\end{array}$$(2)

where n _p  = 4, is the number of nodes, N _i denotes the bilinear shape functions associated with node i and q _i is the vector of nodal degrees of freedom at node i. In this context, u, v and w are the displacements in x, y, and z direction, respectively, while θ _x , θ _y , and θ _z represent the rotations about the corresponding axes. In the kinematic definitions of flat shell elements, the contribution of the drilling degree of freedom (θ _z ) is excluded. As a result, the bilinear shape function matrix, $\widetilde{{\mathit{N}}_{\mathit{i}}}\in {\mathbb{R}}^{6\times 6}$ becomes singular, requiring careful attention during stiffness calculations.

The element mass matrix m ^e  ∈ ℝ^24 × 24 for the e-th element is expressed as

$$\begin{array}{c}\hfill {\mathit{m}}^e={\int }_{-1}^1{\int }_{-1}^1{\mathit{N}}^{\top }\mathit{\rho }\mathit{N}det(\tilde{\mathit{J}})\mathrm{d}\xi \mathrm{d}\eta ,\end{array}$$(3)

where

$$\begin{array}{c}\hfill \mathit{\rho }=\left[\begin{array}{ccc}\varrho h{\mathit{I}}_{3\times 3}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \frac{1}{12}\varrho h^3{\mathit{I}}_{2\times 2}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \frac{1}{12}f\varrho h^3\end{array}\right]\in {\mathbb{R}}^{6\times 6}·\end{array}$$(4)

Here 𝜚, h, and I represent the mass density, element thickness and identity matrices, respectively. The shear correction factor (drill factor) f is taken as 5/6 in this formulation [10]. Following equation (2), the shape function matrix N  ∈ ℝ^6 × 24 consists of four concatenated bilinear shape function matrices, $\widetilde{{\mathit{N}}_{\mathit{i}}}$, which interpolate the 24 nodal degrees of freedom into the rotational and translational element displacements. The Jacobian matrix $\tilde{\mathit{J}}\in {\mathbb{R}}^{2\times 2}$ maps the isoparametric ξ–η coordinates to the global Cartesian x–y coordinates.

When formulating the stiffness matrix for a flat shell element, the element is treated as a combination of a membrane element-based on in-plane deformation and a plate element-based on out-of-plane bending deformation [17, 18]. This formulation enables the final element stiffness to capture membrane, bending, and transverse shear stiffness components. The stiffness matrix of the flat shell element, ${\stackrel{\mathbf{~}}{\mathit{k}}}^e\in {\mathbb{R}}^{24\times 24}$ is expressed as

$$\begin{array}{c}\hfill {\tilde{\mathit{k}}}^e={\int }_{-1}^1{\int }_{-1}^1{\tilde{\mathit{B}}}^{\top }\mathit{E}\tilde{\mathit{B}}det(\tilde{\mathit{J}})\mathrm{d}\xi \mathrm{d}\eta ,\end{array}$$(5)

computed as the summation of the stiffness contributions from both the membrane and plate elements. Here, $\tilde{\mathit{B}}\in {\mathbb{R}}^{8\times 24}$ represents the strain-displacement matrix, which incorporates contributions from both the membrane and plate elements. Additionally, $\mathit{E}\in {\mathbb{R}}^{8\times 8}$ denotes the orthotropic elastic material matrix under a linear plane stress condition as follows. Wood is an orthotropic material with unique and independent properties along three mutually perpendicular directions, its mechanical characteristics are defined by Young’s moduli E_L, E_R, and E_T (longitudinal, radial, and tangential directions), Poisson’s ratios ν_LR, ν_LT, and ν_RT, and shear moduli G_LR, G_LT, and G_RT. The material matrix reads

$$\begin{array}{c}\hfill \mathit{E}=\left[\begin{array}{ccc}{\mathit{E}}_m& 0& 0\\ 0& {\mathit{E}}_b& 0\\ 0& 0& {\mathit{E}}_s\end{array}\right],\end{array}$$(6)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{E}}_m& =h{\left[\begin{array}{ccc}\frac{1}{E_L}& -\frac{{\nu }_{\mathit{LR}}}{E_L}& 0\\ -\frac{{\nu }_{\mathit{LR}}}{E_L}& \frac{1}{E_R}& 0\\ 0& 0& \frac{1}{G_{\mathit{LR}}}\end{array}\right]}^{-1},\hfill \\ \hfill {\mathit{E}}_b& =\frac{h^2}{12}{\mathit{E}}_m\text{and}{\mathit{E}}_s=hf{\left[\begin{array}{cc}\frac{1}{G_{\mathit{LT}}}& 0\\ 0& \frac{1}{G_{\mathit{RT}}}\end{array}\right]}^{-1}·\hfill \end{array}\end{array}$$(7)

Here, ${\mathit{E}}_m\in {\mathbb{R}}^{3\times 3}$ represents the contribution of plane stress from the membrane element, while ${\mathit{E}}_b\in {\mathbb{R}}^{3\times 3}$ and ${\mathit{E}}_s\in {\mathbb{R}}^{2\times 2}$ describe the contributions of bending stress and transverse shear stress from the plate element, respectively. Due to the orthotropic nature of the material, the element stiffness matrix is initially computed in the local coordinate frame. Subsequently, using an appropriate transformation matrix that maps the material orientation to the global coordinate frame, the final element stiffness matrix is obtained. Details on the computation of the transformation matrix can be found in [19].

The contributions of the membrane and plate elements to the overall flat shell element stiffness matrix are represented by the following sub-matrix corresponding to element node i ∈ {1, 2, 3, 4}

$$\begin{array}{c}\hfill {({\stackrel{\mathbf{~}}{\mathit{k}}}^e)}^{i,i}=\left[\begin{array}{ccc}{({\tilde{\mathit{k}}}_m^e)}^{i,i}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {({\tilde{\mathit{k}}}_b^e)}^{i,i}+{({\tilde{\mathit{k}}}_s^e)}^{i,i}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \end{array}\right]\in {\mathbb{R}}^{6\times 6},\end{array}$$(8)

where ${({\tilde{\mathit{k}}}_m^e)}^{i,i}\in {\mathbb{R}}^{2\times 2}$ is the membrane element stiffness contribution and ${({\tilde{\mathit{k}}}_b^e)}^{i,i}+{({\tilde{\mathit{k}}}_s^e)}^{i,i}\in {\mathbb{R}}^{3\times 3}$ is the plate element stiffness contribution.

Gauss-quadrature numerical integration is typically recommended for evaluating equations (3) and (5). For equation (5) both full integration, using four Gaussian quadrature points, and reduced integration, employing a single Gaussian quadrature point, are viable options [14]. Each method, however, presents inherent challenges. Full integration can result in overly stiff behavior and locking phenomena. These issues are avoided using the Enhanced Assumed Strain (EAS) method [20] and the Assumed Natural Strain (ANS) method [16, 21], which address membrane locking and transverse shear locking, respectively. Despite this, full integration demands high computational resources, making it computationally expensive. In contrast, reduced integration may introduce zero energy modes, which are resolved using hourglass stabilization techniques [22]. Furthermore, the reduced integration approach has demonstrated superior efficiency for parameterized modeling [23], primarily due to its lower computational demands and reduced memory requirements. Accordingly, reduced integration-based element stiffness formulations are selected for implementation in this study.

To address matrix singularity issues in the stiffness matrix described in equation (5), an artificial drilling stiffness ${\mathit{k}}_d^e\in {\mathbb{R}}^{24\times 24}$ is often included. This artificial stiffness is applied to the matrix entries corresponding to the drilling-degrees of freedom, θ _{z i} . In this study, the implementation of the drilling stiffness follows the approach outlined in [24] and employs reduced integration. As a result, the final element stiffness matrix is formulated as

$$\begin{array}{c}\hfill {\mathit{k}}^e={\tilde{\mathit{k}}}^e+{\mathit{k}}_d^e·\end{array}$$(9)

Finally, assembling the element matrices yield the system matrices as

$$\begin{array}{c}\hfill \mathit{M}=\bigcup _{e=1}^{n_e}{\mathit{m}}^e\in {\mathbb{R}}^{N\times N}\text{and}\mathit{K}=\bigcup _{e=1}^{n_e}{\mathit{k}}^e\in {\mathbb{R}}^{N\times N}·\end{array}$$(10)

2.2 Shell element parameterization

In this study, the aforementioned general shell element formulation is extended by integrating element-level parameterization. This approach offers a distinct advantage because it confines finite element pre-processing and model order reduction to the offline phase. Consequently, these computationally intensive steps are not repeated during each optimization iteration in the online phase, leading to substantial savings in computational time.

A geometric parameterization of the shell elements introduced in [23] is applied. In this approach, the geometry of the shell element is defined by seven parameters that represent the coordinates of the element nodes. This parameterization enables the representation of all types of shell element geometries in a planar case. The element parameter vector is given by

$$\begin{array}{c}\hfill {\mathit{p}}^{\mathit{e}}={\left[\begin{array}{ccccccc}a& b& c& d& e& g& h\end{array}\right]}^{\top }·\end{array}$$(11)

A visual representation of the element geometric parameters is given in Figure 1. The seventh parameter h, not shown in the figure, represents the thickness of the element. Here, node &#Xtextcircled;1 is taken as the origin of the two-dimensional x–y Cartesian plane. This parameterization allows to represent all kind of shell element geometries in planar mode.

Figure 1. Geometrically parameterized shell element.

With this element parameterization, equations (3) and (9) become parameter-dependent and can be rewritten in an affine parameter dependency form, resulting in

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{m}}^e({\mathit{p}}^{\mathit{e}})& ={\int }_{-1}^1{\int }_{-1}^1{\mathit{N}}^{\top }\mathit{\rho }({\mathit{p}}^{\mathit{e}})\mathit{N}det\left(\tilde{\mathit{J}}({\mathit{p}}^{\mathit{e}})\right)\mathrm{d}\xi \mathrm{d}\eta \hfill \\ \hfill & =\sum _{i=1}^I{\overline{w}}_i^e({\mathit{p}}^e){\mathit{m}}_i^e,\hfill \\ \hfill {\mathit{k}}^e({\mathit{p}}^{\mathit{e}})& ={\int }_{-1}^1{\int }_{-1}^1{\mathit{B}}^{\top }({\mathit{p}}^{\mathit{e}})\mathit{E}({\mathit{p}}^{\mathit{e}})\mathit{B}({\mathit{p}}^{\mathit{e}})det\left(\tilde{\mathit{J}}({\mathit{p}}^{\mathit{e}})\right)\mathrm{d}\xi \mathrm{d}\eta \hfill \\ \hfill & =\sum _{j=1}^J{\overline{w}}_j^e({\mathit{p}}^e){\mathit{k}}_j^e,\hfill \end{array}\end{array}$$(12)

where ${\overline{w}}_{i,j}^e({\mathit{p}}^e)$ are the i, j-th weight functions of the e-th element expressed in term of element parameter vector p ^e and m _i ^e , k _j ^e  ∈ ℝ^24 × 24 are constant matrices. These constant matrices are precomputed in the offline phase, while the weight functions are computed in the online phase. Here, B represents the strain-displacement matrix, which integrates the membrane and plate formulations along with the formulation for the drilling degrees of freedom. The exact derivation of the parameter dependency is omitted here, but an example of a detailed derivation of parameter dependency for a solid brick element can be found in [25].

The implemented hourglass-stabilized reduced integration in stiffness calculations results in a comparatively lower number of $\tilde{\mathit{J}}({\mathit{p}}^{\mathit{e}})$ summations in the affine representation, thereby reducing memory demands during parametric calculations.

2.3 Guitar soundboard parameterization

In this study, the geometric parameterization of the soundboard follows the approach outlined in [26, 27], where the following rational polynomials in polar coordinates are utilized to define the outline of the soundboard

$$\begin{array}{c}\hfill r(\theta )=L\frac{c_1+c_3\theta +c_5{\theta }^2+c_7{\theta }^3+c_9{\theta }^4}{1+c_2\theta +c_4{\theta }^2+c_6{\theta }^3+c_8{\theta }^4}·\end{array}$$(13)

The radius r(θ) defines the contour as a function of θ. The parameters c _i , i = {1, 2, …, 9} determine the shape of the outline, while L represents the scale length. In addition, the plate thickness h and the sound hole radius R are also treated as independent parameters. These parameters provide sensible trade-off between the range of producible soundboard shapes and the number of necessary independent parameters. Collectively, this results in d = 12 independent geometric parameters, represented by

$$\begin{array}{c}\hfill \mathit{p}={\left[\begin{array}{cccccccc}{c}_1& c_2& \dots & c_8& c_9& L& h& R\end{array}\right]}^{\top }\in {\mathbb{R}}^d·\end{array}$$(14)

To generate various soundboard samples for the evaluation during shape optimization, the parameter values in vector p are utilized as global parameters within the defined parameter space. However, these global parameters must be mapped to the element level during system matrix evaluation. Achieving this requires an appropriate mesh topology relationship between the soundboard’s geometric shape and the element geometry. Moreover, the total number of degrees of freedom varies across different parametric configurations, indicating that the model is not inherently parametric at the FE level. Therefore, mesh preservation becomes crucial. Also, preserving identical element counts and node connectivity across all sampled soundboard models, aligned with a predefined reference soundboard, eliminates the need for repeated discretization during each model evaluation in the optimization process, thereby significantly enhancing computational efficiency. Additionally, this mesh preservation provides comprehensive comparison of soundboard mode shapes across different parameter settings.

In this study, we adopt the mesh preservation approach proposed in [9], which employs a barycentric mapping-based mesh morphing technique [28, 29]. This method enables to achieve a parametric model at the FE level, ensuring a consistent number of degrees of freedom across all soundboard models. With this mesh morphing technique and global parameterization, the local element parameters p ^e can be expressed in terms of the global parameter vector p . As a result, this formulation allows the element matrices to be represented as functions of the global parameter vector. For instance, the parameterized element stiffness matrix in equation (12) can be expressed in affine form as

$$\begin{array}{c}\hfill {\mathit{k}}^e(\mathit{p})=\sum _{j=1}^J{\overline{w}}_j^e({\mathit{p}}^e(\mathit{p})){\mathit{k}}_j^e=\sum _{j=1}^J{w}_j^e(\mathit{p}){\mathit{k}}_j^e·\end{array}$$(15)

Here, w _j ^e ( p ) is the j-th weight functions of e-th element expressed in terms of global parameter vector p . However, significant geometric variations can cause excessively distorted mesh distributions, yielding unusable finite element models. Therefore, moderate geometric variations were considered to avoid these potential concerns. With these relations, the parameterized global stiffness matrix can be assembled as

$$\begin{array}{c}\hfill \mathit{K}(\mathit{p})=\bigcup _{e=1}^{n_e}{\mathit{k}}^e(\mathit{p})·\end{array}$$(16)

2.4 Parametric model order reduction

For a complex system like the guitar, to obtain a satisfactory approximation of the solution in FE modeling, a fine spatial discretization is necessary. However, this increases the number of degrees of freedom N, leading to long computational time. This is crucial, particularly when multiple model evaluations are required for numerical optimization. This drawback can be overcome using PMOR techniques. The core concept of this approach is to develop an approximate reduced-order model that retains essential model information in a defined parameter space while being significantly more computationally efficient than the original model.

In this study, PMOR with interpolatory projection-based model [30] is employed. Using this approach, equation (1) is transformed into an input-output system of equations as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{M}\mathbf{(}\mathit{p}\mathbf{)}\mathit{q}+\mathit{K}(\mathit{p})\mathit{q}& =\mathit{D}\mathbf{(}\mathit{p}\mathbf{)}\mathit{u}\text{and}\hfill \\ \hfill \mathit{y}& =\mathit{C}\mathbf{(}\mathit{p}\mathbf{)}\mathit{q},\hfill \end{array}\end{array}$$(17)

where f ( p )= D ( p )  u is substituted, introducing the system input vector u  ∈ ℝ ^l , which contains the system inputs, which are distributed across the degrees of freedom using the input matrix D ( p )∈ℝ ^N × l . The output vector y  ∈ ℝ ^k is determined by the output matrix C ( p )∈ℝ ^k × N , which extracts relevant degrees of freedom from the vector q . The transfer function between the system inputs and outputs is expressed as

$$\begin{array}{c}\hfill \mathit{H}(s,\mathit{p})=\mathit{C}(\mathit{p}){(s^2\mathit{M}(\mathit{p})+\mathit{K}(\mathit{p}))}^{-1}\mathit{D}(\mathit{p}),\end{array}$$(18)

in Laplace domain, where s ∈ ℂ. In the case of s = i ω, the transfer function H (i ω,  p )∈ℝ ^k × l becomes the frequency response function. The full order model vector q is approximated as a reduced-order model vector q _r  ∈ ℝ ⁿ through a projection matrix V  ∈ ℝ ^N × n in lower-dimensional subspace,

$$\begin{array}{c}\hfill \mathit{q}\approx \mathit{V}{\mathit{q}}_r\end{array}$$(19)

with n <  < N. Substituting equation (19) into the full order model equation (17) and left-multiplying the transpose of another projection matrix W  ∈ ℝ ^N × n , the reduced-order model equation is obtained as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{{\mathit{M}}_r(\mathit{p})}{\underbrace{{\mathit{W}}^{\top }\mathit{M}(\mathit{p})\mathit{V}}}{\mathit{q}}_r+\underset{{\mathit{K}}_r(\mathit{p})}{\underbrace{{\mathit{W}}^{\top }\mathit{K}(\mathit{p})\mathit{V}}}{\mathit{q}}_r& =\underset{{\mathit{D}}_r(\mathit{p})}{\underbrace{{\mathit{W}}^{\top }\mathit{D}(\mathit{p})}}\mathit{u},\hfill \\ \hfill {\mathit{y}}_r& =\underset{{\mathit{C}}_r(\mathit{p})}{\underbrace{\mathit{C}(\mathit{p})\mathit{V}}}{\mathit{q}}_r,\hfill \end{array}\end{array}$$(20)

where matrices M _r ( p )∈ℝ ^n × n , K _r ( p )∈ℝ ^n × n , D _r ( p )∈ℝ ^n × l and C _r ( p )∈ℝ ^k × n represent the reduced-order mass, stiffness, input, and output matrices, respectively. According to [31] usually and also here W  =  V is chosen.

To determine an appropriate global projection matrix V , several local projection matrices ${\tilde{\mathit{V}}}_i({\sigma }_i,{\mathit{p}}_i)\in {\mathbb{R}}^{N\times n_i}$ were computed as

$$\begin{array}{c}\hfill \text{span}\left({\tilde{\mathit{V}}}_i({\sigma }_i,{\mathit{p}}_i)\right)=\text{span}\left({({\sigma }_i^2\mathit{M}({\mathit{p}}_i)+\mathit{K}({\mathit{p}}_i))}^{-1}\mathit{D}(\mathit{p})\right)\end{array}$$(21)

discrete parameter expansion points p _i and frequency expansion points σ _i . These projection matrices are constructed using moment matching based on the Krylov-subspace method, ensuring minimal projection error [31, 32]. The local projection matrices are then concatenated into

$$\begin{array}{c}\hfill \tilde{\mathit{V}}=[{\tilde{\mathit{V}}}_1({\sigma }_1,{\mathit{p}}_1),{\tilde{\mathit{V}}}_2({\sigma }_2,{\mathit{p}}_2),\dots ]\in {\mathbb{R}}^{N\times \tilde{n}},\end{array}$$(22)

to form global bases for the projection, where $\tilde{n}=\sum n_i>n$. To eliminate linearly dependent columns in $\stackrel{\mathbf{~}}{\mathit{V}}$, a truncated singular value decomposition (SVD) is applied, resulting in the reduced global projection matrix V . Applications of this PMOR approach can be found in [23, 25, 33].

3 Optimization

A geometrically optimized soundboard entails designing a completely new soundboard configuration that closely replicates the tonal characteristics of a predefined reference soundboard. In this work, as a measure for the tonal characteristics we use eigenfrequencies and eigenmodes of the soundboard. To evaluate the model parameter deviation between parameterized soundboard and reference soundboard, the objective function

$$\begin{array}{c}\hfill J(\mathit{p})=\sum _{m=1}^M[{\left(\frac{f_m^{\text{ref}}-f_m(\mathit{p})}{f_m^{\text{ref}}}\right)}^2+{\left(\frac{1-\sqrt{{\text{MAC}}_{\mathit{mm}}(\mathit{p})}}{\sqrt{{\text{MAC}}_{\mathit{mm}}(\mathit{p})}}\right)}^2]\end{array}$$(23)

is employed as proposed in [9] with the optimal solution

$$\begin{array}{c}\hfill {\mathit{p}}^{\ast }=arg\underset{\mathit{p}\in {\mathbb{R}}^d}{min}J(\mathit{p})·\end{array}$$(24)

This function is applied to the first M non-rigid body modes, where f _m ^ref denotes the reference soundboard’s frequency, and f _m ( p ) represents the parametric soundboard’s frequency for the m-th mode, expressed as a function of the global parameter vector p defined in equation (14). Furthermore, the Model Assurance Criterion (MAC) [34] matrix is used to compare the eigenmode shapes. The diagonal MAC values MAC _{m m} ( p ), compare the m-th eigenmodes between the reference soundboard and the parametric soundboard.

In this study, two soundboard designs (reference and initial) modeled with two different wood species are considered. The goal is to determine the optimal geometric design for the initial soundboard, SB^init, through shape optimization, such that its tonal characteristics closely resembles those of the reference soundboard, SB^ref.

Meanwhile, free boundary conditions are assumed for model analysis. The soundboard meshing is done in Abaqus, while all the other calculations are done in Matlab. Subsequently, the first M = 16 non-rigid body modes are considered for objective function evaluation. The minimization of objective J(p) is performed using Matlab’s fmincon function with the interior-point algorithm. The optimization process is considered as converged when the first-order optimality criterion ∥ΔJ∥< 10⁻¹⁰ is satisfied. Following [9], to choose suitable parameter configurations to be used as starting point for the optimizations, the objective functions has been evaluated for an extensive amount of parameter points in the parameter space. The 10 best-performing parameter points are then selected as initial values for the optimization algorithm. Subsequently, the final solution of equation (24) is identified as the one corresponding to the lowest J(p) value among the 10 optimizations performed. Nevertheless, despite exploring various parameter starting points, it cannot be ruled out that the presented solution represents a local minimum rather than a global one.

4 Verification

Following the model description, a verification is conducted to validate the parametric reduced-order shell element model and the geometric optimization model using a reference guitar soundboard (SB₁^ref) and a different guitar soundboard (SB₁^init). Here, both soundboards are modeled using red cedar wood with material parameters specified in Table 1. The goal is to identify the appropriate geometric design starting from SB₁^init that exhibits the tonal characteristics of SB₁^ref as similar as possible. Since the material properties remain consistent while the initial geometric designs differ between SB₁^ref and SB₁^init, the defined optimization process should yield the reference soundboard geometry as the optimal soundboard geometry, SB₁^opt, since it is the global minimum.

First, the implemented parameterized affine full-order FE shell model is validated against Abaqus S4R shell elements using 1000 soundboard samples, achieving an absolute average error of 0.16% across the first 30 non-rigid body mode frequencies, with an average MAC value of 0.989. Subsequently, using the same soundboard parameter sets, the implemented parameterized reduced-order FE model is tested against its full-order affine representation for the same 30 modes. In this context, the same spatial positions are selected as input and output such that D  =  C ^⊤ in equation (17). Three distinct locations are chosen on the soundboard, deliberately avoiding nodal points, and only the translational degree of freedom in the out-of-plane direction is considered for both input and output. The accuracy of the reduced-order FE model is confirmed with an absolute average error of 0.001% for the mode frequencies and an average MAC value of 0.998.

Then, the geometric optimization model is validated as follows. SB₁ ^ref is selected based on the outline contour of the Martin D28 acoustic guitar soundboard, while the Martin 0018 acoustic guitar soundboard outline contour is considered as the initial design, SB₁ ^init. To determine the corresponding geometric parameters of the two soundboards, a regression technique is applied using equation (13). Both soundboards are discretized with an identical element count of n _e  = 9966 and nodal connectivity, resulting in a total of N = 61068 degrees of freedom. As in the previous case, n _{i n p}  = 3 identical input-output points are selected. Here, n _f  = 30 uniformly distributed frequency expansion points are selected between 10 Hz and 500 Hz, as the eigenfrequencies of interest lie within this range, and n _p  = 50 parameter expansion points are chosen within the relevant parameter space for optimization. After computing the local projection matrices as defined in equation (21) and concatenating them into the global matrix as in equation (22), the matrix $\tilde{\mathit{V}}\in {\mathbb{R}}^{4500\times 61068}$ is formed with 4500 = n _{i n p}  ⋅ n _p  ⋅ n _f. A truncated SVD is then applied, retaining the left singular vectors associated with singular values greater than 4 ⋅ 10⁻¹⁴, resulting in the final global projection matrix V ∈ ℝ^{981 × 61068}. The reduced-order system thus yields n = 981 degrees of freedom, representing a significant reduction in system size. This enables considerably faster iterative eigenvalue analyses during optimization. The reduced model achieves an average computational time of 0.02 s for evaluating the first 16 non-rigid body modes in Matlab offering a 75-fold speed-up compared to full model analysis in Abaqus, thereby significantly reducing computational cost.

Figure 2 illustrates the comparative analysis of the first 16 eigenfrequencies and the corresponding MAC matrix for the eigenmodes of both SB₁^ref and SB₁^init. Notably, the mode shape distributions of the two soundboards appears to be very similar with the exception of modes 8 and 9, which are swapped. The average MAC value on the diagonal is ${\overline{\text{MAC}}}_{\mathit{mm}}=0.845$. These differences in modal parameters are attributed solely to the geometric dissimilarity between the two soundboards.

Figure 2. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).

The purpose of this verification process is to determine whether the final optimized design starting from SB₁^init, denoted as SB₁^opt obtained from p^*, aligns with the geometric parameters of SB₁^ref and produces similar tonal characteristics. An initial evaluation of 100000 soundboard samples is conducted using the geometric parameters of SB₁^init as nominal values to determine a suitable starting point. In this case, a wide parameter space is considered to investigate the existence of any local minimum, rather than directly targeting the known global minimum. Therefore, a large number of samples are evaluated to identify 10 best starting points for the optimization that represent this broad parameter space. This also enables verification of the uniqueness of the optimal solution under a widely distributed set of parameter points. At the end of the optimization, out of the 10 starting points, only four parameter sets yielded results close to the optimum, with parameter values nearly identical to those of SB₁^ref, which represents the known global minimum here. In these cases, approximately 230 modal evaluations were performed. All other starting points failed to produce a satisfactory distribution of modal parameters matching the reference. Table 2 shows the parameter values of the reference plate SB₁^ref, the initial parameter values of SB₁^init used as nominal values for the parameter space, the parameter values of the optimized plate SB₁^opt and the error in percentage between SB₁^opt and SB₁^ref.

Based on these parameter values, Figure 3 displays the comparative analysis of modal parameters after optimization. The optimized soundboard design has a good agreement with the reference modal parameter values. The average MAC value on the diagonal has improved to ${\overline{\text{MAC}}}_{\mathit{mm}}=0.998$.

Figure 3. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).

Figure 4 illustrates the relative error of the initial and optimized soundboard frequencies, compared to the reference soundboard frequencies. The cent is a unit of measurement used to quantify musical intervals based on a logarithmic scale. Typically, cents are used to express musical intervals in acoustics, assess intonation, or compare the audible difference between two frequencies. According to the observed data, the initial average relative error of the eigenfrequencies over the 16 modes is $\overline{|{\epsilon }^{\text{rel}}|}=17.20\%=327\text{cent}$, which is reduced to $\overline{|{\epsilon }^{\text{rel}}|}=0.127\%=2\text{cent}$ after optimization, indicating almost identical musical intervals across the entire frequency range.

Figure 4. Geometrically parameterized shell element.

5 Sensitivity analysis and parameter selection

This study considers a 12-dimensional independent global geometric parameter space for soundboard parameterization. The high dimensionality of this parameter vector increases computational complexity during the optimization procedure. Furthermore, among the 12 independent global geometric parameters in equation (14), certain parameters may have a negligible effect on the modal parameters. To address this, a sensitivity analysis is conducted to identify and exclude geometric parameters with minimal impact from subsequent evaluations, thereby further reducing computational cost during optimization.

The sensitivity analysis is performed using 50000 uniformly distributed soundboard samples generated via a quasi-random Sobol sequence [36], considering all 12 geometric parameters. For the correlation coefficient matrix, 12 geometric parameters and the first 16 nonzero eigenfrequencies of the samples are considered. The correlation coefficient of two random variables x and y is a measure of their linear dependence. If each variable has Λ scalar observations, then the correlation coefficient is defined as

$$\begin{array}{c}\hfill \mathrm{\Gamma }(x,y)=\frac{1}{\mathrm{\Lambda }-1}\sum _{i=1}^{\mathrm{\Lambda }}\left(\frac{x_i-{\mu }_x}{{\sigma }_x}\right)\left(\frac{y_i-{\mu }_y}{{\sigma }_y}\right)\end{array}$$(25)

with the mean μ and standard deviation σ [37].

To have a meaningful matrix of correlation coefficients, a mode shape tracking method is used instead of relying solely on the eigenmodes’ order of appearance. In this mode tracking process, eigenmode shape evaluations are conducted by comparing predefined eigenmodes (A, B, C, …) and sample soundboard’s eigenmodes using the MAC. Then, the selection and elimination of corresponding sample soundboard’s eigenfrequencies are done based on the heuristic minimum MAC value of 0.6. By effectively identifying and pairing eigenfrequencies with their respective mode shapes, this approach ensures a meaningful representation in the calculation of the correlation matrix, even when modes appear in a different order due to variations in the underlying parameters or configurations.

Figure 5, on the top, illustrates the matrix of the correlation coefficient between the 12 independent geometric parameters and the predefined 16 eigenfrequencies. On the bottom in Figure 5, the correlation coefficient values of corresponding geometric parameters ranked and scaled by the 2-norm for all frequencies are shown. The parameters L, h, c₈, and c₉ exhibit the highest sensitivity to the eigenmode frequencies. Additionally, several other parameters exhibit notable effect on eigenmodes. However, a reduced parameter space provides a more efficient outcome in the optimization process. Thus, by considering a reasonably broad number of parameters, the initial 12 independent parameters are reduced to 7, thereby significantly decreasing the computational cost of optimization. The new geometric parameter vector is

Figure 5. Correlation coefficient of 12 geometric parameters and predefined eigenmodes (top). Correlation coefficient ranking scaled by the 2-norm over predefined eigenmodes (bottom).

$$\begin{array}{c}\hfill \mathit{p}={\left[\begin{array}{ccccccc}{c}_1& c_5& c_7& c_8& c_9& L& h\end{array}\right]}^{\top }\in {\mathbb{R}}^7·\end{array}$$(26)

6 Geometric optimization

This section addresses the geometric shape optimization of a soundboard, where the reference and initial soundboard designs are identical but modeled with two different wood species. Using the verified shape optimization method, the goal is to determine the optimal parameter values that yield the appropriate geometric modifications from the initial soundboard design. The seven-dimensional parameter space defined in equation (26) is employed for this geometric optimization. All other parameters are kept fixed at the initial soundboard’s nominal values throughout the optimization process.

In this analysis, the Martin 0018 classical acoustic guitar’s soundboard is chosen as the reference soundboard (SB₂ ^ref) with red cedar selected as the wood species. White cedar was used for the other design, SB₂ ^init, which has an identical initial geometry to the reference soundboard and undergoes geometric shape optimization. The material parameters are also listed in Table 1.

Both soundboard designs are discretized with an identical element count of 10010, resulting in N = 61368 degrees of freedom. For the reduced-order model the same number of frequency and parameter expansion points, as well as the same input-output locations used in the verification, are considered here as well, resulting in $\tilde{\mathit{V}}\in {\mathbb{R}}^{4500\times 61368}$. By applying truncated SVD, retaining the left singular vectors corresponding to singular values greater than 4.5 ⋅ 10⁻¹⁴, the final global projection matrix V ∈ ℝ^{1369 × 61368} is obtained resulting a reduced-order model with n = 1369 degrees of freedom. This reduced-order model speeds up the eigenvalue analysis by 25 times for the first 16 non-rigid body modes, compared to the full order Abaqus model. Figure 6 presents the comparative analysis of reference and initial modal parameters, highlighting significant differences in both eigenfrequencies and eigenmodes. These variations are entirely attributed to the differences in the selected material properties.

Figure 6. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).

For this analysis, 150000 soundboard samples within the seven-dimensional parameter space are initially evaluated, using the geometric parameters of SB₂ ^init as nominal values. A significantly wider parameter space is considered here compared to the previous verification process, as no prior information about the global or local minima is available. Additionally, selecting starting points from this broader space enables the investigation of whether entirely different and widely distributed parameter sets can converge to the same unique or nearly optimal solution. Subsequently, optimization is performed using the best initial parameter set to obtain p ^* for SB₂ ^init. Out of the 10 starting points, only one yielded an acceptable optimal solution with an improved guitar soundboard shape after 346 modal evaluations. Two other parameter sets also converged to distinct optima; however, the resulting soundboard shapes were suboptimal. In these cases, approximately 216 modal evaluations were performed. All other starting points failed to produce a satisfactory distribution of modal parameters that matched the reference. However, it cannot be definitively concluded that the presented optimal design corresponds to the global minimum rather than a local minimum.

On the bottom in Figure 7 displays the geometric comparison between selected optimized soundboard design SB₂^opt and reference soundboard design SB₂^ref, while on the top and middle in Figure 7 illustrate the modal parameter comparison. The eigenfrequencies of the optimized soundboard closely match those of the reference soundboard. All the eigenmode shapes of SB₂^opt are in the correct order compared to SB₂^ref, with the average MAC value on the diagonal being ${\overline{\text{MAC}}}_{\mathit{mm}}=0.972$. It is clearly visible that the soundboard geometry is changed while the modal parameters is now very similar.

Figure 7. Comparison between reference and optimized design: eigenfrequencies (top), eigenmodes (middle), and geometry (bottom).

Figure 8 illustrates the relative eigenfrequency errors across the first 16 modes for initial and optimized cases. The initial average relative error of the eigenfrequencies is $\overline{|{\epsilon }^{\text{rel}}|}=17.81\%=340\text{cent}$. After the optimization process, this has improved to the average relative error of $\overline{|{\epsilon }^{\text{rel}}|}=1.78\%=31\text{cent}$, equivalent to a music interval close to unison. However, the third mode eigenfrequency shows the largest relative error of −6.41%, which corresponds to the cent value slightly grater than the minor second (one semitone).

Figure 8. Relative eigenfrequency errors, in percentage and cent units, between initial (SB2 init) and optimized (SB2 opt) soundboards relative to the reference (SB2 ref) soundboard.

7 Summary

This article presents a validated novel methodology for replicating the vibrational characteristics of a guitar soundboard through a parameterized reduced-order shell element model in the context of geometric shape optimization. By evaluating and comparing the first 16 non-rigid body mode shapes and their corresponding frequencies, the utilized objective function effectively identifies the optimal geometric design of a soundboard that aligns with the desired vibrational characteristics. The integration of the mesh morphing method enables the mapping of different global parameter settings into element-level parameters while preserving an identical mesh topology, thereby facilitating direct comparisons of mode shapes alongside their frequencies. Furthermore, the developed PMOR approach significantly reduces the computational cost associated with iterative model analyses during optimization. The integration of such a wide range of advanced techniques in computational luthiery constitutes a novel methodology, representing a significant step forward in enabling more accurate, efficient, and systematic soundboard shape optimization. While the model reduction technique significantly reduces computational cost on standard hardware, the optimization process is still constrained by the sequential nature of Matlab’s fmincon interior-point algorithm. Although the parameter space is inherently parallelizable, the core solver does not support native GPU or distributed acceleration. In this work, partial parallelization has already been employed for objective function evaluations, yielding measurable performance improvements. Future work may explore alternative optimization frameworks or customized implementations that enable full parallelism, both at the evaluation and algorithmic level, allowing better exploitation of modern hardware such as GPUs or computing clusters to further enhance scalability.

The verification results demonstrate that the implemented parameterized reduced-order shell element model and employed geometrical shape optimization process reliably converge to the global optimum, achieving the same modal parameters. This approach enables the identification of appropriate soundboard geometries to achieve predefined tonal characteristics while also compensating material variability. Additionally, the geometric shape optimization of a soundboard using two different wood species results in a new design that closely matches the targeted modal parameters.

Since the applied geometric parameter range in this work is sufficient to yield acceptable results, a broader parameter space can be explored in future work. Moreover, additional parameter constraints could be introduced to account for manufacturability and ergonomic considerations. The logical next step involves experimental validation of the findings by fabricating two soundboards with identical tonal characteristics, constructed from two different wood species. Furthermore, this methodology could be extended by incorporating appropriate bracing patterns, enabling further geometric optimization by considering variations in brace thickness, lengths, and orientations.

The proposed procedure is applicable to a wide range of string instruments and potentially beyond. From the perspective of musical instrument manufacturing, the findings of this contribution offer significant value. Moreover, this method provides an efficient and scientific approach to guitar design, moving beyond the traditional intuition-based practices. However, consideration of radiation efficiency and psychoacoustic factors is essential for a complete acoustic characterization. Future work will need to extend this approach by incorporating acoustic radiation modeling and perceptual evaluation, enabling a more comprehensive design strategy. Although possessing a similar vibrational modal spectrum alone does not guarantee identical acoustic output, this approach nevertheless brings us one step closer to enabling musicians to experience the timeless sound of historically significant instruments, granting access to the pinnacle of musical craftsmanship previously reserved for a privileged few.

Funding

The authors received no funding to complete this research.

Conicts of interest

The authors declare no conflict of interest.

Data availability statement

Data are available on request from the authors.

Author contribution statement

Tharindu Danushka Nandalal: Methodology, software, data curation, visualization, investigation, writing – original draft; Pierfrancesco Cillo: Conceptualization, methodology, software, writing – review and editing, supervision; Pascal Ziegler: Supervision, funding acquisition, writing – review and editing; Peter Eberhard: Supervision, project administration, funding acquisition, writing – review and editing.

References

All Tables

All Figures

	Figure 1. Geometrically parameterized shell element.
In the text

	Figure 2. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).
In the text

	Figure 3. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).
In the text

	Figure 4. Geometrically parameterized shell element.
In the text

	Figure 5. Correlation coefficient of 12 geometric parameters and predefined eigenmodes (top). Correlation coefficient ranking scaled by the 2-norm over predefined eigenmodes (bottom).
In the text

	Figure 6. Comparison between reference and initial design: eigenfrequencies (top), eigenmodes (bottom).
In the text

	Figure 7. Comparison between reference and optimized design: eigenfrequencies (top), eigenmodes (middle), and geometry (bottom).
In the text

	Figure 8. Relative eigenfrequency errors, in percentage and cent units, between initial (SB2 init) and optimized (SB2 opt) soundboards relative to the reference (SB2 ref) soundboard.
In the text